Changeset 2435 for LMDZ5/branches/testing/libf/phylmd/cosp/cosp.F90

Timestamp:

Jan 28, 2016, 5:02:13 PM (9 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2396:2434 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/cosp/cosp.F90 (modified) (24 diffs)

LMDZ5/branches/testing/libf/phylmd/cosp/cosp.F90

@@ -23 +23 @@
 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-!!#include "cosp_defs.h"
+#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
+  USE MOD_COSP_MODIS_SIMULATOR
   IMPLICIT NONE
 
@@ -37 +36 @@
 !--------------------- SUBROUTINE COSP ---------------------------
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
-
+!#ifdef RTTOV
+!SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
+!#else
+SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
+!#endif
   ! Arguments
   integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
@@ -50 +52 @@
   type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
   type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
+  type(cosp_modis),intent(inout)   :: modis   ! Output from MODIS simulator
+!#ifdef RTTOV
+!  type(cosp_rttov),intent(inout)   :: rttov   ! Output from RTTOV
+!#endif
   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
@@ -59 +65 @@
   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 :: i,Ni
   integer,dimension(2) :: ix,iy
   logical :: reff_zero
-  real :: minv,maxv
   real :: maxp,minp
-  integer,dimension(:),save,  allocatable :: & ! Dimensions nPoints
+  integer,dimension(:),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
@@ -78 +81 @@
   type(cosp_sglidar) :: sglidar_it
   type(cosp_isccp)   :: isccp_it
+  type(cosp_modis)   :: modis_it
   type(cosp_misr)    :: misr_it
+!#ifdef RTTOV
+!  type(cosp_rttov)   :: rttov_it
+!#endif
   type(cosp_radarstats) :: stradar_it
   type(cosp_lidarstats) :: stlidar_it
-  
-  logical,save :: first_cosp=.TRUE.
-!$OMP THREADPRIVATE(first_cosp)
-  
-  !++++++++++ Dimensions ++++++++++++
+
+!++++++++++ Dimensions ++++++++++++
   Npoints  = gbx%Npoints
   Nlevels  = gbx%Nlevels
   Nhydro   = gbx%Nhydro
+
+!++++++++++ Depth of model layers ++++++++++++
+  do i=1,Nlevels-1
+    gbx%dlev(:,i) = gbx%zlev_half(:,i+1) - gbx%zlev_half(:,i)
+  enddo
+  gbx%dlev(:,Nlevels) = 2.0*(gbx%zlev(:,Nlevels) - gbx%zlev_half(:,Nlevels))
 
 !++++++++++ Apply sanity checks to inputs ++++++++++
@@ -129 +139 @@
       !     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
@@ -170 +172 @@
   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   
+   allocate(seed(Npoints))
+   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
-    
+   minp = minval(gbx%psfc)
+   maxp = maxval(gbx%psfc)
+   if (Npoints .gt. 1) seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
+   ! Below it's how it was done in the original implementation of the ISCCP simulator. 
+   ! The one above is better for offline data, when you may have packed data 
+   ! that subsamples the decimal fraction of the surface pressure. 
+!    if (Npoints .gt. 1) seed=(gbx%psfc-int(gbx%psfc))*1000000 
+
+  
    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)
+!#ifdef RTTOV
+!        call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
+!#else
+        call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
+!#endif
    else ! Several iterations to save memory
         Niter = gbx%Npoints/gbx%Npoints_it ! Integer division
@@ -205 +202 @@
             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, &
+                call construct_cosp_gridbox(gbx%time,gbx%time_bnds,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, &
@@ -219 +217 @@
                 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_modis(cfg, Ni, modis_it)
                 call construct_cosp_misr(cfg,Ni,misr_it)
+!#ifdef RTTOV
+!                call construct_cosp_rttov(Ni,gbx%nchan,rttov_it)
+!#endif
                 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)
@@ -229 +231 @@
                 call free_cosp_sglidar(sglidar_it)
                 call free_cosp_isccp(isccp_it)
+                call free_cosp_modis(modis_it)
                 call free_cosp_misr(misr_it)
+!#ifdef RTTOV
+!                call free_cosp_rttov(rttov_it)
+!#endif
                 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, &
+                call construct_cosp_gridbox(gbx%time,gbx%time_bnds,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, &
@@ -250 +257 @@
                 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_modis(cfg,Ni, modis_it)
                 call construct_cosp_misr(cfg,Ni,misr_it)
+!#ifdef RTTOV 
+!                call construct_cosp_rttov(Ni,gbx%nchan,rttov_it) 
+!#endif 
                 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)
@@ -263 +274 @@
             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%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis,modis_it)
             if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr,misr_it)
+!#ifdef RTTOV 
+!            if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov,rttov_it) 
+!#endif
             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
+!#ifdef RTTOV
+!            call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
+!                           sglidar_it,isccp_it,misr_it,modis_it,rttov_it,stradar_it,stlidar_it)
+!#else
             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)
-            
+                           sglidar_it,isccp_it,misr_it,modis_it,stradar_it,stlidar_it)
+!#endif
             ! --- Copy results to output structures ---
-!             call cosp_gridbox_cphp(gbx_it,gbx)
             ix=(/1,Ni/)
             iy=(/i_first,i_last/)
@@ -278 +295 @@
             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%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis_it,modis)
             if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr_it,misr)
+!#ifdef RTTOV 
+!            if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov_it,rttov) 
+!#endif 
             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)
@@ -289 +310 @@
         call free_cosp_sglidar(sglidar_it)
         call free_cosp_isccp(isccp_it)
+        call free_cosp_modis(modis_it)
         call free_cosp_misr(misr_it)
+!#ifdef RTTOV 
+!        call free_cosp_rttov(rttov_it) 
+!#endif
         call free_cosp_radarstats(stradar_it)
         call free_cosp_lidarstats(stlidar_it)
    endif
+   deallocate(seed)
 
     
@@ -300 +326 @@
 !--------------------- SUBROUTINE COSP_ITER ----------------------
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
-
+!#ifdef RTTOV
+!SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
+!#else
+SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
+!#endif
   ! Arguments
   integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
@@ -313 +342 @@
   type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
   type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
+  type(cosp_modis),intent(inout)   :: modis   ! Output from MODIS simulator
+!#ifdef RTTOV
+!  type(cosp_rttov),intent(inout)   :: rttov   ! Output from RTTOV
+!#endif
   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
@@ -321 +354 @@
   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
+  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(:,:),pointer :: column_prec_out ! Array with one column of prec_frac
+  integer :: 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.
@@ -332 +364 @@
   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
+  real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric density
   type(cosp_sghydro) :: sghydro   ! Subgrid info for hydrometeors en each iteration
 
@@ -342 +374 @@
   Nhydro   = gbx%Nhydro
     
-   
   !++++++++++ Climate/NWP mode ++++++++++  
   if (Ncolumns > 1) then
@@ -411 +442 @@
        ! 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)
@@ -420 +451 @@
                 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)
+
+                sghydro%Np(:,k,:,I_LSCLIQ)     = gbx%Np(:,:,I_LSCLIQ)
+                sghydro%Np(:,k,:,I_LSCICE)     = gbx%Np(:,:,I_LSCICE)
+
             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 
+                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%Np(:,k,:,I_CVCLIQ)     = gbx%Np(:,:,I_CVCLIQ)
+                sghydro%Np(:,k,:,I_CVCICE)     = gbx%Np(:,:,I_CVCICE)
+
+            end where
+            column_prec_out => sgx%prec_frac(:,k,:)
+            where ((column_prec_out == 1) .or. (column_prec_out == 3) )  !++++ LS precip ++++
+                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)
+
+                sghydro%Np(:,k,:,I_LSRAIN)     = gbx%Np(:,:,I_LSRAIN)
+                sghydro%Np(:,k,:,I_LSSNOW)     = gbx%Np(:,:,I_LSSNOW)
+                sghydro%Np(:,k,:,I_LSGRPL)     = gbx%Np(:,:,I_LSGRPL)
+            elsewhere ((column_prec_out == 2) .or. (column_prec_out == 3)) !++++ CONV precip ++++
+                sghydro%Reff(:,k,:,I_CVRAIN) = gbx%Reff(:,:,I_CVRAIN)
+                sghydro%Reff(:,k,:,I_CVSNOW) = gbx%Reff(:,:,I_CVSNOW)
+
+                sghydro%Np(:,k,:,I_CVRAIN)     = gbx%Np(:,:,I_CVRAIN)
+                sghydro%Np(:,k,:,I_CVSNOW)     = gbx%Np(:,:,I_CVSNOW)
+            end where
             !--------- Precip -------
             if (.not. gbx%use_precipitation_fluxes) then
@@ -442 +492 @@
                     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) 
+                    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
@@ -486 +536 @@
         
         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
+        
+#ifdef MMF_V3p5_TWO_MOMENT
+
+        write(*,*) 'Precipitation Flux to Mixing Ratio conversion not (yet?) supported ', &
+               'for MMF3.5 Two Moment Microphysics'
+        stop
+#else
+            ! Density
+            allocate(rho(Npoints,Nlevels))
+            I_HYDRO = I_LSRAIN
+            call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
+                    n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
+                    g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
+                    gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
+                    gbx%rain_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
+            I_HYDRO = I_LSSNOW
+            call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
+                    n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
+                    g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
+                    gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
+                    gbx%snow_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
+            I_HYDRO = I_CVRAIN
+            call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,2., &
+                    n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
+                    g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
+                    gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
+                    gbx%rain_cv,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
+            I_HYDRO = I_CVSNOW
+            call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,2., &
+                    n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
+                    g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
+                    gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
+                    gbx%snow_cv,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
+            I_HYDRO = I_LSGRPL
+            call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
+                    n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
+                    g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
+                    gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
+                    gbx%grpl_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
+            if(allocated(rho)) deallocate(rho)
+#endif
+
+        endif
    !++++++++++ CRM mode ++++++++++
    else
+      call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
       sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
       sghydro%Reff(:,1,:,:) = gbx%Reff
+      sghydro%Np(:,1,:,:) = gbx%Np      ! added by Roj with Quickbeam V3.0
+      
       !--------- Clouds -------
       where ((gbx%dtau_s > 0.0))
@@ -502 +592 @@
    endif ! Ncolumns > 1
   
-   
    !++++++++++ Simulator ++++++++++
-    call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,stradar,stlidar)
+!#ifdef RTTOV
+!    call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
+!#else
+    call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
+!#endif
 
     ! Deallocate subgrid arrays