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