cosp.F90

Go to the documentation of this file.

! (c) British Crown Copyright 2008, the Met Office.
! All rights reserved.
! 
! Redistribution and use in source and binary forms, with or without modification, are permitted 
! provided that the following conditions are met:
! 
!     * Redistributions of source code must retain the above copyright notice, this list 
!       of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright notice, this list
!       of conditions and the following disclaimer in the documentation and/or other materials 
!       provided with the distribution.
!     * Neither the name of the Met Office nor the names of its contributors may be used 
!       to endorse or promote products derived from this software without specific prior written 
!       permission.
! 
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

!!#include "cosp_defs.h"
MODULE mod_cosp
  USE mod_cosp_types
  USE mod_cosp_simulator
  USE mod_phys_lmdz_para
  USE mod_grid_phy_lmdz
  IMPLICIT NONE

CONTAINS


!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!--------------------- SUBROUTINE COSP ---------------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)

  ! Arguments
  integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
  integer,intent(in) :: Ncolumns
  type(cosp_config),intent(in) :: cfg   ! Configuration options
  type(cosp_vgrid),intent(in) :: vgrid   ! Information on vertical grid of stats
  type(cosp_gridbox),intent(inout) :: gbx
  type(cosp_subgrid),intent(inout) :: sgx   ! Subgrid info
  type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
  type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
  type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
  type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
  type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
  type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator

  ! Local variables 
  integer :: Npoints   ! Number of gridpoints
  integer :: Nlevels   ! Number of levels
  integer :: Nhydro    ! Number of hydrometeors
  integer :: Niter     ! Number of calls to cosp_simulator
  integer :: i_first,i_last ! First and last gridbox to be processed in each iteration
  integer :: i,j,k,Ni
  integer,dimension(2) :: ix,iy
  logical :: reff_zero
  real :: minv,maxv
  real :: maxp,minp
  integer,dimension(:),save,  allocatable :: & ! Dimensions nPoints
                  seed    !  It is recommended that the seed is set to a different value for each model
                          !  gridbox it is called on, as it is possible that the choice of the same 
                          !  seed value every time may introduce some statistical bias in the results, 
                          !  particularly for low values of NCOL.
!$OMP THREADPRIVATE(seed)
  real,dimension(:),allocatable :: rseed    !  It is recommended that the seed is set to a different value for each model
  ! Types used in one iteration
  type(cosp_gridbox) :: gbx_it
  type(cosp_subgrid) :: sgx_it
  type(cosp_vgrid)   :: vgrid_it
  type(cosp_sgradar) :: sgradar_it
  type(cosp_sglidar) :: sglidar_it
  type(cosp_isccp)   :: isccp_it
  type(cosp_misr)    :: misr_it
  type(cosp_radarstats) :: stradar_it
  type(cosp_lidarstats) :: stlidar_it
  
  logical,save :: first_cosp=.true.
!$OMP THREADPRIVATE(first_cosp)
  
  !++++++++++ Dimensions ++++++++++++
  npoints  = gbx%Npoints
  nlevels  = gbx%Nlevels
  nhydro   = gbx%Nhydro

!++++++++++ Apply sanity checks to inputs ++++++++++
!  call cosp_check_input('longitude',gbx%longitude,min_val=0.0,max_val=360.0)
  call cosp_check_input('longitude',gbx%longitude,min_val=-180.0,max_val=180.0)
  call cosp_check_input('latitude',gbx%latitude,min_val=-90.0,max_val=90.0)
  call cosp_check_input('dlev',gbx%dlev,min_val=0.0)
  call cosp_check_input('p',gbx%p,min_val=0.0)
  call cosp_check_input('ph',gbx%ph,min_val=0.0)
  call cosp_check_input('T',gbx%T,min_val=0.0)
  call cosp_check_input('q',gbx%q,min_val=0.0)
  call cosp_check_input('sh',gbx%sh,min_val=0.0)
  call cosp_check_input('dtau_s',gbx%dtau_s,min_val=0.0)
  call cosp_check_input('dtau_c',gbx%dtau_c,min_val=0.0)
  call cosp_check_input('dem_s',gbx%dem_s,min_val=0.0,max_val=1.0)
  call cosp_check_input('dem_c',gbx%dem_c,min_val=0.0,max_val=1.0)
  ! Point information (Npoints)
  call cosp_check_input('land',gbx%land,min_val=0.0,max_val=1.0)
  call cosp_check_input('psfc',gbx%psfc,min_val=0.0)
  call cosp_check_input('sunlit',gbx%sunlit,min_val=0.0,max_val=1.0)
  call cosp_check_input('skt',gbx%skt,min_val=0.0)
  ! TOTAL and CONV cloud fraction for SCOPS
  call cosp_check_input('tca',gbx%tca,min_val=0.0,max_val=1.0)
  call cosp_check_input('cca',gbx%cca,min_val=0.0,max_val=1.0)
  ! Precipitation fluxes on model levels
  call cosp_check_input('rain_ls',gbx%rain_ls,min_val=0.0)
  call cosp_check_input('rain_cv',gbx%rain_cv,min_val=0.0)
  call cosp_check_input('snow_ls',gbx%snow_ls,min_val=0.0)
  call cosp_check_input('snow_cv',gbx%snow_cv,min_val=0.0)
  call cosp_check_input('grpl_ls',gbx%grpl_ls,min_val=0.0)
  ! Hydrometeors concentration and distribution parameters
  call cosp_check_input('mr_hydro',gbx%mr_hydro,min_val=0.0)
  ! Effective radius [m]. (Npoints,Nlevels,Nhydro)
  call cosp_check_input('Reff',gbx%Reff,min_val=0.0)
  reff_zero=.true.
  if (any(gbx%Reff > 1.e-8)) then
     reff_zero=.false.
      ! reff_zero == .false.
      !     and gbx%use_reff == .true.   Reff use in radar and lidar
      !     and reff_zero    == .false.  Reff use in lidar and set to 0 for radar
  endif
!  if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
!        print *, '---------- COSP ERROR ------------'
!        print *, ''
!        print *, 'use_reff==.true. but Reff is always zero'
!        print *, ''
!        print *, '----------------------------------'
!        stop
!  endif
  if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
        gbx%Reff = default_lidar_reff
        print *, '---------- COSP WARNING ------------'
        print *, ''
        print *, 'Using default Reff in lidar simulations'
        print *, ''
        print *, '----------------------------------'
  endif
  
  ! Aerosols concentration and distribution parameters
  call cosp_check_input('conc_aero',gbx%conc_aero,min_val=0.0)
  ! Checks for CRM mode
  if (ncolumns == 1) then
     if (gbx%use_precipitation_fluxes) then
        print *, '---------- COSP ERROR ------------'
        print *, ''
        print *, 'Use of precipitation fluxes not supported in CRM mode (Ncolumns=1)'
        print *, ''
        print *, '----------------------------------'
        stop
     endif
     if ((maxval(gbx%dtau_c) > 0.0).or.(maxval(gbx%dem_c) > 0.0)) then
        print *, '---------- COSP ERROR ------------'
        print *, ''
        print *, ' dtau_c > 0.0 or dem_c > 0.0. In CRM mode (Ncolumns=1), '
        print *, ' the optical depth (emmisivity) of all clouds must be '
        print *, ' passed through dtau_s (dem_s)'
        print *, ''
        print *, '----------------------------------'
        stop
     endif
  endif

  if (first_cosp) then   
   ! We base the seed in the decimal part of the surface pressure.
     allocate(seed(npoints))

     allocate(rseed(klon_glo))
     CALL gather(gbx%psfc,rseed)
     call bcast(rseed)
!   seed = int(gbx%psfc) ! This is to avoid division by zero when Npoints = 1   
      ! Roj Oct/2008 ... Note: seed value of 0 caused me some problems + I want to 
      ! randomize for each call to COSP even when Npoints ==1
     minp = minval(rseed)
     maxp = maxval(rseed)
   
     if (npoints .gt. 1) THEN
       seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
     else
       seed=int(gbx%psfc-minp)
     endif

     deallocate(rseed)
     first_cosp=.false.
   endif
    
   if (gbx%Npoints_it >= gbx%Npoints) then ! One iteration gbx%Npoints
        call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
   else ! Several iterations to save memory
        niter = gbx%Npoints/gbx%Npoints_it ! Integer division
        if (niter*gbx%Npoints_it < gbx%Npoints) niter = niter + 1
        do i=1,niter
            i_first = (i-1)*gbx%Npoints_it + 1
            i_last  = i_first + gbx%Npoints_it - 1
            i_last  = min(i_last,gbx%Npoints)
            ni = i_last - i_first + 1
            if (i == 1) then
                ! Allocate types for all but last iteration
                call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
                                            gbx%do_ray,gbx%melt_lay,gbx%k2,ni,nlevels,ncolumns,n_hydro,gbx%Nprmts_max_hydro, &
                                            gbx%Naero,gbx%Nprmts_max_aero,ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
                                            gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
                                            gbx%use_precipitation_fluxes,gbx%use_reff, &
                                            gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
                                            gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
                                            gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
                                            gbx_it)
                call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
                call construct_cosp_subgrid(ni, ncolumns, nlevels, sgx_it)
                call construct_cosp_sgradar(cfg,ni,ncolumns,nlevels,n_hydro,sgradar_it)
                call construct_cosp_sglidar(cfg,ni,ncolumns,nlevels,n_hydro,parasol_nrefl,sglidar_it)
                call construct_cosp_isccp(cfg,ni,ncolumns,nlevels,isccp_it)
                call construct_cosp_misr(cfg,ni,misr_it)
                call construct_cosp_radarstats(cfg,ni,ncolumns,vgrid%Nlvgrid,n_hydro,stradar_it)
                call construct_cosp_lidarstats(cfg,ni,ncolumns,vgrid%Nlvgrid,n_hydro,parasol_nrefl,stlidar_it)
            elseif (i == niter) then ! last iteration
                call free_cosp_gridbox(gbx_it,.true.)
                call free_cosp_subgrid(sgx_it)
                call free_cosp_vgrid(vgrid_it)
                call free_cosp_sgradar(sgradar_it)
                call free_cosp_sglidar(sglidar_it)
                call free_cosp_isccp(isccp_it)
                call free_cosp_misr(misr_it)
                call free_cosp_radarstats(stradar_it)
                call free_cosp_lidarstats(stlidar_it)
                ! Allocate types for iterations
                call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
                                            gbx%do_ray,gbx%melt_lay,gbx%k2,ni,nlevels,ncolumns,n_hydro,gbx%Nprmts_max_hydro, &
                                            gbx%Naero,gbx%Nprmts_max_aero,ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
                                            gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
                                            gbx%use_precipitation_fluxes,gbx%use_reff, &
                                            gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
                                            gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
                                            gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
                                            gbx_it)
                ! --- Copy arrays without Npoints as dimension ---
                gbx_it%dist_prmts_hydro = gbx%dist_prmts_hydro
                gbx_it%dist_type_aero   = gbx_it%dist_type_aero
                call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
                call construct_cosp_subgrid(ni, ncolumns, nlevels, sgx_it)
                call construct_cosp_sgradar(cfg,ni,ncolumns,nlevels,n_hydro,sgradar_it)
                call construct_cosp_sglidar(cfg,ni,ncolumns,nlevels,n_hydro,parasol_nrefl,sglidar_it)
                call construct_cosp_isccp(cfg,ni,ncolumns,nlevels,isccp_it)
                call construct_cosp_misr(cfg,ni,misr_it)
                call construct_cosp_radarstats(cfg,ni,ncolumns,vgrid%Nlvgrid,n_hydro,stradar_it)
                call construct_cosp_lidarstats(cfg,ni,ncolumns,vgrid%Nlvgrid,n_hydro,parasol_nrefl,stlidar_it)
            endif
            ! --- Copy sections of arrays with Npoints as dimension ---
            ix=(/i_first,i_last/)
            iy=(/1,ni/)
            call cosp_gridbox_cpsection(ix,iy,gbx,gbx_it)
              ! These serve as initialisation of *_it types
            call cosp_subgrid_cpsection(ix,iy,sgx,sgx_it)
            if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar,sgradar_it)
            if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar,sglidar_it)
            if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp,isccp_it)
            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr,misr_it)
            if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar,stradar_it)
            if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar,stlidar_it)
            print *,'---------ix: ',ix
            call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
                           sglidar_it,isccp_it,misr_it,stradar_it,stlidar_it)
            
            ! --- Copy results to output structures ---
!             call cosp_gridbox_cphp(gbx_it,gbx)
            ix=(/1,ni/)
            iy=(/i_first,i_last/)
            call cosp_subgrid_cpsection(ix,iy,sgx_it,sgx)
            if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar_it,sgradar)
            if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar_it,sglidar)
            if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp_it,isccp)
            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr_it,misr)
            if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar_it,stradar)
            if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar_it,stlidar)
        enddo
        ! Deallocate types
        call free_cosp_gridbox(gbx_it,.true.)
        call free_cosp_subgrid(sgx_it)
        call free_cosp_vgrid(vgrid_it)
        call free_cosp_sgradar(sgradar_it)
        call free_cosp_sglidar(sglidar_it)
        call free_cosp_isccp(isccp_it)
        call free_cosp_misr(misr_it)
        call free_cosp_radarstats(stradar_it)
        call free_cosp_lidarstats(stlidar_it)
   endif

    
END SUBROUTINE cosp

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!--------------------- SUBROUTINE COSP_ITER ----------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)

  ! Arguments
  integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
  integer,dimension(:),intent(in) :: seed
  type(cosp_config),intent(in) :: cfg   ! Configuration options
  type(cosp_vgrid),intent(in) :: vgrid   ! Information on vertical grid of stats
  type(cosp_gridbox),intent(inout) :: gbx
  type(cosp_subgrid),intent(inout) :: sgx   ! Subgrid info
  type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
  type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
  type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
  type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
  type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
  type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator

  ! Local variables 
  integer :: Npoints   ! Number of gridpoints
  integer :: Ncolumns  ! Number of subcolumns
  integer :: Nlevels   ! Number of levels
  integer :: Nhydro    ! Number of hydrometeors
  integer :: Niter     ! Number of calls to cosp_simulator
  integer :: i,j,k
  integer :: I_HYDRO
  real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
  integer,parameter :: scops_debug=0    !  set to non-zero value to print out inputs for debugging in SCOPS
  
  real,dimension(:, :),allocatable :: cca_scops,ls_p_rate,cv_p_rate, &
                     tca_scops ! Cloud cover in each model level (HORIZONTAL gridbox fraction) of total cloud.
                               ! Levels are from TOA to SURFACE. (nPoints, nLev)
  real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
                                                                     ! Levels are from SURFACE to TOA
  real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric dens
  type(cosp_sghydro) :: sghydro   ! Subgrid info for hydrometeors en each iteration

  
  !++++++++++ Dimensions ++++++++++++
  npoints  = gbx%Npoints
  ncolumns = gbx%Ncolumns
  nlevels  = gbx%Nlevels
  nhydro   = gbx%Nhydro
    
   
  !++++++++++ Climate/NWP mode ++++++++++  
  if (ncolumns > 1) then
        !++++++++++ Subgrid sampling ++++++++++
        ! Allocate arrays before calling SCOPS
        allocate(frac_ls(npoints,nlevels),frac_cv(npoints,nlevels),prec_ls(npoints,nlevels),prec_cv(npoints,nlevels))
        allocate(tca_scops(npoints,nlevels),cca_scops(npoints,nlevels), &
                ls_p_rate(npoints,nlevels),cv_p_rate(npoints,nlevels))
        ! Initialize to zero
        frac_ls=0.0
        prec_ls=0.0
        frac_cv=0.0
        prec_cv=0.0
        ! Cloud fractions for SCOPS from TOA to SFC
        tca_scops = gbx%tca(:,nlevels:1:-1)
        cca_scops = gbx%cca(:,nlevels:1:-1)
        
        ! Call to SCOPS
        ! strat and conv arrays are passed with levels from TOA to SURFACE.
        call scops(npoints,nlevels,ncolumns,seed,tca_scops,cca_scops,overlap,sgx%frac_out,scops_debug)
        
        ! temporarily use prec_ls/cv to transfer information about precipitation flux into prec_scops
        if(gbx%use_precipitation_fluxes) then
            ls_p_rate(:,nlevels:1:-1)=gbx%rain_ls(:,1:nlevels)+gbx%snow_ls(:,1:nlevels)+gbx%grpl_ls(:,1:nlevels)
            cv_p_rate(:,nlevels:1:-1)=gbx%rain_cv(:,1:nlevels)+gbx%snow_cv(:,1:nlevels)
        else
            ls_p_rate(:,nlevels:1:-1)=gbx%mr_hydro(:,1:nlevels,i_lsrain)+ &
                                      gbx%mr_hydro(:,1:nlevels,i_lssnow)+ &
                                      gbx%mr_hydro(:,1:nlevels,i_lsgrpl)
            cv_p_rate(:,nlevels:1:-1)=gbx%mr_hydro(:,1:nlevels,i_cvrain)+ &
                                      gbx%mr_hydro(:,1:nlevels,i_cvsnow)
        endif
        
        call prec_scops(npoints,nlevels,ncolumns,ls_p_rate,cv_p_rate,sgx%frac_out,sgx%prec_frac)
        
        ! Precipitation fraction
        do j=1,npoints,1
        do k=1,nlevels,1
            do i=1,ncolumns,1
                if (sgx%frac_out (j,i,nlevels+1-k) == i_lsc) frac_ls(j,k)=frac_ls(j,k)+1.
                if (sgx%frac_out (j,i,nlevels+1-k) == i_cvc) frac_cv(j,k)=frac_cv(j,k)+1.
                if (sgx%prec_frac(j,i,nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
                if (sgx%prec_frac(j,i,nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
                if (sgx%prec_frac(j,i,nlevels+1-k) .eq. 3) then
                    prec_cv(j,k)=prec_cv(j,k)+1.
                    prec_ls(j,k)=prec_ls(j,k)+1.
                endif
            enddo  !i
            frac_ls(j,k)=frac_ls(j,k)/ncolumns
            frac_cv(j,k)=frac_cv(j,k)/ncolumns
            prec_ls(j,k)=prec_ls(j,k)/ncolumns
            prec_cv(j,k)=prec_cv(j,k)/ncolumns
        enddo  !k
        enddo  !j
        
         ! Levels from SURFACE to TOA.
        if (npoints*ncolumns*nlevels < 10000) then
            sgx%frac_out(1:npoints,:,1:nlevels)  = sgx%frac_out(1:npoints,:,nlevels:1:-1)
            sgx%prec_frac(1:npoints,:,1:nlevels) = sgx%prec_frac(1:npoints,:,nlevels:1:-1)
        else
            ! This is done within a loop (unvectorized) over nPoints to save memory
            do j=1,npoints
                sgx%frac_out(j,:,1:nlevels)  = sgx%frac_out(j,:,nlevels:1:-1)
                sgx%prec_frac(j,:,1:nlevels) = sgx%prec_frac(j,:,nlevels:1:-1)
            enddo
        endif
       
       ! Deallocate arrays that will no longer be used
        deallocate(tca_scops,cca_scops,ls_p_rate,cv_p_rate)
         
        ! Populate the subgrid arrays
        call construct_cosp_sghydro(npoints,ncolumns,nlevels,nhydro,sghydro)
        do k=1,ncolumns
            !--------- Mixing ratios for clouds and Reff for Clouds and precip -------
            column_frac_out => sgx%frac_out(:,k,:)
            where (column_frac_out == i_lsc)     !+++++++++++ LS clouds ++++++++
                sghydro%mr_hydro(:,k,:,i_lscliq) = gbx%mr_hydro(:,:,i_lscliq)
                sghydro%mr_hydro(:,k,:,i_lscice) = gbx%mr_hydro(:,:,i_lscice)
                
                sghydro%Reff(:,k,:,i_lscliq)     = gbx%Reff(:,:,i_lscliq)
                sghydro%Reff(:,k,:,i_lscice)     = gbx%Reff(:,:,i_lscice)
                sghydro%Reff(:,k,:,i_lsrain)     = gbx%Reff(:,:,i_lsrain)
                sghydro%Reff(:,k,:,i_lssnow)     = gbx%Reff(:,:,i_lssnow)
                sghydro%Reff(:,k,:,i_lsgrpl)     = gbx%Reff(:,:,i_lsgrpl)
            elsewhere (column_frac_out == i_cvc) !+++++++++++ CONV clouds ++++++++
                sghydro%mr_hydro(:,k,:,i_cvcliq) = gbx%mr_hydro(:,:,i_cvcliq) 
                sghydro%mr_hydro(:,k,:,i_cvcice) = gbx%mr_hydro(:,:,i_cvcice) 
                
                sghydro%Reff(:,k,:,i_cvcliq)     = gbx%Reff(:,:,i_cvcliq) 
                sghydro%Reff(:,k,:,i_cvcice)     = gbx%Reff(:,:,i_cvcice) 
                sghydro%Reff(:,k,:,i_cvrain)     = gbx%Reff(:,:,i_cvrain) 
                sghydro%Reff(:,k,:,i_cvsnow)     = gbx%Reff(:,:,i_cvsnow) 
            end where 
            !--------- Precip -------
            if (.not. gbx%use_precipitation_fluxes) then
                where (column_frac_out == i_lsc)  !+++++++++++ LS Precipitation ++++++++
                    sghydro%mr_hydro(:,k,:,i_lsrain) = gbx%mr_hydro(:,:,i_lsrain)
                    sghydro%mr_hydro(:,k,:,i_lssnow) = gbx%mr_hydro(:,:,i_lssnow)
                    sghydro%mr_hydro(:,k,:,i_lsgrpl) = gbx%mr_hydro(:,:,i_lsgrpl)
                elsewhere (column_frac_out == i_cvc) !+++++++++++ CONV Precipitation ++++++++
                    sghydro%mr_hydro(:,k,:,i_cvrain) = gbx%mr_hydro(:,:,i_cvrain) 
                    sghydro%mr_hydro(:,k,:,i_cvsnow) = gbx%mr_hydro(:,:,i_cvsnow) 
                end where 
            endif
        enddo
        ! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values
        do k=1,nlevels
            do j=1,npoints
                !--------- Clouds -------
                if (frac_ls(j,k) .ne. 0.) then
                    sghydro%mr_hydro(j,:,k,i_lscliq) = sghydro%mr_hydro(j,:,k,i_lscliq)/frac_ls(j,k)
                    sghydro%mr_hydro(j,:,k,i_lscice) = sghydro%mr_hydro(j,:,k,i_lscice)/frac_ls(j,k)
                endif
                if (frac_cv(j,k) .ne. 0.) then
                    sghydro%mr_hydro(j,:,k,i_cvcliq) = sghydro%mr_hydro(j,:,k,i_cvcliq)/frac_cv(j,k)
                    sghydro%mr_hydro(j,:,k,i_cvcice) = sghydro%mr_hydro(j,:,k,i_cvcice)/frac_cv(j,k)
                endif
                !--------- Precip -------
                if (gbx%use_precipitation_fluxes) then
                    if (prec_ls(j,k) .ne. 0.) then
                        gbx%rain_ls(j,k) = gbx%rain_ls(j,k)/prec_ls(j,k)
                        gbx%snow_ls(j,k) = gbx%snow_ls(j,k)/prec_ls(j,k)
                        gbx%grpl_ls(j,k) = gbx%grpl_ls(j,k)/prec_ls(j,k)
                    endif
                    if (prec_cv(j,k) .ne. 0.) then
                        gbx%rain_cv(j,k) = gbx%rain_cv(j,k)/prec_cv(j,k)
                        gbx%snow_cv(j,k) = gbx%snow_cv(j,k)/prec_cv(j,k)
                    endif
                else
                    if (prec_ls(j,k) .ne. 0.) then
                        sghydro%mr_hydro(j,:,k,i_lsrain) = sghydro%mr_hydro(j,:,k,i_lsrain)/prec_ls(j,k)
                        sghydro%mr_hydro(j,:,k,i_lssnow) = sghydro%mr_hydro(j,:,k,i_lssnow)/prec_ls(j,k)
                        sghydro%mr_hydro(j,:,k,i_lsgrpl) = sghydro%mr_hydro(j,:,k,i_lsgrpl)/prec_ls(j,k)
                    endif
                    if (prec_cv(j,k) .ne. 0.) then
                        sghydro%mr_hydro(j,:,k,i_cvrain) = sghydro%mr_hydro(j,:,k,i_cvrain)/prec_cv(j,k)
                        sghydro%mr_hydro(j,:,k,i_cvsnow) = sghydro%mr_hydro(j,:,k,i_cvsnow)/prec_cv(j,k)
                    endif
                endif  
            enddo !k
        enddo !j
        deallocate(frac_ls,prec_ls,frac_cv,prec_cv)
        
        if (gbx%use_precipitation_fluxes) then
            ! convert precipitation flux into mixing ratio
            call pf_to_mr(npoints,nlevels,ncolumns,gbx%rain_ls,gbx%snow_ls,gbx%grpl_ls, &
                        gbx%rain_cv,gbx%snow_cv,sgx%prec_frac,gbx%p,gbx%T, &
                        sghydro%mr_hydro(:,:,:,i_lsrain),sghydro%mr_hydro(:,:,:,i_lssnow),sghydro%mr_hydro(:,:,:,i_lsgrpl), &
                        sghydro%mr_hydro(:,:,:,i_cvrain),sghydro%mr_hydro(:,:,:,i_cvsnow))
       endif
   !++++++++++ CRM mode ++++++++++
   else
      sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
      sghydro%Reff(:,1,:,:) = gbx%Reff
      !--------- Clouds -------
      where ((gbx%dtau_s > 0.0))
             sgx%frac_out(:,1,:) = 1  ! Subgrid cloud array. Dimensions (Npoints,Ncolumns,Nlevels)
      endwhere
   endif ! Ncolumns > 1
  
   
   !++++++++++ Simulator ++++++++++
    call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,stradar,stlidar)

    ! Deallocate subgrid arrays
    call free_cosp_sghydro(sghydro)
END SUBROUTINE cosp_iter

END MODULE mod_cosp