LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - stats_clubb_utilities.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 1143 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 10 | 0.0 %

              Line data    Source code
           1             : !-----------------------------------------------------------------------
           2             : !  $Id$
           3             : !===============================================================================
           4             : module stats_clubb_utilities
           5             :
           6             :   implicit none
           7             :
           8             :   private ! Set Default Scope
           9             :
          10             :   public :: stats_init, stats_begin_timestep, stats_end_timestep, &
          11             :     stats_accumulate, stats_finalize, stats_accumulate_hydromet, &
          12             :     stats_accumulate_lh_tend
          13             :
          14             :   private :: stats_zero, stats_avg, stats_check_num_samples
          15             :
          16             :   contains
          17             :
          18             :   !-----------------------------------------------------------------------
          19           0 :   subroutine stats_init( iunit, fname_prefix, fdir, l_stats_in, &
          20             :                          stats_fmt_in, stats_tsamp_in, stats_tout_in, fnamelist, &
          21             :                          hydromet_dim, sclr_dim, edsclr_dim, sclr_tol, &
          22           0 :                          hydromet_list, l_mix_rat_hm, &
          23           0 :                          nzmax, nlon, nlat, gzt, gzm, nnrad_zt, &
          24           0 :                          grad_zt, nnrad_zm, grad_zm, day, month, year, &
          25           0 :                          lon_vals, lat_vals, time_current, delt, l_silhs_out_in, &
          26             :                          clubb_params, &
          27             :                          l_uv_nudge, &
          28             :                          l_tke_aniso, &
          29             :                          l_standard_term_ta, &
          30             :                          stats_metadata, &
          31             :                          stats_zt, stats_zm, stats_sfc, &
          32             :                          stats_lh_zt, stats_lh_sfc, &
          33             :                          stats_rad_zt, stats_rad_zm )
          34             :     !
          35             :     ! Description:
          36             :     !   Initializes the statistics saving functionality of the CLUBB model.
          37             :     !
          38             :     ! References:
          39             :     !   None
          40             :     !-----------------------------------------------------------------------
          41             :
          42             :     use stats_variables, only: &
          43             :         stats_metadata_type
          44             :
          45             :     use parameter_indices, only: &
          46             :         nparams
          47             :
          48             :     use clubb_precision, only: &
          49             :         time_precision, & ! Constant(s)
          50             :         core_rknd
          51             :
          52             :     use output_grads, only: &
          53             :         open_grads ! Procedure
          54             :
          55             : #ifdef NETCDF
          56             :     use output_netcdf, only: &
          57             :         open_netcdf_for_writing, & ! Procedure
          58             :         first_write
          59             : #endif
          60             :
          61             :     use stats_zm_module, only: &
          62             :         nvarmax_zm, & ! Constant(s)
          63             :         stats_init_zm ! Procedure(s)
          64             :
          65             :     use stats_zt_module, only: &
          66             :         nvarmax_zt, & ! Constant(s)
          67             :         stats_init_zt ! Procedure(s)
          68             :
          69             :     use stats_lh_zt_module, only: &
          70             :         nvarmax_lh_zt, & ! Constant(s)
          71             :         stats_init_lh_zt ! Procedure(s)
          72             :
          73             :     use stats_lh_sfc_module, only: &
          74             :         nvarmax_lh_sfc, & ! Constant(s)
          75             :         stats_init_lh_sfc ! Procedure(s)
          76             :
          77             :     use stats_rad_zt_module, only: &
          78             :         nvarmax_rad_zt, & ! Constant(s)
          79             :         stats_init_rad_zt ! Procedure(s)
          80             :
          81             :     use stats_rad_zm_module, only: &
          82             :         nvarmax_rad_zm, & ! Constant(s)
          83             :         stats_init_rad_zm ! Procedure(s)
          84             :
          85             :     use stats_sfc_module, only: &
          86             :         nvarmax_sfc, & ! Constant(s)
          87             :         stats_init_sfc ! Procedure(s)
          88             :
          89             :     use constants_clubb, only: &
          90             :         fstdout, fstderr, var_length ! Constants
          91             :
          92             :     use error_code, only: &
          93             :         clubb_at_least_debug_level, &   ! Procedure
          94             :         err_code, &                     ! Error Indicator
          95             :         clubb_fatal_error               ! Constant
          96             :
          97             :     use stats_type, only: stats ! Type
          98             :
          99             :     implicit none
         100             :
         101             :     ! Local Constants
         102             :     integer, parameter :: &
         103             :       silhs_num_importance_categories = 8
         104             :
         105             :     ! Input Variables
         106             :     integer, intent(in) :: iunit  ! File unit for fnamelist
         107             :
         108             :     character(len=*), intent(in) ::  &
         109             :       fname_prefix, & ! Start of the stats filenames
         110             :       fdir            ! Directory to output to
         111             :
         112             :     logical, intent(in) :: &
         113             :       l_stats_in      ! Stats on? T/F
         114             :
         115             :     character(len=*), intent(in) :: &
         116             :       stats_fmt_in    ! Format of the stats file output
         117             :
         118             :     real( kind = core_rknd ), intent(in) ::  &
         119             :       stats_tsamp_in,  & ! Sampling interval   [s]
         120             :       stats_tout_in      ! Output interval     [s]
         121             :
         122             :     character(len=*), intent(in) :: &
         123             :       fnamelist          ! Filename holding the &statsnl
         124             :
         125             :     integer, intent(in) :: &
         126             :       nlon, & ! Number of points in the X direction [-]
         127             :       nlat, & ! Number of points in the Y direction [-]
         128             :       nzmax   ! Grid points in the vertical         [-]
         129             :
         130             :     real( kind = core_rknd ), intent(in), dimension(nzmax) ::  &
         131             :       gzt, gzm  ! Thermodynamic and momentum levels           [m]
         132             :
         133             :     integer, intent(in) :: &
         134             :       nnrad_zt,     & ! Grid points in the radiation grid [count]
         135             :       hydromet_dim, &
         136             :       sclr_dim,     &
         137             :       edsclr_dim
         138             :
         139             :     real( kind = core_rknd ), dimension(sclr_dim), intent(in) :: &
         140             :       sclr_tol
         141             :
         142             :     character(len=10), dimension(hydromet_dim), intent(in) :: &
         143             :       hydromet_list
         144             :
         145             :     logical, dimension(hydromet_dim), intent(in) :: &
         146             :       l_mix_rat_hm   ! if true, then the quantity is a hydrometeor mixing ratio
         147             :
         148             :     real( kind = core_rknd ), intent(in), dimension(nnrad_zt) :: grad_zt ! Radiation levels [m]
         149             :
         150             :     integer, intent(in) :: nnrad_zm ! Grid points in the radiation grid [count]
         151             :
         152             :     real( kind = core_rknd ), intent(in), dimension(nnrad_zm) :: grad_zm ! Radiation levels [m]
         153             :
         154             :     integer, intent(in) :: day, month, year  ! Time of year
         155             :
         156             :     real( kind = core_rknd ), dimension(nlon), intent(in) ::  &
         157             :       lon_vals  ! Longitude values [Degrees E]
         158             :
         159             :     real( kind = core_rknd ), dimension(nlat), intent(in) ::  &
         160             :       lat_vals  ! Latitude values  [Degrees N]
         161             :
         162             :     real( kind = time_precision ), intent(in) ::  &
         163             :       time_current ! Model time                         [s]
         164             :
         165             :     real( kind = core_rknd ), intent(in) ::  &
         166             :       delt         ! Timestep (dt_main in CLUBB)         [s]
         167             :
         168             :     logical, intent(in) :: &
         169             :       l_silhs_out_in  ! Whether to output SILHS files (stats_lh_zt, stats_lh_sfc)  [boolean]
         170             :
         171             :     real( kind = core_rknd ), dimension(nparams), intent(in) :: &
         172             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
         173             :
         174             :     logical, intent(in) :: &
         175             :       l_uv_nudge,         & ! For wind speed nudging
         176             :       l_tke_aniso,        & ! For anisotropic turbulent kinetic energy, i.e. TKE = 1/2
         177             :                             ! (u'^2 + v'^2 + w'^2)
         178             :       l_standard_term_ta    ! Use the standard discretization for the turbulent advection terms.
         179             :                             ! Setting to .false. means that a_1 and a_3 are pulled outside of the
         180             :                             ! derivative in advance_wp2_wp3_module.F90 and in
         181             :                             ! advance_xp2_xpyp_module.F90.
         182             :
         183             :     type (stats_metadata_type), intent(inout) :: &
         184             :       stats_metadata
         185             :
         186             :     type (stats), target, intent(inout) :: &
         187             :       stats_zt, &
         188             :       stats_zm, &
         189             :       stats_sfc, &
         190             :       stats_lh_zt, &
         191             :       stats_lh_sfc, &
         192             :       stats_rad_zt, &
         193             :       stats_rad_zm
         194             :
         195             :     ! Local Variables
         196             :     logical :: l_error
         197             :
         198             :     character(len=200) :: fname
         199             :
         200             :     integer :: ivar, ntot, read_status
         201             :
         202             :     ! Namelist Variables
         203             :
         204             :     character(len=10) :: stats_fmt  ! File storage convention
         205             :
         206             :     character(len=var_length), dimension(nvarmax_zt) ::  &
         207             :       vars_zt  ! Variables on the thermodynamic levels
         208             :
         209             :     character(len=var_length), dimension(nvarmax_lh_zt) ::  &
         210             :       vars_lh_zt  ! Latin Hypercube variables on the thermodynamic levels
         211             :
         212             :     character(len=var_length), dimension(nvarmax_lh_sfc) ::  &
         213             :       vars_lh_sfc  ! Latin Hypercube variables at the surface
         214             :
         215             :     character(len=var_length), dimension(nvarmax_zm) ::  &
         216             :       vars_zm  ! Variables on the momentum levels
         217             :
         218             :     character(len=var_length), dimension(nvarmax_rad_zt) ::  &
         219             :       vars_rad_zt  ! Variables on the radiation levels
         220             :
         221             :     character(len=var_length), dimension(nvarmax_rad_zm) ::  &
         222             :       vars_rad_zm  ! Variables on the radiation levels
         223             :
         224             :     character(len=var_length), dimension(nvarmax_sfc) ::  &
         225             :       vars_sfc ! Variables at the model surface
         226             :
         227             :     namelist /statsnl/ &
         228             :       vars_zt, &
         229             :       vars_zm, &
         230             :       vars_lh_zt, &
         231             :       vars_lh_sfc, &
         232             :       vars_rad_zt, &
         233             :       vars_rad_zm, &
         234             :       vars_sfc
         235             :
         236             :     ! ---- Begin Code ----
         237             :
         238             :     ! Initialize
         239           0 :     l_error = .false.
         240             :
         241             :     ! Set stats_variables variables with inputs from calling subroutine
         242           0 :     stats_metadata%l_stats = l_stats_in
         243             :
         244           0 :     stats_metadata%stats_tsamp = stats_tsamp_in
         245           0 :     stats_metadata%stats_tout  = stats_tout_in
         246           0 :     stats_fmt   = trim( stats_fmt_in )
         247           0 :     stats_metadata%l_silhs_out = l_silhs_out_in
         248             :
         249           0 :     if ( .not. stats_metadata%l_stats ) then
         250           0 :       stats_metadata%l_stats_samp  = .false.
         251           0 :       stats_metadata%l_stats_last  = .false.
         252           0 :       return
         253             :     end if
         254             :
         255             :     ! Initialize namelist variables
         256             :
         257           0 :     vars_zt  = ''
         258           0 :     vars_zm  = ''
         259           0 :     vars_lh_zt = ''
         260           0 :     vars_lh_sfc = ''
         261           0 :     vars_rad_zt = ''
         262           0 :     vars_rad_zm = ''
         263           0 :     vars_sfc = ''
         264             :
         265             :     ! Reads list of variables that should be output to GrADS/NetCDF (namelist &statsnl)
         266             :
         267           0 :     open(unit=iunit, file=fnamelist)
         268           0 :     read(unit=iunit, nml=statsnl, iostat=read_status, end=100)
         269           0 :     if ( read_status /= 0 ) then
         270           0 :       if ( read_status > 0 ) then
         271           0 :         write(fstderr,*) "Error reading stats namelist in file ",  &
         272           0 :                          trim( fnamelist )
         273             :       else ! Read status < 0
         274           0 :         write(fstderr,*) "End of file marker reached while reading stats namelist in file ", &
         275           0 :           trim( fnamelist )
         276             :       end if
         277           0 :       write(fstderr,*) "One cause is having more statistical variables ",  &
         278           0 :                        "listed in the namelist for var_zt, var_zm, or ",  &
         279           0 :                        "var_sfc than allowed by nvarmax_zt, nvarmax_zm, ",  &
         280           0 :                        "or nvarmax_sfc, respectively."
         281           0 :       write(fstderr,*) "Maximum variables allowed for var_zt = ", nvarmax_zt
         282           0 :       write(fstderr,*) "Maximum variables allowed for var_zm = ", nvarmax_zm
         283           0 :       write(fstderr,*) "Maximum variables allowed for var_rad_zt = ", nvarmax_rad_zt
         284           0 :       write(fstderr,*) "Maximum variables allowed for var_rad_zm = ", nvarmax_rad_zm
         285           0 :       write(fstderr,*) "Maximum variables allowed for var_sfc = ", nvarmax_sfc
         286           0 :       write(fstderr,*) "stats_init: Error reading stats namelist."
         287           0 :       err_code = clubb_fatal_error
         288           0 :       close(unit=iunit)
         289           0 :       return
         290             :     end if ! read_status /= 0
         291             :
         292           0 :     close(unit=iunit)
         293             :
         294           0 :     if ( clubb_at_least_debug_level( 1 ) ) then
         295           0 :       write(fstdout,*) "--------------------------------------------------"
         296             :
         297           0 :       write(fstdout,*) "Statistics"
         298             :
         299           0 :       write(fstdout,*) "--------------------------------------------------"
         300           0 :       write(fstdout,*) "vars_zt = "
         301           0 :       ivar = 1
         302           0 :       do while ( vars_zt(ivar) /= '' )
         303           0 :         write(fstdout,*) vars_zt(ivar)
         304           0 :         ivar = ivar + 1
         305             :       end do
         306             :
         307           0 :       write(fstdout,*) "vars_zm = "
         308           0 :       ivar = 1
         309           0 :       do while ( vars_zm(ivar) /= '' )
         310           0 :         write(fstdout,*) vars_zm(ivar)
         311           0 :         ivar = ivar + 1
         312             :       end do
         313             :
         314           0 :       if ( stats_metadata%l_silhs_out ) then
         315           0 :         write(fstdout,*) "vars_lh_zt = "
         316           0 :         ivar = 1
         317           0 :         do while ( vars_lh_zt(ivar) /= '' )
         318           0 :           write(fstdout,*) vars_lh_zt(ivar)
         319           0 :           ivar = ivar + 1
         320             :         end do
         321             :
         322           0 :         write(fstdout,*) "vars_lh_sfc = "
         323           0 :         ivar = 1
         324           0 :         do while ( vars_lh_sfc(ivar) /= '' )
         325           0 :           write(fstdout,*) vars_lh_sfc(ivar)
         326           0 :           ivar = ivar + 1
         327             :         end do
         328             :       end if ! l_silhs_out
         329             :
         330           0 :       if ( stats_metadata%l_output_rad_files ) then
         331           0 :         write(fstdout,*) "vars_rad_zt = "
         332           0 :         ivar = 1
         333           0 :         do while ( vars_rad_zt(ivar) /= '' )
         334           0 :           write(fstdout,*) vars_rad_zt(ivar)
         335           0 :           ivar = ivar + 1
         336             :         end do
         337             :
         338           0 :         write(fstdout,*) "vars_rad_zm = "
         339           0 :         ivar = 1
         340           0 :         do while ( vars_rad_zm(ivar) /= '' )
         341           0 :           write(fstdout,*) vars_rad_zm(ivar)
         342           0 :           ivar = ivar + 1
         343             :         end do
         344             :       end if ! l_output_rad_files
         345             :
         346           0 :       write(fstdout,*) "vars_sfc = "
         347           0 :       ivar = 1
         348           0 :       do while ( vars_sfc(ivar) /= '' )
         349           0 :         write(fstdout,*) vars_sfc(ivar)
         350           0 :         ivar = ivar + 1
         351             :       end do
         352             :
         353           0 :       write(fstdout,*) "--------------------------------------------------"
         354             :     end if ! clubb_at_least_debug_level 1
         355             :
         356             :     ! Determine file names for GrADS or NetCDF files
         357           0 :     stats_metadata%fname_zt  = trim( fname_prefix )//"_zt"
         358           0 :     stats_metadata%fname_zm  = trim( fname_prefix )//"_zm"
         359           0 :     stats_metadata%fname_lh_zt  = trim( fname_prefix )//"_lh_zt"
         360           0 :     stats_metadata%fname_lh_sfc  = trim( fname_prefix )//"_lh_sfc"
         361           0 :     stats_metadata%fname_rad_zt  = trim( fname_prefix )//"_rad_zt"
         362           0 :     stats_metadata%fname_rad_zm  = trim( fname_prefix )//"_rad_zm"
         363           0 :     stats_metadata%fname_sfc = trim( fname_prefix )//"_sfc"
         364             :
         365             :     ! Parse the file type for stats output.  Currently only GrADS and
         366             :     ! netCDF > version 3.5 are supported by this code.
         367           0 :     select case ( trim( stats_fmt ) )
         368             :     case ( "GrADS", "grads", "gr" )
         369           0 :       stats_metadata%l_netcdf = .false.
         370           0 :       stats_metadata%l_grads  = .true.
         371             :
         372             :     case ( "NetCDF", "netcdf", "nc" )
         373           0 :       stats_metadata%l_netcdf = .true.
         374           0 :       stats_metadata%l_grads  = .false.
         375             :
         376             :     case default
         377           0 :       write(fstderr,*) "In module stats_clubb_utilities subroutine stats_init: "
         378           0 :       write(fstderr,*) "Invalid stats output format "//trim( stats_fmt )
         379           0 :       err_code = clubb_fatal_error
         380           0 :       return
         381             :
         382             :     end select
         383             :
         384             :     ! Check sampling and output frequencies
         385             :
         386             :     ! The model time step length, delt (which is dt_main), should multiply
         387             :     ! evenly into the statistical sampling time step length, stats_tsamp.
         388           0 :     if ( abs( stats_metadata%stats_tsamp/delt &
         389             :               - real( floor( stats_metadata%stats_tsamp/delt ), kind=core_rknd ) )  &
         390             :            > 1.e-8_core_rknd) then
         391           0 :       l_error = .true.  ! This will cause the run to stop.
         392           0 :       write(fstderr,*) 'Error:  stats_tsamp should be an even multiple of ',  &
         393           0 :                        'delt (which is dt_main).  Check the appropriate ',  &
         394           0 :                        'model.in file.'
         395           0 :       write(fstderr,*) 'stats_tsamp = ', stats_metadata%stats_tsamp
         396           0 :       write(fstderr,*) 'delt = ', delt
         397             :     end if
         398             :
         399             :     ! The statistical sampling time step length, stats_tsamp, should multiply
         400             :     ! evenly into the statistical output time step length, stats_tout.
         401           0 :     if ( abs( stats_metadata%stats_tout/stats_metadata%stats_tsamp &
         402             :            - real( floor( stats_metadata%stats_tout    &
         403             :                           / stats_metadata%stats_tsamp ), kind=core_rknd) ) &
         404             :          > 1.e-8_core_rknd) then
         405           0 :       l_error = .true.  ! This will cause the run to stop.
         406           0 :       write(fstderr,*) 'Error:  stats_tout should be an even multiple of ',  &
         407           0 :                        'stats_tsamp.  Check the appropriate model.in file.'
         408           0 :       write(fstderr,*) 'stats_tout = ', stats_metadata%stats_tout
         409           0 :       write(fstderr,*) 'stats_tsamp = ', stats_metadata%stats_tsamp
         410             :     end if
         411             :
         412             :     ! Initialize zt (mass points)
         413             :
         414           0 :     ivar = 1
         415           0 :     do while ( ichar(vars_zt(ivar)(1:1)) /= 0  &
         416             :                .and. len_trim(vars_zt(ivar)) /= 0 &
         417           0 :                .and. ivar <= nvarmax_zt )
         418           0 :       ivar = ivar + 1
         419             :     end do
         420           0 :     ntot = ivar - 1
         421             :
         422           0 :     if ( any( vars_zt == "hm_i" ) ) then
         423             :        ! Correct for number of variables found under "hm_i".
         424             :        ! Subtract "hm_i" from the number of zt statistical variables.
         425           0 :        ntot = ntot - 1
         426             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         427             :        ! to the number of zt statistical variables.
         428           0 :        ntot = ntot + 2 * hydromet_dim
         429             :     endif
         430           0 :     if ( any( vars_zt == "mu_hm_i" ) ) then
         431             :        ! Correct for number of variables found under "mu_hm_i".
         432             :        ! Subtract "mu_hm_i" from the number of zt statistical variables.
         433           0 :        ntot = ntot - 1
         434             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         435             :        ! to the number of zt statistical variables.
         436           0 :        ntot = ntot + 2 * hydromet_dim
         437             :     endif
         438           0 :     if ( any( vars_zt == "mu_Ncn_i" ) ) then
         439             :        ! Correct for number of variables found under "mu_Ncn_i".
         440             :        ! Subtract "mu_Ncn_i" from the number of zt statistical variables.
         441           0 :        ntot = ntot - 1
         442             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         443             :        ! statistical variables.
         444           0 :        ntot = ntot + 2
         445             :     endif
         446           0 :     if ( any( vars_zt == "mu_hm_i_n" ) ) then
         447             :        ! Correct for number of variables found under "mu_hm_i_n".
         448             :        ! Subtract "mu_hm_i_n" from the number of zt statistical variables.
         449           0 :        ntot = ntot - 1
         450             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         451             :        ! to the number of zt statistical variables.
         452           0 :        ntot = ntot + 2 * hydromet_dim
         453             :     endif
         454           0 :     if ( any( vars_zt == "mu_Ncn_i_n" ) ) then
         455             :        ! Correct for number of variables found under "mu_Ncn_i_n".
         456             :        ! Subtract "mu_Ncn_i_n" from the number of zt statistical variables.
         457           0 :        ntot = ntot - 1
         458             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         459             :        ! statistical variables.
         460           0 :        ntot = ntot + 2
         461             :     endif
         462           0 :     if ( any( vars_zt == "sigma_hm_i" ) ) then
         463             :        ! Correct for number of variables found under "sigma_hm_i".
         464             :        ! Subtract "sigma_hm_i" from the number of zt statistical variables.
         465           0 :        ntot = ntot - 1
         466             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         467             :        ! to the number of zt statistical variables.
         468           0 :        ntot = ntot + 2 * hydromet_dim
         469             :     endif
         470           0 :     if ( any( vars_zt == "sigma_Ncn_i" ) ) then
         471             :        ! Correct for number of variables found under "sigma_Ncn_i".
         472             :        ! Subtract "sigma_Ncn_i" from the number of zt statistical variables.
         473           0 :        ntot = ntot - 1
         474             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         475             :        ! statistical variables.
         476           0 :        ntot = ntot + 2
         477             :     endif
         478           0 :     if ( any( vars_zt == "sigma_hm_i_n" ) ) then
         479             :        ! Correct for number of variables found under "sigma_hm_i_n".
         480             :        ! Subtract "sigma_hm_i_n" from the number of zt statistical variables.
         481           0 :        ntot = ntot - 1
         482             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         483             :        ! to the number of zt statistical variables.
         484           0 :        ntot = ntot + 2 * hydromet_dim
         485             :     endif
         486           0 :     if ( any( vars_zt == "sigma_Ncn_i_n" ) ) then
         487             :        ! Correct for number of variables found under "sigma_Ncn_i_n".
         488             :        ! Subtract "sigma_Ncn_i_n" from the number of zt statistical variables.
         489           0 :        ntot = ntot - 1
         490             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         491             :        ! statistical variables.
         492           0 :        ntot = ntot + 2
         493             :     endif
         494             :
         495           0 :     if ( any( vars_zt == "corr_w_hm_i" ) ) then
         496             :        ! Correct for number of variables found under "corr_w_hm_i".
         497             :        ! Subtract "corr_w_hm_i" from the number of zt statistical variables.
         498           0 :        ntot = ntot - 1
         499             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         500             :        ! to the number of zt statistical variables.
         501           0 :        ntot = ntot + 2 * hydromet_dim
         502             :     endif
         503           0 :     if ( any( vars_zt == "corr_w_Ncn_i" ) ) then
         504             :        ! Correct for number of variables found under "corr_w_Ncn_i".
         505             :        ! Subtract "corr_w_Ncn_i" from the number of zt statistical variables.
         506           0 :        ntot = ntot - 1
         507             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         508             :        ! statistical variables.
         509           0 :        ntot = ntot + 2
         510             :     endif
         511           0 :     if ( any( vars_zt == "corr_chi_hm_i" ) ) then
         512             :        ! Correct for number of variables found under "corr_chi_hm_i".
         513             :        ! Subtract "corr_chi_hm_i" from the number of zt statistical variables.
         514           0 :        ntot = ntot - 1
         515             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         516             :        ! to the number of zt statistical variables.
         517           0 :        ntot = ntot + 2 * hydromet_dim
         518             :     endif
         519           0 :     if ( any( vars_zt == "corr_chi_Ncn_i" ) ) then
         520             :        ! Correct for number of variables found under "corr_chi_Ncn_i".
         521             :        ! Subtract "corr_chi_Ncn_i" from the number of zt statistical variables.
         522           0 :        ntot = ntot - 1
         523             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         524             :        ! statistical variables.
         525           0 :        ntot = ntot + 2
         526             :     endif
         527           0 :     if ( any( vars_zt == "corr_eta_hm_i" ) ) then
         528             :        ! Correct for number of variables found under "corr_eta_hm_i".
         529             :        ! Subtract "corr_eta_hm_i" from the number of zt statistical variables.
         530           0 :        ntot = ntot - 1
         531             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         532             :        ! to the number of zt statistical variables.
         533           0 :        ntot = ntot + 2 * hydromet_dim
         534             :     endif
         535           0 :     if ( any( vars_zt == "corr_eta_Ncn_i" ) ) then
         536             :        ! Correct for number of variables found under "corr_eta_Ncn_i".
         537             :        ! Subtract "corr_eta_Ncn_i" from the number of zt statistical variables.
         538           0 :        ntot = ntot - 1
         539             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         540             :        ! statistical variables.
         541           0 :        ntot = ntot + 2
         542             :     endif
         543           0 :     if ( any( vars_zt == "corr_Ncn_hm_i" ) ) then
         544             :        ! Correct for number of variables found under "corr_Ncn_hm_i".
         545             :        ! Subtract "corr_Ncn_hm_i" from the number of zt statistical variables.
         546           0 :        ntot = ntot - 1
         547             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         548             :        ! to the number of zt statistical variables.
         549           0 :        ntot = ntot + 2 * hydromet_dim
         550             :     endif
         551           0 :     if ( any( vars_zt == "corr_hmx_hmy_i" ) ) then
         552             :        ! Correct for number of variables found under "corr_hmx_hmy_i".
         553             :        ! Subtract "corr_hmx_hmy_i" from the number of zt statistical variables.
         554           0 :        ntot = ntot - 1
         555             :        ! Add 2 (1st PDF component and 2nd PDF component) multipled by the
         556             :        ! number of correlations of two hydrometeors, which is found by:
         557             :        ! (1/2) * hydromet_dim * ( hydromet_dim - 1 );
         558             :        ! to the number of zt statistical variables.
         559           0 :        ntot = ntot + hydromet_dim * ( hydromet_dim - 1 )
         560             :     endif
         561             :
         562           0 :     if ( any( vars_zt == "corr_w_hm_i_n" ) ) then
         563             :        ! Correct for number of variables found under "corr_w_hm_i_n".
         564             :        ! Subtract "corr_w_hm_i_n" from the number of zt statistical variables.
         565           0 :        ntot = ntot - 1
         566             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         567             :        ! to the number of zt statistical variables.
         568           0 :        ntot = ntot + 2 * hydromet_dim
         569             :     endif
         570           0 :     if ( any( vars_zt == "corr_w_Ncn_i_n" ) ) then
         571             :        ! Correct for number of variables found under "corr_w_Ncn_i_n".
         572             :        ! Subtract "corr_w_Ncn_i_n" from the number of zt statistical variables.
         573           0 :        ntot = ntot - 1
         574             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         575             :        ! statistical variables.
         576           0 :        ntot = ntot + 2
         577             :     endif
         578           0 :     if ( any( vars_zt == "corr_chi_hm_i_n" ) ) then
         579             :        ! Correct for number of variables found under "corr_chi_hm_i_n".
         580             :        ! Subtract "corr_chi_hm_i_n" from the number of zt statistical variables.
         581           0 :        ntot = ntot - 1
         582             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         583             :        ! to the number of zt statistical variables.
         584           0 :        ntot = ntot + 2 * hydromet_dim
         585             :     endif
         586           0 :     if ( any( vars_zt == "corr_chi_Ncn_i_n" ) ) then
         587             :        ! Correct for number of variables found under "corr_chi_Ncn_i_n".
         588             :        ! Subtract "corr_chi_Ncn_i_n" from the number of zt statistical variables.
         589           0 :        ntot = ntot - 1
         590             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         591             :        ! statistical variables.
         592           0 :        ntot = ntot + 2
         593             :     endif
         594           0 :     if ( any( vars_zt == "corr_eta_hm_i_n" ) ) then
         595             :        ! Correct for number of variables found under "corr_eta_hm_i_n".
         596             :        ! Subtract "corr_eta_hm_i_n" from the number of zt statistical variables.
         597           0 :        ntot = ntot - 1
         598             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         599             :        ! to the number of zt statistical variables.
         600           0 :        ntot = ntot + 2 * hydromet_dim
         601             :     endif
         602           0 :     if ( any( vars_zt == "corr_eta_Ncn_i_n" ) ) then
         603             :        ! Correct for number of variables found under "corr_eta_Ncn_i_n".
         604             :        ! Subtract "corr_eta_Ncn_i_n" from the number of zt statistical variables.
         605           0 :        ntot = ntot - 1
         606             :        ! Add 2 (1st PDF component and 2nd PDF component) to the number of zt
         607             :        ! statistical variables.
         608           0 :        ntot = ntot + 2
         609             :     endif
         610           0 :     if ( any( vars_zt == "corr_Ncn_hm_i_n" ) ) then
         611             :        ! Correct for number of variables found under "corr_Ncn_hm_i_n".
         612             :        ! Subtract "corr_Ncn_hm_i_n" from the number of zt statistical variables.
         613           0 :        ntot = ntot - 1
         614             :        ! Add 2 (1st PDF component and 2nd PDF component) for each hydrometeor
         615             :        ! to the number of zt statistical variables.
         616           0 :        ntot = ntot + 2 * hydromet_dim
         617             :     endif
         618           0 :     if ( any( vars_zt == "corr_hmx_hmy_i_n" ) ) then
         619             :        ! Correct for number of variables found under "corr_hmx_hmy_i_n".
         620             :        ! Subtract "corr_hmx_hmy_i_n" from the number of zt statistical variables.
         621           0 :        ntot = ntot - 1
         622             :        ! Add 2 (1st PDF component and 2nd PDF component) multipled by the
         623             :        ! number of normal space correlations of two hydrometeors, which is
         624             :        ! found by:  (1/2) * hydromet_dim * ( hydromet_dim - 1 );
         625             :        ! to the number of zt statistical variables.
         626           0 :        ntot = ntot + hydromet_dim * ( hydromet_dim - 1 )
         627             :     endif
         628             :
         629           0 :     if ( any( vars_zt == "hmp2_zt" ) ) then
         630             :        ! Correct for number of variables found under "hmp2_zt".
         631             :        ! Subtract "hmp2_zt" from the number of zt statistical variables.
         632           0 :        ntot = ntot - 1
         633             :        ! Add 1 for each hydrometeor to the number of zt statistical variables.
         634           0 :        ntot = ntot + hydromet_dim
         635             :     endif
         636             :
         637           0 :     if ( any( vars_zt == "wp2hmp" ) ) then
         638             :        ! Correct for number of variables found under "wp2hmp".
         639             :        ! Subtract "wp2hmp" from the number of zt statistical variables.
         640           0 :        ntot = ntot - 1
         641             :        ! Add 1 for each hydrometeor to the number of zt statistical variables.
         642           0 :        ntot = ntot + hydromet_dim
         643             :     endif
         644             :
         645           0 :     if ( any( vars_zt == "sclrm" ) ) then
         646             :        ! Correct for number of variables found under "sclrm".
         647             :        ! Subtract "sclrm" from the number of zt statistical variables.
         648           0 :        ntot = ntot - 1
         649             :        ! Add 1 for each scalar to the number of zt statistical variables.
         650           0 :        ntot = ntot + sclr_dim
         651             :     endif
         652             :
         653           0 :     if ( any( vars_zt == "sclrm_f" ) ) then
         654             :        ! Correct for number of variables found under "sclrm_f".
         655             :        ! Subtract "sclrm_f" from the number of zt statistical variables.
         656           0 :        ntot = ntot - 1
         657             :        ! Add 1 for each scalar to the number of zt statistical variables.
         658           0 :        ntot = ntot + sclr_dim
         659             :     endif
         660             :
         661           0 :     if ( any( vars_zt == "edsclrm" ) ) then
         662             :        ! Correct for number of variables found under "edsclrm".
         663             :        ! Subtract "edsclrm" from the number of zt statistical variables.
         664           0 :        ntot = ntot - 1
         665             :        ! Add 1 for each scalar to the number of zt statistical variables.
         666           0 :        ntot = ntot + edsclr_dim
         667             :     endif
         668             :
         669           0 :     if ( any( vars_zt == "edsclrm_f" ) ) then
         670             :        ! Correct for number of variables found under "edsclrm_f".
         671             :        ! Subtract "edsclrm_f" from the number of zt statistical variables.
         672           0 :        ntot = ntot - 1
         673             :        ! Add 1 for each scalar to the number of zt statistical variables.
         674           0 :        ntot = ntot + edsclr_dim
         675             :     endif
         676             :
         677           0 :     if ( ntot >= nvarmax_zt ) then
         678           0 :       write(fstderr,*) "There are more statistical variables listed in ",  &
         679           0 :                        "vars_zt than allowed for by nvarmax_zt."
         680           0 :       write(fstderr,*) "Check the number of variables listed for vars_zt ",  &
         681           0 :                        "in the stats namelist, or change nvarmax_zt."
         682           0 :       write(fstderr,*) "nvarmax_zt = ", nvarmax_zt
         683           0 :       write(fstderr,*) "number of variables in vars_zt = ", ntot
         684           0 :       write(fstderr,*) "stats_init:  number of zt statistical variables exceeds limit"
         685           0 :       err_code = clubb_fatal_error
         686           0 :       return
         687             :     end if
         688             :
         689           0 :     stats_zt%num_output_fields = ntot
         690           0 :     stats_zt%kk = nzmax
         691           0 :     stats_zt%ii = nlon
         692           0 :     stats_zt%jj = nlat
         693             :
         694           0 :     allocate( stats_zt%z( stats_zt%kk ) )
         695           0 :     stats_zt%z = gzt
         696             :
         697           0 :     allocate( stats_zt%accum_field_values( stats_zt%ii, stats_zt%jj, &
         698           0 :       stats_zt%kk, stats_zt%num_output_fields ) )
         699           0 :     allocate( stats_zt%accum_num_samples( stats_zt%ii, stats_zt%jj, &
         700           0 :       stats_zt%kk, stats_zt%num_output_fields ) )
         701           0 :     allocate( stats_zt%l_in_update( stats_zt%ii, stats_zt%jj, stats_zt%kk, &
         702           0 :       stats_zt%num_output_fields ) )
         703             :     call stats_zero( stats_zt%ii, stats_zt%jj, stats_zt%kk, stats_zt%num_output_fields, & ! In
         704           0 :       stats_zt%accum_field_values, stats_zt%accum_num_samples, stats_zt%l_in_update ) ! Out
         705             :
         706           0 :     allocate( stats_zt%file%grid_avg_var( stats_zt%num_output_fields ) )
         707           0 :     allocate( stats_zt%file%z( stats_zt%kk ) )
         708             :
         709           0 :     fname = trim( stats_metadata%fname_zt )
         710             :
         711             :     ! Default initialization for array indices for zt
         712             :
         713             :     call stats_init_zt( hydromet_dim, sclr_dim, edsclr_dim, & ! intent(in)
         714             :                         hydromet_list, l_mix_rat_hm,        & ! intent(in)
         715             :                         vars_zt,                            & ! intent(in)
         716             :                         l_error,                            & ! intent(inout)
         717           0 :                         stats_metadata, stats_zt )            ! intent(inout)
         718             :
         719           0 :     if ( stats_metadata%l_grads ) then
         720             :
         721             :       ! Open GrADS file
         722             :       call open_grads( iunit, fdir, fname,  &  ! intent(in)
         723             :                        1, stats_zt%kk, nlat, nlon, stats_zt%z, & ! intent(in)
         724             :                        day, month, year, lat_vals, lon_vals, &  ! intent(in)
         725             :                        time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
         726             :                        stats_metadata, stats_zt%num_output_fields, & ! intent(in)
         727           0 :                        stats_zt%file ) ! intent(inout)
         728             :
         729             :     else ! Open NetCDF file
         730             : #ifdef NETCDF
         731             :       call open_netcdf_for_writing( nlat, nlon, fdir, fname, 1, stats_zt%kk, stats_zt%z, &  ! In
         732             :                         day, month, year, lat_vals, lon_vals, &  ! In
         733             :                         time_current, stats_metadata%stats_tout, stats_zt%num_output_fields, &  ! In
         734             :                         stats_zt%file ) ! InOut
         735             :
         736             :       ! Finalize the variable definitions
         737             :       call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
         738             :                         l_uv_nudge, & ! intent(in)
         739             :                         l_tke_aniso, & ! intent(in)
         740             :                         l_standard_term_ta, & ! intent(in)
         741             :                         stats_zt%file ) ! intent(inout)
         742             :
         743             :       if ( err_code == clubb_fatal_error ) then
         744             :         write(fstderr,*) "Fatal error setting up stats_zt"
         745             :         return
         746             :       end if
         747             : #else
         748           0 :       error stop "This CLUBB program was not compiled with netCDF support."
         749             : #endif
         750             :
         751             :     end if
         752             :
         753             :
         754             :     ! Setup output file for stats_lh_zt (Latin Hypercube stats)
         755             :
         756           0 :     if ( stats_metadata%l_silhs_out ) then
         757             :
         758             :       ivar = 1
         759           0 :       do while ( ichar(vars_lh_zt(ivar)(1:1)) /= 0  &
         760             :                  .and. len_trim(vars_lh_zt(ivar)) /= 0 &
         761           0 :                  .and. ivar <= nvarmax_lh_zt )
         762           0 :         ivar = ivar + 1
         763             :       end do
         764           0 :       ntot = ivar - 1
         765           0 :       if ( any( vars_lh_zt == "silhs_variance_category" ) ) then
         766             :         ! Correct for number of variables found under "silhs_variance_category".
         767             :         ! Subtract "silhs_variance_category" from the number of lh_zt statistical
         768             :         ! variables.
         769           0 :         ntot = ntot - 1
         770             :         ! Add 1 for each SILHS category to the number of lh_zt statistical variables
         771           0 :         ntot = ntot + silhs_num_importance_categories
         772             :       end if
         773             :
         774           0 :       if ( any( vars_lh_zt == "lh_samp_frac_category" ) ) then
         775             :         ! Correct for number of variables found under "lh_samp_frac_category".
         776             :         ! Subtract "lh_samp_frac_category" from the number of lh_zt statistical
         777             :         ! variables.
         778           0 :         ntot = ntot - 1
         779             :         ! Add 1 for each SILHS category to the number of lh_zt statistical variables
         780           0 :         ntot = ntot + silhs_num_importance_categories
         781             :       end if
         782             :
         783           0 :       if ( ntot == nvarmax_lh_zt ) then
         784           0 :         write(fstderr,*) "There are more statistical variables listed in ",  &
         785           0 :                          "vars_lh_zt than allowed for by nvarmax_lh_zt."
         786           0 :         write(fstderr,*) "Check the number of variables listed for vars_lh_zt ",  &
         787           0 :                          "in the stats namelist, or change nvarmax_lh_zt."
         788           0 :         write(fstderr,*) "nvarmax_lh_zt = ", nvarmax_lh_zt
         789           0 :         write(fstderr,*) "number of variables in vars_lh_zt = ", ntot
         790           0 :         write(fstderr,*) "stats_init:  number of lh_zt statistical variables exceeds limit"
         791           0 :         err_code = clubb_fatal_error
         792           0 :         return
         793             :       end if
         794             :
         795           0 :       stats_lh_zt%num_output_fields = ntot
         796           0 :       stats_lh_zt%kk = nzmax
         797           0 :       stats_lh_zt%ii = nlon
         798           0 :       stats_lh_zt%jj = nlat
         799             :
         800           0 :       allocate( stats_lh_zt%z( stats_lh_zt%kk ) )
         801           0 :       stats_lh_zt%z = gzt
         802             :
         803           0 :       allocate( stats_lh_zt%accum_field_values( stats_lh_zt%ii, stats_lh_zt%jj, &
         804           0 :         stats_lh_zt%kk, stats_lh_zt%num_output_fields ) )
         805           0 :       allocate( stats_lh_zt%accum_num_samples( stats_lh_zt%ii, stats_lh_zt%jj, &
         806           0 :         stats_lh_zt%kk, stats_lh_zt%num_output_fields ) )
         807           0 :       allocate( stats_lh_zt%l_in_update( stats_lh_zt%ii, stats_lh_zt%jj, stats_lh_zt%kk, &
         808           0 :         stats_lh_zt%num_output_fields ) )
         809             :       call stats_zero( stats_lh_zt%ii, stats_lh_zt%jj, stats_lh_zt%kk, & ! intent(in)
         810             :         stats_lh_zt%num_output_fields, & ! intent(in)
         811             :         stats_lh_zt%accum_field_values, stats_lh_zt%accum_num_samples, & ! intent(out)
         812           0 :         stats_lh_zt%l_in_update ) ! intent(out)
         813             :
         814           0 :       allocate( stats_lh_zt%file%grid_avg_var( stats_lh_zt%num_output_fields ) )
         815           0 :       allocate( stats_lh_zt%file%z( stats_lh_zt%kk ) )
         816             :
         817             :       call stats_init_lh_zt( vars_lh_zt,                    & ! intent(in)
         818             :                              l_error,                       & ! intent(inout)
         819           0 :                              stats_metadata, stats_lh_zt )    ! intent(inout)
         820             :
         821             :
         822           0 :       fname = trim( stats_metadata%fname_lh_zt )
         823             :
         824           0 :       if ( stats_metadata%l_grads ) then
         825             :
         826             :         ! Open GrADS file
         827             :         call open_grads( iunit, fdir, fname,  & ! intent(in)
         828             :                          1, stats_lh_zt%kk, nlat, nlon, stats_lh_zt%z, & ! intent(in)
         829             :                          day, month, year, lat_vals, lon_vals, &  ! intent(in)
         830             :                          time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
         831             :                          stats_metadata, stats_lh_zt%num_output_fields, & ! intent(in)
         832           0 :                          stats_lh_zt%file ) ! intent(inout)
         833             :
         834             :       else ! Open NetCDF file
         835             : #ifdef NETCDF
         836             :         call open_netcdf_for_writing( nlat, nlon, fdir, fname, 1, stats_lh_zt%kk, &  ! In
         837             :                           stats_lh_zt%z, day, month, year, lat_vals, lon_vals, &  ! In
         838             :                           time_current, stats_metadata%stats_tout, stats_lh_zt%num_output_fields, &  ! In
         839             :                           stats_lh_zt%file ) ! InOut
         840             :
         841             :         ! Finalize the variable definitions
         842             :         call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
         843             :                           l_uv_nudge, & ! intent(in)
         844             :                           l_tke_aniso, & ! intent(in)
         845             :                           l_standard_term_ta, & ! intent(in)
         846             :                           stats_lh_zt%file ) ! intent(inout)
         847             :
         848             :         if ( err_code == clubb_fatal_error ) then
         849             :           write(fstderr,*) "Fatal error setting up stats_lh_zt"
         850             :           return
         851             :         end if
         852             : #else
         853           0 :         error stop "This CLUBB program was not compiled with netCDF support."
         854             : #endif
         855             :
         856             :       end if
         857             :
         858           0 :       ivar = 1
         859           0 :       do while ( ichar(vars_lh_sfc(ivar)(1:1)) /= 0  &
         860             :                  .and. len_trim(vars_lh_sfc(ivar)) /= 0 &
         861           0 :                  .and. ivar <= nvarmax_lh_sfc )
         862           0 :         ivar = ivar + 1
         863             :       end do
         864           0 :       ntot = ivar - 1
         865           0 :       if ( ntot == nvarmax_lh_sfc ) then
         866           0 :         write(fstderr,*) "There are more statistical variables listed in ",  &
         867           0 :                          "vars_lh_sfc than allowed for by nvarmax_lh_sfc."
         868           0 :         write(fstderr,*) "Check the number of variables listed for vars_lh_sfc ",  &
         869           0 :                          "in the stats namelist, or change nvarmax_lh_sfc."
         870           0 :         write(fstderr,*) "nvarmax_lh_sfc = ", nvarmax_lh_sfc
         871           0 :         write(fstderr,*) "number of variables in vars_lh_sfc = ", ntot
         872           0 :         write(fstderr,*) "stats_init:  number of lh_sfc statistical variables exceeds limit"
         873           0 :         err_code = clubb_fatal_error
         874           0 :         return
         875             :       end if
         876             :
         877           0 :       stats_lh_sfc%num_output_fields = ntot
         878           0 :       stats_lh_sfc%kk = 1
         879           0 :       stats_lh_sfc%ii = nlon
         880           0 :       stats_lh_sfc%jj = nlat
         881             :
         882           0 :       allocate( stats_lh_sfc%z( stats_lh_sfc%kk ) )
         883           0 :       stats_lh_sfc%z = gzm(1)
         884             :
         885           0 :       allocate( stats_lh_sfc%accum_field_values( stats_lh_sfc%ii, stats_lh_sfc%jj, &
         886           0 :         stats_lh_sfc%kk, stats_lh_sfc%num_output_fields ) )
         887           0 :       allocate( stats_lh_sfc%accum_num_samples( stats_lh_sfc%ii, stats_lh_sfc%jj, &
         888           0 :         stats_lh_sfc%kk, stats_lh_sfc%num_output_fields ) )
         889           0 :       allocate( stats_lh_sfc%l_in_update( stats_lh_sfc%ii, stats_lh_sfc%jj, &
         890           0 :         stats_lh_sfc%kk, stats_lh_sfc%num_output_fields ) )
         891             :
         892             :       call stats_zero( stats_lh_sfc%ii, stats_lh_sfc%jj, stats_lh_sfc%kk, & ! intent(in)
         893             :           stats_lh_sfc%num_output_fields, & ! intent(in)
         894             :           stats_lh_sfc%accum_field_values, & ! intent(out)
         895           0 :           stats_lh_sfc%accum_num_samples, stats_lh_sfc%l_in_update ) ! intent(out)
         896             :
         897           0 :       allocate( stats_lh_sfc%file%grid_avg_var( stats_lh_sfc%num_output_fields ) )
         898           0 :       allocate( stats_lh_sfc%file%z( stats_lh_sfc%kk ) )
         899             :
         900             :       call stats_init_lh_sfc( vars_lh_sfc,                    & ! intent(in)
         901             :                               l_error,                        & ! intent(inout)
         902           0 :                               stats_metadata, stats_lh_sfc )    ! intent(inout)
         903             :
         904           0 :       fname = trim( stats_metadata%fname_lh_sfc )
         905             :
         906           0 :       if ( stats_metadata%l_grads ) then
         907             :
         908             :         ! Open GrADS file
         909             :         call open_grads( iunit, fdir, fname,  & ! intent(in)
         910             :                          1, stats_lh_sfc%kk, nlat, nlon, stats_lh_sfc%z, & ! intent(in)
         911             :                          day, month, year, lat_vals, lon_vals, &  ! intent(in)
         912             :                          time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
         913             :                          stats_metadata, stats_lh_sfc%num_output_fields, & ! intent(in)
         914           0 :                          stats_lh_sfc%file ) ! intent(inout)
         915             :
         916             :       else ! Open NetCDF file
         917             : #ifdef NETCDF
         918             :         call open_netcdf_for_writing( nlat, nlon, fdir, fname, 1, stats_lh_sfc%kk, &  ! In
         919             :                           stats_lh_sfc%z, day, month, year, lat_vals, lon_vals, &  ! In
         920             :                           time_current, stats_metadata%stats_tout, stats_lh_sfc%num_output_fields, &  ! In
         921             :                           stats_lh_sfc%file ) ! InOut
         922             :
         923             :         ! Finalize the variable definitions
         924             :         call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
         925             :                           l_uv_nudge, & ! intent(in)
         926             :                           l_tke_aniso, & ! intent(in)
         927             :                           l_standard_term_ta, & ! intent(in)
         928             :                           stats_lh_sfc%file ) ! intent(inout)
         929             :
         930             :         if ( err_code == clubb_fatal_error ) then
         931             :           write(fstderr,*) "Fatal error setting up stats_lh_sfc"
         932             :           return
         933             :         end if
         934             :
         935             :         if ( err_code == clubb_fatal_error ) return
         936             : #else
         937           0 :         error stop "This CLUBB program was not compiled with netCDF support."
         938             : #endif
         939             :
         940             :       end if
         941             :
         942             :     end if ! l_silhs_out
         943             :
         944             :     ! Initialize stats_zm (momentum points)
         945             :
         946           0 :     ivar = 1
         947           0 :     do while ( ichar(vars_zm(ivar)(1:1)) /= 0  &
         948             :                .and. len_trim(vars_zm(ivar)) /= 0 &
         949           0 :                .and. ivar <= nvarmax_zm )
         950           0 :       ivar = ivar + 1
         951             :     end do
         952           0 :     ntot = ivar - 1
         953             :
         954           0 :     if ( any( vars_zm == "hydrometp2" ) ) then
         955             :        ! Correct for number of variables found under "hydrometp2".
         956             :        ! Subtract "hydrometp2" from the number of zm statistical variables.
         957           0 :        ntot = ntot - 1
         958             :        ! Add 1 for each hydrometeor to the number of zm statistical variables.
         959           0 :        ntot = ntot + hydromet_dim
         960             :     endif
         961             :
         962           0 :     if ( any( vars_zm == "wphydrometp" ) ) then
         963             :        ! Correct for number of variables found under "wphydrometp".
         964             :        ! Subtract "wphydrometp" from the number of zm statistical variables.
         965           0 :        ntot = ntot - 1
         966             :        ! Add 1 for each hydrometeor to the number of zm statistical variables.
         967           0 :        ntot = ntot + hydromet_dim
         968             :     endif
         969             :
         970           0 :     if ( any( vars_zm == "rtphmp" ) ) then
         971             :        ! Correct for number of variables found under "rtphmp".
         972             :        ! Subtract "rtphmp" from the number of zm statistical variables.
         973           0 :        ntot = ntot - 1
         974             :        ! Add 1 for each hydrometeor to the number of zm statistical variables.
         975           0 :        ntot = ntot + hydromet_dim
         976             :     endif
         977             :
         978           0 :     if ( any( vars_zm == "thlphmp" ) ) then
         979             :        ! Correct for number of variables found under "thlphmp".
         980             :        ! Subtract "thlphmp" from the number of zm statistical variables.
         981           0 :        ntot = ntot - 1
         982             :        ! Add 1 for each hydrometeor to the number of zm statistical variables.
         983           0 :        ntot = ntot + hydromet_dim
         984             :     endif
         985             :
         986           0 :     if ( any( vars_zm == "hmxphmyp" ) ) then
         987             :        ! Correct for number of variables found under "hmxphmyp".
         988             :        ! Subtract "hmxphmyp" from the number of zm statistical variables.
         989           0 :        ntot = ntot - 1
         990             :        ! Add the number of overall covariances of two hydrometeors, which is
         991             :        ! found by:  (1/2) * hydromet_dim * ( hydromet_dim - 1 );
         992             :        ! to the number of zm statistical variables.
         993           0 :        ntot = ntot + hydromet_dim * ( hydromet_dim - 1 ) / 2
         994             :     endif
         995             :
         996           0 :     if ( any( vars_zm == "K_hm" ) ) then
         997             :        ! Correct for number of variables found under "K_hm".
         998             :        ! Subtract "K_hm" from the number of zm statistical variables.
         999           0 :        ntot = ntot - 1
        1000             :        ! Add 1 for each hydrometeor to the number of zm statistical variables.
        1001           0 :        ntot = ntot + hydromet_dim
        1002             :     endif
        1003             :
        1004           0 :     if ( any( vars_zm == "sclrprtp" ) ) then
        1005             :        ! Correct for number of variables found under "sclrprtp".
        1006             :        ! Subtract "sclrprtp" from the number of zm statistical variables.
        1007           0 :        ntot = ntot - 1
        1008             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1009           0 :        ntot = ntot + sclr_dim
        1010             :     endif
        1011             :
        1012           0 :     if ( any( vars_zm == "sclrp2" ) ) then
        1013             :        ! Correct for number of variables found under "sclrp2".
        1014             :        ! Subtract "sclrp2" from the number of zm statistical variables.
        1015           0 :        ntot = ntot - 1
        1016             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1017           0 :        ntot = ntot + sclr_dim
        1018             :     endif
        1019             :
        1020             :
        1021           0 :     if ( any( vars_zm == "sclrpthvp" ) ) then
        1022             :        ! Correct for number of variables found under "sclrpthvp".
        1023             :        ! Subtract "sclrpthvp" from the number of zm statistical variables.
        1024           0 :        ntot = ntot - 1
        1025             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1026           0 :        ntot = ntot + sclr_dim
        1027             :     endif
        1028             :
        1029             :
        1030           0 :     if ( any( vars_zm == "sclrpthlp" ) ) then
        1031             :        ! Correct for number of variables found under "sclrpthlp".
        1032             :        ! Subtract "sclrpthlp" from the number of zm statistical variables.
        1033           0 :        ntot = ntot - 1
        1034             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1035           0 :        ntot = ntot + sclr_dim
        1036             :     endif
        1037             :
        1038             :
        1039           0 :     if ( any( vars_zm == "sclrprcp" ) ) then
        1040             :        ! Correct for number of variables found under "sclrprcp".
        1041             :        ! Subtract "sclrprcp" from the number of zm statistical variables.
        1042           0 :        ntot = ntot - 1
        1043             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1044           0 :        ntot = ntot + sclr_dim
        1045             :     endif
        1046             :
        1047             :
        1048           0 :     if ( any( vars_zm == "wpsclrp" ) ) then
        1049             :        ! Correct for number of variables found under "wpsclrp".
        1050             :        ! Subtract "wpsclrp" from the number of zm statistical variables.
        1051           0 :        ntot = ntot - 1
        1052             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1053           0 :        ntot = ntot + sclr_dim
        1054             :     endif
        1055             :
        1056             :
        1057           0 :     if ( any( vars_zm == "wpsclrp2" ) ) then
        1058             :        ! Correct for number of variables found under "wpsclrp2".
        1059             :        ! Subtract "wpsclrp2" from the number of zm statistical variables.
        1060           0 :        ntot = ntot - 1
        1061             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1062           0 :        ntot = ntot + sclr_dim
        1063             :     endif
        1064             :
        1065             :
        1066           0 :     if ( any( vars_zm == "wp2sclrp" ) ) then
        1067             :        ! Correct for number of variables found under "wp2sclrp".
        1068             :        ! Subtract "wp2sclrp" from the number of zm statistical variables.
        1069           0 :        ntot = ntot - 1
        1070             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1071           0 :        ntot = ntot + sclr_dim
        1072             :     endif
        1073             :
        1074             :
        1075           0 :     if ( any( vars_zm == "wpsclrprtp" ) ) then
        1076             :        ! Correct for number of variables found under "wpsclrprtp".
        1077             :        ! Subtract "wpsclrprtp" from the number of zm statistical variables.
        1078           0 :        ntot = ntot - 1
        1079             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1080           0 :        ntot = ntot + sclr_dim
        1081             :     endif
        1082             :
        1083             :
        1084           0 :     if ( any( vars_zm == "wpsclrpthlp" ) ) then
        1085             :        ! Correct for number of variables found under "wpsclrpthlp".
        1086             :        ! Subtract "wpsclrpthlp" from the number of zm statistical variables.
        1087           0 :        ntot = ntot - 1
        1088             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1089           0 :        ntot = ntot + sclr_dim
        1090             :     endif
        1091             :
        1092             :
        1093           0 :     if ( any( vars_zm == "wpedsclrp" ) ) then
        1094             :        ! Correct for number of variables found under "wpedsclrp".
        1095             :        ! Subtract "wpedsclrp" from the number of zm statistical variables.
        1096           0 :        ntot = ntot - 1
        1097             :        ! Add 1 for each scalar to the number of zm statistical variables.
        1098           0 :        ntot = ntot + edsclr_dim
        1099             :     endif
        1100             :
        1101             :
        1102             :
        1103           0 :     if ( ntot == nvarmax_zm ) then
        1104           0 :       write(fstderr,*) "There are more statistical variables listed in ",  &
        1105           0 :                        "vars_zm than allowed for by nvarmax_zm."
        1106           0 :       write(fstderr,*) "Check the number of variables listed for vars_zm ",  &
        1107           0 :                        "in the stats namelist, or change nvarmax_zm."
        1108           0 :       write(fstderr,*) "nvarmax_zm = ", nvarmax_zm
        1109           0 :       write(fstderr,*) "number of variables in vars_zm = ", ntot
        1110           0 :       write(fstderr,*) "stats_init:  number of zm statistical variables exceeds limit"
        1111           0 :       err_code = clubb_fatal_error
        1112           0 :       return
        1113             :     end if
        1114             :
        1115           0 :     stats_zm%num_output_fields = ntot
        1116           0 :     stats_zm%kk = nzmax
        1117           0 :     stats_zm%ii = nlon
        1118           0 :     stats_zm%jj = nlat
        1119             :
        1120           0 :     allocate( stats_zm%z( stats_zm%kk ) )
        1121           0 :     stats_zm%z = gzm
        1122             :
        1123           0 :     allocate( stats_zm%accum_field_values( stats_zm%ii, stats_zm%jj, &
        1124           0 :       stats_zm%kk, stats_zm%num_output_fields ) )
        1125           0 :     allocate( stats_zm%accum_num_samples( stats_zm%ii, stats_zm%jj, &
        1126           0 :       stats_zm%kk, stats_zm%num_output_fields ) )
        1127           0 :     allocate( stats_zm%l_in_update( stats_zm%ii, stats_zm%jj, stats_zm%kk, &
        1128           0 :       stats_zm%num_output_fields ) )
        1129             :
        1130             :     call stats_zero( stats_zm%ii, stats_zm%jj, stats_zm%kk, stats_zm%num_output_fields, & ! In
        1131           0 :       stats_zm%accum_field_values, stats_zm%accum_num_samples, stats_zm%l_in_update ) ! intent(out)
        1132             :
        1133           0 :     allocate( stats_zm%file%grid_avg_var( stats_zm%num_output_fields ) )
        1134           0 :     allocate( stats_zm%file%z( stats_zm%kk ) )
        1135             :
        1136             :     call stats_init_zm( hydromet_dim, sclr_dim, edsclr_dim, & ! intent(in)
        1137             :                         hydromet_list, l_mix_rat_hm,        & ! intent(in)
        1138             :                         vars_zm,                            & ! intent(in)
        1139             :                         l_error,                            & ! intent(inout)
        1140           0 :                         stats_metadata, stats_zm )            ! intent(inout)
        1141             :
        1142           0 :     fname = trim( stats_metadata%fname_zm )
        1143           0 :     if ( stats_metadata%l_grads ) then
        1144             :
        1145             :       ! Open GrADS files
        1146             :       call open_grads( iunit, fdir, fname,  & ! intent(in)
        1147             :                        1, stats_zm%kk, nlat, nlon, stats_zm%z, & ! intent(in)
        1148             :                        day, month, year, lat_vals, lon_vals, & ! intent(in)
        1149             :                        time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
        1150             :                        stats_metadata, stats_zm%num_output_fields, & ! intent(in)
        1151           0 :                        stats_zm%file ) ! intent(inout)
        1152             :
        1153             :     else ! Open NetCDF file
        1154             : #ifdef NETCDF
        1155             :       call open_netcdf_for_writing( nlat, nlon, fdir, fname, 1, stats_zm%kk, stats_zm%z, &  ! In
        1156             :                         day, month, year, lat_vals, lon_vals, &  ! In
        1157             :                         time_current, stats_metadata%stats_tout, stats_zm%num_output_fields, &  ! In
        1158             :                         stats_zm%file ) ! InOut
        1159             :
        1160             :       ! Finalize the variable definitions
        1161             :       call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
        1162             :                         l_uv_nudge, & ! intent(in)
        1163             :                         l_tke_aniso, & ! intent(in)
        1164             :                         l_standard_term_ta, & ! intent(in)
        1165             :                         stats_zm%file ) ! intent(inout)
        1166             :
        1167             :       if ( err_code == clubb_fatal_error ) then
        1168             :         write(fstderr,*) "Fatal error setting up stats_zm"
        1169             :         return
        1170             :       end if
        1171             : #else
        1172           0 :       error stop "This CLUBB program was not compiled with netCDF support."
        1173             : #endif
        1174             :     end if
        1175             :
        1176             :     ! Initialize stats_rad_zt (radiation points)
        1177             :
        1178           0 :     if ( stats_metadata%l_output_rad_files ) then
        1179             :
        1180             :       ivar = 1
        1181           0 :       do while ( ichar(vars_rad_zt(ivar)(1:1)) /= 0  &
        1182             :                  .and. len_trim(vars_rad_zt(ivar)) /= 0 &
        1183           0 :                  .and. ivar <= nvarmax_rad_zt )
        1184           0 :         ivar = ivar + 1
        1185             :       end do
        1186           0 :       ntot = ivar - 1
        1187           0 :       if ( ntot == nvarmax_rad_zt ) then
        1188           0 :         write(fstderr,*) "There are more statistical variables listed in ",  &
        1189           0 :                          "vars_rad_zt than allowed for by nvarmax_rad_zt."
        1190           0 :         write(fstderr,*) "Check the number of variables listed for vars_rad_zt ",  &
        1191           0 :                          "in the stats namelist, or change nvarmax_rad_zt."
        1192           0 :         write(fstderr,*) "nvarmax_rad_zt = ", nvarmax_rad_zt
        1193           0 :         write(fstderr,*) "number of variables in vars_rad_zt = ", ntot
        1194           0 :         write(fstderr,*) "stats_init:  number of rad_zt statistical variables exceeds limit"
        1195           0 :         err_code = clubb_fatal_error
        1196           0 :         return
        1197             :       end if
        1198             :
        1199           0 :       stats_rad_zt%num_output_fields = ntot
        1200           0 :       stats_rad_zt%kk = nnrad_zt
        1201           0 :       stats_rad_zt%ii = nlon
        1202           0 :       stats_rad_zt%jj = nlat
        1203           0 :       allocate( stats_rad_zt%z( stats_rad_zt%kk ) )
        1204           0 :       stats_rad_zt%z = grad_zt
        1205             :
        1206           0 :       allocate( stats_rad_zt%accum_field_values( stats_rad_zt%ii, stats_rad_zt%jj, &
        1207           0 :         stats_rad_zt%kk, stats_rad_zt%num_output_fields ) )
        1208           0 :       allocate( stats_rad_zt%accum_num_samples( stats_rad_zt%ii, stats_rad_zt%jj, &
        1209           0 :         stats_rad_zt%kk, stats_rad_zt%num_output_fields ) )
        1210           0 :       allocate( stats_rad_zt%l_in_update( stats_rad_zt%ii, stats_rad_zt%jj, &
        1211           0 :         stats_rad_zt%kk, stats_rad_zt%num_output_fields ) )
        1212             :
        1213             :       call stats_zero( stats_rad_zt%ii, stats_rad_zt%jj, stats_rad_zt%kk, & ! intent(in)
        1214             :                        stats_rad_zt%num_output_fields, & ! intent(in)
        1215             :                        stats_rad_zt%accum_field_values, & ! intent(out)
        1216           0 :                        stats_rad_zt%accum_num_samples, stats_rad_zt%l_in_update )! intent(out)
        1217             :
        1218           0 :       allocate( stats_rad_zt%file%grid_avg_var( stats_rad_zt%num_output_fields ) )
        1219           0 :       allocate( stats_rad_zt%file%z( stats_rad_zt%kk ) )
        1220             :
        1221             :       call stats_init_rad_zt( vars_rad_zt,                    & ! intent(in)
        1222             :                               l_error,                        & ! intent(inout)
        1223           0 :                               stats_metadata, stats_rad_zt )    ! intent(inout)
        1224             :
        1225           0 :       fname = trim( stats_metadata%fname_rad_zt )
        1226           0 :       if ( stats_metadata%l_grads ) then
        1227             :
        1228             :         ! Open GrADS files
        1229             :         call open_grads( iunit, fdir, fname,  & ! intent(in)
        1230             :                          1, stats_rad_zt%kk, nlat, nlon, stats_rad_zt%z, & ! intent(in)
        1231             :                          day, month, year, lat_vals, lon_vals, &
        1232             :                          time_current+real(stats_metadata%stats_tout, kind=time_precision), & !intent(in)
        1233             :                          stats_metadata, stats_rad_zt%num_output_fields, & ! intent(in)
        1234           0 :                          stats_rad_zt%file ) ! intent(inout)
        1235             :
        1236             :       else ! Open NetCDF file
        1237             : #ifdef NETCDF
        1238             :         call open_netcdf_for_writing( nlat, nlon, fdir, fname,  & ! intent(in)
        1239             :                           1, stats_rad_zt%kk, stats_rad_zt%z, & ! intent(in)
        1240             :                           day, month, year, lat_vals, lon_vals, & ! intent(in)
        1241             :                           time_current, stats_metadata%stats_tout, & ! intent(in)
        1242             :                           stats_rad_zt%num_output_fields, & ! intent(in)
        1243             :                           stats_rad_zt%file ) ! intent(inout)
        1244             :
        1245             :         ! Finalize the variable definitions
        1246             :         call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
        1247             :                           l_uv_nudge, & ! intent(in)
        1248             :                           l_tke_aniso, & ! intent(in)
        1249             :                           l_standard_term_ta, & ! intent(in)
        1250             :                           stats_rad_zt%file ) ! intent(inout)
        1251             :
        1252             :         if ( err_code == clubb_fatal_error ) then
        1253             :           write(fstderr,*) "Fatal error setting up stats_rad_zt"
        1254             :           return
        1255             :         end if
        1256             : #else
        1257           0 :         error stop "This CLUBB program was not compiled with netCDF support."
        1258             : #endif
        1259             :       end if
        1260             :
        1261             :       ! Initialize stats_rad_zm (radiation points)
        1262             :
        1263           0 :       ivar = 1
        1264           0 :       do while ( ichar(vars_rad_zm(ivar)(1:1)) /= 0  &
        1265             :                  .and. len_trim(vars_rad_zm(ivar)) /= 0 &
        1266           0 :                  .and. ivar <= nvarmax_rad_zm )
        1267           0 :         ivar = ivar + 1
        1268             :       end do
        1269           0 :       ntot = ivar - 1
        1270           0 :       if ( ntot == nvarmax_rad_zm ) then
        1271           0 :         write(fstderr,*) "There are more statistical variables listed in ",  &
        1272           0 :                          "vars_rad_zm than allowed for by nvarmax_rad_zm."
        1273           0 :         write(fstderr,*) "Check the number of variables listed for vars_rad_zm ",  &
        1274           0 :                          "in the stats namelist, or change nvarmax_rad_zm."
        1275           0 :         write(fstderr,*) "nvarmax_rad_zm = ", nvarmax_rad_zm
        1276           0 :         write(fstderr,*) "number of variables in vars_rad_zm = ", ntot
        1277           0 :         write(fstderr,*) "stats_init:  number of rad_zm statistical variables exceeds limit"
        1278           0 :         err_code = clubb_fatal_error
        1279           0 :         return
        1280             :       end if
        1281             :
        1282           0 :       stats_rad_zm%num_output_fields = ntot
        1283           0 :       stats_rad_zm%kk = nnrad_zm
        1284           0 :       stats_rad_zm%ii = nlon
        1285           0 :       stats_rad_zm%jj = nlat
        1286             :
        1287           0 :       allocate( stats_rad_zm%z( stats_rad_zm%kk ) )
        1288           0 :       stats_rad_zm%z = grad_zm
        1289             :
        1290           0 :       allocate( stats_rad_zm%accum_field_values( stats_rad_zm%ii, stats_rad_zm%jj, &
        1291           0 :         stats_rad_zm%kk, stats_rad_zm%num_output_fields ) )
        1292           0 :       allocate( stats_rad_zm%accum_num_samples( stats_rad_zm%ii, stats_rad_zm%jj, &
        1293           0 :         stats_rad_zm%kk, stats_rad_zm%num_output_fields ) )
        1294           0 :       allocate( stats_rad_zm%l_in_update( stats_rad_zm%ii, stats_rad_zm%jj, &
        1295           0 :         stats_rad_zm%kk, stats_rad_zm%num_output_fields ) )
        1296             :
        1297             :       call stats_zero( stats_rad_zm%ii, stats_rad_zm%jj, stats_rad_zm%kk, & ! intent(in)
        1298             :                        stats_rad_zm%num_output_fields, & ! intent(in)
        1299             :                        stats_rad_zm%accum_field_values, & ! intent(out)
        1300           0 :                        stats_rad_zm%accum_num_samples, stats_rad_zm%l_in_update ) ! intent(out)
        1301             :
        1302           0 :       allocate( stats_rad_zm%file%grid_avg_var( stats_rad_zm%num_output_fields ) )
        1303           0 :       allocate( stats_rad_zm%file%z( stats_rad_zm%kk ) )
        1304             :
        1305             :       call stats_init_rad_zm( vars_rad_zm,                    & ! intent(in)
        1306             :                               l_error,                        & ! intent(inout)
        1307           0 :                               stats_metadata, stats_rad_zm )    ! intent(inout)
        1308             :
        1309           0 :       fname = trim( stats_metadata%fname_rad_zm )
        1310           0 :       if ( stats_metadata%l_grads ) then
        1311             :
        1312             :         ! Open GrADS files
        1313             :         call open_grads( iunit, fdir, fname,  & ! intent(in)
        1314             :                          1, stats_rad_zm%kk, nlat, nlon, stats_rad_zm%z, & ! intent(in)
        1315             :                          day, month, year, lat_vals, lon_vals, &
        1316             :                          time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
        1317             :                          stats_metadata, stats_rad_zm%num_output_fields, & ! intent(in)
        1318           0 :                          stats_rad_zm%file ) ! intent(inout)
        1319             :
        1320             :       else ! Open NetCDF file
        1321             : #ifdef NETCDF
        1322             :         call open_netcdf_for_writing( nlat, nlon, fdir, fname,  & ! intent(in)
        1323             :                           1, stats_rad_zm%kk, stats_rad_zm%z, & ! intent(in)
        1324             :                           day, month, year, lat_vals, lon_vals, & ! intent(in)
        1325             :                           time_current, stats_metadata%stats_tout, & ! intent(in)
        1326             :                           stats_rad_zm%num_output_fields, & ! intent(in)
        1327             :                           stats_rad_zm%file ) ! intent(inout)
        1328             :
        1329             :         ! Finalize the variable definitions
        1330             :         call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
        1331             :                           l_uv_nudge, & ! intent(in)
        1332             :                           l_tke_aniso, & ! intent(in)
        1333             :                           l_standard_term_ta, & ! intent(in)
        1334             :                           stats_rad_zm%file ) ! intent(inout)
        1335             :
        1336             :         if ( err_code == clubb_fatal_error ) then
        1337             :           write(fstderr,*) "Fatal error setting up stats_rad_zm"
        1338             :           return
        1339             :         end if
        1340             :
        1341             :         if ( err_code == clubb_fatal_error ) return
        1342             : #else
        1343           0 :         error stop "This CLUBB program was not compiled with netCDF support."
        1344             : #endif
        1345             :       end if
        1346             :
        1347             :     end if ! l_output_rad_files
        1348             :
        1349             :
        1350             :     ! Initialize stats_sfc (surface point)
        1351             :
        1352           0 :     ivar = 1
        1353           0 :     do while ( ichar(vars_sfc(ivar)(1:1)) /= 0  &
        1354             :                .and. len_trim(vars_sfc(ivar)) /= 0 &
        1355           0 :                .and. ivar <= nvarmax_sfc )
        1356           0 :       ivar = ivar + 1
        1357             :     end do
        1358           0 :     ntot = ivar - 1
        1359           0 :     if ( ntot == nvarmax_sfc ) then
        1360           0 :       write(fstderr,*) "There are more statistical variables listed in ",  &
        1361           0 :                        "vars_sfc than allowed for by nvarmax_sfc."
        1362           0 :       write(fstderr,*) "Check the number of variables listed for vars_sfc ",  &
        1363           0 :                        "in the stats namelist, or change nvarmax_sfc."
        1364           0 :       write(fstderr,*) "nvarmax_sfc = ", nvarmax_sfc
        1365           0 :       write(fstderr,*) "number of variables in vars_sfc = ", ntot
        1366           0 :       write(fstderr,*) "stats_init:  number of sfc statistical variables exceeds limit"
        1367           0 :       err_code = clubb_fatal_error
        1368           0 :       return
        1369             :
        1370             :     end if
        1371             :
        1372           0 :     stats_sfc%num_output_fields = ntot
        1373           0 :     stats_sfc%kk = 1
        1374           0 :     stats_sfc%ii = nlon
        1375           0 :     stats_sfc%jj = nlat
        1376             :
        1377           0 :     allocate( stats_sfc%z( stats_sfc%kk ) )
        1378           0 :     stats_sfc%z = gzm(1)
        1379             :
        1380           0 :     allocate( stats_sfc%accum_field_values( stats_sfc%ii, stats_sfc%jj, &
        1381           0 :       stats_sfc%kk, stats_sfc%num_output_fields ) )
        1382           0 :     allocate( stats_sfc%accum_num_samples( stats_sfc%ii, stats_sfc%jj, &
        1383           0 :       stats_sfc%kk, stats_sfc%num_output_fields ) )
        1384           0 :     allocate( stats_sfc%l_in_update( stats_sfc%ii, stats_sfc%jj, &
        1385           0 :       stats_sfc%kk, stats_sfc%num_output_fields ) )
        1386             :
        1387             :     call stats_zero( stats_sfc%ii, stats_sfc%jj, stats_sfc%kk, stats_sfc%num_output_fields, & ! In
        1388           0 :       stats_sfc%accum_field_values, stats_sfc%accum_num_samples, stats_sfc%l_in_update ) ! out
        1389             :
        1390           0 :     allocate( stats_sfc%file%grid_avg_var( stats_sfc%num_output_fields ) )
        1391           0 :     allocate( stats_sfc%file%z( stats_sfc%kk ) )
        1392             :
        1393             :     call stats_init_sfc( vars_sfc,                    & ! intent(in)
        1394             :                          l_error,                     & ! intent(inout)
        1395           0 :                          stats_metadata, stats_sfc )    ! intent(inout)
        1396             :
        1397           0 :     fname = trim( stats_metadata%fname_sfc )
        1398             :
        1399           0 :     if ( stats_metadata%l_grads ) then
        1400             :
        1401             :       ! Open GrADS files
        1402             :       call open_grads( iunit, fdir, fname,  & ! intent(in)
        1403             :                        1, stats_sfc%kk, nlat, nlon, stats_sfc%z, & ! intent(in)
        1404             :                        day, month, year, lat_vals, lon_vals, & ! intent(in)
        1405             :                        time_current+real(stats_metadata%stats_tout,kind=time_precision), & !intent(in)
        1406             :                        stats_metadata, stats_sfc%num_output_fields, & ! intent(in)
        1407           0 :                        stats_sfc%file ) ! intent(inout)
        1408             :
        1409             :     else ! Open NetCDF files
        1410             : #ifdef NETCDF
        1411             :       call open_netcdf_for_writing( nlat, nlon, fdir, fname, 1, stats_sfc%kk, stats_sfc%z, &  ! In
        1412             :                         day, month, year, lat_vals, lon_vals, &  ! In
        1413             :                         time_current, stats_metadata%stats_tout, stats_sfc%num_output_fields, &  ! In
        1414             :                         stats_sfc%file ) ! InOut
        1415             :
        1416             :       ! Finalize the variable definitions
        1417             :       call first_write( clubb_params, sclr_dim, sclr_tol, & ! intent(in)
        1418             :                         l_uv_nudge, & ! intent(in)
        1419             :                         l_tke_aniso, & ! intent(in)
        1420             :                         l_standard_term_ta, & ! intent(in)
        1421             :                         stats_sfc%file ) ! intent(inout)
        1422             :
        1423             :       if ( err_code == clubb_fatal_error ) then
        1424             :         write(fstderr,*) "Fatal error setting up stats_sfc"
        1425             :         return
        1426             :       end if
        1427             : #else
        1428           0 :       error stop "This CLUBB program was not compiled with netCDF support."
        1429             : #endif
        1430             :     end if
        1431             :
        1432             :     ! Check for errors
        1433             :
        1434           0 :     if ( l_error ) then
        1435           0 :       write(fstderr,*) 'stats_init:  errors found'
        1436           0 :       err_code = clubb_fatal_error
        1437           0 :       return
        1438             :     endif
        1439             :
        1440             :     return
        1441             :
        1442             :     ! If namelist was not found in input file, turn off statistics
        1443             :
        1444             :     100 continue
        1445           0 :     write(fstderr,*) 'Error with statsnl, statistics is turned off'
        1446           0 :     stats_metadata%l_stats       = .false.
        1447           0 :     stats_metadata%l_stats_samp  = .false.
        1448           0 :     stats_metadata%l_stats_last  = .false.
        1449             :
        1450           0 :     if ( err_code == clubb_fatal_error ) error stop
        1451             :
        1452             :     return
        1453             :
        1454             :   end subroutine stats_init
        1455             :
        1456             :   !-----------------------------------------------------------------------
        1457           0 :   subroutine stats_zero( ii, jj, kk, nn, &
        1458           0 :                          x, n, l_in_update )
        1459             :
        1460             :     ! Description:
        1461             :     !   Initialize stats to zero
        1462             :     ! References:
        1463             :     !   None
        1464             :     !-----------------------------------------------------------------------
        1465             :     use clubb_precision, only: &
        1466             :         stat_rknd,   & ! Variable(s)
        1467             :         stat_nknd
        1468             :
        1469             :     implicit none
        1470             :
        1471             :     ! Input Variable(s)
        1472             :     integer, intent(in) :: ii, jj, kk, nn
        1473             :
        1474             :     ! Output Variable(s)
        1475             :     real(kind=stat_rknd), dimension(ii,jj,kk,nn), intent(out)    :: x
        1476             :     integer(kind=stat_nknd), dimension(ii,jj,kk,nn), intent(out) :: n
        1477             :     logical, dimension(ii,jj,kk,nn), intent(out) :: l_in_update
        1478             :
        1479             :     ! Zero out arrays
        1480             :
        1481           0 :     if ( nn > 0 ) then
        1482           0 :       x(:,:,:,:) = 0.0_stat_rknd
        1483           0 :       n(:,:,:,:) = 0_stat_nknd
        1484           0 :       l_in_update(:,:,:,:) = .false.
        1485             :     end if
        1486             :
        1487           0 :     return
        1488             :   end subroutine stats_zero
        1489             :
        1490             :   !-----------------------------------------------------------------------
        1491           0 :   subroutine stats_avg( ii, jj, kk, nn, n, &
        1492           0 :                         x )
        1493             :
        1494             :     ! Description:
        1495             :     !   Compute the average of stats fields
        1496             :     ! References:
        1497             :     !   None
        1498             :     !-----------------------------------------------------------------------
        1499             :     use clubb_precision, only: &
        1500             :         stat_rknd,   & ! Variable(s)
        1501             :         stat_nknd
        1502             :
        1503             :     use stat_file_module, only: &
        1504             :         clubb_i, clubb_j ! Variable(s)
        1505             :
        1506             :     implicit none
        1507             :
        1508             :     ! External
        1509             :     intrinsic :: real
        1510             :
        1511             :     ! Input Variable(s)
        1512             :     integer, intent(in) :: &
        1513             :       ii, & ! Number of points in X (i.e. latitude) dimension
        1514             :       jj, & ! Number of points in Y (i.e. longitude) dimension
        1515             :       kk, & ! Number of levels in vertical (i.e. Z) dimension
        1516             :       nn    ! Number of variables being output to disk (e.g. cloud_frac, rain rate, etc.)
        1517             :
        1518             :     integer(kind=stat_nknd), dimension(ii,jj,kk,nn), intent(in) :: &
        1519             :       n ! n is the number of samples for each of the nn fields
        1520             :         ! and each of the kk vertical levels
        1521             :
        1522             :     ! Output Variable(s)
        1523             :     real(kind=stat_rknd), dimension(ii,jj,kk,nn), intent(inout) :: &
        1524             :       x ! The variable x contains the cumulative sums of n sample values of each of
        1525             :         ! the nn output fields (e.g. the sum of the sampled rain rate values)
        1526             :
        1527             :     ! ---- Begin Code ----
        1528             :
        1529             :     ! Compute averages
        1530           0 :     where ( n(1,1,1:kk,1:nn) > 0 )
        1531           0 :       x(clubb_i,clubb_j,1:kk,1:nn) = x(clubb_i,clubb_j,1:kk,1:nn) &
        1532             :          / real( n(clubb_i,clubb_j,1:kk,1:nn), kind=stat_rknd )
        1533             :     end where
        1534             :
        1535           0 :     return
        1536             :   end subroutine stats_avg
        1537             :
        1538             :   !-----------------------------------------------------------------------
        1539           0 :   subroutine stats_begin_timestep( itime, stats_nsamp, stats_nout, &
        1540             :                                    stats_metadata )
        1541             :
        1542             :     !     Description:
        1543             :     !       Given the elapsed time, set flags determining specifics such as
        1544             :     !       if this time set should be sampled or if this is the first or
        1545             :     !       last time step.
        1546             :     !-----------------------------------------------------------------------
        1547             :
        1548             :     use stats_variables, only: &
        1549             :       stats_metadata_type
        1550             :
        1551             :     implicit none
        1552             :
        1553             :     ! External
        1554             :     intrinsic :: mod
        1555             :
        1556             :     ! Input Variable(s)
        1557             :     integer, intent(in) ::  &
        1558             :       itime, &       ! Elapsed model time       [timestep]
        1559             :       stats_nsamp, & ! Stats sampling interval  [timestep]
        1560             :       stats_nout     ! Stats output interval    [timestep]
        1561             :
        1562             :     type (stats_metadata_type), intent(inout) :: &
        1563             :       stats_metadata
        1564             :
        1565           0 :     if ( .not. stats_metadata%l_stats ) return
        1566             :
        1567             :     ! Only sample time steps that are multiples of "stats_tsamp"
        1568             :     ! in a case's "model.in" file to shorten length of run
        1569           0 :     if ( mod( itime, stats_nsamp ) == 0 ) then
        1570           0 :       stats_metadata%l_stats_samp = .true.
        1571             :     else
        1572           0 :       stats_metadata%l_stats_samp = .false.
        1573             :     end if
        1574             :
        1575             :     ! Indicates the end of the sampling time period. Signals to start writing to the file
        1576           0 :     if ( mod( itime, stats_nout ) == 0 ) then
        1577           0 :       stats_metadata%l_stats_last = .true.
        1578             :     else
        1579           0 :       stats_metadata%l_stats_last = .false.
        1580             :     end if
        1581             :
        1582             :     return
        1583             :
        1584             :   end subroutine stats_begin_timestep
        1585             :
        1586             :   !-----------------------------------------------------------------------
        1587           0 :   subroutine stats_end_timestep( clubb_params, stats_metadata,  & ! intent(in)
        1588             :                                  stats_zt, stats_zm, stats_sfc, & ! intent(inout)
        1589             :                                  stats_lh_zt, stats_lh_sfc,     & ! intent(inout)
        1590             :                                  stats_rad_zt, stats_rad_zm     & ! intent(inout)
        1591             : #ifdef NETCDF
        1592             :                                  , l_uv_nudge, &
        1593             :                                  l_tke_aniso, &
        1594             :                                  l_standard_term_ta &
        1595             : #endif
        1596             :                                   )
        1597             :
        1598             :     ! Description:
        1599             :     !   Called when the stats timestep has ended. This subroutine
        1600             :     !   is responsible for calling statistics to be written to the output
        1601             :     !   format.
        1602             :     !
        1603             :     ! References:
        1604             :     !   None
        1605             :     !-----------------------------------------------------------------------
        1606             :
        1607             :     use clubb_precision, only: &
        1608             :         core_rknd    ! Variable(s)
        1609             :
        1610             :     use constants_clubb, only: &
        1611             :         fstderr ! Constant(s)
        1612             :
        1613             :     use stats_variables, only: &
        1614             :         stats_metadata_type
        1615             :
        1616             :     use output_grads, only: &
        1617             :         write_grads ! Procedure(s)
        1618             :
        1619             :     use stat_file_module, only: &
        1620             :         clubb_i, & ! Variable(s)
        1621             :         clubb_j
        1622             :
        1623             :     use parameter_indices, only: &
        1624             :         nparams    ! Variable(s)
        1625             :
        1626             : #ifdef NETCDF
        1627             :     use output_netcdf, only: &
        1628             :         write_netcdf ! Procedure(s)
        1629             : #endif
        1630             :
        1631             :     use error_code, only : &
        1632             :         err_code, &         ! Error Indicator
        1633             :         clubb_fatal_error   ! Constant
        1634             :
        1635             :     use stats_type, only: stats ! Type
        1636             :
        1637             :     implicit none
        1638             :
        1639             :     real( kind = core_rknd ), dimension(nparams), intent(in) :: &
        1640             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
        1641             :
        1642             :     type (stats_metadata_type), intent(in) :: &
        1643             :       stats_metadata
        1644             :
        1645             :     type (stats), target, intent(inout) :: &
        1646             :       stats_zt, &
        1647             :       stats_zm, &
        1648             :       stats_sfc, &
        1649             :       stats_lh_zt, &
        1650             :       stats_lh_sfc, &
        1651             :       stats_rad_zt, &
        1652             :       stats_rad_zm
        1653             :
        1654             :     ! External
        1655             :     intrinsic :: floor
        1656             :
        1657             : #ifdef NETCDF
        1658             :     ! Input Variables
        1659             :     logical, intent(in) :: &
        1660             :       l_uv_nudge,         & ! For wind speed nudging.
        1661             :       l_tke_aniso,        & ! For anisotropic turbulent kinetic energy, i.e.
        1662             :                             ! TKE = 1/2 (u'^2 + v'^2 + w'^2)
        1663             :       l_standard_term_ta    ! Use the standard discretization for the turbulent advection terms.
        1664             :                             ! Setting to .false. means that a_1 and a_3 are pulled outside of the
        1665             :                             ! derivative in advance_wp2_wp3_module.F90 and in
        1666             :                             ! advance_xp2_xpyp_module.F90.
        1667             : #endif
        1668             :
        1669             :     ! Local Variables
        1670             :
        1671             :     logical :: l_error
        1672             :
        1673             :     ! ---- Begin Code ----
        1674             :
        1675             :     ! Check if it is time to write to file
        1676             :
        1677           0 :     if ( .not. stats_metadata%l_stats_last ) return
        1678             :
        1679             :     ! Initialize
        1680           0 :     l_error = .false.
        1681             :
        1682             :     call stats_check_num_samples( stats_zt, stats_metadata,  & ! intent(in)
        1683           0 :                                   l_error )                    ! intent(inout)
        1684             :     call stats_check_num_samples( stats_zm, stats_metadata,  & ! intent(in)
        1685           0 :                                   l_error )                    ! intent(inout)
        1686             :     call stats_check_num_samples( stats_sfc, stats_metadata, & ! intent(in)
        1687           0 :                                   l_error )                    ! intent(inout)
        1688           0 :     if ( stats_metadata%l_silhs_out ) then
        1689             :       call stats_check_num_samples( stats_lh_zt, stats_metadata,  & ! intent(in)
        1690           0 :                                     l_error )                       ! intent(inout)
        1691             :       call stats_check_num_samples( stats_lh_sfc, stats_metadata, & ! intent(in)
        1692           0 :                                     l_error )                       ! intent(inout)
        1693             :     end if
        1694           0 :     if ( stats_metadata%l_output_rad_files ) then
        1695             :       call stats_check_num_samples( stats_rad_zt, stats_metadata, & ! intent(in)
        1696           0 :                                     l_error )                       ! intent(inout)
        1697             :       call stats_check_num_samples( stats_rad_zm, stats_metadata, & ! intent(in)
        1698           0 :                                     l_error )                       ! intent(inout)
        1699             :     end if
        1700             :
        1701             :     ! Return if errors are found.
        1702           0 :     if ( l_error ) then
        1703           0 :       write(fstderr,*) 'Possible statistical sampling error'
        1704           0 :       write(fstderr,*) 'For details, set debug_level to a value of at ',  &
        1705           0 :                        'least 1 in the appropriate model.in file.'
        1706           0 :       write(fstderr,*) 'stats_end_timestep:  error(s) found'
        1707           0 :       err_code = clubb_fatal_error
        1708           0 :       return
        1709             :     end if ! l_error
        1710             :
        1711             :     ! Compute averages
        1712             :     call stats_avg( stats_zt%ii, stats_zt%jj, stats_zt%kk, stats_zt%num_output_fields, & ! In
        1713             :                     stats_zt%accum_num_samples, & ! intent(in)
        1714           0 :                     stats_zt%accum_field_values ) ! intent(inout)
        1715             :     call stats_avg( stats_zm%ii, stats_zm%jj, stats_zm%kk, stats_zm%num_output_fields, & ! In
        1716             :                     stats_zm%accum_num_samples, & ! intent(in)
        1717           0 :                     stats_zm%accum_field_values ) ! intent(inout)
        1718           0 :     if ( stats_metadata%l_silhs_out ) then
        1719             :       call stats_avg( stats_lh_zt%ii, stats_lh_zt%jj, stats_lh_zt%kk, & ! intent(in)
        1720             :          stats_lh_zt%num_output_fields, stats_lh_zt%accum_num_samples, & ! intent(in)
        1721           0 :          stats_lh_zt%accum_field_values ) ! intent(inout)
        1722             :       call stats_avg( stats_lh_sfc%ii, stats_lh_sfc%jj, stats_lh_sfc%kk, & ! intent(in)
        1723             :         stats_lh_sfc%num_output_fields, stats_lh_sfc%accum_num_samples, & ! intent(in)
        1724           0 :         stats_lh_sfc%accum_field_values ) ! intent(inout)
        1725             :     end if
        1726           0 :     if ( stats_metadata%l_output_rad_files ) then
        1727             :       call stats_avg( stats_rad_zt%ii, stats_rad_zt%jj, stats_rad_zt%kk, & ! intent(in)
        1728             :         stats_rad_zt%num_output_fields, & ! intent(in)
        1729             :         stats_rad_zt%accum_num_samples, & ! intent(in)
        1730           0 :         stats_rad_zt%accum_field_values ) ! intent(inout)
        1731             :       call stats_avg( stats_rad_zm%ii, stats_rad_zm%jj, stats_rad_zm%kk, & ! intent(in)
        1732             :         stats_rad_zm%num_output_fields, & ! intent(in)
        1733             :         stats_rad_zm%accum_num_samples, & ! intent(in)
        1734           0 :         stats_rad_zm%accum_field_values ) ! intent(inout)
        1735             :     end if
        1736             :     call stats_avg( stats_sfc%ii, stats_sfc%jj, stats_sfc%kk, stats_sfc%num_output_fields, & ! In
        1737             :         stats_sfc%accum_num_samples, & ! intent(in)
        1738           0 :         stats_sfc%accum_field_values ) ! intent(inout)
        1739             :
        1740             :     ! Only write to the file and zero out the stats fields if we've reach the horizontal
        1741             :     ! limits of the domain (this is always true in the single-column case because it's 1x1).
        1742           0 :     if ( clubb_i == stats_zt%ii .and. clubb_j == stats_zt%jj ) then
        1743             :       ! Write to file
        1744           0 :       if ( stats_metadata%l_grads ) then
        1745           0 :         call write_grads( stats_zt%file  ) ! intent(inout)
        1746           0 :         call write_grads( stats_zm%file  ) ! intent(inout)
        1747           0 :         if ( stats_metadata%l_silhs_out ) then
        1748           0 :           call write_grads( stats_lh_zt%file  ) ! intent(inout)
        1749           0 :           call write_grads( stats_lh_sfc%file  ) ! intent(inout)
        1750             :         end if
        1751           0 :         if ( stats_metadata%l_output_rad_files ) then
        1752           0 :           call write_grads( stats_rad_zt%file  ) ! intent(inout)
        1753           0 :           call write_grads( stats_rad_zm%file  ) ! intent(inout)
        1754             :         end if
        1755           0 :         call write_grads( stats_sfc%file  ) ! intent(inout)
        1756             :       else ! l_netcdf
        1757             :
        1758             : #ifdef NETCDF
        1759             :         call write_netcdf( stats_zt%file  ) ! intent(inout)
        1760             :
        1761             :         call write_netcdf( stats_zm%file  ) ! intent(inout)
        1762             :
        1763             :         if ( stats_metadata%l_silhs_out ) then
        1764             :
        1765             :           call write_netcdf( stats_lh_zt%file  ) ! intent(inout)
        1766             :
        1767             :           call write_netcdf( stats_lh_sfc%file  ) ! intent(inout)
        1768             :
        1769             :         end if
        1770             :         if ( stats_metadata%l_output_rad_files ) then
        1771             :
        1772             :           call write_netcdf( stats_rad_zt%file  ) ! intent(inout)
        1773             :
        1774             :           call write_netcdf( stats_rad_zm%file  ) ! intent(inout)
        1775             :
        1776             :         end if
        1777             :
        1778             :         call write_netcdf( stats_sfc%file  ) ! intent(inout)
        1779             :
        1780             :         if ( err_code == clubb_fatal_error ) return
        1781             : #else
        1782           0 :         error stop "This program was not compiled with netCDF support"
        1783             : #endif /* NETCDF */
        1784             :       end if ! l_grads
        1785             :
        1786             :       ! Reset sample fields
        1787             :       call stats_zero( stats_zt%ii, stats_zt%jj, stats_zt%kk, stats_zt%num_output_fields, & ! In
        1788           0 :       stats_zt%accum_field_values, stats_zt%accum_num_samples, stats_zt%l_in_update ) ! out
        1789             :       call stats_zero( stats_zm%ii, stats_zm%jj, stats_zm%kk, stats_zm%num_output_fields, & ! In
        1790           0 :         stats_zm%accum_field_values, stats_zm%accum_num_samples, stats_zm%l_in_update ) ! Out
        1791           0 :       if ( stats_metadata%l_silhs_out ) then
        1792             :         call stats_zero( stats_lh_zt%ii, stats_lh_zt%jj, stats_lh_zt%kk, & ! intent(in)
        1793             :           stats_lh_zt%num_output_fields, & ! intent(in)
        1794             :           stats_lh_zt%accum_field_values, & ! intent(out)
        1795           0 :           stats_lh_zt%accum_num_samples, stats_lh_zt%l_in_update ) ! intent(out)
        1796             :         call stats_zero( stats_lh_sfc%ii, stats_lh_sfc%jj, stats_lh_sfc%kk, & ! intent(in)
        1797             :           stats_lh_sfc%num_output_fields, & ! intent(in)
        1798             :           stats_lh_sfc%accum_field_values, & ! intent(out)
        1799           0 :           stats_lh_sfc%accum_num_samples, stats_lh_sfc%l_in_update ) ! intent(out)
        1800             :       end if
        1801           0 :       if ( stats_metadata%l_output_rad_files ) then
        1802             :         call stats_zero( stats_rad_zt%ii, stats_rad_zt%jj, stats_rad_zt%kk, & ! intent(in)
        1803             :           stats_rad_zt%num_output_fields, & ! intent(in)
        1804             :           stats_rad_zt%accum_field_values, & ! intent(out)
        1805           0 :           stats_rad_zt%accum_num_samples, stats_rad_zt%l_in_update ) ! intent(out)
        1806             :         call stats_zero( stats_rad_zt%ii, stats_rad_zt%jj, stats_rad_zm%kk, & ! intent(in)
        1807             :           stats_rad_zm%num_output_fields, & ! intent(in)
        1808             :           stats_rad_zm%accum_field_values, & ! intent(out)
        1809           0 :           stats_rad_zm%accum_num_samples, stats_rad_zm%l_in_update ) ! intent(out)
        1810             :       end if
        1811             :       call stats_zero( stats_sfc%ii, stats_sfc%jj, stats_sfc%kk, stats_sfc%num_output_fields, & !IN
        1812             :         stats_sfc%accum_field_values, & ! intent(out)
        1813           0 :         stats_sfc%accum_num_samples, stats_sfc%l_in_update ) ! intent(out)
        1814             :
        1815             :     end if ! clubb_i = stats_zt%ii .and. clubb_j == stats_zt%jj
        1816             :
        1817             :
        1818             :     return
        1819             :   end subroutine stats_end_timestep
        1820             :
        1821             :   !----------------------------------------------------------------------
        1822           0 :   subroutine stats_accumulate( &
        1823             :                      nz, sclr_dim, edsclr_dim, &
        1824           0 :                      invrs_dzm, zt, dzm, dzt, dt, &
        1825           0 :                      um, vm, upwp, vpwp, up2, vp2, &
        1826           0 :                      thlm, rtm, wprtp, wpthlp, &
        1827           0 :                      wp2, wp3, rtp2, rtp3, thlp2, thlp3, rtpthlp, &
        1828           0 :                      wpthvp, wp2thvp, rtpthvp, thlpthvp, &
        1829           0 :                      p_in_Pa, exner, rho, rho_zm, &
        1830           0 :                      rho_ds_zm, rho_ds_zt, thv_ds_zm, thv_ds_zt, &
        1831           0 :                      wm_zt, wm_zm, rcm, wprcp, rc_coef, rc_coef_zm, &
        1832           0 :                      rcm_zm, rtm_zm, thlm_zm, cloud_frac, ice_supersat_frac, &
        1833           0 :                      cloud_frac_zm, ice_supersat_frac_zm, rcm_in_layer, &
        1834           0 :                      cloud_cover, rcm_supersat_adj, sigma_sqd_w, &
        1835           0 :                      thvm, ug, vg, Lscale, wpthlp2, wp2thlp, wprtp2, wp2rtp, &
        1836           0 :                      Lscale_up, Lscale_down, tau_zt, Kh_zt, wp2rcp, &
        1837           0 :                      wprtpthlp, sigma_sqd_w_zt, rsat, wp2_zt, thlp2_zt, &
        1838           0 :                      wpthlp_zt, wprtp_zt, rtp2_zt, rtpthlp_zt, up2_zt, &
        1839           0 :                      vp2_zt, upwp_zt, vpwp_zt, wpup2, wpvp2, &
        1840           0 :                      wp2up2, wp2vp2, wp4, &
        1841           0 :                      tau_zm, Kh_zm, thlprcp, &
        1842           0 :                      rtprcp, rcp2, em, a3_coef, a3_coef_zt, &
        1843           0 :                      wp3_zm, wp3_on_wp2, wp3_on_wp2_zt, Skw_velocity, &
        1844           0 :                      w_up_in_cloud, w_down_in_cloud, cloudy_updraft_frac, &
        1845           0 :                      cloudy_downdraft_frac, pdf_params, pdf_params_zm, &
        1846           0 :                      sclrm, sclrp2, sclrprtp, sclrpthlp, sclrm_forcing, sclrpthvp, &
        1847           0 :                      wpsclrp, sclrprcp, wp2sclrp, wpsclrp2, wpsclrprtp, &
        1848           0 :                      wpsclrpthlp, wpedsclrp, edsclrm, edsclrm_forcing, &
        1849             :                      saturation_formula, &
        1850             :                      stats_metadata, &
        1851             :                      stats_zt, stats_zm, stats_sfc )
        1852             :
        1853             :     ! Description:
        1854             :     !   Accumulate those stats variables that are preserved in CLUBB from timestep to
        1855             :     !   timestep, but not those stats that are not, (e.g. budget terms, longwave and
        1856             :     !   shortwave components, etc.)
        1857             :     !
        1858             :     ! References:
        1859             :     !   None
        1860             :     !----------------------------------------------------------------------
        1861             :
        1862             :     use constants_clubb, only: &
        1863             :         cloud_frac_min, &  ! Constant
        1864             :         eps
        1865             :
        1866             :     use pdf_utilities, only: &
        1867             :         compute_variance_binormal    ! Procedure
        1868             :
        1869             :     use stats_variables, only: &
        1870             :         stats_metadata_type
        1871             :
        1872             :     use grid_class, only: &
        1873             :         zt2zm ! Procedure(s)
        1874             :
        1875             :     use pdf_parameter_module, only: &
        1876             :         pdf_parameter ! Type
        1877             :
        1878             :     use T_in_K_module, only: &
        1879             :         thlm2T_in_K ! Procedure
        1880             :
        1881             :     use constants_clubb, only: &
        1882             :         rc_tol, fstderr    ! Constant(s)
        1883             :
        1884             :     use stats_type_utilities, only: &
        1885             :         stat_update_var,  & ! Procedure(s)
        1886             :         stat_update_var_pt
        1887             :
        1888             :     use advance_helper_module, only: &
        1889             :         vertical_avg, &     ! Procedure(s)
        1890             :         vertical_integral
        1891             :
        1892             :     use interpolation, only: &
        1893             :         lin_interpolate_two_points             ! Procedure
        1894             :
        1895             :     use saturation, only: &
        1896             :         sat_mixrat_ice ! Procedure
        1897             :
        1898             :     use clubb_precision, only: &
        1899             :         core_rknd ! Variable(s)
        1900             :
        1901             :     use stats_type, only: stats ! Type
        1902             :
        1903             :     implicit none
        1904             :
        1905             :     ! Input Variable(s)
        1906             :     integer, intent(in) :: &
        1907             :       nz, &
        1908             :       sclr_dim, &
        1909             :       edsclr_dim
        1910             :
        1911             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1912             :       invrs_dzm, & ! The inverse spacing between thermodynamic grid
        1913             :                    ! levels; centered over momentum grid levels.
        1914             :       zt,        & ! Thermo grid
        1915             :       dzm,       & ! Spacing between thermodynamic grid levels; centered over
        1916             :                    ! momentum grid levels
        1917             :       dzt          ! Spcaing between momentum grid levels; centered over
        1918             :                    ! thermodynamic grid levels
        1919             :
        1920             :     real( kind = core_rknd ), intent(in) ::  &
        1921             :       dt           ! Model timestep                        [s]
        1922             :
        1923             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1924             :       um,       & ! u wind (thermodynamic levels)          [m/s]
        1925             :       vm,       & ! v wind (thermodynamic levels)          [m/s]
        1926             :       upwp,     & ! vertical u momentum flux (m levs.)     [m^2/s^2]
        1927             :       vpwp,     & ! vertical v momentum flux (m levs.)     [m^2/s^2]
        1928             :       up2,      & ! < u'^2 > (momentum levels)             [m^2/s^2]
        1929             :       vp2,      & ! < v'^2 > (momentum levels)             [m^2/s^2]
        1930             :       thlm,     & ! liquid potential temperature (t levs.) [K]
        1931             :       rtm,      & ! total water mixing ratio (t levs.)     [kg/kg]
        1932             :       wprtp,    & ! < w' r_t' > (momentum levels)          [m/s kg/kg]
        1933             :       wpthlp,   & ! < w' th_l' > (momentum levels)         [m/s K]
        1934             :       wp2,      & ! < w'^2 > (momentum levels)             [m^2/s^2]
        1935             :       wp3,      & ! < w'^3 > (thermodynamic levels)        [m^3/s^3]
        1936             :       rtp2,     & ! < r_t'^2 > (momentum levels)           [(kg/kg)^2]
        1937             :       rtp3,     & ! < r_t'^3 > (thermodynamic levels)      [(kg/kg)^3]
        1938             :       thlp2,    & ! < th_l'^2 > (momentum levels)          [K^2]
        1939             :       thlp3,    & ! < th_l'^3 > (thermodynamic levels)     [K^3]
        1940             :       rtpthlp,  & ! < r_t' th_l' > (momentum levels)       [kg/kg K]
        1941             :       wpthvp,   & ! < w' th_v' > (momentum levels)         [kg/kg K]
        1942             :       wp2thvp,  & ! < w'^2 th_v' > (thermodynamic levels)  [m^2/s^2 K]
        1943             :       rtpthvp,  & ! < r_t' th_v' > (momentum levels)       [kg/kg K]
        1944             :       thlpthvp    ! < th_l' th_v' > (momentum levels)      [K^2]
        1945             :
        1946             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1947             :       p_in_Pa,      & ! Pressure (Pa) on thermodynamic points    [Pa]
        1948             :       exner,        & ! Exner function = ( p / p0 ) ** kappa     [-]
        1949             :       rho,          & ! Density (thermodynamic levels)           [kg/m^3]
        1950             :       rho_zm,       & ! Density on momentum levels               [kg/m^3]
        1951             :       rho_ds_zm,    & ! Dry, static density (momentum levels)    [kg/m^3]
        1952             :       rho_ds_zt,    & ! Dry, static density (thermo. levs.)      [kg/m^3]
        1953             :       thv_ds_zm,    & ! Dry, base-state theta_v (momentum levs.) [K]
        1954             :       thv_ds_zt,    & ! Dry, base-state theta_v (thermo. levs.)  [K]
        1955             :       wm_zt,        & ! w on thermodynamic levels                [m/s]
        1956             :       wm_zm           ! w on momentum levels                     [m/s]
        1957             :
        1958             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1959             :       rcm_zm,               & ! Cloud water mixing ratio on m levs.      [kg/kg]
        1960             :       rtm_zm,               & ! Total water mixing ratio on m levs.      [kg/kg]
        1961             :       thlm_zm,              & ! Liquid potential temperature on m levs.  [K]
        1962             :       rcm,                  & ! Cloud water mixing ratio (t levs.)       [kg/kg]
        1963             :       wprcp,                & ! < w' r_c' > (momentum levels)            [m/s kg/kg]
        1964             :       rc_coef,              & ! Coefficient of X'r_c' (t-levs.)      [K/(kg/kg)]
        1965             :       rc_coef_zm,           & ! Coefficient of X'r_c' on m-levs.     [K/(kg/kg)]
        1966             :       cloud_frac,           & ! Cloud fraction (thermodynamic levels)    [-]
        1967             :       ice_supersat_frac,    & ! Ice cloud fracion (thermodynamic levels) [-]
        1968             :       cloud_frac_zm,        & ! Cloud fraction on zm levels              [-]
        1969             :       ice_supersat_frac_zm, & ! Ice cloud fraction on zm levels          [-]
        1970             :       rcm_in_layer,         & ! Cloud water mixing ratio in cloud layer  [kg/kg]
        1971             :       cloud_cover,          & ! Cloud cover                              [-]
        1972             :       rcm_supersat_adj        ! rcm adjustment due to supersaturation    [kg/kg]
        1973             :
        1974             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1975             :       sigma_sqd_w    ! PDF width parameter (momentum levels)    [-]
        1976             :
        1977             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1978             :         thvm,           & ! Virtual potential temperature        [K]
        1979             :         ug,             & ! u geostrophic wind                   [m/s]
        1980             :         vg,             & ! v geostrophic wind                   [m/s]
        1981             :         Lscale,         & ! Length scale                         [m]
        1982             :         wpthlp2,        & ! w'thl'^2                             [m K^2/s]
        1983             :         wp2thlp,        & ! w'^2 thl'                            [m^2 K/s^2]
        1984             :         wprtp2,         & ! w'rt'^2                              [m/s kg^2/kg^2]
        1985             :         wp2rtp,         & ! w'^2rt'                              [m^2/s^2 kg/kg]
        1986             :         Lscale_up,      & ! Length scale (upwards component)     [m]
        1987             :         Lscale_down,    & ! Length scale (downwards component)   [m]
        1988             :         tau_zt,         & ! Eddy diss. time scale; thermo. levs. [s]
        1989             :         Kh_zt,          & ! Eddy diff. coef. on thermo. levels   [m^2/s]
        1990             :         wp2rcp,         & ! w'^2 rc'                             [m^2/s^2 kg/kg]
        1991             :         wprtpthlp,      & ! w'rt'thl'                            [m/s kg/kg K]
        1992             :         sigma_sqd_w_zt, & ! PDF width parameter (thermo. levels) [-]
        1993             :         rsat              ! Saturation mixing ratio              [kg/kg]
        1994             :
        1995             :     real( kind = core_rknd ), intent(in), dimension(nz) :: &
        1996             :       wp2_zt,                & ! w'^2 on thermo. grid                  [m^2/s^2]
        1997             :       thlp2_zt,              & ! thl'^2 on thermo. grid                [K^2]
        1998             :       wpthlp_zt,             & ! w'thl' on thermo. grid                [m K/s]
        1999             :       wprtp_zt,              & ! w'rt' on thermo. grid                 [m kg/(kg s)]
        2000             :       rtp2_zt,               & ! rt'^2 on therm. grid                  [(kg/kg)^2]
        2001             :       rtpthlp_zt,            & ! rt'thl' on thermo. grid               [kg K/kg]
        2002             :       up2_zt,                & ! u'^2 on thermo. grid                  [m^2/s^2]
        2003             :       vp2_zt,                & ! v'^2 on thermo. grid                  [m^2/s^2]
        2004             :       upwp_zt,               & ! u'w' on thermo. grid                  [m^2/s^2]
        2005             :       vpwp_zt,               & ! v'w' on thermo. grid                  [m^2/s^2]
        2006             :       wpup2,                 & ! w'u'^2 (thermodynamic levels)         [m^3/s^3]
        2007             :       wpvp2,                 & ! w'v'^2 (thermodynamic levels)         [m^3/s^3]
        2008             :       wp2up2,                & ! < w'^2u'^2 > (momentum levels)        [m^4/s^4]
        2009             :       wp2vp2,                & ! < w'^2v'^2 > (momentum levels)        [m^4/s^4]
        2010             :       wp4,                   & ! < w'^4 > (momentum levels)            [m^4/s^4]
        2011             :       tau_zm,                & ! Eddy diss. time scale; momentum levs. [s]
        2012             :       Kh_zm,                 & ! Eddy diff. coef. on momentum levels   [m^2/s]
        2013             :       thlprcp,               & ! thl'rc'                               [K kg/kg]
        2014             :       rtprcp,                & ! rt'rc'                                [kg^2/kg^2]
        2015             :       rcp2,                  & ! rc'^2                                 [kg^2/kg^2]
        2016             :       em,                    & ! Turbulent Kinetic Energy (TKE)        [m^2/s^2]
        2017             :       a3_coef,               & ! The a3 coefficient from CLUBB eqns    [-]
        2018             :       a3_coef_zt,            & ! The a3 coef. interp. to the zt grid   [-]
        2019             :       wp3_zm,                & ! w'^3 interpolated to momentum levels  [m^3/s^3]
        2020             :       wp3_on_wp2,            & ! w'^3 / w'^2 on the zm grid            [m/s]
        2021             :       wp3_on_wp2_zt,         & ! w'^3 / w'^2 on the zt grid            [m/s]
        2022             :       Skw_velocity,          & ! Skewness velocity                     [m/s]
        2023             :       w_up_in_cloud,         & ! Mean cloudy updraft speed             [m/s]
        2024             :       w_down_in_cloud,       & ! Mean cloudy downdraft speed           [m/s]
        2025             :       cloudy_updraft_frac,   & ! Cloudy updraft fraction               [-]
        2026             :       cloudy_downdraft_frac    ! Cloudy downdraft fraction             [-]
        2027             :
        2028             :     type(pdf_parameter), intent(in) :: &
        2029             :       pdf_params,    & ! PDF parameters (thermodynamic levels)    [units vary]
        2030             :       pdf_params_zm    ! PDF parameters on momentum levels        [units vary]
        2031             :
        2032             :     real( kind = core_rknd ), intent(in), dimension(nz,sclr_dim) :: &
        2033             :       sclrm,           & ! High-order passive scalar            [units vary]
        2034             :       sclrp2,          & ! High-order passive scalar variance   [units^2]
        2035             :       sclrprtp,        & ! High-order passive scalar covariance [units kg/kg]
        2036             :       sclrpthlp,       & ! High-order passive scalar covariance [units K]
        2037             :       sclrm_forcing,   & ! Large-scale forcing of scalar        [units/s]
        2038             :       sclrpthvp,       & ! High-order passive scalar covariance [units K]
        2039             :       wpsclrp            ! w'sclr'                              [units m/s]
        2040             :
        2041             :     real( kind = core_rknd ), intent(in), dimension(nz,sclr_dim) :: &
        2042             :       sclrprcp,    & ! sclr'rc'     [units vary]
        2043             :       wp2sclrp,    & ! w'^2 sclr'   [units vary]
        2044             :       wpsclrp2,    & ! w'sclr'^2    [units vary]
        2045             :       wpsclrprtp,  & ! w'sclr'rt'   [units vary]
        2046             :       wpsclrpthlp    ! w'sclr'thl'  [units vary]
        2047             :
        2048             :     real( kind = core_rknd ), intent(in), dimension(nz,edsclr_dim) :: &
        2049             :       wpedsclrp,       & ! w'edsclr'                        [units vary]
        2050             :       edsclrm,         & ! Eddy-diff passive scalar         [units vary]
        2051             :       edsclrm_forcing    ! Large-scale forcing of edscalar  [units vary]
        2052             :
        2053             :     integer, intent(in) :: &
        2054             :       saturation_formula
        2055             :
        2056             :     type (stats_metadata_type), intent(in) :: &
        2057             :       stats_metadata
        2058             :
        2059             :     type (stats), target, intent(inout) :: &
        2060             :       stats_zt, &
        2061             :       stats_zm, &
        2062             :       stats_sfc
        2063             :
        2064             :     ! Local Variables
        2065             :
        2066             :     integer :: isclr, k
        2067             :     integer :: grid_level = 1  ! grid level for stats where there is only one sensible level (eg timeseries)
        2068             :
        2069             :     real( kind = core_rknd ), dimension(nz) :: &
        2070           0 :       T_in_K,      &  ! Absolute temperature         [K]
        2071           0 :       rsati,       &  ! Saturation w.r.t ice         [kg/kg]
        2072           0 :       shear,       &  ! Wind shear production term   [m^2/s^3]
        2073           0 :       chi,         &  ! Mellor's 's'                 [kg/kg]
        2074           0 :       chip2,         &  ! Variance of Mellor's 's'     [kg/kg]
        2075           0 :       rcm_in_cloud    ! rcm in cloud                 [kg/kg]
        2076             :
        2077             :     real( kind = core_rknd ) :: xtmp
        2078             :
        2079             :     ! ---- Begin Code ----
        2080             :
        2081             :     ! Sample fields
        2082             :
        2083           0 :     if ( stats_metadata%l_stats_samp ) then
        2084             :
        2085             :       ! stats_zt variables
        2086             :
        2087             :
        2088           0 :       if ( stats_metadata%iT_in_K > 0 .or. stats_metadata%irsati > 0 ) then
        2089           0 :         T_in_K = thlm2T_in_K( nz, thlm, exner, rcm )
        2090             :       else
        2091           0 :         T_in_K = -999._core_rknd
        2092             :       end if
        2093             :
        2094             :       call stat_update_var( stats_metadata%iT_in_K, T_in_K, & ! intent(in)
        2095           0 :                             stats_zt ) ! intent(inout)
        2096             :
        2097             :       call stat_update_var( stats_metadata%ithlm, thlm, & ! intent(in)
        2098           0 :                             stats_zt ) ! intent(inout)
        2099             :       call stat_update_var( stats_metadata%ithvm, thvm, & ! intent(in)
        2100           0 :                             stats_zt ) ! intent(inout)
        2101             :       call stat_update_var( stats_metadata%irtm, rtm, & ! intent(in)
        2102           0 :                             stats_zt ) ! intent(inout)
        2103             :       call stat_update_var( stats_metadata%ircm, rcm, & ! intent(in)
        2104           0 :                             stats_zt ) ! intent(inout)
        2105             :       call stat_update_var( stats_metadata%ium, um, & ! intent(in)
        2106           0 :                             stats_zt ) ! intent(inout)
        2107             :       call stat_update_var( stats_metadata%ivm, vm, & ! intent(in)
        2108           0 :                             stats_zt ) ! intent(inout)
        2109             :       call stat_update_var( stats_metadata%iwm_zt, wm_zt, & ! intent(in)
        2110           0 :                             stats_zt ) ! intent(inout)
        2111             :       call stat_update_var( stats_metadata%iwm_zm, wm_zm, & ! intent(in)
        2112           0 :                             stats_zm ) ! intent(inout)
        2113             :       call stat_update_var( stats_metadata%iug, ug, & ! intent(in)
        2114           0 :                             stats_zt ) ! intent(inout)
        2115             :       call stat_update_var( stats_metadata%ivg, vg, & ! intent(in)
        2116           0 :                             stats_zt ) ! intent(inout)
        2117             :       call stat_update_var( stats_metadata%icloud_frac, cloud_frac, & ! intent(in)
        2118           0 :                             stats_zt ) ! intent(inout)
        2119             :       call stat_update_var( stats_metadata%iice_supersat_frac, ice_supersat_frac, & ! intent(in)
        2120           0 :                             stats_zt) ! intent(inout)
        2121             :       call stat_update_var( stats_metadata%ircm_in_layer, rcm_in_layer, & ! intent(in)
        2122           0 :                             stats_zt ) ! intent(inout)
        2123             :       call stat_update_var( stats_metadata%icloud_cover, cloud_cover, & ! intent(in)
        2124           0 :                             stats_zt ) ! intent(inout)
        2125             :       call stat_update_var( stats_metadata%ircm_supersat_adj, rcm_supersat_adj, & ! intent(in)
        2126           0 :                             stats_zt ) ! intent(inout)
        2127             :       call stat_update_var( stats_metadata%ip_in_Pa, p_in_Pa, & ! intent(in)
        2128           0 :                             stats_zt ) ! intent(inout)
        2129             :       call stat_update_var( stats_metadata%iexner, exner, & ! intent(in)
        2130           0 :                             stats_zt ) ! intent(inout)
        2131             :       call stat_update_var( stats_metadata%irho_ds_zt, rho_ds_zt, & ! intent(in)
        2132           0 :                             stats_zt ) ! intent(inout)
        2133             :       call stat_update_var( stats_metadata%ithv_ds_zt, thv_ds_zt, & ! intent(in)
        2134           0 :                             stats_zt ) ! intent(inout)
        2135             :       call stat_update_var( stats_metadata%iLscale, Lscale, & ! intent(in)
        2136           0 :                             stats_zt ) ! intent(inout)
        2137             :       call stat_update_var( stats_metadata%iwpup2, wpup2, & ! intent(in)
        2138           0 :                             stats_zt ) ! intent(inout)
        2139             :       call stat_update_var( stats_metadata%iwpvp2, wpvp2, & ! intent(in)
        2140           0 :                             stats_zt ) ! intent(inout)
        2141             :       call stat_update_var( stats_metadata%iwp3, wp3, & ! intent(in)
        2142           0 :                             stats_zt ) ! intent(inout)
        2143             :       call stat_update_var( stats_metadata%iwpthlp2, wpthlp2, & ! intent(in)
        2144           0 :                             stats_zt ) ! intent(inout)
        2145             :       call stat_update_var( stats_metadata%iwp2thlp, wp2thlp, & ! intent(in)
        2146           0 :                             stats_zt ) ! intent(inout)
        2147             :       call stat_update_var( stats_metadata%iwprtp2, wprtp2, & ! intent(in)
        2148           0 :                             stats_zt ) ! intent(inout)
        2149             :       call stat_update_var( stats_metadata%iwp2rtp, wp2rtp, & ! intent(in)
        2150           0 :                             stats_zt ) ! intent(inout)
        2151             :       call stat_update_var( stats_metadata%iLscale_up, Lscale_up, & ! intent(in)
        2152           0 :                             stats_zt ) ! intent(inout)
        2153             :       call stat_update_var( stats_metadata%iLscale_down, Lscale_down, & ! intent(in)
        2154           0 :                             stats_zt ) ! intent(inout)
        2155             :       call stat_update_var( stats_metadata%itau_zt, tau_zt, & ! intent(in)
        2156           0 :                             stats_zt ) ! intent(inout)
        2157             :       call stat_update_var( stats_metadata%iKh_zt, Kh_zt, & ! intent(in)
        2158           0 :                             stats_zt ) ! intent(inout)
        2159             :       call stat_update_var( stats_metadata%iwp2thvp, wp2thvp, & ! intent(in)
        2160           0 :                             stats_zt ) ! intent(inout)
        2161             :       call stat_update_var( stats_metadata%iwp2rcp, wp2rcp, & ! intent(in)
        2162           0 :                             stats_zt ) ! intent(inout)
        2163             :       call stat_update_var( stats_metadata%iw_up_in_cloud, w_up_in_cloud, & ! intent(in)
        2164           0 :                             stats_zt ) ! intent(inout)
        2165             :       call stat_update_var( stats_metadata%iw_down_in_cloud, w_down_in_cloud, & ! intent(in)
        2166           0 :                             stats_zt ) ! intent(inout)
        2167             :       call stat_update_var( stats_metadata%icld_updr_frac, cloudy_updraft_frac, & ! intent(in)
        2168           0 :                             stats_zt ) ! intent(inout)
        2169             :       call stat_update_var( stats_metadata%icld_downdr_frac, cloudy_downdraft_frac, & ! intent(in)
        2170           0 :                             stats_zt ) ! intent(inout)
        2171             :       call stat_update_var( stats_metadata%iwprtpthlp, wprtpthlp, & ! intent(in)
        2172           0 :                             stats_zt ) ! intent(inout)
        2173             :       call stat_update_var( stats_metadata%irc_coef, rc_coef, & ! intent(in)
        2174           0 :                             stats_zt ) ! intent(inout)
        2175             :       call stat_update_var( stats_metadata%isigma_sqd_w_zt, sigma_sqd_w_zt, & ! intent(in)
        2176           0 :                             stats_zt ) ! intent(inout)
        2177             :       call stat_update_var( stats_metadata%irho, rho, & ! intent(in)
        2178           0 :                             stats_zt ) ! intent(inout)
        2179             :       call stat_update_var( stats_metadata%irsat, rsat, & ! intent(in)
        2180           0 :                             stats_zt ) ! intent(inout)
        2181           0 :       if ( stats_metadata%irsati > 0 ) then
        2182           0 :         do k = 1, nz
        2183           0 :           rsati(k) = sat_mixrat_ice( p_in_Pa(k), T_in_K(k), saturation_formula )
        2184             :         end do
        2185             :         call stat_update_var( stats_metadata%irsati, rsati, & ! intent(in)
        2186           0 :                               stats_zt ) ! intent(inout)
        2187             :       end if
        2188             :
        2189           0 :       call stat_update_var( stats_metadata%imixt_frac, pdf_params%mixt_frac(1,:), & ! intent(in)
        2190           0 :                             stats_zt ) ! intent(inout)
        2191           0 :       call stat_update_var( stats_metadata%iw_1, pdf_params%w_1(1,:), & ! intent(in)
        2192           0 :                             stats_zt ) ! intent(inout)
        2193           0 :       call stat_update_var( stats_metadata%iw_2, pdf_params%w_2(1,:), & ! intent(in)
        2194           0 :                             stats_zt ) ! intent(inout)
        2195           0 :       call stat_update_var( stats_metadata%ivarnce_w_1, pdf_params%varnce_w_1(1,:), & ! intent(in)
        2196           0 :                             stats_zt ) ! intent(inout)
        2197           0 :       call stat_update_var( stats_metadata%ivarnce_w_2, pdf_params%varnce_w_2(1,:), & ! intent(in)
        2198           0 :                             stats_zt ) ! intent(inout)
        2199           0 :       call stat_update_var( stats_metadata%ithl_1, pdf_params%thl_1(1,:), & ! intent(in)
        2200           0 :                             stats_zt ) ! intent(inout)
        2201           0 :       call stat_update_var( stats_metadata%ithl_2, pdf_params%thl_2(1,:), & ! intent(in)
        2202           0 :                             stats_zt ) ! intent(inout)
        2203           0 :       call stat_update_var( stats_metadata%ivarnce_thl_1, pdf_params%varnce_thl_1(1,:), & ! intent(in)
        2204           0 :                             stats_zt ) ! intent(inout)
        2205           0 :       call stat_update_var( stats_metadata%ivarnce_thl_2, pdf_params%varnce_thl_2(1,:), & ! intent(in)
        2206           0 :                             stats_zt ) ! intent(inout)
        2207           0 :       call stat_update_var( stats_metadata%irt_1, pdf_params%rt_1(1,:), & ! intent(in)
        2208           0 :                             stats_zt ) ! intent(inout)
        2209           0 :       call stat_update_var( stats_metadata%irt_2, pdf_params%rt_2(1,:), & ! intent(in)
        2210           0 :                             stats_zt ) ! intent(inout)
        2211           0 :       call stat_update_var( stats_metadata%ivarnce_rt_1, pdf_params%varnce_rt_1(1,:), & ! intent(in)
        2212           0 :                             stats_zt ) ! intent(inout)
        2213           0 :       call stat_update_var( stats_metadata%ivarnce_rt_2, pdf_params%varnce_rt_2(1,:), & ! intent(in)
        2214           0 :                             stats_zt ) ! intent(inout )
        2215           0 :       call stat_update_var( stats_metadata%irc_1, pdf_params%rc_1(1,:), & ! intent(in)
        2216           0 :                             stats_zt ) ! intent(inout)
        2217           0 :       call stat_update_var( stats_metadata%irc_2, pdf_params%rc_2(1,:), & ! intent(in)
        2218           0 :                             stats_zt ) ! intent(inout)
        2219           0 :       call stat_update_var( stats_metadata%irsatl_1, pdf_params%rsatl_1(1,:), & ! intent(in)
        2220           0 :                             stats_zt ) ! intent(inout)
        2221           0 :       call stat_update_var( stats_metadata%irsatl_2, pdf_params%rsatl_2(1,:), & ! intent(in)
        2222           0 :                             stats_zt ) ! intent(inout)
        2223           0 :       call stat_update_var( stats_metadata%icloud_frac_1, pdf_params%cloud_frac_1(1,:), & ! intent(in)
        2224           0 :                             stats_zt ) ! intent(inout)
        2225           0 :       call stat_update_var( stats_metadata%icloud_frac_2, pdf_params%cloud_frac_2(1,:), & ! intent(in)
        2226           0 :                             stats_zt ) ! intent(inout)
        2227           0 :       call stat_update_var( stats_metadata%ichi_1, pdf_params%chi_1(1,:), & ! intent(in)
        2228           0 :                             stats_zt ) ! intent(inout)
        2229           0 :       call stat_update_var( stats_metadata%ichi_2, pdf_params%chi_2(1,:), & ! intent(in)
        2230           0 :                             stats_zt ) ! intent(inout)
        2231           0 :       call stat_update_var( stats_metadata%istdev_chi_1, pdf_params%stdev_chi_1(1,:), & ! intent(in)
        2232           0 :                             stats_zt ) ! intent(inout)
        2233           0 :       call stat_update_var( stats_metadata%istdev_chi_2, pdf_params%stdev_chi_2(1,:), & ! intent(in)
        2234           0 :                             stats_zt ) ! intent(inout)
        2235           0 :       call stat_update_var( stats_metadata%istdev_eta_1, pdf_params%stdev_eta_1(1,:), & ! intent(in)
        2236           0 :                             stats_zt ) ! intent(inout)
        2237           0 :       call stat_update_var( stats_metadata%istdev_eta_2, pdf_params%stdev_eta_2(1,:), & ! intent(in)
        2238           0 :                             stats_zt ) ! intent(inout)
        2239           0 :       call stat_update_var( stats_metadata%icovar_chi_eta_1, pdf_params%covar_chi_eta_1(1,:), & ! intent(in)
        2240           0 :                             stats_zt ) ! intent(inout)
        2241           0 :       call stat_update_var( stats_metadata%icovar_chi_eta_2, pdf_params%covar_chi_eta_2(1,:), & ! intent(in)
        2242           0 :                             stats_zt ) ! intent(inout)
        2243           0 :       call stat_update_var( stats_metadata%icorr_w_chi_1, pdf_params%corr_w_chi_1(1,:), & ! intent(in)
        2244           0 :                             stats_zt ) ! intent(inout)
        2245           0 :       call stat_update_var( stats_metadata%icorr_w_chi_2, pdf_params%corr_w_chi_2(1,:), & ! intent(in)
        2246           0 :                             stats_zt ) ! intent(inout)
        2247           0 :       call stat_update_var( stats_metadata%icorr_w_eta_1, pdf_params%corr_w_eta_1(1,:), & ! intent(in)
        2248           0 :                             stats_zt ) ! intent(inout)
        2249           0 :       call stat_update_var( stats_metadata%icorr_w_eta_2, pdf_params%corr_w_eta_2(1,:), & ! intent(in)
        2250           0 :                             stats_zt ) ! intent(inout)
        2251           0 :       call stat_update_var( stats_metadata%icorr_chi_eta_1, pdf_params%corr_chi_eta_1(1,:), & ! intent(in)
        2252           0 :                             stats_zt ) ! intent(inout)
        2253           0 :       call stat_update_var( stats_metadata%icorr_chi_eta_2, pdf_params%corr_chi_eta_2(1,:), & ! intent(in)
        2254           0 :                             stats_zt ) ! intent(inout)
        2255           0 :       call stat_update_var( stats_metadata%icorr_w_rt_1, pdf_params%corr_w_rt_1(1,:), & ! intent(in)
        2256           0 :                             stats_zt ) ! intent(inout)
        2257           0 :       call stat_update_var( stats_metadata%icorr_w_rt_2, pdf_params%corr_w_rt_2(1,:), & ! intent(in)
        2258           0 :                             stats_zt ) ! intent(inout)
        2259           0 :       call stat_update_var( stats_metadata%icorr_w_thl_1, pdf_params%corr_w_thl_1(1,:), & ! intent(in)
        2260           0 :                             stats_zt ) ! intent(inout)
        2261           0 :       call stat_update_var( stats_metadata%icorr_w_thl_2, pdf_params%corr_w_thl_2(1,:), & ! intent(in)
        2262           0 :                             stats_zt ) ! intent(inout)
        2263           0 :       call stat_update_var( stats_metadata%icorr_rt_thl_1, pdf_params%corr_rt_thl_1(1,:), & ! intent(in)
        2264           0 :                             stats_zt ) ! intent(inout)
        2265           0 :       call stat_update_var( stats_metadata%icorr_rt_thl_2, pdf_params%corr_rt_thl_2(1,:), & ! intent(in)
        2266           0 :                             stats_zt ) ! intent(inout)
        2267           0 :       call stat_update_var( stats_metadata%icrt_1, pdf_params%crt_1(1,:), & ! intent(in)
        2268           0 :                             stats_zt ) ! intent(inout)
        2269           0 :       call stat_update_var( stats_metadata%icrt_2, pdf_params%crt_2(1,:), & ! intent(in)
        2270           0 :                             stats_zt ) ! intent(inout)
        2271           0 :       call stat_update_var( stats_metadata%icthl_1, pdf_params%cthl_1(1,:), & ! intent(in)
        2272           0 :                             stats_zt ) ! intent(inout)
        2273           0 :       call stat_update_var( stats_metadata%icthl_2, pdf_params%cthl_2(1,:), & ! intent(in)
        2274           0 :                             stats_zt ) ! intent(inout)
        2275             :       call stat_update_var( stats_metadata%iwp2_zt, wp2_zt, & ! intent(in)
        2276           0 :                             stats_zt ) ! intent(inout)
        2277             :       call stat_update_var( stats_metadata%ithlp2_zt, thlp2_zt, & ! intent(in)
        2278           0 :                             stats_zt ) ! intent(inout)
        2279             :       call stat_update_var( stats_metadata%ithlp3, thlp3, & ! intent(in)
        2280           0 :                             stats_zt ) ! intent(inout)
        2281             :       call stat_update_var( stats_metadata%iwpthlp_zt, wpthlp_zt, & ! intent(in)
        2282           0 :                             stats_zt ) ! intent(inout)
        2283             :       call stat_update_var( stats_metadata%iwprtp_zt, wprtp_zt, & ! intent(in)
        2284           0 :                             stats_zt ) ! intent(inout)
        2285             :       call stat_update_var( stats_metadata%irtp2_zt, rtp2_zt, & ! intent(in)
        2286           0 :                             stats_zt ) ! intent(inout)
        2287             :       call stat_update_var( stats_metadata%irtp3, rtp3, & ! intent(in)
        2288           0 :                             stats_zt ) ! intent(inout)
        2289             :       call stat_update_var( stats_metadata%irtpthlp_zt, rtpthlp_zt, & ! intent(in)
        2290           0 :                             stats_zt ) ! intent(inout)
        2291             :       call stat_update_var( stats_metadata%iup2_zt, up2_zt, & ! intent(in)
        2292           0 :                             stats_zt ) ! intent(inout)
        2293             :       call stat_update_var( stats_metadata%ivp2_zt, vp2_zt, & ! intent(in)
        2294           0 :                             stats_zt ) ! intent(inout)
        2295             :       call stat_update_var( stats_metadata%iupwp_zt, upwp_zt, & ! intent(in)
        2296           0 :                             stats_zt ) ! intent(inout)
        2297             :       call stat_update_var( stats_metadata%ivpwp_zt, vpwp_zt, & ! intent(in)
        2298           0 :                             stats_zt ) ! intent(inout)
        2299             :       call stat_update_var( stats_metadata%ia3_coef_zt, a3_coef_zt, & ! intent(in)
        2300           0 :                             stats_zt ) ! intent(inout)
        2301             :       call stat_update_var( stats_metadata%iwp3_on_wp2_zt, wp3_on_wp2_zt, & ! intent(in)
        2302           0 :                             stats_zt ) ! intent(inout)
        2303             :
        2304           0 :       if ( stats_metadata%ichi > 0 ) then
        2305             :         ! Determine 's' from Mellor (1977) (extended liquid water)
        2306           0 :         chi(:) = pdf_params%mixt_frac(1,:) * pdf_params%chi_1(1,:) &
        2307           0 :                     + (1.0_core_rknd-pdf_params%mixt_frac(1,:)) * pdf_params%chi_2(1,:)
        2308             :         call stat_update_var( stats_metadata%ichi, chi, & ! intent(in)
        2309           0 :                              stats_zt ) ! intent(inout)
        2310             :       end if
        2311             :
        2312             :       ! Calculate variance of chi
        2313           0 :       if ( stats_metadata%ichip2 > 0 ) then
        2314           0 :         chip2 = compute_variance_binormal( chi, pdf_params%chi_1(1,:), pdf_params%chi_2(1,:), &
        2315           0 :                                          pdf_params%stdev_chi_1(1,:), pdf_params%stdev_chi_2(1,:), &
        2316           0 :                                          pdf_params%mixt_frac(1,:) )
        2317             :         call stat_update_var( stats_metadata%ichip2, chip2, & ! intent(in)
        2318           0 :                               stats_zt ) ! intent(inout)
        2319             :       end if
        2320             :
        2321           0 :       if ( sclr_dim > 0 ) then
        2322           0 :         do isclr=1, sclr_dim
        2323           0 :           call stat_update_var( stats_metadata%isclrm(isclr), sclrm(:,isclr), & ! intent(in)
        2324           0 :                                 stats_zt ) ! intent(inout)
        2325           0 :           call stat_update_var( stats_metadata%isclrm_f(isclr), sclrm_forcing(:,isclr),  & ! intent(in)
        2326           0 :                                 stats_zt ) ! intent(inout)
        2327             :         end do
        2328             :       end if
        2329             :
        2330           0 :       if ( edsclr_dim > 0 ) then
        2331           0 :         do isclr = 1, edsclr_dim
        2332           0 :           call stat_update_var( stats_metadata%iedsclrm(isclr), edsclrm(:,isclr), & ! intent(in)
        2333           0 :                                 stats_zt ) ! intent(inout)
        2334           0 :           call stat_update_var( stats_metadata%iedsclrm_f(isclr), edsclrm_forcing(:,isclr), & ! intent(in)
        2335           0 :                                 stats_zt ) ! intent(inout)
        2336             :         end do
        2337             :       end if
        2338             :
        2339             :       ! Calculate rcm in cloud
        2340           0 :       if ( stats_metadata%ircm_in_cloud > 0 ) then
        2341           0 :         where ( cloud_frac(:) > cloud_frac_min )
        2342             :             rcm_in_cloud(:) = rcm / cloud_frac
        2343             :         elsewhere
        2344             :             rcm_in_cloud(:) = rcm
        2345             :         endwhere
        2346             :
        2347             :         call stat_update_var( stats_metadata%ircm_in_cloud, rcm_in_cloud, & ! intent(in)
        2348           0 :                               stats_zt ) ! intent(inout)
        2349             :       end if
        2350             :
        2351             :       ! stats_zm variables
        2352             :
        2353             :       call stat_update_var( stats_metadata%iwp2, wp2, & ! intent(in)
        2354           0 :                             stats_zm ) ! intent(inout)
        2355             :       call stat_update_var( stats_metadata%iwp3_zm, wp3_zm, & ! intent(in)
        2356           0 :                             stats_zm ) ! intent(inout)
        2357             :       call stat_update_var( stats_metadata%irtp2, rtp2, & ! intent(in)
        2358           0 :                             stats_zm ) ! intent(inout)
        2359             :       call stat_update_var( stats_metadata%ithlp2, thlp2, & ! intent(in)
        2360           0 :                             stats_zm ) ! intent(inout)
        2361             :       call stat_update_var( stats_metadata%irtpthlp, rtpthlp, & ! intent(in)
        2362           0 :                             stats_zm ) ! intent(inout)
        2363             :       call stat_update_var( stats_metadata%iwprtp, wprtp, & ! intent(in)
        2364           0 :                             stats_zm ) ! intent(inout)
        2365             :       call stat_update_var( stats_metadata%iwpthlp, wpthlp, & ! intent(in)
        2366           0 :                             stats_zm ) ! intent(inout)
        2367             :       call stat_update_var( stats_metadata%iwp2up2, wp2up2, & ! intent(in)
        2368           0 :                             stats_zm ) ! intent(inout)
        2369             :       call stat_update_var( stats_metadata%iwp2vp2, wp2vp2, &  ! intent(in)
        2370           0 :                             stats_zm ) ! intent(inout)
        2371             :       call stat_update_var( stats_metadata%iwp4, wp4, & ! intent(in)
        2372           0 :                             stats_zm ) ! intent(inout)
        2373             :       call stat_update_var( stats_metadata%iwpthvp, wpthvp, & ! intent(in)
        2374           0 :                             stats_zm ) ! intent(inout)
        2375             :       call stat_update_var( stats_metadata%irtpthvp, rtpthvp, & ! intent(in)
        2376           0 :                             stats_zm ) ! intent(inout)
        2377             :       call stat_update_var( stats_metadata%ithlpthvp, thlpthvp, & ! intent(in)
        2378           0 :                             stats_zm ) ! intent(inout)
        2379             :       call stat_update_var( stats_metadata%itau_zm, tau_zm, & ! intent(in)
        2380           0 :                             stats_zm ) ! intent(inout)
        2381             :       call stat_update_var( stats_metadata%iKh_zm, Kh_zm, & ! intent(in)
        2382           0 :                             stats_zm ) ! intent(inout)
        2383             :       call stat_update_var( stats_metadata%iwprcp, wprcp, & ! intent(in)
        2384           0 :                             stats_zm ) ! intent(inout)
        2385             :       call stat_update_var( stats_metadata%irc_coef_zm, rc_coef_zm, & ! intent(in)
        2386           0 :                             stats_zm ) ! intent(inout)
        2387             :       call stat_update_var( stats_metadata%ithlprcp, thlprcp, & ! intent(in)
        2388           0 :                             stats_zm ) ! intent(inout)
        2389             :       call stat_update_var( stats_metadata%irtprcp, rtprcp, & ! intent(in)
        2390           0 :                             stats_zm ) ! intent(inout)
        2391             :       call stat_update_var( stats_metadata%ircp2, rcp2, & ! intent(in)
        2392           0 :                             stats_zm ) ! intent(inout)
        2393             :       call stat_update_var( stats_metadata%iupwp, upwp, & ! intent(in)
        2394           0 :                             stats_zm ) ! intent(inout)
        2395             :       call stat_update_var( stats_metadata%ivpwp, vpwp, & ! intent(in)
        2396           0 :                             stats_zm ) ! intent(inout)
        2397             :       call stat_update_var( stats_metadata%ivp2, vp2, & ! intent(in)
        2398           0 :                             stats_zm ) ! intent(inout)
        2399             :       call stat_update_var( stats_metadata%iup2, up2, & ! intent(in)
        2400           0 :                             stats_zm ) ! intent(inout)
        2401             :       call stat_update_var( stats_metadata%irho_zm, rho_zm, & ! intent(in)
        2402           0 :                             stats_zm ) ! intent(inout)
        2403             :       call stat_update_var( stats_metadata%isigma_sqd_w, sigma_sqd_w, & ! intent(in)
        2404           0 :                             stats_zm ) ! intent(inout)
        2405             :       call stat_update_var( stats_metadata%irho_ds_zm, rho_ds_zm, & ! intent(in)
        2406           0 :                             stats_zm ) ! intent(inout)
        2407             :       call stat_update_var( stats_metadata%ithv_ds_zm, thv_ds_zm, & ! intent(in)
        2408           0 :                             stats_zm ) ! intent(inout)
        2409             :       call stat_update_var( stats_metadata%iem, em, & ! intent(in)
        2410           0 :                             stats_zm ) ! intent(inout)
        2411             :       call stat_update_var( stats_metadata%iSkw_velocity, Skw_velocity, & ! intent(in)
        2412           0 :                             stats_zm ) ! intent(inout)
        2413             :       call stat_update_var( stats_metadata%ia3_coef, a3_coef, & ! intent(in)
        2414           0 :                             stats_zm ) ! intent(inout)
        2415             :       call stat_update_var( stats_metadata%iwp3_on_wp2, wp3_on_wp2, & ! intent(in)
        2416           0 :                             stats_zm ) ! intent(inout)
        2417             :       call stat_update_var( stats_metadata%iwp3_on_wp2_cfl_num, wp3_on_wp2 * dt / dzm, & ! intent(in)
        2418           0 :                             stats_zm ) ! intent(inout)
        2419             :
        2420             :       call stat_update_var( stats_metadata%icloud_frac_zm, cloud_frac_zm, & ! intent(in)
        2421           0 :                             stats_zm ) ! intent(inout)
        2422             :       call stat_update_var( stats_metadata%iice_supersat_frac_zm, ice_supersat_frac_zm, & ! intent(in)
        2423           0 :                             stats_zm ) ! intent(inout)
        2424             :       call stat_update_var( stats_metadata%ircm_zm, rcm_zm, & ! intent(in)
        2425           0 :                             stats_zm ) ! intent(inout)
        2426             :       call stat_update_var( stats_metadata%irtm_zm, rtm_zm, & ! intent(in)
        2427           0 :                             stats_zm ) ! intent(inout)
        2428             :       call stat_update_var( stats_metadata%ithlm_zm, thlm_zm, & ! intent(in)
        2429           0 :                             stats_zm ) ! intent(inout)
        2430           0 :       call stat_update_var( stats_metadata%iw_1_zm, pdf_params_zm%w_1(1,:), & ! intent(in)
        2431           0 :                             stats_zm ) ! intent(inout)
        2432           0 :       call stat_update_var( stats_metadata%iw_2_zm, pdf_params_zm%w_2(1,:), & ! intent(in)
        2433           0 :                             stats_zm ) ! intent(inout)
        2434           0 :       call stat_update_var( stats_metadata%ivarnce_w_1_zm, pdf_params_zm%varnce_w_1(1,:), & ! intent(in)
        2435           0 :                             stats_zm ) ! intent(inout)
        2436           0 :       call stat_update_var( stats_metadata%ivarnce_w_2_zm, pdf_params_zm%varnce_w_2(1,:), & ! intent(in)
        2437           0 :                             stats_zm ) ! intent(inout)
        2438           0 :       call stat_update_var( stats_metadata%imixt_frac_zm, pdf_params_zm%mixt_frac(1,:), & ! intent(in)
        2439           0 :                             stats_zm ) ! intent(inout)
        2440             :
        2441           0 :       if ( sclr_dim > 0 ) then
        2442           0 :         do isclr=1, sclr_dim
        2443           0 :           call stat_update_var( stats_metadata%isclrp2(isclr), sclrp2(:,isclr), & ! intent(in)
        2444           0 :                                 stats_zm ) ! intent(inout)
        2445           0 :           call stat_update_var( stats_metadata%isclrprtp(isclr), sclrprtp(:,isclr), & ! intent(in)
        2446           0 :                                 stats_zm ) ! intent(inout)
        2447           0 :           call stat_update_var( stats_metadata%isclrpthvp(isclr), sclrpthvp(:,isclr), & ! intent(in)
        2448           0 :                                 stats_zm ) ! intent(inout)
        2449           0 :           call stat_update_var( stats_metadata%isclrpthlp(isclr), sclrpthlp(:,isclr), & ! intent(in)
        2450           0 :                                  stats_zm ) ! intent(inout)
        2451           0 :           call stat_update_var( stats_metadata%isclrprcp(isclr), sclrprcp(:,isclr), & ! intent(in)
        2452           0 :                                 stats_zm ) ! intent(inout)
        2453           0 :           call stat_update_var( stats_metadata%iwpsclrp(isclr), wpsclrp(:,isclr), & ! intent(in)
        2454           0 :                                stats_zm ) ! intent(inout)
        2455           0 :           call stat_update_var( stats_metadata%iwp2sclrp(isclr), wp2sclrp(:,isclr), & ! intent(in)
        2456           0 :                                 stats_zm ) ! intent(inout)
        2457           0 :           call stat_update_var( stats_metadata%iwpsclrp2(isclr), wpsclrp2(:,isclr), & ! intent(in)
        2458           0 :                                 stats_zm ) ! intent(inout)
        2459           0 :           call stat_update_var( stats_metadata%iwpsclrprtp(isclr), wpsclrprtp(:,isclr), & ! intent(in)
        2460           0 :                                 stats_zm ) ! intent(inout)
        2461           0 :           call stat_update_var( stats_metadata%iwpsclrpthlp(isclr), wpsclrpthlp(:,isclr), & ! intent(in)
        2462           0 :                                 stats_zm ) ! intent(inout)
        2463             :         end do
        2464             :       end if
        2465           0 :       if ( edsclr_dim > 0 ) then
        2466           0 :         do isclr = 1, edsclr_dim
        2467           0 :           call stat_update_var( stats_metadata%iwpedsclrp(isclr), wpedsclrp(:,isclr), & ! intent(in)
        2468           0 :                                 stats_zm ) ! intent(inout)
        2469             :         end do
        2470             :       end if
        2471             :
        2472             :       ! Calculate shear production
        2473           0 :       if ( stats_metadata%ishear > 0 ) then
        2474           0 :         do k = 1, nz-1, 1
        2475           0 :           shear(k) = - upwp(k) * ( um(k+1) - um(k) ) * invrs_dzm(k)  &
        2476           0 :                      - vpwp(k) * ( vm(k+1) - vm(k) ) * invrs_dzm(k)
        2477             :         enddo
        2478           0 :         shear(nz) = 0.0_core_rknd
        2479             :       end if
        2480             :       call stat_update_var( stats_metadata%ishear, shear, & ! intent(in)
        2481           0 :                             stats_zm ) ! intent(inout)
        2482             :
        2483             :       ! stats_sfc variables
        2484             :
        2485             :       ! Cloud cover
        2486           0 :       call stat_update_var_pt( stats_metadata%icc, grid_level, maxval( cloud_frac(1:nz) ), & ! intent(in)
        2487           0 :                                stats_sfc ) ! intent(inout)
        2488             :
        2489             :       ! Cloud base
        2490           0 :       if ( stats_metadata%iz_cloud_base > 0 ) then
        2491             :
        2492             :         k = 1
        2493           0 :         do while ( rcm(k) < rc_tol .and. k < nz )
        2494           0 :           k = k + 1
        2495             :         enddo
        2496             :
        2497           0 :         if ( k > 1 .and. k < nz) then
        2498             :
        2499             :           ! Use linear interpolation to find the exact height of the
        2500             :           ! rc_tol kg/kg level.  Brian.
        2501             :           call stat_update_var_pt( stats_metadata%iz_cloud_base, grid_level, & ! intent(in)
        2502             :                                    lin_interpolate_two_points( rc_tol, rcm(k), & ! intent(in)
        2503           0 :                                    rcm(k-1), zt(k), zt(k-1) ), & ! intent(in)
        2504           0 :                                    stats_sfc ) ! intent(inout)
        2505             :
        2506             :         else
        2507             :
        2508             :           ! Set the cloud base output to -10m, if it's clear.
        2509             :           ! Known magic number
        2510             :           call stat_update_var_pt( stats_metadata%iz_cloud_base, grid_level, -10.0_core_rknd , & ! intent(in)
        2511           0 :                                    stats_sfc ) ! intent(inout)
        2512             :
        2513             :         end if ! k > 1 and k < nz
        2514             :
        2515             :       end if ! stats_metadata%iz_cloud_base > 0
        2516             :
        2517             :       ! Liquid Water Path
        2518           0 :       if ( stats_metadata%ilwp > 0 ) then
        2519             :
        2520             :         xtmp &
        2521             :         = vertical_integral &
        2522           0 :                ( (nz - 2 + 1), rho_ds_zt(2:nz), &
        2523           0 :                  rcm(2:nz), dzt(2:nz) )
        2524             :
        2525             :         call stat_update_var_pt( stats_metadata%ilwp, grid_level, xtmp, & ! intent(in)
        2526           0 :                                  stats_sfc ) ! intent(inout)
        2527             :
        2528             :       end if
        2529             :
        2530             :       ! Vapor Water Path (Precipitable Water)
        2531           0 :       if ( stats_metadata%ivwp > 0 ) then
        2532             :
        2533             :         xtmp &
        2534             :         = vertical_integral &
        2535           0 :                ( (nz - 2 + 1), rho_ds_zt(2:nz), &
        2536           0 :                  ( rtm(2:nz) - rcm(2:nz) ), dzt(2:nz) )
        2537             :
        2538             :         call stat_update_var_pt( stats_metadata%ivwp, grid_level, xtmp, & ! intent(in)
        2539           0 :                                  stats_sfc ) ! intent(inout)
        2540             :
        2541             :       end if
        2542             :
        2543             :
        2544             :       ! Vertical average of thermodynamic level variables.
        2545             :
        2546             :       ! Find the vertical average of thermodynamic level variables, averaged from
        2547             :       ! level 2 (the first thermodynamic level above model surface) through
        2548             :       ! level nz (the top of the model).  Use the vertical averaging function
        2549             :       ! found in advance_helper_module.F90.
        2550             :
        2551             :       ! Vertical average of thlm.
        2552             :       call stat_update_var_pt( stats_metadata%ithlm_vert_avg, grid_level,  & ! intent(in)
        2553           0 :            vertical_avg( (nz-2+1), rho_ds_zt(2:nz), & ! intent(in)
        2554             :                          thlm(2:nz), dzt(2:nz) ), & ! intent(in)
        2555           0 :                                stats_sfc ) ! intent(inout)
        2556             :
        2557             :       ! Vertical average of rtm.
        2558             :       call stat_update_var_pt( stats_metadata%irtm_vert_avg, grid_level,  & ! intent(in)
        2559             :            vertical_avg( (nz-2+1), rho_ds_zt(2:nz), & ! intent(in)
        2560             :                          rtm(2:nz), dzt(2:nz) ), & ! intent(in)
        2561           0 :                                stats_sfc ) ! intent(inout)
        2562             :
        2563             :       ! Vertical average of um.
        2564             :       call stat_update_var_pt( stats_metadata%ium_vert_avg, grid_level,  & ! intent(in)
        2565             :            vertical_avg( (nz-2+1), rho_ds_zt(2:nz), & ! intent(in)
        2566             :                          um(2:nz), dzt(2:nz) ), & ! intent(in)
        2567           0 :                                stats_sfc ) ! intent(inout)
        2568             :
        2569             :       ! Vertical average of vm.
        2570             :       call stat_update_var_pt( stats_metadata%ivm_vert_avg, grid_level,  & ! intent(in)
        2571             :            vertical_avg( (nz-2+1), rho_ds_zt(2:nz), & ! intent(in)
        2572             :                          vm(2:nz), dzt(2:nz) ), & ! intent(in)
        2573           0 :                                stats_sfc ) ! intent(inout)
        2574             :
        2575             :       ! Vertical average of momentum level variables.
        2576             :
        2577             :       ! Find the vertical average of momentum level variables, averaged over the
        2578             :       ! entire vertical profile (level 1 through level nz).  Use the vertical
        2579             :       ! averaging function found in advance_helper_module.F90.
        2580             :
        2581             :       ! Vertical average of wp2.
        2582             :       call stat_update_var_pt( stats_metadata%iwp2_vert_avg, grid_level,  & ! intent(in)
        2583             :            vertical_avg( (nz-1+1), rho_ds_zm(1:nz), & ! intent(in)
        2584             :                          wp2(1:nz), dzm(1:nz) ), & ! intent(in)
        2585           0 :                                stats_sfc ) ! intent(inout)
        2586             :
        2587             :       ! Vertical average of up2.
        2588             :       call stat_update_var_pt( stats_metadata%iup2_vert_avg, grid_level,  & ! intent(in)
        2589             :            vertical_avg( (nz-1+1), rho_ds_zm(1:nz), & ! intent(in)
        2590             :                          up2(1:nz), dzm(1:nz) ), & ! intent(in)
        2591           0 :                                stats_sfc ) ! intent(inout)
        2592             :
        2593             :       ! Vertical average of vp2.
        2594             :       call stat_update_var_pt( stats_metadata%ivp2_vert_avg, grid_level,  & ! intent(in)
        2595             :            vertical_avg( (nz-1+1), rho_ds_zm(1:nz), & ! intent(in)
        2596             :                          vp2(1:nz), dzm(1:nz) ), & ! intent(in)
        2597           0 :                                stats_sfc ) ! intent(inout)
        2598             :
        2599             :       ! Vertical average of rtp2.
        2600             :       call stat_update_var_pt( stats_metadata%irtp2_vert_avg, grid_level,  & ! intent(in)
        2601             :            vertical_avg( (nz-1+1), rho_ds_zm(1:nz), & ! intent(in)
        2602             :                          rtp2(1:nz), dzm(1:nz) ), & ! intent(in)
        2603           0 :                                stats_sfc ) ! intent(inout)
        2604             :
        2605             :       ! Vertical average of thlp2.
        2606             :       call stat_update_var_pt( stats_metadata%ithlp2_vert_avg, grid_level,  & ! intent(in)
        2607             :            vertical_avg( (nz-1+1), rho_ds_zm(1:nz), & ! intent(in)
        2608             :                          thlp2(1:nz), dzm(1:nz) ), & ! intent(in)
        2609           0 :                                stats_sfc ) ! intent(inout)
        2610             :
        2611             :
        2612           0 :       if (stats_metadata%itot_vartn_normlzd_rtm > 0) then
        2613           0 :         if (abs(rtm(nz) - rtm(1)) < eps) then
        2614           0 :           write(fstderr, *) "Warning: tot_vartn_normlzd_rtm tried to divide by zero denominator ", &
        2615           0 :                             "(surface level value was equal to top level value)"
        2616           0 :           xtmp = -999_core_rknd  ! workaround to signify zero denominator
        2617             :         else
        2618           0 :           xtmp = sum(abs(rtm(2 : nz) - rtm(1 : nz-1)) / abs(rtm(nz) - rtm(1)))
        2619             :         end if
        2620             :
        2621             :         call stat_update_var_pt( stats_metadata%itot_vartn_normlzd_rtm, grid_level, xtmp, & ! intent(in)
        2622           0 :                                  stats_sfc ) ! intent(inout)
        2623             :       end if
        2624             :
        2625           0 :       if (stats_metadata%itot_vartn_normlzd_thlm > 0) then
        2626           0 :         if (abs(thlm(nz) - thlm(1)) < eps) then
        2627           0 :           write(fstderr, *) "Warning: tot_vartn_normlzd_thlm tried to divide by zero denominator ", &
        2628           0 :                             "(surface level value was equal to top level value)"
        2629           0 :           xtmp = -999_core_rknd  ! workaround to signify zero denominator
        2630             :         else
        2631           0 :           xtmp = sum(abs(thlm(2 : nz) - thlm(1 : nz-1)) / abs(thlm(nz) - thlm(1)))
        2632             :         end if
        2633             :
        2634             :         call stat_update_var_pt( stats_metadata%itot_vartn_normlzd_thlm, grid_level, xtmp, & ! intent(in)
        2635           0 :                                  stats_sfc ) ! intent(inout)
        2636             :       end if
        2637             :
        2638           0 :       if (stats_metadata%itot_vartn_normlzd_wprtp > 0) then
        2639           0 :         if (abs(wprtp(nz) - wprtp(1)) < eps) then
        2640           0 :           write(fstderr, *) "Warning: tot_vartn_normlzd_wprtp tried to divide by zero denominator ", &
        2641           0 :                             "(surface level value was equal to top level value)"
        2642           0 :           xtmp = -999_core_rknd  ! workaround to signify zero denominator
        2643             :         else
        2644           0 :           xtmp = sum(abs(wprtp(2 : nz) - wprtp(1 : nz-1)) / abs(wprtp(nz) - wprtp(1)))
        2645             :         end if
        2646             :
        2647             :         call stat_update_var_pt( stats_metadata%itot_vartn_normlzd_wprtp, grid_level, xtmp, & ! intent(in)
        2648           0 :                                  stats_sfc ) ! intent(inout)
        2649             :       end if
        2650             :     end if ! stats_metadata%l_stats_samp
        2651             :
        2652             :
        2653           0 :     return
        2654             :   end subroutine stats_accumulate
        2655             : !------------------------------------------------------------------------------
        2656           0 :   subroutine stats_accumulate_hydromet( gr, hydromet_dim, hm_metadata, & ! intent(in)
        2657           0 :                                         hydromet, rho_ds_zt,              & ! intent(in)
        2658             :                                         stats_metadata,                   & ! intent(in)
        2659             :                                         stats_zt, stats_sfc )               ! intent(inout)
        2660             : ! Description:
        2661             : !   Compute stats related the hydrometeors
        2662             :
        2663             : ! References:
        2664             : !   None
        2665             : !------------------------------------------------------------------------------
        2666             :
        2667             :     use grid_class, only: &
        2668             :         grid ! Type
        2669             :
        2670             :     use corr_varnce_module, only: &
        2671             :         hm_metadata_type
        2672             :
        2673             :     use advance_helper_module, only: &
        2674             :         vertical_integral ! Procedure(s)
        2675             :
        2676             :     use stats_type_utilities, only: &
        2677             :         stat_update_var, & ! Procedure(s)
        2678             :         stat_update_var_pt
        2679             :
        2680             :     use clubb_precision, only: &
        2681             :         core_rknd ! Variable(s)
        2682             :
        2683             :     use stats_type, only: &
        2684             :         stats ! Type
        2685             :
        2686             :     use stats_variables, only: &
        2687             :         stats_metadata_type
        2688             :
        2689             :     implicit none
        2690             :
        2691             :     type (grid), target, intent(in) :: &
        2692             :       gr
        2693             :
        2694             :     integer, intent(in) :: &
        2695             :       hydromet_dim
        2696             :
        2697             :     type (hm_metadata_type), intent(in) :: &
        2698             :       hm_metadata
        2699             :
        2700             :     type (stats_metadata_type), intent(in) :: &
        2701             :       stats_metadata
        2702             :
        2703             :     type (stats), target, intent(inout) :: &
        2704             :       stats_zt, &
        2705             :       stats_sfc
        2706             :
        2707             :     ! Input Variables
        2708             :     real( kind = core_rknd ), dimension(gr%nz,hydromet_dim), intent(in) :: &
        2709             :       hydromet ! All hydrometeors except for rcm        [units vary]
        2710             :
        2711             :     real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
        2712             :       rho_ds_zt ! Dry, static density (thermo. levs.)      [kg/m^3]
        2713             :
        2714             :     ! Local Variables
        2715             :     real(kind=core_rknd) :: xtmp
        2716             :
        2717             :     integer :: grid_level = 1
        2718             :
        2719             :     ! ---- Begin Code ----
        2720             :
        2721           0 :     if ( stats_metadata%l_stats_samp ) then
        2722             :
        2723           0 :       if ( hm_metadata%iirr > 0 ) then
        2724             :         call stat_update_var( stats_metadata%irrm, hydromet(:,hm_metadata%iirr), & ! intent(in)
        2725           0 :                               stats_zt ) ! intent(inout)
        2726             :       end if
        2727             :
        2728           0 :       if ( hm_metadata%iirs > 0 ) then
        2729             :         call stat_update_var( stats_metadata%irsm, hydromet(:,hm_metadata%iirs), & ! intent(in)
        2730           0 :                               stats_zt ) ! intent(inout)
        2731             :       end if
        2732             :
        2733           0 :       if ( hm_metadata%iiri > 0 ) then
        2734             :         call stat_update_var( stats_metadata%irim, hydromet(:,hm_metadata%iiri), & ! intent(in)
        2735           0 :                               stats_zt ) ! intent(inout)
        2736             :       end if
        2737             :
        2738           0 :       if ( hm_metadata%iirg > 0 ) then
        2739             :         call stat_update_var( stats_metadata%irgm,  &  ! intent(in)
        2740             :                               hydromet(:,hm_metadata%iirg), & ! intent(in)
        2741           0 :                               stats_zt ) ! intent(inout)
        2742             :       end if
        2743             :
        2744           0 :       if ( hm_metadata%iiNi > 0 ) then
        2745             :         call stat_update_var( stats_metadata%iNim, hydromet(:,hm_metadata%iiNi), & ! intent(in)
        2746           0 :                               stats_zt ) ! intent(inout)
        2747             :       end if
        2748             :
        2749           0 :       if ( hm_metadata%iiNr > 0 ) then
        2750             :         call stat_update_var( stats_metadata%iNrm, hydromet(:,hm_metadata%iiNr), & ! intent(in)
        2751           0 :                               stats_zt ) ! intent(inout)
        2752             :       end if
        2753             :
        2754           0 :       if ( hm_metadata%iiNs > 0 ) then
        2755             :         call stat_update_var( stats_metadata%iNsm, hydromet(:,hm_metadata%iiNs), & ! intent(in)
        2756           0 :                               stats_zt ) ! intent(inout)
        2757             :       end if
        2758             :
        2759           0 :       if ( hm_metadata%iiNg > 0 ) then
        2760             :         call stat_update_var( stats_metadata%iNgm, hydromet(:,hm_metadata%iiNg), & ! intent(in)
        2761           0 :                               stats_zt ) ! intent(inout)
        2762             :       end if
        2763             :
        2764             :       ! Snow Water Path
        2765           0 :       if ( stats_metadata%iswp > 0 .and. hm_metadata%iirs > 0 ) then
        2766             :
        2767             :         ! Calculate snow water path
        2768             :         xtmp &
        2769             :         = vertical_integral &
        2770           0 :                ( (gr%nz - 2 + 1), rho_ds_zt(2:gr%nz), &
        2771           0 :                  hydromet(2:gr%nz,hm_metadata%iirs), gr%dzt(1,2:gr%nz) )
        2772             :
        2773             :         call stat_update_var_pt( stats_metadata%iswp, grid_level, xtmp, & ! intent(in)
        2774           0 :                                  stats_sfc ) ! intent(inout)
        2775             :
        2776             :       end if ! stats_metadata%iswp > 0 .and. iirs > 0
        2777             :
        2778             :       ! Ice Water Path
        2779           0 :       if ( stats_metadata%iiwp > 0 .and. hm_metadata%iiri > 0 ) then
        2780             :
        2781             :         xtmp &
        2782             :         = vertical_integral &
        2783           0 :                ( (gr%nz - 2 + 1), rho_ds_zt(2:gr%nz), &
        2784           0 :                  hydromet(2:gr%nz,hm_metadata%iiri), gr%dzt(1,2:gr%nz) )
        2785             :
        2786             :         call stat_update_var_pt( stats_metadata%iiwp, grid_level, xtmp, & ! intent(in)
        2787           0 :                                  stats_sfc ) ! intent(inout)
        2788             :
        2789             :       end if
        2790             :
        2791             :       ! Rain Water Path
        2792           0 :       if ( stats_metadata%irwp > 0 .and. hm_metadata%iirr > 0 ) then
        2793             :
        2794             :         xtmp &
        2795             :         = vertical_integral &
        2796           0 :                ( (gr%nz - 2 + 1), rho_ds_zt(2:gr%nz), &
        2797           0 :                  hydromet(2:gr%nz,hm_metadata%iirr), gr%dzt(1,2:gr%nz) )
        2798             :
        2799             :         call stat_update_var_pt( stats_metadata%irwp, grid_level, xtmp, & ! intent(in)
        2800           0 :                                  stats_sfc ) ! intent(inout)
        2801             :
        2802             :       end if ! stats_metadata%irwp > 0 .and. stats_metadata%irrm > 0
        2803             :     end if ! stats_metadata%l_stats_samp
        2804             :
        2805           0 :     return
        2806             :   end subroutine stats_accumulate_hydromet
        2807             : !------------------------------------------------------------------------------
        2808           0 :   subroutine stats_accumulate_lh_tend( gr, hydromet_dim, hm_metadata, &
        2809           0 :                                        lh_hydromet_mc, lh_Ncm_mc, &
        2810           0 :                                        lh_thlm_mc, lh_rvm_mc, lh_rcm_mc, &
        2811           0 :                                        lh_AKm, AKm, AKstd, AKstd_cld, &
        2812           0 :                                        lh_rcm_avg, AKm_rcm, AKm_rcc, &
        2813             :                                        stats_metadata, &
        2814             :                                        stats_lh_zt )
        2815             :
        2816             : ! Description:
        2817             : !   Compute stats for the tendency of latin hypercube sample points.
        2818             :
        2819             : ! References:
        2820             : !   None
        2821             : !------------------------------------------------------------------------------
        2822             :
        2823             :     use grid_class, only: &
        2824             :         grid ! Type
        2825             :
        2826             :     use stats_type_utilities, only: &
        2827             :         stat_update_var ! Procedure(s)
        2828             :
        2829             :     use stats_variables, only: &
        2830             :         stats_metadata_type
        2831             :
        2832             :     use clubb_precision, only: &
        2833             :         core_rknd ! Variable(s)
        2834             :
        2835             :     use corr_varnce_module, only: &
        2836             :       hm_metadata_type
        2837             :
        2838             :     use stats_type, only: &
        2839             :       stats ! Type
        2840             :
        2841             :     implicit none
        2842             :
        2843             :     !----------------------- Input Variables -----------------------
        2844             :     type (grid), target, intent(in) :: &
        2845             :       gr
        2846             :
        2847             :     integer, intent(in) :: &
        2848             :       hydromet_dim
        2849             :
        2850             :     type (hm_metadata_type), intent(in) :: &
        2851             :       hm_metadata
        2852             :
        2853             :     real( kind = core_rknd ), dimension(gr%nz,hydromet_dim), intent(in) :: &
        2854             :       lh_hydromet_mc ! Tendency of hydrometeors except for rvm/rcm  [units vary]
        2855             :
        2856             :     real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
        2857             :       lh_Ncm_mc,  & ! Tendency of cloud droplet concentration  [num/kg/s]
        2858             :       lh_thlm_mc, & ! Tendency of liquid potential temperature [kg/kg/s]
        2859             :       lh_rcm_mc,  & ! Tendency of cloud water                  [kg/kg/s]
        2860             :       lh_rvm_mc     ! Tendency of vapor                        [kg/kg/s]
        2861             :
        2862             :     real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
        2863             :       lh_AKm,     & ! Kessler ac estimate                 [kg/kg/s]
        2864             :       AKm,        & ! Exact Kessler ac                    [kg/kg/s]
        2865             :       AKstd,      & ! St dev of exact Kessler ac          [kg/kg/s]
        2866             :       AKstd_cld,  & ! Stdev of exact w/in cloud ac        [kg/kg/s]
        2867             :       lh_rcm_avg, & ! Monte Carlo rcm estimate            [kg/kg]
        2868             :       AKm_rcm,    & ! Kessler ac based on rcm             [kg/kg/s]
        2869             :       AKm_rcc       ! Kessler ac based on rcm/cloud_frac  [kg/kg/s]
        2870             :
        2871             :     type (stats_metadata_type), intent(in) :: &
        2872             :       stats_metadata
        2873             :
        2874             :     !----------------------- InOut Variables -----------------------
        2875             :     type (stats), target, intent(inout) :: &
        2876             :       stats_lh_zt
        2877             :
        2878             :     !----------------------- Begin Code -----------------------
        2879             :
        2880           0 :     if ( stats_metadata%l_stats_samp ) then
        2881             :
        2882             :       call stat_update_var( stats_metadata%ilh_thlm_mc, lh_thlm_mc, & ! intent(in)
        2883           0 :                             stats_lh_zt ) ! intent(inout)
        2884             :       call stat_update_var( stats_metadata%ilh_rcm_mc, lh_rcm_mc, & ! intent(in)
        2885           0 :                             stats_lh_zt ) ! intent(inout)
        2886             :       call stat_update_var( stats_metadata%ilh_rvm_mc, lh_rvm_mc, & ! intent(in)
        2887           0 :                             stats_lh_zt ) ! intent(inout)
        2888             :
        2889             :       call stat_update_var( stats_metadata%ilh_Ncm_mc, lh_Ncm_mc, & ! intent(in)
        2890           0 :                             stats_lh_zt ) ! intent(inout)
        2891             :
        2892           0 :       if ( hm_metadata%iirr > 0 ) then
        2893             :         call stat_update_var( stats_metadata%ilh_rrm_mc, lh_hydromet_mc(:,hm_metadata%iirr), & ! intent(in)
        2894           0 :                               stats_lh_zt ) ! intent(inout)
        2895             :       end if
        2896             :
        2897           0 :       if ( hm_metadata%iirs > 0 ) then
        2898             :         call stat_update_var( stats_metadata%ilh_rsm_mc,              & ! intent(in)
        2899             :                               lh_hydromet_mc(:,hm_metadata%iirs),  & ! intent(in)
        2900           0 :                               stats_lh_zt )                             ! intent(inout)
        2901             :       end if
        2902             :
        2903           0 :       if ( hm_metadata%iiri > 0 ) then
        2904             :         call stat_update_var( stats_metadata%ilh_rim_mc,              & ! intent(in)
        2905             :                               lh_hydromet_mc(:,hm_metadata%iiri),  & ! intent(in)
        2906           0 :                               stats_lh_zt )                             ! intent(inout)
        2907             :       end if
        2908             :
        2909           0 :       if ( hm_metadata%iirg > 0 ) then
        2910             :         call stat_update_var( stats_metadata%ilh_rgm_mc,              & ! intent(in)
        2911             :                               lh_hydromet_mc(:,hm_metadata%iirg),  & ! intent(in)
        2912           0 :                               stats_lh_zt )                             ! intent(inout)
        2913             :       end if
        2914             :
        2915           0 :       if ( hm_metadata%iiNi > 0 ) then
        2916             :         call stat_update_var( stats_metadata%ilh_Nim_mc,              & ! intent(in)
        2917             :                               lh_hydromet_mc(:,hm_metadata%iiNi),  & ! intent(in)
        2918           0 :                               stats_lh_zt )                             ! intent(inout)
        2919             :       end if
        2920             :
        2921           0 :       if ( hm_metadata%iiNr > 0 ) then
        2922             :         call stat_update_var( stats_metadata%ilh_Nrm_mc,              & ! intent(in)
        2923             :                               lh_hydromet_mc(:,hm_metadata%iiNr),  & ! intent(in)
        2924           0 :                               stats_lh_zt )                             ! intent(inout)
        2925             :       end if
        2926             :
        2927           0 :       if ( hm_metadata%iiNs > 0 ) then
        2928             :         call stat_update_var( stats_metadata%ilh_Nsm_mc,              & ! intent(in)
        2929             :                               lh_hydromet_mc(:,hm_metadata%iiNs),  & ! intent(in)
        2930           0 :                               stats_lh_zt )                             ! intent(inout)
        2931             :       end if
        2932             :
        2933           0 :       if ( hm_metadata%iiNg > 0 ) then
        2934             :         call stat_update_var( stats_metadata%ilh_Ngm_mc,              & ! intent(in)
        2935             :                               lh_hydromet_mc(:,hm_metadata%iiNg),  & ! intent(in)
        2936           0 :                               stats_lh_zt )                             ! intent(inout)
        2937             :       end if
        2938             :
        2939             :       call stat_update_var( stats_metadata%iAKm, AKm, & ! intent(in)
        2940           0 :                             stats_lh_zt )               ! intent(inout)
        2941             :
        2942             :       call stat_update_var( stats_metadata%ilh_AKm, lh_AKm, & ! intent(in)
        2943           0 :                             stats_lh_zt)                      ! intent(inout)
        2944             :
        2945             :       call stat_update_var( stats_metadata%ilh_rcm_avg, lh_rcm_avg, & ! intent(in)
        2946           0 :                             stats_lh_zt )                             ! intent(inout)
        2947             :
        2948             :       call stat_update_var( stats_metadata%iAKstd, AKstd, & ! intent(in)
        2949           0 :                             stats_lh_zt )                   ! intent(inout)
        2950             :
        2951             :       call stat_update_var( stats_metadata%iAKstd_cld, AKstd_cld, & ! intent(in)
        2952           0 :                             stats_lh_zt )                       ! intent(inout)
        2953             :
        2954             :       call stat_update_var( stats_metadata%iAKm_rcm, AKm_rcm, & ! intent(in)
        2955           0 :                             stats_lh_zt)                        ! intent(inout)
        2956             :
        2957             :       call stat_update_var( stats_metadata%iAKm_rcc, AKm_rcc, & ! intent(in)
        2958           0 :                             stats_lh_zt )                       ! intent(inout)
        2959             :
        2960             :     end if ! stats_metadata%l_stats_samp
        2961             :
        2962           0 :     return
        2963             :   end subroutine stats_accumulate_lh_tend
        2964             :
        2965             :   !-----------------------------------------------------------------------
        2966           0 :   subroutine stats_finalize( stats_metadata, &
        2967             :                              stats_zt, stats_zm, stats_sfc, &
        2968             :                              stats_lh_zt, stats_lh_sfc, &
        2969             :                              stats_rad_zt, stats_rad_zm )
        2970             :
        2971             :     !     Description:
        2972             :     !     Close NetCDF files and deallocate scratch space and
        2973             :     !     stats file structures.
        2974             :     !-----------------------------------------------------------------------
        2975             :
        2976             :     use stats_variables, only: &
        2977             :         stats_metadata_type
        2978             :
        2979             : #ifdef NETCDF
        2980             :     use output_netcdf, only:  &
        2981             :         close_netcdf ! Procedure
        2982             : #endif
        2983             :
        2984             :     use stats_type, only: stats ! Type
        2985             :
        2986             :     implicit none
        2987             :
        2988             :     type (stats), target, intent(inout) :: &
        2989             :       stats_zt, &
        2990             :       stats_zm, &
        2991             :       stats_sfc, &
        2992             :       stats_lh_zt, &
        2993             :       stats_lh_sfc, &
        2994             :       stats_rad_zt, &
        2995             :       stats_rad_zm
        2996             :
        2997             :     type (stats_metadata_type), intent(inout) :: &
        2998             :       stats_metadata
        2999             :
        3000           0 :     if ( stats_metadata%l_stats .and. stats_metadata%l_netcdf ) then
        3001             : #ifdef NETCDF
        3002             :       call close_netcdf( stats_zt%file ) ! intent(inout)
        3003             :       call close_netcdf( stats_lh_zt%file ) ! intent(inout)
        3004             :       call close_netcdf( stats_lh_sfc%file ) ! intent(inout)
        3005             :       call close_netcdf( stats_zm%file ) ! intent(inout)
        3006             :       call close_netcdf( stats_rad_zt%file ) ! intent(inout)
        3007             :       call close_netcdf( stats_rad_zm%file ) ! intent(inout)
        3008             :       call close_netcdf( stats_sfc%file ) ! intent(inout)
        3009             : #else
        3010           0 :       error stop "This program was not compiled with netCDF support"
        3011             : #endif
        3012             :     end if
        3013             :
        3014           0 :     if ( stats_metadata%l_stats ) then
        3015             :       ! De-allocate all stats_zt variables
        3016           0 :       if (allocated(stats_zt%z)) then
        3017           0 :         deallocate( stats_zt%z )
        3018             :
        3019           0 :         deallocate( stats_zt%accum_field_values )
        3020           0 :         deallocate( stats_zt%accum_num_samples )
        3021           0 :         deallocate( stats_zt%l_in_update )
        3022             :
        3023           0 :         deallocate( stats_zt%file%grid_avg_var )
        3024           0 :         deallocate( stats_zt%file%z )
        3025             :
        3026             :         ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3027           0 :         if ( allocated( stats_zt%file%lat_vals ) ) then
        3028           0 :           deallocate( stats_zt%file%lat_vals )
        3029             :         end if
        3030           0 :         if ( allocated( stats_zt%file%lon_vals ) ) then
        3031           0 :           deallocate( stats_zt%file%lon_vals )
        3032             :         end if
        3033             :
        3034             :       end if
        3035             :
        3036           0 :       if ( stats_metadata%l_silhs_out .and. allocated(stats_lh_zt%z) ) then
        3037             :         ! De-allocate all stats_lh_zt variables
        3038           0 :         deallocate( stats_lh_zt%z )
        3039             :
        3040           0 :         deallocate( stats_lh_zt%accum_field_values )
        3041           0 :         deallocate( stats_lh_zt%accum_num_samples )
        3042           0 :         deallocate( stats_lh_zt%l_in_update )
        3043             :
        3044           0 :         deallocate( stats_lh_zt%file%grid_avg_var )
        3045           0 :         deallocate( stats_lh_zt%file%z )
        3046             :
        3047             :         ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3048           0 :         if ( allocated(stats_lh_zt%file%lat_vals ) ) then
        3049           0 :           deallocate( stats_lh_zt%file%lat_vals )
        3050             :         end if
        3051           0 :         if ( allocated(stats_lh_zt%file%lon_vals ) ) then
        3052           0 :           deallocate( stats_lh_zt%file%lon_vals )
        3053             :         end if
        3054             :
        3055             :         ! De-allocate all stats_lh_sfc variables
        3056           0 :         deallocate( stats_lh_sfc%z )
        3057             :
        3058           0 :         deallocate( stats_lh_sfc%accum_field_values )
        3059           0 :         deallocate( stats_lh_sfc%accum_num_samples )
        3060           0 :         deallocate( stats_lh_sfc%l_in_update )
        3061             :
        3062           0 :         deallocate( stats_lh_sfc%file%grid_avg_var )
        3063           0 :         deallocate( stats_lh_sfc%file%z )
        3064             :
        3065             :         ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3066           0 :         if ( allocated( stats_lh_sfc%file%lat_vals ) ) then
        3067           0 :           deallocate( stats_lh_sfc%file%lat_vals )
        3068             :         end if
        3069           0 :         if ( allocated( stats_lh_sfc%file%lon_vals ) ) then
        3070           0 :           deallocate( stats_lh_sfc%file%lon_vals )
        3071             :         end if
        3072             :       end if ! l_silhs_out
        3073             :
        3074             :       ! De-allocate all stats_zm variables
        3075           0 :       if (allocated(stats_zm%z)) then
        3076           0 :         deallocate( stats_zm%z )
        3077             :
        3078           0 :         deallocate( stats_zm%accum_field_values )
        3079           0 :         deallocate( stats_zm%accum_num_samples )
        3080           0 :         deallocate( stats_zm%l_in_update )
        3081             :
        3082           0 :         deallocate( stats_zm%file%grid_avg_var )
        3083           0 :         deallocate( stats_zm%file%z )
        3084             :
        3085             :         ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3086           0 :         if ( allocated( stats_zm%file%lat_vals ) ) then
        3087           0 :           deallocate( stats_zm%file%lat_vals )
        3088             :         end if
        3089           0 :         if ( allocated( stats_zm%file%lon_vals ) ) then
        3090           0 :           deallocate( stats_zm%file%lon_vals )
        3091             :         end if
        3092             :
        3093             :       end if
        3094             :
        3095           0 :       if ( stats_metadata%l_output_rad_files ) then
        3096             :         ! De-allocate all stats_rad_zt variables
        3097           0 :         if (allocated(stats_rad_zt%z)) then
        3098           0 :           deallocate( stats_rad_zt%z )
        3099             :
        3100           0 :           deallocate( stats_rad_zt%accum_field_values )
        3101           0 :           deallocate( stats_rad_zt%accum_num_samples )
        3102           0 :           deallocate( stats_rad_zt%l_in_update )
        3103             :
        3104           0 :           deallocate( stats_rad_zt%file%grid_avg_var )
        3105           0 :           deallocate( stats_rad_zt%file%z )
        3106             :
        3107             :           ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3108           0 :           if ( allocated( stats_rad_zt%file%lat_vals ) ) then
        3109           0 :             deallocate( stats_rad_zt%file%lat_vals )
        3110             :           end if
        3111           0 :           if ( allocated( stats_rad_zt%file%lon_vals ) ) then
        3112           0 :             deallocate( stats_rad_zt%file%lon_vals )
        3113             :           end if
        3114             :
        3115             :           ! De-allocate all stats_rad_zm variables
        3116           0 :           deallocate( stats_rad_zm%z )
        3117             :
        3118           0 :           deallocate( stats_rad_zm%accum_field_values )
        3119           0 :           deallocate( stats_rad_zm%accum_num_samples )
        3120           0 :           deallocate( stats_rad_zm%l_in_update )
        3121             :
        3122           0 :           deallocate( stats_rad_zm%file%grid_avg_var )
        3123           0 :           deallocate( stats_rad_zm%file%z )
        3124             :
        3125             :           ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3126           0 :           if ( allocated( stats_rad_zm%file%lat_vals ) ) then
        3127           0 :             deallocate( stats_rad_zm%file%lat_vals )
        3128             :           end if
        3129           0 :           if ( allocated( stats_rad_zm%file%lon_vals ) ) then
        3130           0 :             deallocate( stats_rad_zm%file%lon_vals )
        3131             :           end if
        3132             :
        3133             :         end if
        3134             :
        3135             :       end if ! l_output_rad_files
        3136             :
        3137             :       ! De-allocate all stats_sfc variables
        3138           0 :       if (allocated(stats_sfc%z)) then
        3139           0 :         deallocate( stats_sfc%z )
        3140             :
        3141           0 :         deallocate( stats_sfc%accum_field_values )
        3142           0 :         deallocate( stats_sfc%accum_num_samples )
        3143           0 :         deallocate( stats_sfc%l_in_update )
        3144             :
        3145           0 :         deallocate( stats_sfc%file%grid_avg_var )
        3146           0 :         deallocate( stats_sfc%file%z )
        3147             :       end if
        3148             :
        3149             :       ! Check if pointer is allocated to prevent crash in netcdf (ticket 765)
        3150           0 :       if ( allocated( stats_sfc%file%lat_vals ) ) then
        3151           0 :         deallocate( stats_sfc%file%lat_vals )
        3152             :       end if
        3153           0 :       if ( allocated( stats_sfc%file%lon_vals ) ) then
        3154           0 :         deallocate( stats_sfc%file%lon_vals )
        3155             :       end if
        3156             :
        3157             :       ! De-allocate scalar indices
        3158           0 :       if (allocated(stats_metadata%isclrm)) then
        3159           0 :         deallocate( stats_metadata%isclrm )
        3160           0 :         deallocate( stats_metadata%isclrm_f )
        3161           0 :         deallocate( stats_metadata%iedsclrm )
        3162           0 :         deallocate( stats_metadata%iedsclrm_f )
        3163           0 :         deallocate( stats_metadata%isclrprtp )
        3164           0 :         deallocate( stats_metadata%isclrp2 )
        3165           0 :         deallocate( stats_metadata%isclrpthvp )
        3166           0 :         deallocate( stats_metadata%isclrpthlp )
        3167           0 :         deallocate( stats_metadata%isclrprcp )
        3168           0 :         deallocate( stats_metadata%iwpsclrp )
        3169           0 :         deallocate( stats_metadata%iwp2sclrp )
        3170           0 :         deallocate( stats_metadata%iwpsclrp2 )
        3171           0 :         deallocate( stats_metadata%iwpsclrprtp )
        3172           0 :         deallocate( stats_metadata%iwpsclrpthlp )
        3173           0 :         deallocate( stats_metadata%iwpedsclrp )
        3174             :       end if
        3175             :
        3176             :       ! De-allocate hyderometeor statistical variables
        3177           0 :       if (allocated(stats_metadata%ihm_1)) then
        3178           0 :         deallocate( stats_metadata%ihm_1 )
        3179           0 :         deallocate( stats_metadata%ihm_2 )
        3180           0 :         deallocate( stats_metadata%imu_hm_1 )
        3181           0 :         deallocate( stats_metadata%imu_hm_2 )
        3182           0 :         deallocate( stats_metadata%imu_hm_1_n )
        3183           0 :         deallocate( stats_metadata%imu_hm_2_n )
        3184           0 :         deallocate( stats_metadata%isigma_hm_1 )
        3185           0 :         deallocate( stats_metadata%isigma_hm_2 )
        3186           0 :         deallocate( stats_metadata%isigma_hm_1_n )
        3187           0 :         deallocate( stats_metadata%isigma_hm_2_n )
        3188           0 :         deallocate( stats_metadata%icorr_w_hm_1 )
        3189           0 :         deallocate( stats_metadata%icorr_w_hm_2 )
        3190           0 :         deallocate( stats_metadata%icorr_chi_hm_1 )
        3191           0 :         deallocate( stats_metadata%icorr_chi_hm_2 )
        3192           0 :         deallocate( stats_metadata%icorr_eta_hm_1 )
        3193           0 :         deallocate( stats_metadata%icorr_eta_hm_2 )
        3194           0 :         deallocate( stats_metadata%icorr_Ncn_hm_1 )
        3195           0 :         deallocate( stats_metadata%icorr_Ncn_hm_2 )
        3196           0 :         deallocate( stats_metadata%icorr_hmx_hmy_1 )
        3197           0 :         deallocate( stats_metadata%icorr_hmx_hmy_2 )
        3198           0 :         deallocate( stats_metadata%icorr_w_hm_1_n )
        3199           0 :         deallocate( stats_metadata%icorr_w_hm_2_n )
        3200           0 :         deallocate( stats_metadata%icorr_chi_hm_1_n )
        3201           0 :         deallocate( stats_metadata%icorr_chi_hm_2_n )
        3202           0 :         deallocate( stats_metadata%icorr_eta_hm_1_n )
        3203           0 :         deallocate( stats_metadata%icorr_eta_hm_2_n )
        3204           0 :         deallocate( stats_metadata%icorr_Ncn_hm_1_n )
        3205           0 :         deallocate( stats_metadata%icorr_Ncn_hm_2_n )
        3206           0 :         deallocate( stats_metadata%icorr_hmx_hmy_1_n )
        3207           0 :         deallocate( stats_metadata%icorr_hmx_hmy_2_n )
        3208           0 :         deallocate( stats_metadata%ihmp2_zt )
        3209           0 :         deallocate( stats_metadata%iwp2hmp )
        3210           0 :         deallocate( stats_metadata%ihydrometp2 )
        3211           0 :         deallocate( stats_metadata%iwphydrometp )
        3212           0 :         deallocate( stats_metadata%irtphmp )
        3213           0 :         deallocate( stats_metadata%ithlphmp )
        3214           0 :         deallocate( stats_metadata%ihmxphmyp )
        3215           0 :         deallocate( stats_metadata%iK_hm )
        3216             :       end if
        3217             :
        3218           0 :       if ( allocated( stats_metadata%isilhs_variance_category ) ) then
        3219           0 :         deallocate( stats_metadata%isilhs_variance_category )
        3220           0 :         deallocate( stats_metadata%ilh_samp_frac_category )
        3221             :       end if
        3222             :
        3223             :     end if ! l_stats
        3224             :
        3225           0 :     return
        3226             :   end subroutine stats_finalize
        3227             :
        3228             : !===============================================================================
        3229             :
        3230             : !-----------------------------------------------------------------------
        3231           0 : subroutine stats_check_num_samples( stats_grid, stats_metadata, &
        3232             :                                     l_error )
        3233             :
        3234             : ! Description:
        3235             : !   Ensures that each variable in a stats grid is sampled the correct
        3236             : !   number of times.
        3237             : ! References:
        3238             : !   None
        3239             : !-----------------------------------------------------------------------
        3240             :
        3241             :     use constants_clubb, only: &
        3242             :         fstderr ! Constant
        3243             :
        3244             :     use stats_type, only: &
        3245             :         stats ! Type
        3246             :
        3247             :     use stats_variables, only: &
        3248             :         stats_metadata_type
        3249             :
        3250             :     use error_code, only: &
        3251             :         clubb_at_least_debug_level   ! Procedure
        3252             :
        3253             :     implicit none
        3254             :
        3255             :   ! Input Variables
        3256             :     type (stats), intent(in) :: &
        3257             :       stats_grid               ! Grid type              [grid]
        3258             :
        3259             :     type (stats_metadata_type), intent(in) :: &
        3260             :       stats_metadata
        3261             :
        3262             :   ! Input/Output Variables
        3263             :     logical, intent(inout) :: &
        3264             :       l_error                  ! Indicates an error     [boolean]
        3265             :
        3266             :   ! Local Variables
        3267             :     integer :: ivar, kvar      ! Loop variable          [index]
        3268             :
        3269             :     logical :: l_proper_sample
        3270             :
        3271             : !-----------------------------------------------------------------------
        3272             :
        3273             :   !----- Begin Code -----
        3274             :
        3275             :   ! Look for errors by checking the number of sampling points
        3276             :   ! for each variable in the statistics grid at each vertical level.
        3277           0 :   do ivar = 1, stats_grid%num_output_fields
        3278           0 :     do kvar = 1, stats_grid%kk
        3279             :
        3280           0 :       l_proper_sample = ( stats_grid%accum_num_samples(1,1,kvar,ivar) == 0 .or. &
        3281             :                           stats_grid%accum_num_samples(1,1,kvar,ivar) &
        3282             :                           == floor(  stats_metadata%stats_tout  &
        3283           0 :                                    / stats_metadata%stats_tsamp ) )
        3284             :
        3285           0 :       if ( .not. l_proper_sample ) then
        3286             :
        3287           0 :         l_error = .true.  ! This will stop the run
        3288             :
        3289           0 :         if ( clubb_at_least_debug_level( 1 ) ) then
        3290           0 :           write(fstderr,*) 'Possible sampling error for variable ',  &
        3291           0 :                            trim(stats_grid%file%grid_avg_var(ivar)%name), ' in stats_grid ',  &
        3292           0 :                            'at k = ', kvar,  &
        3293           0 :                            '; stats_grid%accum_num_samples(',kvar,',',ivar,') = ', &
        3294           0 :                             stats_grid%accum_num_samples(1,1,kvar,ivar)
        3295             :         end if ! clubb_at_lest_debug_level 1
        3296             :
        3297             :
        3298             :       end if ! .not. l_proper_sample
        3299             :
        3300             :     end do ! kvar = 1 .. stats_grid%kk
        3301             :   end do ! ivar = 1 .. stats_grid%num_output_fields
        3302             :
        3303           0 :   return
        3304             : end subroutine stats_check_num_samples
        3305             : !-----------------------------------------------------------------------
        3306             :
        3307             : end module stats_clubb_utilities

Generated by: LCOV version 1.14