moistconvection.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

            ds1(i) = etagdt(i)*(sbh(i,k+1) - sc(i))*rpdel(i,k+1)
            dq1(i) = etagdt(i)*(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
            ds2(i) = (etagdt(i)*(sc(i) - sbh(i,k+1)) +  &
                     hlat*grav*cldwtr(i) - beta(i)*etagdt(i)*(sc(i) - sbh(i,k)))*rpdel(i,k)
            dq2(i) = (etagdt(i)*(shc(i) - shbh(i,k+1)) - grav*rnwtr(i) - beta(i)* &
                     etagdt(i)*(shc(i) - shbh(i,k)))*rpdel(i,k)
            ds3(i) = beta(i)*(etagdt(i)*(sc(i) - sbh(i,k)) - hlat*grav*cldwtr(i))* &
                     rpdel(i,k-1)
            dq3(i) = beta(i)*etagdt(i)*(shc(i) - shbh(i,k))*rpdel(i,k-1)
!
! Isolate convective fluxes for later diagnostics
!
            fslkp = eta(i)*(sc(i) - sbh(i,k+1))
            fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - hlat*cldwtr(i)*rtdt)
            fqlkp = eta(i)*(shc(i) - shbh(i,k+1))
            fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
!
! Update thermodynamic profile (update sb, hb, & hbs later)
!
            tb (i,k+1) = tb(i,k+1)  + ds1(i)*rcp
            tb (i,k  ) = tb(i,k  )  + ds2(i)*rcp
            tb (i,k-1) = tb(i,k-1)  + ds3(i)*rcp
            shb(i,k+1) = shb(i,k+1) + dq1(i)
            shb(i,k  ) = shb(i,k  ) + dq2(i)
            shb(i,k-1) = shb(i,k-1) + dq3(i)
!
! ** Update diagnostic information for final budget **
! Tracking precipitation, temperature & specific humidity tendencies,
! rainout term, convective mass flux, convective liquid
! water static energy flux, and convective total water flux
! The variable afac makes the necessary adjustment to the
! diagnostic fluxes to account for adjustment time scale based on
! how relaxation time scale is to be applied (column vs. triplet)
!
            prec(i)    = prec(i) + (rnwtr(i)/rhoh2o)*adjfac
!
! The following variables have units of "units"/second
!
            cmfdt (i,k+1) = cmfdt (i,k+1) + ds1(i)*rtdt*adjfac
            cmfdt (i,k  ) = cmfdt (i,k  ) + ds2(i)*rtdt*adjfac
            cmfdt (i,k-1) = cmfdt (i,k-1) + ds3(i)*rtdt*adjfac
            cmfdq (i,k+1) = cmfdq (i,k+1) + dq1(i)*rtdt*adjfac
            cmfdq (i,k  ) = cmfdq (i,k  ) + dq2(i)*rtdt*adjfac
            cmfdq (i,k-1) = cmfdq (i,k-1) + dq3(i)*rtdt*adjfac
            qc    (i,k  ) = (grav*rnwtr(i)*rpdel(i,k))*rtdt*adjfac
            cmfdqr(i,k  ) = cmfdqr(i,k  ) + qc(i,k)
            cmfmc (i,k+1) = cmfmc (i,k+1) + eta(i)*adjfac
            cmfmc (i,k  ) = cmfmc (i,k  ) + beta(i)*eta(i)*adjfac
!
! The following variables have units of w/m**2
!
            cmfsl (i,k+1) = cmfsl (i,k+1) + fslkp*adjfac
            cmfsl (i,k  ) = cmfsl (i,k  ) + fslkm*adjfac
            cmflq (i,k+1) = cmflq (i,k+1) + hlat*fqlkp*adjfac
            cmflq (i,k  ) = cmflq (i,k  ) + hlat*fqlkm*adjfac
30          continue
!
! Next, convectively modify passive constituents
! For now, when applying relaxation time scale to thermal fields after
! entire column has undergone convective overturning, constituents will
! be mixed using a "relaxed" value of the mass flux determined above
! Although this will be inconsistant with the treatment of the thermal
! fields, it's computationally much cheaper, no more-or-less justifiable,
! and consistent with how the history tape mass fluxes would be used in
! an off-line mode (i.e., using an off-line transport model)
!
            do 50 m=2,pcnst+pnats    ! note: indexing assumes water is first field
!CDIR$ IVDEP
               do 40 ii=1,len1
                  i = indx1(ii)
!
! If any of the reported values of the constituent is negative in
! the three adjacent levels, nothing will be done to the profile
!
                  if ((dq(i,k+1,m) < 0.0) .or. (dq(i,k,m) < 0.0) .or. (dq(i,k-1,m) < 0.0)) go to 40
!
! Specify constituent interface values (linear interpolation)
!
                  cmrh(i,k  ) = 0.5*(dq(i,k-1,m) + dq(i,k  ,m))
                  cmrh(i,k+1) = 0.5*(dq(i,k  ,m) + dq(i,k+1,m))
!
! Specify perturbation properties of constituents in PBL
!
                  pblhgt = max(pblht(i),1.0_r8)
                  if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0 ) then
                     fac1 = max(0.0_r8,1.0-zm(i,k+1)/pblhgt)
                     cmrc(i) = dq(i,k+1,m) + qpert(i,m)*fac1
                  else
                     cmrc(i) = dq(i,k+1,m)
                  end if
!
! Determine fluxes, flux divergence => changes due to convection
! Logic must be included to avoid producing negative values. A bit
! messy since there are no a priori assumptions about profiles.
! Tendency is modified (reduced) when pending disaster detected.
!
                  botflx   = etagdt(i)*(cmrc(i) - cmrh(i,k+1))*adjfac
                  topflx   = beta(i)*etagdt(i)*(cmrc(i)-cmrh(i,k))*adjfac
                  dcmr1(i) = -botflx*rpdel(i,k+1)
                  efac1    = 1.0
                  efac2    = 1.0
                  efac3    = 1.0
!
                  if (dq(i,k+1,m)+dcmr1(i) < 0.0) then
                     efac1 = max(tiny,abs(dq(i,k+1,m)/dcmr1(i)) - eps)
                  end if
!
                  if (efac1 == tiny .or. efac1 > 1.0) efac1 = 0.0
                  dcmr1(i) = -efac1*botflx*rpdel(i,k+1)
                  dcmr2(i) = (efac1*botflx - topflx)*rpdel(i,k)
!
                  if (dq(i,k,m)+dcmr2(i) < 0.0) then
                     efac2 = max(tiny,abs(dq(i,k  ,m)/dcmr2(i)) - eps)
                  end if
!
                  if (efac2 == tiny .or. efac2 > 1.0) efac2 = 0.0
                  dcmr2(i) = (efac1*botflx - efac2*topflx)*rpdel(i,k)
                  dcmr3(i) = efac2*topflx*rpdel(i,k-1)
!
                  if (dq(i,k-1,m)+dcmr3(i) < 0.0) then
                     efac3 = max(tiny,abs(dq(i,k-1,m)/dcmr3(i)) - eps)
                  end if
!
                  if (efac3 == tiny .or. efac3 > 1.0) efac3 = 0.0
                  efac3    = min(efac2,efac3)
                  dcmr2(i) = (efac1*botflx - efac3*topflx)*rpdel(i,k)
                  dcmr3(i) = efac3*topflx*rpdel(i,k-1)
!
                  dq(i,k+1,m) = dq(i,k+1,m) + dcmr1(i)
                  dq(i,k  ,m) = dq(i,k  ,m) + dcmr2(i)
                  dq(i,k-1,m) = dq(i,k-1,m) + dcmr3(i)
40                continue
50                continue                ! end of m=2,pcnst+pnats loop
!
! Constituent modifications complete
!
                  if (k == limcnv+1) go to 60
!
! Complete update of thermodynamic structure at integer levels
! gather/scatter points that need new values of shbs and gamma
!
                  do ii=1,len1
                     i = indx1(ii)
                     vtemp1(ii     ) = tb(i,k)
                     vtemp1(ii+len1) = tb(i,k-1)
                     vtemp2(ii     ) = pmid(i,k)
                     vtemp2(ii+len1) = pmid(i,k-1)
                  end do
                  call vqsatd (vtemp1  ,vtemp2  ,estemp  ,vtemp3  , vtemp4  , &
                               2*len1   )    ! using estemp as extra long vector
!CDIR$ IVDEP
                  do ii=1,len1
                     i = indx1(ii)
                     shbs(i,k  ) = vtemp3(ii     )
                     shbs(i,k-1) = vtemp3(ii+len1)
                     gam(i,k  ) = vtemp4(ii     )
                     gam(i,k-1) = vtemp4(ii+len1)
                     sb (i,k  ) = sb(i,k  ) + ds2(i)
                     sb (i,k-1) = sb(i,k-1) + ds3(i)
                     hb (i,k  ) = sb(i,k  ) + hlat*shb(i,k  )
                     hb (i,k-1) = sb(i,k-1) + hlat*shb(i,k-1)
                     hbs(i,k  ) = sb(i,k  ) + hlat*shbs(i,k  )
                     hbs(i,k-1) = sb(i,k-1) + hlat*shbs(i,k-1)
                  end do
!
! Update thermodynamic information at half (i.e., interface) levels
!
!CDIR$ IVDEP
                  do ii=1,len1
                     i = indx1(ii)
                     sbh (i,k) = 0.5*(sb(i,k) + sb(i,k-1))
                     shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
                     hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
                     sbh (i,k-1) = 0.5*(sb(i,k-1) + sb(i,k-2))
                     shbh(i,k-1) = qhalf(shb(i,k-2),shb(i,k-1),shbs(i,k-2),shbs(i,k-1))
                     hbh (i,k-1) = sbh(i,k-1) + hlat*shbh(i,k-1)
                  end do
!
#if ( defined DIAGNS )
!#######################################################################
!#                                                                     #
!#    this update necessary, only if debugging diagnostics requested   #
!#                                                                     #
                  do i=1,ncol
                     if (lat(i) == jloc .and. nstep >= nsloc) then
                        call vqsatd(tb(i,k+1),pmid(i,k+1),es,shbs(i,k+1),gam(i,k+1), &
                                    1)
                        sb (i,k+1) = sb(i,k+1) + ds1(i)
                        hb (i,k+1) = sb(i,k+1) + hlat*shb(i,k+1)
                        hbs(i,k+1) = sb(i,k+1) + hlat*shbs(i,k+1)
                        kpp = k + 2
                        if (k+1 == pver) kpp = k + 1
                        do kp=k+1,kpp
                           kpm1 = kp-1
                           sbh(i,kp)  = 0.5*(sb(i,kpm1) + sb(i,kp))
                           shbh(i,kp) = qhalf(shb(i,kpm1),shb(i,kp),shbs(i,kpm1),shbs(i,kp))
                           hbh(i,kp)  = sbh(i,kp) + hlat*shbh(i,kp)
                        end do
                     end if
                  end do
!#                                                                     #
!#          diagnostic output                                          #
!#                                                                     #
                  do i=1,ncol
                    if ((lat(i)== jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
                     write(6, 8060) k
                     fac = grav*864.
                     do n=limcnv,pver
                        rh  = shb(i,n)/shbs(i,n)
                        write(6,8020)shbh(i,n),sbh(i,n),hbh(i,n),fac*cmfmc(i,n), &
                                     cmfsl(i,n), cmflq(i,n)
!--------------write(6, 8050)
!--------------write(6, 8030) fac*cmfmc(i,n),cmfsl(i,n), cmflq(i,n)
                        write(6, 8040) tb(i,n),shb(i,n),rh,sb(i,n),hb(i,n), &
                                       hbs(i,n), ztodt*cmfdt(i,n),ztodt*cmfdq(i,n), &
	                               ztodt*cmfdqr(i,n)
                     end do
                     write(6, 8000) prec(i)
                    end if
                  end do
!#                                                                     #
!#                                                                     #
!#######################################################################
#endif
!
! Ensure that dzcld is reset if convective mass flux zero
! specify the current vertical extent of the convective activity
! top of convective layer determined by size of overshoot param.
!
60                do i=1,ncol
                     etagt0 = eta(i).gt.0.0
                     if ( .not. etagt0) dzcld(i) = 0.0
                     if (etagt0 .and. beta(i) > betamn) then
                        ktp = k-1
                     else
                        ktp = k
                     end if
                     if (etagt0) then
                        rk=k
                        rktp=ktp
                        cnt(i) = min(cnt(i),rktp)
                        cnb(i) = max(cnb(i),rk)
                     end if
                  end do
70                continue                  ! end of k loop
!
! ** apply final thermodynamic tendencies **
!
!**BAB don't update input profiles
!!$                  do k=limcnv,pver
!!$                     do i=1,ncol
!!$                        t (i,k) = t (i,k) + cmfdt(i,k)*ztodt
!!$                        q(i,k,1) = q(i,k,1) + cmfdq(i,k)*ztodt
!!$                     end do
!!$                  end do
! Set output q tendencies 
      dq(:ncol,:,1 ) = cmfdq(:ncol,:)
      dq(:ncol,:,2:) = (dq(:ncol,:,2:) - q(:ncol,:,2:))/ztodt
!
! Kludge to prevent cnb-cnt from being zero (in the event
! someone decides that they want to divide by this quantity)
!
                  do i=1,ncol
                     if (cnb(i) /= 0.0 .and. cnb(i) == cnt(i)) then
                        cnt(i) = cnt(i) - 1.0
                     end if
                  end do
!
                  do i=1,ncol
                     precc(i) = prec(i)*rtdt
                  end do
!
#if ( defined DIAGNS )
!#######################################################################
!#                                                                     #
!#    we're done ... show final result if debug diagnostics requested  #
!#                                                                     #
                  do i=1,ncol
                    if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
                     fac = grav*864.
                     write(6, 8010)
                     do k=limcnv,pver
                        rh = shb(i,k)/shbs(i,k)
                        write(6, 8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k), &
                                       cmfsl(i,k), cmflq(i,k)
                        write(6, 8040) tb(i,k),shb(i,k),rh   ,sb(i,k),hb(i,k), &
                                       hbs(i,k), ztodt*cmfdt(i,k),ztodt*cmfdq(i,k), &
                                       ztodt*cmfdqr(i,k)
                     end do
                     write(6, 8000) prec(i)
!#                                                                     #
!#       approximate vertical integral of moist static energy and      #
!#       total preciptable water after adjustment and output changes   #
!#                                                                     #
                     hsum2 = 0.0
                     qsum2 = 0.0
                     do k=limcnv,pver
                        hsum2 = hsum2 + pdel(i,k)*rgrav*hb(i,k)
                        qsum2 = qsum2 + pdel(i,k)*rgrav*shb(i,k)
                     end do
!#                                                                     #
                     write (6,8070) hsum1, hsum2, abs(hsum2-hsum1)/hsum2, &
                                    qsum1, qsum2, abs(qsum2-qsum1)/qsum2
                    end if
                  enddo
!#                                                                     #
!#                                                                     #
!#######################################################################
#endif
                  return                 ! we're all done ... return to calling procedure
#if ( defined DIAGNS )
!
! Formats
!
8000              format(///,10x,'PREC = ',3pf12.6,/)
8010              format('1**        TB      SHB      RH       SB', &
                        '       HB      HBS      CAH      CAM       PRECC ', &
                        '     ETA      FSL       FLQ     **', /)
8020              format(' ----- ',     9x,3p,f7.3,2x,2p,     9x,-3p,f7.3,2x, &
                        f7.3, 37x, 0p,2x,f8.2,  0p,2x,f8.2,2x,f8.2, ' ----- ')
8030              format(' ----- ',  0p,82x,f8.2,  0p,2x,f8.2,2x,f8.2, &
                         ' ----- ')
8040              format(' - - - ',f7.3,2x,3p,f7.3,2x,2p,f7.3,2x,-3p,f7.3,2x, &
                        f7.3, 2x,f8.3,2x,0p,f7.3,3p,2x,f7.3,2x,f7.3,30x, &
                         ' - - - ')
8050              format(' ----- ',110x,' ----- ')
8060              format('1 K =>',  i4,/, &
                           '           TB      SHB      RH       SB', &
                           '       HB      HBS      CAH      CAM       PREC ', &
                           '     ETA      FSL       FLQ', /)
8070              format(' VERTICALLY INTEGRATED MOIST STATIC ENERGY BEFORE, AFTER', &
                        ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/, &
                        ' VERTICALLY INTEGRATED MOISTURE            BEFORE, AFTER' &,
                        ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/)
8080              format(' BETA, ETA => ', 1p,2e12.3)
8090              format (' k+1, sc, shc, hc => ', 1x, i2, 1p, 3e12.4)
#endif
!
end subroutine cmfmca
end module moistconvection