zm_conv.f90

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

!--------------------------Data statements------------------------------
!
   logical momentm
   data momentm/.FALSE./
!
! Set internal variable "msg" (convection limit) to "limcnv-1"
!
   msg = limcnv - 1
!
! initialize necessary arrays.
! zero out variables not used in cam
!
   qtg(:,:) = 0.
   ttg(:,:) = 0.
   mcon(:,:) = 0.
   do i = 1,ncol
      paprc(i) = 0.
      paprs(i) = 0.
      pcpc(i) = 0
   end do
   psdiss = 0.
   psheat = 0.
   psevap = 0.
   psrain = 0.
!      jlatpr = 32
!
! initialize convective tendencies
!
   do k = 1,pver
      do i = 1,ncol
         dqdt(i,k) = 0.
         dsdt(i,k) = 0.
         dudt(i,k) = 0.
         dvdt(i,k) = 0.
!            deltaq(i,k) = qh(i,k,1)
!            deltat(i,k) = t(i,k)
         pcpck(i,k) = 0.
         pflx(i,k) = 0.
         pflxg(i,k) = 0.
         cme(i,k) = 0.
!++bee
         cmfdqr(i,k) = 0.
!--bee
         zdu(i,k) = 0.
         ql(i,k) = 0.
         qlg(i,k) = 0.
      end do
   end do
   do i = 1,ncol
      pflx(i,pverp) = 0
      pflxg(i,pverp) = 0
   end do
   if (.not.momentm) then
      do k = 1,pver
         do i = 1,ncol
            u(i,k) = 0.
            v(i,k) = 0.
         end do
      end do
   end if
!
   do i = 1,ncol
      pblt(i) = pver
      pcpr(i) = 0.
      pcps(i) = 0.
      dsubcld(i) = 0.
      sumq(i) = 0.
!         sumdt(i) = 0.
!         sumdq(i) = 0.
      pcpdh(i) = rgrav
      jctop(i) = pver
      jcbot(i) = 1
   end do
!
! calculate local pressure (mbs) and height (m) for both interface
! and mid-layer locations.
!
   do i = 1,ncol
      zs(i) = geos(i)*rgrav
      pf(i,pver+1) = paph(i,pver+1)*0.01
      zf(i,pver+1) = zi(i,pver+1) + zs(i)
   end do
   do k = 1,pver
      do i = 1,ncol
         p(i,k) = pap(i,k)*0.01
         pf(i,k) = paph(i,k)*0.01
         z(i,k) = zm(i,k) + zs(i)
         zf(i,k) = zi(i,k) + zs(i)
      end do
   end do
!
   do k = pver - 1,msg + 1,-1
      do i = 1,ncol
         if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5) pblt(i) = k
      end do
   end do
!
! store incoming specific humidity field for subsequent calculation
! of precipitation (through change in storage).
! convert from specific humidity (bounded by qmin) to mixing ratio.
! define dry static energy (normalized by cp).
!
   do k = 1,pver
      do i = 1,ncol
         qeff = max(qh(i,k,1),qmin)
         q(i,k) = qeff
         s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
         tp(i,k)=0.0
         shat(i,k) = s(i,k)
         qhat(i,k) = q(i,k)
         dp(i,k) = dpp(i,k)*0.01
         qg(i,k) = q(i,k)
         tg(i,k) = t(i,k)
         pg(i,k) = p(i,k)
         zg(i,k) = z(i,k)
         sg(i,k) = s(i,k)
         tpg(i,k) = tp(i,k)
         zfg(i,k) = zf(i,k)
         qstpg(i,k) = q(i,k)
         ug(i,k) = u(i,k)
         vg(i,k) = v(i,k)
         dlg(i,k) = 0
         dlf(i,k) = 0
      end do
   end do
   do i = 1,ncol
      zfg(i,pver+1) = zf(i,pver+1)
      capeg(i) = 0.
      lclg(i) = 1
      lelg(i) = pver
      maxg(i) = 1
      tlg(i) = 400.
      dsubcld(i) = 0.
   end do
!
! evaluate covective available potential energy (cape).
!
   call buoyan(lchnk   ,ncol    , &
               q       ,t       ,p       ,z       ,pf       , &
               tp      ,qstp    ,tl      ,rl      ,cape     , &
               pblt    ,lcl     ,lel     ,lon     ,maxi     , &
               rgas    ,grav    ,cpres   ,msg     ,nstep    , &
               tpert   )
!
! determine whether grid points will undergo some deep convection
! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
! (require cape.gt. 0 and lel<lcl as minimum conditions).
!
   lengath = 0
   do i=1,ncol
      if (cape(i) > capelmt) then
         lengath = lengath + 1
         index(lengath) = i
      end if
   end do

   if (lengath.eq.0) return
!CDIR$ IVDEP
   do ii=1,lengath
      i=index(ii)
      ideep(ii)=i
   end do
!c        jyes = 0
!c        jno = nlon - 1 + 2
!c        do il = 1,nlon
!c           if (cape(il).gt.capelmt) then
!c              jyes = jyes + 1
!c              ideep(jyes) = il
!c           else
!c              jno = jno - 1
!c              ideep(jno) = il
!c           end if
!c        end do
!c        lengath = jyes
!c        if (lengath.eq.0) return
!
! obtain gathered arrays necessary for ensuing calculations.
!
   do k = 1,pver
      do i = 1,lengath
         dp(i,k) = 0.01*dpp(ideep(i),k)
         qg(i,k) = q(ideep(i),k)
         tg(i,k) = t(ideep(i),k)
         pg(i,k) = p(ideep(i),k)
         zg(i,k) = z(ideep(i),k)
         sg(i,k) = s(ideep(i),k)
         tpg(i,k) = tp(ideep(i),k)
         zfg(i,k) = zf(ideep(i),k)
         qstpg(i,k) = qstp(ideep(i),k)
         ug(i,k) = u(ideep(i),k)
         vg(i,k) = v(ideep(i),k)
      end do
   end do
!
   do i = 1,lengath
      zfg(i,pver+1) = zf(ideep(i),pver+1)
   end do
   do i = 1,lengath
      capeg(i) = cape(ideep(i))
      lclg(i) = lcl(ideep(i))
      lelg(i) = lel(ideep(i))
      maxg(i) = maxi(ideep(i))
      tlg(i) = tl(ideep(i))
   end do
!
! calculate sub-cloud layer pressure "thickness" for use in
! closure and tendency routines.
!
   do k = msg + 1,pver
      do i = 1,lengath
         if (k >= maxg(i)) then
            dsubcld(i) = dsubcld(i) + dp(i,k)
         end if
      end do
   end do
!
! define array of factors (alpha) which defines interfacial
! values, as well as interfacial values for (q,s) used in
! subsequent routines.
!
   do k = msg + 2,pver
      do i = 1,lengath
!            alpha(i,k) = 0.5
         sdifr = 0.
         qdifr = 0.
         if (sg(i,k) > 0. .or. sg(i,k-1) > 0.) &
            sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
         if (qg(i,k) > 0. .or. qg(i,k-1) > 0.) &
            qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
         if (sdifr > 1.E-6) then
            shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k))
         else
            shat(i,k) = 0.5* (sg(i,k)+sg(i,k-1))
         end if
         if (qdifr > 1.E-6) then
            qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k))
         else
            qhat(i,k) = 0.5* (qg(i,k)+qg(i,k-1))
         end if
      end do
   end do
!
! obtain cloud properties.
!
   call cldprp(lchnk   , &
               qg      ,tg      ,ug      ,vg      ,pg      , &
               zg      ,sg      ,mu      ,eu      ,du      , &
               md      ,ed      ,sd      ,qd      ,mc      , &
               qu      ,su      ,zfg     ,qs      ,hmn     , &
               hsat    ,shat    ,qlg     ,totpcp  ,totevp  , &
               cmeg    ,maxg    ,lelg    ,jt      ,jlcl    , &
               maxg    ,j0      ,jd      ,rl      ,lengath , &
               rgas    ,grav    ,cpres   ,msg     ,nstep   , &
                        pflxg   ,evpg    ,cug     ,mu2     , &
               eu2     ,du2     ,md2     ,ed2     ,cmfdqrg , &
               limcnv  )
!
! convert detrainment from units of "1/m" to "1/mb".
!
   do k = msg + 1,pver
      do i = 1,lengath
         du(i,k) = du(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         eu(i,k) = eu(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         ed(i,k) = ed(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         cug(i,k) = cug(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
!++bee, bug fix from /home/pjr/ccm2/omega0.10.1/fspj01/.
         cmeg(i,k) = cmeg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         cmfdqrg(i,k) = cmfdqrg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
!--bee
         evpg(i,k) = evpg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         du2(i,k) = du2(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         eu2(i,k) = eu2(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
         ed2(i,k) = ed2(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
      end do
   end do

   call closure(lchnk   , &
                qg      ,tg      ,pg      ,zg      ,sg      , &
                tpg     ,qs      ,qu      ,su      ,mc      , &
                du      ,mu      ,md      ,qd      ,sd      , &
                qhat    ,shat    ,dp      ,qstpg   ,zfg     , &
                qlg     ,dsubcld ,mb      ,capeg   ,tlg     , &
                lclg    ,lelg    ,jt      ,maxg    ,1       , &
                lengath ,rgas    ,grav    ,cpres   ,rl      , &
                msg     ,capelmt ,nstep   )
!
! limit cloud base mass flux to theoretical upper bound.
!
   do i=1,lengath
      mumax(i) = 0
   end do
   do k=msg + 2,pver
      do i=1,lengath
         mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
      end do
   end do
   do i=1,lengath
      if (mumax(i) > 0.) then
         mb(i) = min(mb(i),0.5/(delt*mumax(i)))
      else
         mb(i) = 0.
      endif
   end do
   do k=msg+1,pver
      do i=1,lengath
         mu(i,k) = mu(i,k)*mb(i)
         md(i,k) = md(i,k)*mb(i)
         mc(i,k) = mc(i,k)*mb(i)
         du(i,k) = du(i,k)*mb(i)
         eu(i,k) = eu(i,k)*mb(i)
         ed(i,k) = ed(i,k)*mb(i)
         cmeg(i,k) = cmeg(i,k)*mb(i)
!++bee
         cmfdqrg(i,k) = cmfdqrg(i,k)*mb(i)
!--bee
         cug(i,k) = cug(i,k)*mb(i)
         evpg(i,k) = evpg(i,k)*mb(i)
         pflxg(i,k+1) = pflxg(i,k+1)*mb(i)*100./grav
         mu2(i,k) = mu2(i,k)*mb(i)
         md2(i,k) = md2(i,k)*mb(i)
         du2(i,k) = du2(i,k)*mb(i)
         eu2(i,k) = eu2(i,k)*mb(i)
         ed2(i,k) = ed2(i,k)*mb(i)
      end do
   end do
   do i = 1,lengath
!
! totpcp from cldprp has the dimension of kg/kg, here it is
! converted to kg/(m^2*s), the precipitation rate
!
      totpcp(i) = totpcp(i)*mb(i)*100./grav
      totevp(i) = totevp(i)*mb(i)*100./grav
   end do
!
! compute temperature and moisture changes due to convection.
!
   call q1q2_pjr(lchnk   , &
                 dqdt    ,dsdt    ,qg      ,qs      ,qu      , &
                 su      ,du      ,qhat    ,shat    ,dp      , &
                 mu      ,md      ,sd      ,qd      ,qlg     , &
                 dsubcld ,jt      ,maxg    ,1       ,lengath , &
                 cpres   ,rl      ,msg     ,nstep   ,          &
                 dlg     ,evpg    ,cug     )
!
! compute momentum changes due to convection, if desired (i.e
! if logical switch set).
!
!
!      if(momentm)
!       then
!        call moment(dudt,dvdt,du,alpha,dp,ed,eu,mc,md,mu,
!     1             pg,qd,qu,qhat,sd,su,shat,ud,vd,tg,ug,vg,zg,zfg,
!     2             dsubcld,maxg,jd,jt,rl,
!     3             msg,2.*delt,grav,cpres,rgas,pver,1,lengath,nlond,lat(1))
!      endif
!
! move the convective transport outside of conv_cam.
!
! gather back temperature and mixing ratio.
!
   do k = msg + 1,pver
      do i = 1,lengath
         psdiss = psdiss + (dudt(i,k)*u(ideep(i),k)+dvdt(i,k)*v(ideep(i),k))* &
                  dpp(ideep(i),k)/grav
!
! q is updated to compute net precip, and then reset to old value.
! the last line is overwritten. so the input basic variables, i.e.
! q, t, u and v are updated to include convective increments.
! (5/2/95)
!
         q(ideep(i),k) = q(ideep(i),k) + 2.*delt*dqdt(i,k)
!**BAB don't update t (real state variable)
!!$         t(ideep(i),k) = t(ideep(i),k) + 2.*delt*dsdt(i,k)
         u(ideep(i),k) = u(ideep(i),k) + 2.*delt*dudt(i,k)
         v(ideep(i),k) = v(ideep(i),k) + 2.*delt*dvdt(i,k)