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: 2025-01-13 21:54:50 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