!############################################################################### !############################################################################### ! Code for runoff routing. !############################################################################### !############################################################################### ! subroutine route_drive ! Driver routine for runoff routing. SUBROUTINE route_drive( land_pts,sub_surf_roff,surf_roff,roffAccumLand ) USE route_mod, ONLY : & ! imported scalar parameters routeTypeTRIP & ! imported scalars with intent(in) ,routeTimestep,routeType & ! imported scalars with intent(inout) ,routeCount,routeStep USE timeConst, ONLY : & iSecInDay !------------------------------------------------------------------------------- IMPLICIT NONE ! Array arguments with intent(in) INTEGER, INTENT(in) :: land_pts ! number of land points ! Array arguments with intent(in) REAL, INTENT(in) :: sub_surf_roff(land_pts) ! Sub-surface runoff (kg m-2 s-1) REAL, INTENT(in) :: surf_roff(land_pts) ! Surface runoff (kg m-2 s-1) ! Array arguments with intent(inout) REAL, INTENT(inout) :: roffAccumLand(land_pts) ! average runoff (production) ! rate between calls to routing (kg m-2 s-1) !------------------------------------------------------------------------------- IF ( routeCount == 0 ) THEN SELECT CASE ( routeType ) CASE ( routeTypeTRIP ) ! Reset accumulated runoff to zero. roffAccumLand(:) = 0.0 CASE default WRITE(*,*)'ERROR: route_drive: do not recognise routeType=',TRIM(routeType) STOP END SELECT ENDIF ! Increment counters. routeStep = routeStep + 1 routeCount = routeCount + 1 !------------------------------------------------------------------------------- ! Accumulate runoff. !------------------------------------------------------------------------------- SELECT CASE ( routeType ) CASE ( routeTypeTRIP ) ! Accumulate the total runoff. roffAccumLand(:) = roffAccumLand(:) + ( surf_roff(:) + sub_surf_roff(:) ) & / REAL(routeTimeStep) CASE default WRITE(*,*)'ERROR: route_drive: do not recognise routeType=',TRIM(routeType) STOP END SELECT !------------------------------------------------------------------------------- ! Decide if routing model is to be called on this timestep and call. !------------------------------------------------------------------------------- IF ( routeStep >= routeTimeStep ) THEN !------------------------------------------------------------------------------- ! Routing is to be called. !------------------------------------------------------------------------------- SELECT CASE ( routeType ) CASE ( routeTypeTRIP ) !------------------------------------------------------------------------------- ! Correct average runoff rate if did not get expected number of increments. ! This deals with awkward cases - e.g. when spin-up interupts. !------------------------------------------------------------------------------- IF ( routeCount /= routeTimeStep ) THEN IF ( routeCount > 0 ) THEN roffAccumLand(:) = roffAccumLand(:) * REAL(routeTimeStep) & / REAL(routeCount) ELSE ! At present I'm not sure if we'd ever get here! But just in case... WRITE(*,*) 'ERROR: route_drive div by routeCount=0' STOP ENDIF ENDIF !------------------------------------------------------------------------------- ! Call TRIP routing driver. !------------------------------------------------------------------------------- CALL route_drive_trip( roffAccumLand ) CASE default !------------------------------------------------------------------------------- ! Default case for routeType. !------------------------------------------------------------------------------- WRITE(*,*)'ERROR: route_drive: do not recognise routeType=',TRIM(routeType) STOP END SELECT !------------------------------------------------------------------------------- ! Reset counters after a call to routing. !------------------------------------------------------------------------------- routeStep = 0 routeCount = 0 ! ELSE ! Routing is not done on this timestep. ! Set routing diagnostics to missing data. Prognostics will retain current value. ! select case ( routeType ) ! case ( routeTypeTRIP ) ! if ( outPosRFlow > 0 ) outValExtra(outPosRFlow:outPosRFlow+npRoute-1) = undefOut ! if ( outPosRRun > 0 ) outValExtra(outPosRRun:outPosRRun+npRoute-1) = undefOut ! case default ! write(*,*)'ERROR: route_drive: do not recognise routeType=',trim(routeType) ! stop ! end select ENDIF ! routeStep END SUBROUTINE route_drive !############################################################################### !############################################################################### !############################################################################### !############################################################################### ! subroutine route_drive_trip ! Driver routine for runoff routing by the TRIP (linear) model. SUBROUTINE route_drive_trip( roffAccumLand ) !------------------------------------------------------------------------------- ! Modules used: USE ancil_info, ONLY : & ! imported scalars with intent(in) land_pts,nx=>row_length,ny=>rows & ! imported arrays with intent(in) ,land_index USE earth_utils, ONLY : & ! imported procdures earthArea USE grid_utils, ONLY : & ! imported procedures getXYpos USE offline_diag, ONLY : & ! imported scalars with intent(in) useRflowDiag,useRrunDiag & ! imported arrays with intent(inout) ,rflowDiag,rrunDiag USE route_mod, ONLY : & ! imported scalars with intent(in) npRoute,nxRoute,nyRoute,routeDlat,routeDlon,routeLat1,routeLon1 & ,routeRegLatLon,routeRegrid & ! imported arrays with intent(in) ,routeIndex & ! imported arrays with intent(inout) ,routeNext USE time_loc, ONLY : & ! imported scalars with intent(in) nxGrid,nyGrid,regDlat,regDlon,regLat1,regLon1 & ! imported arrays with intent(in) ,latitude,longitude !------------------------------------------------------------------------------- IMPLICIT NONE ! Array arguments with intent(in) REAL, INTENT(in) :: roffAccumLand(land_pts) ! average rate of runoff since ! last routing call (kg m-2 s-1) ! Local scalar variables. INTEGER :: i ! work INTEGER :: ip ! loop counter INTEGER :: ix ! work INTEGER :: iy ! work INTEGER :: maxLin ! an array size ! Local array variables. INTEGER :: x(npRoute) ! x index on the routing grid of active routing points INTEGER :: xNext(npRoute) ! x location of the next downstream point ! <0 indicates next point is off grid INTEGER :: y(npRoute) ! y index on the routing grid of active routing points INTEGER :: yNext(npRoute) ! y location of the next downstream point ! If xNext(i)<0, yNext(i) gives the direction of the ! next point, in terms of index in flowDirSet. REAL :: outFlow(nxRoute,nyRoute) ! rate of channel flow leaving gridbox (kg s-1) REAL :: routeArea(nxRoute,nyRoute) ! area of each gridbox on routing grid (m2) REAL :: routeLat(nxRoute,nyRoute) ! latitude of points on routing grid (degrees) REAL :: routeLength(nxRoute,nyRoute) ! distance to next downstream gridpoint (m) ! This is only calculated at "active" routing points. REAL :: routeLon(nxRoute,nyRoute) ! longitude of points on routing grid (degrees) REAL :: runoff(nxRoute,nyRoute) ! average rate of runoff, on routing grid (kg s-1) ! During work stages, units may be kg m-2 s-1. REAL :: roffAccumGrid(nxGrid,nyGrid) ! average rate of runoff since last routing ! call (kg m-2 s-1), on a grid that ! is a superset of (can be equal to) the model grid LOGICAL :: routeMask(nxRoute,nyRoute) ! TRUE at land points on routing grid, ! else FALSE !------------------------------------------------------------------------------- ! Curent code can only deal with "regular" lat/lon grids. !------------------------------------------------------------------------------- IF ( .NOT. routeRegLatLon ) THEN WRITE(*,*)'ERROR: route_drive_trip: code only exists for regular lat/lon grids.' WRITE(*,*)'Sorry. Maybe you''d like to write some new code...?' STOP ENDIF !------------------------------------------------------------------------------- ! Calculate location on routing grid of each land point, and find the point ! immediately downstream. Also set a mask to show land points on the routing ! grid. !------------------------------------------------------------------------------- ! Initialise. xNext(:) = 0 yNext(:) = 0 routeMask(:,:) = .FALSE. DO ip=1,npRoute ! Get location of this point. CALL getXYPos( routeIndex(ip),nxRoute,nyRoute,x(ip),y(ip) ) ! Set mask to TRUE at this land point. routeMask(x(ip),y(ip)) = .TRUE. ! Get location of next downstream point. IF ( routeNext(ip) > 0 ) THEN ! Downstream point is on the grid. CALL getXYPos( routeNext(ip),nxRoute,nyRoute,xNext(ip),yNext(ip) ) ELSEIF ( routeNext(ip) == -9 ) THEN ! There is no defined flow direction, and so no downstream point. xNext(ip) = -9 yNext(ip) = -9 ELSE ! Downstream point is off edge of grid. ! routeNext has been set to -1*index of direction in flowDirSet. xNext(ip) = -1 * routeNext(ip) yNext(ip) = xNext(ip) ENDIF ENDDO !------------------------------------------------------------------------------- ! Calculate latitude and longitude of each gridpoint on routing grid. !------------------------------------------------------------------------------- DO iy=1,nyRoute DO ix=1,nxRoute routeLat(ix,iy) = routeLat1 + (iy-1)*routeDlat routeLon(ix,iy) = routeLon1 + (ix-1)*routeDlon ENDDO ENDDO !------------------------------------------------------------------------------- ! Calculate distance between grid points. !------------------------------------------------------------------------------- CALL getRouteLen( x,xNext,y,yNext,routeLat,routeLon,routeLength ) !------------------------------------------------------------------------------- ! Calculate area of each routing gridbox. !------------------------------------------------------------------------------- DO iy=1,nyRoute DO ix=1,nxRoute routeArea(ix,iy) = earthArea( & routeLat(ix,iy)-0.5*routeDlat,routeLat(ix,iy)+0.5*routeDlat & ,routeLon(ix,iy)-0.5*routeDlon,routeLon(ix,iy)+0.5*routeDlon ) * 1.0e6 ENDDO ENDDO !------------------------------------------------------------------------------- ! Put accumulated runoff onto a grid. ! If the "main" model grid is 2-D, this target grid is the 2-D grid. ! If the "main" grid is a vector (in offline applications of JULES this is ! possible if points from a larger grid have been compressed - e.g. land points ! selected from a larger grid.), the target grid is the larger grid, across ! which the points are to be scattered. ! ! Arguably this would be better done at a higher level, since it's unlikely that ! the requirement to grid the runoff will be unique to TRIP. But it's here for ! now. An advantage of being here is that we only grid the fields we need for ! this particular model (if done higher, we would have to test which fields to ! grid for which model - not impossible!). ! ! Note that the regridded runoff may be non-zero at TRIP sea points. ! !------------------------------------------------------------------------------- roffAccumGrid(:,:) = 0.0 IF ( nxGrid /= nx .OR. nyGrid /= ny ) THEN ! The "main" model grid is not a regular lat/lon grid, but is a subset of ! such a regular grid. ! We cannot use land_index to scatter the points across this grid. Instead, ! use lat/lon to locate in the larger grid of shape (nxGrid,nyGrid). ! Note that this code will not currently work with grids that are just not ! a subset of a "regular" grid - e.g. rotated grid - since then nxGrid and ! nyGrid have not been given useful values.... DO i=1,land_pts ip = land_index(i) ix = NINT( (longitude(ip,1) - regLon1) / regDlon ) + 1 iy = NINT( (latitude(ip,1) - regLat1) / regDlat ) + 1 roffAccumGrid(ix,iy) = roffAccumLand(ip) ENDDO ELSE ! The "main" model grid is a regular, lat/lon grid. ! land_index can be used to scatter the points across this grid. DO ip=1,land_pts CALL getXYPos( land_index(ip),nx,ny,ix,iy ) roffAccumGrid(ix,iy) = roffAccumLand(ip) ENDDO ENDIF !------------------------------------------------------------------------------- ! Get runoff from land grid onto routing grid, and convert from flux density ! (kg m-2 s-1) to flux (kg s-1). !------------------------------------------------------------------------------- runoff(:,:) = 0.0 IF ( .NOT. routeRegrid ) THEN !------------------------------------------------------------------------------- ! "Main" and routing grids are identical (including points in same order). ! No need to regrid, simply convert units. !------------------------------------------------------------------------------- runoff(:,:) = roffAccumGrid(:,:) * routeArea(:,:) ELSE !------------------------------------------------------------------------------- ! Regrid from land to routing grid (and convert to flux). !------------------------------------------------------------------------------- ! Calculate an array size used in regridding. maxLin = ( nx + nxRoute )*( ny + nyRoute ) CALL route_regrid( maxLin,roffAccumGrid,runoff ) ! Units of runoff are currently those of roffAccumGrid *(kg m-2 s-1). ! Convert from flux density to flux. runoff(:,:) = runoff(:,:) * routeArea(:,:) ENDIF !------------------------------------------------------------------------------- ! Call the routing routine. !------------------------------------------------------------------------------- CALL route_trip( routeLength,runoff,outFlow ) !------------------------------------------------------------------------------- ! Add runoff on TRIP sea points to the river flow variable. ! Regridding may have resulted in some land runoff appearing in TRIP sea ! gridboxes, so we need to account for this water. Even without regridding, ! a land point that is sea in TRIP will not have been added to flow. !------------------------------------------------------------------------------- WHERE( .NOT. routeMask(:,:) ) outFlow(:,:) = runoff(:,:) ! do iy=1,nyRoute ! do ix=1,nxRoute ! if ( .NOT. routeMask(ix,iy) ) outFlow(ix,iy) = runoff(ix,iy) ! enddo ! enddo !------------------------------------------------------------------------------- ! Save any extra diagnostics that are requested for offline application. !------------------------------------------------------------------------------- IF ( useRflowDiag ) rFlowDiag(:,:) = outFlow(:,:) IF ( useRrunDiag ) rRunDiag(:,:) = runoff(:,:) END SUBROUTINE route_drive_trip !############################################################################### !############################################################################### !############################################################################### ! subroutine getRouteLen ! Driver routine that calls procedures that calculates distance between gridpoints. SUBROUTINE getRouteLen( x,xNext,y,yNext,routeLat,routeLon,length ) USE earth_utils, ONLY : & ! imported procedures giveLength USE route_mod, ONLY : & ! imported scalars with intent(in) npRoute,nxRoute,nyRoute,routeDlat,routeDlon & ! imported arrays with intent(in) ,flowDirDelta,flowDirSet USE time_loc, ONLY : & ! imported arrays with intent(in) latitude,longitude !------------------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------------------- ! Array arguemnts with intent(in) !------------------------------------------------------------------------------- INTEGER, INTENT(in) :: & x(npRoute) &! x index on the routing grid of active routing points ,xNext(npRoute) &! x location of the next downstream point ! <0 indicates next point is off grid and gives the ! direction of the next point, in terms of index in ! flowDirSet. ,y(npRoute) &! y index on the routing grid of active routing points ,yNext(npRoute) ! y location of the next downstream point ! <0 interpreted as for xNext. REAL, INTENT(in) :: & routeLat(nxRoute,nyRoute) &! latitude of points on routing grid (degrees) ,routeLon(nxRoute,nyRoute) ! longitude of points on routing grid (degrees) !------------------------------------------------------------------------------- ! Array arguments with intent(out) !------------------------------------------------------------------------------- REAL, INTENT(out) :: & length(nxRoute,nyRoute) ! distance between each gridpoint and downstream gridpoint (m) !------------------------------------------------------------------------------- ! Local scalar variables. !------------------------------------------------------------------------------- INTEGER :: & ip,ix,iy &! work ,jx,jy ! work REAL :: & dx,dy ! work !------------------------------------------------------------------------------- ! Only need to calculate length at land points on routing grid. DO ip=1,npRoute !------------------------------------------------------------------------------- ! Get coords of this point and any point immediately downstream. !------------------------------------------------------------------------------- ix = x(ip) iy = y(ip) jx = xNext(ip) jy = yNext(ip) IF ( jx > 0 ) THEN !------------------------------------------------------------------------------- ! There is an outflow direction, and the downstream point is on the grid. !------------------------------------------------------------------------------- length(ix,iy) = giveLength( routeLat(ix,iy),routeLon(ix,iy) & ,routeLat(jx,jy),routeLon(jx,jy) ) * 1000.0 ELSEIF ( jx == -9 ) THEN !------------------------------------------------------------------------------- ! Outflow to sea, or no outflow (pit or depression). ! This is relatively common, hence treatment separately from flow across grid edge. ! Use distance to an adjacent point at same latitude. !------------------------------------------------------------------------------- length(ix,iy) = giveLength( routeLat(ix,iy),routeLon(ix,iy) & ,routeLat(ix,iy),routeLon(ix,iy)+routeDlon ) * 1000.0 ELSE !------------------------------------------------------------------------------- ! Outflow across the edge of the grid. We would still like to calculate a length, ! so that, as far as possible, a regional implementation gives the same results ! as a global implementation, even for these edge points. ! NB This works for outflow from grid, but if a regional run is missing an ! inflow from outside the grid, answers will generally be different in global ! and regional implementations. ! The lat/lon of the downstream point ! is not available in these cases, so we ASSUME a "regular" lat/lon grid. ! jx is -1 * index in flowDirSet. !------------------------------------------------------------------------------- ! Get number of gridboxes in each direction to the downstream location. dx = REAL( flowDirDelta(ABS(jx),1) ) dy = REAL( flowDirDelta(ABS(jx),2) ) length(ix,iy) = giveLength( routeLat(ix,iy),routeLon(ix,iy) & ,routeLat(ix,iy)+dy*routeDlat & ,routeLon(ix,iy)+dx*routeDlon ) * 1000.0 ENDIF ENDDO END SUBROUTINE getRouteLen !############################################################################### !############################################################################### !############################################################################### ! subroutine route_trip ! Calculate river outflow and update channel storage for TRIP routing model. ! Based on UM routine outflow1. ! At present there is no regridding and therefore no need to account for water ! that may have been allocated to e.g. a sea point during regridding. SUBROUTINE route_trip( routeLength,runoff,outFlow ) USE grid_utils, ONLY : & ! imported procedures getXYpos USE prognostics, ONLY : & ! imported arrays with intent(inout) routeStore USE route_mod, ONLY : & ! imported scalars with intent(in) npRoute,nxRoute,nyRoute,routeMeander,routeSpeed,routeTimeStep & ! imported arrays with intent(in) ,routeIndex,routeNext,routeOrder USE time_loc, ONLY : & ! imported scalars with intent(in) timeStep !------------------------------------------------------------------------------- IMPLICIT NONE INTEGER :: &! local scalars i,ip,ix,iy,jx,jy ! work/loop counters REAL, INTENT(in) :: &! in arrays routeLength(nxRoute,nyRoute) &! distance between gridpoints (m) ,runoff(nxRoute,nyRoute) ! rate of runoff generation in each gridbox (kg s-1) REAL, INTENT(out) :: &! out arrays outFlow(nxRoute,nyRoute) ! rate of channel flow leaving gridbox (kg s-1) REAL :: &! local scalars coeff &! coefficient in the routing model (s-1) ,dt &! timestep of routing model (s) ,storeOld ! channel storage (kg) REAL :: &! local arrays inflow(nxRoute,nyRoute) ! rate of channel flow entering gridbox (kg s-1) !------------------------------------------------------------------------------- dt = REAL(routeTimeStep) * timeStep ! Initialise inflow with runoff generated over each gridbox. inflow(:,:) = runoff(:,:) ! Initialise outlow (for clarity at non-routing points). outFlow(:,:) = 0.0 ! Loop over active routing points. DO i=1,npRoute ! Get index (location in routing vector) of the point to consider. ip = routeOrder(i) ! Get location of this point in routing grid. CALL getXyPos( routeIndex(ip),nxRoute,nyRoute,ix,iy ) ! Calculate the coefficient "c" of the model. ! c=u/(d*r), where u is effective flow speed, d is distance between gridpoints, ! and r is meander ratio. coeff = routeSpeed / ( routeLength(ix,iy) * routeMeander ) ! Save value of channel storage at start of timestep. storeOld = routeStore(ix,iy) ! Calculate channel storage at end of timestep. ! Eqn.4 of Oki et al, 1999, J.Met.Soc.Japan, 77, 235-255. routeStore(ix,iy) = storeOld * EXP(-(coeff*dt)) & + (1.0 - EXP(-(coeff*dt))) * inflow(ix,iy) / coeff ! Calculate outflow as inflow minus change in storage. outflow(ix,iy) = inflow(ix,iy) + (storeOld-routeStore(ix,iy)) / dt ! Add outflow to inflow of next downstream point. IF ( routeNext(ip) > 0 ) THEN ! Get location in grid of next downstream point. CALL getXYPos( routeNext(ip),nxRoute,nyRoute,jx,jy ) inflow(jx,jy) = inflow(jx,jy) + outflow(ix,iy) ENDIF ENDDO ! points END SUBROUTINE route_trip !############################################################################### !############################################################################### !############################################################################### ! subroutine route_regrid ! Handles regridding of runoff from "main" to routing grids. SUBROUTINE route_regrid( maxLin,roffAccumGrid,runoff ) !------------------------------------------------------------------------------- ! Description: ! Regrids runoff from a source grid to a target grid (the routing grid). ! Both grids must be regular in latitude and longitude. ! !------------------------------------------------------------------------------- ! Modules used: USE route_mod, ONLY : & ! imported scalars with intent(in) nxRoute,nyRoute,routeDlat,routeDlon,routeLat1,routeLon1 USE time_loc, ONLY : & ! imported scalars with intent(in) nxGrid,nyGrid,regDlat,regDlon,regLat1,regLon1 !------------------------------------------------------------------------------- IMPLICIT NONE ! Scalar arguments with intent(in) INTEGER, INTENT(in) :: maxLin ! a vector length ! Array arguments with intent(in) REAL, intent(in) :: roffAccumGrid(nxGrid,nyGrid) ! runoff rate on model grid (kg m-2 s-1) ! Array arguments with intent(out) REAL, intent(out) :: runoff(nxRoute,nyRoute) ! runoff rate on routing grid (kg m-2 s-1) !------------------------------------------------------------------------------- ! Local scalar variables. INTEGER :: adjust ! mode (option number) of area averaging INTEGER :: icode ! exit code from subroutines INTEGER :: ix ! loop counter INTEGER :: iy ! loop counter INTEGER :: maxL ! number of entries in output from pre_areaver LOGICAL :: cyclic_srce ! TRUE source (model) grid is cyclic in x LOGICAL :: cyclic_targ ! TRUE target (model) grid is cyclic in x LOGICAL :: invert_srce ! TRUE source (model) grid runs N to S ! FALSE source grid runs S to N LOGICAL :: spherical ! TRUE coordinates are lat/lon on a sphere ! FALSE Cartesian axes. LOGICAL :: want ! Value of masks at locations for which data are required CHARACTER(len=80) :: cmessage ! error message from subroutines !------------------------------------------------------------------------------- ! Local array variables. INTEGER :: count_targ(nxRoute,nyRoute) ! number of model gridboxes that ! contribute to each routing gridbox INTEGER :: base_targ(nxRoute,nyRoute) ! base (starting) index for each ! routing gridbox (i.e. location of first element in lists) INTEGER :: index_srce(maxLin) ! list of source (model) gridboxes ! that contribute to each routing gridbox REAL :: sourceLat(nyGrid+1) ! latitudes of edges of source (model) gridboxes. ! First value is the S edge of first gridbox, all other values are the ! N edge of each gridbox. REAL :: sourceLon(nxGrid+1) ! longitudes of edges of source (model) gridboxes ! First value is the W edge of first gridbox, all other values are the ! E edge of each gridbox. REAL :: targetLat(nyRoute+1) ! latitudes of edges of target (mrouting) gridboxes. ! First value is the S edge of first gridbox, all other values are the ! N edge of each gridbox. REAL :: targetLon(nxRoute+1) ! longitudes of edges of target (routing) gridboxes ! First value is the W edge of first gridbox, all other values are the ! E edge of each gridbox. REAL :: adjust_targ(nxRoute,nyRoute) ! adjustment factors (not used) REAL :: weight(maxLin) ! lists of weights for each source (model) gridbox LOGICAL :: mask_srce(nxGrid,nyGrid) LOGICAL :: mask_targ(nxRoute,nyRoute) !------------------------------------------------------------------------------- ! Set values. spherical = .TRUE. ! Calculations are for lat/lon coordinates. adjust = 0 ! "normal" area averaging maxL = maxlIn invert_srce = .FALSE. ! model grid runs S to N ! Set coordinates of edges of model gridboxes. DO ix=1,nxGrid+1 sourceLon(ix) = regLon1 + (REAL(ix)-0.5) * regDlon ENDDO DO iy=1,nyGrid+1 sourceLat(iy) = regLat1 + (REAL(iy)-0.5) * regDlat ENDDO ! Set coordinates of edges of routing gridboxes. DO ix=1,nxRoute+1 targetLon(ix) = routeLon1 + (REAL(ix)-0.5) * routeDlon ENDDO DO iy=1,nyRoute+1 targetLat(iy) = routeLat1 + (REAL(iy)-0.5) * routeDlat ENDDO ! Decide if grids are cyclic in longitude. cyclic_srce = .FALSE. cyclic_targ = .FALSE. IF ( REAL(nxGrid)*regDlon > 359.9 ) cyclic_srce = .TRUE. IF ( REAL(nxRoute)*routeDlon > 359.9 ) cyclic_targ = .TRUE. ! Set masks to indicate that all points in both grids are to be used. want = .TRUE. mask_srce(:,:) = want mask_targ(:,:) = want ! Call setup routing for averaging. CALL PRE_AREAVER( nxGrid,sourceLon,nyGrid,sourceLat & ,cyclic_srce,nxGrid,want,mask_srce & ,nxRoute,targetLon,nyRoute,targetLat & ,cyclic_targ,spherical & ,maxL,count_targ,base_targ,index_srce,weight,icode,cmessage ) ! Call averaging routine. CALL DO_AREAVER( nxGrid,nyGrid,nxGrid & ,invert_srce,roffAccumGrid,nxRoute,nyRoute,count_targ & ,base_targ,nxRoute,want,mask_targ,index_srce,weight,adjust & ,runoff,adjust_targ,icode,cmessage ) END SUBROUTINE route_regrid !############################################################################### !############################################################################### !###############################################################################