rtmmod.f90

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

        call histslf ('QCHOCNR ', qchocn2(begpatch:endpatch))
        call histslf ('QCHANR  ', qchan2(begpatch:endpatch))
        RETURN
     endif

  endif

! --------------------------------------------------------------------
! RTM runoff - only master processor performs runoff computation
! --------------------------------------------------------------------

  if (masterproc) then

! Map from land model grid to RTM grid (intepolate to 1/2 degree resolution)

!$OMP PARALLEL DO PRIVATE (jr,ir,n,is,js,wt)
     do jr = 1,rtmlat
        do ir = 1,numlon_r(jr)
           totrunin_r(ir,jr) = 0.
           do n = 1, novr_s2r(ir,jr)
              if (wovr_s2r(ir,jr,n) > 0.) then
                 is = iovr_s2r(ir,jr,n)
                 js = jovr_s2r(ir,jr,n)
                 wt = wovr_s2r(ir,jr,n)
                 totrunin_r(ir,jr) = totrunin_r(ir,jr) + wt*totruninxy(is,js)*landfrac(is,js)
              end if
           end do
        end do
     end do

! Determine fluxes on 1/2 degree grid

     call Rtm 

#if (defined COUP_CSM)

! Determine ocean runoff vector to send to coupler 
     
     do n = 1,size(ocnrof_vec)
        i = ocnrof_iindx(n)
        j = ocnrof_jindx(n)
        ocnrof_vec(n) = flxocn_r(i,j)/(area_r(i,j)*1000.) ! units of kg/m^2/s
     end do

#endif

! Determine ocean runoff and total runoff on land model grid and
! compute global input runoff on rtm grid, global ocean runoff on
! rtm grid and global change in storage on rtm grid

     do js = 1, lsmlat
        do is =1, numlon(js)
           flxout_s(is,js)  = 0.
           flxocn_s(is,js)  = 0.
        end do
     end do
     do jr = 1,rtmlat
        do ir = 1,numlon_r(jr)
           do n = 1, novr_s2r(ir,jr)
              if (wovr_s2r(ir,jr,n) > 0.) then
                 is = iovr_s2r(ir,jr,n)
                 js = jovr_s2r(ir,jr,n)
                 wt = wovr_s2r(ir,jr,n)
                 flxocn_s(is,js)  = flxocn_s(is,js) + wt*flxocn_r(ir,jr)
                 flxout_s(is,js)  = max(flxout_s(is,js), wt*flxlnd_r(ir,jr))
              end if
           end do
           runrtm(ir,jr) = totrunin_r(ir,jr)*1000.*area_r(ir,jr)
           volrtm(ir,jr) = dvolrdt_r(ir,jr)*1000.*area_r(ir,jr)
        end do
     end do
     
! Determine global quantities and increment global counter

!$OMP PARALLEL DO PRIVATE (js,is)
     do js = 1,lsmlat
        do is = 1,numlon(js)
           precxy(is,js) = precxy(is,js)*area(is,js)*1000.*landfrac(is,js)
           evapxy(is,js) = evapxy(is,js)*area(is,js)*1000.*landfrac(is,js)
           totruninxy(is,js) = totruninxy(is,js)*area(is,js)*1000.*landfrac(is,js)
        end do
     end do

     prec_sum   = sum(precxy)
     evap_sum   = sum(evapxy)
     runlnd_sum = sum(totruninxy)
     runrtm_sum = sum(runrtm)
     volrtm_sum = sum(volrtm)
     ocnrtm_sum = sum(flxocn_r)

     prec_global   = prec_global   + prec_sum
     evap_global   = evap_global   + evap_sum
     runlnd_global = runlnd_global + runlnd_sum
     runrtm_global = runrtm_global + runrtm_sum
     volrtm_global = volrtm_global + volrtm_sum
     ocnrtm_global = ocnrtm_global + ocnrtm_sum

     ncount_global = ncount_global + 1

! Print out diagnostics if appropriate

     write(6,*)
     write(6,F41)'water inst   '
     write(6,F21)'water inst   '
     write(6,F22)'water inst   ',get_nstep(), prec_sum, evap_sum, &
          runlnd_sum, runrtm_sum, volrtm_sum, ocnrtm_sum
     write(6,*)
     
     call get_curr_date(yrnew, mon, day, ncsec)
     ncdate = yrnew*10000 + mon*100 + day
     if (yrnew /= yrold) then
        prec_global   = prec_global/ncount_global
        evap_global   = evap_global/ncount_global
        runlnd_global = runlnd_global/ncount_global
        runrtm_global = runrtm_global/ncount_global
        volrtm_global = volrtm_global/ncount_global
        ocnrtm_global = ocnrtm_global/ncount_global
        ncount_global = 0
        write(6,*)
        write(6,F40)'water tavg   '
        write(6,F21)'water tavg   '
        write(6,F22)'water tavg   ',ncdate, prec_global, evap_global,&
             runlnd_global, runrtm_global, volrtm_global, ocnrtm_global
        write(6,*)
     endif
     yrold = yrnew

! Convert gridded output to vector form 

     call xy2v (lsmlon, lsmlat, flxocn_s, 1, numpatch, qchocn2)
     call xy2v (lsmlon, lsmlat, flxout_s, 1, numpatch, qchan2 )

  endif  !end of if-masterproc block

! Add to history file

#if (defined SPMD)
!
! need to do the following since all other calls to histslf are done from
! begpatch to endpatch in other modules and then an mpi-gather is done in histwrt
!
  call compute_mpigs_patch(1, numrecv, numsendv, displsv)
  if (masterproc) scatter1d(:) = qchocn2(:)
  call mpi_scatterv (scatter1d, numsendv, displsv, mpir8, &
                     qchocn2(begpatch), numrecv , mpir8  , 0, mpicom, ier)
  if (masterproc) scatter1d(:) = qchan2(:)
  call mpi_scatterv (scatter1d, numsendv, displsv, mpir8, &
                     qchan2(begpatch) , numrecv , mpir8  , 0, mpicom, ier)
#endif

  call histslf ('QCHOCNR ', qchocn2(begpatch:endpatch))
  call histslf ('QCHANR  ', qchan2(begpatch:endpatch))
     
  return
end subroutine Rtmriverflux

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

subroutine Rtm 

!----------------------------------------------------------------------- 
! 
! Purpose: 
! River routing model (based on U. Texas code)
! 
! Method: 
! inputs are totrunin_r, flxocn_r, flxlnd_r
! output is dvolrdt_r
! 
! Author: Sam Levis
! 
!-----------------------------------------------------------------------

! ------------------------ local variables ---------------------------
  integer  :: i, j                        !loop indices
  logical  :: lfirst = .true.             !indicates first time through
  real(r8) :: buffern(rtmlon)             !temp buffer
  real(r8) :: buffers(rtmlon)             !temp buffer
  real(r8) :: dvolrdt                     !change in storage (m3/s)
  real(r8) :: sumdvolr(rtmlat)            !global sum (m3/s)
  real(r8) :: sumrunof(rtmlat)            !global sum (m3/s)
  real(r8) :: sumdvolr_tot                !global sum (m3/s)
  real(r8) :: sumrunof_tot                !global sum (m3/s)
  real(r8), parameter :: effvel = 0.35    !effective velocity (m/s)
!------------------------------------------------------------------
             
! At the beginning of a run calculate initialize fluxout 

  if (lfirst) then
     do j=1,rtmlat
        do i=1,rtmlon
           flxocn_r(i,j) = 0.  !runoff on ocean cell
           flxlnd_r(i,j) = 0.  !runoff on land  cell
           fluxout(i,j) = 0.
           if (mask_r(i,j)==1) then
              fluxout(i,j) = volr(i,j) * effvel/ddist(i,j)
              fluxout(i,j) = min(fluxout(i,j), volr(i,j) / delt_rtm)
           endif
        enddo
     enddo
     lfirst = .false.
  end if

! Determine fluxout at extended points and at southern and northern outer lats

  fluxout(rtmlon+1,rtmlat+1) = fluxout(1,1)
  fluxout(rtmlon+1,0)        = fluxout(1,rtmlat)
  fluxout(0,0)               = fluxout(rtmlon,rtmlat)
  fluxout(0,rtmlat+1)        = fluxout(rtmlon,1)

  do i=1,rtmlon
     fluxout(i,0)        = fluxout(i,rtmlat)
     fluxout(i,rtmlat+1) = fluxout(i,1)
     buffern(i)          = fluxout(i,rtmlat)
     buffers(i)          = fluxout(i,1)
  enddo
  do j=1,rtmlat
     fluxout(0,j)        = fluxout(rtmlon,j)
     fluxout(rtmlon+1,j) = fluxout(1,j)
  enddo
  do i=0,rtmlon+1
     fluxout(i,0)        = buffern(mod(i+rtmlon/2-1,rtmlon)+1)
     fluxout(i,rtmlat+1) = buffers(mod(i+rtmlon/2-1,rtmlon)+1)
  enddo

! Determine cell-to-cell transport - calculate sfluxin

!$OMP PARALLEL DO PRIVATE (I,J)
  do j=1,rtmlat
     do i=1,rtmlon
        sfluxin(i,j) = 0.
        if (rdirc(i  ,j-1)==1) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i  ,j-1)
        if (rdirc(i-1,j-1)==2) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i-1,j-1)
        if (rdirc(i-1,j  )==3) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i-1,j  )
        if (rdirc(i-1,j+1)==4) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i-1,j+1)
        if (rdirc(i  ,j+1)==5) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i  ,j+1)
        if (rdirc(i+1,j+1)==6) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i+1,j+1)
        if (rdirc(i+1,j  )==7) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i+1,j  )
        if (rdirc(i+1,j-1)==8) sfluxin(i,j) = sfluxin(i  ,j) + fluxout(i+1,j-1)
     enddo
  enddo

! Loops above and below must remain separate because fluxout is updated below

  sumdvolr(:) = 0.
  sumrunof(:) = 0.
!$OMP PARALLEL DO PRIVATE (I,J,DVOLRDT)
  do j=1,rtmlat
     do i=1,rtmlon

! calculate change in cell storage volume change units for totrunin from kg/m2s==mm/s -> m3/s
        
        dvolrdt = sfluxin(i,j) - fluxout(i,j) + 0.001*totrunin_r(i,j)*rivarea(i,j)

! calculate flux out of a cell:
! land: do not permit change in cell storage volume greater than volume present
! make up for the difference with an extraction term (eg from aquifers)
! ocean: do not permit negative change in cell storage volume,
! because at ocean points cell storage volume equals zero
! water balance check (in mm/s), convert runinxy from mm/s to m/s (* 1.e-3) 
! and land model area from km**2 to m**2 (* 1.e6)

        if (mask_r(i,j) == 1) then         ! land points
           volr(i,j)    = volr(i,j) + dvolrdt*delt_rtm
           fluxout(i,j) = volr(i,j) * effvel/ddist(i,j)
           fluxout(i,j) = min(fluxout(i,j), volr(i,j) / delt_rtm)
           flxlnd_r(i,j) = fluxout(i,j)
           flxocn_r(i,j) = 0.
        else                               ! ocean points
           flxlnd_r(i,j) = 0.
           flxocn_r(i,j) = dvolrdt
        endif
        sumdvolr(j) = sumdvolr(j) + dvolrdt
        sumrunof(j) = sumrunof(j) + totrunin_r(i,j)*1000.*area_r(i,j)
        dvolrdt_r(i,j) = 1000.*dvolrdt/rivarea(i,j)

     enddo
  enddo

! Global water balance calculation and error check

  sumdvolr_tot = 0.
  sumrunof_tot = 0.
  do j = 1,rtmlat
     sumdvolr_tot = sumdvolr_tot + sumdvolr(j)
     sumrunof_tot = sumrunof_tot + sumrunof(j)
  end do
  if (abs((sumdvolr_tot-sumrunof_tot)/sumrunof_tot) > 0.01) then
     write(6,*) 'RTM Error: sumdvolr= ',sumdvolr_tot,&
          ' not equal to sumrunof= ',sumrunof_tot
     call endrun
  end if
  
  return
end subroutine Rtm

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

#endif

end module RtmMod