#include <misc.h>#include <preproc.h>module RtmMod#if (defined RTM) !----------------------------------------------------------------------- ! ! Purpose: ! River Routing Model! Contains routines: Rtmgridini, Rtmlandini, Rtmriverflux, Rtm! ! Method: ! (U. of Texas River Transport Model) !! Author: Sam Levis! !-----------------------------------------------------------------------! $Id: RtmMod.F90,v 1.10.6.5 2002/04/27 15:38:55 erik Exp $!-----------------------------------------------------------------------  use precision  use clm_varpar, only : lsmlon, lsmlat, rtmlon, rtmlat   implicit none! RTM grid info  integer , private :: numlon_r(rtmlat)               !number of lon points at each lat  real(r8), private, dimension(4) :: rtmedge = (/ 90., 180., -90., -180. /)  !N,E,S,W edges of rtm grid  real(r8), public, allocatable :: latixy_r(:,:)      !rtm latitudes  of grid cells (degrees)         real(r8), public, allocatable :: longxy_r(:,:)      !rtm longitudes of grid cells (degrees)         real(r8), public :: area_r(rtmlon,rtmlat)           !rtm gridcell area (km^2)            integer , public :: mask_r(rtmlon,rtmlat)           !rtm landmask (land=1,ocean=0)! land model to RTM mapping. for each rtm grid cell:  integer , private :: mxovr_s2r                      !max number of overlapping cells  integer , private :: novr_s2r(rtmlon,rtmlat)        !number    of overlapping cells  integer , private, allocatable :: iovr_s2r(:,:,:)   !lon index of overlapping cells  integer , private, allocatable :: jovr_s2r(:,:,:)   !lat index of overlapping cells  real(r8), private, allocatable :: wovr_s2r(:,:,:)   !weight    of overlapping cells! RTM runoff for coupled communication  integer , public, allocatable :: ocnrof_iindx(:)    !rtm longitude index of ocean runoff point  integer , public, allocatable :: ocnrof_jindx(:)    !rtm latitude index of ocean runoff point  real(r8), public, allocatable :: ocnrof_vec(:)      !rtm runoff vector (1/2 deg grid, kg/m^2/s)! RTM history file variables  real(r8), public, allocatable :: qchan2(:)          !river (channel) flow (m**3 H2O /s)  real(r8), public, allocatable :: qchocn2(:)         !river (channel) flow into ocean (m**3/s)! time averaging for rtm calculatino  real(r8), public, allocatable :: totrunin_ave(:)    !time averaged vector of input fluxes  real(r8), public, allocatable :: prec_ave(:)        !time averaged vector of precipitation  real(r8), public, allocatable :: evap_ave(:)        !time averaged vector of evaporation  real(r8), public :: delt_rtm                        !rtm time step  integer , public :: ncount_rtm                      !number of time samples to average over! fluxes  integer , private :: rdirc(0:rtmlon+1,0:rtmlat+1)   !rtm river flow direction (0-8)  real(r8), private :: fluxout(0:rtmlon+1,0:rtmlat+1) !water flux out of cell (m^3/s)  real(r8), private :: ddist(rtmlon,rtmlat)           !downstream distance (m)  real(r8), private :: rivarea(rtmlon,rtmlat)         !cell area (m^2)  real(r8), public  :: volr(rtmlon,rtmlat)            !water volume in cell (m^3)  real(r8), private, allocatable :: latsh(:)          !southern edge of cells at rtm grid       real(r8), private, allocatable :: lonwh(:,:)        !western  edge of cells at rtm grid     ! inputs to RTM at 1/2 degree resolution  real(r8), private :: totrunin_r(rtmlon,rtmlat)      !surface runoff (mm/s)! outputs returned from RTM at 1/2 degree resolution  real(r8), private :: flxlnd_r(rtmlon,rtmlat)        !river flux (m**3/s)  real(r8), private :: flxocn_r(rtmlon,rtmlat)        !river flux to the ocean (m**3/s)  real(r8), private :: dvolrdt_r(rtmlon,rtmlat)       !change in storage (mm/s)  real(r8), private :: volrtm(rtmlon,rtmlat)          !change in storage (m**3/s)  real(r8), private :: runrtm(rtmlon,rtmlat)          !input runoff on rtm grid (m**3/s)! RTM water flux into cell  real(r8), private :: sfluxin(rtmlon,rtmlat)         !water flux into cell (m3/s)! global averaging  character(len=*),parameter :: F40="('(diag) ',a17,'    date  ', &  &     '   prec        evap        runoff(lnd)   runoff(rtm) dvoldt(rtm) runoff-ocn(rtm)  (m^3/sec)')"  character(len=*),parameter :: F41="('(diag) ',a17,'   nstep  ', &  &     '   prec        evap        runoff(lnd)   runoff(rtm) dvoldt(rtm) runoff-ocn(rtm)  (m^3/sec)')"  character(len=*),parameter :: F21="('(diag) ',a17,' ----------------------', &  &     7('----------'))"  character(len=*),parameter :: F22="('(diag) ',a17,i8,6(d13.4))"  real(r8) prec_global            !total precipitation (m^3/sec)   real(r8) evap_global            !total evaporation (m^3/sec)  real(r8) runlnd_global          !total input runoff on land grid (m^3/sec)  real(r8) runrtm_global          !total input runoff on rtm grid (m^3/sec)  real(r8) ocnrtm_global          !total ocean runoff on rtm grid (m^3/sec)  real(r8) volrtm_global          !total change in storage on rtm (m^3/sec)  integer  ncount_global          !global counter   integer  yrold                  !old year  SAVE!=======================================================================CONTAINS!=======================================================================subroutine Rtmgridini!----------------------------------------------------------------------- ! ! Purpose: ! Initialize RTM grid and land mask (U. of Texas River Transport Model)! ! Method: ! ! Author: Sam Levis! !-----------------------------------------------------------------------  use spmdMod   , only : masterproc  use areaMod   , only : celledge, cellarea     use clm_varctl, only : frivinp_rtm  use clm_varcon, only : re  use shr_const_mod, only: SHR_CONST_PI! ------------------------ local variables ---------------------------  integer  :: ioff(0:8) = (/0,0,1,1,1,0,-1,-1,-1/) !calc dist as in hydra  integer  :: joff(0:8) = (/0,1,1,0,-1,-1,-1,0,1/) !of grid cell down stream  integer  :: i,j,k,n                       !loop indices  integer  :: i2,j2                         !downstream i and j  real(r8) :: deg2rad                       !pi/180  real(r8) :: dx                            !lon dist. between grid cells (m)  real(r8) :: dy                            !lat dist. between grid cells (m)  real(r8) :: dist(rtmlon,rtmlat)           !dist. of the grid cell down stream (m)   real(r8) :: tempg(rtmlon,rtmlat)          !temporary buffer  integer  :: tempgp(0:rtmlon+1,0:rtmlat+1) !temporary buffer  ! --------------------------------------------------------------------    if (masterproc) then! --------------------------------------------------------------------! Useful constants and initial values! --------------------------------------------------------------------     write(6,*)'Columns in RTM = ',rtmlon     write(6,*)'Rows in RTM    = ',rtmlat     allocate(latixy_r(rtmlon,rtmlat))     allocate(longxy_r(rtmlon,rtmlat))     allocate(latsh(rtmlat+1))        !southern edge of cells at rtm grid          allocate(lonwh(rtmlon+1,rtmlat)) !western  edge of cells at rtm grid               deg2rad = SHR_CONST_PI / 180.     volr = 0.! --------------------------------------------------------------------! Open and read input data (river direction file)! rtm operates from south to north and from the dateline! --------------------------------------------------------------------     open (1,file=frivinp_rtm)     write(6,*)'opened river direction data'     do j = 1,rtmlat        numlon_r(j) = 0        do i = 1,rtmlon           read(1,*) latixy_r(i,j),longxy_r(i,j),tempg(i,j)           if (longxy_r(i,j) /= 1.e36) numlon_r(j) = numlon_r(j) + 1           tempgp(i,j) = nint(tempg(i,j))        enddo     enddo     close(1)     write(6,*)'closed river direction data'     write(6,*)     ! --------------------------------------------------------------------! Determine RTM celledges, areas and interpolation masks! --------------------------------------------------------------------     call celledge (rtmlat    , rtmlon    , numlon_r  , longxy_r  , &                    latixy_r  , rtmedge(1), rtmedge(2), rtmedge(3), &                    rtmedge(4), latsh     , lonwh     )     call cellarea (rtmlat    , rtmlon    , numlon_r  , latsh     , lonwh , &                    rtmedge(1), rtmedge(2), rtmedge(3), rtmedge(4), area_r) ! --------------------------------------------------------------------! Determine rtm mask, downstream distance and area! --------------------------------------------------------------------! determine rtm ocn/land mask     do i=1,rtmlon        tempgp(i,0)        = tempgp(mod(i+rtmlon/2-1,rtmlon)+1,1)        tempgp(i,rtmlat+1) = tempgp(mod(i+rtmlon/2-1,rtmlon)+1,rtmlat)        if (tempgp(i,0)        /= 0) tempgp(i,0)        = mod(tempgp(i,0)       +4-1,8)+1        if (tempgp(i,rtmlat+1) /= 0) tempgp(i,rtmlat+1) = mod(tempgp(i,rtmlat+1)+4-1,8)+1     enddo     do j=0,rtmlat+1        tempgp(0,j) =tempgp(rtmlon,j)        tempgp(rtmlon+1,j)=tempgp(1,j)     enddo          do j=0,rtmlat+1        do i=0,rtmlon+1           rdirc(i,j)=tempgp(i,j)        enddo     enddo          do j=1,rtmlat        do i=1,rtmlon           if (rdirc(i,j) == 0) then              mask_r(i,j) = 0           else              mask_r(i,j) = 1           end if        enddo     enddo! determine downstream distance - instead of reading a distance file ! calculate the downstream distance as in hydra     do j=1,rtmlat        do i=1,rtmlon           i2 = i + ioff(tempgp(i,j))           j2 = j + joff(tempgp(i,j))           if (i2 == 0) i2 = 2                 !avoids i2 out of bounds in the following           if (i2 == rtmlon+1) i2 = rtmlon-1   !avoids i2 out of bounds in the following             dy = deg2rad * abs(latixy_r(i,j)-latixy_r(i2,j2)) * re*1000.           dx = deg2rad * abs(longxy_r(i,j)-longxy_r(i2,j2)) * re*1000. &                *0.5*(cos(latixy_r(i,j)*deg2rad)+cos(latixy_r(i2,j2)*deg2rad))           dist(i,j) = sqrt(dx*dx + dy*dy)           ddist(i,j) = dist(i,j)           rivarea(i,j)=1.e6 * area_r(i,j)     !convert into m**2        enddo     enddo    endif  ! end of if-masterproc blockend subroutine Rtmgridini!=======================================================================subroutine Rtmlandini!----------------------------------------------------------------------- ! ! Purpose: ! Initialize RTM-land interpolation weights (U. of Texas River Transport Model)! and variables related to runoff time averaging! ! Method: ! ! Author: Mariana Vertenstein! !-----------------------------------------------------------------------  use spmdMod     , only : masterproc  use areaMod     , only : areaini              use clm_varsur  , only : numlon, area, lats, lonw, landmask   use clm_varmap  , only : numpatch  use time_manager, only : get_curr_date! ------------------------ local variables ---------------------------  integer  :: i,j,k,n                    !loop indices  integer  :: is,js                      !land model grid indices  real(r8) :: maskone_s(lsmlon,lsmlat)   !dummy field: see below                   real(r8) :: maskone_r(rtmlon,rtmlat)   !dummy field: see below                   integer  :: ocnrof_mask(rtmlon,rtmlat) !rtm mask for ocean points with possible nonzero runoff  integer  :: ocnrof_num                 !number of valid ocean points with possible nonzero runoff  integer  :: yrnew                      !year (0, ...)  integer  :: mon                        !month (1, ..., 12)   integer  :: day                        !day of month (1, ..., 31)  integer  :: ncsec                      !seconds of current date! --------------------------------------------------------------------  if (masterproc) then! --------------------------------------------------------------------! The following section allows RTM and land model to coexist at different! horizontal resolutions! --------------------------------------------------------------------     write(6,*)     write(6,*) 'Initializing area-averaging interpolation for RTM.....'! To find fraction of each land model grid cell that is land based on rtm grid.! For this purpose, want all rtm grid cells to contribute to grid cell ! average on land model grid, i.e., all cells used regardless of whether land ! or ocean. Do this by setting [maskone_s] = 1 ! [maskone_s] = 1 means all grid cells on land model grid, regardless of whether! land or ocean, will contribute to rtm grid.     do j = 1, lsmlat        do i = 1, numlon(j)