LCOV - code coverage report
Current view: | top level - physics/cam - rk_stratiform.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 17 | 409 | 4.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 6 | 16.7 %

              Line data    Source code
           1             : module rk_stratiform
           2             :
           3             : !-------------------------------------------------------------------------------------------------------
           4             : !
           5             : ! Provides the CAM interface to the Rasch and Kristjansson (RK)
           6             : ! prognostic cloud microphysics, and the cam4 macrophysics.
           7             : !
           8             : !-------------------------------------------------------------------------------------------------------
           9             :
          10             : use shr_kind_mod,   only: r8=>shr_kind_r8
          11             : use ppgrid,         only: pcols, pver, pverp
          12             : use physconst,      only: gravit, latvap, latice
          13             : use phys_control,   only: phys_getopts
          14             : use constituents,   only: pcnst
          15             : use spmd_utils,     only: masterproc
          16             :
          17             : use cam_logfile,    only: iulog
          18             : use cam_abortutils, only: endrun
          19             : use perf_mod
          20             :
          21             : implicit none
          22             : private
          23             : save
          24             :
          25             : public :: rk_stratiform_register, rk_stratiform_init_cnst, rk_stratiform_implements_cnst
          26             : public :: rk_stratiform_init
          27             : public :: rk_stratiform_tend
          28             : public :: rk_stratiform_readnl
          29             :
          30             : ! Physics buffer indices
          31             : integer  ::  landm_idx          = 0
          32             :
          33             : integer  ::  qcwat_idx          = 0
          34             : integer  ::  lcwat_idx          = 0
          35             : integer  ::  tcwat_idx          = 0
          36             :
          37             : integer  ::  cld_idx            = 0
          38             : integer  ::  ast_idx            = 0
          39             : integer  ::  concld_idx         = 0
          40             : integer  ::  fice_idx           = 0
          41             :
          42             : integer  ::  qme_idx            = 0
          43             : integer  ::  prain_idx          = 0
          44             : integer  ::  nevapr_idx         = 0
          45             :
          46             : integer  ::  wsedl_idx          = 0
          47             :
          48             : integer  ::  rei_idx            = 0
          49             : integer  ::  rel_idx            = 0
          50             :
          51             : integer  ::  shfrc_idx          = 0
          52             : integer  ::  cmfmc_sh_idx       = 0
          53             :
          54             : integer  ::  prec_str_idx       = 0
          55             : integer  ::  snow_str_idx       = 0
          56             : integer  ::  prec_sed_idx       = 0
          57             : integer  ::  snow_sed_idx       = 0
          58             : integer  ::  prec_pcw_idx       = 0
          59             : integer  ::  snow_pcw_idx       = 0
          60             :
          61             : integer  ::  ls_flxprc_idx      = 0
          62             : integer  ::  ls_flxsnw_idx      = 0
          63             :
          64             : integer, parameter :: ncnst = 2                    ! Number of constituents
          65             : character(len=8), dimension(ncnst), parameter &    ! Constituent names
          66             :                    :: cnst_names = (/'CLDLIQ', 'CLDICE'/)
          67             : logical            :: use_shfrc                       ! Local copy of flag from convect_shallow_use_shfrc
          68             :
          69             : logical            :: do_cnst = .false. ! True when this module has registered constituents.
          70             :
          71             : integer :: &
          72             :    ixcldliq,     &! cloud liquid amount index
          73             :    ixcldice       ! cloud ice amount index
          74             :
          75             : real(r8), parameter :: unset_r8 = huge(1.0_r8)
          76             : logical :: do_psrhmin
          77             :
          78             : !===============================================================================
          79             : contains
          80             : !===============================================================================
          81        1536 :   subroutine rk_stratiform_readnl(nlfile)
          82             :
          83             :    use namelist_utils,  only: find_group_name
          84             :    use units,           only: getunit, freeunit
          85             :    use cldwat,          only: cldwat_init
          86             :    use mpishorthand
          87             :
          88             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          89             :
          90             :    ! Local variables
          91             :    integer :: unitn, ierr
          92             :    character(len=*), parameter :: subname = 'rk_stratiform_readnl'
          93             :
          94             :    ! Namelist variables
          95             :    real(r8) :: rk_strat_icritw  = unset_r8    !   icritw  = threshold for autoconversion of warm ice
          96             :    real(r8) :: rk_strat_icritc  = unset_r8    !   icritc  = threshold for autoconversion of cold ice
          97             :    real(r8) :: rk_strat_conke   = unset_r8    !   conke   = tunable constant for evaporation of precip
          98             :    real(r8) :: rk_strat_r3lcrit = unset_r8    !   r3lcrit = critical radius where liq conversion begins
          99             :    real(r8) :: rk_strat_polstrat_rhmin = unset_r8 ! condensation threadhold in polar stratosphere
         100             :
         101             :    namelist /rk_stratiform_nl/ rk_strat_icritw, rk_strat_icritc, rk_strat_conke, rk_strat_r3lcrit
         102             :    namelist /rk_stratiform_nl/ rk_strat_polstrat_rhmin
         103             :
         104             :    !-----------------------------------------------------------------------------
         105             :
         106        1536 :    if (masterproc) then
         107           2 :       unitn = getunit()
         108           2 :       open( unitn, file=trim(nlfile), status='old' )
         109           2 :       call find_group_name(unitn, 'rk_stratiform_nl', status=ierr)
         110           2 :       if (ierr == 0) then
         111           0 :          read(unitn, rk_stratiform_nl, iostat=ierr)
         112           0 :          if (ierr /= 0) then
         113           0 :             call endrun(subname // ':: ERROR reading namelist')
         114             :          end if
         115             :       end if
         116           2 :       close(unitn)
         117           2 :       call freeunit(unitn)
         118             :    end if
         119             :
         120             : #ifdef SPMD
         121             :    ! Broadcast namelist variables
         122        1536 :    call mpibcast(rk_strat_icritw,  1, mpir8, 0, mpicom)
         123        1536 :    call mpibcast(rk_strat_icritc,  1, mpir8, 0, mpicom)
         124        1536 :    call mpibcast(rk_strat_conke,   1, mpir8, 0, mpicom)
         125        1536 :    call mpibcast(rk_strat_r3lcrit, 1, mpir8, 0, mpicom)
         126        1536 :    call mpibcast(rk_strat_polstrat_rhmin, 1, mpir8, 0, mpicom)
         127             : #endif
         128             :
         129        1536 :    do_psrhmin = rk_strat_polstrat_rhmin .ne. unset_r8
         130             :
         131             :    call cldwat_init(rk_strat_icritw,rk_strat_icritc,rk_strat_conke,&
         132        1536 :                     rk_strat_r3lcrit,rk_strat_polstrat_rhmin,do_psrhmin)
         133             :
         134        1536 : end subroutine rk_stratiform_readnl
         135             :
         136           0 : subroutine rk_stratiform_register
         137             :
         138             :    !---------------------------------------------------------------------- !
         139             :    !                                                                       !
         140             :    ! Register the constituents (cloud liquid and cloud ice) and the fields !
         141             :    ! in the physics buffer.                                                !
         142             :    !                                                                       !
         143             :    !---------------------------------------------------------------------- !
         144             :
         145        1536 :    use constituents, only: cnst_add, pcnst
         146             :    use physconst,    only: mwh2o, cpair
         147             :
         148             :    use physics_buffer, only : pbuf_add_field, dtype_r8, dyn_time_lvls
         149             :
         150             :    !-----------------------------------------------------------------------
         151             :
         152             :    ! Take note of the fact that we are registering constituents.
         153           0 :    do_cnst = .true.
         154             :
         155             :    ! Register cloud water and save indices.
         156             :    call cnst_add(cnst_names(1), mwh2o, cpair, 0._r8, ixcldliq, &
         157           0 :       longname='Grid box averaged cloud liquid amount', is_convtran1=.true.)
         158             :    call cnst_add(cnst_names(2), mwh2o, cpair, 0._r8, ixcldice, &
         159           0 :       longname='Grid box averaged cloud ice amount', is_convtran1=.true.)
         160             :
         161           0 :    call pbuf_add_field('QCWAT',  'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), qcwat_idx)
         162           0 :    call pbuf_add_field('LCWAT',  'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), lcwat_idx)
         163           0 :    call pbuf_add_field('TCWAT',  'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), tcwat_idx)
         164             :
         165           0 :    call pbuf_add_field('CLD',    'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), cld_idx)
         166           0 :    call pbuf_add_field('AST',    'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), ast_idx)
         167           0 :    call pbuf_add_field('CONCLD', 'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), concld_idx)
         168             :
         169           0 :    call pbuf_add_field('FICE',   'physpkg', dtype_r8, (/pcols,pver/), fice_idx)
         170             :
         171           0 :    call pbuf_add_field('QME',       'physpkg', dtype_r8, (/pcols,pver/), qme_idx)
         172           0 :    call pbuf_add_field('PRAIN',     'physpkg', dtype_r8, (/pcols,pver/), prain_idx)
         173           0 :    call pbuf_add_field('NEVAPR',    'physpkg', dtype_r8, (/pcols,pver/), nevapr_idx)
         174             :
         175           0 :    call pbuf_add_field('WSEDL',     'physpkg', dtype_r8, (/pcols,pver/), wsedl_idx)
         176             :
         177           0 :    call pbuf_add_field('REI',       'physpkg', dtype_r8, (/pcols,pver/), rei_idx)
         178           0 :    call pbuf_add_field('REL',       'physpkg', dtype_r8, (/pcols,pver/), rel_idx)
         179             :
         180           0 :    call pbuf_add_field('LS_FLXPRC', 'physpkg', dtype_r8, (/pcols,pverp/), ls_flxprc_idx)
         181           0 :    call pbuf_add_field('LS_FLXSNW', 'physpkg', dtype_r8, (/pcols,pverp/), ls_flxsnw_idx)
         182             :
         183           0 : end subroutine rk_stratiform_register
         184             :
         185             : !===============================================================================
         186             :
         187           0 : function rk_stratiform_implements_cnst(name)
         188             :
         189             :   !----------------------------------------------------------------------------- !
         190             :   !                                                                              !
         191             :   ! Return true if specified constituent is implemented by this package          !
         192             :   !                                                                              !
         193             :   !----------------------------------------------------------------------------- !
         194             :
         195             :    character(len=*), intent(in) :: name      ! constituent name
         196             :    logical :: rk_stratiform_implements_cnst     ! return value
         197             :
         198             :    !-----------------------------------------------------------------------
         199             :
         200           0 :    rk_stratiform_implements_cnst = (do_cnst .and. any(name == cnst_names))
         201             :
         202           0 : end function rk_stratiform_implements_cnst
         203             :
         204             : !===============================================================================
         205             :
         206           0 : subroutine rk_stratiform_init_cnst(name, latvals, lonvals, mask, q)
         207             :
         208             :    !----------------------------------------------------------------------- !
         209             :    !                                                                        !
         210             :    ! Initialize the cloud water mixing ratios (liquid and ice), if they are !
         211             :    ! not read from the initial file                                         !
         212             :    !                                                                        !
         213             :    !----------------------------------------------------------------------- !
         214             :
         215             :    character(len=*), intent(in)  :: name       ! constituent name
         216             :    real(r8),         intent(in)  :: latvals(:) ! lat in degrees (ncol)
         217             :    real(r8),         intent(in)  :: lonvals(:) ! lon in degrees (ncol)
         218             :    logical,          intent(in)  :: mask(:)    ! Only initialize where .true.
         219             :    real(r8),         intent(out) :: q(:,:)     ! kg tracer/kg dry air (gcol, plev
         220             :    !-----------------------------------------------------------------------
         221             :    integer :: k
         222             :
         223           0 :    if (any(name == cnst_names)) then
         224           0 :      do k = 1, size(q, 2)
         225           0 :        where(mask)
         226           0 :          q(:, k) = 0.0_r8
         227             :        end where
         228             :      end do
         229             :    end if
         230             :
         231           0 : end subroutine rk_stratiform_init_cnst
         232             :
         233             : !===============================================================================
         234             :
         235           0 : subroutine rk_stratiform_init()
         236             :
         237             :    !-------------------------------------------- !
         238             :    !                                             !
         239             :    ! Initialize the cloud water parameterization !
         240             :    !                                             !
         241             :    !-------------------------------------------- !
         242             :
         243             :    use physics_buffer,  only: physics_buffer_desc, pbuf_get_index
         244             :    use constituents,    only: cnst_get_ind, cnst_name, cnst_longname, sflxnam, apcnst, bpcnst
         245             :    use cam_history,     only: addfld, add_default, horiz_only
         246             :    use convect_shallow, only: convect_shallow_use_shfrc
         247             :    use phys_control,    only: cam_physpkg_is
         248             :    use physconst,       only: tmelt, rhodair, rh2o
         249             :    use cldwat,          only: inimc
         250             :
         251             :    integer :: m, mm
         252             :    logical :: history_amwg         ! output the variables used by the AMWG diag package
         253             :    logical :: history_aerosol      ! Output the MAM aerosol tendencies
         254             :    logical :: history_budget       ! Output tendencies and state variables for CAM4
         255             :                                    ! temperature, water vapor, cloud ice and cloud
         256             :                                    ! liquid budgets.
         257             :    integer :: history_budget_histfile_num ! output history file number for budget fields
         258             :    !-----------------------------------------------------------------------
         259             :
         260             :    call phys_getopts( history_aerosol_out        = history_aerosol      , &
         261             :                       history_amwg_out   = history_amwg                 , &
         262             :                       history_budget_out         = history_budget       , &
         263           0 :                       history_budget_histfile_num_out = history_budget_histfile_num)
         264             :
         265           0 :    landm_idx = pbuf_get_index("LANDM")
         266             :
         267             :    ! Find out whether shfrc from convect_shallow will be used in cldfrc
         268           0 :    if( convect_shallow_use_shfrc() ) then
         269           0 :       use_shfrc = .true.
         270           0 :       shfrc_idx = pbuf_get_index('shfrc')
         271             :    else
         272           0 :       use_shfrc = .false.
         273             :    endif
         274             :
         275             :    ! Register history variables
         276             :
         277           0 :    do m = 1, ncnst
         278           0 :       call cnst_get_ind( cnst_names(m), mm )
         279           0 :       call addfld( cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg',   cnst_longname(mm)                    )
         280           0 :       call addfld( sflxnam  (mm), horiz_only,  'A', 'kg/m2/s', trim(cnst_name(mm))//' surface flux' )
         281           0 :       if (history_amwg) then
         282           0 :          call add_default( cnst_name(mm), 1, ' ' )
         283           0 :          call add_default( sflxnam  (mm), 1, ' ' )
         284             :       endif
         285             :    enddo
         286             :
         287           0 :    call addfld (apcnst(ixcldliq), (/ 'lev' /), 'A', 'kg/kg', trim(cnst_name(ixcldliq))//' after physics'  )
         288           0 :    call addfld (apcnst(ixcldice), (/ 'lev' /), 'A', 'kg/kg', trim(cnst_name(ixcldice))//' after physics'  )
         289           0 :    call addfld (bpcnst(ixcldliq), (/ 'lev' /), 'A', 'kg/kg', trim(cnst_name(ixcldliq))//' before physics' )
         290           0 :    call addfld (bpcnst(ixcldice), (/ 'lev' /), 'A', 'kg/kg', trim(cnst_name(ixcldice))//' before physics' )
         291             :
         292           0 :    if( history_budget) then
         293           0 :       call add_default (cnst_name(ixcldliq), history_budget_histfile_num, ' ')
         294           0 :       call add_default (cnst_name(ixcldice), history_budget_histfile_num, ' ')
         295           0 :       call add_default (apcnst   (ixcldliq), history_budget_histfile_num, ' ')
         296           0 :       call add_default (apcnst   (ixcldice), history_budget_histfile_num, ' ')
         297           0 :       call add_default (bpcnst   (ixcldliq), history_budget_histfile_num, ' ')
         298           0 :       call add_default (bpcnst   (ixcldice), history_budget_histfile_num, ' ')
         299             :    end if
         300             :
         301           0 :    call addfld ('FWAUT',      (/ 'lev' /),   'A', 'fraction', 'Relative importance of liquid autoconversion'            )
         302           0 :    call addfld ('FSAUT',      (/ 'lev' /),   'A', 'fraction', 'Relative importance of ice autoconversion'               )
         303           0 :    call addfld ('FRACW',      (/ 'lev' /),   'A', 'fraction', 'Relative importance of rain accreting liquid'            )
         304           0 :    call addfld ('FSACW',      (/ 'lev' /),   'A', 'fraction', 'Relative importance of snow accreting liquid'            )
         305           0 :    call addfld ('FSACI',      (/ 'lev' /),   'A', 'fraction', 'Relative importance of snow accreting ice'               )
         306           0 :    call addfld ('CME',        (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of cond-evap within the cloud'                      )
         307           0 :    call addfld ('CMEICE',     (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of cond-evap of ice within the cloud'               )
         308           0 :    call addfld ('CMELIQ',     (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of cond-evap of liq within the cloud'               )
         309           0 :    call addfld ('ICE2PR',     (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of conversion of ice to precip'                     )
         310           0 :    call addfld ('LIQ2PR',     (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of conversion of liq to precip'                     )
         311           0 :    call addfld ('ZMDLF',      (/ 'lev' /),   'A', 'kg/kg/s' , 'Detrained liquid water from ZM convection'               )
         312           0 :    call addfld ('SHDLF',      (/ 'lev' /),   'A', 'kg/kg/s' , 'Detrained liquid water from shallow convection'          )
         313             :
         314           0 :    call addfld ('PRODPREC',   (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of conversion of condensate to precip'              )
         315           0 :    call addfld ('EVAPPREC',   (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of evaporation of falling precip'                   )
         316           0 :    call addfld ('EVAPSNOW',   (/ 'lev' /),   'A', 'kg/kg/s' , 'Rate of evaporation of falling snow'                     )
         317           0 :    call addfld ('HPROGCLD',   (/ 'lev' /),   'A', 'W/kg'    , 'Heating from prognostic clouds'                          )
         318           0 :    call addfld ('HCME',       (/ 'lev' /),   'A', 'W/kg'    , 'Heating from cond-evap within the cloud'                 )
         319           0 :    call addfld ('HEVAP',      (/ 'lev' /),   'A', 'W/kg'    , 'Heating from evaporation of falling precip'              )
         320           0 :    call addfld ('HFREEZ',     (/ 'lev' /),   'A', 'W/kg'    , 'Heating rate due to freezing of precip'                  )
         321           0 :    call addfld ('HMELT',      (/ 'lev' /),   'A', 'W/kg'    , 'Heating from snow melt'                                  )
         322           0 :    call addfld ('HREPART',    (/ 'lev' /),   'A', 'W/kg'    , 'Heating from cloud ice/liquid repartitioning'            )
         323           0 :    call addfld ('REPARTICE',  (/ 'lev' /),   'A', 'kg/kg/s' , 'Cloud ice tendency from cloud ice/liquid repartitioning' )
         324           0 :    call addfld ('REPARTLIQ',  (/ 'lev' /),   'A', 'kg/kg/s' , 'Cloud liq tendency from cloud ice/liquid repartitioning' )
         325           0 :    call addfld ('FICE',       (/ 'lev' /),   'A', 'fraction', 'Fractional ice content within cloud'                     )
         326           0 :    call addfld ('ICWMR',      (/ 'lev' /),   'A', 'kg/kg'   , 'Prognostic in-cloud water mixing ratio'                  )
         327           0 :    call addfld ('ICIMR',      (/ 'lev' /),   'A', 'kg/kg'   , 'Prognostic in-cloud ice mixing ratio'                    )
         328           0 :    call addfld ('PCSNOW',     horiz_only   , 'A', 'm/s'     , 'Snow fall from prognostic clouds'                        )
         329             :
         330           0 :    call addfld ('DQSED',      (/ 'lev' /),   'A', 'kg/kg/s' , 'Water vapor tendency from cloud sedimentation'           )
         331           0 :    call addfld ('DLSED',      (/ 'lev' /),   'A', 'kg/kg/s' , 'Cloud liquid tendency from sedimentation'                )
         332           0 :    call addfld ('DISED',      (/ 'lev' /),   'A', 'kg/kg/s' , 'Cloud ice tendency from sedimentation'                   )
         333           0 :    call addfld ('HSED',       (/ 'lev' /),   'A', 'W/kg'    , 'Heating from cloud sediment evaporation'                 )
         334           0 :    call addfld ('SNOWSED',    horiz_only,    'A', 'm/s'     , 'Snow from cloud ice sedimentation'                       )
         335           0 :    call addfld ('RAINSED',    horiz_only,    'A', 'm/s'     , 'Rain from cloud liquid sedimentation'                    )
         336           0 :    call addfld ('PRECSED',    horiz_only,    'A', 'm/s'     , 'Precipitation from cloud sedimentation'                  )
         337             :
         338             :
         339           0 :    call addfld ('CNVCLD',     horiz_only,    'A', 'fraction', 'Vertically integrated convective cloud amount'           )
         340           0 :    call addfld ('CLDST',      (/ 'lev' /),   'A', 'fraction', 'Stratus cloud fraction'                                  )
         341           0 :    call addfld ('CONCLD',     (/ 'lev' /),   'A', 'fraction', 'Convective cloud cover'                                  )
         342             :
         343           0 :    call addfld ('AST',        (/ 'lev' /),   'A','fraction' , 'Stratus cloud fraction'                                  )
         344           0 :    call addfld ('LIQCLDF',    (/ 'lev' /),   'A', 'fraction', 'Stratus Liquid cloud fraction'                           )
         345           0 :    call addfld ('ICECLDF',    (/ 'lev' /),   'A', 'fraction', 'Stratus ICE cloud fraction'                              )
         346           0 :    call addfld ('IWC',        (/ 'lev' /),   'A', 'kg/m3'   , 'Grid box average ice water content'                      )
         347           0 :    call addfld ('LWC',        (/ 'lev' /),   'A', 'kg/m3'   , 'Grid box average liquid water content'                   )
         348           0 :    call addfld ('ICWNC',      (/ 'lev' /),   'A', 'm-3'     , 'Prognostic in-cloud water number conc'                   )
         349           0 :    call addfld ('ICINC',      (/ 'lev' /),   'A', 'm-3'     , 'Prognostic in-cloud ice number conc'                     )
         350           0 :    call addfld ('REL',        (/ 'lev' /),   'A', 'micron'  , 'effective liquid drop radius'                            )
         351           0 :    call addfld ('REI',        (/ 'lev' /),   'A', 'micron'  , 'effective ice particle radius'                           )
         352             :
         353           0 :    if ( history_budget ) then
         354             :
         355           0 :       call add_default ('EVAPSNOW ', history_budget_histfile_num, ' ')
         356           0 :       call add_default ('EVAPPREC ', history_budget_histfile_num, ' ')
         357           0 :       call add_default ('CMELIQ   ', history_budget_histfile_num, ' ')
         358             :
         359           0 :       if( cam_physpkg_is('cam4') ) then
         360             :
         361           0 :          call add_default ('ZMDLF    ', history_budget_histfile_num, ' ')
         362           0 :          call add_default ('CME      ', history_budget_histfile_num, ' ')
         363           0 :          call add_default ('DQSED    ', history_budget_histfile_num, ' ')
         364           0 :          call add_default ('DISED    ', history_budget_histfile_num, ' ')
         365           0 :          call add_default ('DLSED    ', history_budget_histfile_num, ' ')
         366           0 :          call add_default ('HSED     ', history_budget_histfile_num, ' ')
         367           0 :          call add_default ('CMEICE   ', history_budget_histfile_num, ' ')
         368           0 :          call add_default ('LIQ2PR   ', history_budget_histfile_num, ' ')
         369           0 :          call add_default ('ICE2PR   ', history_budget_histfile_num, ' ')
         370           0 :          call add_default ('HCME     ', history_budget_histfile_num, ' ')
         371           0 :          call add_default ('HEVAP    ', history_budget_histfile_num, ' ')
         372           0 :          call add_default ('HFREEZ   ', history_budget_histfile_num, ' ')
         373           0 :          call add_default ('HMELT    ', history_budget_histfile_num, ' ')
         374           0 :          call add_default ('HREPART  ', history_budget_histfile_num, ' ')
         375           0 :          call add_default ('HPROGCLD ', history_budget_histfile_num, ' ')
         376           0 :          call add_default ('REPARTLIQ', history_budget_histfile_num, ' ')
         377           0 :          call add_default ('REPARTICE', history_budget_histfile_num, ' ')
         378             :
         379             :       end if
         380             :
         381             :    end if
         382             :
         383           0 :    if (history_amwg) then
         384           0 :       call add_default ('ICWMR', 1, ' ')
         385           0 :       call add_default ('ICIMR', 1, ' ')
         386           0 :       call add_default ('CONCLD  ', 1, ' ')
         387           0 :       call add_default ('FICE    ', 1, ' ')
         388             :    endif
         389             :
         390             :    ! History Variables for COSP/CFMIP
         391           0 :    call addfld ('LS_FLXPRC', (/ 'ilev' /), 'A', 'kg/m2/s', 'ls stratiform gbm interface rain+snow flux')
         392           0 :    call addfld ('LS_FLXSNW', (/ 'ilev' /), 'A', 'kg/m2/s', 'ls stratiform gbm interface snow flux')
         393           0 :    call addfld ('PRACWO',    (/ 'lev' /),  'A', '1/s', 'Accretion of cloud water by rain')
         394           0 :    call addfld ('PSACWO',    (/ 'lev' /),  'A', '1/s', 'Accretion of cloud water by snow')
         395           0 :    call addfld ('PSACIO',    (/ 'lev' /),  'A', '1/s', 'Accretion of cloud ice by snow')
         396             :
         397           0 :    call addfld ('CLDLIQSTR', (/ 'lev' /),  'A', 'kg/kg', 'Stratiform CLDLIQ')
         398           0 :    call addfld ('CLDICESTR', (/ 'lev' /),  'A', 'kg/kg', 'Stratiform CLDICE')
         399           0 :    call addfld ('CLDLIQCON', (/ 'lev' /),  'A', 'kg/kg', 'Convective CLDLIQ')
         400           0 :    call addfld ('CLDICECON', (/ 'lev' /),  'A', 'kg/kg', 'Convective CLDICE')
         401             :
         402           0 :    cmfmc_sh_idx = pbuf_get_index('CMFMC_SH')
         403           0 :    prec_str_idx = pbuf_get_index('PREC_STR')
         404           0 :    snow_str_idx = pbuf_get_index('SNOW_STR')
         405           0 :    prec_pcw_idx = pbuf_get_index('PREC_PCW')
         406           0 :    snow_pcw_idx = pbuf_get_index('SNOW_PCW')
         407           0 :    prec_sed_idx = pbuf_get_index('PREC_SED')
         408           0 :    snow_sed_idx = pbuf_get_index('SNOW_SED')
         409             :
         410             :    ! Initialize cldwat with constants.
         411           0 :    call inimc(tmelt, rhodair/1000.0_r8, gravit, rh2o)
         412             :
         413           0 : end subroutine rk_stratiform_init
         414             :
         415             : !===============================================================================
         416             :
         417           0 : subroutine rk_stratiform_tend( &
         418             :    state, ptend_all, pbuf, dtime, icefrac, &
         419             :    landfrac, ocnfrac, snowh, dlf,   &
         420             :    dlf2, rliq, cmfmc, ts,                  &
         421             :    sst, zdu)
         422             :
         423             :    !-------------------------------------------------------- !
         424             :    !                                                         !
         425             :    ! Interface to sedimentation, detrain, cloud fraction and !
         426             :    !        cloud macro - microphysics subroutines           !
         427             :    !                                                         !
         428             :    !-------------------------------------------------------- !
         429             :
         430           0 :    use cloud_fraction,   only: cldfrc, cldfrc_fice
         431             :    use physics_types,    only: physics_state, physics_ptend
         432             :    use physics_types,    only: physics_ptend_init, physics_update
         433             :    use physics_types,    only: physics_ptend_sum,  physics_state_copy
         434             :    use physics_types,    only: physics_state_dealloc
         435             :    use cam_history,      only: outfld
         436             :    use physics_buffer,   only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
         437             :    use pkg_cld_sediment, only: cld_sediment_vel, cld_sediment_tend
         438             :    use cldwat,           only: pcond
         439             :    use pkg_cldoptics,    only: cldefr
         440             :    use phys_control,     only: cam_physpkg_is
         441             :    use tropopause,       only: tropopause_find_cam
         442             :    use phys_grid,        only: get_rlat_all_p
         443             :    use physconst,        only: pi
         444             :
         445             :    ! Arguments
         446             :    type(physics_state), intent(in)    :: state       ! State variables
         447             :    type(physics_ptend), intent(out)   :: ptend_all   ! Package tendencies
         448             :    type(physics_buffer_desc), pointer :: pbuf(:)
         449             :
         450             :    real(r8), intent(in)  :: dtime                    ! Timestep
         451             :    real(r8), intent(in)  :: icefrac (pcols)          ! Sea ice fraction (fraction)
         452             :    real(r8), intent(in)  :: landfrac(pcols)          ! Land fraction (fraction)
         453             :    real(r8), intent(in)  :: ocnfrac (pcols)          ! Ocean fraction (fraction)
         454             :    real(r8), intent(in)  :: snowh(pcols)             ! Snow depth over land, water equivalent (m)
         455             :
         456             :    real(r8), intent(in)  :: dlf(pcols,pver)          ! Detrained water from convection schemes
         457             :    real(r8), intent(in)  :: dlf2(pcols,pver)         ! Detrained water from shallow convection scheme
         458             :    real(r8), intent(in)  :: rliq(pcols)              ! Vertical integral of liquid not yet in q(ixcldliq)
         459             :    real(r8), intent(in)  :: cmfmc(pcols,pverp)       ! Deep + Shallow Convective mass flux [ kg /s/m^2 ]
         460             :
         461             :    real(r8), intent(in)  :: ts(pcols)                ! Surface temperature
         462             :    real(r8), intent(in)  :: sst(pcols)               ! Sea surface temperature
         463             :    real(r8), intent(in)  :: zdu(pcols,pver)          ! Detrainment rate from deep convection
         464             :
         465             :   ! Local variables
         466             :
         467           0 :    type(physics_state)   :: state1                   ! Local copy of the state variable
         468           0 :    type(physics_ptend)   :: ptend_loc                ! Package tendencies
         469             :
         470             :    integer :: i, k, m
         471             :    integer :: lchnk                                  ! Chunk identifier
         472             :    integer :: ncol                                   ! Number of atmospheric columns
         473             :    integer :: itim_old
         474             :
         475             :    ! Physics buffer fields
         476           0 :    real(r8), pointer :: landm(:)             ! Land fraction ramped over water
         477             :
         478           0 :    real(r8), pointer :: prec_str(:)          ! [Total] Sfc flux of precip from stratiform [ m/s ]
         479           0 :    real(r8), pointer :: snow_str(:)          ! [Total] Sfc flux of snow from stratiform   [ m/s ]
         480           0 :    real(r8), pointer :: prec_sed(:)          ! Surface flux of total cloud water from sedimentation
         481           0 :    real(r8), pointer :: snow_sed(:)          ! Surface flux of cloud ice from sedimentation
         482           0 :    real(r8), pointer :: prec_pcw(:)          ! Sfc flux of precip from microphysics [ m/s ]
         483           0 :    real(r8), pointer :: snow_pcw(:)          ! Sfc flux of snow from microphysics [ m/s ]
         484             :
         485           0 :    real(r8), pointer, dimension(:,:) :: qcwat        ! Cloud water old q
         486           0 :    real(r8), pointer, dimension(:,:) :: tcwat        ! Cloud water old temperature
         487           0 :    real(r8), pointer, dimension(:,:) :: lcwat        ! Cloud liquid water old q
         488           0 :    real(r8), pointer, dimension(:,:) :: cld          ! Total cloud fraction
         489           0 :    real(r8), pointer, dimension(:,:) :: fice         ! Cloud ice/water partitioning ratio.
         490           0 :    real(r8), pointer, dimension(:,:) :: ast          ! Relative humidity cloud fraction
         491           0 :    real(r8), pointer, dimension(:,:) :: concld       ! Convective cloud fraction
         492           0 :    real(r8), pointer, dimension(:,:) :: qme          ! rate of cond-evap of condensate (1/s)
         493           0 :    real(r8), pointer, dimension(:,:) :: prain        ! Total precipitation (rain + snow)
         494           0 :    real(r8), pointer, dimension(:,:) :: nevapr       ! Evaporation of total precipitation (rain + snow)
         495           0 :    real(r8), pointer, dimension(:,:) :: rel          ! Liquid effective drop radius (microns)
         496           0 :    real(r8), pointer, dimension(:,:) :: rei          ! Ice effective drop size (microns)
         497           0 :    real(r8), pointer, dimension(:,:) :: wsedl        ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
         498           0 :    real(r8), pointer, dimension(:,:) :: shfrc        ! Cloud fraction from shallow convection scheme
         499           0 :    real(r8), pointer, dimension(:,:) :: cmfmc_sh     ! Shallow convective mass flux (pcols,pverp) [ kg/s/m^2 ]
         500             :
         501             :    real(r8), target :: shfrc_local(pcols,pver)
         502             :
         503             :    ! physics buffer fields for COSP simulator (RK only)
         504           0 :    real(r8), pointer, dimension(:,:) :: rkflxprc     ! RK grid-box mean flux_large_scale_cloud_rain+snow at interfaces (kg/m2/s)
         505           0 :    real(r8), pointer, dimension(:,:) :: rkflxsnw     ! RK grid-box mean flux_large_scale_cloud_snow at interfaces (kg/m2/s)
         506             :
         507             :    ! Local variables for stratiform_sediment
         508             :    real(r8) :: rain(pcols)                           ! Surface flux of cloud liquid
         509             :    real(r8) :: pvliq(pcols,pverp)                    ! Vertical velocity of cloud liquid drops (Pa/s)
         510             :    real(r8) :: pvice(pcols,pverp)                    ! Vertical velocity of cloud ice particles (Pa/s)
         511             :
         512             :    ! Local variables for cldfrc
         513             :
         514             :    real(r8) :: cldst(pcols,pver)                       ! Stratus cloud fraction
         515             :    real(r8) :: rhcloud(pcols,pver)                     ! Relative humidity cloud (last timestep)
         516             :    real(r8) :: rhcloud2(pcols,pver)                    ! Relative humidity cloud (perturbation)
         517             :    real(r8) :: clc(pcols)                              ! Column convective cloud amount
         518             :    real(r8) :: relhum(pcols,pver)                      ! RH, output to determine drh/da
         519             :    real(r8) :: rhu00(pcols,pver)
         520             :    real(r8) :: rhu002(pcols,pver)                      ! Same as rhu00 but for perturbed rh
         521             :    real(r8) :: rhdfda(pcols,pver)
         522             :    real(r8) :: cld2(pcols,pver)                        ! Same as cld but for perturbed rh
         523             :    real(r8) :: concld2(pcols,pver)                     ! Same as concld but for perturbed rh
         524             :    real(r8) :: cldst2(pcols,pver)                      ! Same as cldst but for perturbed rh
         525             :    real(r8) :: relhum2(pcols,pver)                     ! RH after  perturbation
         526             :    real(r8) :: icecldf(pcols,pver)                     ! Ice cloud fraction
         527             :    real(r8) :: liqcldf(pcols,pver)                     ! Liquid cloud fraction (combined into cloud)
         528             :    real(r8) :: icecldf_out(pcols,pver)                 ! Ice cloud fraction
         529             :    real(r8) :: liqcldf_out(pcols,pver)                 ! Liquid cloud fraction (combined into cloud)
         530             :    real(r8) :: icecldf2(pcols,pver)                    ! Ice cloud fraction
         531             :    real(r8) :: liqcldf2(pcols,pver)                    ! Liquid cloud fraction (combined into cloud)
         532             :
         533             :    ! Local variables for microphysics
         534             :
         535             :    real(r8) :: rdtime                                  ! 1./dtime
         536             :    real(r8) :: qtend(pcols,pver)                       ! Moisture tendencies
         537             :    real(r8) :: ttend(pcols,pver)                       ! Temperature tendencies
         538             :    real(r8) :: ltend(pcols,pver)                       ! Cloud liquid water tendencies
         539             :    real(r8) :: evapheat(pcols,pver)                    ! Heating rate due to evaporation of precip
         540             :    real(r8) :: evapsnow(pcols,pver)                    ! Local evaporation of snow
         541             :    real(r8) :: prfzheat(pcols,pver)                    ! Heating rate due to freezing of precip (W/kg)
         542             :    real(r8) :: meltheat(pcols,pver)                    ! Heating rate due to phase change of precip
         543             :    real(r8) :: cmeheat (pcols,pver)                    ! Heating rate due to phase change of precip
         544             :    real(r8) :: prodsnow(pcols,pver)                    ! Local production of snow
         545             :    real(r8) :: totcw(pcols,pver)                       ! Total cloud water mixing ratio
         546             :    real(r8) :: fsnow(pcols,pver)                       ! Fractional snow production
         547             :    real(r8) :: repartht(pcols,pver)                    ! Heating rate due to phase repartition of input precip
         548             :    real(r8) :: icimr(pcols,pver)                       ! In cloud ice mixing ratio
         549             :    real(r8) :: icwmr(pcols,pver)                       ! In cloud water mixing ratio
         550             :    real(r8) :: fwaut(pcols,pver)
         551             :    real(r8) :: fsaut(pcols,pver)
         552             :    real(r8) :: fracw(pcols,pver)
         553             :    real(r8) :: fsacw(pcols,pver)
         554             :    real(r8) :: fsaci(pcols,pver)
         555             :    real(r8) :: cmeice(pcols,pver)                      ! Rate of cond-evap of ice within the cloud
         556             :    real(r8) :: cmeliq(pcols,pver)                      ! Rate of cond-evap of liq within the cloud
         557             :    real(r8) :: ice2pr(pcols,pver)                      ! Rate of conversion of ice to precip
         558             :    real(r8) :: liq2pr(pcols,pver)                      ! Rate of conversion of liquid to precip
         559             :    real(r8) :: liq2snow(pcols,pver)                    ! Rate of conversion of liquid to snow
         560             :
         561             :    ! Local variables for CFMIP calculations
         562             :    real(r8) :: mr_lsliq(pcols,pver)                     ! mixing_ratio_large_scale_cloud_liquid (kg/kg)
         563             :    real(r8) :: mr_lsice(pcols,pver)                     ! mixing_ratio_large_scale_cloud_ice (kg/kg)
         564             :    real(r8) :: mr_ccliq(pcols,pver)                     ! mixing_ratio_convective_cloud_liquid (kg/kg)
         565             :    real(r8) :: mr_ccice(pcols,pver)                     ! mixing_ratio_convective_cloud_ice (kg/kg)
         566             :
         567             :    real(r8) :: pracwo(pcols,pver)                       ! RK accretion of cloud water by rain (1/s)
         568             :    real(r8) :: psacwo(pcols,pver)                       ! RK accretion of cloud water by snow (1/s)
         569             :    real(r8) :: psacio(pcols,pver)                       ! RK accretion of cloud ice by snow (1/s)
         570             :
         571             :    real(r8) :: iwc(pcols,pver)                         ! Grid box average ice water content
         572             :    real(r8) :: lwc(pcols,pver)                         ! Grid box average liquid water content
         573             :
         574             :    logical  :: lq(pcnst)
         575             :    integer  :: troplev(pcols)
         576             :    real(r8) :: rlat(pcols)
         577             :    real(r8) :: dlat(pcols)
         578             :    real(r8), parameter :: rad2deg = 180._r8/pi
         579             :
         580             :    ! ======================================================================
         581             :
         582           0 :    lchnk = state%lchnk
         583           0 :    ncol  = state%ncol
         584             :
         585           0 :    call physics_state_copy(state,state1)             ! Copy state to local state1.
         586             :
         587             :    ! Associate pointers with physics buffer fields
         588             :
         589           0 :    call pbuf_get_field(pbuf, landm_idx,   landm)
         590             :
         591           0 :    itim_old = pbuf_old_tim_idx()
         592           0 :    call pbuf_get_field(pbuf, qcwat_idx,   qcwat,   start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         593           0 :    call pbuf_get_field(pbuf, tcwat_idx,   tcwat,   start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         594           0 :    call pbuf_get_field(pbuf, lcwat_idx,   lcwat,   start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         595             :
         596           0 :    call pbuf_get_field(pbuf, cld_idx,     cld,     start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         597           0 :    call pbuf_get_field(pbuf, concld_idx,  concld,  start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         598           0 :    call pbuf_get_field(pbuf, ast_idx,     ast,     start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         599             :
         600           0 :    call pbuf_get_field(pbuf, fice_idx,    fice)
         601             :
         602           0 :    call pbuf_get_field(pbuf, cmfmc_sh_idx, cmfmc_sh)
         603             :
         604           0 :    call pbuf_get_field(pbuf, prec_str_idx, prec_str)
         605           0 :    call pbuf_get_field(pbuf, snow_str_idx, snow_str)
         606           0 :    call pbuf_get_field(pbuf, prec_sed_idx, prec_sed)
         607           0 :    call pbuf_get_field(pbuf, snow_sed_idx, snow_sed)
         608           0 :    call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw)
         609           0 :    call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw)
         610             :
         611           0 :    call pbuf_get_field(pbuf, qme_idx,    qme )
         612           0 :    call pbuf_get_field(pbuf, prain_idx,  prain)
         613           0 :    call pbuf_get_field(pbuf, nevapr_idx, nevapr)
         614             :
         615           0 :    call pbuf_get_field(pbuf, rel_idx,    rel)
         616           0 :    call pbuf_get_field(pbuf, rei_idx,    rei)
         617             :
         618           0 :    call pbuf_get_field(pbuf, wsedl_idx,  wsedl)
         619             :
         620             :    ! check that qcwat and tcwat were initialized; if not then do it now.
         621           0 :    if (qcwat(1,1) == huge(1._r8)) then
         622           0 :       qcwat(:ncol,:) = state%q(:ncol,:,1)
         623             :    end if
         624           0 :    if (tcwat(1,1) == huge(1._r8)) then
         625           0 :       tcwat(:ncol,:) = state%t(:ncol,:)
         626             :    end if
         627             :
         628           0 :    if ( do_psrhmin ) then
         629             :       !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         630           0 :       troplev(:) = 0
         631             :       !REMOVECAM_END
         632           0 :       call tropopause_find_cam(state, troplev)
         633           0 :       call get_rlat_all_p(lchnk,ncol,rlat)
         634           0 :       dlat(:ncol) = rlat(:ncol)*rad2deg
         635             :    endif
         636             :
         637             :    ! ------------- !
         638             :    ! Sedimentation !
         639             :    ! ------------- !
         640             :
         641             :    ! Allow the cloud liquid drops and ice particles to sediment.
         642             :    ! This is done before adding convectively detrained cloud water,
         643             :    ! because the phase of the detrained water is unknown.
         644             :
         645           0 :    call t_startf('stratiform_sediment')
         646             :
         647             :    call cld_sediment_vel( ncol,                                                           &
         648             :                           icefrac, landfrac, ocnfrac, state1%pmid, state1%pdel, state1%t, &
         649             :                           cld, state1%q(:,:,ixcldliq), state1%q(:,:,ixcldice),            &
         650           0 :                           pvliq, pvice, landm, snowh )
         651             :
         652           0 :    wsedl(:ncol,:pver) = pvliq(:ncol,:pver)/gravit/(state1%pmid(:ncol,:pver)/(287.15_r8*state1%t(:ncol,:pver)))
         653             :
         654           0 :    lq(:)        = .FALSE.
         655           0 :    lq(1)        = .TRUE.
         656           0 :    lq(ixcldice) = .TRUE.
         657           0 :    lq(ixcldliq) = .TRUE.
         658           0 :    call physics_ptend_init(ptend_loc, state%psetcols, 'pcwsediment', ls=.true., lq=lq)! Initialize local ptend type
         659             :
         660             :    call cld_sediment_tend( ncol, dtime ,                                                             &
         661             :                            state1%pint, state1%pmid, state1%pdel, state1%t,                          &
         662             :                            cld, state1%q(:,:,ixcldliq), state1%q(:,:,ixcldice), pvliq, pvice,        &
         663             :                            ptend_loc%q(:,:,ixcldliq), ptend_loc%q(:,:,ixcldice), ptend_loc%q(:,:,1), &
         664           0 :                            ptend_loc%s, rain, snow_sed )
         665             :
         666             :    ! Convert rain and snow fluxes at the surface from [kg/m2/s] to [m/s]
         667             :    ! Compute total precipitation flux at the surface in [m/s]
         668             :
         669           0 :    snow_sed(:ncol) = snow_sed(:ncol)/1000._r8
         670           0 :    rain(:ncol)     = rain(:ncol)/1000._r8
         671           0 :    prec_sed(:ncol) = rain(:ncol) + snow_sed(:ncol)
         672             :
         673             :    ! Record history variables
         674           0 :    call outfld( 'DQSED'   ,ptend_loc%q(:,:,1)       , pcols,lchnk )
         675           0 :    call outfld( 'DISED'   ,ptend_loc%q(:,:,ixcldice), pcols,lchnk )
         676           0 :    call outfld( 'DLSED'   ,ptend_loc%q(:,:,ixcldliq), pcols,lchnk )
         677           0 :    call outfld( 'HSED'    ,ptend_loc%s              , pcols,lchnk )
         678           0 :    call outfld( 'PRECSED' ,prec_sed                 , pcols,lchnk )
         679           0 :    call outfld( 'SNOWSED' ,snow_sed                 , pcols,lchnk )
         680           0 :    call outfld( 'RAINSED' ,rain                     , pcols,lchnk )
         681             :
         682             :    ! Add tendency from this process to tend from other processes here
         683           0 :    call physics_ptend_init(ptend_all, state%psetcols, 'stratiform')
         684           0 :    call physics_ptend_sum( ptend_loc, ptend_all, ncol )
         685             :
         686             :    ! Update physics state type state1 with ptend_loc
         687           0 :    call physics_update( state1, ptend_loc, dtime )
         688             :
         689           0 :    call t_stopf('stratiform_sediment')
         690             :
         691             :    ! Accumulate prec and snow flux at the surface [ m/s ]
         692           0 :    prec_str(:ncol) = prec_sed(:ncol)
         693           0 :    snow_str(:ncol) = snow_sed(:ncol)
         694             :
         695             :    ! ----------------------------------------------------------------------------- !
         696             :    ! Detrainment of convective condensate into the environment or stratiform cloud !
         697             :    ! ----------------------------------------------------------------------------- !
         698             :
         699             :    ! Put all of the detraining cloud water from convection into the large scale cloud.
         700             :    ! It all goes in liquid for the moment.
         701             :    ! Strictly speaking, this approach is detraining all the cconvective water into
         702             :    ! the environment, not the large-scale cloud.
         703             :
         704           0 :    lq(:)        = .FALSE.
         705           0 :    lq(ixcldliq) = .TRUE.
         706           0 :    call physics_ptend_init( ptend_loc, state1%psetcols, 'pcwdetrain', lq=lq)
         707             :
         708           0 :    do k = 1, pver
         709           0 :       do i = 1, state1%ncol
         710           0 :          ptend_loc%q(i,k,ixcldliq) = dlf(i,k)
         711             :       end do
         712             :    end do
         713             :
         714           0 :    call outfld( 'ZMDLF', dlf, pcols, lchnk )
         715           0 :    call outfld( 'SHDLF', dlf2, pcols, lchnk )
         716             :
         717             :    ! Add hie detrainment tendency to tend from the other prior processes
         718             :
         719           0 :    call physics_ptend_sum( ptend_loc, ptend_all, ncol )
         720           0 :    call physics_update( state1, ptend_loc, dtime )
         721             :
         722             :    ! Accumulate prec and snow, reserved liquid has now been used.
         723             :
         724           0 :    prec_str(:ncol) = prec_str(:ncol) - rliq(:ncol)  ! ( snow contribution is zero )
         725             :
         726             :    ! -------------------------------------- !
         727             :    ! Computation of Various Cloud Fractions !
         728             :    ! -------------------------------------- !
         729             :
         730             :    ! ----------------------------------------------------------------------------- !
         731             :    ! Treatment of cloud fraction in CAM4 and CAM5 differs                          !
         732             :    ! (1) CAM4                                                                      !
         733             :    !     . Cumulus AMT = Deep    Cumulus AMT ( empirical fcn of mass flux ) +      !
         734             :    !                     Shallow Cumulus AMT ( empirical fcn of mass flux )        !
         735             :    !     . Stratus AMT = max( RH stratus AMT, Stability Stratus AMT )              !
         736             :    !     . Cumulus and Stratus are 'minimally' overlapped without hierarchy.       !
         737             :    !     . Cumulus LWC,IWC is assumed to be the same as Stratus LWC,IWC            !
         738             :    ! (2) CAM5                                                                      !
         739             :    !     . Cumulus AMT = Deep    Cumulus AMT ( empirical fcn of mass flux ) +      !
         740             :    !                     Shallow Cumulus AMT ( internally fcn of mass flux and w ) !
         741             :    !     . Stratus AMT = fcn of environmental-mean RH ( no Stability Stratus )     !
         742             :    !     . Cumulus and Stratus are non-overlapped with higher priority on Cumulus  !
         743             :    !     . Cumulus ( both Deep and Shallow ) has its own LWC and IWC.              !
         744             :    ! ----------------------------------------------------------------------------- !
         745             :
         746           0 :    if( use_shfrc ) then
         747           0 :        call pbuf_get_field(pbuf, shfrc_idx, shfrc )
         748             :    else
         749           0 :        shfrc=>shfrc_local
         750           0 :        shfrc(:,:) = 0._r8
         751             :    endif
         752             :
         753             :    ! Stratus ('ast' = max(alst,aist)) and total cloud fraction ('cld = ast + concld')
         754             :    ! will be computed using this updated 'concld' in the stratiform macrophysics
         755             :    ! scheme (mmacro_pcond) later below.
         756             :
         757           0 :    call t_startf("cldfrc")
         758             :    call cldfrc( lchnk, ncol, pbuf,                                  &
         759             :                 state1%pmid, state1%t, state1%q(:,:,1), state1%omega, state1%phis, &
         760             :                 shfrc, use_shfrc,                                                  &
         761             :                 cld, rhcloud, clc, state1%pdel,                                    &
         762             :                 cmfmc, cmfmc_sh, landfrac,snowh, concld, cldst,                    &
         763             :                 ts, sst, state1%pint(:,pverp), zdu, ocnfrac, rhu00,                &
         764             :                 state1%q(:,:,ixcldice), icecldf, liqcldf,                          &
         765           0 :                 relhum, 0 )
         766             :
         767             :    ! Re-calculate cloud with perturbed rh add call cldfrc to estimate rhdfda.
         768             :
         769             :    call cldfrc( lchnk, ncol, pbuf,                                  &
         770             :                 state1%pmid, state1%t, state1%q(:,:,1), state1%omega, state1%phis, &
         771             :                 shfrc, use_shfrc,                                                  &
         772             :                 cld2, rhcloud2, clc, state1%pdel,                                  &
         773             :                 cmfmc, cmfmc_sh, landfrac, snowh, concld2, cldst2,                 &
         774             :                 ts, sst, state1%pint(:,pverp), zdu, ocnfrac, rhu002,               &
         775             :                 state1%q(:,:,ixcldice), icecldf2, liqcldf2,                        &
         776           0 :                 relhum2, 1 )
         777             :
         778           0 :    call t_stopf("cldfrc")
         779             :
         780           0 :    rhu00(:ncol,1) = 2.0_r8
         781           0 :    do k = 1, pver
         782           0 :       do i = 1, ncol
         783           0 :          if( relhum(i,k) < rhu00(i,k) ) then
         784           0 :             rhdfda(i,k) = 0.0_r8
         785           0 :          elseif( relhum(i,k) >= 1.0_r8 ) then
         786           0 :             rhdfda(i,k) = 0.0_r8
         787             :          else
         788             :             ! Under certain circumstances, rh+ cause cld not to changed
         789             :             ! when at an upper limit, or w/ strong subsidence
         790           0 :             if( ( cld2(i,k) - cld(i,k) ) < 1.e-4_r8 ) then
         791           0 :                rhdfda(i,k) = 0.01_r8*relhum(i,k)*1.e+4_r8
         792             :             else
         793           0 :                rhdfda(i,k) = 0.01_r8*relhum(i,k)/(cld2(i,k)-cld(i,k))
         794             :             endif
         795             :          endif
         796             :       enddo
         797             :    enddo
         798             :
         799             :    ! ---------------------------------------------- !
         800             :    ! Stratiform Cloud Macrophysics and Microphysics !
         801             :    ! ---------------------------------------------- !
         802             :
         803           0 :    call t_startf('stratiform_microphys')
         804             :
         805           0 :    rdtime = 1._r8/dtime
         806             :
         807             :    ! Define fractional amount of stratus condensate and precipitation in ice phase.
         808             :    ! This uses a ramp ( -30 ~ -10 for fice, -5 ~ 0 for fsnow ).
         809             :    ! The ramp within convective cloud may be different
         810             :
         811             : !REMOVECAM - no longer need these when CAM is retired and pcols no longer exists
         812           0 :    fice(:,:) = 0._r8
         813           0 :    fsnow(:,:) = 0._r8
         814             : !REMOVECAM_END
         815           0 :    call cldfrc_fice(ncol, state1%t(1:ncol,:), fice(1:ncol,:), fsnow(1:ncol,:))
         816             :
         817             :    ! Perform repartitioning of stratiform condensate.
         818             :    ! Corresponding heating tendency will be added later.
         819             :
         820           0 :    lq(:)        = .FALSE.
         821           0 :    lq(ixcldice) = .true.
         822           0 :    lq(ixcldliq) = .true.
         823           0 :    call physics_ptend_init( ptend_loc, state1%psetcols, 'cldwat-repartition', lq=lq )
         824             :
         825           0 :    totcw(:ncol,:pver)     = state1%q(:ncol,:pver,ixcldice) + state1%q(:ncol,:pver,ixcldliq)
         826           0 :    repartht(:ncol,:pver)  = state1%q(:ncol,:pver,ixcldice)
         827           0 :    ptend_loc%q(:ncol,:pver,ixcldice) = rdtime * ( totcw(:ncol,:pver)*fice(:ncol,:pver)          - state1%q(:ncol,:pver,ixcldice) )
         828           0 :    ptend_loc%q(:ncol,:pver,ixcldliq) = rdtime * ( totcw(:ncol,:pver)*(1.0_r8-fice(:ncol,:pver)) - state1%q(:ncol,:pver,ixcldliq) )
         829             :
         830           0 :    call outfld( 'REPARTICE', ptend_loc%q(:,:,ixcldice), pcols, lchnk )
         831           0 :    call outfld( 'REPARTLIQ', ptend_loc%q(:,:,ixcldliq), pcols, lchnk )
         832             :
         833           0 :    call physics_ptend_sum( ptend_loc, ptend_all, ncol )
         834           0 :    call physics_update( state1, ptend_loc, dtime )
         835             :
         836             :    ! Determine repartition heating from change in cloud ice.
         837             :
         838           0 :    repartht(:ncol,:pver) = (latice/dtime) * ( state1%q(:ncol,:pver,ixcldice) - repartht(:ncol,:pver) )
         839             :
         840             :    ! Non-micro and non-macrophysical external advective forcings to compute net condensation rate.
         841             :    ! Note that advective forcing of condensate is aggregated into liquid phase.
         842             :
         843           0 :    qtend(:ncol,:pver) = ( state1%q(:ncol,:pver,1) - qcwat(:ncol,:pver) ) * rdtime
         844           0 :    ttend(:ncol,:pver) = ( state1%t(:ncol,:pver)   - tcwat(:ncol,:pver) ) * rdtime
         845           0 :    ltend(:ncol,:pver) = ( totcw   (:ncol,:pver)   - lcwat(:ncol,:pver) ) * rdtime
         846             :
         847             :    ! Compute Stratiform Macro-Microphysical Tendencies
         848             :
         849             :    ! Add rain and snow fluxes as output variables from pcond, and into physics buffer
         850           0 :    call pbuf_get_field(pbuf, ls_flxprc_idx, rkflxprc)
         851           0 :    call pbuf_get_field(pbuf, ls_flxsnw_idx, rkflxsnw)
         852             :
         853           0 :    call t_startf('pcond')
         854             :    call pcond( lchnk, ncol, troplev, dlat,                                 &
         855           0 :                state1%t, ttend, state1%q(1,1,1), qtend, state1%omega,      &
         856             :                totcw, state1%pmid , state1%pdel, cld, fice, fsnow,         &
         857             :                qme, prain, prodsnow, nevapr, evapsnow, evapheat, prfzheat, &
         858             :                meltheat, prec_pcw, snow_pcw, dtime, fwaut,                 &
         859             :                fsaut, fracw, fsacw, fsaci, ltend,                          &
         860             :                rhdfda, rhu00, landm, icefrac, state1%zi, ice2pr, liq2pr,   &
         861           0 :                liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio )
         862           0 :    call t_stopf('pcond')
         863             :
         864           0 :    lq(:)        = .FALSE.
         865           0 :    lq(1)        = .true.
         866           0 :    lq(ixcldice) = .true.
         867           0 :    lq(ixcldliq) = .true.
         868           0 :    call physics_ptend_init( ptend_loc, state1%psetcols, 'cldwat', ls=.true., lq=lq)
         869             :
         870           0 :    do k = 1, pver
         871           0 :       do i = 1, ncol
         872           0 :          ptend_loc%s(i,k)          =   qme(i,k)*( latvap + latice*fice(i,k) ) + &
         873           0 :                                        evapheat(i,k) + prfzheat(i,k) + meltheat(i,k) + repartht(i,k)
         874           0 :          ptend_loc%q(i,k,1)        = - qme(i,k) + nevapr(i,k)
         875           0 :          ptend_loc%q(i,k,ixcldice) =   qme(i,k)*fice(i,k)         - ice2pr(i,k)
         876           0 :          ptend_loc%q(i,k,ixcldliq) =   qme(i,k)*(1._r8-fice(i,k)) - liq2pr(i,k)
         877             :       end do
         878             :    end do
         879             :
         880           0 :    do k = 1, pver
         881           0 :       do i = 1, ncol
         882           0 :          ast(i,k)   = cld(i,k)
         883           0 :          icimr(i,k) = (state1%q(i,k,ixcldice) + dtime*ptend_loc%q(i,k,ixcldice)) / max(0.01_r8,ast(i,k))
         884           0 :          icwmr(i,k) = (state1%q(i,k,ixcldliq) + dtime*ptend_loc%q(i,k,ixcldliq)) / max(0.01_r8,ast(i,k))
         885             :       end do
         886             :    end do
         887             :
         888             :    ! Convert precipitation from [ kg/m2 ] to [ m/s ]
         889           0 :    snow_pcw(:ncol) = snow_pcw(:ncol)/1000._r8
         890           0 :    prec_pcw(:ncol) = prec_pcw(:ncol)/1000._r8
         891             :
         892           0 :    do k = 1, pver
         893           0 :       do i = 1, ncol
         894           0 :          cmeheat(i,k) = qme(i,k) * ( latvap + latice*fice(i,k) )
         895           0 :          cmeice (i,k) = qme(i,k) *   fice(i,k)
         896           0 :          cmeliq (i,k) = qme(i,k) * ( 1._r8 - fice(i,k) )
         897             :       end do
         898             :    end do
         899             :
         900             :    ! Record history variables
         901             :
         902           0 :    call outfld( 'FWAUT'   , fwaut,       pcols, lchnk )
         903           0 :    call outfld( 'FSAUT'   , fsaut,       pcols, lchnk )
         904           0 :    call outfld( 'FRACW'   , fracw,       pcols, lchnk )
         905           0 :    call outfld( 'FSACW'   , fsacw,       pcols, lchnk )
         906           0 :    call outfld( 'FSACI'   , fsaci,       pcols, lchnk )
         907             :
         908           0 :    call outfld( 'PCSNOW'  , snow_pcw,    pcols, lchnk )
         909           0 :    call outfld( 'FICE'    , fice,        pcols, lchnk )
         910           0 :    call outfld( 'CMEICE'  , cmeice,      pcols, lchnk )
         911           0 :    call outfld( 'CMELIQ'  , cmeliq,      pcols, lchnk )
         912           0 :    call outfld( 'ICE2PR'  , ice2pr,      pcols, lchnk )
         913           0 :    call outfld( 'LIQ2PR'  , liq2pr,      pcols, lchnk )
         914           0 :    call outfld( 'HPROGCLD', ptend_loc%s, pcols, lchnk )
         915           0 :    call outfld( 'HEVAP   ', evapheat,    pcols, lchnk )
         916           0 :    call outfld( 'HMELT'   , meltheat,    pcols, lchnk )
         917           0 :    call outfld( 'HCME'    , cmeheat ,    pcols, lchnk )
         918           0 :    call outfld( 'HFREEZ'  , prfzheat,    pcols, lchnk )
         919           0 :    call outfld( 'HREPART' , repartht,    pcols, lchnk )
         920           0 :    call outfld('LS_FLXPRC', rkflxprc,    pcols, lchnk )
         921           0 :    call outfld('LS_FLXSNW', rkflxsnw,    pcols, lchnk )
         922           0 :    call outfld('PRACWO'   , pracwo,      pcols, lchnk )
         923           0 :    call outfld('PSACWO'   , psacwo,      pcols, lchnk )
         924           0 :    call outfld('PSACIO'   , psacio,      pcols, lchnk )
         925             :
         926             :    ! initialize local variables
         927           0 :    mr_ccliq(1:ncol,1:pver) = 0._r8
         928           0 :    mr_ccice(1:ncol,1:pver) = 0._r8
         929           0 :    mr_lsliq(1:ncol,1:pver) = 0._r8
         930           0 :    mr_lsice(1:ncol,1:pver) = 0._r8
         931             :
         932           0 :    do k=1,pver
         933           0 :       do i=1,ncol
         934           0 :          if (cld(i,k) .gt. 0._r8) then
         935           0 :             mr_ccliq(i,k) = (state%q(i,k,ixcldliq)/cld(i,k))*concld(i,k)
         936           0 :             mr_ccice(i,k) = (state%q(i,k,ixcldice)/cld(i,k))*concld(i,k)
         937           0 :             mr_lsliq(i,k) = (state%q(i,k,ixcldliq)/cld(i,k))*(cld(i,k)-concld(i,k))
         938           0 :             mr_lsice(i,k) = (state%q(i,k,ixcldice)/cld(i,k))*(cld(i,k)-concld(i,k))
         939             :          else
         940           0 :             mr_ccliq(i,k) = 0._r8
         941           0 :             mr_ccice(i,k) = 0._r8
         942           0 :             mr_lsliq(i,k) = 0._r8
         943           0 :             mr_lsice(i,k) = 0._r8
         944             :          end if
         945             :       end do
         946             :    end do
         947             :
         948           0 :    call outfld( 'CLDLIQSTR  ', mr_lsliq,    pcols, lchnk )
         949           0 :    call outfld( 'CLDICESTR  ', mr_lsice,    pcols, lchnk )
         950           0 :    call outfld( 'CLDLIQCON  ', mr_ccliq,    pcols, lchnk )
         951           0 :    call outfld( 'CLDICECON  ', mr_ccice,    pcols, lchnk )
         952             :
         953             :    ! ------------------------------- !
         954             :    ! Update microphysical tendencies !
         955             :    ! ------------------------------- !
         956             :
         957           0 :    call physics_ptend_sum( ptend_loc, ptend_all, ncol )
         958           0 :    call physics_update( state1, ptend_loc, dtime )
         959             :
         960           0 :    call t_startf("cldfrc")
         961             :    call cldfrc( lchnk, ncol, pbuf,                                  &
         962             :                 state1%pmid, state1%t, state1%q(:,:,1), state1%omega, state1%phis, &
         963             :                 shfrc, use_shfrc,                                                  &
         964             :                 cld, rhcloud, clc, state1%pdel,                                    &
         965             :                 cmfmc, cmfmc_sh, landfrac, snowh, concld, cldst,                   &
         966             :                 ts, sst, state1%pint(:,pverp), zdu, ocnfrac, rhu00,                &
         967             :                 state1%q(:,:,ixcldice), icecldf, liqcldf,                          &
         968           0 :                 relhum, 0 )
         969           0 :    call t_stopf("cldfrc")
         970             :
         971           0 :    call outfld( 'CONCLD  ', concld, pcols, lchnk )
         972           0 :    call outfld( 'CLDST   ', cldst,  pcols, lchnk )
         973           0 :    call outfld( 'CNVCLD  ', clc,    pcols, lchnk )
         974           0 :    call outfld( 'AST',      ast,    pcols, lchnk )
         975             :
         976           0 :    do k = 1, pver
         977           0 :       do i = 1, ncol
         978           0 :          iwc(i,k)   = state1%q(i,k,ixcldice)*state1%pmid(i,k)/(287.15_r8*state1%t(i,k))
         979           0 :          lwc(i,k)   = state1%q(i,k,ixcldliq)*state1%pmid(i,k)/(287.15_r8*state1%t(i,k))
         980           0 :          icimr(i,k) = state1%q(i,k,ixcldice) / max(0.01_r8,rhcloud(i,k))
         981           0 :          icwmr(i,k) = state1%q(i,k,ixcldliq) / max(0.01_r8,rhcloud(i,k))
         982             :       end do
         983             :    end do
         984             :
         985           0 :    call outfld( 'IWC'      , iwc,         pcols, lchnk )
         986           0 :    call outfld( 'LWC'      , lwc,         pcols, lchnk )
         987           0 :    call outfld( 'ICIMR'    , icimr,       pcols, lchnk )
         988           0 :    call outfld( 'ICWMR'    , icwmr,       pcols, lchnk )
         989           0 :    call outfld( 'CME'      , qme,         pcols, lchnk )
         990           0 :    call outfld( 'PRODPREC' , prain,       pcols, lchnk )
         991           0 :    call outfld( 'EVAPPREC' , nevapr,      pcols, lchnk )
         992           0 :    call outfld( 'EVAPSNOW' , evapsnow,    pcols, lchnk )
         993             :
         994           0 :    call t_stopf('stratiform_microphys')
         995             :
         996           0 :    prec_str(:ncol) = prec_str(:ncol) + prec_pcw(:ncol)
         997           0 :    snow_str(:ncol) = snow_str(:ncol) + snow_pcw(:ncol)
         998             :
         999             :    ! Save variables for use in the macrophysics at the next time step
        1000             :
        1001           0 :    do k = 1, pver
        1002           0 :       qcwat(:ncol,k) = state1%q(:ncol,k,1)
        1003           0 :       tcwat(:ncol,k) = state1%t(:ncol,k)
        1004           0 :       lcwat(:ncol,k) = state1%q(:ncol,k,ixcldice) + state1%q(:ncol,k,ixcldliq)
        1005             :    end do
        1006             :
        1007             :    ! Cloud water and ice particle sizes, saved in physics buffer for radiation
        1008             :
        1009           0 :    call cldefr( lchnk, ncol, landfrac, state1%t, rel, rei, state1%ps, state1%pmid, landm, icefrac, snowh )
        1010             :
        1011           0 :    call outfld('REL', rel, pcols, lchnk)
        1012           0 :    call outfld('REI', rei, pcols, lchnk)
        1013             :
        1014           0 :    call physics_state_dealloc(state1)
        1015             :
        1016           0 : end subroutine rk_stratiform_tend
        1017             :
        1018             :   !============================================================================ !
        1019             :   !                                                                             !
        1020             :   !============================================================================ !
        1021             :
        1022             :    subroutine debug_microphys_1(state1,ptend,i,k, &
        1023             :         dtime,qme,fice,snow_pcw,prec_pcw, &
        1024             :         prain,nevapr,prodsnow, evapsnow, &
        1025             :         ice2pr,liq2pr,liq2snow)
        1026             :
        1027           0 :      use physics_types, only: physics_state, physics_ptend
        1028             :      use physconst,     only: tmelt
        1029             :
        1030             :      implicit none
        1031             :
        1032             :      integer, intent(in) :: i,k
        1033             :      type(physics_state), intent(in) :: state1   ! local copy of the state variable
        1034             :      type(physics_ptend), intent(in) :: ptend  ! local copy of the ptend variable
        1035             :      real(r8), intent(in)  :: dtime                ! timestep
        1036             :      real(r8), intent(in) :: qme(pcols,pver)          ! local condensation - evaporation of cloud water
        1037             :
        1038             :      real(r8), intent(in) :: prain(pcols,pver)          ! local production of precipitation
        1039             :      real(r8), intent(in) :: nevapr(pcols,pver)          ! local evaporation of precipitation
        1040             :      real(r8), intent(in) :: prodsnow(pcols,pver)          ! local production of snow
        1041             :      real(r8), intent(in) :: evapsnow(pcols,pver)          ! local evaporation of snow
        1042             :      real(r8), intent(in) :: ice2pr(pcols,pver)   ! rate of conversion of ice to precip
        1043             :      real(r8), intent(in) :: liq2pr(pcols,pver)   ! rate of conversion of liquid to precip
        1044             :      real(r8), intent(in) :: liq2snow(pcols,pver)   ! rate of conversion of liquid to snow
        1045             :      real(r8), intent(in) :: fice    (pcols,pver)          ! Fractional ice content within cloud
        1046             :      real(r8), intent(in) :: snow_pcw(pcols)
        1047             :      real(r8), intent(in) :: prec_pcw(pcols)
        1048             :
        1049             :      real(r8) hs1, qv1, ql1, qi1, qs1, qr1, fice2, pr1, w1, w2, w3, fliq, res
        1050             :      real(r8) w4, wl, wv, wi, wlf, wvf, wif, qif, qlf, qvf
        1051             :
        1052             :      pr1 = 0
        1053             :      hs1 = 0
        1054             :      qv1 = 0
        1055             :      ql1 = 0
        1056             :      qi1 = 0
        1057             :      qs1 = 0
        1058             :      qr1 = 0
        1059             :      w1 = 0
        1060             :      wl = 0
        1061             :      wv = 0
        1062             :      wi = 0
        1063             :      wlf = 0
        1064             :      wvf = 0
        1065             :      wif = 0
        1066             :
        1067             :
        1068             :      write(iulog,*)
        1069             :      write(iulog,*) ' input state, t, q, l, i ', k, state1%t(i,k), state1%q(i,k,1), state1%q(i,k,ixcldliq),  state1%q(i,k,ixcldice)
        1070             :      write(iulog,*) ' rain, snow, total from components before accumulation ', qr1, qs1, qr1+qs1
        1071             :      write(iulog,*) ' total precip before accumulation                      ', k, pr1
        1072             :
        1073             :      wv = wv + state1%q(i,k,1       )*state1%pdel(i,k)/gravit
        1074             :      wl = wl + state1%q(i,k,ixcldliq)*state1%pdel(i,k)/gravit
        1075             :      wi = wi + state1%q(i,k,ixcldice)*state1%pdel(i,k)/gravit
        1076             :
        1077             :      qvf = state1%q(i,k,1) + ptend%q(i,k,1)*dtime
        1078             :      qlf = state1%q(i,k,ixcldliq) + ptend%q(i,k,ixcldliq)*dtime
        1079             :      qif = state1%q(i,k,ixcldice) + ptend%q(i,k,ixcldice)*dtime
        1080             :
        1081             :      if (qvf.lt.0._r8) then
        1082             :         write(iulog,*) ' qvf is negative *******', qvf
        1083             :      endif
        1084             :      if (qlf.lt.0._r8) then
        1085             :         write(iulog,*) ' qlf is negative *******', qlf
        1086             :      endif
        1087             :      if (qif.lt.0._r8) then
        1088             :         write(iulog,*) ' qif is negative *******', qif
        1089             :      endif
        1090             :      write(iulog,*) ' qvf, qlf, qif ', qvf, qlf, qif
        1091             :
        1092             :      wvf = wvf + qvf*state1%pdel(i,k)/gravit
        1093             :      wlf = wlf + qlf*state1%pdel(i,k)/gravit
        1094             :      wif = wif + qif*state1%pdel(i,k)/gravit
        1095             :
        1096             :      hs1 = hs1 + ptend%s(i,k)*state1%pdel(i,k)/gravit
        1097             :      pr1 = pr1 + state1%pdel(i,k)/gravit*(prain(i,k)-nevapr(i,k))
        1098             :      qv1 = qv1 - (qme(i,k)-nevapr(i,k))*state1%pdel(i,k)/gravit    ! vdot
        1099             :      w1  = w1  + (qme(i,k)-prain(i,k))*state1%pdel(i,k)/gravit    ! cdot
        1100             :      qi1 = qi1 + ((qme(i,k))*fice(i,k)        -ice2pr(i,k) )*state1%pdel(i,k)/gravit   ! idot
        1101             :      ql1 = ql1 + ((qme(i,k))*(1._r8-fice(i,k))-liq2pr(i,k) )*state1%pdel(i,k)/gravit   ! ldot
        1102             :
        1103             :      qr1 = qr1 &
        1104             :           + ( liq2pr(i,k)-liq2snow(i,k)   &     ! production of rain
        1105             :           -(nevapr(i,k)-evapsnow(i,k)) &     ! rain evaporation
        1106             :           )*state1%pdel(i,k)/gravit
        1107             :      qs1 = qs1 &
        1108             :           + ( ice2pr(i,k) + liq2snow(i,k) &     ! production of snow.Note last term has phase change
        1109             :           -evapsnow(i,k)               &     ! snow evaporation
        1110             :           )*state1%pdel(i,k)/gravit
        1111             :
        1112             :      if (state1%t(i,k).gt.tmelt) then
        1113             :         qr1 = qr1 + qs1
        1114             :         qs1 = 0._r8
        1115             :      endif
        1116             :      write(iulog,*) ' rain, snow, total after accumulation ', qr1, qs1, qr1+qs1
        1117             :      write(iulog,*) ' total precip after accumulation      ', k, pr1
        1118             :      write(iulog,*)
        1119             :      write(iulog,*) ' layer prain, nevapr, pdel ', prain(i,k), nevapr(i,k), state1%pdel(i,k)
        1120             :      write(iulog,*) ' layer prodsnow, ice2pr+liq2snow ', prodsnow(i,k), ice2pr(i,k)+liq2snow(i,k)
        1121             :      write(iulog,*) ' layer prain-prodsnow, liq2pr-liq2snow ', prain(i,k)-prodsnow(i,k), liq2pr(i,k)-liq2snow(i,k)
        1122             :      write(iulog,*) ' layer evapsnow, evaprain ', k, evapsnow(i,k), nevapr(i,k)-evapsnow(i,k)
        1123             :      write(iulog,*) ' layer ice2pr, liq2pr, liq2snow ', ice2pr(i,k), liq2pr(i,k), liq2snow(i,k)
        1124             :      write(iulog,*) ' layer ice2pr+liq2pr, prain ', ice2pr(i,k)+liq2pr(i,k), prain(i,k)
        1125             :      write(iulog,*)
        1126             :      write(iulog,*) ' qv1 vapor removed from col after accum  (vdot)   ', k, qv1
        1127             :      write(iulog,*) ' - (precip produced - vapor removed) after accum  ', k, -pr1-qv1
        1128             :      write(iulog,*) ' condensate produce after accum                   ', k, w1
        1129             :      write(iulog,*) ' liq+ice tends accum                              ', k, ql1+qi1
        1130             :      write(iulog,*) ' change in total water after accum                ', k, qv1+ql1+qi1
        1131             :      write(iulog,*) ' imbalance in colum after accum                   ', k, qs1+qr1+qv1+ql1+qi1
        1132             :      write(iulog,*) ' fice at this lev ', fice(i,k)
        1133             :      write(iulog,*)
        1134             :
        1135             :      res = abs((qs1+qr1+qv1+ql1+qi1)/max(abs(qv1),abs(ql1),abs(qi1),abs(qs1),abs(qr1),1.e-36_r8))
        1136             :      write(iulog,*) ' relative residual in column method 1             ', k, res
        1137             :
        1138             :      write(iulog,*) ' relative residual in column method 2             ',&
        1139             :          k, abs((qs1+qr1+qv1+ql1+qi1)/max(abs(qv1+ql1+qi1),1.e-36_r8))
        1140             :      !            if (abs((qs1+qr1+qv1+ql1+qi1)/(qs1+qr1+1.e-36)).gt.1.e-14) then
        1141             :      if (res.gt.1.e-14_r8) then
        1142             :         call endrun ('STRATIFORM_TEND')
        1143             :      endif
        1144             :
        1145             :      !             w3  = qme(i,k) * (latvap + latice*fice(i,k)) &
        1146             :      !               + evapheat(i,k) + prfzheat(i,k) + meltheat(i,k)
        1147             :
        1148             :      res = qs1+qr1-pr1
        1149             :      w4 = max(abs(qs1),abs(qr1),abs(pr1))
        1150             :      if (w4.gt.0._r8)  then
        1151             :         if (res/w4.gt.1.e-14_r8) then
        1152             :            write(iulog,*) ' imbalance in precips calculated two ways '
        1153             :            write(iulog,*) ' res/w4, pr1, qr1, qs1, qr1+qs1 ', &
        1154             :                 res/w4, pr1, qr1, qs1, qr1+qs1
        1155             :            !                   call endrun()
        1156             :         endif
        1157             :      endif
        1158             :      if (k.eq.pver) then
        1159             :         write(iulog,*) ' pcond returned precip, rain and snow rates ', prec_pcw(i), prec_pcw(i)-snow_pcw(i), snow_pcw(i)
        1160             :         write(iulog,*) ' I calculate ', pr1, qr1, qs1
        1161             :         !               call endrun
        1162             :         write(iulog,*) ' byrons water check ', wv+wl+wi-pr1*dtime, wvf+wlf+wif
        1163             :      endif
        1164             :      write(iulog,*)
        1165             :
        1166             :
        1167             :    end subroutine debug_microphys_1
        1168             :
        1169             :   !============================================================================ !
        1170             :   !                                                                             !
        1171             :   !============================================================================ !
        1172             :
        1173             :    subroutine debug_microphys_2(state1,&
        1174             :         snow_pcw,fsaut,fsacw ,fsaci, meltheat)
        1175             :
        1176             :      use ppgrid,        only: pver
        1177             :      use physconst,     only: tmelt
        1178             :      use physics_types, only: physics_state
        1179             :
        1180             :      implicit none
        1181             :
        1182             :      type(physics_state), intent(in) :: state1   ! local copy of the state variable
        1183             :      real(r8), intent(in) :: snow_pcw(pcols)
        1184             :      real(r8), intent(in) :: fsaut(pcols,pver)
        1185             :      real(r8), intent(in) :: fsacw(pcols,pver)
        1186             :      real(r8), intent(in) :: fsaci(pcols,pver)
        1187             :      real(r8), intent(in) :: meltheat(pcols,pver)          ! heating rate due to phase change of precip
        1188             :
        1189             :
        1190             :      integer  i,ncol,lchnk
        1191             :
        1192             :
        1193             :      ncol = state1%ncol
        1194             :      lchnk = state1%lchnk
        1195             :
        1196             :      do i = 1,ncol
        1197             :         if (snow_pcw(i) .gt. 0.01_r8/8.64e4_r8  .and.  state1%t(i,pver) .gt. tmelt) then
        1198             :            write(iulog,*) ' stratiform: snow, temp, ', i, lchnk, &
        1199             :                 snow_pcw(i), state1%t(i,pver)
        1200             :            write(iulog,*) ' t ', state1%t(i,:)
        1201             :            write(iulog,*) ' fsaut ', fsaut(i,:)
        1202             :            write(iulog,*) ' fsacw ', fsacw(i,:)
        1203             :            write(iulog,*) ' fsaci ', fsaci(i,:)
        1204             :            write(iulog,*) ' meltheat ', meltheat(i,:)
        1205             :            call endrun ('STRATIFORM_TEND')
        1206             :         endif
        1207             :
        1208             :         if (snow_pcw(i)*8.64e4_r8 .lt. -1.e-5_r8) then
        1209             :            write(iulog,*) ' neg snow ', snow_pcw(i)*8.64e4_r8
        1210             :            write(iulog,*) ' stratiform: snow_pcw, temp, ', i, lchnk, &
        1211             :                 snow_pcw(i), state1%t(i,pver)
        1212             :            write(iulog,*) ' t ', state1%t(i,:)
        1213             :            write(iulog,*) ' fsaut ', fsaut(i,:)
        1214             :            write(iulog,*) ' fsacw ', fsacw(i,:)
        1215             :            write(iulog,*) ' fsaci ', fsaci(i,:)
        1216             :            write(iulog,*) ' meltheat ', meltheat(i,:)
        1217             :            call endrun ('STRATIFORM_TEND')
        1218             :         endif
        1219             :      end do
        1220             :
        1221             :    end subroutine debug_microphys_2
        1222             :
        1223             :   end module rk_stratiform

Generated by: LCOV version 1.14