LCOV - code coverage report
Current view: | top level - physics/cam - clubb_mf.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 20 | 283 | 7.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 8 | 12.5 %

              Line data    Source code
           1             : module clubb_mf
           2             :
           3             : ! =============================================================================== !
           4             : ! Mass-flux module for use with CLUBB                                             !
           5             : ! Together (CLUBB+MF) they comprise a eddy-diffusivity mass-flux approach (EDMF)  !
           6             : ! =============================================================================== !
           7             :
           8             :   use shr_kind_mod,  only: r8=>shr_kind_r8
           9             :   use spmd_utils,    only: masterproc
          10             :   use cam_logfile,   only: iulog
          11             :
          12             :   implicit none
          13             :   private
          14             :   save
          15             :
          16             :   public :: integrate_mf, &
          17             :             clubb_mf_readnl, &
          18             :             do_clubb_mf, &
          19             :             do_clubb_mf_diag
          20             :
          21             :   real(r8) :: clubb_mf_L0   = 0._r8
          22             :   real(r8) :: clubb_mf_ent0 = 0._r8
          23             :   integer  :: clubb_mf_nup  = 0
          24             :   logical, protected :: do_clubb_mf = .false.
          25             :   logical, protected :: do_clubb_mf_diag = .false.
          26             :
          27             :   contains
          28             :
          29        1536 :   subroutine clubb_mf_readnl(nlfile)
          30             :
          31             :   ! =============================================================================== !
          32             :   ! MF namelists                                                                    !
          33             :   ! =============================================================================== !
          34             :
          35             :     use namelist_utils,  only: find_group_name
          36             :     use cam_abortutils,  only: endrun
          37             :     use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_integer, mpi_logical
          38             :
          39             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          40             :
          41             :     character(len=*), parameter :: sub = 'clubb_mf_readnl'
          42             :
          43             :     integer :: iunit, read_status, ierr
          44             :
          45             :
          46             :     namelist /clubb_mf_nl/ clubb_mf_L0, clubb_mf_ent0, clubb_mf_nup, do_clubb_mf, do_clubb_mf_diag
          47             :
          48        1536 :     if (masterproc) then
          49           2 :       open( newunit=iunit, file=trim(nlfile), status='old' )
          50           2 :       call find_group_name(iunit, 'clubb_mf_nl', status=read_status)
          51           2 :       if (read_status == 0) then
          52           2 :          read(iunit, clubb_mf_nl, iostat=ierr)
          53           2 :          if (ierr /= 0) then
          54           0 :             call endrun('clubb_mf_readnl: ERROR reading namelist')
          55             :          end if
          56             :       end if
          57           2 :       close(iunit)
          58             :     end if
          59             :
          60        1536 :     call mpi_bcast(clubb_mf_L0,   1, mpi_real8,   mstrid, mpicom, ierr)
          61        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_L0")
          62        1536 :     call mpi_bcast(clubb_mf_ent0, 1, mpi_real8,   mstrid, mpicom, ierr)
          63        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_ent0")
          64        1536 :     call mpi_bcast(clubb_mf_nup,  1, mpi_integer, mstrid, mpicom, ierr)
          65        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_nup")
          66        1536 :     call mpi_bcast(do_clubb_mf,      1, mpi_logical, mstrid, mpicom, ierr)
          67        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_clubb_mf")
          68        1536 :     call mpi_bcast(do_clubb_mf_diag, 1, mpi_logical, mstrid, mpicom, ierr)
          69        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_clubb_mf_diag")
          70             :
          71        1536 :     if ((.not. do_clubb_mf) .and. do_clubb_mf_diag ) then
          72           0 :        call endrun('clubb_mf_readnl: Error - cannot turn on do_clubb_mf_diag without also turning on do_clubb_mf')
          73             :     end if
          74             :
          75             :
          76        1536 :   end subroutine clubb_mf_readnl
          77             :
          78           0 :   subroutine integrate_mf( nz,      dzt,     zm,      p_zm,      iexner_zm,         & ! input
          79           0 :                                                       p_zt,      iexner_zt,         & ! input
          80           0 :                            u,       v,       thl,     qt,        thv,               & ! input
          81           0 :                                              thl_zm,  qt_zm,                        & ! input
          82             :                                              wthl,    wqt,       pblh,              & ! input
          83           0 :                            dry_a,   moist_a,                                        & ! output - plume diagnostics
          84           0 :                            dry_w,   moist_w,                                        & ! output - plume diagnostics
          85           0 :                            dry_qt,  moist_qt,                                       & ! output - plume diagnostics
          86           0 :                            dry_thl, moist_thl,                                      & ! output - plume diagnostics
          87           0 :                            dry_u,   moist_u,                                        & ! output - plume diagnostics
          88           0 :                            dry_v,   moist_v,                                        & ! output - plume diagnostics
          89           0 :                                     moist_qc,                                       & ! output - plume diagnostics
          90           0 :                            ae,      aw,                                             & ! output - diagnosed fluxes BEFORE mean field update
          91           0 :                            awthl,   awqt,                                           & ! output - diagnosed fluxes BEFORE mean field update
          92           0 :                            awql,    awqi,                                           & ! output - diagnosed fluxes BEFORE mean field update
          93           0 :                            awu,     awv,                                            & ! output - diagnosed fluxes BEFORE mean field update
          94           0 :                            thlflx,  qtflx )                                           ! output - variables needed for solver
          95             :
          96             :   ! ================================================================================= !
          97             :   ! Mass-flux algorithm                                                               !
          98             :   !                                                                                   !
          99             :   ! Provides rtm and thl fluxes due to mass flux ensemble,                            !
         100             :   ! which are fed into the mixed explicit/implicit clubb solver as explicit terms     !
         101             :   !                                                                                   !
         102             :   ! Variables needed for solver                                                       !
         103             :   ! ae = sum_i (1-a_i)                                                                !
         104             :   ! aw3 = sum (a_i w_i)                                                               !
         105             :   ! aws3 = sum (a_i w_i*s_i); s=thl*cp                                                !
         106             :   ! aws3,awqv3,awql3,awqi3,awu3,awv3 similar as above except for different variables  !
         107             :   !                                                                                   !
         108             :   ! Mass flux variables are computed on edges (i.e. momentum grid):                   !
         109             :   ! upa,upw,upqt,...                                                                  !
         110             :   ! dry_a,moist_a,dry_w,moist_w, ...                                                  !
         111             :   !                                                                                   !
         112             :   ! In CLUBB (unlike CAM) nlevs of momentum grid = nlevs of thermodynamic grid,       !
         113             :   ! due to a subsurface thermodynamic layer. To avoid confusion, below the variables  !
         114             :   ! are grouped by the grid they are on.                                              !
         115             :   !                                                                                   !
         116             :   ! *note that state on the lowest thermo level is equal to state on the lowest       !
         117             :   ! momentum level due to state_zt(1) = state_zt(2), and lowest momentum level        !
         118             :   ! is a weighted combination of the lowest two thermodynamic levels.                 !
         119             :   !                                                                                   !
         120             :   ! ---------------------------------Authors----------------------------------------  !
         121             :   ! Marcin Kurowski, JPL                                                              !
         122             :   ! Modified heavily by Mikael Witte, UCLA/JPL for implementation in CESM2/E3SM       !
         123             :   ! Additional modifications by Adam Herrington, NCAR                                 !
         124             :   ! ================================================================================= !
         125             :
         126             :      use physconst,          only: rair, cpair, gravit, latvap, latice, zvir
         127             :
         128             :      integer,  intent(in)                :: nz
         129             :      real(r8), dimension(nz), intent(in) :: u,      v,            & ! thermodynamic grid
         130             :                                             thl,    thv,          & ! thermodynamic grid
         131             :                                             qt,                   & ! thermodynamic grid
         132             :                                             dzt,                  & ! thermodynamic grid
         133             :                                             p_zt,   iexner_zt,    & ! thermodynamic grid
         134             :                                             thl_zm,               & ! momentum grid
         135             :                                             qt_zm,                & ! momentum grid
         136             :                                             zm,                   & ! momentum grid
         137             :                                             p_zm,   iexner_zm       ! momentum grid
         138             :
         139             :      real(r8), intent(in)                :: wthl,wqt
         140             :      real(r8), intent(inout)             :: pblh
         141             :
         142             :      real(r8),dimension(nz), intent(out) :: dry_a,   moist_a,     & ! momentum grid
         143             :                                             dry_w,   moist_w,     & ! momentum grid
         144             :                                             dry_qt,  moist_qt,    & ! momentum grid
         145             :                                             dry_thl, moist_thl,   & ! momentum grid
         146             :                                             dry_u,   moist_u,     & ! momentum grid
         147             :                                             dry_v,   moist_v,     & ! momentum grid
         148             :                                                      moist_qc,    & ! momentum grid
         149             :                                             ae,      aw,          & ! momentum grid
         150             :                                             awthl,   awqt,        & ! momentum grid
         151             :                                             awql,    awqi,        & ! momentum grid
         152             :                                             awu,     awv,         & ! momentum grid
         153             :                                             thlflx,  qtflx          ! momentum grid
         154             :
         155             :      ! =============================================================================== !
         156             :      ! INTERNAL VARIABLES
         157             :      !
         158             :      ! sums over all plumes
         159             :      real(r8), dimension(nz)              :: moist_th, dry_th,         & ! momentum grid
         160           0 :                                              awqv,     awth              ! momentum grid
         161             :      !
         162             :      ! updraft properties
         163           0 :      real(r8), dimension(nz,clubb_mf_nup) :: upw,      upa,            & ! momentum grid
         164           0 :                                              upqt,     upqc,           & ! momentum grid
         165           0 :                                              upqv,     upqs,           & ! momentum grid
         166           0 :                                              upql,     upqi,           & ! momentum grid
         167           0 :                                              upth,     upthv,          & ! momentum grid
         168           0 :                                                        upthl,          & ! momentum grid
         169           0 :                                              upu,      upv               ! momentum grid
         170             :      !
         171             :      ! entrainment profiles
         172           0 :      real(r8), dimension(nz,clubb_mf_nup) :: ent,      entf              ! thermodynamic grid
         173           0 :      integer,  dimension(nz,clubb_mf_nup) :: enti                        ! thermodynamic grid
         174             :      !
         175             :      ! other variables
         176             :      integer                              :: k,i,ih
         177           0 :      real(r8), dimension(clubb_mf_nup)    :: zcb
         178             :      real(r8)                             :: zcb_unset,                &
         179             :                                              wthv,                     &
         180             :                                              wstar,  qstar,   thvstar, &
         181             :                                              sigmaw, sigmaqt, sigmathv,&
         182             :                                                      wmin,    wmax,    &
         183             :                                              wlv,    wtv,     wp,      &
         184             :                                              B,                        & ! thermodynamic grid
         185             :                                              entexp, entexpu, entw,    & ! thermodynamic grid
         186             :                                              thln,   thvn,    thn,     & ! momentum grid
         187             :                                              qtn,    qsn,              & ! momentum grid
         188             :                                              qcn,    qln,     qin,     & ! momentum grid
         189             :                                              un,     vn,      wn2,     & ! momentum grid
         190             :                                              lmixn,                    & ! momentum grid
         191             :                                              supqt,  supthl              ! thermodynamic grid
         192             :      !
         193             :      ! parameters defining initial conditions for updrafts
         194             :      real(r8),parameter                   :: pwmin = 1.5_r8,           &
         195             :                                              pwmax = 3._r8
         196             :
         197             :      !
         198             :      ! alpha, z-scores after Suselj etal 2019
         199             :      real(r8),parameter                   :: alphw   = 0.572_r8,       &
         200             :                                              alphqt  = 2.890_r8,       &
         201             :                                              alphthv = 2.890_r8
         202             :      !
         203             :      ! w' covariance after Suselj etal 2019
         204             :      real(r8),parameter                   :: cwqt  = 0.32_r8,          &
         205             :                                              cwthv = 0.58_r8
         206             :      !
         207             :      ! virtual mass coefficients for w-eqn after Suselj etal 2019
         208             :      real(r8),parameter                   :: wa = 1.0_r8,              &
         209             :                                              wb = 1.5_r8
         210             :      !
         211             :      ! min values to avoid singularities
         212             :      real(r8),parameter                   :: wstarmin = 1.e-3_r8,      &
         213             :                                              pblhmin  = 100._r8
         214             :      !
         215             :      ! to condensate or not to condensate
         216             :      logical                              :: do_condensation = .true.
         217             :      !
         218             :      ! to precip or not to precip
         219             :      logical                              :: do_precip = .false.
         220             :      !
         221             :      ! to debug flag (overides stochastic entrainment)
         222             :      logical                              :: debug  = .false.
         223             :      real(r8),parameter                   :: fixent = 0.004_r8
         224             :
         225             :      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         226             :      !!!!!!!!!!!!!!!!!!!!!! BEGIN CODE !!!!!!!!!!!!!!!!!!!!!!!
         227             :      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         228             :
         229             :      ! INITIALIZE OUTPUT VARIABLES
         230             :      ! set updraft properties to zero
         231           0 :      dry_a     = 0._r8
         232           0 :      moist_a   = 0._r8
         233           0 :      dry_w     = 0._r8
         234           0 :      moist_w   = 0._r8
         235           0 :      dry_qt    = 0._r8
         236           0 :      moist_qt  = 0._r8
         237           0 :      dry_thl   = 0._r8
         238           0 :      moist_thl = 0._r8
         239           0 :      dry_u     = 0._r8
         240           0 :      moist_u   = 0._r8
         241           0 :      dry_v     = 0._r8
         242           0 :      moist_v   = 0._r8
         243           0 :      moist_qc  = 0._r8
         244             :      ! outputs - variables needed for solver
         245           0 :      aw        = 0._r8
         246           0 :      awthl     = 0._r8
         247           0 :      awqt      = 0._r8
         248           0 :      awqv      = 0._r8
         249           0 :      awql      = 0._r8
         250           0 :      awqi      = 0._r8
         251           0 :      awu       = 0._r8
         252           0 :      awv       = 0._r8
         253           0 :      thlflx    = 0._r8
         254           0 :      qtflx     = 0._r8
         255             :
         256           0 :      ent       = 0._r8
         257           0 :      entf      = 0._r8
         258           0 :      enti      = 0
         259             :
         260             :      ! this is the environmental area - by default 1.
         261           0 :      ae = 1._r8
         262             :
         263             :      ! START MAIN COMPUTATION
         264           0 :      upw   = 0._r8
         265           0 :      upthl = 0._r8
         266           0 :      upthv = 0._r8
         267           0 :      upqt  = 0._r8
         268           0 :      upa   = 0._r8
         269           0 :      upu   = 0._r8
         270           0 :      upv   = 0._r8
         271           0 :      upqc  = 0._r8
         272           0 :      upth  = 0._r8
         273           0 :      upql  = 0._r8
         274           0 :      upqi  = 0._r8
         275           0 :      upqv  = 0._r8
         276           0 :      upqs  = 0._r8
         277             :
         278             :      ! unique identifier
         279             :      zcb_unset = 9999999._r8
         280           0 :      zcb       = zcb_unset
         281             :
         282           0 :      pblh = max(pblh,pblhmin)
         283           0 :      wthv = wthl+zvir*thv(1)*wqt
         284             :
         285             :      ! if surface buoyancy is positive then do mass-flux
         286           0 :      if ( wthv > 0._r8 ) then
         287             :
         288           0 :        if (debug) then
         289             :          ! overide stochastic entrainment with fixent
         290           0 :          ent(:,:) = fixent
         291             :        else
         292             :
         293             :          ! get entrainment coefficient, dz/L0
         294           0 :          do i=1,clubb_mf_nup
         295           0 :            do k=1,nz
         296           0 :              entf(k,i) = dzt(k) / clubb_mf_L0
         297             :            enddo
         298             :          enddo
         299             :
         300             :          ! get poisson, P(dz/L0)
         301           0 :          call poisson( nz, clubb_mf_nup, entf, enti, u(2:5))
         302             :
         303             :          ! get entrainment, ent=ent0/dz*P(dz/L0)
         304           0 :          do i=1,clubb_mf_nup
         305           0 :            do k=1,nz
         306           0 :              ent(k,i) = real( enti(k,i))*clubb_mf_ent0/dzt(k)
         307             :            enddo
         308             :          enddo
         309             :
         310             :        end if
         311             :
         312             :        ! get surface conditions
         313           0 :        wstar   = max( wstarmin, (gravit/thv(1)*wthv*pblh)**(1._r8/3._r8) )
         314           0 :        qstar   = wqt / wstar
         315           0 :        thvstar = wthv / wstar
         316             :
         317           0 :        sigmaw   = alphw * wstar
         318           0 :        sigmaqt  = alphqt * abs(qstar)
         319           0 :        sigmathv = alphthv * abs(thvstar)
         320             :
         321           0 :        wmin = sigmaw * pwmin
         322           0 :        wmax = sigmaw * pwmax
         323             :
         324           0 :        do i=1,clubb_mf_nup
         325             :
         326           0 :          wlv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * (real(i-1, r8))
         327           0 :          wtv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * real(i,r8)
         328             :
         329           0 :          upw(1,i) = 0.5_r8 * (wlv+wtv)
         330             :          upa(1,i) = 0.5_r8 * erf( wtv/(sqrt(2._r8)*sigmaw) ) &
         331           0 :                     - 0.5_r8 * erf( wlv/(sqrt(2._r8)*sigmaw) )
         332             :
         333           0 :          upu(1,i) = u(1)
         334           0 :          upv(1,i) = v(1)
         335             :
         336           0 :          upqt(1,i)  = qt(1)  + cwqt * upw(1,i) * sigmaqt/sigmaw
         337           0 :          upthv(1,i) = thv(1) + cwthv * upw(1,i) * sigmathv/sigmaw
         338           0 :          upthl(1,i) = upthv(1,i) / (1._r8+zvir*upqt(1,i))
         339             :
         340             :          ! get cloud, lowest momentum level
         341           0 :          if (do_condensation) then
         342             :            call condensation_mf(upqt(1,i), upthl(1,i), p_zm(1), iexner_zm(1), &
         343           0 :                                 thvn, qcn, thn, qln, qin, qsn, lmixn)
         344           0 :            upthv(1,i) = thvn
         345           0 :            upqc(1,i)  = qcn
         346           0 :            upql(1,i)  = qln
         347           0 :            upqi(1,i)  = qin
         348           0 :            upqs(1,i)  = qsn
         349           0 :            if (qcn > 0._r8) zcb(i) = zm(1)
         350             :          else
         351             :            ! assume no cldliq
         352           0 :            upqc(1,i)  = 0._r8
         353             :          end if
         354             :
         355             :        enddo
         356             :
         357             :        ! get updraft properties
         358           0 :        do i=1,clubb_mf_nup
         359           0 :          do k=1,nz-1
         360             :
         361             :            ! get microphysics, autoconversion
         362           0 :            if (do_precip .and. upqc(k,i) > 0._r8) then
         363           0 :              call precip_mf(upqs(k,i),upqt(k,i),upw(k,i),dzt(k+1),zm(k+1)-zcb(i),supqt)
         364             :
         365           0 :              supthl = -1._r8*lmixn*supqt*iexner_zt(k+1)/cpair
         366             :            else
         367           0 :              supqt  = 0._r8
         368           0 :              supthl = 0._r8
         369             :            end if
         370             :
         371             :            ! integrate updraft
         372           0 :            entexp  = exp(-ent(k+1,i)*dzt(k+1))
         373           0 :            entexpu = exp(-ent(k+1,i)*dzt(k+1)/3._r8)
         374             :
         375           0 :            qtn  = qt(k+1) *(1._r8-entexp ) + upqt (k,i)*entexp + supqt
         376           0 :            thln = thl(k+1)*(1._r8-entexp ) + upthl(k,i)*entexp + supthl
         377           0 :            un   = u(k+1)  *(1._r8-entexpu) + upu  (k,i)*entexpu
         378           0 :            vn   = v(k+1)  *(1._r8-entexpu) + upv  (k,i)*entexpu
         379             :
         380             :            ! get cloud, momentum levels
         381           0 :            if (do_condensation) then
         382             :              call condensation_mf(qtn, thln, p_zm(k+1), iexner_zm(k+1), &
         383           0 :                                   thvn, qcn, thn, qln, qin, qsn, lmixn)
         384           0 :              if (zcb(i).eq.zcb_unset .and. qcn > 0._r8) zcb(i) = zm(k+1)
         385             :            else
         386           0 :              thvn = thln*(1._r8+zvir*qtn)
         387             :            end if
         388             :
         389             :            ! get buoyancy
         390           0 :            B=gravit*(0.5_r8*(thvn + upthv(k,i))/thv(k+1)-1._r8)
         391           0 :            if (debug) then
         392           0 :              if ( masterproc ) then
         393           0 :                write(iulog,*) "B(k,i), k, i ", B, k, i
         394             :              end if
         395             :            end if
         396             :
         397             :            ! get wn^2
         398           0 :            wp = wb*ent(k+1,i)
         399           0 :            if (wp==0._r8) then
         400           0 :              wn2 = upw(k,i)**2._r8+2._r8*wa*B*dzt(k+1)
         401             :            else
         402           0 :              entw = exp(-2._r8*wp*dzt(k+1))
         403           0 :              wn2 = entw*upw(k,i)**2._r8+wa*B/(wb*ent(k+1,i))*(1._r8-entw)
         404             :            end if
         405             :
         406           0 :            if (wn2>0._r8) then
         407           0 :              upw(k+1,i)   = sqrt(wn2)
         408           0 :              upthv(k+1,i) = thvn
         409           0 :              upthl(k+1,i) = thln
         410           0 :              upqt(k+1,i)  = qtn
         411           0 :              upqc(k+1,i)  = qcn
         412           0 :              upqs(k+1,i)  = qsn
         413           0 :              upu(k+1,i)   = un
         414           0 :              upv(k+1,i)   = vn
         415           0 :              upa(k+1,i)   = upa(k,i)
         416           0 :              upql(k+1,i)  = qln
         417           0 :              upqi(k+1,i)  = qin
         418           0 :              upqv(k+1,i)  = qtn - qcn
         419             :            else
         420             :              exit
         421             :            end if
         422             :          enddo
         423             :        enddo
         424             :
         425             :        ! writing updraft properties for output
         426           0 :        do k=1,nz
         427             :
         428             :          ! first sum over all i-updrafts
         429           0 :          do i=1,clubb_mf_nup
         430           0 :            if (upqc(k,i)>0._r8) then
         431           0 :              moist_a(k)   = moist_a(k)   + upa(k,i)
         432           0 :              moist_w(k)   = moist_w(k)   + upa(k,i)*upw(k,i)
         433           0 :              moist_qt(k)  = moist_qt(k)  + upa(k,i)*upqt(k,i)
         434           0 :              moist_thl(k) = moist_thl(k) + upa(k,i)*upthl(k,i)
         435           0 :              moist_u(k)   = moist_u(k)   + upa(k,i)*upu(k,i)
         436           0 :              moist_v(k)   = moist_v(k)   + upa(k,i)*upv(k,i)
         437           0 :              moist_qc(k)  = moist_qc(k)  + upa(k,i)*upqc(k,i)
         438             :            else
         439           0 :              dry_a(k)     = dry_a(k)     + upa(k,i)
         440           0 :              dry_w(k)     = dry_w(k)     + upa(k,i)*upw(k,i)
         441           0 :              dry_qt(k)    = dry_qt(k)    + upa(k,i)*upqt(k,i)
         442           0 :              dry_thl(k)   = dry_thl(k)   + upa(k,i)*upthl(k,i)
         443           0 :              dry_u(k)     = dry_u(k)     + upa(k,i)*upu(k,i)
         444           0 :              dry_v(k)     = dry_v(k)     + upa(k,i)*upv(k,i)
         445             :            endif
         446             :          enddo
         447             :
         448           0 :          if ( dry_a(k) > 0._r8 ) then
         449           0 :            dry_w(k)   = dry_w(k)   / dry_a(k)
         450           0 :            dry_qt(k)  = dry_qt(k)  / dry_a(k)
         451           0 :            dry_thl(k) = dry_thl(k) / dry_a(k)
         452           0 :            dry_u(k)   = dry_u(k)   / dry_a(k)
         453           0 :            dry_v(k)   = dry_v(k)   / dry_a(k)
         454             :          else
         455           0 :            dry_w(k)   = 0._r8
         456           0 :            dry_qt(k)  = 0._r8
         457           0 :            dry_thl(k) = 0._r8
         458           0 :            dry_u(k)   = 0._r8
         459           0 :            dry_v(k)   = 0._r8
         460             :          endif
         461             :
         462           0 :          if ( moist_a(k) > 0._r8 ) then
         463           0 :            moist_w(k)   = moist_w(k)   / moist_a(k)
         464           0 :            moist_qt(k)  = moist_qt(k)  / moist_a(k)
         465           0 :            moist_thl(k) = moist_thl(k) / moist_a(k)
         466           0 :            moist_u(k)   = moist_u(k)   / moist_a(k)
         467           0 :            moist_v(k)   = moist_v(k)   / moist_a(k)
         468           0 :            moist_qc(k)  = moist_qc(k)  / moist_a(k)
         469             :          else
         470           0 :            moist_w(k)   = 0._r8
         471           0 :            moist_qt(k)  = 0._r8
         472           0 :            moist_thl(k) = 0._r8
         473           0 :            moist_u(k)   = 0._r8
         474           0 :            moist_v(k)   = 0._r8
         475           0 :            moist_qc(k)  = 0._r8
         476             :          endif
         477             :
         478             :        enddo
         479             :
         480           0 :        do k=1,nz
         481           0 :          do i=1,clubb_mf_nup
         482           0 :            ae  (k) = ae  (k) - upa(k,i)
         483           0 :            aw  (k) = aw  (k) + upa(k,i)*upw(k,i)
         484           0 :            awu (k) = awu (k) + upa(k,i)*upw(k,i)*upu(k,i)
         485           0 :            awv (k) = awv (k) + upa(k,i)*upw(k,i)*upv(k,i)
         486           0 :            awthl(k)= awthl(k)+ upa(k,i)*upw(k,i)*upthl(k,i)
         487           0 :            awqt(k) = awqt(k) + upa(k,i)*upw(k,i)*upqt(k,i)
         488           0 :            awqv(k) = awqv(k) + upa(k,i)*upw(k,i)*upqv(k,i)
         489           0 :            awql(k) = awql(k) + upa(k,i)*upw(k,i)*upql(k,i)
         490           0 :            awqi(k) = awqi(k) + upa(k,i)*upw(k,i)*upqi(k,i)
         491             :          enddo
         492             :        enddo
         493             :
         494           0 :        do k=1,nz
         495           0 :          thlflx(k)= awthl(k) - aw(k)*thl_zm(k)
         496           0 :          qtflx(k)= awqt(k) - aw(k)*qt_zm(k)
         497             :        enddo
         498             :
         499             :      end if  ! ( wthv > 0.0 )
         500             :
         501           0 :   end subroutine integrate_mf
         502             :
         503           0 :   subroutine condensation_mf( qt, thl, p, iex, thv, qc, th, ql, qi, qs, lmix )
         504             :   ! =============================================================================== !
         505             :   ! zero or one condensation for edmf: calculates thv and qc                        !
         506             :   ! =============================================================================== !
         507             :      use physconst,          only: cpair, zvir, h2otrip
         508             :      use wv_saturation,      only : qsat
         509             :
         510             :      real(r8),intent(in) :: qt,thl,p,iex
         511             :      real(r8),intent(out):: thv,qc,th,ql,qi,qs,lmix
         512             :
         513             :      !local variables
         514             :      integer  :: niter,i
         515             :      real(r8) :: diff,t,qstmp,qcold,es,wf
         516             :
         517             :      ! max number of iterations
         518           0 :      niter=50
         519             :      ! minimum difference
         520           0 :      diff=2.e-5_r8
         521             :
         522           0 :      qc=0._r8
         523           0 :      t=thl/iex
         524             :
         525             :      !by definition:
         526             :      ! T   = Th*Exner, Exner=(p/p0)^(R/cp)   (1)
         527             :      ! Thl = Th - L/cp*ql/Exner              (2)
         528             :      !so:
         529             :      ! Th  = Thl + L/cp*ql/Exner             (3)
         530             :      ! T   = Th*Exner=(Thl+L/cp*ql/Exner)*Exner    (4)
         531             :      !     = Thl*Exner + L/cp*ql
         532           0 :      do i=1,niter
         533           0 :        wf = get_watf(t)
         534           0 :        t = thl/iex+get_alhl(wf)/cpair*qc   !as in (4)
         535             :
         536             :        ! qsat, p is in pascal (check!)
         537           0 :        call qsat(t,p,es,qstmp)
         538           0 :        qcold = qc
         539           0 :        qc = max(0.5_r8*qc+0.5_r8*(qt-qstmp),0._r8)
         540           0 :        if (abs(qc-qcold)<diff) exit
         541             :      enddo
         542             :
         543           0 :      wf = get_watf(t)
         544           0 :      t = thl/iex+get_alhl(wf)/cpair*qc
         545             :
         546           0 :      call qsat(t,p,es,qs)
         547           0 :      qc = max(qt-qs,0._r8)
         548           0 :      thv = (thl+get_alhl(wf)/cpair*iex*qc)*(1._r8+zvir*(qt-qc)-qc)
         549           0 :      lmix = get_alhl(wf)
         550           0 :      th = t*iex
         551           0 :      qi = qc*(1._r8-wf)
         552           0 :      ql = qc*wf
         553             :
         554             :      contains
         555             :
         556           0 :      function get_watf(t)
         557             :        real(r8)            :: t,get_watf,tc
         558             :        real(r8), parameter :: &
         559             :                               tmax=-10._r8, &
         560             :                               tmin=-40._r8
         561             :
         562           0 :        tc=t-h2otrip
         563             :
         564           0 :        if (tc>tmax) then
         565             :          get_watf=1._r8
         566           0 :        else if (tc<tmin) then
         567             :          get_watf=0._r8
         568             :        else
         569           0 :          get_watf=(tc-tmin)/(tmax-tmin);
         570             :        end if
         571             :
         572           0 :      end function get_watf
         573             :
         574             :
         575           0 :      function get_alhl(wf)
         576             :      !latent heat of the mixture based on water fraction
         577             :        use physconst,        only : latvap , latice
         578             :        real(r8) :: get_alhl,wf
         579             :
         580           0 :        get_alhl = wf*latvap+(1._r8-wf)*(latvap+latice)
         581             :
         582           0 :      end function get_alhl
         583             :
         584             :   end subroutine condensation_mf
         585             :
         586           0 :   subroutine precip_mf(qs,qt,w,dz,dzcld,Supqt)
         587             :   !**********************************************************************
         588             :   ! Precipitation microphysics
         589             :   ! By Adam Herrington, after Kay Suselj
         590             :   !**********************************************************************
         591             :
         592             :        real(r8),intent(in)  :: qs,qt,w,dz,dzcld
         593             :        real(r8),intent(out) :: Supqt
         594             :        !
         595             :        ! local vars
         596             :        real(r8)            :: tauwgt, tau,       & ! time-scale vars
         597             :                               qstar                ! excess cloud liquid
         598             :
         599             :        real(r8),parameter  :: tau0  = 15._r8,    & ! base time-scale
         600             :                               zmin  = 300._r8,   & ! small cloud thick
         601             :                               zmax  = 3000._r8,  & ! large cloud thick
         602             :                               qcmin = 0.00125_r8   ! supersat threshold
         603             :
         604           0 :        qstar = qs+qcmin
         605             :
         606           0 :        if (qt > qstar) then
         607             :          ! get precip efficiency
         608           0 :          tauwgt = (dzcld-zmin)/(zmax-zmin)
         609           0 :          tauwgt = min(max(tauwgt,0._r8),1._r8)
         610           0 :          tau    = tauwgt/tau0
         611             :
         612             :          ! get source for updraft
         613           0 :          Supqt = (qstar-qt)*(1._r8 - exp(-1._r8*tau*dz/w))
         614             :        else
         615           0 :          Supqt = 0._r8
         616             :        end if
         617             :
         618           0 :   end subroutine precip_mf
         619             :
         620           0 :   subroutine poisson(nz,nup,lambda,poi,state)
         621             :   !**********************************************************************
         622             :   ! Set a unique (but reproduceble) seed for the kiss RNG
         623             :   ! Call Poisson deviate
         624             :   ! By Adam Herrington
         625             :   !**********************************************************************
         626             :    use shr_RandNum_mod, only: ShrKissRandGen
         627             :
         628             :        integer,                     intent(in)  :: nz,nup
         629             :        real(r8), dimension(4),      intent(in)  :: state
         630             :        real(r8), dimension(nz,nup), intent(in)  :: lambda
         631             :        integer,  dimension(nz,nup), intent(out) :: poi
         632             :        integer,  dimension(1,4)                 :: tmpseed
         633             :        integer                                  :: i,j
         634           0 :        type(ShrKissRandGen)                     :: kiss_gen
         635             :
         636             :        ! Compute seed
         637           0 :        tmpseed(1,1) = int((state(1) - int(state(1))) * 1000000000._r8)
         638           0 :        tmpseed(1,2) = int((state(2) - int(state(2))) * 1000000000._r8)
         639           0 :        tmpseed(1,3) = int((state(3) - int(state(3))) * 1000000000._r8)
         640           0 :        tmpseed(1,4) = int((state(4) - int(state(4))) * 1000000000._r8)
         641             :
         642             :        ! Set seed
         643           0 :        kiss_gen = ShrKissRandGen(tmpseed)
         644             :
         645           0 :        do i=1,nz
         646           0 :          do j=1,nup
         647           0 :            call knuth(kiss_gen,lambda(i,j),poi(i,j))
         648             :          enddo
         649             :        enddo
         650             :
         651           0 :   end subroutine poisson
         652             :
         653           0 :   subroutine knuth(kiss_gen,lambda,kout)
         654             :   !**********************************************************************
         655             :   ! Discrete random poisson from Knuth
         656             :   ! The Art of Computer Programming, v2, 137-138
         657             :   ! By Adam Herrington
         658             :   !**********************************************************************
         659             :    use shr_RandNum_mod, only: ShrKissRandGen
         660             :
         661             :        type(ShrKissRandGen), intent(inout) :: kiss_gen
         662             :        real(r8),             intent(in)    :: lambda
         663             :        integer,              intent(out)   :: kout
         664             :
         665             :        ! Local variables
         666             :        real(r8), dimension(1,1) :: tmpuni
         667             :        real(r8)                 :: puni, explam
         668             :        integer                  :: k
         669             :
         670           0 :        k = 0
         671           0 :        explam = exp(-1._r8*lambda)
         672           0 :        puni = 1._r8
         673           0 :        do while (puni > explam)
         674           0 :          k = k + 1
         675           0 :          call kiss_gen%random(tmpuni)
         676           0 :          puni = puni*tmpuni(1,1)
         677             :        end do
         678           0 :        kout = k - 1
         679             :
         680           0 :   end subroutine knuth
         681             :
         682             : end module clubb_mf

Generated by: LCOV version 1.14