rtmmod.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 910 行 · 第 1/3 页

           maskone_s(i,j) = 1.
        end do
     end do

! [maskone_r] = 1 means all the rtm grid is land. Used as dummy
! variable so code will not abort with false, non-valid error check

     do j = 1, rtmlat
        do i = 1, numlon_r(j)
           maskone_r(i,j) = 1.
        end do
     end do
     
! --------------------------------------------------------------------
! Map weights from land model grid to rtm grid
! --------------------------------------------------------------------

     write(6,*) 'Initializing land model -> rtm interpolation .....'

! For each rtm grid cell: get lat [jovr_s2r] and lon [iovr_s2r] indices 
! and weights [wovr_s2r] of overlapping atm grid cells 

     call mkmxovr (lsmlon, lsmlat, numlon  , lonw , lats , &
                   rtmlon, rtmlat, numlon_r, lonwh, latsh, &
                   mxovr_s2r     , novr_s2r)

     allocate(iovr_s2r(rtmlon,rtmlat,mxovr_s2r))
     allocate(jovr_s2r(rtmlon,rtmlat,mxovr_s2r))
     allocate(wovr_s2r(rtmlon,rtmlat,mxovr_s2r))
     
     call areaini (lsmlon   , lsmlat  , numlon  , lonw , lats , area  , maskone_s,  &
                   rtmlon   , rtmlat  , numlon_r, lonwh, latsh, area_r, maskone_r,  &
                   mxovr_s2r, novr_s2r, iovr_s2r, jovr_s2r, wovr_s2r)

     write(6,*) 'Successfully made land model -> rtm interpolation'
     write(6,*)

#if (defined COUP_CSM)

! --------------------------------------------------------------------
! Determine which ocean cells might have runoff values. 
! --------------------------------------------------------------------

! First loop over all ocean points and determine which are at the 
! end of rivers by examining if any neighboring points are land and 
! if that land neighbor points into this ocean point. Next loop over all
! ocean points and determine which overlap with at least one land cell.

     ocnrof_num = 0
     ocnrof_mask(:,:) = 0
     do j=1,rtmlat
        do i=1,rtmlon
           if (mask_r(i,j) == 0) then
              if (rdirc(i  ,j-1)==1) ocnrof_mask(i,j) = 1
              if (rdirc(i-1,j-1)==2) ocnrof_mask(i,j) = 1
              if (rdirc(i-1,j  )==3) ocnrof_mask(i,j) = 1
              if (rdirc(i-1,j+1)==4) ocnrof_mask(i,j) = 1
              if (rdirc(i  ,j+1)==5) ocnrof_mask(i,j) = 1
              if (rdirc(i+1,j+1)==6) ocnrof_mask(i,j) = 1
              if (rdirc(i+1,j  )==7) ocnrof_mask(i,j) = 1
              if (rdirc(i+1,j-1)==8) ocnrof_mask(i,j) = 1
              if (ocnrof_mask(i,j) == 0) then
                 do n=1,novr_s2r(i,j)
                    is = iovr_s2r(i,j,n)
                    js = jovr_s2r(i,j,n)
                    if (landmask(is,js)==1 .and. wovr_s2r(i,j,n)>0.) then
                       ocnrof_mask(i,j) = 1
                    end if
	         end do 
              endif
           endif
           if (ocnrof_mask(i,j) == 1) ocnrof_num = ocnrof_num +1
        enddo
     enddo

! allocate ocean runoff vector and indices and determine indices
! need to reset ocnrof_num to 0 and do the counting again because need to first
! first count to allocate vector and must now count to actually determine indices

     allocate(ocnrof_vec  (ocnrof_num)) 
     allocate(ocnrof_iindx(ocnrof_num)) 
     allocate(ocnrof_jindx(ocnrof_num))
     
     ocnrof_num = 0
     do j=1,rtmlat
        do i=1,rtmlon
           if (ocnrof_mask(i,j) == 1) then
              ocnrof_num = ocnrof_num + 1
              ocnrof_iindx(ocnrof_num) = i
              ocnrof_jindx(ocnrof_num) = j
              ocnrof_vec(ocnrof_num) = 0.
           endif
        end do
     enddo
     
#endif

! Deallocate memory for rtm grid  - needed to be done here because
! rtm grid information had to be sent to coupler between calls to 
! Rtmgridini and Rtmlandini

     deallocate(latixy_r)
     deallocate(longxy_r)
     deallocate(latsh)
     deallocate(lonwh) 
     
  endif ! end of if-masterproc block

! Initialize rtm time averaging variables
! Upon restart, the following variables will get new values 
! from the restart file - these values are only valid for
! initial runs

  ncount_rtm    = 0

  ncount_global = 0

  prec_global   = 0.  
  evap_global   = 0.
  runlnd_global = 0.
  runrtm_global = 0.
  volrtm_global = 0.
  ocnrtm_global = 0.

  call get_curr_date(yrold, mon, day, ncsec)

  allocate (totrunin_ave(numpatch)); totrunin_ave(:) = 0.

  allocate (prec_ave(numpatch)); prec_ave(:) = 0.

  allocate (evap_ave(numpatch)); evap_ave(:) = 0.

  allocate (qchan2(numpatch)); qchan2(:) = 0.

  allocate (qchocn2(numpatch)); qchocn2(:) = 0.

  return
end subroutine Rtmlandini

!=======================================================================

subroutine Rtmriverflux ()

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Interface with RTM river routing model
! 
! Method: 
! 
! Author: Sam Levis
! 
!-----------------------------------------------------------------------

  use precision	
  use clm_varder	
  use clm_varpar  , only : lsmlon, lsmlat
  use clm_varmap  , only : begpatch, endpatch, numpatch
  use clm_varsur  , only : numlon, area, landfrac 
  use clm_varctl  , only : rtm_nsteps
  use histFileMod , only : histslf
#if (defined SPMD)
  use spmdMod     , only : masterproc, npes, compute_mpigs_patch, iam
  use mpishorthand, only : mpir8, mpilog, mpicom
#else
  use spmdMod     , only : masterproc
#endif
  use time_manager, only : get_step_size, get_curr_date, get_nstep
  implicit none

! ----------------------- local  variables ---------------------------
! misc variables
  integer  :: io,jo,ir,jr,is,js         !mapping indices
  integer  :: k,n,i,j	                !indices
  real(r8) :: totrunin(numpatch)        !input runoff
  real(r8) :: prec(numpatch)            !input precipitation
  real(r8) :: evap(numpatch)            !input evaporation
#if (defined SPMD)
  real(r8) :: gather1d(numpatch)        !MPI temporary
  real(r8) :: scatter1d(numpatch)       !MPI temporary
#endif
  real(r8) :: wt                        !weight

! inputs to RTM at land model resolution
  real(r8) :: totruninxy(lsmlon,lsmlat) !surface runoff (mm H2O /s)
  real(r8) :: precxy(lsmlon,lsmlat)     !precipitation (mm H2O /s)
  real(r8) :: evapxy(lsmlon,lsmlat)     !evaporation (mm H2O /s)

! outputs returned from RTM converted to land model resolution 
  real(r8) :: flxout_s(lsmlon,lsmlat)   !river flow (m**3)
  real(r8) :: flxocn_s(lsmlon,lsmlat)   !flow into ocean (m**3)

! global balance
  integer  :: yrnew         !year (0, ...)
  integer  :: mon           !month (1, ..., 12) 
  integer  :: day           !day of month (1, ..., 31)
  integer  :: ncsec         !seconds of current date  
  integer  :: ncdate        !current date   

  real(r8) :: prec_sum      !total precipitation (m^3/sec) 
  real(r8) :: evap_sum      !total evaporation (m^3/sec)
  real(r8) :: runlnd_sum    !total input runoff on land grid (m^3/sec)
  real(r8) :: runrtm_sum    !total input runoff on rtm grid (m^3/sec)
  real(r8) :: ocnrtm_sum    !total ocean runoff on rtm grid (m^3/sec)
  real(r8) :: volrtm_sum    !total change in storage on rtm (m^3/sec)

#if (defined SPMD)
  integer :: numsendv(0:npes-1)   !vector of items to be sent
  integer :: numrecvv(0:npes-1)   !vector of items to be received  
  integer :: displsv(0:npes-1)    !displacement vector
  integer :: numsend              !number of items to be sent
  integer :: numrecv              !number of items to be received
  integer :: ier                  !MPI error status
#endif
  integer :: nstep                !time step index
! --------------------------------------------------------------------

! --------------------------------------------------------------------
! RTM inputs 
! --------------------------------------------------------------------

! Make gridded representation of runoff from subgrid patch data
! total surface runoff = surface runoff on soils 
!                      + runoff on glaciers, wetlands, and lakes (P-E) 

!$OMP PARALLEL DO PRIVATE (k)
  do k = begpatch, endpatch
     totrunin(k) = clm(k)%qflx_surf + clm(k)%qflx_qrgwl + clm(k)%qflx_drain
     prec(k)     = clm(k)%forc_rain + clm(k)%forc_snow
     evap(k)     = clm(k)%qflx_evap_tot
  end do

! --------------------------------------------------------------------
! Average fluxes for RTM calculation if appropriate
! --------------------------------------------------------------------

! RTM averaging is not done

  if (rtm_nsteps <= 1) then

#if (defined SPMD)
     call compute_mpigs_patch(1, numsend, numrecvv, displsv)
     call mpi_gatherv (totrunin(begpatch), numsend , mpir8, &
          gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
     if (masterproc) totrunin(:) = gather1d(:)
     call mpi_gatherv (prec(begpatch), numsend , mpir8, &
          gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
     if (masterproc) prec(:) = gather1d(:)
     call mpi_gatherv (evap(begpatch), numsend , mpir8, &
          gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
     if (masterproc) evap(:) = gather1d(:)
#endif           
     if (masterproc) then
        call v2xy (totrunin, 0._r8, totruninxy)
        call v2xy (prec    , 0._r8, precxy)
        call v2xy (evap    , 0._r8, evapxy)
        delt_rtm = get_step_size()
     endif

! RTM averaging is done only done by master processor - however
! all SPMD processe will continue below

  else

     do k = begpatch, endpatch
        totrunin_ave(k) = totrunin_ave(k) + totrunin(k)
        prec_ave(k) = prec_ave(k) + prec(k)
        evap_ave(k) = evap_ave(k) + evap(k)
     end do
     if (masterproc) then
        ncount_rtm = ncount_rtm + 1     
     endif
     nstep = get_nstep()
     if ((mod(nstep,rtm_nsteps)==0) .and. (nstep>1)) then
#if (defined SPMD)
        call compute_mpigs_patch(1, numsend, numrecvv, displsv)
        call mpi_gatherv (totrunin_ave(begpatch), numsend , mpir8, &
             gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
        if (masterproc) totrunin_ave(:) = gather1d(:)
        call mpi_gatherv (prec_ave(begpatch), numsend , mpir8, &
             gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
        if (masterproc) prec_ave(:) = gather1d(:)
        call mpi_gatherv (evap_ave(begpatch), numsend , mpir8, &
             gather1d, numrecvv, displsv, mpir8, 0, mpicom, ier)
        if (masterproc) evap_ave(:) = gather1d(:)
#endif           
        if (masterproc) then
           do k = 1,numpatch
              totrunin_ave(k) = totrunin_ave(k)/ncount_rtm
              prec_ave(k) = prec_ave(k)/ncount_rtm
              evap_ave(k) = evap_ave(k)/ncount_rtm
           end do
           call v2xy(totrunin_ave, 0._r8, totruninxy)
           call v2xy(prec_ave    , 0._r8, precxy)
           call v2xy(evap_ave    , 0._r8, evapxy)
           delt_rtm = ncount_rtm*get_step_size()   !compute delt for rtm
           ncount_rtm = 0                          !reset counter to 0
        endif
        do k = begpatch,endpatch
           totrunin_ave(k) = 0.                    !reset averager
           prec_ave(k) = 0.
           evap_ave(k) = 0.
        end do
     else