LCOV - code coverage report
Current view: | top level - physics/cam - cospsimulator_intr.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 5 | 8 | 62.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 4 | 50.0 %

              Line data    Source code
           1             : module cospsimulator_intr
           2             :   ! ######################################################################################
           3             :   ! Purpose: CAM interface to
           4             :   !         Name:         CFMIP Observational Simulator Package Version 2 (COSP2)
           5             :   !         What:         Simulate ISCCP/CloudSat/CALIPSO/MISR/MODIS cloud products from
           6             :   !                       GCM inputs
           7             :   !         Version:      v2.1.4 (August 2019)
           8             :   !         Authors:      Dustin Swales (dustin.swales@noaa.gov)
           9             :   !
          10             :   ! Modifications:
          11             :   !
          12             :   ! ######################################################################################
          13             :   use shr_kind_mod,         only: r8 => shr_kind_r8
          14             :   use spmd_utils,           only: masterproc
          15             :   use ppgrid,               only: pcols, pver, pverp, begchunk, endchunk
          16             :   use ref_pres,             only: ktop => trop_cloud_top_lev
          17             :   use perf_mod,             only: t_startf, t_stopf
          18             :   use cam_abortutils,       only: endrun, handle_allocate_error
          19             :   use phys_control,         only: cam_physpkg_is
          20             :   use cam_logfile,          only: iulog
          21             : #ifdef USE_COSP
          22             :   use quickbeam,            only: radar_cfg
          23             :   use mod_quickbeam_optics, only: size_distribution
          24             :   use mod_cosp,             only: cosp_outputs, cosp_optical_inputs, cosp_column_inputs
          25             :   use mod_cosp_config,      only: pres_binCenters, pres_binEdges, tau_binCenters,      &
          26             :        tau_binEdges, cloudsat_binCenters, cloudsat_binEdges, calipso_binCenters,       &
          27             :        calipso_binEdges, misr_histHgtCenters, misr_histHgtEdges,  PARASOL_SZA,         &
          28             :        R_UNDEF, PARASOL_NREFL, LIDAR_NCAT,SR_BINS, N_HYDRO, RTTOV_MAX_CHANNELS,        &
          29             :        numMISRHgtBins, CLOUDSAT_DBZE_BINS, LIDAR_NTEMP, calipso_histBsct,              &
          30             :        numMODISTauBins, numMODISPresBins, numMODISReffIceBins, numMODISReffLiqBins,    &
          31             :        numISCCPTauBins, numISCCPPresBins, numMISRTauBins, reffICE_binEdges,            &
          32             :        reffICE_binCenters, reffLIQ_binEdges, reffLIQ_binCenters, LIDAR_NTYPE,          &
          33             :        nCloudsatPrecipClass, &
          34             :        nsza_cosp         => PARASOL_NREFL,       &
          35             :        nprs_cosp         => npres,               &
          36             :        ntau_cosp         => ntau,                &
          37             :        ntau_cosp_modis   => ntau,                &
          38             :        nsr_cosp          => SR_BINS,             &
          39             :        nhtmisr_cosp      => numMISRHgtBins,      &
          40             :        nhydro            => N_HYDRO, &
          41             :        cloudsat_preclvl
          42             :     use mod_cosp_stats,       only: cosp_change_vertical_grid
          43             : #endif
          44             :   implicit none
          45             :   private
          46             :   save
          47             :
          48             :   ! Public functions/subroutines
          49             :   public :: &
          50             :        cospsimulator_intr_readnl,  &
          51             :        cospsimulator_intr_register,&
          52             :        cospsimulator_intr_init,    &
          53             :        cospsimulator_intr_run
          54             :
          55             :   ! ######################################################################################
          56             :   ! Public declarations
          57             :   ! ######################################################################################
          58             :   ! Whether to do COSP calcs and I/O, default is false. If docosp is specified in
          59             :   ! the atm_in namelist, this value is overwritten and cosp is run
          60             :   logical, public, protected :: docosp = .false.
          61             :
          62             :   ! Frequency at which cosp is called, every cosp_nradsteps radiation timestep
          63             :   integer, public, protected :: cosp_nradsteps = 1
          64             :
          65             : #ifdef USE_COSP
          66             :
          67             :   ! ######################################################################################
          68             :   ! Local declarations
          69             :   ! ######################################################################################
          70             :   integer ::  &
          71             :        nlay,        &     ! Number of CAM layers used by COSP.
          72             :        nlayp,       &     ! Number of CAM layer interfaces used by COSP.
          73             :        nscol_cosp,  &     ! Number of subcolumns, allow namelist input to set.
          74             :        nht_cosp           ! Number of height for COSP radar and calipso simulator outputs.
          75             :                           !  *set to 40 if csat_vgrid=.true., else set to Nlr*
          76             :
          77             :
          78             :   ! ######################################################################################
          79             :   ! Bin-boundaries for mixed dimensions. Calculated in cospsetupvales OR in cosp_config.F90
          80             :   ! ######################################################################################
          81             :   real(r8), target :: prsmid_cosp(nprs_cosp)            ! pressure midpoints of COSP ISCCP output
          82             :   real(r8), target :: prslim_cosp(2,nprs_cosp)
          83             :   real(r8), target :: taumid_cosp(ntau_cosp)            ! optical depth midpoints of COSP ISCCP output
          84             :   real(r8), target :: taulim_cosp(2,ntau_cosp)
          85             :   real(r8), target :: srmid_cosp(nsr_cosp)              ! sr midpoints of COSP lidar output
          86             :   real(r8), target :: srlim_cosp(2,nsr_cosp)
          87             :   real(r8), target :: sza_cosp(nsza_cosp)
          88             :   real(r8), target :: dbzemid_cosp(CLOUDSAT_DBZE_BINS)          ! dbze midpoints of COSP radar output
          89             :   real(r8), target :: dbzelim_cosp(2,CLOUDSAT_DBZE_BINS)
          90             :   real(r8), target :: htmisrmid_cosp(nhtmisr_cosp)      ! htmisr midpoints of COSP misr simulator output
          91             :   real(r8), target :: htmisrlim_cosp(2,nhtmisr_cosp)
          92             :   real(r8), target :: taumid_cosp_modis(ntau_cosp_modis)! optical depth midpoints of COSP MODIS output
          93             :   real(r8), target :: taulim_cosp_modis(2,ntau_cosp_modis)
          94             :   real(r8), target :: reffICE_binEdges_cosp(2,numMODISReffIceBins)
          95             :   real(r8), target :: reffLIQ_binEdges_cosp(2,numMODISReffLiqBins)
          96             :   real(r8), target :: reffICE_binCenters_cosp(numMODISReffIceBins)
          97             :   real(r8), target :: reffLIQ_binCenters_cosp(numMODISReffLiqBins)
          98             :
          99             :   integer  :: prstau_cosp(nprs_cosp*ntau_cosp)             ! ISCCP mixed output dimension index
         100             :   integer  :: prstau_cosp_modis(nprs_cosp*ntau_cosp_modis) ! MODIS mixed output dimension index
         101             :   integer  :: htmisrtau_cosp(nhtmisr_cosp*ntau_cosp)       ! MISR mixed output dimension index
         102             :   real(r8) :: prstau_prsmid_cosp(nprs_cosp*ntau_cosp)
         103             :   real(r8) :: prstau_taumid_cosp(nprs_cosp*ntau_cosp)
         104             :   real(r8) :: prstau_prsmid_cosp_modis(nprs_cosp*ntau_cosp_modis)
         105             :   real(r8) :: prstau_taumid_cosp_modis(nprs_cosp*ntau_cosp_modis)
         106             :   real(r8) :: htmisrtau_htmisrmid_cosp(nhtmisr_cosp*ntau_cosp)
         107             :   real(r8) :: htmisrtau_taumid_cosp(nhtmisr_cosp*ntau_cosp)
         108             :   real(r8),allocatable         :: htmlmid_cosp(:)          ! Model level height midpoints for output (nlay)
         109             :   real(r8),allocatable, public :: htdbze_dbzemid_cosp(:)   ! (nht_cosp*CLOUDSAT_DBZE_BINS)
         110             :   real(r8),allocatable, target :: htlim_cosp(:,:)          ! height limits for COSP outputs (nht_cosp+1)
         111             :   real(r8),allocatable, target :: htmid_cosp(:)            ! height midpoints of COSP radar/lidar output (nht_cosp)
         112             :   real(r8),allocatable         :: htlim_cosp_1d(:)         ! height limits for COSP outputs (nht_cosp+1)
         113             :   real(r8),allocatable         :: htdbze_htmid_cosp(:)     ! (nht_cosp*CLOUDSAT_DBZE_BINS)
         114             :   real(r8),allocatable         :: htsr_htmid_cosp(:)       ! (nht_cosp*nsr_cosp)
         115             :   real(r8),allocatable         :: htsr_srmid_cosp(:)       ! (nht_cosp*nsr_cosp)
         116             :   real(r8),allocatable         :: htmlscol_htmlmid_cosp(:) ! (nlay*nscol_cosp)
         117             :   real(r8),allocatable         :: htmlscol_scol_cosp(:)    ! (nlay*nscol_cosp)
         118             :   integer, allocatable, target :: scol_cosp(:)             ! sub-column number (nscol_cosp)
         119             :   integer, allocatable         :: htdbze_cosp(:)           ! radar CFAD mixed output dimension index (nht_cosp*CLOUDSAT_DBZE_BINS)
         120             :   integer, allocatable         :: htsr_cosp(:)             ! lidar CFAD mixed output dimension index (nht_cosp*nsr_cosp)
         121             :   integer, allocatable         :: htmlscol_cosp(:)         ! html-subcolumn mixed output dimension index (nlay*nscol_cosp)
         122             :
         123             :   ! ######################################################################################
         124             :   ! Default CAM namelist settings
         125             :   ! ######################################################################################
         126             :   logical :: cosp_amwg             = .false.
         127             :   logical :: cosp_lite             = .false.
         128             :   logical :: cosp_passive          = .false.
         129             :   logical :: cosp_active           = .false.
         130             :   logical :: cosp_isccp            = .false.
         131             :   logical :: cosp_lradar_sim       = .false.
         132             :   logical :: cosp_llidar_sim       = .false.
         133             :   logical :: cosp_lisccp_sim       = .false.
         134             :   logical :: cosp_lmisr_sim        = .false.
         135             :   logical :: cosp_lmodis_sim       = .false.
         136             :   logical :: cosp_histfile_aux     = .false.
         137             :   logical :: cosp_lfrac_out        = .false.
         138             :   logical :: cosp_runall           = .false.
         139             :   integer :: cosp_ncolumns         = 50
         140             :   integer :: cosp_histfile_num     =  1
         141             :   integer :: cosp_histfile_aux_num = -1
         142             :
         143             :   ! COSP
         144             :   logical :: lradar_sim       = .false.
         145             :   logical :: llidar_sim       = .false.
         146             :   logical :: lparasol_sim     = .false.
         147             :   logical :: lgrLidar532      = .false.
         148             :   logical :: latlid           = .false.
         149             :   logical :: lisccp_sim       = .false.
         150             :   logical :: lmisr_sim        = .false.
         151             :   logical :: lmodis_sim       = .false.
         152             :   logical :: lrttov_sim       = .false.
         153             :   logical :: lfrac_out        = .false.
         154             :
         155             :   ! ######################################################################################
         156             :   ! COSP parameters
         157             :   ! ######################################################################################
         158             :   integer, parameter :: Npoints_it = 10000       ! Max # gridpoints to be processed in one iteration (10,000)
         159             :   integer :: ncolumns = 50                       ! Number of subcolumns in SCOPS (50)
         160             :   integer :: nlr = 40                            ! Number of levels in statistical outputs
         161             :                                                  ! (only used if USE_VGRID=.true.)  (40)
         162             :   logical :: use_vgrid = .true.                  ! Use fixed vertical grid for outputs?
         163             :                                                  ! (if .true. then define # of levels with nlr)  (.true.)
         164             :   logical :: csat_vgrid = .true.                 ! CloudSat vertical grid?
         165             :
         166             :   ! Variables for COSP input related to radar simulator
         167             :   real(r8) :: radar_freq = 94.0_r8               ! CloudSat radar frequency (GHz) (94.0)
         168             :   integer :: surface_radar = 0                   ! surface=1, spaceborne=0 (0)
         169             :   integer :: use_gas_abs = 1                     ! include gaseous absorption? yes=1,no=0 (1)
         170             :   integer :: do_ray = 0                          ! calculate/output Rayleigh refl=1, not=0 (0)
         171             :   real(r8) :: k2 = -1                            ! |K|^2, -1=use frequency dependent default (-1)
         172             :
         173             :   ! Variables for COSP input related to lidar simulator
         174             :   integer, parameter :: Nprmts_max_hydro = 12    ! Max # params for hydrometeor size distributions (12)
         175             :   integer, parameter :: Naero = 1                ! Number of aerosol species (Not used) (1)
         176             :   integer, parameter :: Nprmts_max_aero = 1      ! Max # params for aerosol size distributions (not used) (1)
         177             :   integer :: lidar_ice_type = 0                  ! Ice particle shape in lidar calculations
         178             :                                                  ! (0=ice-spheres ; 1=ice-non-spherical) (0)
         179             :   integer, parameter :: overlap = 3              ! overlap type: 1=max, 2=rand, 3=max/rand (3)
         180             :
         181             :   ! Variables for COSP input related to ISCCP simulator
         182             :   integer :: isccp_topheight = 1                 ! 1 = adjust top height using both a computed infrared
         183             :                                                  ! brightness temperature and the visible
         184             :                                                  ! optical depth to adjust cloud top pressure.
         185             :                                                  ! Note that this calculation is most appropriate to compare
         186             :                                                  ! to ISCCP data during sunlit hours.
         187             :                                                  ! 2 = do not adjust top height, that is cloud top pressure
         188             :                                                  ! is the actual cloud top pressure in the model
         189             :                                                  ! 3 = adjust top height using only the computed infrared
         190             :                                                  ! brightness temperature. Note that this calculation is most
         191             :                                                  ! appropriate to compare to ISCCP IR only algortihm (i.e.
         192             :                                                  ! you can compare to nighttime ISCCP data with this option) (1)
         193             :   integer :: isccp_topheight_direction = 2       ! direction for finding atmosphere pressure level with
         194             :                                                  ! interpolated temperature equal to the radiance
         195             :                                                  ! determined cloud-top temperature
         196             :                                                  ! 1 = find the *lowest* altitude (highest pressure) level
         197             :                                                  ! with interpolated temperature
         198             :                                                  ! equal to the radiance determined cloud-top temperature
         199             :                                                  ! 2 = find the *highest* altitude (lowest pressure) level
         200             :                                                  ! with interpolated temperature
         201             :                                                  ! equal to the radiance determined cloud-top temperature
         202             :                                                  ! ONLY APPLICABLE IF top_height EQUALS 1 or 3
         203             :                                                  ! 1 = default setting in COSP v1.1, matches all versions of
         204             :                                                  ! ISCCP simulator with versions numbers 3.5.1 and lower
         205             :                                                  ! 2 = default setting in COSP v1.3. default since V4.0 of ISCCP simulator
         206             :
         207             :   ! ######################################################################################
         208             :   ! Other variables
         209             :   ! ######################################################################################
         210             :   logical,allocatable :: first_run_cosp(:)      !.true. if run_cosp has been populated (allocatable->begchunk:endchunk)
         211             :   logical,allocatable :: run_cosp(:,:)          !.true. if cosp should be run by column and
         212             :                                                 !       chunk (allocatable->1:pcols,begchunk:endchunk)
         213             :   ! pbuf indices
         214             :   integer :: cld_idx, concld_idx, lsreffrain_idx, lsreffsnow_idx, cvreffliq_idx
         215             :   integer :: cvreffice_idx
         216             :   integer :: gb_totcldliqmr_idx, gb_totcldicemr_idx
         217             :   integer :: dpflxprc_idx
         218             :   integer :: dpflxsnw_idx, shflxprc_idx, shflxsnw_idx, lsflxprc_idx, lsflxsnw_idx
         219             :   integer :: rei_idx, rel_idx
         220             :
         221             :   ! ######################################################################################
         222             :   ! Declarations specific to COSP2
         223             :   ! ######################################################################################
         224             :   type(radar_cfg)              :: rcfg_cloudsat ! Radar configuration (Cloudsat)
         225             :   type(radar_cfg), allocatable :: rcfg_cs(:)    ! chunked version of rcfg_cloudsat
         226             :   type(size_distribution)              :: sd       ! Size distribution used by radar simulator
         227             :   type(size_distribution), allocatable :: sd_cs(:) ! chunked version of sd
         228             :   character(len=64)         :: cloudsat_micro_scheme = 'MMF_v3.5_single_moment'
         229             :
         230             :   integer,parameter :: &
         231             :        I_LSCLIQ = 1, & ! Large-scale (stratiform) liquid
         232             :        I_LSCICE = 2, & ! Large-scale (stratiform) ice
         233             :        I_LSRAIN = 3, & ! Large-scale (stratiform) rain
         234             :        I_LSSNOW = 4, & ! Large-scale (stratiform) snow
         235             :        I_CVCLIQ = 5, & ! Convective liquid
         236             :        I_CVCICE = 6, & ! Convective ice
         237             :        I_CVRAIN = 7, & ! Convective rain
         238             :        I_CVSNOW = 8, & ! Convective snow
         239             :        I_LSGRPL = 9    ! Large-scale (stratiform) groupel
         240             :
         241             :   ! Stratiform and convective clouds in frac_out (scops output).
         242             :   integer, parameter :: &
         243             :        I_LSC = 1, & ! Large-scale clouds
         244             :        I_CVC = 2    ! Convective clouds
         245             :
         246             :   ! Microphysical settings for the precipitation flux to mixing ratio conversion
         247             :   real(r8),parameter,dimension(nhydro) :: &
         248             :                  !  LSL     LSI         LSR         LSS       CVL     CVI        CVR          CVS          LSG
         249             :        N_ax    = (/-1._r8, -1._r8,     8.e6_r8,     3.e6_r8, -1._r8, -1._r8,     8.e6_r8,     3.e6_r8,     4.e6_r8/),&
         250             :        N_bx    = (/-1._r8, -1._r8,      0.0_r8,      0.0_r8, -1._r8, -1._r8,      0.0_r8,      0.0_r8,      0.0_r8/),&
         251             :        alpha_x = (/-1._r8, -1._r8,      0.0_r8,      0.0_r8, -1._r8, -1._r8,      0.0_r8,      0.0_r8,      0.0_r8/),&
         252             :        c_x     = (/-1._r8, -1._r8,    842.0_r8,     4.84_r8, -1._r8, -1._r8,    842.0_r8,     4.84_r8,     94.5_r8/),&
         253             :        d_x     = (/-1._r8, -1._r8,      0.8_r8,     0.25_r8, -1._r8, -1._r8,      0.8_r8,     0.25_r8,      0.5_r8/),&
         254             :        g_x     = (/-1._r8, -1._r8,      0.5_r8,      0.5_r8, -1._r8, -1._r8,      0.5_r8,      0.5_r8,      0.5_r8/),&
         255             :        a_x     = (/-1._r8, -1._r8,    524.0_r8,    52.36_r8, -1._r8, -1._r8,    524.0_r8,    52.36_r8,   209.44_r8/),&
         256             :        b_x     = (/-1._r8, -1._r8,      3.0_r8,      3.0_r8, -1._r8, -1._r8,      3.0_r8,      3.0_r8,      3.0_r8/),&
         257             :        gamma_1 = (/-1._r8, -1._r8, 17.83725_r8, 8.284701_r8, -1._r8, -1._r8, 17.83725_r8, 8.284701_r8, 11.63230_r8/),&
         258             :        gamma_2 = (/-1._r8, -1._r8,      6.0_r8,      6.0_r8, -1._r8, -1._r8,      6.0_r8,      6.0_r8,      6.0_r8/),&
         259             :        gamma_3 = (/-1._r8, -1._r8,      2.0_r8,      2.0_r8, -1._r8, -1._r8,      2.0_r8,      2.0_r8,      2.0_r8/),&
         260             :        gamma_4 = (/-1._r8, -1._r8,      6.0_r8,      6.0_r8, -1._r8, -1._r8,      6.0_r8,      6.0_r8,      6.0_r8/)
         261             : #endif
         262             :
         263             : CONTAINS
         264             :
         265             :   ! ######################################################################################
         266             :   ! SUBROUTINE cospsimulator_intr_readnl
         267             :   ! ######################################################################################
         268        1536 :   subroutine cospsimulator_intr_readnl(nlfile)
         269             :     use namelist_utils,  only: find_group_name
         270             :     use units,           only: getunit, freeunit
         271             : #ifdef SPMD
         272             :     use mpishorthand,    only: mpicom, mpilog, mpiint
         273             : #endif
         274             :
         275             :     character(len=*), intent(in) :: nlfile  ! file containing namelist input  (nlfile=atm_in)
         276             :
         277             :     ! Local variables
         278             :     integer :: unitn, ierr
         279             :     character(len=*), parameter :: subname = 'cospsimulator_intr_readnl'
         280             :
         281             : #ifdef USE_COSP
         282             :     namelist /cospsimulator_nl/ docosp, cosp_ncolumns, cosp_nradsteps,           &
         283             :        cosp_amwg, cosp_lite, cosp_passive, cosp_active, cosp_isccp, cosp_runall, &
         284             :        cosp_lfrac_out, cosp_lradar_sim, cosp_llidar_sim, cosp_lisccp_sim,        &
         285             :        cosp_lmisr_sim, cosp_lmodis_sim,                                          &
         286             :        cosp_histfile_num, cosp_histfile_aux, cosp_histfile_aux_num
         287             :
         288             :     !! read in the namelist
         289             :     if (masterproc) then
         290             :        unitn = getunit()
         291             :        open( unitn, file=trim(nlfile), status='old' )
         292             :        call find_group_name(unitn, 'cospsimulator_nl', status=ierr)
         293             :        if (ierr == 0) then
         294             :           read(unitn, cospsimulator_nl, iostat=ierr)
         295             :           if (ierr /= 0) then
         296             :              call endrun(subname // ':: ERROR reading namelist')
         297             :           end if
         298             :        end if
         299             :        close(unitn)
         300             :        call freeunit(unitn)
         301             :     end if
         302             :
         303             : #ifdef SPMD
         304             :     ! Broadcast namelist variables
         305             :     call mpibcast(docosp,               1,  mpilog, 0, mpicom)
         306             :     call mpibcast(cosp_amwg,            1,  mpilog, 0, mpicom)
         307             :     call mpibcast(cosp_lite,            1,  mpilog, 0, mpicom)
         308             :     call mpibcast(cosp_passive,         1,  mpilog, 0, mpicom)
         309             :     call mpibcast(cosp_active,          1,  mpilog, 0, mpicom)
         310             :     call mpibcast(cosp_isccp,           1,  mpilog, 0, mpicom)
         311             :     call mpibcast(cosp_runall,          1,  mpilog, 0, mpicom)
         312             :     call mpibcast(cosp_lfrac_out,       1,  mpilog, 0, mpicom)
         313             :     call mpibcast(cosp_lradar_sim,      1,  mpilog, 0, mpicom)
         314             :     call mpibcast(cosp_llidar_sim,      1,  mpilog, 0, mpicom)
         315             :     call mpibcast(cosp_lisccp_sim,      1,  mpilog, 0, mpicom)
         316             :     call mpibcast(cosp_lmisr_sim,       1,  mpilog, 0, mpicom)
         317             :     call mpibcast(cosp_lmodis_sim,      1,  mpilog, 0, mpicom)
         318             :     call mpibcast(cosp_ncolumns,        1,  mpiint, 0, mpicom)
         319             :     call mpibcast(cosp_histfile_num,    1,  mpiint, 0, mpicom)
         320             :     call mpibcast(cosp_histfile_aux_num,1,  mpiint, 0, mpicom)
         321             :     call mpibcast(cosp_histfile_aux,    1,  mpilog, 0, mpicom)
         322             :     call mpibcast(cosp_nradsteps,       1,  mpiint, 0, mpicom)
         323             : #endif
         324             :
         325             :     if (cosp_lfrac_out) then
         326             :        lfrac_out = .true.
         327             :     end if
         328             :     if (cosp_lradar_sim) then
         329             :        lradar_sim = .true.
         330             :     end if
         331             :     if (cosp_llidar_sim) then
         332             :        llidar_sim = .true.
         333             :        lparasol_sim = .true.
         334             :     end if
         335             :     if (cosp_lisccp_sim) then
         336             :        lisccp_sim = .true.
         337             :     end if
         338             :     if (cosp_lmisr_sim) then
         339             :        lmisr_sim = .true.
         340             :     end if
         341             :     if (cosp_lmodis_sim) then
         342             :        lmodis_sim = .true.
         343             :     end if
         344             :
         345             :     if (cosp_histfile_aux .and. cosp_histfile_aux_num == -1) then
         346             :        cosp_histfile_aux_num = cosp_histfile_num
         347             :     end if
         348             :
         349             :     if (cosp_lite) then
         350             :        llidar_sim = .true.
         351             :        lparasol_sim = .true.
         352             :        lisccp_sim = .true.
         353             :        lmisr_sim = .true.
         354             :        lmodis_sim = .true.
         355             :        cosp_ncolumns = 10
         356             :        cosp_nradsteps = 3
         357             :     end if
         358             :
         359             :     if (cosp_passive) then
         360             :        lisccp_sim = .true.
         361             :        lmisr_sim = .true.
         362             :        lmodis_sim = .true.
         363             :        cosp_ncolumns = 10
         364             :        cosp_nradsteps = 3
         365             :     end if
         366             :
         367             :     if (cosp_active) then
         368             :        lradar_sim = .true.
         369             :        llidar_sim = .true.
         370             :        lparasol_sim = .true.
         371             :        cosp_ncolumns = 10
         372             :        cosp_nradsteps = 3
         373             :     end if
         374             :
         375             :     if (cosp_isccp) then
         376             :        lisccp_sim = .true.
         377             :        cosp_ncolumns = 10
         378             :        cosp_nradsteps = 3
         379             :     end if
         380             :
         381             :     if (cosp_runall) then
         382             :        lradar_sim = .true.
         383             :        llidar_sim = .true.
         384             :        lparasol_sim = .true.
         385             :        lisccp_sim = .true.
         386             :        lmisr_sim = .true.
         387             :        lmodis_sim = .true.
         388             :        lfrac_out = .true.
         389             :     end if
         390             :
         391             :     !! if no simulators are turned on at all and docosp is, set cosp_amwg = .true.
         392             :     if((docosp) .and. (.not.lradar_sim) .and. (.not.llidar_sim) .and. (.not.lisccp_sim) .and. &
         393             :          (.not.lmisr_sim) .and. (.not.lmodis_sim)) then
         394             :        cosp_amwg = .true.
         395             :     end if
         396             :     if (cosp_amwg) then
         397             :        lradar_sim = .true.
         398             :        llidar_sim = .true.
         399             :        lparasol_sim = .true.
         400             :        lisccp_sim = .true.
         401             :        lmisr_sim = .true.
         402             :        lmodis_sim = .true.
         403             :        cosp_ncolumns = 10
         404             :        cosp_nradsteps = 3
         405             :     end if
         406             :
         407             :     ! Set number of sub-columns, from namelist
         408             :     ncolumns   = cosp_ncolumns
         409             :     nscol_cosp = cosp_ncolumns
         410             :
         411             :     if (masterproc) then
         412             :        if (docosp) then
         413             :           write(iulog,*)'COSP configuration:'
         414             :           write(iulog,*)'  Number of COSP subcolumns                = ', cosp_ncolumns
         415             :           write(iulog,*)'  COSP frequency in radiation steps        = ', cosp_nradsteps
         416             :           write(iulog,*)'  Enable radar simulator                   = ', lradar_sim
         417             :           write(iulog,*)'  Enable calipso simulator                 = ', llidar_sim
         418             :           write(iulog,*)'  Enable ISCCP simulator                   = ', lisccp_sim
         419             :           write(iulog,*)'  Enable MISR simulator                    = ', lmisr_sim
         420             :           write(iulog,*)'  Enable MODIS simulator                   = ', lmodis_sim
         421             :           write(iulog,*)'  RADAR_SIM microphysics scheme            = ', trim(cloudsat_micro_scheme)
         422             :           write(iulog,*)'  Write COSP output to history file        = ', cosp_histfile_num
         423             :           write(iulog,*)'  Write COSP input fields                  = ', cosp_histfile_aux
         424             :           write(iulog,*)'  Write COSP input fields to history file  = ', cosp_histfile_aux_num
         425             :           write(iulog,*)'  Write COSP subcolumn fields              = ', lfrac_out
         426             :        else
         427             :           write(iulog,*)'COSP not enabled'
         428             :        end if
         429             :     end if
         430             : #endif
         431        1536 :   end subroutine cospsimulator_intr_readnl
         432             :
         433             :   ! ######################################################################################
         434             :   ! SUBROUTINE cospsimulator_intr_register
         435             :   ! ######################################################################################
         436        1536 :   subroutine cospsimulator_intr_register()
         437             :
         438             :     ! The coordinate variables used for COSP output are defined here.  This
         439             :     ! needs to be done before the call to read_restart_history in order for
         440             :     ! restarts to work.
         441             :
         442             :     use cam_history_support, only: add_hist_coord
         443             :     !---------------------------------------------------------------------------
         444             :
         445             : #ifdef USE_COSP
         446             :     ! Set number of levels used by COSP to the number of levels used by
         447             :     ! CAM's cloud macro/microphysics parameterizations.
         448             :     nlay = pver - ktop + 1
         449             :     nlayp = nlay + 1
         450             :
         451             :     ! Set COSP coordinate arrays
         452             :     call setcosp2values()
         453             :
         454             :     ! Define coordinate variables for COSP outputs.
         455             :     if (lisccp_sim .or. lmodis_sim) then
         456             :        call add_hist_coord('cosp_prs', nprs_cosp, 'COSP Mean ISCCP pressure',  &
         457             :             'hPa', prsmid_cosp, bounds_name='cosp_prs_bnds', bounds=prslim_cosp)
         458             :     end if
         459             :
         460             :     if (lisccp_sim .or. lmisr_sim) then
         461             :        call add_hist_coord('cosp_tau', ntau_cosp,                              &
         462             :             'COSP Mean ISCCP optical depth', '1', taumid_cosp,                 &
         463             :             bounds_name='cosp_tau_bnds', bounds=taulim_cosp)
         464             :     end if
         465             :
         466             :     if (lisccp_sim .or. llidar_sim .or. lradar_sim .or. lmisr_sim) then
         467             :        call add_hist_coord('cosp_scol', nscol_cosp, 'COSP subcolumn',          &
         468             :             values=scol_cosp)
         469             :     end if
         470             :
         471             :     if (llidar_sim .or. lradar_sim) then
         472             :        call add_hist_coord('cosp_ht', nht_cosp,                                &
         473             :             'COSP Mean Height for calipso and radar simulator outputs', 'm',   &
         474             :             htmid_cosp, bounds_name='cosp_ht_bnds', bounds=htlim_cosp,         &
         475             :             vertical_coord=.true.)
         476             :     end if
         477             :
         478             :     if (llidar_sim) then
         479             :        call add_hist_coord('cosp_sr', nsr_cosp,                                &
         480             :             'COSP Mean Scattering Ratio for calipso simulator CFAD output', '1', &
         481             :             srmid_cosp, bounds_name='cosp_sr_bnds', bounds=srlim_cosp)
         482             :     end if
         483             :
         484             :     if (llidar_sim) then
         485             :        call add_hist_coord('cosp_sza', nsza_cosp, 'COSP Parasol SZA',          &
         486             :             'degrees', sza_cosp)
         487             :     end if
         488             :
         489             :     if (lradar_sim) then
         490             :        call add_hist_coord('cosp_dbze', CLOUDSAT_DBZE_BINS,                    &
         491             :             'COSP Mean dBZe for radar simulator CFAD output', 'dBZ',           &
         492             :             dbzemid_cosp, bounds_name='cosp_dbze_bnds', bounds=dbzelim_cosp)
         493             :     end if
         494             :
         495             :     if (lmisr_sim) then
         496             :        call add_hist_coord('cosp_htmisr', nhtmisr_cosp, 'COSP MISR height', &
         497             :             'km', htmisrmid_cosp,                                           &
         498             :             bounds_name='cosp_htmisr_bnds', bounds=htmisrlim_cosp)
         499             :     end if
         500             :
         501             :     if (lmodis_sim) then
         502             :        call add_hist_coord('cosp_tau_modis', ntau_cosp_modis,                  &
         503             :             'COSP Mean MODIS optical depth', '1', taumid_cosp_modis,           &
         504             :             bounds_name='cosp_tau_modis_bnds', bounds=taulim_cosp_modis)
         505             :        call add_hist_coord('cosp_reffice',numMODISReffIceBins,                 &
         506             :             'COSP Mean MODIS effective radius (ice)', 'microns', reffICE_binCenters_cosp, &
         507             :             bounds_name='cosp_reffice_bnds',bounds=reffICE_binEdges_cosp)
         508             :        call add_hist_coord('cosp_reffliq',numMODISReffLiqBins,                 &
         509             :             'COSP Mean MODIS effective radius (liquid)', 'microns', reffLIQ_binCenters_cosp, &
         510             :             bounds_name='cosp_reffliq_bnds',bounds=reffLIQ_binEdges_cosp)
         511             :     end if
         512             :
         513             : #endif
         514        1536 :   end subroutine cospsimulator_intr_register
         515             :
         516             :   ! ######################################################################################
         517             :   ! SUBROUTINE cospsimulator_intr_init
         518             :   ! ######################################################################################
         519           0 :   subroutine cospsimulator_intr_init()
         520             :
         521             : #ifdef USE_COSP
         522             :
         523             :     use cam_history,         only: addfld, add_default, horiz_only
         524             :     use physics_buffer,      only: pbuf_get_index
         525             :
         526             :     integer :: i, ierr, istat
         527             :     character(len=*), parameter :: sub = 'cospsimulator_intr_init'
         528             :     !---------------------------------------------------------------------------
         529             :
         530             :     ! The COSP init method (setcosp2values) was run from cospsimulator_intr_register in order to add
         531             :     ! the history coordinate variables earlier as needed for the restart time sequencing.
         532             :
         533             :     ! ISCCP OUTPUTS
         534             :     if (lisccp_sim) then
         535             :        call addfld('FISCCP1_COSP', (/'cosp_tau','cosp_prs'/), 'A', 'percent', &
         536             :             'Grid-box fraction covered by each ISCCP D level cloud type', flag_xyfill=.true., fill_value=R_UNDEF)
         537             :        call addfld('CLDTOT_ISCCP', horiz_only, 'A', 'percent', &
         538             :             'Total Cloud Fraction Calculated by the ISCCP Simulator ',flag_xyfill=.true., fill_value=R_UNDEF)
         539             :        call addfld('MEANCLDALB_ISCCP', horiz_only, 'A', '1', &
         540             :             'Mean cloud albedo*CLDTOT_ISCCP', flag_xyfill=.true., fill_value=R_UNDEF)
         541             :        call addfld('MEANPTOP_ISCCP', horiz_only, 'A', 'Pa', &
         542             :             'Mean cloud top pressure*CLDTOT_ISCCP',flag_xyfill=.true., fill_value=R_UNDEF)
         543             :        call addfld('MEANTAU_ISCCP', horiz_only, 'A', '1', &
         544             :             'Mean optical thickness*CLDTOT_ISCCP',flag_xyfill=.true., fill_value=R_UNDEF)
         545             :        call addfld('MEANTB_ISCCP', horiz_only, 'A', 'K', &
         546             :             'Mean Infrared Tb from ISCCP simulator',flag_xyfill=.true., fill_value=R_UNDEF)
         547             :        call addfld('MEANTBCLR_ISCCP', horiz_only, 'A', 'K', &
         548             :             'Mean Clear-sky Infrared Tb from ISCCP simulator', flag_xyfill=.true., fill_value=R_UNDEF)
         549             :        call addfld('TAU_ISCCP', (/'cosp_scol'/), 'I', '1', &
         550             :             'Optical Depth in each Subcolumn', flag_xyfill=.true., fill_value=R_UNDEF)
         551             :        call addfld('CLDPTOP_ISCCP', (/'cosp_scol'/), 'I', 'Pa', &
         552             :             'Cloud Top Pressure in each Subcolumn', flag_xyfill=.true., fill_value=R_UNDEF)
         553             :
         554             :        call add_default('FISCCP1_COSP',cosp_histfile_num,' ')
         555             :        call add_default('CLDTOT_ISCCP',cosp_histfile_num,' ')
         556             :        call add_default('MEANCLDALB_ISCCP',cosp_histfile_num,' ')
         557             :        call add_default('MEANPTOP_ISCCP',cosp_histfile_num,' ')
         558             :        call add_default('MEANTAU_ISCCP',cosp_histfile_num,' ')
         559             :        call add_default('MEANTB_ISCCP',cosp_histfile_num,' ')
         560             :        call add_default('MEANTBCLR_ISCCP',cosp_histfile_num,' ')
         561             :
         562             :     end if
         563             :
         564             :     ! CALIPSO SIMULATOR OUTPUTS
         565             :     if (llidar_sim) then
         566             :        call addfld('CLDLOW_CAL', horiz_only, 'A', 'percent', &
         567             :             'Calipso Low-level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         568             :        call addfld('CLDMED_CAL', horiz_only, 'A', 'percent', &
         569             :             'Calipso Mid-level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         570             :        call addfld('CLDHGH_CAL', horiz_only, 'A', 'percent', &
         571             :             'Calipso High-level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         572             :        call addfld('CLDTOT_CAL', horiz_only, 'A', 'percent', &
         573             :             'Calipso Total Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         574             :        call addfld('CLD_CAL', (/'cosp_ht'/), 'A', 'percent', &
         575             :             'Calipso Cloud Fraction (532 nm)', flag_xyfill=.true., fill_value=R_UNDEF)
         576             :        call addfld('RFL_PARASOL', (/'cosp_sza'/), 'A', 'fraction', &
         577             :             'PARASOL-like mono-directional reflectance ', flag_xyfill=.true., fill_value=R_UNDEF)
         578             :        call addfld('CFAD_SR532_CAL', (/'cosp_sr','cosp_ht'/), 'A', 'fraction', &
         579             :             'Calipso Scattering Ratio CFAD (532 nm)', flag_xyfill=.true., fill_value=R_UNDEF)
         580             :        call addfld('MOL532_CAL', (/'trop_pref'/), 'A', 'm-1 sr-1', &
         581             :             'Calipso Molecular Backscatter (532 nm) ', flag_xyfill=.true., fill_value=R_UNDEF)
         582             :        call addfld('ATB532_CAL', (/'cosp_scol','trop_pref'/), 'I', 'no_unit_log10(x)', &
         583             :             'Calipso Attenuated Total Backscatter (532 nm) in each Subcolumn', flag_xyfill=.true., fill_value=R_UNDEF)
         584             :        call addfld('CLD_CAL_LIQ', (/'cosp_ht'/), 'A', 'percent', &
         585             :             'Calipso Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         586             :        call addfld('CLD_CAL_ICE', (/'cosp_ht'/), 'A', 'percent', &
         587             :             'Calipso Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         588             :        call addfld('CLD_CAL_UN', (/'cosp_ht'/), 'A', 'percent',  &
         589             :             'Calipso Undefined-Phase Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         590             :        call addfld('CLD_CAL_TMP', (/'cosp_ht'/), 'A', 'K', &
         591             :             'Calipso Cloud Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         592             :        call addfld('CLD_CAL_TMPLIQ', (/'cosp_ht'/), 'A', 'K', &
         593             :             'Calipso Liquid Cloud Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         594             :        call addfld('CLD_CAL_TMPICE', (/'cosp_ht'/), 'A', 'K', &
         595             :             'Calipso Ice Cloud Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         596             :        call addfld('CLD_CAL_TMPUN', (/'cosp_ht'/), 'A', 'K', &
         597             :             'Calipso Undefined-Phase Cloud Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         598             :        call addfld('CLDTOT_CAL_ICE', horiz_only, 'A', 'percent', &
         599             :             'Calipso Total Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         600             :        call addfld('CLDTOT_CAL_LIQ', horiz_only, 'A', 'percent', &
         601             :             'Calipso Total Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         602             :        call addfld('CLDTOT_CAL_UN', horiz_only, 'A', 'percent', &
         603             :             'Calipso Total Undefined-Phase Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         604             :        call addfld('CLDHGH_CAL_ICE', horiz_only, 'A', 'percent', &
         605             :             'Calipso High-level Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         606             :        call addfld('CLDHGH_CAL_LIQ', horiz_only, 'A', 'percent', &
         607             :             'Calipso High-level Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         608             :        call addfld('CLDHGH_CAL_UN', horiz_only, 'A', 'percent', &
         609             :             'Calipso High-level Undefined-Phase Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         610             :        call addfld('CLDMED_CAL_ICE', horiz_only, 'A', 'percent', &
         611             :             'Calipso Mid-level Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         612             :        call addfld('CLDMED_CAL_LIQ', horiz_only, 'A', 'percent', &
         613             :             'Calipso Mid-level Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         614             :        call addfld('CLDMED_CAL_UN', horiz_only, 'A', 'percent', &
         615             :             'Calipso Mid-level Undefined-Phase Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         616             :        call addfld('CLDLOW_CAL_ICE', horiz_only, 'A', 'percent', &
         617             :             'Calipso Low-level Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         618             :        call addfld('CLDLOW_CAL_LIQ', horiz_only, 'A', 'percent', &
         619             :             'Calipso Low-level Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         620             :        call addfld('CLDLOW_CAL_UN', horiz_only, 'A', 'percent', &
         621             :             'Calipso Low-level Undefined-Phase Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         622             :
         623             :        call add_default('CLDLOW_CAL',cosp_histfile_num,' ')
         624             :        call add_default('CLDMED_CAL',cosp_histfile_num,' ')
         625             :        call add_default('CLDHGH_CAL',cosp_histfile_num,' ')
         626             :        call add_default('CLDTOT_CAL',cosp_histfile_num,' ')
         627             :        call add_default('CLD_CAL',cosp_histfile_num,' ')
         628             :        call add_default('RFL_PARASOL',cosp_histfile_num,' ')
         629             :        call add_default('CFAD_SR532_CAL',cosp_histfile_num,' ')
         630             :        call add_default('CLD_CAL_LIQ',cosp_histfile_num,' ')
         631             :        call add_default('CLD_CAL_ICE',cosp_histfile_num,' ')
         632             :        call add_default('CLD_CAL_UN',cosp_histfile_num,' ')
         633             :        call add_default('CLDTOT_CAL_ICE',cosp_histfile_num,' ')
         634             :        call add_default('CLDTOT_CAL_LIQ',cosp_histfile_num,' ')
         635             :        call add_default('CLDTOT_CAL_UN',cosp_histfile_num,' ')
         636             :        call add_default('CLDHGH_CAL_ICE',cosp_histfile_num,' ')
         637             :        call add_default('CLDHGH_CAL_LIQ',cosp_histfile_num,' ')
         638             :        call add_default('CLDHGH_CAL_UN',cosp_histfile_num,' ')
         639             :        call add_default('CLDMED_CAL_ICE',cosp_histfile_num,' ')
         640             :        call add_default('CLDMED_CAL_LIQ',cosp_histfile_num,' ')
         641             :        call add_default('CLDMED_CAL_UN',cosp_histfile_num,' ')
         642             :        call add_default('CLDLOW_CAL_ICE',cosp_histfile_num,' ')
         643             :        call add_default('CLDLOW_CAL_LIQ',cosp_histfile_num,' ')
         644             :        call add_default('CLDLOW_CAL_UN',cosp_histfile_num,' ')
         645             :
         646             :        if ((.not.cosp_amwg) .and. (.not.cosp_lite) .and. (.not.cosp_passive) .and. (.not.cosp_active) &
         647             :             .and. (.not.cosp_isccp)) then
         648             :           call add_default('MOL532_CAL',cosp_histfile_num,' ')
         649             :        end if
         650             :     end if
         651             :
         652             :     ! RADAR SIMULATOR OUTPUTS
         653             :     allocate(sd_cs(begchunk:endchunk), rcfg_cs(begchunk:endchunk), stat=istat)
         654             :     call handle_allocate_error(istat, sub, 'sd_cs,rcfg_cs')
         655             :     if (lradar_sim) then
         656             :
         657             :        do i = begchunk, endchunk
         658             :           sd_cs(i)   = sd
         659             :           rcfg_cs(i) = rcfg_cloudsat
         660             :        end do
         661             :
         662             :        call addfld('CFAD_DBZE94_CS',(/'cosp_dbze','cosp_ht  '/), 'A', 'fraction', &
         663             :             'Radar Reflectivity Factor CFAD (94 GHz)', flag_xyfill=.true., fill_value=R_UNDEF)
         664             :        call addfld('CLD_CAL_NOTCS', (/'cosp_ht'/), 'A', 'percent', &
         665             :             'Cloud occurrence seen by CALIPSO but not CloudSat ', flag_xyfill=.true., fill_value=R_UNDEF)
         666             :        call addfld('CLDTOT_CALCS', horiz_only, 'A', 'percent', &
         667             :             'Calipso and Radar Total Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         668             :        call addfld('CLDTOT_CS', horiz_only, 'A', 'percent', &
         669             :             'Radar total cloud amount', flag_xyfill=.true., fill_value=R_UNDEF)
         670             :        call addfld('CLDTOT_CS2', horiz_only, 'A', 'percent', &
         671             :             'Radar total cloud amount without the data for the first kilometer above surface ', &
         672             :             flag_xyfill=.true., fill_value=R_UNDEF)
         673             :        call addfld('DBZE_CS', (/'cosp_scol','trop_pref'/), 'I', 'dBZe', &
         674             :             'Radar dBZe (94 GHz) in each Subcolumn', flag_xyfill=.true., fill_value=R_UNDEF)
         675             :
         676             :        ! Cloudsat near-sfc precipitation diagnostics
         677             :        call addfld('CS_NOPRECIP',  horiz_only, 'A', '1', &
         678             :             'CloudSat No Rain Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         679             :        call addfld('CS_RAINPOSS',  horiz_only, 'A', '1', &
         680             :             'Cloudsat Rain Possible Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         681             :        call addfld('CS_RAINPROB',  horiz_only, 'A', '1', &
         682             :             'CloudSat Rain Probable Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         683             :        call addfld('CS_RAINCERT',  horiz_only, 'A', '1', &
         684             :             'CloudSat Rain Certain Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         685             :        call addfld('CS_SNOWPOSS',  horiz_only, 'A', '1', &
         686             :             'CloudSat Snow Possible Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         687             :        call addfld('CS_SNOWCERT',  horiz_only, 'A', '1', &
         688             :             'CloudSat Snow Certain Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         689             :        call addfld('CS_MIXPOSS',   horiz_only, 'A', '1', &
         690             :             'CloudSat Mixed Possible Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         691             :        call addfld('CS_MIXCERT',   horiz_only, 'A', '1', &
         692             :             'CloudSat Mixed Certain Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         693             :        call addfld('CS_RAINHARD',  horiz_only, 'A', '1', &
         694             :             'CloudSat Heavy Rain Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         695             :        call addfld('CS_UN',        horiz_only, 'A', '1', &
         696             :             'CloudSat Unclassified Precipitation Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         697             :        call addfld('CS_PIA',       horiz_only, 'A', 'dBZ', &
         698             :             'CloudSat Radar Path Integrated Attenuation', flag_xyfill=.true., fill_value=R_UNDEF)
         699             :
         700             :        call add_default('CFAD_DBZE94_CS',cosp_histfile_num,' ')
         701             :        call add_default('CLD_CAL_NOTCS', cosp_histfile_num,' ')
         702             :        call add_default('CLDTOT_CALCS',  cosp_histfile_num,' ')
         703             :        call add_default('CLDTOT_CS',     cosp_histfile_num,' ')
         704             :        call add_default('CLDTOT_CS2',    cosp_histfile_num,' ')
         705             :        call add_default('CS_NOPRECIP',   cosp_histfile_num,' ')
         706             :        call add_default('CS_RAINPOSS',   cosp_histfile_num,' ')
         707             :        call add_default('CS_RAINPROB',   cosp_histfile_num,' ')
         708             :        call add_default('CS_RAINCERT',   cosp_histfile_num,' ')
         709             :        call add_default('CS_SNOWPOSS',   cosp_histfile_num,' ')
         710             :        call add_default('CS_SNOWCERT',   cosp_histfile_num,' ')
         711             :        call add_default('CS_MIXPOSS',    cosp_histfile_num,' ')
         712             :        call add_default('CS_MIXCERT',    cosp_histfile_num,' ')
         713             :        call add_default('CS_RAINHARD',   cosp_histfile_num,' ')
         714             :        call add_default('CS_UN',         cosp_histfile_num,' ')
         715             :        call add_default('CS_PIA',        cosp_histfile_num,' ')
         716             :     end if
         717             :
         718             :     ! MISR SIMULATOR OUTPUTS
         719             :     if (lmisr_sim) then
         720             :        call addfld('CLD_MISR', (/'cosp_tau   ','cosp_htmisr'/), 'A', 'percent', &
         721             :             'Cloud Fraction from MISR Simulator', flag_xyfill=.true., fill_value=R_UNDEF)
         722             :
         723             :        call add_default('CLD_MISR',cosp_histfile_num,' ')
         724             :     end if
         725             :
         726             :     ! MODIS OUTPUT
         727             :     if (lmodis_sim) then
         728             :        call addfld('CLTMODIS', horiz_only, 'A', '%', &
         729             :             'MODIS Total Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         730             :        call addfld('CLWMODIS', horiz_only, 'A', '%', &
         731             :             'MODIS Liquid Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         732             :        call addfld('CLIMODIS', horiz_only, 'A', '%', &
         733             :             'MODIS Ice Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         734             :        call addfld('CLHMODIS', horiz_only, 'A', '%', &
         735             :             'MODIS High Level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         736             :        call addfld('CLMMODIS', horiz_only, 'A', '%', &
         737             :             'MODIS Mid Level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         738             :        call addfld('CLLMODIS', horiz_only, 'A', '%', &
         739             :             'MODIS Low Level Cloud Fraction', flag_xyfill=.true., fill_value=R_UNDEF)
         740             :        call addfld('TAUTMODIS', horiz_only, 'A', '1', &
         741             :             'MODIS Total Cloud Optical Thickness*CLTMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         742             :        call addfld('TAUWMODIS', horiz_only, 'A', '1', &
         743             :             'MODIS Liquid Cloud Optical Thickness*CLWMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         744             :        call addfld('TAUIMODIS', horiz_only, 'A', '1', &
         745             :             'MODIS Ice Cloud Optical Thickness*CLIMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         746             :        call addfld('TAUTLOGMODIS', horiz_only, 'A', '1', &
         747             :             'MODIS Total Cloud Optical Thickness (Log10 Mean)*CLTMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         748             :        call addfld('TAUWLOGMODIS', horiz_only, 'A', '1', &
         749             :             'MODIS Liquid Cloud Optical Thickness (Log10 Mean)*CLWMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         750             :        call addfld('TAUILOGMODIS', horiz_only, 'A', '1', &
         751             :             'MODIS Ice Cloud Optical Thickness (Log10 Mean)*CLIMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         752             :        call addfld('REFFCLWMODIS', horiz_only, 'A', 'm', &
         753             :             'MODIS Liquid Cloud Particle Size*CLWMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         754             :        call addfld('REFFCLIMODIS', horiz_only, 'A', 'm', &
         755             :             'MODIS Ice Cloud Particle Size*CLIMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         756             :        call addfld('PCTMODIS', horiz_only, 'A', 'Pa', &
         757             :             'MODIS Cloud Top Pressure*CLTMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         758             :        call addfld('LWPMODIS', horiz_only, 'A', 'kg m-2', &
         759             :             'MODIS Cloud Liquid Water Path*CLWMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         760             :        call addfld('IWPMODIS', horiz_only, 'A', 'kg m-2', &
         761             :             'MODIS Cloud Ice Water Path*CLIMODIS', flag_xyfill=.true., fill_value=R_UNDEF)
         762             :        call addfld('CLMODIS', (/'cosp_tau_modis','cosp_prs      '/), 'A', '%', &
         763             :             'MODIS Cloud Area Fraction (tau-pressure histogram)', flag_xyfill=.true., fill_value=R_UNDEF)
         764             :        call addfld('CLRIMODIS', (/'cosp_tau_modis','cosp_reffice  '/), 'A', '%', &
         765             :             'MODIS Cloud Area Fraction (tau-reffice histogram)', flag_xyfill=.true., fill_value=R_UNDEF)
         766             :        call addfld('CLRLMODIS', (/'cosp_tau_modis','cosp_reffliq  '/), 'A', '%', &
         767             :             'MODIS Cloud Area Fraction (tau-reffliq histogram)', flag_xyfill=.true., fill_value=R_UNDEF)
         768             :
         769             :        call add_default('CLTMODIS',cosp_histfile_num,' ')
         770             :        call add_default('CLWMODIS',cosp_histfile_num,' ')
         771             :        call add_default('CLIMODIS',cosp_histfile_num,' ')
         772             :        call add_default('CLHMODIS',cosp_histfile_num,' ')
         773             :        call add_default('CLMMODIS',cosp_histfile_num,' ')
         774             :        call add_default('CLLMODIS',cosp_histfile_num,' ')
         775             :        call add_default('TAUTMODIS',cosp_histfile_num,' ')
         776             :        call add_default('TAUWMODIS',cosp_histfile_num,' ')
         777             :        call add_default('TAUIMODIS',cosp_histfile_num,' ')
         778             :        call add_default('TAUTLOGMODIS',cosp_histfile_num,' ')
         779             :        call add_default('TAUWLOGMODIS',cosp_histfile_num,' ')
         780             :        call add_default('TAUILOGMODIS',cosp_histfile_num,' ')
         781             :        call add_default('REFFCLWMODIS',cosp_histfile_num,' ')
         782             :        call add_default('REFFCLIMODIS',cosp_histfile_num,' ')
         783             :        call add_default('PCTMODIS',cosp_histfile_num,' ')
         784             :        call add_default('LWPMODIS',cosp_histfile_num,' ')
         785             :        call add_default('IWPMODIS',cosp_histfile_num,' ')
         786             :        call add_default('CLMODIS',cosp_histfile_num,' ')
         787             :        call add_default('CLRIMODIS',cosp_histfile_num,' ')
         788             :        call add_default('CLRLMODIS',cosp_histfile_num,' ')
         789             :     end if
         790             :
         791             :     ! SUB-COLUMN OUTPUT
         792             :     if (lfrac_out) then
         793             :        call addfld('SCOPS_OUT', (/'cosp_scol','trop_pref'/), 'I', '0=nocld,1=strcld,2=cnvcld', &
         794             :             'SCOPS Subcolumn output', flag_xyfill=.true., fill_value=R_UNDEF)
         795             :
         796             :        call add_default('SCOPS_OUT',cosp_histfile_num,' ')
         797             :
         798             :        if (lisccp_sim) then
         799             :           call add_default('TAU_ISCCP',cosp_histfile_num,' ')
         800             :           call add_default('CLDPTOP_ISCCP',cosp_histfile_num,' ')
         801             :        end if
         802             :
         803             :        if (llidar_sim) then
         804             :           call add_default('ATB532_CAL',cosp_histfile_num,' ')
         805             :        end if
         806             :
         807             :        if (lradar_sim) then
         808             :           call add_default('DBZE_CS',cosp_histfile_num,' ')
         809             :        end if
         810             :     end if
         811             :
         812             :     !! ADDFLD, ADD_DEFAULT, OUTFLD CALLS FOR COSP OUTPUTS IF RUNNING COSP OFF-LINE
         813             :     if (cosp_histfile_aux) then
         814             :        call addfld ('PS_COSP',         horiz_only,            'I','Pa', &
         815             :           'COSP Surface Pressure', flag_xyfill=.true., fill_value=R_UNDEF)
         816             :        call addfld ('TS_COSP',         horiz_only,            'I','K',  &
         817             :           'COSP Skin Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         818             :        call addfld ('P_COSP',          (/            'trop_pref'/), 'I','Pa', &
         819             :           'COSP Pressure (layer midpoint)', flag_xyfill=.true., fill_value=R_UNDEF)
         820             :        call addfld ('PH_COSP',         (/            'trop_prefi'/), 'I','Pa', &
         821             :           'COSP Pressure (layer interface)', flag_xyfill=.true., fill_value=R_UNDEF)
         822             :        call addfld ('ZLEV_COSP',       (/            'trop_pref'/), 'I','m',  &
         823             :           'COSP Height (layer midpoint)', flag_xyfill=.true., fill_value=R_UNDEF)
         824             :        call addfld ('ZLEV_HALF_COSP',  (/            'trop_prefi'/), 'I','m',  &
         825             :           'COSP Height (layer interface)', flag_xyfill=.true., fill_value=R_UNDEF)
         826             :        call addfld ('T_COSP',          (/            'trop_pref'/), 'I','K',  &
         827             :           'COSP Temperature', flag_xyfill=.true., fill_value=R_UNDEF)
         828             :        call addfld ('Q_COSP',          (/            'trop_pref'/), 'I','percent', &
         829             :           'COSP Specific Humidity', flag_xyfill=.true., fill_value=R_UNDEF)
         830             :
         831             :        call addfld ('TAU_067',         (/'cosp_scol','trop_pref'/), 'I','1', &
         832             :           'Subcolumn 0.67micron optical depth', flag_xyfill=.true., fill_value=R_UNDEF)
         833             :        call addfld ('EMISS_11',        (/'cosp_scol','trop_pref'/), 'I','1', &
         834             :           'Subcolumn 11micron emissivity', flag_xyfill=.true., fill_value=R_UNDEF)
         835             :        call addfld ('MODIS_fracliq',   (/'cosp_scol','trop_pref'/), 'I','1', &
         836             :           'Fraction of tau from liquid water', flag_xyfill=.true., fill_value=R_UNDEF)
         837             :        call addfld ('MODIS_asym',      (/'cosp_scol','trop_pref'/), 'I','1', &
         838             :           'Asymmetry parameter (MODIS)', flag_xyfill=.true., fill_value=R_UNDEF)
         839             :        call addfld ('MODIS_ssa',       (/'cosp_scol','trop_pref'/), 'I','1', &
         840             :           'Single-scattering albedo (MODIS)', flag_xyfill=.true., fill_value=R_UNDEF)
         841             :        call addfld ('CS_z_vol',        (/'cosp_scol','trop_pref'/), 'I','1', &
         842             :           'Effective reflectivity factor (CLOUDSAT)', flag_xyfill=.true., fill_value=R_UNDEF)
         843             :        call addfld ('CS_kr_vol',       (/'cosp_scol','trop_pref'/), 'I','1', &
         844             :           'Attenuation coefficient (hydro) (CLOUDSAT)', flag_xyfill=.true., fill_value=R_UNDEF)
         845             :        call addfld ('CS_g_vol',        (/'cosp_scol','trop_pref'/), 'I','1', &
         846             :           'Attenuation coefficient (gases) (CLOUDSAT)', flag_xyfill=.true., fill_value=R_UNDEF)
         847             :
         848             :        call add_default('PS_COSP',         cosp_histfile_aux_num,' ')
         849             :        call add_default('TS_COSP',         cosp_histfile_aux_num,' ')
         850             :        call add_default('P_COSP',          cosp_histfile_aux_num,' ')
         851             :        call add_default('PH_COSP',         cosp_histfile_aux_num,' ')
         852             :        call add_default('ZLEV_COSP',       cosp_histfile_aux_num,' ')
         853             :        call add_default('ZLEV_HALF_COSP',  cosp_histfile_aux_num,' ')
         854             :        call add_default('T_COSP',          cosp_histfile_aux_num,' ')
         855             :        call add_default('Q_COSP',          cosp_histfile_aux_num,' ')
         856             :        call add_default('TAU_067',         cosp_histfile_aux_num,' ')
         857             :        call add_default('EMISS_11',        cosp_histfile_aux_num,' ')
         858             :        call add_default('MODIS_fracliq',   cosp_histfile_aux_num,' ')
         859             :        call add_default('MODIS_asym',      cosp_histfile_aux_num,' ')
         860             :        call add_default('MODIS_ssa',       cosp_histfile_aux_num,' ')
         861             :        call add_default('CS_z_vol',        cosp_histfile_aux_num,' ')
         862             :        call add_default('CS_kr_vol',       cosp_histfile_aux_num,' ')
         863             :        call add_default('CS_g_vol',        cosp_histfile_aux_num,' ')
         864             :     end if
         865             :
         866             :     rei_idx        = pbuf_get_index('REI')
         867             :     rel_idx        = pbuf_get_index('REL')
         868             :     cld_idx        = pbuf_get_index('CLD')
         869             :     concld_idx     = pbuf_get_index('CONCLD')
         870             :     lsreffrain_idx = pbuf_get_index('LS_REFFRAIN')
         871             :     lsreffsnow_idx = pbuf_get_index('LS_REFFSNOW')
         872             :     cvreffliq_idx  = pbuf_get_index('CV_REFFLIQ')
         873             :     cvreffice_idx  = pbuf_get_index('CV_REFFICE')
         874             :     gb_totcldliqmr_idx = pbuf_get_index('GB_TOTCLDLIQMR') ! grid box total cloud liquid water mr (kg/kg)
         875             :     gb_totcldicemr_idx = pbuf_get_index('GB_TOTCLDICEMR') ! grid box total cloud ice water mr (kg/kg)
         876             :     dpflxprc_idx   = pbuf_get_index('DP_FLXPRC')
         877             :     dpflxsnw_idx   = pbuf_get_index('DP_FLXSNW')
         878             :     shflxprc_idx   = pbuf_get_index('SH_FLXPRC', errcode=ierr)
         879             :     shflxsnw_idx   = pbuf_get_index('SH_FLXSNW', errcode=ierr)
         880             :     lsflxprc_idx   = pbuf_get_index('LS_FLXPRC')
         881             :     lsflxsnw_idx   = pbuf_get_index('LS_FLXSNW')
         882             :
         883             :     allocate(first_run_cosp(begchunk:endchunk), run_cosp(1:pcols,begchunk:endchunk), &
         884             :        stat=istat)
         885             :     call handle_allocate_error(istat, sub, '*run_cosp')
         886             :     first_run_cosp(begchunk:endchunk)=.true.
         887             :     run_cosp(1:pcols,begchunk:endchunk)=.false.
         888             :
         889             : #endif
         890        1536 :   end subroutine cospsimulator_intr_init
         891             :
         892             :   ! ######################################################################################
         893             :   ! SUBROUTINE setcosp2values
         894             :   ! ######################################################################################
         895             : #ifdef USE_COSP
         896             :   subroutine setcosp2values()
         897             :     use mod_cosp,             only: cosp_init
         898             :     use mod_cosp_config,      only: vgrid_zl, vgrid_zu, vgrid_z
         899             :     use mod_quickbeam_optics, only: hydro_class_init, quickbeam_optics_init
         900             :
         901             :     ! Local
         902             :     logical :: ldouble=.false.
         903             :     logical :: lsingle=.true. ! Default is to use single moment
         904             :     integer :: k
         905             :     integer :: istat
         906             :     character(len=*), parameter :: sub = 'setcosp2values'
         907             :     !--------------------------------------------------------------------------------------
         908             :
         909             :     prsmid_cosp  = pres_binCenters
         910             :     prslim_cosp  = pres_binEdges
         911             :     taumid_cosp  = tau_binCenters
         912             :     taulim_cosp  = tau_binEdges
         913             :     srmid_cosp   = calipso_binCenters
         914             :     srlim_cosp   = calipso_binEdges
         915             :     sza_cosp     = parasol_sza
         916             :     dbzemid_cosp = cloudsat_binCenters
         917             :     dbzelim_cosp = cloudsat_binEdges
         918             :     htmisrmid_cosp = misr_histHgtCenters
         919             :     htmisrlim_cosp = misr_histHgtEdges
         920             :     taumid_cosp_modis = tau_binCenters
         921             :     taulim_cosp_modis = tau_binEdges
         922             :     reffICE_binCenters_cosp = reffICE_binCenters
         923             :     reffICE_binEdges_cosp   = reffICE_binEdges
         924             :     reffLIQ_binCenters_cosp = reffLIQ_binCenters
         925             :     reffLIQ_binEdges_cosp   = reffLIQ_binEdges
         926             :
         927             :     ! Initialize the distributional parameters for hydrometeors in radar simulator. In COSPv1.4, this was declared in
         928             :     ! cosp_defs.f.
         929             :     if (cloudsat_micro_scheme == 'MMF_v3.5_two_moment')  then
         930             :        ldouble = .true.
         931             :        lsingle = .false.
         932             :     endif
         933             :     call hydro_class_init(lsingle,ldouble,sd)
         934             :     call quickbeam_optics_init()
         935             :
         936             :     ! DS2017: The setting up of the vertical grid for regridding the CALIPSO and Cloudsat products is
         937             :     !         now done in cosp_init, but these fields are stored in cosp_config.F90.
         938             :     !         Additionally all static fields used by the individual simulators are set up by calls
         939             :     !         to _init functions in cosp_init.
         940             :     ! DS2019: Add logicals, default=.false., for new Lidar simuldators (Earthcare (atlid) and ground-based
         941             :     !         lidar at 532nm)
         942             :     call COSP_INIT(Lisccp_sim, Lmodis_sim, Lmisr_sim, Lradar_sim, Llidar_sim, LgrLidar532, &
         943             :          Latlid, Lparasol_sim, Lrttov_sim, radar_freq, k2, use_gas_abs, do_ray,              &
         944             :          isccp_topheight, isccp_topheight_direction, surface_radar, rcfg_cloudsat,           &
         945             :          use_vgrid, csat_vgrid, Nlr, nlay, cloudsat_micro_scheme)
         946             :
         947             :     if (use_vgrid) then      !! using fixed vertical grid
         948             :        if (csat_vgrid) then
         949             :           nht_cosp = 40
         950             :        else
         951             :           nht_cosp = Nlr
         952             :        endif
         953             :     endif
         954             :
         955             :     ! DJS2017: In COSP2, most of the bin boundaries, centers, and edges are declared in src/cosp_config.F90.
         956             :     !          Above I just assign them accordingly in the USE statement. Other bin bounds needed by CAM
         957             :     !          are calculated here.
         958             :
         959             :     allocate( &
         960             :        htmlmid_cosp(nlay), &
         961             :        htdbze_dbzemid_cosp(nht_cosp*CLOUDSAT_DBZE_BINS), &
         962             :        htlim_cosp(2,nht_cosp), &
         963             :        htmid_cosp(nht_cosp), &
         964             :        htlim_cosp_1d(nht_cosp+1), &
         965             :        htdbze_htmid_cosp(nht_cosp*CLOUDSAT_DBZE_BINS), &
         966             :        htsr_htmid_cosp(nht_cosp*nsr_cosp), &
         967             :        htsr_srmid_cosp(nht_cosp*nsr_cosp), &
         968             :        htmlscol_htmlmid_cosp(nlay*nscol_cosp), &
         969             :        htmlscol_scol_cosp(nlay*nscol_cosp), &
         970             :        scol_cosp(nscol_cosp), &
         971             :        htdbze_cosp(nht_cosp*CLOUDSAT_DBZE_BINS), &
         972             :        htsr_cosp(nht_cosp*nsr_cosp), &
         973             :        htmlscol_cosp(nlay*nscol_cosp), stat=istat)
         974             :     call handle_allocate_error(istat, sub, 'htmlmid_cosp,..,htmlscol_cosp')
         975             :
         976             :     ! DJS2017: Just pull from cosp_config
         977             :     if (use_vgrid) then
         978             :        htlim_cosp_1d(1)            = vgrid_zu(1)
         979             :        htlim_cosp_1d(2:nht_cosp+1) = vgrid_zl
         980             :     endif
         981             :     htmid_cosp      = vgrid_z
         982             :     htlim_cosp(1,:) = vgrid_zu
         983             :     htlim_cosp(2,:) = vgrid_zl
         984             :
         985             :     scol_cosp(:) = (/(k,k=1,nscol_cosp)/)
         986             :
         987             :     !  Just using an index here, model height is a prognostic variable
         988             :     htmlmid_cosp(:) = (/(k,k=1,nlay)/)
         989             :
         990             :     ! assign mixed dimensions an integer index for cam_history.F90
         991             :     do k=1,nprs_cosp*ntau_cosp
         992             :        prstau_cosp(k) = k
         993             :     end do
         994             :     do k=1,nprs_cosp*ntau_cosp_modis
         995             :        prstau_cosp_modis(k) = k
         996             :     end do
         997             :     do k=1,nht_cosp*CLOUDSAT_DBZE_BINS
         998             :        htdbze_cosp(k) = k
         999             :     end do
        1000             :     do k=1,nht_cosp*nsr_cosp
        1001             :        htsr_cosp(k) = k
        1002             :     end do
        1003             :     do k=1,nlay*nscol_cosp
        1004             :        htmlscol_cosp(k) = k
        1005             :     end do
        1006             :     do k=1,nhtmisr_cosp*ntau_cosp
        1007             :        htmisrtau_cosp(k) = k
        1008             :     end do
        1009             :
        1010             :     ! next, assign collapsed reference vectors for cam_history.F90
        1011             :     ! convention for saving output = prs1,tau1 ... prs1,tau7 ; prs2,tau1 ... prs2,tau7 etc.
        1012             :     ! actual output is specified in cospsimulator_intr_init.
        1013             :     do k=1,nprs_cosp
        1014             :        prstau_taumid_cosp(ntau_cosp*(k-1)+1:k*ntau_cosp)=taumid_cosp(1:ntau_cosp)
        1015             :        prstau_prsmid_cosp(ntau_cosp*(k-1)+1:k*ntau_cosp)=prsmid_cosp(k)
        1016             :        prstau_taumid_cosp_modis(ntau_cosp_modis*(k-1)+1:k*ntau_cosp_modis)=taumid_cosp_modis(1:ntau_cosp_modis)
        1017             :        prstau_prsmid_cosp_modis(ntau_cosp_modis*(k-1)+1:k*ntau_cosp_modis)=prsmid_cosp(k)
        1018             :     enddo
        1019             :
        1020             :     do k=1,nht_cosp
        1021             :        htdbze_dbzemid_cosp(CLOUDSAT_DBZE_BINS*(k-1)+1:k*CLOUDSAT_DBZE_BINS)=dbzemid_cosp(1:CLOUDSAT_DBZE_BINS)
        1022             :        htdbze_htmid_cosp(CLOUDSAT_DBZE_BINS*(k-1)+1:k*CLOUDSAT_DBZE_BINS)=htmid_cosp(k)
        1023             :     enddo
        1024             :
        1025             :     do k=1,nht_cosp
        1026             :        htsr_srmid_cosp(nsr_cosp*(k-1)+1:k*nsr_cosp)=srmid_cosp(1:nsr_cosp)
        1027             :        htsr_htmid_cosp(nsr_cosp*(k-1)+1:k*nsr_cosp)=htmid_cosp(k)
        1028             :     enddo
        1029             :
        1030             :     do k=1,nlay
        1031             :        htmlscol_scol_cosp(nscol_cosp*(k-1)+1:k*nscol_cosp)=scol_cosp(1:nscol_cosp)
        1032             :        htmlscol_htmlmid_cosp(nscol_cosp*(k-1)+1:k*nscol_cosp)=htmlmid_cosp(k)
        1033             :     enddo
        1034             :
        1035             :     do k=1,nhtmisr_cosp
        1036             :        htmisrtau_taumid_cosp(ntau_cosp*(k-1)+1:k*ntau_cosp)=taumid_cosp(1:ntau_cosp)
        1037             :        htmisrtau_htmisrmid_cosp(ntau_cosp*(k-1)+1:k*ntau_cosp)=htmisrmid_cosp(k)
        1038             :     enddo
        1039             :
        1040             :   end subroutine setcosp2values
        1041             : #endif
        1042             :
        1043             :   ! ######################################################################################
        1044             :   ! SUBROUTINE cospsimulator_intr_run
        1045             :   ! ######################################################################################
        1046           0 :   subroutine cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs,     &
        1047             :                                     cld_swtau_in, snow_tau_in, snow_emis_in)
        1048             :
        1049             :     use physics_types,        only: physics_state
        1050             :     use physics_buffer,       only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
        1051             :     use camsrfexch,           only: cam_in_t
        1052             :     use constituents,         only: cnst_get_ind
        1053             :     use rad_constituents,     only: rad_cnst_get_gas
        1054             :     use interpolate_data,     only: lininterp_init,lininterp,lininterp_finish,interp_type
        1055             :     use physconst,            only: pi, inverse_gravit => rga
        1056             :     use cam_history,          only: outfld,hist_fld_col_active
        1057             :     use cam_history_support,  only: max_fieldname_len
        1058             :
        1059             : #ifdef USE_COSP
        1060             :     use mod_cosp_config,      only: R_UNDEF,parasol_nrefl, Nlvgrid
        1061             :     use mod_cosp,             only: cosp_simulator
        1062             :     use mod_quickbeam_optics, only: size_distribution
        1063             : #endif
        1064             :
        1065             :     ! ######################################################################################
        1066             :     ! Inputs
        1067             :     ! ######################################################################################
        1068             :     type(physics_state), intent(in),target  :: state
        1069             :     type(physics_buffer_desc),      pointer :: pbuf(:)
        1070             :     type(cam_in_t),      intent(in)         :: cam_in
        1071             :     real(r8), intent(in) :: emis(pcols,pver)                  ! cloud longwave emissivity
        1072             :     real(r8), intent(in) :: coszrs(pcols)                     ! cosine solar zenith angle (to tell if day or night)
        1073             :     real(r8), intent(in),optional :: cld_swtau_in(pcols,pver) ! RRTM cld_swtau_in, read in using this variable
        1074             :     real(r8), intent(in),optional :: snow_tau_in(pcols,pver)  ! RRTM grid-box mean SW snow optical depth, used for CAM5 simulations
        1075             :     real(r8), intent(in),optional :: snow_emis_in(pcols,pver) ! RRTM grid-box mean LW snow optical depth, used for CAM5 simulations
        1076             :
        1077             : #ifdef USE_COSP
        1078             :     ! ######################################################################################
        1079             :     ! Local variables
        1080             :     ! ######################################################################################
        1081             :     integer :: lchnk                             ! chunk identifier
        1082             :     integer :: ncol                              ! number of active atmospheric columns
        1083             :     integer :: i, k, kk
        1084             :     integer :: itim_old
        1085             :     integer :: ip, it
        1086             :     integer :: ipt
        1087             :     integer :: ih, ihd, ihs, ihsc, ihm, ihmt, ihml
        1088             :     integer :: isc
        1089             :     integer :: is
        1090             :     integer :: id
        1091             :
        1092             :     real(r8), parameter :: rad2deg = 180._r8/pi
        1093             :
        1094             :     ! Microphysics variables
        1095             :     integer :: ixcldliq                                   ! cloud liquid amount index for state%q
        1096             :     integer :: ixcldice                                   ! cloud ice amount index
        1097             :
        1098             :     ! COSP-related local vars
        1099             :     type(cosp_outputs)        :: cospOUT                  ! COSP simulator outputs
        1100             :     type(cosp_optical_inputs) :: cospIN                   ! COSP optical (or derived?) fields needed by simulators
        1101             :     type(cosp_column_inputs)  :: cospstateIN              ! COSP model fields needed by simulators
        1102             :
        1103             :     ! COSP input variables that depend on CAM
        1104             :     integer :: Npoints                                    ! Number of gridpoints COSP will process
        1105             :     real(r8), parameter :: emsfc_lw = 0.99_r8             ! longwave emissivity of surface at 10.5 microns
        1106             :
        1107             :     ! Local vars related to calculations to go from CAM input to COSP input
        1108             :     ! cosp convective value includes both deep and shallow convection
        1109             :     real(r8), allocatable :: &
        1110             :        zmid(:,:), &           ! layer midpoint height asl (m)
        1111             :        zint(:,:), &           ! layer interface height asl (m)
        1112             :        surf_hgt(:), &         ! surface height (m)
        1113             :        landmask(:), &         ! landmask (0 or 1)
        1114             :        mr_ccliq(:,:), &       ! mixing_ratio_convective_cloud_liquid (kg/kg)
        1115             :        mr_ccice(:,:), &       ! mixing_ratio_convective_cloud_ice (kg/kg)
        1116             :        mr_lsliq(:,:), &       ! mixing_ratio_large_scale_cloud_liquid (kg/kg)
        1117             :        mr_lsice(:,:), &       ! mixing_ratio_large_scale_cloud_ice (kg/kg)
        1118             :        rain_cv(:,:), &        ! interface flux_convective_cloud_rain (kg m^-2 s^-1)
        1119             :        snow_cv(:,:), &        ! interface flux_convective_cloud_snow (kg m^-2 s^-1)
        1120             :        rain_cv_interp(:,:), & ! midpoint flux_convective_cloud_rain (kg m^-2 s^-1)
        1121             :        snow_cv_interp(:,:), & ! midpoint flux_convective_cloud_snow (kg m^-2 s^-1)
        1122             :        rain_ls_interp(:,:), & ! midpoint ls rain flux (kg m^-2 s^-1)
        1123             :        snow_ls_interp(:,:), & ! midpoint ls snow flux
        1124             :        grpl_ls_interp(:,:), & ! midpoint ls grp flux, set to 0
        1125             :        reff_cosp(:,:,:), &    ! effective radius for cosp input
        1126             :        dtau_s(:,:), &         ! Optical depth of stratiform cloud at 0.67 um
        1127             :        dtau_c(:,:), &         ! Optical depth of convective cloud at 0.67 um
        1128             :        dtau_s_snow(:,:), &    ! Grid-box mean Optical depth of stratiform snow at 0.67 um
        1129             :        dem_s(:,:), &          ! Longwave emis of stratiform cloud at 10.5 um
        1130             :        dem_c(:,:), &          ! Longwave emis of convective cloud at 10.5 um
        1131             :        dem_s_snow(:,:)        ! Grid-box mean Optical depth of stratiform snow at 10.5 um
        1132             :
        1133             :     integer :: cam_sunlit(pcols) ! cam_sunlit - Sunlit flag(1-sunlit/0-dark).
        1134             :     integer :: nSunLit           ! Number of sunlit (not sunlit) scenes.
        1135             :
        1136             :     ! ######################################################################################
        1137             :     ! Simulator output info
        1138             :     ! ######################################################################################
        1139             :     integer, parameter :: nf_radar=17                    ! number of radar outputs
        1140             :     integer, parameter :: nf_calipso=28                  ! number of calipso outputs
        1141             :     integer, parameter :: nf_isccp=9                     ! number of isccp outputs
        1142             :     integer, parameter :: nf_misr=1                      ! number of misr outputs
        1143             :     integer, parameter :: nf_modis=20                    ! number of modis outputs
        1144             :
        1145             :     ! Cloudsat outputs
        1146             :     character(len=max_fieldname_len),dimension(nf_radar),parameter ::          &
        1147             :          fname_radar = (/'CFAD_DBZE94_CS', 'CLD_CAL_NOTCS ', 'DBZE_CS       ', &
        1148             :                          'CLDTOT_CALCS  ', 'CLDTOT_CS     ', 'CLDTOT_CS2    ', &
        1149             :                          'CS_NOPRECIP   ', 'CS_RAINPOSS   ', 'CS_RAINPROB   ', &
        1150             :                          'CS_RAINCERT   ', 'CS_SNOWPOSS   ', 'CS_SNOWCERT   ', &
        1151             :                          'CS_MIXPOSS    ', 'CS_MIXCERT    ', 'CS_RAINHARD   ', &
        1152             :                          'CS_UN         ', 'CS_PIA        '/)
        1153             :
        1154             :     ! CALIPSO outputs
        1155             :     character(len=max_fieldname_len),dimension(nf_calipso),parameter :: &
        1156             :          fname_calipso=(/'CLDLOW_CAL     ','CLDMED_CAL     ','CLDHGH_CAL     ','CLDTOT_CAL     ','CLD_CAL        ',&
        1157             :                          'RFL_PARASOL    ','CFAD_SR532_CAL ','ATB532_CAL     ','MOL532_CAL     ','CLD_CAL_LIQ    ',&
        1158             :                          'CLD_CAL_ICE    ','CLD_CAL_UN     ','CLD_CAL_TMP    ','CLD_CAL_TMPLIQ ','CLD_CAL_TMPICE ',&
        1159             :                          'CLD_CAL_TMPUN  ','CLDTOT_CAL_ICE ','CLDTOT_CAL_LIQ ','CLDTOT_CAL_UN  ','CLDHGH_CAL_ICE ',&
        1160             :                          'CLDHGH_CAL_LIQ ','CLDHGH_CAL_UN  ','CLDMED_CAL_ICE ','CLDMED_CAL_LIQ ','CLDMED_CAL_UN  ',&
        1161             :                          'CLDLOW_CAL_ICE ','CLDLOW_CAL_LIQ ','CLDLOW_CAL_UN  '/)
        1162             :     ! ISCCP outputs
        1163             :     character(len=max_fieldname_len),dimension(nf_isccp),parameter :: &
        1164             :          fname_isccp=(/'FISCCP1_COSP    ','CLDTOT_ISCCP    ','MEANCLDALB_ISCCP',&
        1165             :                        'MEANPTOP_ISCCP  ','TAU_ISCCP       ','CLDPTOP_ISCCP   ','MEANTAU_ISCCP   ',&
        1166             :                        'MEANTB_ISCCP    ','MEANTBCLR_ISCCP '/)
        1167             :     ! MISR outputs
        1168             :     character(len=max_fieldname_len),dimension(nf_misr),parameter :: &
        1169             :          fname_misr=(/'CLD_MISR '/)
        1170             :     ! MODIS outputs
        1171             :     character(len=max_fieldname_len),dimension(nf_modis) :: &
        1172             :          fname_modis=(/'CLTMODIS    ','CLWMODIS    ','CLIMODIS    ','CLHMODIS    ','CLMMODIS    ',&
        1173             :                        'CLLMODIS    ','TAUTMODIS   ','TAUWMODIS   ','TAUIMODIS   ','TAUTLOGMODIS',&
        1174             :                        'TAUWLOGMODIS','TAUILOGMODIS','REFFCLWMODIS','REFFCLIMODIS',&
        1175             :                        'PCTMODIS    ','LWPMODIS    ','IWPMODIS    ','CLMODIS     ','CLRIMODIS   ',&
        1176             :                        'CLRLMODIS   '/)
        1177             :
        1178             :     logical :: run_radar(nf_radar,pcols)                 ! logical telling you if you should run radar simulator
        1179             :     logical :: run_calipso(nf_calipso,pcols)             ! logical telling you if you should run calipso simulator
        1180             :     logical :: run_isccp(nf_isccp,pcols)                 ! logical telling you if you should run isccp simulator
        1181             :     logical :: run_misr(nf_misr,pcols)                   ! logical telling you if you should run misr simulator
        1182             :     logical :: run_modis(nf_modis,pcols)                 ! logical telling you if you should run modis simulator
        1183             :
        1184             :     ! CAM pointers to get variables from radiation interface (get from rad_cnst_get_gas)
        1185             :     real(r8), pointer, dimension(:,:) :: q               ! specific humidity (kg/kg)
        1186             :     real(r8), pointer, dimension(:,:) :: o3              ! Mass mixing ratio 03
        1187             :
        1188             :     ! CAM pointers to get variables from the physics buffer
        1189             :     real(r8), pointer, dimension(:,:) :: cld             ! cloud fraction, tca - total_cloud_amount (0-1)
        1190             :     real(r8), pointer, dimension(:,:) :: concld          ! concld fraction, cca - convective_cloud_amount (0-1)
        1191             :     real(r8), pointer, dimension(:,:) :: rel             ! liquid effective drop radius (microns)
        1192             :     real(r8), pointer, dimension(:,:) :: rei             ! ice effective drop size (microns)
        1193             :     real(r8), pointer, dimension(:,:) :: ls_reffrain     ! rain effective drop radius (microns)
        1194             :     real(r8), pointer, dimension(:,:) :: ls_reffsnow     ! snow effective drop size (microns)
        1195             :     real(r8), pointer, dimension(:,:) :: cv_reffliq      ! convective cld liq effective drop radius (microns)
        1196             :     real(r8), pointer, dimension(:,:) :: cv_reffice      ! convective cld ice effective drop size (microns)
        1197             :
        1198             :     !! precip flux pointers
        1199             :     real(r8), target, dimension(pcols,pverp) :: zero_ifc ! zero array for interface fields not in the pbuf
        1200             :     real(r8), pointer, dimension(:,:) :: dp_flxprc       ! deep interface gbm flux_convective_cloud_rain+snow (kg m^-2 s^-1)
        1201             :     real(r8), pointer, dimension(:,:) :: dp_flxsnw       ! deep interface gbm flux_convective_cloud_snow (kg m^-2 s^-1)
        1202             :     real(r8), pointer, dimension(:,:) :: sh_flxprc       ! shallow interface gbm flux_convective_cloud_rain+snow (kg m^-2 s^-1)
        1203             :     real(r8), pointer, dimension(:,:) :: sh_flxsnw       ! shallow interface gbm flux_convective_cloud_snow (kg m^-2 s^-1)
        1204             :     real(r8), pointer, dimension(:,:) :: ls_flxprc       ! stratiform interface gbm flux_cloud_rain+snow (kg m^-2 s^-1)
        1205             :     real(r8), pointer, dimension(:,:) :: ls_flxsnw       ! stratiform interface gbm flux_cloud_snow (kg m^-2 s^-1)
        1206             :
        1207             :     !! grid box total cloud mixing ratio (large-scale + convective)
        1208             :     real(r8), pointer, dimension(:,:) :: totg_liq       ! gbm total cloud liquid water (kg/kg)
        1209             :     real(r8), pointer, dimension(:,:) :: totg_ice       ! gbm total cloud ice water (kg/kg)
        1210             :
        1211             :     ! Output CAM variables
        1212             :     ! Multiple "mdims" are collapsed because CAM history buffers only support one mdim.
        1213             :     ! MIXED DIMS: ntau_cosp*nprs_cosp, CLOUDSAT_DBZE_BINS*nht_cosp, nsr_cosp*nht_cosp, nscol_cosp*nlay,
        1214             :     !             ntau_cosp*nhtmisr_cosp
        1215             :     real(r8) :: clisccp2(pcols,ntau_cosp,nprs_cosp)
        1216             :     real(r8) :: cfad_dbze94(pcols,CLOUDSAT_DBZE_BINS,nht_cosp)
        1217             :     real(r8) :: cfad_lidarsr532(pcols,nsr_cosp,nht_cosp)
        1218             :     real(r8) :: dbze94(pcols,nscol_cosp,nlay)
        1219             :     real(r8) :: atb532(pcols,nscol_cosp,nlay)
        1220             :     real(r8) :: clMISR(pcols,ntau_cosp,nhtmisr_cosp)
        1221             :     real(r8) :: frac_out(pcols,nscol_cosp,nlay)
        1222             :     real(r8) :: cldtot_isccp(pcols)
        1223             :     real(r8) :: meancldalb_isccp(pcols)
        1224             :     real(r8) :: meanptop_isccp(pcols)
        1225             :     real(r8) :: cldlow_cal(pcols)
        1226             :     real(r8) :: cldmed_cal(pcols)
        1227             :     real(r8) :: cldhgh_cal(pcols)
        1228             :     real(r8) :: cldtot_cal(pcols)
        1229             :     real(r8) :: cldtot_cal_ice(pcols)
        1230             :     real(r8) :: cldtot_cal_liq(pcols)
        1231             :     real(r8) :: cldtot_cal_un(pcols)
        1232             :     real(r8) :: cldhgh_cal_ice(pcols)
        1233             :     real(r8) :: cldhgh_cal_liq(pcols)
        1234             :     real(r8) :: cldhgh_cal_un(pcols)
        1235             :     real(r8) :: cldmed_cal_ice(pcols)
        1236             :     real(r8) :: cldmed_cal_liq(pcols)
        1237             :     real(r8) :: cldmed_cal_un(pcols)
        1238             :     real(r8) :: cldlow_cal_ice(pcols)
        1239             :     real(r8) :: cldlow_cal_liq(pcols)
        1240             :     real(r8) :: cldlow_cal_un(pcols)
        1241             :     real(r8) :: cld_cal(pcols,nht_cosp)
        1242             :     real(r8) :: cld_cal_liq(pcols,nht_cosp)
        1243             :     real(r8) :: cld_cal_ice(pcols,nht_cosp)
        1244             :     real(r8) :: cld_cal_un(pcols,nht_cosp)
        1245             :     real(r8) :: cld_cal_tmp(pcols,nht_cosp)
        1246             :     real(r8) :: cld_cal_tmpliq(pcols,nht_cosp)
        1247             :     real(r8) :: cld_cal_tmpice(pcols,nht_cosp)
        1248             :     real(r8) :: cld_cal_tmpun(pcols,nht_cosp)
        1249             :     real(r8) :: cfad_dbze94_cs(pcols,nht_cosp*CLOUDSAT_DBZE_BINS)
        1250             :     real(r8) :: cfad_sr532_cal(pcols,nht_cosp*nsr_cosp)
        1251             :     real(r8) :: tau_isccp(pcols,nscol_cosp)
        1252             :     real(r8) :: cldptop_isccp(pcols,nscol_cosp)
        1253             :     real(r8) :: meantau_isccp(pcols)
        1254             :     real(r8) :: meantb_isccp(pcols)
        1255             :     real(r8) :: meantbclr_isccp(pcols)
        1256             :     real(r8) :: dbze_cs(pcols,nlay*nscol_cosp)
        1257             :     real(r8) :: cldtot_calcs(pcols)
        1258             :     real(r8) :: cldtot_cs(pcols)
        1259             :     real(r8) :: cldtot_cs2(pcols)
        1260             :     real(r8) :: ptcloudsatflag0(pcols)
        1261             :     real(r8) :: ptcloudsatflag1(pcols)
        1262             :     real(r8) :: ptcloudsatflag2(pcols)
        1263             :     real(r8) :: ptcloudsatflag3(pcols)
        1264             :     real(r8) :: ptcloudsatflag4(pcols)
        1265             :     real(r8) :: ptcloudsatflag5(pcols)
        1266             :     real(r8) :: ptcloudsatflag6(pcols)
        1267             :     real(r8) :: ptcloudsatflag7(pcols)
        1268             :     real(r8) :: ptcloudsatflag8(pcols)
        1269             :     real(r8) :: ptcloudsatflag9(pcols)
        1270             :     real(r8) :: cloudsatpia(pcols)
        1271             :     real(r8) :: cld_cal_notcs(pcols,nht_cosp)
        1272             :     real(r8) :: atb532_cal(pcols,nlay*nscol_cosp)
        1273             :     real(r8) :: mol532_cal(pcols,nlay)
        1274             :     real(r8) :: cld_misr(pcols,nhtmisr_cosp*ntau_cosp)
        1275             :     real(r8) :: refl_parasol(pcols,nsza_cosp)
        1276             :     real(r8) :: scops_out(pcols,nlay*nscol_cosp)
        1277             :     real(r8) :: cltmodis(pcols)
        1278             :     real(r8) :: clwmodis(pcols)
        1279             :     real(r8) :: climodis(pcols)
        1280             :     real(r8) :: clhmodis(pcols)
        1281             :     real(r8) :: clmmodis(pcols)
        1282             :     real(r8) :: cllmodis(pcols)
        1283             :     real(r8) :: tautmodis(pcols)
        1284             :     real(r8) :: tauwmodis(pcols)
        1285             :     real(r8) :: tauimodis(pcols)
        1286             :     real(r8) :: tautlogmodis(pcols)
        1287             :     real(r8) :: tauwlogmodis(pcols)
        1288             :     real(r8) :: tauilogmodis(pcols)
        1289             :     real(r8) :: reffclwmodis(pcols)
        1290             :     real(r8) :: reffclimodis(pcols)
        1291             :     real(r8) :: pctmodis(pcols)
        1292             :     real(r8) :: lwpmodis(pcols)
        1293             :     real(r8) :: iwpmodis(pcols)
        1294             :     real(r8) :: clmodis_cam(pcols,ntau_cosp_modis*nprs_cosp)
        1295             :     real(r8) :: clmodis(pcols,ntau_cosp_modis,nprs_cosp)
        1296             :     real(r8) :: clrimodis_cam(pcols,ntau_cosp*numMODISReffIceBins)
        1297             :     real(r8) :: clrimodis(pcols,ntau_cosp,numMODISReffIceBins)
        1298             :     real(r8) :: clrlmodis_cam(pcols,ntau_cosp*numMODISReffLiqBins)
        1299             :     real(r8) :: clrlmodis(pcols,ntau_cosp,numMODISReffLiqBins)
        1300             :     real(r8), dimension(pcols,nlay*nscol_cosp) :: &
        1301             :        tau067_out, emis11_out, fracliq_out, asym34_out, ssa34_out
        1302             :
        1303             :     type(interp_type)  :: interp_wgts
        1304             :     integer, parameter :: extrap_method = 1   ! sets extrapolation method to boundary value (1)
        1305             :
        1306             :     ! COSPv2 stuff
        1307             :     character(len=256),dimension(100) :: cosp_status
        1308             :     integer :: nerror
        1309             :
        1310             :     integer :: istat
        1311             :     character(len=*), parameter :: sub = 'cospsimulator_intr_run'
        1312             :     !--------------------------------------------------------------------------------------
        1313             :
        1314             :     call t_startf("init_and_stuff")
        1315             :     ! ######################################################################################
        1316             :     ! Initialization
        1317             :     ! ######################################################################################
        1318             :
        1319             :     lchnk = state%lchnk    ! chunk ID
        1320             :     ncol  = state%ncol     ! number of columns in the chunk
        1321             :     Npoints = ncol         ! number of COSP gridpoints
        1322             :
        1323             :     zero_ifc = 0._r8
        1324             :
        1325             :     ! Initialize CAM variables as R_UNDEF, important for history files because it will exclude these from averages
        1326             :     ! initialize over all pcols, not just ncol.  missing values needed in chunks where ncol<pcols
        1327             :     clisccp2(1:pcols,1:ntau_cosp,1:nprs_cosp)     = R_UNDEF
        1328             :     cfad_dbze94(1:pcols,1:CLOUDSAT_DBZE_BINS,1:nht_cosp)  = R_UNDEF
        1329             :     cfad_lidarsr532(1:pcols,1:nsr_cosp,1:nht_cosp)= R_UNDEF
        1330             :     dbze94(1:pcols,1:nscol_cosp,1:nlay)           = R_UNDEF
        1331             :     atb532(1:pcols,1:nscol_cosp,1:nlay)           = R_UNDEF
        1332             :     clMISR(1:pcols,ntau_cosp,1:nhtmisr_cosp)      = R_UNDEF
        1333             :     frac_out(1:pcols,1:nscol_cosp,1:nlay)         = R_UNDEF
        1334             :
        1335             :     ! (all CAM output variables. including collapsed variables)
        1336             :     cldtot_isccp(1:pcols)                            = R_UNDEF
        1337             :     meancldalb_isccp(1:pcols)                        = R_UNDEF
        1338             :     meanptop_isccp(1:pcols)                          = R_UNDEF
        1339             :     cldlow_cal(1:pcols)                              = R_UNDEF
        1340             :     cldmed_cal(1:pcols)                              = R_UNDEF
        1341             :     cldhgh_cal(1:pcols)                              = R_UNDEF
        1342             :     cldtot_cal(1:pcols)                              = R_UNDEF
        1343             :     cldtot_cal_ice(1:pcols)                          = R_UNDEF !+cosp1.4
        1344             :     cldtot_cal_liq(1:pcols)                          = R_UNDEF
        1345             :     cldtot_cal_un(1:pcols)                           = R_UNDEF
        1346             :     cldhgh_cal_ice(1:pcols)                          = R_UNDEF
        1347             :     cldhgh_cal_liq(1:pcols)                          = R_UNDEF
        1348             :     cldhgh_cal_un(1:pcols)                           = R_UNDEF
        1349             :     cldmed_cal_ice(1:pcols)                          = R_UNDEF
        1350             :     cldmed_cal_liq(1:pcols)                          = R_UNDEF
        1351             :     cldmed_cal_un(1:pcols)                           = R_UNDEF
        1352             :     cldlow_cal_liq(1:pcols)                          = R_UNDEF
        1353             :     cldlow_cal_ice(1:pcols)                          = R_UNDEF
        1354             :     cldlow_cal_un(1:pcols)                           = R_UNDEF !+cosp1.4
        1355             :     cld_cal(1:pcols,1:nht_cosp)                      = R_UNDEF
        1356             :     cld_cal_liq(1:pcols,1:nht_cosp)                  = R_UNDEF !+cosp1.4
        1357             :     cld_cal_ice(1:pcols,1:nht_cosp)                  = R_UNDEF
        1358             :     cld_cal_un(1:pcols,1:nht_cosp)                   = R_UNDEF
        1359             :     cld_cal_tmp(1:pcols,1:nht_cosp)                  = R_UNDEF
        1360             :     cld_cal_tmpliq(1:pcols,1:nht_cosp)               = R_UNDEF
        1361             :     cld_cal_tmpice(1:pcols,1:nht_cosp)               = R_UNDEF
        1362             :     cld_cal_tmpun(1:pcols,1:nht_cosp)                = R_UNDEF
        1363             :     cfad_dbze94_cs(1:pcols,1:nht_cosp*CLOUDSAT_DBZE_BINS)    = R_UNDEF
        1364             :     cfad_sr532_cal(1:pcols,1:nht_cosp*nsr_cosp)      = R_UNDEF
        1365             :     tau_isccp(1:pcols,1:nscol_cosp)                  = R_UNDEF
        1366             :     cldptop_isccp(1:pcols,1:nscol_cosp)              = R_UNDEF
        1367             :     meantau_isccp(1:pcols)                           = R_UNDEF
        1368             :     meantb_isccp(1:pcols)                            = R_UNDEF
        1369             :     meantbclr_isccp(1:pcols)                         = R_UNDEF
        1370             :     dbze_cs(1:pcols,1:nlay*nscol_cosp)               = R_UNDEF
        1371             :     ptcloudsatflag0(1:pcols)                         = R_UNDEF
        1372             :     ptcloudsatflag1(1:pcols)                         = R_UNDEF
        1373             :     ptcloudsatflag2(1:pcols)                         = R_UNDEF
        1374             :     ptcloudsatflag3(1:pcols)                         = R_UNDEF
        1375             :     ptcloudsatflag4(1:pcols)                         = R_UNDEF
        1376             :     ptcloudsatflag5(1:pcols)                         = R_UNDEF
        1377             :     ptcloudsatflag6(1:pcols)                         = R_UNDEF
        1378             :     ptcloudsatflag7(1:pcols)                         = R_UNDEF
        1379             :     ptcloudsatflag8(1:pcols)                         = R_UNDEF
        1380             :     ptcloudsatflag9(1:pcols)                         = R_UNDEF
        1381             :     cloudsatpia(1:pcols)                             = R_UNDEF
        1382             :     cldtot_calcs(1:pcols)                            = R_UNDEF
        1383             :     cldtot_cs(1:pcols)                               = R_UNDEF
        1384             :     cldtot_cs2(1:pcols)                              = R_UNDEF
        1385             :     cld_cal_notcs(1:pcols,1:nht_cosp)                = R_UNDEF
        1386             :     atb532_cal(1:pcols,1:nlay*nscol_cosp)            = R_UNDEF
        1387             :     mol532_cal(1:pcols,1:nlay)                       = R_UNDEF
        1388             :     cld_misr(1:pcols,1:nhtmisr_cosp*ntau_cosp)       = R_UNDEF
        1389             :     refl_parasol(1:pcols,1:nsza_cosp)                = R_UNDEF
        1390             :     scops_out(1:pcols,1:nlay*nscol_cosp)             = R_UNDEF
        1391             :     cltmodis(1:pcols)                                = R_UNDEF
        1392             :     clwmodis(1:pcols)                                = R_UNDEF
        1393             :     climodis(1:pcols)                                = R_UNDEF
        1394             :     clhmodis(1:pcols)                                = R_UNDEF
        1395             :     clmmodis(1:pcols)                                = R_UNDEF
        1396             :     cllmodis(1:pcols)                                = R_UNDEF
        1397             :     tautmodis(1:pcols)                               = R_UNDEF
        1398             :     tauwmodis(1:pcols)                               = R_UNDEF
        1399             :     tauimodis(1:pcols)                               = R_UNDEF
        1400             :     tautlogmodis(1:pcols)                            = R_UNDEF
        1401             :     tauwlogmodis(1:pcols)                            = R_UNDEF
        1402             :     tauilogmodis(1:pcols)                            = R_UNDEF
        1403             :     reffclwmodis(1:pcols)                            = R_UNDEF
        1404             :     reffclimodis(1:pcols)                            = R_UNDEF
        1405             :     pctmodis(1:pcols)                                = R_UNDEF
        1406             :     lwpmodis(1:pcols)                                = R_UNDEF
        1407             :     iwpmodis(1:pcols)                                = R_UNDEF
        1408             :     clmodis_cam(1:pcols,1:ntau_cosp_modis*nprs_cosp) = R_UNDEF
        1409             :     clmodis(1:pcols,1:ntau_cosp_modis,1:nprs_cosp)   = R_UNDEF
        1410             :     clrimodis_cam(1:pcols,1:ntau_cosp_modis*numMODISReffIceBins) = R_UNDEF ! +cosp2
        1411             :     clrimodis(1:pcols,1:ntau_cosp_modis,1:numMODISReffIceBins)   = R_UNDEF ! +cosp2
        1412             :     clrlmodis_cam(1:pcols,1:ntau_cosp_modis*numMODISReffLiqBins) = R_UNDEF ! +cosp2
        1413             :     clrlmodis(1:pcols,1:ntau_cosp_modis,1:numMODISReffLiqBins)   = R_UNDEF ! +cosp2
        1414             :     tau067_out(1:pcols,1:nlay*nscol_cosp)            = R_UNDEF ! +cosp2
        1415             :     emis11_out(1:pcols,1:nlay*nscol_cosp)            = R_UNDEF ! +cosp2
        1416             :     asym34_out(1:pcols,1:nlay*nscol_cosp)            = R_UNDEF ! +cosp2
        1417             :     ssa34_out(1:pcols,1:nlay*nscol_cosp)             = R_UNDEF ! +cosp2
        1418             :     fracLiq_out(1:pcols,1:nlay*nscol_cosp)           = R_UNDEF ! +cosp2
        1419             :
        1420             :     ! ######################################################################################
        1421             :     ! DECIDE WHICH COLUMNS YOU ARE GOING TO RUN COSP ON....
        1422             :     ! ######################################################################################
        1423             :
        1424             :     !! run_cosp is set for each column in each chunk in the first timestep of the run
        1425             :     !! hist_fld_col_active in cam_history.F90 is used to decide if you need to run cosp.
        1426             :     if (first_run_cosp(lchnk)) then
        1427             :        !! initalize to run logicals as false
        1428             :        run_cosp(1:ncol,lchnk)=.false.
        1429             :        run_radar(1:nf_radar,1:ncol)=.false.
        1430             :        run_calipso(1:nf_calipso,1:ncol)=.false.
        1431             :        run_isccp(1:nf_isccp,1:ncol)=.false.
        1432             :        run_misr(1:nf_misr,1:ncol)=.false.
        1433             :        run_modis(1:nf_modis,1:ncol)=.false.
        1434             :
        1435             :        if (lradar_sim) then
        1436             :           do i=1,nf_radar
        1437             :              run_radar(i,1:pcols)=hist_fld_col_active(fname_radar(i),lchnk,pcols)
        1438             :           end do
        1439             :        end if
        1440             :        if (llidar_sim) then
        1441             :           do i=1,nf_calipso
        1442             :              run_calipso(i,1:pcols)=hist_fld_col_active(fname_calipso(i),lchnk,pcols)
        1443             :           end do
        1444             :        end if
        1445             :        if (lisccp_sim) then
        1446             :           do i=1,nf_isccp
        1447             :              run_isccp(i,1:pcols)=hist_fld_col_active(fname_isccp(i),lchnk,pcols)
        1448             :           end do
        1449             :        end if
        1450             :        if (lmisr_sim) then
        1451             :           do i=1,nf_misr
        1452             :              run_misr(i,1:pcols)=hist_fld_col_active(fname_misr(i),lchnk,pcols)
        1453             :           end do
        1454             :        end if
        1455             :        if (lmodis_sim) then
        1456             :           do i=1,nf_modis
        1457             :              run_modis(i,1:pcols)=hist_fld_col_active(fname_modis(i),lchnk,pcols)
        1458             :           end do
        1459             :        end if
        1460             :
        1461             :        do i=1,ncol
        1462             :           if ((any(run_radar(:,i))) .or. (any(run_calipso(:,i))) .or. (any(run_isccp(:,i))) &
        1463             :                .or. (any(run_misr(:,i))) .or. (any(run_modis(:,i)))) then
        1464             :              run_cosp(i,lchnk)=.true.
        1465             :           end if
        1466             :        end do
        1467             :
        1468             :        first_run_cosp(lchnk)=.false.
        1469             :     endif
        1470             :
        1471             :     ! ######################################################################################
        1472             :     ! GET CAM GEOPHYSICAL VARIABLES NEEDED FOR COSP INPUT
        1473             :     ! ######################################################################################
        1474             :     ! state variables (prognostic variables, see physics_types.F90)
        1475             :     ! state%lat   ! lat (radians)
        1476             :     ! state%lon   ! lon (radians)
        1477             :     ! state%t     ! temperature (K)
        1478             :     ! state%ps    ! surface pressure (Pa)
        1479             :     ! state%pint  ! p - p_in_full_levels (Pa)
        1480             :     ! state%pmid  ! ph - p_in_half_levels (Pa)
        1481             :     ! state%zm    ! geopotential height above surface at midpoints (m), pver
        1482             :     ! state%zi    ! geopotential height above surface at interfaces (m), pverp
        1483             :     ! state%phis  ! surface geopotential (m2/s2)
        1484             :     ! NOTE: The state variables state%q(:,:,ixcldliq)/state%q(:,:,ixcldice) are grid-box
        1485             :     ! quantities for the stratiform clouds only.  stratiform water * stratiform cloud fraction
        1486             :     ! state%q(:,:,ixcldliq) ! for CAM4: cldliq = stratiform incld water content * total cloud fraction
        1487             :     ! state%q(:,:,ixcldice) ! for CAM4: cldice = stratiform incld ice content * total cloud fraction
        1488             :
        1489             :     ! advected constiutent indices
        1490             :     call cnst_get_ind('CLDLIQ', ixcldliq)
        1491             :     call cnst_get_ind('CLDICE', ixcldice)
        1492             :
        1493             :     ! radiative constituents (prognostic or data)
        1494             :     call rad_cnst_get_gas(0,'H2O', state, pbuf,  q)
        1495             :     call rad_cnst_get_gas(0,'O3',  state, pbuf,  o3)
        1496             :
        1497             :     ! fields from physics buffer
        1498             :     itim_old = pbuf_old_tim_idx()
        1499             :     call pbuf_get_field(pbuf, cld_idx,    cld,    start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
        1500             :     call pbuf_get_field(pbuf, concld_idx, concld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
        1501             :     call pbuf_get_field(pbuf, rel_idx, rel  )
        1502             :     call pbuf_get_field(pbuf, rei_idx, rei)
        1503             :     call pbuf_get_field(pbuf, lsreffrain_idx, ls_reffrain  )
        1504             :     call pbuf_get_field(pbuf, lsreffsnow_idx, ls_reffsnow  )
        1505             :     call pbuf_get_field(pbuf, cvreffliq_idx,  cv_reffliq   )
        1506             :     call pbuf_get_field(pbuf, cvreffice_idx,  cv_reffice   )
        1507             :
        1508             :     !! grid box total cloud mixing ratios
        1509             :     call pbuf_get_field(pbuf, gb_totcldliqmr_idx, totg_liq)
        1510             :     call pbuf_get_field(pbuf, gb_totcldicemr_idx, totg_ice)
        1511             :
        1512             :     !! precipitation fluxes
        1513             :     call pbuf_get_field(pbuf, dpflxprc_idx, dp_flxprc  )
        1514             :     call pbuf_get_field(pbuf, dpflxsnw_idx, dp_flxsnw  )
        1515             :     if (shflxprc_idx > 0) then
        1516             :        call pbuf_get_field(pbuf, shflxprc_idx, sh_flxprc  )
        1517             :     else
        1518             :        sh_flxprc => zero_ifc
        1519             :     end if
        1520             :     if (shflxsnw_idx > 0) then
        1521             :        call pbuf_get_field(pbuf, shflxsnw_idx, sh_flxsnw  )
        1522             :     else
        1523             :        sh_flxsnw => zero_ifc
        1524             :     end if
        1525             :     call pbuf_get_field(pbuf, lsflxprc_idx, ls_flxprc  )
        1526             :     call pbuf_get_field(pbuf, lsflxsnw_idx, ls_flxsnw  )
        1527             :
        1528             :     !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        1529             :     ! CALCULATE COSP INPUT VARIABLES FROM CAM VARIABLES, done for all columns within chunk
        1530             :     !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        1531             :
        1532             :     ! These arrays are dimensioned to only include active columns (ncol), and the number
        1533             :     ! of layers (nlay) and layer interfaces (nlayp) operated on by COSP.
        1534             :     allocate( &
        1535             :        zmid(ncol,nlay), &
        1536             :        zint(ncol,nlayp), &
        1537             :        surf_hgt(ncol),   &
        1538             :        landmask(ncol),   &
        1539             :        mr_ccliq(ncol,nlay), &
        1540             :        mr_ccice(ncol,nlay), &
        1541             :        mr_lsliq(ncol,nlay), &
        1542             :        mr_lsice(ncol,nlay), &
        1543             :        rain_cv(ncol,nlayp), &
        1544             :        snow_cv(ncol,nlayp), &
        1545             :        rain_cv_interp(ncol,nlay), &
        1546             :        snow_cv_interp(ncol,nlay), &
        1547             :        rain_ls_interp(ncol,nlay), &
        1548             :        snow_ls_interp(ncol,nlay), &
        1549             :        grpl_ls_interp(ncol,nlay), &
        1550             :        reff_cosp(ncol,nlay,nhydro), &
        1551             :        dtau_s(ncol,nlay), &
        1552             :        dtau_c(ncol,nlay), &
        1553             :        dtau_s_snow(ncol,nlay), &
        1554             :        dem_s(ncol,nlay), &
        1555             :        dem_c(ncol,nlay), &
        1556             :        dem_s_snow(ncol,nlay), stat=istat)
        1557             :     call handle_allocate_error(istat, sub, 'zmid,..,dem_s_snow')
        1558             :
        1559             :     ! add surface height (surface geopotential/gravity) to convert CAM heights based on
        1560             :     ! geopotential above surface into height above sea level
        1561             :     surf_hgt = state%phis(:ncol)*inverse_gravit
        1562             :     do k = 1, nlay
        1563             :        zmid(:,k) = state%zm(:ncol,ktop+k-1) + surf_hgt
        1564             :        zint(:,k) = state%zi(:ncol,ktop+k-1) + surf_hgt
        1565             :     end do
        1566             :     zint(:,nlayp) = surf_hgt
        1567             :
        1568             :     landmask = 0._r8
        1569             :     do i = 1, ncol
        1570             :        if (cam_in%landfrac(i) > 0.01_r8) landmask(i)= 1
        1571             :     end do
        1572             :
        1573             :     ! Add together deep and shallow convection precipitation fluxes.
        1574             :     ! Note: sh_flxprc and dp_flxprc variables are rain+snow
        1575             :     rain_cv = (sh_flxprc(:ncol,ktop:pverp) - sh_flxsnw(:ncol,ktop:pverp)) + &
        1576             :               (dp_flxprc(:ncol,ktop:pverp) - dp_flxsnw(:ncol,ktop:pverp))
        1577             :     snow_cv = sh_flxsnw(:ncol,ktop:pverp) + dp_flxsnw(:ncol,ktop:pverp)
        1578             :
        1579             :     ! interpolate interface precip fluxes to mid points
        1580             :     do i = 1, ncol
        1581             :        ! find weights
        1582             :        call lininterp_init(state%zi(i,ktop:pverp), nlayp, state%zm(i,ktop:pver), nlay, &
        1583             :                            extrap_method, interp_wgts)
        1584             :        ! interpolate  lininterp(arrin, nin, arrout, nout, interp_wgts)
        1585             :        call lininterp(rain_cv(i,:), nlayp, rain_cv_interp(i,:), nlay, interp_wgts)
        1586             :        call lininterp(snow_cv(i,:), nlayp, snow_cv_interp(i,:), nlay, interp_wgts)
        1587             :        call lininterp(ls_flxprc(i,ktop:pverp), nlayp, rain_ls_interp(i,:), nlay, interp_wgts)
        1588             :        call lininterp(ls_flxsnw(i,ktop:pverp), nlayp, snow_ls_interp(i,:), nlay, interp_wgts)
        1589             :        call lininterp_finish(interp_wgts)
        1590             :        !! ls_flxprc is for rain+snow, find rain_ls_interp by subtracting off snow_ls_interp
        1591             :        rain_ls_interp(i,:) = rain_ls_interp(i,:) - snow_ls_interp(i,:)
        1592             :     end do
        1593             :
        1594             :     !! Make sure interpolated values are not less than 0
        1595             :     do k = 1, nlay
        1596             :        do i = 1, ncol
        1597             :           if (rain_ls_interp(i,k) < 0._r8) then
        1598             :              rain_ls_interp(i,k) = 0._r8
        1599             :           end if
        1600             :           if (snow_ls_interp(i,k) < 0._r8) then
        1601             :              snow_ls_interp(i,k) = 0._r8
        1602             :           end if
        1603             :           if (rain_cv_interp(i,k) < 0._r8) then
        1604             :              rain_cv_interp(i,k) = 0._r8
        1605             :           end if
        1606             :           if (snow_cv_interp(i,k) < 0._r8) then
        1607             :              snow_cv_interp(i,k) = 0._r8
        1608             :           end if
        1609             :        end do
        1610             :     end do
        1611             :
        1612             :     grpl_ls_interp = 0._r8
        1613             :
        1614             :     ! subroutine subsample_and_optics provides separate arguments to pass
        1615             :     ! the large scale and convective cloud condensate.  Below the grid box
        1616             :     ! total cloud water mixing ratios are passed in the arrays for the
        1617             :     ! large scale contributions and the arrays for the convective
        1618             :     ! contributions are set to zero.  This is consistent with the treatment
        1619             :     ! of cloud water by the radiation code.
        1620             :     mr_ccliq = 0._r8
        1621             :     mr_ccice = 0._r8
        1622             :     mr_lsliq = 0._r8
        1623             :     mr_lsice = 0._r8
        1624             :     do k = 1, nlay
        1625             :        kk = ktop + k -1
        1626             :        do i = 1, ncol
        1627             :           if (cld(i,k) > 0._r8) then
        1628             :              mr_lsliq(i,k) = totg_liq(i,kk)
        1629             :              mr_lsice(i,k) = totg_ice(i,kk)
        1630             :           end if
        1631             :        end do
        1632             :     end do
        1633             :
        1634             :     !! The specification of reff_cosp now follows e-mail discussion with Yuying in January 2011.
        1635             :     !! The values from the physics buffer are in microns... convert to meters for COSP.
        1636             :     reff_cosp(:,:,I_LSCLIQ) = rel(:ncol,ktop:pver)*1.e-6_r8
        1637             :     reff_cosp(:,:,I_LSCICE) = rei(:ncol,ktop:pver)*1.e-6_r8
        1638             :     reff_cosp(:,:,I_LSRAIN) = ls_reffrain(:ncol,ktop:pver)*1.e-6_r8
        1639             :     reff_cosp(:,:,I_LSSNOW) = ls_reffsnow(:ncol,ktop:pver)*1.e-6_r8
        1640             :     reff_cosp(:,:,I_CVCLIQ) = cv_reffliq(:ncol,ktop:pver)*1.e-6_r8
        1641             :     reff_cosp(:,:,I_CVCICE) = cv_reffice(:ncol,ktop:pver)*1.e-6_r8
        1642             :     reff_cosp(:,:,I_CVRAIN) = ls_reffrain(:ncol,ktop:pver)*1.e-6_r8  !! same as stratiform per Andrew
        1643             :     reff_cosp(:,:,I_CVSNOW) = ls_reffsnow(:ncol,ktop:pver)*1.e-6_r8  !! same as stratiform per Andrew
        1644             :     reff_cosp(:,:,I_LSGRPL) = 0._r8                                  !! using radar default reff
        1645             :
        1646             :     ! assign optical depths and emissivities
        1647             :     ! CAM4 assumes same radiative properties for stratiform and convective clouds,
        1648             :     !   (see ISCCP_CLOUD_TYPES subroutine call in cloudsimulator.F90)
        1649             :     !   Assume CAM5 is doing the same thing based on the ISCCP simulator calls within RRTM's radiation.F90
        1650             :     ! COSP wants in-cloud values.  CAM5 values cld_swtau are in-cloud.
        1651             :     ! snow_tau_in and snow_emis_in are passed without modification to COSP
        1652             :     dtau_s      = cld_swtau_in(:ncol,ktop:pver)
        1653             :     dtau_c      = cld_swtau_in(:ncol,ktop:pver)
        1654             :     dtau_s_snow = snow_tau_in(:ncol,ktop:pver)
        1655             :     dem_s       = emis(:ncol,ktop:pver)
        1656             :     dem_c       = emis(:ncol,ktop:pver)
        1657             :     dem_s_snow  = snow_emis_in(:ncol,ktop:pver)
        1658             :
        1659             :     ! ######################################################################################
        1660             :     ! Compute sunlit flag. If cosp_runall=.true., then run on all points.
        1661             :     ! ######################################################################################
        1662             :     cam_sunlit(:) = 0
        1663             :     if (cosp_runall) then
        1664             :        cam_sunlit(:) = 1
        1665             :        nSunLit   = ncol
        1666             :     else
        1667             :        nSunLit   = 0
        1668             :        do i=1,ncol
        1669             :           if ((coszrs(i) > 0.0_r8) .and. (run_cosp(i,lchnk))) then
        1670             :              cam_sunlit(i) = 1
        1671             :              nSunLit   = nSunLit+1
        1672             :           endif
        1673             :        enddo
        1674             :     endif
        1675             :     call t_stopf("init_and_stuff")
        1676             :
        1677             :     ! ######################################################################################
        1678             :     ! Construct COSP output derived type.
        1679             :     ! ######################################################################################
        1680             :     call t_startf("construct_cosp_outputs")
        1681             :     call construct_cosp_outputs(ncol, nscol_cosp, nlay, Nlvgrid, cospOUT)
        1682             :     call t_stopf("construct_cosp_outputs")
        1683             :
        1684             :     ! ######################################################################################
        1685             :     ! Construct and populate COSP input types
        1686             :     ! ######################################################################################
        1687             :     ! Model state
        1688             :     call t_startf("construct_cospstateIN")
        1689             :
        1690             :     call construct_cospstateIN(ncol, nlay, 0, cospstateIN)
        1691             :
        1692             :     ! convert to degrees.  Lat in range [-90,..,90], Lon in range [0,..,360]
        1693             :     cospstateIN%lat             = state%lat(:ncol)*rad2deg
        1694             :     cospstateIN%lon             = state%lon(:ncol)*rad2deg
        1695             :     cospstateIN%at              = state%t(:ncol,ktop:pver)
        1696             :     cospstateIN%qv              = q(:ncol,ktop:pver)
        1697             :     cospstateIN%o3              = o3(:ncol,ktop:pver)
        1698             :     cospstateIN%sunlit          = cam_sunlit(:ncol)
        1699             :     cospstateIN%skt             = cam_in%ts(:ncol)
        1700             :     cospstateIN%land            = landmask
        1701             :     cospstateIN%pfull           = state%pmid(:ncol,ktop:pver)
        1702             :     cospstateIN%phalf           = state%pint(:ncol,ktop:pverp)
        1703             :     cospstateIN%hgt_matrix      = zmid
        1704             :     cospstateIN%hgt_matrix_half = zint
        1705             :     cospstateIN%surfelev        = surf_hgt
        1706             :     call t_stopf("construct_cospstateIN")
        1707             :
        1708             :     ! Optical inputs
        1709             :     call t_startf("construct_cospIN")
        1710             :     call construct_cospIN(ncol, nscol_cosp, nlay, cospIN)
        1711             :     cospIN%emsfc_lw = emsfc_lw
        1712             :     if (lradar_sim) cospIN%rcfg_cloudsat = rcfg_cs(lchnk)
        1713             :     call t_stopf("construct_cospIN")
        1714             :
        1715             :     call t_startf("subsample_and_optics")
        1716             :     ! The arrays passed here contain only active columns and the limited vertical
        1717             :     ! domain operated on by COSP.  Unsubscripted array arguments have already been
        1718             :     ! allocated to the correct size.  Arrays the size of a CAM chunk (pcol,pver)
        1719             :     ! need to pass the correct section (:ncol,ktop:pver).
        1720             :     call subsample_and_optics( &
        1721             :        ncol, nlay, nscol_cosp, nhydro, overlap, &
        1722             :        lidar_ice_type, sd_cs(lchnk), &
        1723             :        cld(:ncol,ktop:pver), concld(:ncol,ktop:pver), &
        1724             :        rain_ls_interp, snow_ls_interp, grpl_ls_interp, rain_cv_interp, &
        1725             :        snow_cv_interp, mr_lsliq, mr_lsice, mr_ccliq, mr_ccice, &
        1726             :        reff_cosp, dtau_c, dtau_s ,dem_c, dem_s, dtau_s_snow, &
        1727             :        dem_s_snow, state%ps(:ncol), cospstateIN, cospIN)
        1728             :     call t_stopf("subsample_and_optics")
        1729             :
        1730             :     ! ######################################################################################
        1731             :     ! Call COSP
        1732             :     ! ######################################################################################
        1733             :     call t_startf("cosp_simulator")
        1734             :     cosp_status = COSP_SIMULATOR(cospIN, cospstateIN, cospOUT, start_idx=1, stop_idx=ncol,debug=.false.)
        1735             :
        1736             :     ! Check status flags
        1737             :     nerror = 0
        1738             :     do i = 1, ubound(cosp_status, 1)
        1739             :        if (len_trim(cosp_status(i)) > 0) then
        1740             :           write(iulog,*) "cosp_simulator: ERROR: "//trim(cosp_status(i))
        1741             :           nerror = nerror + 1
        1742             :        end if
        1743             :     end do
        1744             :     if (nerror > 0) then
        1745             :        call endrun('cospsimulator_intr_run: error return from cosp_simulator')
        1746             :     end if
        1747             :     call t_stopf("cosp_simulator")
        1748             :
        1749             :     ! ######################################################################################
        1750             :     ! Write COSP inputs to output file for offline use.
        1751             :     ! ######################################################################################
        1752             :     call t_startf("cosp_histfile_aux")
        1753             :     if (cosp_histfile_aux) then
        1754             :        ! 1D outputs
        1755             :        call outfld('PS_COSP',        state%ps(1:ncol),             ncol,lchnk)
        1756             :        call outfld('TS_COSP',        cospstateIN%skt,              ncol,lchnk)
        1757             :
        1758             :        ! 2D outputs
        1759             :        call outfld('P_COSP',         cospstateIN%pfull,            ncol,lchnk)
        1760             :        call outfld('PH_COSP',        cospstateIN%phalf,            ncol,lchnk)
        1761             :        call outfld('ZLEV_COSP',      cospstateIN%hgt_matrix,       ncol,lchnk)
        1762             :        call outfld('ZLEV_HALF_COSP', cospstateIN%hgt_matrix_half,  ncol,lchnk)
        1763             :        call outfld('T_COSP',         cospstateIN%at,               ncol,lchnk)
        1764             :        call outfld('Q_COSP',         cospstateIN%qv,               ncol,lchnk)
        1765             :
        1766             :        ! 3D outputs, but first compress to 2D
        1767             :        do i=1,ncol
        1768             :           do ihml=1,nlay
        1769             :              do isc=1,nscol_cosp
        1770             :                 ihsc = (ihml-1)*nscol_cosp+isc
        1771             :                 tau067_out(i,ihsc)  = cospIN%tau_067(i,isc,ihml)
        1772             :                 emis11_out(i,ihsc)  = cospIN%emiss_11(i,isc,ihml)
        1773             :                 ssa34_out(i,ihsc)   = cospIN%ss_alb(i,isc,ihml)
        1774             :                 asym34_out(i,ihsc)  = cospIN%asym(i,isc,ihml)
        1775             :                 fracLiq_out(i,ihsc) = cospIN%fracLiq(i,isc,ihml)
        1776             :              end do
        1777             :           end do
        1778             :        end do
        1779             :        call outfld('TAU_067',      tau067_out, pcols,lchnk)
        1780             :        call outfld('EMISS_11',     emis11_out, pcols,lchnk)
        1781             :        call outfld('MODIS_asym',   asym34_out, pcols,lchnk)
        1782             :        call outfld('MODIS_ssa',    ssa34_out,  pcols,lchnk)
        1783             :        call outfld('MODIS_fracliq',fracLiq_out,pcols,lchnk)
        1784             :     end if
        1785             :     call t_stopf("cosp_histfile_aux")
        1786             :
        1787             :     ! ######################################################################################
        1788             :     ! Set dark-scenes to fill value. Only done for passive simulators and when cosp_runall=F
        1789             :     ! ######################################################################################
        1790             :     call t_startf("sunlit_passive")
        1791             :     if (.not. cosp_runall) then
        1792             :        ! ISCCP simulator
        1793             :        if (lisccp_sim) then
        1794             :           ! 1D
        1795             :           where(cam_sunlit(1:ncol) .eq. 0)
        1796             :              cospOUT%isccp_totalcldarea(1:ncol)  = R_UNDEF
        1797             :              cospOUT%isccp_meanptop(1:ncol)      = R_UNDEF
        1798             :              cospOUT%isccp_meantaucld(1:ncol)    = R_UNDEF
        1799             :              cospOUT%isccp_meanalbedocld(1:ncol) = R_UNDEF
        1800             :              cospOUT%isccp_meantb(1:ncol)        = R_UNDEF
        1801             :              cospOUT%isccp_meantbclr(1:ncol)     = R_UNDEF
        1802             :           end where
        1803             :           ! 2D
        1804             :           do i=1,nscol_cosp
        1805             :              where (cam_sunlit(1:ncol) .eq. 0)
        1806             :                 cospOUT%isccp_boxtau(1:ncol,i)  = R_UNDEF
        1807             :                 cospOUT%isccp_boxptop(1:ncol,i) = R_UNDEF
        1808             :              end where
        1809             :           enddo
        1810             :           ! 3D
        1811             :           do i=1,nprs_cosp
        1812             :              do k=1,ntau_cosp
        1813             :                 where(cam_sunlit(1:ncol) .eq. 0)
        1814             :                    cospOUT%isccp_fq(1:ncol,k,i) = R_UNDEF
        1815             :                 end where
        1816             :              end do
        1817             :           end do
        1818             :        endif
        1819             :
        1820             :        ! MISR simulator
        1821             :        if (lmisr_sim) then
        1822             :           do i=1,nhtmisr_cosp
        1823             :              do k=1,ntau_cosp
        1824             :                 where(cam_sunlit(1:ncol) .eq. 0)
        1825             :                    cospOUT%misr_fq(1:ncol,k,i) = R_UNDEF
        1826             :                 end where
        1827             :              end do
        1828             :           end do
        1829             :        end if
        1830             :
        1831             :        ! MODIS simulator
        1832             :        if (lmodis_sim) then
        1833             :           ! 1D
        1834             :           where(cam_sunlit(1:ncol) .eq. 0)
        1835             :              cospOUT%modis_Cloud_Fraction_Total_Mean(1:ncol)       = R_UNDEF
        1836             :              cospOUT%modis_Cloud_Fraction_Water_Mean(1:ncol)       = R_UNDEF
        1837             :              cospOUT%modis_Cloud_Fraction_Ice_Mean(1:ncol)         = R_UNDEF
        1838             :              cospOUT%modis_Cloud_Fraction_High_Mean(1:ncol)        = R_UNDEF
        1839             :              cospOUT%modis_Cloud_Fraction_Mid_Mean(1:ncol)         = R_UNDEF
        1840             :              cospOUT%modis_Cloud_Fraction_Low_Mean(1:ncol)         = R_UNDEF
        1841             :              cospOUT%modis_Optical_Thickness_Total_Mean(1:ncol)    = R_UNDEF
        1842             :              cospOUT%modis_Optical_Thickness_Water_Mean(1:ncol)    = R_UNDEF
        1843             :              cospOUT%modis_Optical_Thickness_Ice_Mean(1:ncol)      = R_UNDEF
        1844             :              cospOUT%modis_Optical_Thickness_Total_LogMean(1:ncol) = R_UNDEF
        1845             :              cospOUT%modis_Optical_Thickness_Water_LogMean(1:ncol) = R_UNDEF
        1846             :              cospOUT%modis_Optical_Thickness_Ice_LogMean(1:ncol)   = R_UNDEF
        1847             :              cospOUT%modis_Cloud_Particle_Size_Water_Mean(1:ncol)  = R_UNDEF
        1848             :              cospOUT%modis_Cloud_Particle_Size_Ice_Mean(1:ncol)    = R_UNDEF
        1849             :              cospOUT%modis_Cloud_Top_Pressure_Total_Mean(1:ncol)   = R_UNDEF
        1850             :              cospOUT%modis_Liquid_Water_Path_Mean(1:ncol)          = R_UNDEF
        1851             :              cospOUT%modis_Ice_Water_Path_Mean(1:ncol)             = R_UNDEF
        1852             :           endwhere
        1853             :           ! 3D
        1854             :           do i=1,ntau_cosp_modis
        1855             :              do k=1,nprs_cosp
        1856             :                 where(cam_sunlit(1:ncol) .eq. 0)
        1857             :                    cospOUT%modis_Optical_Thickness_vs_Cloud_Top_Pressure(1:ncol,i,k) = R_UNDEF
        1858             :                 end where
        1859             :              enddo
        1860             :              do k=1,numMODISReffIceBins
        1861             :                 where(cam_sunlit(1:ncol) .eq. 0)
        1862             :                    cospOUT%modis_Optical_Thickness_vs_ReffICE(1:ncol,i,k) = R_UNDEF
        1863             :                 end where
        1864             :              end do
        1865             :              do k=1,numMODISReffLiqBins
        1866             :                 where(cam_sunlit(1:ncol) .eq. 0)
        1867             :                    cospOUT%modis_Optical_Thickness_vs_ReffLIQ(1:ncol,i,k) = R_UNDEF
        1868             :                 end where
        1869             :              enddo
        1870             :           enddo
        1871             :        end if
        1872             :     end if
        1873             :     call t_stopf("sunlit_passive")
        1874             :
        1875             :     ! ######################################################################################
        1876             :     ! Copy COSP outputs to CAM fields.
        1877             :     ! ######################################################################################
        1878             :     call t_startf("output_copying")
        1879             :     if (allocated(cospIN%frac_out)) &
        1880             :          frac_out(1:ncol,1:nscol_cosp,1:nlay) = cospIN%frac_out
        1881             :
        1882             :     ! Cloudsat
        1883             :     if (lradar_sim) then
        1884             :        cfad_dbze94(1:ncol,1:CLOUDSAT_DBZE_BINS,1:nht_cosp) = cospOUT%cloudsat_cfad_ze
        1885             :        dbze94(1:ncol,1:nscol_cosp,1:nlay)    = cospOUT%cloudsat_Ze_tot
        1886             :        cldtot_cs(1:ncol)  = 0._r8
        1887             :        cldtot_cs2(1:ncol) = 0._r8
        1888             :        ! *NOTE* These two fields are joint-simulator products, but in CAM they are controlled
        1889             :        !        by the radar simulator control.
        1890             :        cldtot_calcs(1:ncol) = cospOUT%radar_lidar_tcc
        1891             :        cld_cal_notcs(1:ncol,1:nht_cosp) = cospOUT%lidar_only_freq_cloud
        1892             :
        1893             :        ! Cloudsat near-surface precipitation diagnostics
        1894             :        ptcloudsatflag0(1:ncol) = cospOUT%cloudsat_precip_cover(:,1)
        1895             :        ptcloudsatflag1(1:ncol) = cospOUT%cloudsat_precip_cover(:,2)
        1896             :        ptcloudsatflag2(1:ncol) = cospOUT%cloudsat_precip_cover(:,3)
        1897             :        ptcloudsatflag3(1:ncol) = cospOUT%cloudsat_precip_cover(:,4)
        1898             :        ptcloudsatflag4(1:ncol) = cospOUT%cloudsat_precip_cover(:,5)
        1899             :        ptcloudsatflag5(1:ncol) = cospOUT%cloudsat_precip_cover(:,6)
        1900             :        ptcloudsatflag6(1:ncol) = cospOUT%cloudsat_precip_cover(:,7)
        1901             :        ptcloudsatflag7(1:ncol) = cospOUT%cloudsat_precip_cover(:,8)
        1902             :        ptcloudsatflag8(1:ncol) = cospOUT%cloudsat_precip_cover(:,9)
        1903             :        ptcloudsatflag9(1:ncol) = cospOUT%cloudsat_precip_cover(:,10)
        1904             :        cloudsatpia(1:ncol)     = cospOUT%cloudsat_pia
        1905             :
        1906             :     endif
        1907             :
        1908             :     ! CALIPSO
        1909             :     if (llidar_sim) then
        1910             :        cldlow_cal(1:ncol)                = cospOUT%calipso_cldlayer(:,1)
        1911             :        cldmed_cal(1:ncol)                = cospOUT%calipso_cldlayer(:,2)
        1912             :        cldhgh_cal(1:ncol)                = cospOUT%calipso_cldlayer(:,3)
        1913             :        cldtot_cal(1:ncol)                = cospOUT%calipso_cldlayer(:,4)
        1914             :        cldlow_cal_ice(1:ncol)            = cospOUT%calipso_cldlayerphase(:,1,1)
        1915             :        cldmed_cal_ice(1:ncol)            = cospOUT%calipso_cldlayerphase(:,2,1)
        1916             :        cldhgh_cal_ice(1:ncol)            = cospOUT%calipso_cldlayerphase(:,3,1)
        1917             :        cldtot_cal_ice(1:ncol)            = cospOUT%calipso_cldlayerphase(:,4,1)
        1918             :        cldlow_cal_liq(1:ncol)            = cospOUT%calipso_cldlayerphase(:,1,2)
        1919             :        cldmed_cal_liq(1:ncol)            = cospOUT%calipso_cldlayerphase(:,2,2)
        1920             :        cldhgh_cal_liq(1:ncol)            = cospOUT%calipso_cldlayerphase(:,3,2)
        1921             :        cldtot_cal_liq(1:ncol)            = cospOUT%calipso_cldlayerphase(:,4,2)
        1922             :        cldlow_cal_un(1:ncol)             = cospOUT%calipso_cldlayerphase(:,1,3)
        1923             :        cldmed_cal_un(1:ncol)             = cospOUT%calipso_cldlayerphase(:,2,3)
        1924             :        cldhgh_cal_un(1:ncol)             = cospOUT%calipso_cldlayerphase(:,3,3)
        1925             :        cldtot_cal_un(1:ncol)             = cospOUT%calipso_cldlayerphase(:,4,3)
        1926             :        cld_cal_ice(1:ncol,1:nht_cosp)    = cospOUT%calipso_lidarcldphase(:,:,1)
        1927             :        cld_cal_liq(1:ncol,1:nht_cosp)    = cospOUT%calipso_lidarcldphase(:,:,2)
        1928             :        cld_cal_un(1:ncol,1:nht_cosp)     = cospOUT%calipso_lidarcldphase(:,:,3)
        1929             :        cld_cal_tmp(1:ncol,1:nht_cosp)    = cospOUT%calipso_lidarcldtmp(:,:,1)
        1930             :        cld_cal_tmpliq(1:ncol,1:nht_cosp) = cospOUT%calipso_lidarcldtmp(:,:,2)
        1931             :        cld_cal_tmpice(1:ncol,1:nht_cosp) = cospOUT%calipso_lidarcldtmp(:,:,3)
        1932             :        cld_cal_tmpun(1:ncol,1:nht_cosp)  = cospOUT%calipso_lidarcldtmp(:,:,4)
        1933             :        cld_cal(1:ncol,1:nht_cosp)        = cospOUT%calipso_lidarcld(:,1:nht_cosp)
        1934             :        mol532_cal(1:ncol,1:nlay)         = cospOUT%calipso_beta_mol
        1935             :        atb532(1:ncol,1:nscol_cosp,1:nlay)= cospOUT%calipso_beta_tot
        1936             :        cfad_lidarsr532(1:ncol,1:nsr_cosp,1:nht_cosp) = cospOUT%calipso_cfad_sr(:,:,:)
        1937             :        refl_parasol(1:ncol,1:nsza_cosp)  = cospOUT%parasolGrid_refl
        1938             :     endif
        1939             :
        1940             :     ! ISCCP
        1941             :     if (lisccp_sim) then
        1942             :        clisccp2(1:ncol,1:ntau_cosp,1:nprs_cosp) = cospOUT%isccp_fq
        1943             :        tau_isccp(1:ncol,1:nscol_cosp)           = cospOUT%isccp_boxtau
        1944             :        cldptop_isccp(1:ncol,1:nscol_cosp)       = cospOUT%isccp_boxptop
        1945             :        cldtot_isccp(1:ncol)                     = cospOUT%isccp_totalcldarea
        1946             :        meanptop_isccp(1:ncol)                   = cospOUT%isccp_meanptop
        1947             :        meantau_isccp(1:ncol)                    = cospOUT%isccp_meantaucld
        1948             :        meancldalb_isccp(1:ncol)                 = cospOUT%isccp_meanalbedocld
        1949             :        meantb_isccp(1:ncol)                     = cospOUT%isccp_meantb
        1950             :        meantbclr_isccp(1:ncol)                  = cospOUT%isccp_meantbclr
        1951             :     endif
        1952             :
        1953             :     ! MISR
        1954             :     if (lmisr_sim) then
        1955             :        clMISR(1:ncol,1:ntau_cosp,1:nhtmisr_cosp) = cospOUT%misr_fq
        1956             :     endif
        1957             :
        1958             :     ! MODIS
        1959             :     if (lmodis_sim) then
        1960             :        cltmodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_Total_Mean
        1961             :        clwmodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_Water_Mean
        1962             :        climodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_Ice_Mean
        1963             :        clhmodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_High_Mean
        1964             :        clmmodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_Mid_Mean
        1965             :        cllmodis(1:ncol)     = cospOUT%modis_Cloud_Fraction_Low_Mean
        1966             :        tautmodis(1:ncol)    = cospOUT%modis_Optical_Thickness_Total_Mean
        1967             :        tauwmodis(1:ncol)    = cospOUT%modis_Optical_Thickness_Water_Mean
        1968             :        tauimodis(1:ncol)    = cospOUT%modis_Optical_Thickness_Ice_Mean
        1969             :        tautlogmodis(1:ncol) = cospOUT%modis_Optical_Thickness_Total_LogMean
        1970             :        tauwlogmodis(1:ncol) = cospOUT%modis_Optical_Thickness_Water_LogMean
        1971             :        tauilogmodis(1:ncol) = cospOUT%modis_Optical_Thickness_Ice_LogMean
        1972             :        reffclwmodis(1:ncol) = cospOUT%modis_Cloud_Particle_Size_Water_Mean
        1973             :        reffclimodis(1:ncol) = cospOUT%modis_Cloud_Particle_Size_Ice_Mean
        1974             :        pctmodis(1:ncol)     = cospOUT%modis_Cloud_Top_Pressure_Total_Mean
        1975             :        lwpmodis(1:ncol)     = cospOUT%modis_Liquid_Water_Path_Mean
        1976             :        iwpmodis(1:ncol)     = cospOUT%modis_Ice_Water_Path_Mean
        1977             :        clmodis(1:ncol,1:ntau_cosp_modis,1:nprs_cosp)  = cospOUT%modis_Optical_Thickness_vs_Cloud_Top_Pressure
        1978             :        clrimodis(1:ncol,1:ntau_cosp_modis,1:numMODISReffIceBins) = cospOUT%modis_Optical_Thickness_vs_ReffICE
        1979             :        clrlmodis(1:ncol,1:ntau_cosp_modis,1:numMODISReffLiqBins) = cospOUT%modis_Optical_Thickness_vs_ReffLIQ
        1980             :     endif
        1981             :
        1982             :     ! Use COSP output to populate CAM collapsed output variables
        1983             :     do i=1,ncol
        1984             :        if (lradar_sim) then
        1985             :           do ih=1,nht_cosp
        1986             :              do id=1,CLOUDSAT_DBZE_BINS
        1987             :                 ihd=(ih-1)*CLOUDSAT_DBZE_BINS+id
        1988             :                 cfad_dbze94_cs(i,ihd) = cfad_dbze94(i,id,ih)
        1989             :              end do
        1990             :           end do
        1991             :           do ihml=1,nlay
        1992             :              do isc=1,nscol_cosp
        1993             :                 ihsc=(ihml-1)*nscol_cosp+isc
        1994             :                 dbze_cs(i,ihsc) = dbze94(i,isc,ihml)
        1995             :              end do
        1996             :           end do
        1997             :        endif
        1998             :
        1999             :        if (llidar_sim) then
        2000             :           do ih=1,nht_cosp
        2001             :              do is=1,nsr_cosp
        2002             :                 ihs=(ih-1)*nsr_cosp+is
        2003             :                 cfad_sr532_cal(i,ihs) = cfad_lidarsr532(i,is,ih)
        2004             :              end do
        2005             :           end do
        2006             :           do ihml=1,nlay
        2007             :              do isc=1,nscol_cosp
        2008             :                 ihsc=(ihml-1)*nscol_cosp+isc
        2009             :                 atb532_cal(i,ihsc) = atb532(i,isc,ihml)
        2010             :              end do
        2011             :           end do
        2012             :        endif
        2013             :
        2014             :        if (lmisr_sim) then
        2015             :           do ihm=1,nhtmisr_cosp
        2016             :              do it=1,ntau_cosp
        2017             :                 ihmt=(ihm-1)*ntau_cosp+it
        2018             :                 cld_misr(i,ihmt) = clMISR(i,it,ihm)
        2019             :              end do
        2020             :           end do
        2021             :        endif
        2022             :
        2023             :        if (lmodis_sim) then
        2024             :           do ip=1,nprs_cosp
        2025             :              do it=1,ntau_cosp_modis
        2026             :                 ipt=(ip-1)*ntau_cosp_modis+it
        2027             :                 clmodis_cam(i,ipt) = clmodis(i,it,ip)
        2028             :              end do
        2029             :           end do
        2030             :           do ip=1,numMODISReffIceBins
        2031             :              do it=1,ntau_cosp_modis
        2032             :                 ipt=(ip-1)*ntau_cosp_modis+it
        2033             :                 clrimodis_cam(i,ipt) = clrimodis(i,it,ip)
        2034             :              end do
        2035             :           end do
        2036             :           do ip=1,numMODISReffLiqBins
        2037             :              do it=1,ntau_cosp_modis
        2038             :                 ipt=(ip-1)*ntau_cosp_modis+it
        2039             :                 clrlmodis_cam(i,ipt) = clrlmodis(i,it,ip)
        2040             :              end do
        2041             :           end do
        2042             :        endif
        2043             :
        2044             :        ! Subcolums
        2045             :        do ihml=1,nlay
        2046             :           do isc=1,nscol_cosp
        2047             :              ihsc=(ihml-1)*nscol_cosp+isc
        2048             :              scops_out(i,ihsc) = frac_out(i,isc,ihml)
        2049             :           end do
        2050             :        end do
        2051             :     end do
        2052             :     call t_stopf("output_copying")
        2053             :
        2054             :     ! ######################################################################################
        2055             :     ! Clean up
        2056             :     ! ######################################################################################
        2057             :     call t_startf("destroy_cospIN")
        2058             :     call destroy_cospIN(cospIN)
        2059             :     call t_stopf("destroy_cospIN")
        2060             :     call t_startf("destroy_cospstateIN")
        2061             :     call destroy_cospstateIN(cospstateIN)
        2062             :     call t_stopf("destroy_cospstateIN")
        2063             :     call t_startf("destroy_cospOUT")
        2064             :     call destroy_cosp_outputs(cospOUT)
        2065             :     call t_stopf("destroy_cospOUT")
        2066             :
        2067             :     ! ######################################################################################
        2068             :     ! OUTPUT
        2069             :     ! ######################################################################################
        2070             :     call t_startf("writing_output")
        2071             :     ! ISCCP OUTPUTS
        2072             :     if (lisccp_sim) then
        2073             :        call outfld('FISCCP1_COSP',clisccp2,     pcols,lchnk)
        2074             :        call outfld('CLDTOT_ISCCP',cldtot_isccp, pcols,lchnk)
        2075             :        !! weight meancldalb_isccp by the cloud fraction
        2076             :        !! where there is no isccp cloud fraction, set meancldalb_isccp = R_UNDEF
        2077             :        !! weight meanptop_isccp  by the cloud fraction
        2078             :        !! where there is no isccp cloud fraction, set meanptop_isccp = R_UNDEF
        2079             :        !! weight meantau_isccp by the cloud fraction
        2080             :        !! where there is no isccp cloud fraction, set meantau_isccp = R_UNDEF
        2081             :        where (cldtot_isccp(:ncol) .eq. R_UNDEF)
        2082             :           meancldalb_isccp(:ncol) = R_UNDEF
        2083             :           meanptop_isccp(:ncol)   = R_UNDEF
        2084             :           meantau_isccp(:ncol)    = R_UNDEF
        2085             :        elsewhere
        2086             :           meancldalb_isccp(:ncol) = meancldalb_isccp(:ncol)*cldtot_isccp(:ncol)
        2087             :           meanptop_isccp(:ncol)   = meanptop_isccp(:ncol)*cldtot_isccp(:ncol)
        2088             :           meantau_isccp(:ncol)    = meantau_isccp(:ncol)*cldtot_isccp(:ncol)
        2089             :        end where
        2090             :        call outfld('MEANCLDALB_ISCCP',meancldalb_isccp,pcols,lchnk)
        2091             :        call outfld('MEANPTOP_ISCCP',  meanptop_isccp,  pcols,lchnk)
        2092             :        call outfld('MEANTAU_ISCCP',   meantau_isccp,   pcols,lchnk)
        2093             :        call outfld('MEANTB_ISCCP',    meantb_isccp,    pcols,lchnk)
        2094             :        call outfld('MEANTBCLR_ISCCP', meantbclr_isccp, pcols,lchnk)
        2095             :     end if
        2096             :
        2097             :     ! CALIPSO SIMULATOR OUTPUTS
        2098             :     if (llidar_sim) then
        2099             :        call outfld('CLDLOW_CAL',    cldlow_cal,     pcols,lchnk)
        2100             :        call outfld('CLDMED_CAL',    cldmed_cal,     pcols,lchnk)
        2101             :        call outfld('CLDHGH_CAL',    cldhgh_cal,     pcols,lchnk)
        2102             :        call outfld('CLDTOT_CAL',    cldtot_cal,     pcols,lchnk)
        2103             :        call outfld('CLDTOT_CAL_ICE',cldtot_cal_ice, pcols,lchnk) !+1.4
        2104             :        call outfld('CLDTOT_CAL_LIQ',cldtot_cal_liq, pcols,lchnk)
        2105             :        call outfld('CLDTOT_CAL_UN', cldtot_cal_un,  pcols,lchnk)
        2106             :        call outfld('CLDHGH_CAL_ICE',cldhgh_cal_ice, pcols,lchnk)
        2107             :        call outfld('CLDHGH_CAL_LIQ',cldhgh_cal_liq, pcols,lchnk)
        2108             :        call outfld('CLDHGH_CAL_UN', cldhgh_cal_un,  pcols,lchnk)
        2109             :        call outfld('CLDMED_CAL_ICE',cldmed_cal_ice, pcols,lchnk)
        2110             :        call outfld('CLDMED_CAL_LIQ',cldmed_cal_liq, pcols,lchnk)
        2111             :        call outfld('CLDMED_CAL_UN', cldmed_cal_un,  pcols,lchnk)
        2112             :        call outfld('CLDLOW_CAL_ICE',cldlow_cal_ice, pcols,lchnk)
        2113             :        call outfld('CLDLOW_CAL_LIQ',cldlow_cal_liq, pcols,lchnk)
        2114             :        call outfld('CLDLOW_CAL_UN', cldlow_cal_un,  pcols,lchnk) !+1.4
        2115             :        where (cld_cal(:ncol,:nht_cosp) .eq. R_UNDEF)
        2116             :           !! setting missing values to 0 (clear air).
        2117             :           !! I'm not sure why COSP produces a mix of R_UNDEF and realvalue in the nht_cosp dimension.
        2118             :           cld_cal(:ncol,:nht_cosp) = 0.0_r8
        2119             :        end where
        2120             :        call outfld('CLD_CAL',        cld_cal,       pcols,lchnk)  !! fails check_accum if 'A'
        2121             :        call outfld('MOL532_CAL',     mol532_cal,    pcols,lchnk)
        2122             :
        2123             :        where (cfad_sr532_cal(:ncol,:nht_cosp*nsr_cosp) .eq. R_UNDEF)
        2124             :           !! fails check_accum if this is set... with ht_cosp set relative to sea level, mix of R_UNDEF and realvalue
        2125             :           !!            cfad_sr532_cal(:ncol,:nht_cosp*nsr_cosp) = R_UNDEF
        2126             :           cfad_sr532_cal(:ncol,:nht_cosp*nsr_cosp) = 0.0_r8
        2127             :        end where
        2128             :        call outfld('CFAD_SR532_CAL',cfad_sr532_cal    ,pcols,lchnk)
        2129             :
        2130             :        where (refl_parasol(:ncol,:nsza_cosp) .eq. R_UNDEF)
        2131             :           !! setting missing values to 0 (clear air).
        2132             :           refl_parasol(:ncol,:nsza_cosp) = 0
        2133             :        end where
        2134             :        call outfld('RFL_PARASOL',refl_parasol   ,pcols,lchnk) !!
        2135             :
        2136             :        where (cld_cal_liq(:ncol,:nht_cosp) .eq. R_UNDEF) !+cosp1.4
        2137             :           !! setting missing values to 0 (clear air), likely below sea level
        2138             :           cld_cal_liq(:ncol,:nht_cosp) = 0.0_r8
        2139             :        end where
        2140             :        call outfld('CLD_CAL_LIQ',cld_cal_liq    ,pcols,lchnk)  !!
        2141             :
        2142             :        where (cld_cal_ice(:ncol,:nht_cosp) .eq. R_UNDEF)
        2143             :           !! setting missing values to 0 (clear air), likely below sea level
        2144             :           cld_cal_ice(:ncol,:nht_cosp) = 0.0_r8
        2145             :        end where
        2146             :        call outfld('CLD_CAL_ICE',cld_cal_ice    ,pcols,lchnk)  !!
        2147             :
        2148             :        where (cld_cal_un(:ncol,:nht_cosp) .eq. R_UNDEF)
        2149             :           !! setting missing values to 0 (clear air), likely below sea level
        2150             :           cld_cal_un(:ncol,:nht_cosp) = 0.0_r8
        2151             :        end where
        2152             :        call outfld('CLD_CAL_UN',cld_cal_un    ,pcols,lchnk)  !!
        2153             :
        2154             :        where (cld_cal_tmp(:ncol,:nht_cosp) .eq. R_UNDEF)
        2155             :           !! setting missing values to 0 (clear air), likely below sea level
        2156             :           cld_cal_tmp(:ncol,:nht_cosp) = 0.0_r8
        2157             :        end where
        2158             :        call outfld('CLD_CAL_TMP',cld_cal_tmp    ,pcols,lchnk)  !!
        2159             :
        2160             :        where (cld_cal_tmpliq(:ncol,:nht_cosp) .eq. R_UNDEF)
        2161             :           !! setting missing values to 0 (clear air), likely below sea level
        2162             :           cld_cal_tmpliq(:ncol,:nht_cosp) = 0.0_r8
        2163             :        end where
        2164             :        call outfld('CLD_CAL_TMPLIQ',cld_cal_tmpliq    ,pcols,lchnk)  !!
        2165             :
        2166             :        where (cld_cal_tmpice(:ncol,:nht_cosp) .eq. R_UNDEF)
        2167             :           !! setting missing values to 0 (clear air), likely below sea level
        2168             :           cld_cal_tmpice(:ncol,:nht_cosp) = 0.0_r8
        2169             :        end where
        2170             :        call outfld('CLD_CAL_TMPICE',cld_cal_tmpice    ,pcols,lchnk)  !!
        2171             :
        2172             :        where (cld_cal_tmpun(:ncol,:nht_cosp) .eq. R_UNDEF)
        2173             :           !! setting missing values to 0 (clear air), likely below sea level
        2174             :           cld_cal_tmpun(:ncol,:nht_cosp) = 0.0_r8
        2175             :        end where
        2176             :        call outfld('CLD_CAL_TMPUN',cld_cal_tmpun    ,pcols,lchnk)  !!  !+cosp1.4
        2177             :
        2178             :     end if
        2179             :
        2180             :     ! RADAR SIMULATOR OUTPUTS
        2181             :     if (lradar_sim) then
        2182             :        where (cfad_dbze94_cs(:ncol,:nht_cosp*CLOUDSAT_DBZE_BINS) .eq. R_UNDEF)
        2183             :           !! fails check_accum if this is set... with ht_cosp set relative to sea level, mix of R_UNDEF and realvalue
        2184             :           !           cfad_dbze94_cs(:ncol,:nht_cosp*CLOUDSAT_DBZE_BINS) = R_UNDEF
        2185             :           cfad_dbze94_cs(:ncol,:nht_cosp*CLOUDSAT_DBZE_BINS) = 0.0_r8
        2186             :        end where
        2187             :        call outfld('CFAD_DBZE94_CS',cfad_dbze94_cs,   pcols, lchnk)
        2188             :        call outfld('CLDTOT_CALCS',  cldtot_calcs,     pcols, lchnk)
        2189             :        call outfld('CLDTOT_CS',     cldtot_cs,        pcols, lchnk)
        2190             :        call outfld('CLDTOT_CS2',    cldtot_cs2,       pcols, lchnk)
        2191             :        call outfld('CLD_CAL_NOTCS', cld_cal_notcs,    pcols, lchnk)
        2192             :        call outfld('CS_NOPRECIP',   ptcloudsatflag0,  pcols, lchnk)
        2193             :        call outfld('CS_RAINPOSS',   ptcloudsatflag1,  pcols, lchnk)
        2194             :        call outfld('CS_RAINPROB',   ptcloudsatflag2,  pcols, lchnk)
        2195             :        call outfld('CS_RAINCERT',   ptcloudsatflag3,  pcols, lchnk)
        2196             :        call outfld('CS_SNOWPOSS',   ptcloudsatflag4,  pcols, lchnk)
        2197             :        call outfld('CS_SNOWCERT',   ptcloudsatflag5,  pcols, lchnk)
        2198             :        call outfld('CS_MIXPOSS',    ptcloudsatflag6,  pcols, lchnk)
        2199             :        call outfld('CS_MIXCERT',    ptcloudsatflag7,  pcols, lchnk)
        2200             :        call outfld('CS_RAINHARD',   ptcloudsatflag8,  pcols, lchnk)
        2201             :        call outfld('CS_UN',         ptcloudsatflag9,  pcols, lchnk)
        2202             :        call outfld('CS_PIA',        cloudsatpia,      pcols, lchnk)
        2203             :     end if
        2204             :
        2205             :     ! MISR SIMULATOR OUTPUTS
        2206             :     if (lmisr_sim) then
        2207             :        call outfld('CLD_MISR',cld_misr    ,pcols,lchnk)
        2208             :     end if
        2209             :
        2210             :     ! MODIS SIMULATOR OUTPUTS
        2211             :     if (lmodis_sim) then
        2212             :        call outfld('CLTMODIS',cltmodis    ,pcols,lchnk)
        2213             :        call outfld('CLWMODIS',clwmodis    ,pcols,lchnk)
        2214             :        call outfld('CLIMODIS',climodis    ,pcols,lchnk)
        2215             :        call outfld('CLHMODIS',clhmodis    ,pcols,lchnk)
        2216             :        call outfld('CLMMODIS',clmmodis    ,pcols,lchnk)
        2217             :        call outfld('CLLMODIS',cllmodis    ,pcols,lchnk)
        2218             :
        2219             :        !! where there is no cloud fraction or no retrieval, set to R_UNDEF,
        2220             :        !! otherwise weight retrieval by cloud fraction
        2221             :        where ((cltmodis(:ncol) .eq. R_UNDEF) .or. (tautmodis(:ncol) .eq. R_UNDEF))
        2222             :           tautmodis(:ncol) = R_UNDEF
        2223             :        elsewhere
        2224             :           !! weight by the cloud fraction cltmodis
        2225             :           tautmodis(:ncol) = tautmodis(:ncol)*cltmodis(:ncol)
        2226             :        end where
        2227             :        call outfld('TAUTMODIS',tautmodis    ,pcols,lchnk)
        2228             :
        2229             :        where ((tauwmodis(:ncol) .eq. R_UNDEF) .or. (clwmodis(:ncol) .eq. R_UNDEF))
        2230             :           tauwmodis(:ncol) = R_UNDEF
        2231             :        elsewhere
        2232             :           !! weight by the cloud fraction clwmodis
        2233             :           tauwmodis(:ncol) = tauwmodis(:ncol)*clwmodis(:ncol)
        2234             :        end where
        2235             :        call outfld('TAUWMODIS',tauwmodis    ,pcols,lchnk)
        2236             :
        2237             :        where ((tauimodis(:ncol) .eq. R_UNDEF) .or. (climodis(:ncol) .eq. R_UNDEF))
        2238             :           tauimodis(:ncol) = R_UNDEF
        2239             :        elsewhere
        2240             :           !! weight by the cloud fraction climodis
        2241             :           tauimodis(:ncol) = tauimodis(:ncol)*climodis(:ncol)
        2242             :        end where
        2243             :        call outfld('TAUIMODIS',tauimodis    ,pcols,lchnk)
        2244             :
        2245             :        where ((tautlogmodis(:ncol)  .eq. R_UNDEF) .or. (cltmodis(:ncol) .eq. R_UNDEF))
        2246             :           tautlogmodis(:ncol) = R_UNDEF
        2247             :        elsewhere
        2248             :           !! weight by the cloud fraction cltmodis
        2249             :           tautlogmodis(:ncol) = tautlogmodis(:ncol)*cltmodis(:ncol)
        2250             :        end where
        2251             :        call outfld('TAUTLOGMODIS',tautlogmodis    ,pcols,lchnk)
        2252             :
        2253             :        where ((tauwlogmodis(:ncol)  .eq. R_UNDEF) .or. (clwmodis(:ncol) .eq. R_UNDEF))
        2254             :           tauwlogmodis(:ncol) = R_UNDEF
        2255             :        elsewhere
        2256             :           !! weight by the cloud fraction clwmodis
        2257             :           tauwlogmodis(:ncol) = tauwlogmodis(:ncol)*clwmodis(:ncol)
        2258             :        end where
        2259             :        call outfld('TAUWLOGMODIS',tauwlogmodis    ,pcols,lchnk)
        2260             :
        2261             :        where ((tauilogmodis(:ncol)  .eq. R_UNDEF) .or. (climodis(:ncol) .eq. R_UNDEF))
        2262             :           tauilogmodis(:ncol) = R_UNDEF
        2263             :        elsewhere
        2264             :           !! weight by the cloud fraction climodis
        2265             :           tauilogmodis(:ncol) = tauilogmodis(:ncol)*climodis(:ncol)
        2266             :        end where
        2267             :        call outfld('TAUILOGMODIS',tauilogmodis    ,pcols,lchnk)
        2268             :
        2269             :        where ((reffclwmodis(:ncol)  .eq. R_UNDEF) .or. (clwmodis(:ncol) .eq. R_UNDEF))
        2270             :           reffclwmodis(:ncol) = R_UNDEF
        2271             :        elsewhere
        2272             :           !! weight by the cloud fraction clwmodis
        2273             :           reffclwmodis(:ncol) = reffclwmodis(:ncol)*clwmodis(:ncol)
        2274             :        end where
        2275             :        call outfld('REFFCLWMODIS',reffclwmodis    ,pcols,lchnk)
        2276             :
        2277             :        where ((reffclimodis(:ncol)  .eq. R_UNDEF) .or. (climodis(:ncol) .eq. R_UNDEF))
        2278             :           reffclimodis(:ncol) = R_UNDEF
        2279             :        elsewhere
        2280             :           !! weight by the cloud fraction climodis
        2281             :           reffclimodis(:ncol) = reffclimodis(:ncol)*climodis(:ncol)
        2282             :        end where
        2283             :        call outfld('REFFCLIMODIS',reffclimodis    ,pcols,lchnk)
        2284             :
        2285             :        where ((pctmodis(:ncol)  .eq. R_UNDEF) .or. ( cltmodis(:ncol) .eq. R_UNDEF))
        2286             :           pctmodis(:ncol) = R_UNDEF
        2287             :        elsewhere
        2288             :           !! weight by the cloud fraction cltmodis
        2289             :           pctmodis(:ncol) = pctmodis(:ncol)*cltmodis(:ncol)
        2290             :        end where
        2291             :        call outfld('PCTMODIS',pctmodis    ,pcols,lchnk)
        2292             :
        2293             :        where ((lwpmodis(:ncol)  .eq. R_UNDEF) .or. (clwmodis(:ncol) .eq. R_UNDEF))
        2294             :           lwpmodis(:ncol) = R_UNDEF
        2295             :        elsewhere
        2296             :           !! weight by the cloud fraction clwmodis
        2297             :           lwpmodis(:ncol) = lwpmodis(:ncol)*clwmodis(:ncol)
        2298             :        end where
        2299             :        call outfld('LWPMODIS',lwpmodis    ,pcols,lchnk)
        2300             :
        2301             :        where ((iwpmodis(:ncol)  .eq. R_UNDEF) .or. (climodis(:ncol) .eq. R_UNDEF))
        2302             :           iwpmodis(:ncol) = R_UNDEF
        2303             :        elsewhere
        2304             :           !! weight by the cloud fraction climodis
        2305             :           iwpmodis(:ncol) = iwpmodis(:ncol)*climodis(:ncol)
        2306             :        end where
        2307             :        call outfld('IWPMODIS',iwpmodis    ,pcols,lchnk)
        2308             :
        2309             :        call outfld('CLMODIS',clmodis_cam  ,pcols,lchnk)
        2310             :        call outfld('CLRIMODIS',clrimodis_cam  ,pcols,lchnk)
        2311             :        call outfld('CLRLMODIS',clrlmodis_cam  ,pcols,lchnk)
        2312             :     end if
        2313             :
        2314             :     ! SUB-COLUMN OUTPUT
        2315             :     if (lfrac_out) then
        2316             :        call outfld('SCOPS_OUT',scops_out   ,pcols,lchnk)!!!-1.00000E+30 !! fails check_accum if 'A'
        2317             :        if (lisccp_sim) then
        2318             :           call outfld('TAU_ISCCP',    tau_isccp,    pcols,lchnk) !! fails check_accum if 'A'
        2319             :           call outfld('CLDPTOP_ISCCP',cldptop_isccp,pcols,lchnk) !! fails check_accum if 'A'
        2320             :        end if
        2321             :        if (llidar_sim) then
        2322             :           call outfld('ATB532_CAL',atb532_cal,pcols,lchnk) !! fails check_accum if 'A'
        2323             :        end if
        2324             :        if (lradar_sim) then
        2325             :           call outfld('DBZE_CS',dbze_cs,pcols,lchnk) !! fails check_accum if 'A'
        2326             :        end if
        2327             :     end if
        2328             :     call t_stopf("writing_output")
        2329             : #endif
        2330           0 :   end subroutine cospsimulator_intr_run
        2331             :
        2332             : #ifdef USE_COSP
        2333             :   ! ######################################################################################
        2334             :   ! SUBROUTINE subsample_and_optics
        2335             :   ! ######################################################################################
        2336             :   subroutine subsample_and_optics(nPoints, nLevels, nColumns, nHydro,overlap,            &
        2337             :                                   lidar_ice_type, sd, tca, cca,                          &
        2338             :                                   fl_lsrainIN, fl_lssnowIN, fl_lsgrplIN, fl_ccrainIN,    &
        2339             :                                   fl_ccsnowIN, mr_lsliq, mr_lsice, mr_ccliq, mr_ccice,   &
        2340             :                                   reffIN, dtau_c, dtau_s, dem_c, dem_s, dtau_s_snow,     &
        2341             :                                   dem_s_snow, sfcP, cospstateIN, cospIN)
        2342             :     ! Dependencies
        2343             :     use cosp_kinds,           only: wp
        2344             :     use mod_rng,              only: rng_state, init_rng
        2345             :     use mod_cosp_config,      only: R_UNDEF
        2346             :     use mod_scops,            only: scops
        2347             :     use mod_prec_scops,       only: prec_scops
        2348             :     use mod_cosp_utils,       only: cosp_precip_mxratio
        2349             :     use mod_quickbeam_optics, only: quickbeam_optics, gases
        2350             :     use cosp_optics,          only: cosp_simulator_optics,lidar_optics,modis_optics,    &
        2351             :                                     modis_optics_partition
        2352             :     use mod_cosp_config,      only: Nlvgrid, vgrid_zl, vgrid_zu
        2353             :     use mod_cosp_stats,       only: cosp_change_vertical_grid
        2354             :     ! Inputs
        2355             :     integer,intent(in) :: &
        2356             :          nPoints,      & ! Number of gridpoints
        2357             :          nLevels,      & ! Number of vertical levels
        2358             :          nColumns,     & ! Number of subcolumns
        2359             :          nHydro,       & ! Number pf hydrometeor types
        2360             :          overlap,      & ! Overlap assumption (1/2/3)
        2361             :          lidar_ice_type  ! Ice type assumption used by lidar optics
        2362             :     real(wp),intent(in),dimension(nPoints,nLevels) :: &
        2363             :          tca,          & ! Total cloud amount (0-1)
        2364             :          cca,          & ! Convective cloud amount (0-1)
        2365             :          mr_lsliq,     & ! Mixing ratio (kg/kg)
        2366             :          mr_lsice,     & ! Mixing ratio (kg/kg)
        2367             :          mr_ccliq,     & ! Mixing ratio (kg/kg)
        2368             :          mr_ccice,     & ! Mixing ratio (kg/kg)
        2369             :          dtau_c,       & ! 0.67-micron optical depth (convective)
        2370             :          dtau_s,       & ! 0.67-micron optical depth (stratiform)
        2371             :          dem_c,        & ! 11-micron emissivity (convective)
        2372             :          dem_s,        & ! 11-micron emissivity (stratiform)
        2373             :          fl_lsrainIN,  & ! Precipitation flux
        2374             :          fl_lssnowIN,  & ! Precipitation flux
        2375             :          fl_lsgrplIN,  & ! Precipitation flux
        2376             :          fl_ccrainIN,  & ! Precipitation flux
        2377             :          fl_ccsnowIN     ! Precipitation flux
        2378             :     real(wp),intent(inout),dimension(nPoints,nLevels) :: &
        2379             :          dtau_s_snow,  & ! 0.67-micron optical depth (snow)
        2380             :          dem_s_snow      ! 11-micron emissivity (snow)
        2381             :     real(wp),intent(in),dimension(nPoints,nLevels,nHydro) :: &
        2382             :          reffIN          !
        2383             :     real(wp),intent(in),dimension(nPoints) :: &
        2384             :          sfcP            ! Surface pressure
        2385             :     type(size_distribution),intent(inout) :: &
        2386             :          sd
        2387             :
        2388             :     ! Outputs
        2389             :     type(cosp_optical_inputs),intent(inout) :: cospIN
        2390             :     type(cosp_column_inputs),intent(inout)  :: cospstateIN
        2391             :
        2392             :     ! Local variables
        2393             :     integer :: i, j, k, istat
        2394             :     real(wp),dimension(nPoints,nLevels)      :: column_frac_out,column_prec_out,         &
        2395             :                                                 fl_lsrain,fl_lssnow,fl_lsgrpl,fl_ccrain, &
        2396             :                                                 fl_ccsnow
        2397             :     real(wp),dimension(nPoints,nLevels,nHydro) :: ReffTemp
        2398             :     type(rng_state),allocatable,dimension(:) :: rngs  ! Seeds for random number generator
        2399             :     integer,dimension(:),allocatable         :: seed
        2400             :     real(wp),dimension(:,:),allocatable      :: ls_p_rate,cv_p_rate,frac_ls,frac_cv,     &
        2401             :                                                 prec_ls,prec_cv,g_vol
        2402             :     real(wp),dimension(:,:,:),  allocatable  :: frac_prec,&
        2403             :                                                  MODIS_cloudWater,MODIS_cloudIce,        &
        2404             :                                                  MODIS_watersize,MODIS_iceSize,          &
        2405             :                                                  MODIS_snowSize,MODIS_cloudSnow,         &
        2406             :                                                  MODIS_opticalThicknessLiq,              &
        2407             :                                                  MODIS_opticalThicknessSnow,             &
        2408             :                                                  MODIS_opticalThicknessIce,              &
        2409             :                                                  fracPrecipIce, fracPrecipIce_statGrid
        2410             :     real(wp),dimension(:,:,:,:),allocatable   :: mr_hydro,Reff,Np
        2411             :
        2412             :     character(len=*), parameter :: sub = 'subsample_and_optics'
        2413             :     !--------------------------------------------------------------------------------------
        2414             :
        2415             :     call t_startf("scops")
        2416             :     if (Ncolumns .gt. 1) then
        2417             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2418             :        ! Generate subcolumns for clouds (SCOPS) and precipitation type (PREC_SCOPS)
        2419             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2420             :        ! RNG used for subcolumn generation
        2421             :        allocate(rngs(nPoints), seed(nPoints), stat=istat)
        2422             :        call handle_allocate_error(istat, sub, 'rngs, seed')
        2423             :        seed = int(sfcP)
        2424             :        if (Npoints .gt. 1) seed=(sfcP-int(sfcP))*1000000
        2425             :        call init_rng(rngs, seed)
        2426             :
        2427             :        ! Call scops
        2428             :        call scops(NPoints,Nlevels,Ncolumns,rngs,tca,cca,overlap,cospIN%frac_out,0)
        2429             :        deallocate(seed,rngs)
        2430             :
        2431             :        ! Sum up precipitation rates.
        2432             :        allocate(ls_p_rate(nPoints,nLevels), cv_p_rate(nPoints,Nlevels), stat=istat)
        2433             :        call handle_allocate_error(istat, sub, 'ls_p_rate, cv_p_rate')
        2434             :        ls_p_rate(:,1:nLevels) = fl_lsrainIN + fl_lssnowIN + fl_lsgrplIN
        2435             :        cv_p_rate(:,1:nLevels) = fl_ccrainIN + fl_ccsnowIN
        2436             :
        2437             :        ! Call PREC_SCOPS
        2438             :        allocate(frac_prec(nPoints,nColumns,nLevels), stat=istat)
        2439             :        call handle_allocate_error(istat, sub, 'frac_prec')
        2440             :        call prec_scops(nPoints,nLevels,nColumns,ls_p_rate,cv_p_rate,cospIN%frac_out,frac_prec)
        2441             :        deallocate(ls_p_rate,cv_p_rate)
        2442             :
        2443             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2444             :        ! Compute precipitation fraction in each gridbox
        2445             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2446             :        allocate(frac_ls(nPoints,nLevels),prec_ls(nPoints,nLevels), &
        2447             :                 frac_cv(nPoints,nLevels),prec_cv(nPoints,nLevels), stat=istat)
        2448             :        call handle_allocate_error(istat, sub, 'frac_ls,..,prec_cv')
        2449             :
        2450             :        ! Initialize
        2451             :        frac_ls(1:nPoints,1:nLevels) = 0._wp
        2452             :        prec_ls(1:nPoints,1:nLevels) = 0._wp
        2453             :        frac_cv(1:nPoints,1:nLevels) = 0._wp
        2454             :        prec_cv(1:nPoints,1:nLevels) = 0._wp
        2455             :        do j=1,nPoints
        2456             :           do k=1,nLevels
        2457             :              do i=1,nColumns
        2458             :                 if (cospIN%frac_out(j,i,k)  .eq. 1)  frac_ls(j,k) = frac_ls(j,k)+1._wp
        2459             :                 if (cospIN%frac_out(j,i,k)  .eq. 2)  frac_cv(j,k) = frac_cv(j,k)+1._wp
        2460             :                 if (frac_prec(j,i,k) .eq. 1)         prec_ls(j,k) = prec_ls(j,k)+1._wp
        2461             :                 if (frac_prec(j,i,k) .eq. 2)         prec_cv(j,k) = prec_cv(j,k)+1._wp
        2462             :                 if (frac_prec(j,i,k) .eq. 3)         prec_cv(j,k) = prec_cv(j,k)+1._wp
        2463             :                 if (frac_prec(j,i,k) .eq. 3)         prec_ls(j,k) = prec_ls(j,k)+1._wp
        2464             :              enddo
        2465             :              frac_ls(j,k)=frac_ls(j,k)/nColumns
        2466             :              frac_cv(j,k)=frac_cv(j,k)/nColumns
        2467             :              prec_ls(j,k)=prec_ls(j,k)/nColumns
        2468             :              prec_cv(j,k)=prec_cv(j,k)/nColumns
        2469             :
        2470             :              ! Adjust grid-box mean snow properties to local properties
        2471             :              ! Convert longwave optical depth to longwave emissivity
        2472             :              if (prec_ls(j,k) .ne. 0._r8 .and. dtau_s_snow(j,k) .gt. 0._r8) then
        2473             :                 dtau_s_snow(j,k) = dtau_s_snow(j,k)/prec_ls(j,k)
        2474             :              end if
        2475             :              if (prec_ls(j,k) .ne. 0._r8 .and. dem_s_snow(j,k) .gt. 0._r8) then
        2476             :                 dem_s_snow(j,k) = dem_s_snow(j,k)/prec_ls(j,k)
        2477             :                 dem_s_snow(j,k) = 1._r8 - exp ( -1._r8*dem_s_snow(j,k))
        2478             :              end if !!+JEK
        2479             :           enddo
        2480             :        enddo
        2481             :
        2482             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2483             :        ! Compute mixing ratios, effective radii and precipitation fluxes for clouds
        2484             :        ! and precipitation
        2485             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2486             :        allocate(mr_hydro(nPoints,nColumns,nLevels,nHydro), &
        2487             :                 Reff(nPoints,nColumns,nLevels,nHydro),     &
        2488             :                 Np(nPoints,nColumns,nLevels,nHydro), stat=istat)
        2489             :        call handle_allocate_error(istat, sub, 'mr_hydro,Reff,Np')
        2490             :
        2491             :        ! Initialize
        2492             :        mr_hydro(:,:,:,:) = 0._wp
        2493             :        Reff(:,:,:,:)     = 0._wp
        2494             :        Np(:,:,:,:)       = 0._wp
        2495             :
        2496             :        do k=1,nColumns
        2497             :           ! Subcolumn clouds
        2498             :           column_frac_out = cospIN%frac_out(:,k,:)
        2499             :
        2500             :           ! LS clouds
        2501             :           where (column_frac_out == I_LSC)
        2502             :              mr_hydro(:,k,:,I_LSCLIQ) = mr_lsliq
        2503             :              mr_hydro(:,k,:,I_LSCICE) = mr_lsice
        2504             :              Reff(:,k,:,I_LSCLIQ)     = ReffIN(:,:,I_LSCLIQ)
        2505             :              Reff(:,k,:,I_LSCICE)     = ReffIN(:,:,I_LSCICE)
        2506             :           ! CONV clouds
        2507             :           elsewhere (column_frac_out == I_CVC)
        2508             :              mr_hydro(:,k,:,I_CVCLIQ) = mr_ccliq
        2509             :              mr_hydro(:,k,:,I_CVCICE) = mr_ccice
        2510             :              Reff(:,k,:,I_CVCLIQ)     = ReffIN(:,:,I_CVCLIQ)
        2511             :              Reff(:,k,:,I_CVCICE)     = ReffIN(:,:,I_CVCICE)
        2512             :           end where
        2513             :
        2514             :           ! Subcolumn precipitation
        2515             :           column_prec_out = frac_prec(:,k,:)
        2516             :
        2517             :           ! LS Precipitation
        2518             :           where ((column_prec_out == 1) .or. (column_prec_out == 3) )
        2519             :              Reff(:,k,:,I_LSRAIN) = ReffIN(:,:,I_LSRAIN)
        2520             :              Reff(:,k,:,I_LSSNOW) = ReffIN(:,:,I_LSSNOW)
        2521             :              Reff(:,k,:,I_LSGRPL) = ReffIN(:,:,I_LSGRPL)
        2522             :           ! CONV precipitation
        2523             :           elsewhere ((column_prec_out == 2) .or. (column_prec_out == 3))
        2524             :              Reff(:,k,:,I_CVRAIN) = ReffIN(:,:,I_CVRAIN)
        2525             :              Reff(:,k,:,I_CVSNOW) = ReffIN(:,:,I_CVSNOW)
        2526             :           end where
        2527             :        enddo
        2528             :
        2529             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2530             :        ! Convert the mixing ratio and precipitation fluxes from gridbox mean to
        2531             :        ! the fraction-based values
        2532             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2533             :        do k=1,nLevels
        2534             :           do j=1,nPoints
        2535             :              ! Clouds
        2536             :              if (frac_ls(j,k) .ne. 0._r8) then
        2537             :                 mr_hydro(j,:,k,I_LSCLIQ) = mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
        2538             :                 mr_hydro(j,:,k,I_LSCICE) = mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
        2539             :              endif
        2540             :              if (frac_cv(j,k) .ne. 0._r8) then
        2541             :                 mr_hydro(j,:,k,I_CVCLIQ) = mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
        2542             :                 mr_hydro(j,:,k,I_CVCICE) = mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
        2543             :              endif
        2544             :
        2545             :              ! Precipitation
        2546             :              if (prec_ls(j,k) .ne. 0._r8) then
        2547             :                 fl_lsrain(j,k) = fl_lsrainIN(j,k)/prec_ls(j,k)
        2548             :                 fl_lssnow(j,k) = fl_lssnowIN(j,k)/prec_ls(j,k)
        2549             :                 fl_lsgrpl(j,k) = fl_lsgrplIN(j,k)/prec_ls(j,k)
        2550             :              endif
        2551             :              if (prec_cv(j,k) .ne. 0._r8) then
        2552             :                 fl_ccrain(j,k) = fl_ccrainIN(j,k)/prec_cv(j,k)
        2553             :                 fl_ccsnow(j,k) = fl_ccsnowIN(j,k)/prec_cv(j,k)
        2554             :              endif
        2555             :           enddo
        2556             :        enddo
        2557             :
        2558             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2559             :        ! Convert precipitation fluxes to mixing ratios
        2560             :        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2561             :
        2562             :        ! LS rain
        2563             :        call cosp_precip_mxratio(nPoints, nLevels, nColumns, cospstateIN%pfull,        &
        2564             :             cospstateIN%at, frac_prec, 1._wp, n_ax(I_LSRAIN), n_bx(I_LSRAIN),         &
        2565             :             alpha_x(I_LSRAIN), c_x(I_LSRAIN),   d_x(I_LSRAIN),   g_x(I_LSRAIN),       &
        2566             :             a_x(I_LSRAIN),   b_x(I_LSRAIN),   gamma_1(I_LSRAIN), gamma_2(I_LSRAIN),   &
        2567             :             gamma_3(I_LSRAIN), gamma_4(I_LSRAIN), fl_lsrain,                          &
        2568             :             mr_hydro(:,:,:,I_LSRAIN), Reff(:,:,:,I_LSRAIN))
        2569             :        ! LS snow
        2570             :        call cosp_precip_mxratio(nPoints, nLevels, nColumns, cospstateIN%pfull,        &
        2571             :             cospstateIN%at, frac_prec, 1._wp,  n_ax(I_LSSNOW),  n_bx(I_LSSNOW),       &
        2572             :             alpha_x(I_LSSNOW), c_x(I_LSSNOW),  d_x(I_LSSNOW),  g_x(I_LSSNOW),         &
        2573             :             a_x(I_LSSNOW),   b_x(I_LSSNOW),   gamma_1(I_LSSNOW),  gamma_2(I_LSSNOW),  &
        2574             :             gamma_3(I_LSSNOW), gamma_4(I_LSSNOW), fl_lssnow,                          &
        2575             :             mr_hydro(:,:,:,I_LSSNOW), Reff(:,:,:,I_LSSNOW))
        2576             :        ! CV rain
        2577             :        call cosp_precip_mxratio(nPoints, nLevels, nColumns, cospstateIN%pfull,        &
        2578             :             cospstateIN%at, frac_prec, 2._wp, n_ax(I_CVRAIN),  n_bx(I_CVRAIN),        &
        2579             :             alpha_x(I_CVRAIN), c_x(I_CVRAIN),   d_x(I_CVRAIN),   g_x(I_CVRAIN),       &
        2580             :             a_x(I_CVRAIN),   b_x(I_CVRAIN),   gamma_1(I_CVRAIN), gamma_2(I_CVRAIN),   &
        2581             :             gamma_3(I_CVRAIN), gamma_4(I_CVRAIN), fl_ccrain,                          &
        2582             :             mr_hydro(:,:,:,I_CVRAIN), Reff(:,:,:,I_CVRAIN))
        2583             :        ! CV snow
        2584             :        call cosp_precip_mxratio(nPoints, nLevels, nColumns, cospstateIN%pfull,        &
        2585             :             cospstateIN%at, frac_prec, 2._wp, n_ax(I_CVSNOW),  n_bx(I_CVSNOW),        &
        2586             :             alpha_x(I_CVSNOW),  c_x(I_CVSNOW),   d_x(I_CVSNOW),   g_x(I_CVSNOW),      &
        2587             :             a_x(I_CVSNOW),   b_x(I_CVSNOW),   gamma_1(I_CVSNOW), gamma_2(I_CVSNOW),   &
        2588             :             gamma_3(I_CVSNOW), gamma_4(I_CVSNOW), fl_ccsnow,                          &
        2589             :             mr_hydro(:,:,:,I_CVSNOW), Reff(:,:,:,I_CVSNOW))
        2590             :        ! LS groupel.
        2591             :        call cosp_precip_mxratio(nPoints, nLevels, nColumns, cospstateIN%pfull,        &
        2592             :             cospstateIN%at, frac_prec, 1._wp, n_ax(I_LSGRPL),  n_bx(I_LSGRPL),        &
        2593             :             alpha_x(I_LSGRPL), c_x(I_LSGRPL),   d_x(I_LSGRPL),   g_x(I_LSGRPL),       &
        2594             :             a_x(I_LSGRPL),   b_x(I_LSGRPL),   gamma_1(I_LSGRPL),  gamma_2(I_LSGRPL),  &
        2595             :             gamma_3(I_LSGRPL), gamma_4(I_LSGRPL), fl_lsgrpl,                          &
        2596             :             mr_hydro(:,:,:,I_LSGRPL), Reff(:,:,:,I_LSGRPL))
        2597             :
        2598             :     else
        2599             :        cospIN%frac_out(:,:,:) = 1
        2600             :        allocate(mr_hydro(nPoints, 1,nLevels,nHydro),Reff(nPoints,1,nLevels,nHydro),      &
        2601             :                 Np(nPoints,1,nLevels,nHydro), stat=istat)
        2602             :        call handle_allocate_error(istat, sub, 'mr_hydro,Reff,Np')
        2603             :        mr_hydro(:,1,:,I_LSCLIQ) = mr_lsliq
        2604             :        mr_hydro(:,1,:,I_LSCICE) = mr_lsice
        2605             :        mr_hydro(:,1,:,I_CVCLIQ) = mr_ccliq
        2606             :        mr_hydro(:,1,:,I_CVCICE) = mr_ccice
        2607             :        Reff(:,1,:,:)            = ReffIN
        2608             :     endif
        2609             :     call t_stopf("scops")
        2610             :
        2611             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2612             :     ! CLOUDSAT RADAR OPTICS
        2613             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2614             :     call t_startf("cloudsat_optics")
        2615             :     if (lradar_sim) then
        2616             :        ! Compute gaseous absorption (assume identical for each subcolun)
        2617             :        allocate(g_vol(nPoints,nLevels), stat=istat)
        2618             :        call handle_allocate_error(istat, sub, 'g_vol')
        2619             :        g_vol(:,:)=0._wp
        2620             :        do i = 1, nPoints
        2621             :           do j = 1, nLevels
        2622             :              if (cospIN%rcfg_cloudsat%use_gas_abs == 1 .or. &
        2623             :                 (cospIN%rcfg_cloudsat%use_gas_abs == 2 .and. j == 1)) then
        2624             :                 g_vol(i,j) = gases(cospstateIN%pfull(i,j), cospstateIN%at(i,j),    &
        2625             :                                    cospstateIN%qv(i,j), cospIN%rcfg_cloudsat%freq)
        2626             :              endif
        2627             :              cospIN%g_vol_cloudsat(i,:,j) = g_vol(i,j)
        2628             :           end do
        2629             :        end do
        2630             :
        2631             :        ! Loop over all subcolumns
        2632             :        allocate(fracPrecipIce(nPoints,nColumns,nLevels), stat=istat)
        2633             :        call handle_allocate_error(istat, sub, 'fracPrecipIce')
        2634             :        fracPrecipIce(:,:,:) = 0._wp
        2635             :        do k=1,nColumns
        2636             :           call quickbeam_optics(sd, cospIN%rcfg_cloudsat, nPoints, nLevels, R_UNDEF, &
        2637             :                mr_hydro(:,k,:,1:nHydro)*1000._wp, Reff(:,k,:,1:nHydro)*1.e6_wp,      &
        2638             :                Np(:,k,:,1:nHydro), cospstateIN%pfull, cospstateIN%at,                &
        2639             :                cospstateIN%qv, cospIN%z_vol_cloudsat(1:nPoints,k,:),                 &
        2640             :                cospIN%kr_vol_cloudsat(1:nPoints,k,:))
        2641             :
        2642             :         ! At each model level, what fraction of the precipitation is frozen?
        2643             :           where(mr_hydro(:,k,:,I_LSRAIN) .gt. 0 .or. mr_hydro(:,k,:,I_LSSNOW) .gt. 0 .or. &
        2644             :                 mr_hydro(:,k,:,I_CVRAIN) .gt. 0 .or. mr_hydro(:,k,:,I_CVSNOW) .gt. 0 .or. &
        2645             :                 mr_hydro(:,k,:,I_LSGRPL) .gt. 0)
        2646             :              fracPrecipIce(:,k,:) = (mr_hydro(:,k,:,I_LSSNOW) + mr_hydro(:,k,:,I_CVSNOW) + &
        2647             :                   mr_hydro(:,k,:,I_LSGRPL)) / &
        2648             :                   (mr_hydro(:,k,:,I_LSSNOW) + mr_hydro(:,k,:,I_CVSNOW) + mr_hydro(:,k,:,I_LSGRPL) + &
        2649             :                   mr_hydro(:,k,:,I_LSRAIN)  + mr_hydro(:,k,:,I_CVRAIN))
        2650             :           elsewhere
        2651             :              fracPrecipIce(:,k,:) = 0._wp
        2652             :           endwhere
        2653             :        enddo
        2654             :
        2655             :        ! Regrid frozen fraction to Cloudsat/Calipso statistical grid
        2656             :        allocate(fracPrecipIce_statGrid(nPoints,nColumns,Nlvgrid), stat=istat)
        2657             :        call handle_allocate_error(istat, sub, 'fracPrecipIce_statGrid')
        2658             :        fracPrecipIce_statGrid(:,:,:) = 0._wp
        2659             :        call cosp_change_vertical_grid(Npoints, Ncolumns, Nlevels, cospstateIN%hgt_matrix(:,Nlevels:1:-1), &
        2660             :             cospstateIN%hgt_matrix_half(:,Nlevels:1:-1), fracPrecipIce(:,:,Nlevels:1:-1), Nlvgrid,  &
        2661             :             vgrid_zl(Nlvgrid:1:-1),  vgrid_zu(Nlvgrid:1:-1), fracPrecipIce_statGrid(:,:,Nlvgrid:1:-1))
        2662             :
        2663             :        ! For near-surface diagnostics, we only need the frozen fraction at one layer.
        2664             :        cospIN%fracPrecipIce(:,:) = fracPrecipIce_statGrid(:,:,cloudsat_preclvl)
        2665             :
        2666             :     endif
        2667             :     call t_stopf("cloudsat_optics")
        2668             :
        2669             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2670             :     ! CALIPSO Polarized optics
        2671             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2672             :     call t_startf("calipso_optics")
        2673             :     if (Llidar_sim) then
        2674             :        ReffTemp = ReffIN
        2675             :        call lidar_optics(nPoints,nColumns,nLevels,5,lidar_ice_type,                      &
        2676             :                          mr_hydro(1:nPoints,1:nColumns,1:nLevels,I_LSCLIQ),              &
        2677             :                          mr_hydro(1:nPoints,1:nColumns,1:nLevels,I_LSCICE),              &
        2678             :                          mr_hydro(1:nPoints,1:nColumns,1:nLevels,I_CVCLIQ),              &
        2679             :                          mr_hydro(1:nPoints,1:nColumns,1:nLevels,I_CVCICE),              &
        2680             :                          mr_hydro(1:nPoints,1:nColumns,1:nLevels,I_LSSNOW),              &
        2681             :                          ReffTemp(1:nPoints,1:nLevels,I_LSCLIQ),                         &
        2682             :                          ReffTemp(1:nPoints,1:nLevels,I_LSCICE),                         &
        2683             :                          ReffTemp(1:nPoints,1:nLevels,I_CVCLIQ),                         &
        2684             :                          ReffTemp(1:nPoints,1:nLevels,I_CVCICE),                         &
        2685             :                          ReffTemp(1:nPoints,1:nLevels,I_LSSNOW),                         &
        2686             :                          cospstateIN%pfull(1:nPoints,1:nLevels),                         &
        2687             :                          cospstateIN%phalf(1:nPoints,1:nLevels+1),                       &
        2688             :                          cospstateIN%at(1:nPoints,1:nLevels),                            &
        2689             :                          cospIN%beta_mol_calipso(1:nPoints,1:nLevels),                   &
        2690             :                          cospIN%betatot_calipso(1:nPoints,1:nColumns,1:nLevels),         &
        2691             :                          cospIN%tau_mol_calipso(1:nPoints,1:nLevels),                    &
        2692             :                          cospIN%tautot_calipso(1:nPoints,1:nColumns,1:nLevels),          &
        2693             :                          cospIN%tautot_S_liq(1:nPoints,1:nColumns),                      &
        2694             :                          cospIN%tautot_S_ice(1:nPoints,1:nColumns),                      &
        2695             :                          cospIN%betatot_ice_calipso(1:nPoints,1:nColumns,1:nLevels),     &
        2696             :                          cospIN%betatot_liq_calipso(1:nPoints,1:nColumns,1:nLevels),     &
        2697             :                          cospIN%tautot_ice_calipso(1:nPoints,1:nColumns,1:nLevels),      &
        2698             :                          cospIN%tautot_liq_calipso(1:nPoints,1:nColumns,1:nLevels))
        2699             :     endif
        2700             :     call t_stopf("calipso_optics")
        2701             :
        2702             :
        2703             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2704             :     ! Compute optical fields for passive simulators (i.e. only sunlit points)
        2705             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2706             :
        2707             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2708             :     ! 11 micron emissivity (needed by the ISCCP simulator)
        2709             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2710             :     call t_startf("11micron_emissivity")
        2711             :     if (Lisccp_sim) then
        2712             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,dem_c,dem_s,  &
        2713             :             cospIN%emiss_11)
        2714             :        ! Add in contributions from radiative snow
        2715             :        do j=1,nColumns
        2716             :           where(frac_prec(:,j,:) .eq. 1 .or. frac_prec(:,j,:) .eq. 3)
        2717             :              cospIN%emiss_11(:,j,:) = 1._wp - (1- cospIN%emiss_11(:,j,:))*(1-dem_s_snow)
        2718             :           endwhere
        2719             :        enddo
        2720             :     endif
        2721             :     call t_stopf("11micron_emissivity")
        2722             :
        2723             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2724             :     ! 0.67 micron optical depth (needed by ISCCP, MISR and MODIS simulators)
        2725             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2726             :     call t_startf("067tau")
        2727             :     if (Lisccp_sim .or. Lmisr_sim .or. Lmodis_sim) then
        2728             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,dtau_c,dtau_s,&
        2729             :             cospIN%tau_067)
        2730             :
        2731             :        ! Add in contributions from snow
        2732             :        do j=1,nColumns
        2733             :           where((frac_prec(:,j,:) .eq. 1 .or. frac_prec(:,j,:) .eq. 3) .and. &
        2734             :                Reff(:,j,:,I_LSSNOW) .gt. 0._r8 .and. dtau_s_snow .gt. 0._r8)
        2735             :              cospIN%tau_067(:,j,:)  = cospIN%tau_067(:,j,:)+dtau_s_snow
        2736             :           endwhere
        2737             :        enddo
        2738             :     endif
        2739             :     call t_stopf("067tau")
        2740             :
        2741             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2742             :     ! MODIS optics
        2743             :     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2744             :     call t_startf("modis_optics")
        2745             :     if (lmodis_sim) then
        2746             :        allocate(MODIS_cloudWater(nPoints,nColumns,nLevels),                              &
        2747             :                 MODIS_cloudIce(nPoints,nColumns,nLevels),                                &
        2748             :                 MODIS_cloudSnow(nPoints,nColumns,nLevels),                               &
        2749             :                 MODIS_waterSize(nPoints,nColumns,nLevels),                               &
        2750             :                 MODIS_iceSize(nPoints,nColumns,nLevels),                                 &
        2751             :                 MODIS_snowSize(nPoints,nColumns,nLevels),                                &
        2752             :                 MODIS_opticalThicknessLiq(nPoints,nColumns,nLevels),                     &
        2753             :                 MODIS_opticalThicknessIce(nPoints,nColumns,nLevels),                     &
        2754             :                 MODIS_opticalThicknessSnow(nPoints,nColumns,nLevels), stat=istat)
        2755             :        call handle_allocate_error(istat, sub, 'MODIS_*')
        2756             :
        2757             :        ! Cloud water
        2758             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,              &
        2759             :             mr_hydro(:,:,:,I_CVCLIQ),mr_hydro(:,:,:,I_LSCLIQ),MODIS_cloudWater)
        2760             :        ! Cloud ice
        2761             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,              &
        2762             :             mr_hydro(:,:,:,I_CVCICE),mr_hydro(:,:,:,I_LSCICE),MODIS_cloudIce)
        2763             :        ! Cloud water droplet size
        2764             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,              &
        2765             :             Reff(:,:,:,I_CVCLIQ),Reff(:,:,:,I_LSCLIQ),MODIS_waterSize)
        2766             :        ! Cloud ice crystal size
        2767             :        call cosp_simulator_optics(nPoints,nColumns,nLevels,cospIN%frac_out,              &
        2768             :             Reff(:,:,:,I_CVCICE),Reff(:,:,:,I_LSCICE),MODIS_iceSize)
        2769             :
        2770             :        ! Cloud snow and size
        2771             :        MODIS_snowSize(:,:,:)  = Reff(:,:,:,I_LSSNOW)
        2772             :        do j=1,nColumns
        2773             :           where((frac_prec(:,j,:) .eq. 1 .or. frac_prec(:,j,:) .eq. 3) .and. &
        2774             :                Reff(:,j,:,I_LSSNOW) .gt. 0._r8 .and. dtau_s_snow .gt. 0._r8)
        2775             :              MODIS_cloudSnow(:,j,:) = mr_hydro(:,j,:,I_LSSNOW)
        2776             :              MODIS_snowSize(:,j,:)  = Reff(:,j,:,I_LSSNOW)
        2777             :           elsewhere
        2778             :              MODIS_snowSize(:,j,:)  = 0._wp
        2779             :              MODIS_cloudSnow(:,j,:) = 0._wp
        2780             :           endwhere
        2781             :        enddo
        2782             :
        2783             :        ! Partition optical thickness into liquid and ice parts
        2784             :        call modis_optics_partition(nPoints, nLevels, nColumns, MODIS_cloudWater,     &
        2785             :             MODIS_cloudIce, MODIS_cloudSnow, MODIS_waterSize, MODIS_iceSize,         &
        2786             :             MODIS_snowSize, cospIN%tau_067, MODIS_opticalThicknessLiq,               &
        2787             :             MODIS_opticalThicknessIce, MODIS_opticalThicknessSnow)
        2788             :
        2789             :        ! Compute asymmetry parameter and single scattering albedo
        2790             :        call modis_optics(nPoints, nLevels, nColumns, MODIS_opticalThicknessLiq,      &
        2791             :             MODIS_waterSize*1.0e6_wp, MODIS_opticalThicknessIce,                     &
        2792             :             MODIS_iceSize*1.0e6_wp, MODIS_opticalThicknessSnow,                      &
        2793             :             MODIS_snowSize*1.0e6_wp, cospIN%fracLiq, cospIN%asym, cospIN%ss_alb)
        2794             :
        2795             :     endif ! MODIS simulator optics
        2796             :     call t_stopf("modis_optics")
        2797             :
        2798             :   end subroutine subsample_and_optics
        2799             :
        2800             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2801             :   ! SUBROUTINE construct_cospIN
        2802             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2803             :   subroutine construct_cospIN(npoints,ncolumns,nlevels,y)
        2804             :     ! Inputs
        2805             :     integer,intent(in) :: &
        2806             :          npoints,  & ! Number of horizontal gridpoints
        2807             :          ncolumns, & ! Number of subcolumns
        2808             :          nlevels     ! Number of vertical levels
        2809             :     ! Outputs
        2810             :     type(cosp_optical_inputs),intent(out) :: y
        2811             :
        2812             :     ! local
        2813             :     integer :: istat
        2814             :     character(len=*), parameter :: sub = 'construct_cospIN'
        2815             :     !--------------------------------------------------------------------------------------
        2816             :
        2817             :     ! Dimensions
        2818             :     y%Npoints  = Npoints
        2819             :     y%Ncolumns = Ncolumns
        2820             :     y%Nlevels  = Nlevels
        2821             :     y%Npart    = 4
        2822             :     y%Nrefl    = PARASOL_NREFL
        2823             :
        2824             :     allocate(y%tau_067(            npoints, ncolumns, nlevels),&
        2825             :              y%emiss_11(           npoints, ncolumns, nlevels),&
        2826             :              y%frac_out(           npoints, ncolumns, nlevels),&
        2827             :              y%betatot_calipso(    npoints, ncolumns, nlevels),&
        2828             :              y%betatot_ice_calipso(npoints, ncolumns, nlevels),&
        2829             :              y%fracLiq(            npoints, ncolumns, nlevels),&
        2830             :              y%betatot_liq_calipso(npoints, ncolumns, nlevels),&
        2831             :              y%tautot_calipso(     npoints, ncolumns, nlevels),&
        2832             :              y%tautot_ice_calipso( npoints, ncolumns, nlevels),&
        2833             :              y%tautot_liq_calipso( npoints, ncolumns, nlevels),&
        2834             :              y%z_vol_cloudsat(     npoints, ncolumns, nlevels),&
        2835             :              y%kr_vol_cloudsat(    npoints, ncolumns, nlevels),&
        2836             :              y%g_vol_cloudsat(     npoints, ncolumns, nlevels),&
        2837             :              y%asym(               npoints, ncolumns, nlevels),&
        2838             :              y%ss_alb(             npoints, ncolumns, nlevels),&
        2839             :              y%beta_mol_calipso(   npoints,           nlevels),&
        2840             :              y%tau_mol_calipso(    npoints,           nlevels),&
        2841             :              y%tautot_S_ice(       npoints, ncolumns         ),&
        2842             :              y%tautot_S_liq(       npoints, ncolumns)         ,&
        2843             :              y%fracPrecipIce(npoints,   ncolumns), stat=istat)
        2844             :     call handle_allocate_error(istat, sub, 'tau_067,..,fracPrecipIce')
        2845             :
        2846             :   end subroutine construct_cospIN
        2847             :
        2848             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2849             :   ! SUBROUTINE construct_cospstateIN
        2850             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        2851             :   subroutine construct_cospstateIN(npoints,nlevels,nchan,y)
        2852             :     ! Inputs
        2853             :     integer,intent(in) :: &
        2854             :          npoints, & ! Number of horizontal gridpoints
        2855             :          nlevels, & ! Number of vertical levels
        2856             :          nchan      ! Number of channels
        2857             :     ! Outputs
        2858             :     type(cosp_column_inputs),intent(out) :: y
        2859             :
        2860             :     ! local
        2861             :     integer :: istat
        2862             :     character(len=*), parameter :: sub = 'construct_cospstateIN'
        2863             :     !--------------------------------------------------------------------------------------
        2864             :
        2865             :     allocate( &
        2866             :        y%sunlit(npoints), &
        2867             :        y%at(npoints,nlevels), &
        2868             :        y%pfull(npoints,nlevels), &
        2869             :        y%phalf(npoints,nlevels+1), &
        2870             :        y%qv(npoints,nlevels), &
        2871             :        y%hgt_matrix(npoints,nlevels), &
        2872             :        y%hgt_matrix_half(npoints,nlevels+1), &
        2873             :        y%land(npoints), &
        2874             :        y%skt(npoints), &
        2875             :        y%surfelev(nPoints), &
        2876             :        y%emis_sfc(nchan), &
        2877             :        y%u_sfc(npoints), &
        2878             :        y%v_sfc(npoints), &
        2879             :        y%seaice(npoints), &
        2880             :        y%lat(npoints), &
        2881             :        y%lon(nPoints), &
        2882             :        y%o3(npoints,nlevels), &
        2883             :        y%tca(nPoints,nLevels), &
        2884             :        y%cloudIce(nPoints,nLevels), &
        2885             :        y%cloudLiq(nPoints,nLevels), &
        2886             :        y%fl_rain(nPoints,nLevels), &
        2887             :        y%fl_snow(nPoints,nLevels), stat=istat)
        2888             :     call handle_allocate_error(istat, sub, 'sunlit,..,fl_snow')
        2889             :
        2890             :   end subroutine construct_cospstateIN
        2891             :   ! ######################################################################################
        2892             :   ! SUBROUTINE construct_cosp_outputs
        2893             :   !
        2894             :   ! This subroutine allocates output fields based on input logical flag switches.
        2895             :   ! ######################################################################################
        2896             :   subroutine construct_cosp_outputs(Npoints,Ncolumns,Nlevels,Nlvgrid,x)
        2897             :     ! Inputs
        2898             :     integer,intent(in) :: &
        2899             :          Npoints,         & ! Number of sampled points
        2900             :          Ncolumns,        & ! Number of subgrid columns
        2901             :          Nlevels,         & ! Number of model levels
        2902             :          Nlvgrid            ! Number of levels in L3 stats computation
        2903             :
        2904             :     ! Outputs
        2905             :     type(cosp_outputs),intent(out) :: &
        2906             :          x           ! COSP output structure
        2907             :
        2908             :     ! local
        2909             :     integer :: istat
        2910             :     character(len=*), parameter :: sub = 'construct_cosp_outputs'
        2911             :     !--------------------------------------------------------------------------------------
        2912             :
        2913             :      ! ISCCP simulator outputs
        2914             :     if (lisccp_sim) then
        2915             :        allocate( &
        2916             :           x%isccp_boxtau(Npoints,Ncolumns), &
        2917             :           x%isccp_boxptop(Npoints,Ncolumns), &
        2918             :           x%isccp_fq(Npoints,numISCCPTauBins,numISCCPPresBins), &
        2919             :           x%isccp_totalcldarea(Npoints), &
        2920             :           x%isccp_meanptop(Npoints), &
        2921             :           x%isccp_meantaucld(Npoints), &
        2922             :           x%isccp_meantb(Npoints), &
        2923             :           x%isccp_meantbclr(Npoints), &
        2924             :           x%isccp_meanalbedocld(Npoints), stat=istat)
        2925             :        call handle_allocate_error(istat, sub, 'isccp_*')
        2926             :     endif
        2927             :
        2928             :     ! MISR simulator
        2929             :     if (lmisr_sim) then
        2930             :        allocate( &
        2931             :           x%misr_fq(Npoints,numMISRTauBins,numMISRHgtBins), &
        2932             :           ! *NOTE* These 3 fields are not output, but were part of the v1.4.0 cosp_misr, so
        2933             :           !        they are still computed. Should probably have a logical to control these
        2934             :           !        outputs.
        2935             :           x%misr_dist_model_layertops(Npoints,numMISRHgtBins), &
        2936             :           x%misr_meanztop(Npoints), &
        2937             :           x%misr_cldarea(Npoints), stat=istat)
        2938             :        call handle_allocate_error(istat, sub, 'misr_*')
        2939             :     endif
        2940             :
        2941             :     ! MODIS simulator
        2942             :     if (lmodis_sim) then
        2943             :        allocate( &
        2944             :           x%modis_Cloud_Fraction_Total_Mean(Npoints), &
        2945             :           x%modis_Cloud_Fraction_Water_Mean(Npoints), &
        2946             :           x%modis_Cloud_Fraction_Ice_Mean(Npoints), &
        2947             :           x%modis_Cloud_Fraction_High_Mean(Npoints), &
        2948             :           x%modis_Cloud_Fraction_Mid_Mean(Npoints), &
        2949             :           x%modis_Cloud_Fraction_Low_Mean(Npoints), &
        2950             :           x%modis_Optical_Thickness_Total_Mean(Npoints), &
        2951             :           x%modis_Optical_Thickness_Water_Mean(Npoints), &
        2952             :           x%modis_Optical_Thickness_Ice_Mean(Npoints), &
        2953             :           x%modis_Optical_Thickness_Total_LogMean(Npoints), &
        2954             :           x%modis_Optical_Thickness_Water_LogMean(Npoints), &
        2955             :           x%modis_Optical_Thickness_Ice_LogMean(Npoints), &
        2956             :           x%modis_Cloud_Particle_Size_Water_Mean(Npoints), &
        2957             :           x%modis_Cloud_Particle_Size_Ice_Mean(Npoints), &
        2958             :           x%modis_Cloud_Top_Pressure_Total_Mean(Npoints), &
        2959             :           x%modis_Liquid_Water_Path_Mean(Npoints), &
        2960             :           x%modis_Ice_Water_Path_Mean(Npoints), &
        2961             :           x%modis_Optical_Thickness_vs_Cloud_Top_Pressure(nPoints,numModisTauBins,numMODISPresBins), &
        2962             :           x%modis_Optical_thickness_vs_ReffLIQ(nPoints,numMODISTauBins,numMODISReffLiqBins), &
        2963             :           x%modis_Optical_Thickness_vs_ReffICE(nPoints,numMODISTauBins,numMODISReffIceBins), &
        2964             :           stat=istat)
        2965             :        call handle_allocate_error(istat, sub, 'modis_*')
        2966             :     endif
        2967             :
        2968             :     ! CALIPSO simulator
        2969             :     if (llidar_sim) then
        2970             :        allocate( &
        2971             :           x%calipso_beta_mol(Npoints,Nlevels), &
        2972             :           x%calipso_beta_tot(Npoints,Ncolumns,Nlevels), &
        2973             :           x%calipso_srbval(SR_BINS+1), &
        2974             :           x%calipso_cfad_sr(Npoints,SR_BINS,Nlvgrid), &
        2975             :           x%calipso_betaperp_tot(Npoints,Ncolumns,Nlevels), &
        2976             :           x%calipso_lidarcld(Npoints,Nlvgrid), &
        2977             :           x%calipso_cldlayer(Npoints,LIDAR_NCAT), &
        2978             :           x%calipso_lidarcldphase(Npoints,Nlvgrid,6), &
        2979             :           x%calipso_lidarcldtmp(Npoints,LIDAR_NTEMP,5), &
        2980             :           x%calipso_cldlayerphase(Npoints,LIDAR_NCAT,6), &
        2981             :           x%calipso_tau_tot(Npoints,Ncolumns,Nlevels), &
        2982             :           x%calipso_temp_tot(Npoints,Nlevels), stat=istat)
        2983             :        call handle_allocate_error(istat, sub, 'calipso_*')
        2984             :     endif
        2985             :
        2986             :     ! PARASOL
        2987             :     if (lparasol_sim) then
        2988             :        allocate( &
        2989             :           x%parasolPix_refl(Npoints,Ncolumns,PARASOL_NREFL), &
        2990             :           x%parasolGrid_refl(Npoints,PARASOL_NREFL), stat=istat)
        2991             :        call handle_allocate_error(istat, sub, 'parasol*')
        2992             :     endif
        2993             :
        2994             :     ! Cloudsat simulator
        2995             :     if (lradar_sim) then
        2996             :        allocate( &
        2997             :           x%cloudsat_Ze_tot(Npoints,Ncolumns,Nlevels), &
        2998             :           x%cloudsat_cfad_ze(Npoints,CLOUDSAT_DBZE_BINS,Nlvgrid), &
        2999             :           x%lidar_only_freq_cloud(Npoints,Nlvgrid), &
        3000             :           x%radar_lidar_tcc(Npoints), &
        3001             :           x%cloudsat_precip_cover(Npoints,nCloudsatPrecipClass), &
        3002             :           x%cloudsat_pia(Npoints), stat=istat)
        3003             :        call handle_allocate_error(istat, sub, 'cloudsat*')
        3004             :     endif
        3005             :
        3006             :   end subroutine construct_cosp_outputs
        3007             :
        3008             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3009             :   ! SUBROUTINE destroy_cospIN
        3010             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3011             :   subroutine destroy_cospIN(y)
        3012             :     type(cosp_optical_inputs),intent(inout) :: y
        3013             :
        3014             :     if (allocated(y%tau_067))             deallocate(y%tau_067)
        3015             :     if (allocated(y%emiss_11))            deallocate(y%emiss_11)
        3016             :     if (allocated(y%frac_out))            deallocate(y%frac_out)
        3017             :     if (allocated(y%beta_mol_calipso))    deallocate(y%beta_mol_calipso)
        3018             :     if (allocated(y%tau_mol_calipso))     deallocate(y%tau_mol_calipso)
        3019             :     if (allocated(y%betatot_calipso))     deallocate(y%betatot_calipso)
        3020             :     if (allocated(y%betatot_ice_calipso)) deallocate(y%betatot_ice_calipso)
        3021             :     if (allocated(y%betatot_liq_calipso)) deallocate(y%betatot_liq_calipso)
        3022             :     if (allocated(y%tautot_calipso))      deallocate(y%tautot_calipso)
        3023             :     if (allocated(y%tautot_ice_calipso))  deallocate(y%tautot_ice_calipso)
        3024             :     if (allocated(y%tautot_liq_calipso))  deallocate(y%tautot_liq_calipso)
        3025             :     if (allocated(y%tautot_S_liq))        deallocate(y%tautot_S_liq)
        3026             :     if (allocated(y%tautot_S_ice))        deallocate(y%tautot_S_ice)
        3027             :     if (allocated(y%z_vol_cloudsat))      deallocate(y%z_vol_cloudsat)
        3028             :     if (allocated(y%kr_vol_cloudsat))     deallocate(y%kr_vol_cloudsat)
        3029             :     if (allocated(y%g_vol_cloudsat))      deallocate(y%g_vol_cloudsat)
        3030             :     if (allocated(y%asym))                deallocate(y%asym)
        3031             :     if (allocated(y%ss_alb))              deallocate(y%ss_alb)
        3032             :     if (allocated(y%fracLiq))             deallocate(y%fracLiq)
        3033             :     if (allocated(y%fracPrecipIce))       deallocate(y%fracPrecipIce)
        3034             :   end subroutine destroy_cospIN
        3035             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3036             :   ! SUBROUTINE destroy_cospstateIN
        3037             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3038             :   subroutine destroy_cospstateIN(y)
        3039             :     type(cosp_column_inputs),intent(inout) :: y
        3040             :
        3041             :     if (allocated(y%surfelev))        deallocate(y%surfelev)
        3042             :     if (allocated(y%sunlit))          deallocate(y%sunlit)
        3043             :     if (allocated(y%skt))             deallocate(y%skt)
        3044             :     if (allocated(y%land))            deallocate(y%land)
        3045             :     if (allocated(y%at))              deallocate(y%at)
        3046             :     if (allocated(y%pfull))           deallocate(y%pfull)
        3047             :     if (allocated(y%phalf))           deallocate(y%phalf)
        3048             :     if (allocated(y%qv))              deallocate(y%qv)
        3049             :     if (allocated(y%o3))              deallocate(y%o3)
        3050             :     if (allocated(y%hgt_matrix))      deallocate(y%hgt_matrix)
        3051             :     if (allocated(y%u_sfc))           deallocate(y%u_sfc)
        3052             :     if (allocated(y%v_sfc))           deallocate(y%v_sfc)
        3053             :     if (allocated(y%lat))             deallocate(y%lat)
        3054             :     if (allocated(y%lon))             deallocate(y%lon)
        3055             :     if (allocated(y%emis_sfc))        deallocate(y%emis_sfc)
        3056             :     if (allocated(y%cloudIce))        deallocate(y%cloudIce)
        3057             :     if (allocated(y%cloudLiq))        deallocate(y%cloudLiq)
        3058             :     if (allocated(y%seaice))          deallocate(y%seaice)
        3059             :     if (allocated(y%fl_rain))         deallocate(y%fl_rain)
        3060             :     if (allocated(y%fl_snow))         deallocate(y%fl_snow)
        3061             :     if (allocated(y%tca))             deallocate(y%tca)
        3062             :     if (allocated(y%hgt_matrix_half)) deallocate(y%hgt_matrix_half)
        3063             :
        3064             :   end subroutine destroy_cospstateIN
        3065             :
        3066             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3067             :   ! SUBROUTINE destroy_cosp_outputs
        3068             :   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        3069             :   subroutine destroy_cosp_outputs(y)
        3070             :      type(cosp_outputs),intent(inout) :: y
        3071             :
        3072             :      ! Deallocate and nullify
        3073             :      if (associated(y%calipso_beta_mol))          then
        3074             :         deallocate(y%calipso_beta_mol)
        3075             :         nullify(y%calipso_beta_mol)
        3076             :      endif
        3077             :      if (associated(y%calipso_temp_tot))          then
        3078             :         deallocate(y%calipso_temp_tot)
        3079             :         nullify(y%calipso_temp_tot)
        3080             :      endif
        3081             :      if (associated(y%calipso_betaperp_tot))      then
        3082             :         deallocate(y%calipso_betaperp_tot)
        3083             :         nullify(y%calipso_betaperp_tot)
        3084             :      endif
        3085             :      if (associated(y%calipso_beta_tot))          then
        3086             :         deallocate(y%calipso_beta_tot)
        3087             :         nullify(y%calipso_beta_tot)
        3088             :      endif
        3089             :      if (associated(y%calipso_tau_tot))           then
        3090             :         deallocate(y%calipso_tau_tot)
        3091             :         nullify(y%calipso_tau_tot)
        3092             :      endif
        3093             :      if (associated(y%calipso_lidarcldphase))     then
        3094             :         deallocate(y%calipso_lidarcldphase)
        3095             :         nullify(y%calipso_lidarcldphase)
        3096             :      endif
        3097             :      if (associated(y%calipso_cldlayerphase))     then
        3098             :         deallocate(y%calipso_cldlayerphase)
        3099             :         nullify(y%calipso_cldlayerphase)
        3100             :      endif
        3101             :      if (associated(y%calipso_lidarcldtmp))       then
        3102             :         deallocate(y%calipso_lidarcldtmp)
        3103             :         nullify(y%calipso_lidarcldtmp)
        3104             :      endif
        3105             :      if (associated(y%calipso_cldlayer))          then
        3106             :         deallocate(y%calipso_cldlayer)
        3107             :         nullify(y%calipso_cldlayer)
        3108             :      endif
        3109             :      if (associated(y%calipso_lidarcld))         then
        3110             :         deallocate(y%calipso_lidarcld)
        3111             :         nullify(y%calipso_lidarcld)
        3112             :      endif
        3113             :      if (associated(y%calipso_srbval))            then
        3114             :         deallocate(y%calipso_srbval)
        3115             :         nullify(y%calipso_srbval)
        3116             :      endif
        3117             :      if (associated(y%calipso_cfad_sr))          then
        3118             :         deallocate(y%calipso_cfad_sr)
        3119             :         nullify(y%calipso_cfad_sr)
        3120             :      endif
        3121             :      if (associated(y%parasolPix_refl))           then
        3122             :         deallocate(y%parasolPix_refl)
        3123             :         nullify(y%parasolPix_refl)
        3124             :      endif
        3125             :      if (associated(y%parasolGrid_refl))          then
        3126             :         deallocate(y%parasolGrid_refl)
        3127             :         nullify(y%parasolGrid_refl)
        3128             :      endif
        3129             :      if (associated(y%cloudsat_Ze_tot))           then
        3130             :         deallocate(y%cloudsat_Ze_tot)
        3131             :         nullify(y%cloudsat_Ze_tot)
        3132             :      endif
        3133             :      if (associated(y%cloudsat_precip_cover)) then
        3134             :         deallocate(y%cloudsat_precip_cover)
        3135             :         nullify(y%cloudsat_precip_cover)
        3136             :      endif
        3137             :      if (associated(y%cloudsat_pia)) then
        3138             :         deallocate(y%cloudsat_pia)
        3139             :         nullify(y%cloudsat_pia)
        3140             :      endif
        3141             :      if (associated(y%cloudsat_cfad_ze))          then
        3142             :         deallocate(y%cloudsat_cfad_ze)
        3143             :         nullify(y%cloudsat_cfad_ze)
        3144             :      endif
        3145             :      if (associated(y%radar_lidar_tcc))           then
        3146             :         deallocate(y%radar_lidar_tcc)
        3147             :         nullify(y%radar_lidar_tcc)
        3148             :      endif
        3149             :      if (associated(y%lidar_only_freq_cloud))     then
        3150             :         deallocate(y%lidar_only_freq_cloud)
        3151             :         nullify(y%lidar_only_freq_cloud)
        3152             :      endif
        3153             :      if (associated(y%isccp_totalcldarea))        then
        3154             :         deallocate(y%isccp_totalcldarea)
        3155             :         nullify(y%isccp_totalcldarea)
        3156             :      endif
        3157             :      if (associated(y%isccp_meantb))              then
        3158             :         deallocate(y%isccp_meantb)
        3159             :         nullify(y%isccp_meantb)
        3160             :      endif
        3161             :      if (associated(y%isccp_meantbclr))           then
        3162             :         deallocate(y%isccp_meantbclr)
        3163             :         nullify(y%isccp_meantbclr)
        3164             :      endif
        3165             :      if (associated(y%isccp_meanptop))            then
        3166             :         deallocate(y%isccp_meanptop)
        3167             :         nullify(y%isccp_meanptop)
        3168             :      endif
        3169             :      if (associated(y%isccp_meantaucld))          then
        3170             :         deallocate(y%isccp_meantaucld)
        3171             :         nullify(y%isccp_meantaucld)
        3172             :      endif
        3173             :      if (associated(y%isccp_meanalbedocld))       then
        3174             :         deallocate(y%isccp_meanalbedocld)
        3175             :         nullify(y%isccp_meanalbedocld)
        3176             :      endif
        3177             :      if (associated(y%isccp_boxtau))              then
        3178             :         deallocate(y%isccp_boxtau)
        3179             :         nullify(y%isccp_boxtau)
        3180             :      endif
        3181             :      if (associated(y%isccp_boxptop))             then
        3182             :         deallocate(y%isccp_boxptop)
        3183             :         nullify(y%isccp_boxptop)
        3184             :      endif
        3185             :      if (associated(y%isccp_fq))                  then
        3186             :         deallocate(y%isccp_fq)
        3187             :         nullify(y%isccp_fq)
        3188             :      endif
        3189             :      if (associated(y%misr_fq))                   then
        3190             :         deallocate(y%misr_fq)
        3191             :         nullify(y%misr_fq)
        3192             :      endif
        3193             :      if (associated(y%misr_dist_model_layertops)) then
        3194             :         deallocate(y%misr_dist_model_layertops)
        3195             :         nullify(y%misr_dist_model_layertops)
        3196             :      endif
        3197             :      if (associated(y%misr_meanztop))             then
        3198             :         deallocate(y%misr_meanztop)
        3199             :         nullify(y%misr_meanztop)
        3200             :      endif
        3201             :      if (associated(y%misr_cldarea))              then
        3202             :         deallocate(y%misr_cldarea)
        3203             :         nullify(y%misr_cldarea)
        3204             :      endif
        3205             :      if (associated(y%rttov_tbs))                 then
        3206             :         deallocate(y%rttov_tbs)
        3207             :         nullify(y%rttov_tbs)
        3208             :      endif
        3209             :      if (associated(y%modis_Cloud_Fraction_Total_Mean))                      then
        3210             :         deallocate(y%modis_Cloud_Fraction_Total_Mean)
        3211             :         nullify(y%modis_Cloud_Fraction_Total_Mean)
        3212             :      endif
        3213             :      if (associated(y%modis_Cloud_Fraction_Ice_Mean))                        then
        3214             :         deallocate(y%modis_Cloud_Fraction_Ice_Mean)
        3215             :         nullify(y%modis_Cloud_Fraction_Ice_Mean)
        3216             :      endif
        3217             :      if (associated(y%modis_Cloud_Fraction_Water_Mean))                      then
        3218             :         deallocate(y%modis_Cloud_Fraction_Water_Mean)
        3219             :         nullify(y%modis_Cloud_Fraction_Water_Mean)
        3220             :      endif
        3221             :      if (associated(y%modis_Cloud_Fraction_High_Mean))                       then
        3222             :         deallocate(y%modis_Cloud_Fraction_High_Mean)
        3223             :         nullify(y%modis_Cloud_Fraction_High_Mean)
        3224             :      endif
        3225             :      if (associated(y%modis_Cloud_Fraction_Mid_Mean))                        then
        3226             :         deallocate(y%modis_Cloud_Fraction_Mid_Mean)
        3227             :         nullify(y%modis_Cloud_Fraction_Mid_Mean)
        3228             :      endif
        3229             :      if (associated(y%modis_Cloud_Fraction_Low_Mean))                        then
        3230             :         deallocate(y%modis_Cloud_Fraction_Low_Mean)
        3231             :         nullify(y%modis_Cloud_Fraction_Low_Mean)
        3232             :      endif
        3233             :      if (associated(y%modis_Optical_Thickness_Total_Mean))                   then
        3234             :         deallocate(y%modis_Optical_Thickness_Total_Mean)
        3235             :         nullify(y%modis_Optical_Thickness_Total_Mean)
        3236             :      endif
        3237             :      if (associated(y%modis_Optical_Thickness_Water_Mean))                   then
        3238             :         deallocate(y%modis_Optical_Thickness_Water_Mean)
        3239             :         nullify(y%modis_Optical_Thickness_Water_Mean)
        3240             :      endif
        3241             :      if (associated(y%modis_Optical_Thickness_Ice_Mean))                     then
        3242             :         deallocate(y%modis_Optical_Thickness_Ice_Mean)
        3243             :         nullify(y%modis_Optical_Thickness_Ice_Mean)
        3244             :      endif
        3245             :      if (associated(y%modis_Optical_Thickness_Total_LogMean))                then
        3246             :         deallocate(y%modis_Optical_Thickness_Total_LogMean)
        3247             :         nullify(y%modis_Optical_Thickness_Total_LogMean)
        3248             :      endif
        3249             :      if (associated(y%modis_Optical_Thickness_Water_LogMean))                then
        3250             :         deallocate(y%modis_Optical_Thickness_Water_LogMean)
        3251             :         nullify(y%modis_Optical_Thickness_Water_LogMean)
        3252             :      endif
        3253             :      if (associated(y%modis_Optical_Thickness_Ice_LogMean))                  then
        3254             :         deallocate(y%modis_Optical_Thickness_Ice_LogMean)
        3255             :         nullify(y%modis_Optical_Thickness_Ice_LogMean)
        3256             :      endif
        3257             :      if (associated(y%modis_Cloud_Particle_Size_Water_Mean))                 then
        3258             :         deallocate(y%modis_Cloud_Particle_Size_Water_Mean)
        3259             :         nullify(y%modis_Cloud_Particle_Size_Water_Mean)
        3260             :      endif
        3261             :      if (associated(y%modis_Cloud_Particle_Size_Ice_Mean))                   then
        3262             :         deallocate(y%modis_Cloud_Particle_Size_Ice_Mean)
        3263             :         nullify(y%modis_Cloud_Particle_Size_Ice_Mean)
        3264             :      endif
        3265             :      if (associated(y%modis_Cloud_Top_Pressure_Total_Mean))                  then
        3266             :         deallocate(y%modis_Cloud_Top_Pressure_Total_Mean)
        3267             :         nullify(y%modis_Cloud_Top_Pressure_Total_Mean)
        3268             :      endif
        3269             :      if (associated(y%modis_Liquid_Water_Path_Mean))                         then
        3270             :         deallocate(y%modis_Liquid_Water_Path_Mean)
        3271             :         nullify(y%modis_Liquid_Water_Path_Mean)
        3272             :      endif
        3273             :      if (associated(y%modis_Ice_Water_Path_Mean))                            then
        3274             :         deallocate(y%modis_Ice_Water_Path_Mean)
        3275             :         nullify(y%modis_Ice_Water_Path_Mean)
        3276             :      endif
        3277             :      if (associated(y%modis_Optical_Thickness_vs_Cloud_Top_Pressure))        then
        3278             :         deallocate(y%modis_Optical_Thickness_vs_Cloud_Top_Pressure)
        3279             :         nullify(y%modis_Optical_Thickness_vs_Cloud_Top_Pressure)
        3280             :      endif
        3281             :      if (associated(y%modis_Optical_thickness_vs_ReffLIQ))                   then
        3282             :         deallocate(y%modis_Optical_thickness_vs_ReffLIQ)
        3283             :         nullify(y%modis_Optical_thickness_vs_ReffLIQ)
        3284             :      endif
        3285             :      if (associated(y%modis_Optical_thickness_vs_ReffICE))                   then
        3286             :         deallocate(y%modis_Optical_thickness_vs_ReffICE)
        3287             :         nullify(y%modis_Optical_thickness_vs_ReffICE)
        3288             :      endif
        3289             :      if (associated(y%calipso_cldtype)) then
        3290             :         deallocate(y%calipso_cldtype)
        3291             :         nullify(y%calipso_cldtype)
        3292             :      endif
        3293             :      if (associated(y%calipso_cldtypetemp)) then
        3294             :         deallocate(y%calipso_cldtypetemp)
        3295             :         nullify(y%calipso_cldtypetemp)
        3296             :      endif
        3297             :      if (associated(y%calipso_cldtypemeanz)) then
        3298             :         deallocate(y%calipso_cldtypemeanz)
        3299             :         nullify(y%calipso_cldtypemeanz)
        3300             :      endif
        3301             :      if (associated(y%calipso_cldtypemeanzse)) then
        3302             :         deallocate(y%calipso_cldtypemeanzse)
        3303             :         nullify(y%calipso_cldtypemeanzse)
        3304             :      endif
        3305             :      if (associated(y%calipso_cldthinemis)) then
        3306             :         deallocate(y%calipso_cldthinemis)
        3307             :         nullify(y%calipso_cldthinemis)
        3308             :      endif
        3309             :      if (associated(y%calipso_lidarcldtype)) then
        3310             :         deallocate(y%calipso_lidarcldtype)
        3311             :         nullify(y%calipso_lidarcldtype)
        3312             :      endif
        3313             :
        3314             :    end subroutine destroy_cosp_outputs
        3315             : #endif
        3316             :
        3317             : !#######################################################################
        3318             : end module cospsimulator_intr

Generated by: LCOV version 1.14