source: LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/cosp.F90 @ 5404

Last change on this file since 5404 was 1318, checked in by yann meurdesoif, 15 years ago

YM : Parallelisation COSP MPI+OpenMP

File size: 26.2 KB
Line 
1! (c) British Crown Copyright 2008, the Met Office.
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without modification, are permitted
5! provided that the following conditions are met:
6!
7!     * Redistributions of source code must retain the above copyright notice, this list
8!       of conditions and the following disclaimer.
9!     * Redistributions in binary form must reproduce the above copyright notice, this list
10!       of conditions and the following disclaimer in the documentation and/or other materials
11!       provided with the distribution.
12!     * Neither the name of the Met Office nor the names of its contributors may be used
13!       to endorse or promote products derived from this software without specific prior written
14!       permission.
15!
16! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
17! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
18! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
19! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
22! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
23! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25MODULE MOD_COSP
26  USE MOD_COSP_TYPES
27  USE MOD_COSP_SIMULATOR
28  USE mod_phys_lmdz_para
29  USE mod_grid_phy_lmdz
30  IMPLICIT NONE
31
32CONTAINS
33
34
35!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36!--------------------- SUBROUTINE COSP ---------------------------
37!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
39
40  ! Arguments
41  integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
42  integer,intent(in) :: Ncolumns
43  type(cosp_config),intent(in) :: cfg   ! Configuration options
44  type(cosp_vgrid),intent(in) :: vgrid   ! Information on vertical grid of stats
45  type(cosp_gridbox),intent(inout) :: gbx
46  type(cosp_subgrid),intent(inout) :: sgx   ! Subgrid info
47  type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
48  type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
49  type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
50  type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
51  type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
52  type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
53
54  ! Local variables
55  integer :: Npoints   ! Number of gridpoints
56  integer :: Nlevels   ! Number of levels
57  integer :: Nhydro    ! Number of hydrometeors
58  integer :: Niter     ! Number of calls to cosp_simulator
59  integer :: i_first,i_last ! First and last gridbox to be processed in each iteration
60  integer :: i,j,k,Ni
61  integer,dimension(2) :: ix,iy
62  logical :: reff_zero
63  real :: minv,maxv
64  real :: maxp,minp
65  integer,dimension(:),save,  allocatable :: & ! Dimensions nPoints
66                  seed    !  It is recommended that the seed is set to a different value for each model
67                          !  gridbox it is called on, as it is possible that the choice of the same
68                          !  seed value every time may introduce some statistical bias in the results,
69                          !  particularly for low values of NCOL.
70!$OMP THREADPRIVATE(seed)
71  real,dimension(:),allocatable :: rseed    !  It is recommended that the seed is set to a different value for each model
72  ! Types used in one iteration
73  type(cosp_gridbox) :: gbx_it
74  type(cosp_subgrid) :: sgx_it
75  type(cosp_vgrid)   :: vgrid_it
76  type(cosp_sgradar) :: sgradar_it
77  type(cosp_sglidar) :: sglidar_it
78  type(cosp_isccp)   :: isccp_it
79  type(cosp_misr)    :: misr_it
80  type(cosp_radarstats) :: stradar_it
81  type(cosp_lidarstats) :: stlidar_it
82 
83  logical,save :: first_cosp=.TRUE.
84!$OMP THREADPRIVATE(first_cosp)
85 
86  !++++++++++ Dimensions ++++++++++++
87  Npoints  = gbx%Npoints
88  Nlevels  = gbx%Nlevels
89  Nhydro   = gbx%Nhydro
90
91!++++++++++ Apply sanity checks to inputs ++++++++++
92!  call cosp_check_input('longitude',gbx%longitude,min_val=0.0,max_val=360.0)
93  call cosp_check_input('longitude',gbx%longitude,min_val=-180.0,max_val=180.0)
94  call cosp_check_input('latitude',gbx%latitude,min_val=-90.0,max_val=90.0)
95  call cosp_check_input('dlev',gbx%dlev,min_val=0.0)
96  call cosp_check_input('p',gbx%p,min_val=0.0)
97  call cosp_check_input('ph',gbx%ph,min_val=0.0)
98  call cosp_check_input('T',gbx%T,min_val=0.0)
99  call cosp_check_input('q',gbx%q,min_val=0.0)
100  call cosp_check_input('sh',gbx%sh,min_val=0.0)
101  call cosp_check_input('dtau_s',gbx%dtau_s,min_val=0.0)
102  call cosp_check_input('dtau_c',gbx%dtau_c,min_val=0.0)
103  call cosp_check_input('dem_s',gbx%dem_s,min_val=0.0,max_val=1.0)
104  call cosp_check_input('dem_c',gbx%dem_c,min_val=0.0,max_val=1.0)
105  ! Point information (Npoints)
106  call cosp_check_input('land',gbx%land,min_val=0.0,max_val=1.0)
107  call cosp_check_input('psfc',gbx%psfc,min_val=0.0)
108  call cosp_check_input('sunlit',gbx%sunlit,min_val=0.0,max_val=1.0)
109  call cosp_check_input('skt',gbx%skt,min_val=0.0)
110  ! TOTAL and CONV cloud fraction for SCOPS
111  call cosp_check_input('tca',gbx%tca,min_val=0.0,max_val=1.0)
112  call cosp_check_input('cca',gbx%cca,min_val=0.0,max_val=1.0)
113  ! Precipitation fluxes on model levels
114  call cosp_check_input('rain_ls',gbx%rain_ls,min_val=0.0)
115  call cosp_check_input('rain_cv',gbx%rain_cv,min_val=0.0)
116  call cosp_check_input('snow_ls',gbx%snow_ls,min_val=0.0)
117  call cosp_check_input('snow_cv',gbx%snow_cv,min_val=0.0)
118  call cosp_check_input('grpl_ls',gbx%grpl_ls,min_val=0.0)
119  ! Hydrometeors concentration and distribution parameters
120  call cosp_check_input('mr_hydro',gbx%mr_hydro,min_val=0.0)
121  ! Effective radius [m]. (Npoints,Nlevels,Nhydro)
122  call cosp_check_input('Reff',gbx%Reff,min_val=0.0)
123  reff_zero=.true.
124  if (any(gbx%Reff > 1.e-8)) then
125     reff_zero=.false.
126      ! reff_zero == .false.
127      !     and gbx%use_reff == .true.   Reff use in radar and lidar
128      !     and reff_zero    == .false.  Reff use in lidar and set to 0 for radar
129  endif
130  if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
131        print *, '---------- COSP ERROR ------------'
132        print *, ''
133        print *, 'use_reff==.true. but Reff is always zero'
134        print *, ''
135        print *, '----------------------------------'
136        stop
137  endif
138  if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
139        gbx%Reff = DEFAULT_LIDAR_REFF
140        print *, '---------- COSP WARNING ------------'
141        print *, ''
142        print *, 'Using default Reff in lidar simulations'
143        print *, ''
144        print *, '----------------------------------'
145  endif
146 
147  ! Aerosols concentration and distribution parameters
148  call cosp_check_input('conc_aero',gbx%conc_aero,min_val=0.0)
149  ! Checks for CRM mode
150  if (Ncolumns == 1) then
151     if (gbx%use_precipitation_fluxes) then
152        print *, '---------- COSP ERROR ------------'
153        print *, ''
154        print *, 'Use of precipitation fluxes not supported in CRM mode (Ncolumns=1)'
155        print *, ''
156        print *, '----------------------------------'
157        stop
158     endif
159     if ((maxval(gbx%dtau_c) > 0.0).or.(maxval(gbx%dem_c) > 0.0)) then
160        print *, '---------- COSP ERROR ------------'
161        print *, ''
162        print *, ' dtau_c > 0.0 or dem_c > 0.0. In CRM mode (Ncolumns=1), '
163        print *, ' the optical depth (emmisivity) of all clouds must be '
164        print *, ' passed through dtau_s (dem_s)'
165        print *, ''
166        print *, '----------------------------------'
167        stop
168     endif
169  endif
170
171  if (first_cosp) then   
172   ! We base the seed in the decimal part of the surface pressure.
173     allocate(seed(Npoints))
174
175     allocate(rseed(klon_glo))
176     CALL gather(gbx%psfc,rseed)
177     call bcast(rseed)
178!   seed = int(gbx%psfc) ! This is to avoid division by zero when Npoints = 1   
179      ! Roj Oct/2008 ... Note: seed value of 0 caused me some problems + I want to
180      ! randomize for each call to COSP even when Npoints ==1
181     minp = minval(rseed)
182     maxp = maxval(rseed)
183   
184     if (Npoints .gt. 1) THEN
185       seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
186     else
187       seed=int(gbx%psfc-minp)
188     endif
189
190     deallocate(rseed)
191     first_cosp=.false.
192   endif
193   
194   if (gbx%Npoints_it >= gbx%Npoints) then ! One iteration gbx%Npoints
195        call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
196   else ! Several iterations to save memory
197        Niter = gbx%Npoints/gbx%Npoints_it ! Integer division
198        if (Niter*gbx%Npoints_it < gbx%Npoints) Niter = Niter + 1
199        do i=1,Niter
200            i_first = (i-1)*gbx%Npoints_it + 1
201            i_last  = i_first + gbx%Npoints_it - 1
202            i_last  = min(i_last,gbx%Npoints)
203            Ni = i_last - i_first + 1
204            if (i == 1) then
205                ! Allocate types for all but last iteration
206                call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
207                                            gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels,Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
208                                            gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
209                                            gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
210                                            gbx%use_precipitation_fluxes,gbx%use_reff, &
211                                            gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
212                                            gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
213                                            gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
214                                            gbx_it)
215                call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
216                call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
217                call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
218                call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
219                call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
220                call construct_cosp_misr(cfg,Ni,misr_it)
221                call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
222                call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
223            elseif (i == Niter) then ! last iteration
224                call free_cosp_gridbox(gbx_it,.true.)
225                call free_cosp_subgrid(sgx_it)
226                call free_cosp_vgrid(vgrid_it)
227                call free_cosp_sgradar(sgradar_it)
228                call free_cosp_sglidar(sglidar_it)
229                call free_cosp_isccp(isccp_it)
230                call free_cosp_misr(misr_it)
231                call free_cosp_radarstats(stradar_it)
232                call free_cosp_lidarstats(stlidar_it)
233                ! Allocate types for iterations
234                call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
235                                            gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels,Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
236                                            gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
237                                            gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
238                                            gbx%use_precipitation_fluxes,gbx%use_reff, &
239                                            gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
240                                            gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
241                                            gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
242                                            gbx_it)
243                ! --- Copy arrays without Npoints as dimension ---
244                gbx_it%dist_prmts_hydro = gbx%dist_prmts_hydro
245                gbx_it%dist_type_aero   = gbx_it%dist_type_aero
246                call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
247                call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
248                call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
249                call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
250                call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
251                call construct_cosp_misr(cfg,Ni,misr_it)
252                call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
253                call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
254            endif
255            ! --- Copy sections of arrays with Npoints as dimension ---
256            ix=(/i_first,i_last/)
257            iy=(/1,Ni/)
258            call cosp_gridbox_cpsection(ix,iy,gbx,gbx_it)
259              ! These serve as initialisation of *_it types
260            call cosp_subgrid_cpsection(ix,iy,sgx,sgx_it)
261            if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar,sgradar_it)
262            if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar,sglidar_it)
263            if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp,isccp_it)
264            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr,misr_it)
265            if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar,stradar_it)
266            if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar,stlidar_it)
267            print *,'---------ix: ',ix
268            call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
269                           sglidar_it,isccp_it,misr_it,stradar_it,stlidar_it)
270           
271            ! --- Copy results to output structures ---
272!             call cosp_gridbox_cphp(gbx_it,gbx)
273            ix=(/1,Ni/)
274            iy=(/i_first,i_last/)
275            call cosp_subgrid_cpsection(ix,iy,sgx_it,sgx)
276            if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar_it,sgradar)
277            if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar_it,sglidar)
278            if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp_it,isccp)
279            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr_it,misr)
280            if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar_it,stradar)
281            if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar_it,stlidar)
282        enddo
283        ! Deallocate types
284        call free_cosp_gridbox(gbx_it,.true.)
285        call free_cosp_subgrid(sgx_it)
286        call free_cosp_vgrid(vgrid_it)
287        call free_cosp_sgradar(sgradar_it)
288        call free_cosp_sglidar(sglidar_it)
289        call free_cosp_isccp(isccp_it)
290        call free_cosp_misr(misr_it)
291        call free_cosp_radarstats(stradar_it)
292        call free_cosp_lidarstats(stlidar_it)
293   endif
294
295   
296END SUBROUTINE COSP
297
298!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299!--------------------- SUBROUTINE COSP_ITER ----------------------
300!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
302
303  ! Arguments
304  integer,intent(in) :: overlap !  overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
305  integer,dimension(:),intent(in) :: seed
306  type(cosp_config),intent(in) :: cfg   ! Configuration options
307  type(cosp_vgrid),intent(in) :: vgrid   ! Information on vertical grid of stats
308  type(cosp_gridbox),intent(inout) :: gbx
309  type(cosp_subgrid),intent(inout) :: sgx   ! Subgrid info
310  type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
311  type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
312  type(cosp_isccp),intent(inout)   :: isccp   ! Output from ISCCP simulator
313  type(cosp_misr),intent(inout)    :: misr    ! Output from MISR simulator
314  type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
315  type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
316
317  ! Local variables
318  integer :: Npoints   ! Number of gridpoints
319  integer :: Ncolumns  ! Number of subcolumns
320  integer :: Nlevels   ! Number of levels
321  integer :: Nhydro    ! Number of hydrometeors
322  integer :: Niter     ! Number of calls to cosp_simulator
323  integer :: i,j,k
324  real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
325  integer,parameter :: scops_debug=0    !  set to non-zero value to print out inputs for debugging in SCOPS
326 
327  real,dimension(:, :),allocatable :: cca_scops,ls_p_rate,cv_p_rate, &
328                     tca_scops ! Cloud cover in each model level (HORIZONTAL gridbox fraction) of total cloud.
329                               ! Levels are from TOA to SURFACE. (nPoints, nLev)
330  real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
331                                                                     ! Levels are from SURFACE to TOA
332  type(cosp_sghydro) :: sghydro   ! Subgrid info for hydrometeors en each iteration
333
334 
335  !++++++++++ Dimensions ++++++++++++
336  Npoints  = gbx%Npoints
337  Ncolumns = gbx%Ncolumns
338  Nlevels  = gbx%Nlevels
339  Nhydro   = gbx%Nhydro
340   
341   
342  !++++++++++ Climate/NWP mode ++++++++++ 
343  if (Ncolumns > 1) then
344        !++++++++++ Subgrid sampling ++++++++++
345        ! Allocate arrays before calling SCOPS
346        allocate(frac_ls(Npoints,Nlevels),frac_cv(Npoints,Nlevels),prec_ls(Npoints,Nlevels),prec_cv(Npoints,Nlevels))
347        allocate(tca_scops(Npoints,Nlevels),cca_scops(Npoints,Nlevels), &
348                ls_p_rate(Npoints,Nlevels),cv_p_rate(Npoints,Nlevels))
349        ! Initialize to zero
350        frac_ls=0.0
351        prec_ls=0.0
352        frac_cv=0.0
353        prec_cv=0.0
354        ! Cloud fractions for SCOPS from TOA to SFC
355        tca_scops = gbx%tca(:,Nlevels:1:-1)
356        cca_scops = gbx%cca(:,Nlevels:1:-1)
357       
358        ! Call to SCOPS
359        ! strat and conv arrays are passed with levels from TOA to SURFACE.
360        call scops(Npoints,Nlevels,Ncolumns,seed,tca_scops,cca_scops,overlap,sgx%frac_out,scops_debug)
361       
362        ! temporarily use prec_ls/cv to transfer information about precipitation flux into prec_scops
363        if(gbx%use_precipitation_fluxes) then
364            ls_p_rate(:,Nlevels:1:-1)=gbx%rain_ls(:,1:Nlevels)+gbx%snow_ls(:,1:Nlevels)+gbx%grpl_ls(:,1:Nlevels)
365            cv_p_rate(:,Nlevels:1:-1)=gbx%rain_cv(:,1:Nlevels)+gbx%snow_cv(:,1:Nlevels)
366        else
367            ls_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_LSRAIN)+ &
368                                      gbx%mr_hydro(:,1:Nlevels,I_LSSNOW)+ &
369                                      gbx%mr_hydro(:,1:Nlevels,I_LSGRPL)
370            cv_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_CVRAIN)+ &
371                                      gbx%mr_hydro(:,1:Nlevels,I_CVSNOW)
372        endif
373       
374        call prec_scops(Npoints,Nlevels,Ncolumns,ls_p_rate,cv_p_rate,sgx%frac_out,sgx%prec_frac)
375       
376        ! Precipitation fraction
377        do j=1,Npoints,1
378        do k=1,Nlevels,1
379            do i=1,Ncolumns,1
380                if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 1) frac_ls(j,k)=frac_ls(j,k)+1.
381                if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 2) frac_cv(j,k)=frac_cv(j,k)+1.
382                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
383                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
384                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 3) then
385                    prec_cv(j,k)=prec_cv(j,k)+1.
386                    prec_ls(j,k)=prec_ls(j,k)+1.
387                endif
388            enddo  !i
389            frac_ls(j,k)=frac_ls(j,k)/Ncolumns
390            frac_cv(j,k)=frac_cv(j,k)/Ncolumns
391            prec_ls(j,k)=prec_ls(j,k)/Ncolumns
392            prec_cv(j,k)=prec_cv(j,k)/Ncolumns
393        enddo  !k
394        enddo  !j
395       
396         ! Levels from SURFACE to TOA.
397        if (Npoints*Ncolumns*Nlevels < 10000) then
398            sgx%frac_out(1:Npoints,:,1:Nlevels)  = sgx%frac_out(1:Npoints,:,Nlevels:1:-1)
399            sgx%prec_frac(1:Npoints,:,1:Nlevels) = sgx%prec_frac(1:Npoints,:,Nlevels:1:-1)
400        else
401            ! This is done within a loop (unvectorized) over nPoints to save memory
402            do j=1,Npoints
403                sgx%frac_out(j,:,1:Nlevels)  = sgx%frac_out(j,:,Nlevels:1:-1)
404                sgx%prec_frac(j,:,1:Nlevels) = sgx%prec_frac(j,:,Nlevels:1:-1)
405            enddo
406        endif
407       
408       ! Deallocate arrays that will no longer be used
409        deallocate(tca_scops,cca_scops,ls_p_rate,cv_p_rate)
410         
411        ! Populate the subgrid arrays
412        call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
413        do k=1,Ncolumns
414            !--------- Mixing ratios for clouds and Reff for Clouds and precip -------
415            column_frac_out => sgx%frac_out(:,k,:)
416            where (column_frac_out == 1)     !+++++++++++ LS clouds ++++++++
417                sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
418                sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
419               
420                sghydro%Reff(:,k,:,I_LSCLIQ)     = gbx%Reff(:,:,I_LSCLIQ)
421                sghydro%Reff(:,k,:,I_LSCICE)     = gbx%Reff(:,:,I_LSCICE)
422                sghydro%Reff(:,k,:,I_LSRAIN)     = gbx%Reff(:,:,I_LSRAIN)
423                sghydro%Reff(:,k,:,I_LSSNOW)     = gbx%Reff(:,:,I_LSSNOW)
424                sghydro%Reff(:,k,:,I_LSGRPL)     = gbx%Reff(:,:,I_LSGRPL)
425            elsewhere (column_frac_out == 2) !+++++++++++ CONV clouds ++++++++
426                sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ)
427                sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE)
428               
429                sghydro%Reff(:,k,:,I_CVCLIQ)     = gbx%Reff(:,:,I_CVCLIQ)
430                sghydro%Reff(:,k,:,I_CVCICE)     = gbx%Reff(:,:,I_CVCICE)
431                sghydro%Reff(:,k,:,I_CVRAIN)     = gbx%Reff(:,:,I_CVRAIN)
432                sghydro%Reff(:,k,:,I_CVSNOW)     = gbx%Reff(:,:,I_CVSNOW)
433            end where
434            !--------- Precip -------
435            if (.not. gbx%use_precipitation_fluxes) then
436                where (column_frac_out == 1)  !+++++++++++ LS Precipitation ++++++++
437                    sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN)
438                    sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW)
439                    sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL)
440                elsewhere (column_frac_out == 2) !+++++++++++ CONV Precipitation ++++++++
441                    sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN)
442                    sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW)
443                end where
444            endif
445        enddo
446        ! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values
447        do k=1,Nlevels
448            do j=1,Npoints
449                !--------- Clouds -------
450                if (frac_ls(j,k) .ne. 0.) then
451                    sghydro%mr_hydro(j,:,k,I_LSCLIQ) = sghydro%mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
452                    sghydro%mr_hydro(j,:,k,I_LSCICE) = sghydro%mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
453                endif
454                if (frac_cv(j,k) .ne. 0.) then
455                    sghydro%mr_hydro(j,:,k,I_CVCLIQ) = sghydro%mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
456                    sghydro%mr_hydro(j,:,k,I_CVCICE) = sghydro%mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
457                endif
458                !--------- Precip -------
459                if (gbx%use_precipitation_fluxes) then
460                    if (prec_ls(j,k) .ne. 0.) then
461                        gbx%rain_ls(j,k) = gbx%rain_ls(j,k)/prec_ls(j,k)
462                        gbx%snow_ls(j,k) = gbx%snow_ls(j,k)/prec_ls(j,k)
463                        gbx%grpl_ls(j,k) = gbx%grpl_ls(j,k)/prec_ls(j,k)
464                    endif
465                    if (prec_cv(j,k) .ne. 0.) then
466                        gbx%rain_cv(j,k) = gbx%rain_cv(j,k)/prec_cv(j,k)
467                        gbx%snow_cv(j,k) = gbx%snow_cv(j,k)/prec_cv(j,k)
468                    endif
469                else
470                    if (prec_ls(j,k) .ne. 0.) then
471                        sghydro%mr_hydro(j,:,k,I_LSRAIN) = sghydro%mr_hydro(j,:,k,I_LSRAIN)/prec_ls(j,k)
472                        sghydro%mr_hydro(j,:,k,I_LSSNOW) = sghydro%mr_hydro(j,:,k,I_LSSNOW)/prec_ls(j,k)
473                        sghydro%mr_hydro(j,:,k,I_LSGRPL) = sghydro%mr_hydro(j,:,k,I_LSGRPL)/prec_ls(j,k)
474                    endif
475                    if (prec_cv(j,k) .ne. 0.) then
476                        sghydro%mr_hydro(j,:,k,I_CVRAIN) = sghydro%mr_hydro(j,:,k,I_CVRAIN)/prec_cv(j,k)
477                        sghydro%mr_hydro(j,:,k,I_CVSNOW) = sghydro%mr_hydro(j,:,k,I_CVSNOW)/prec_cv(j,k)
478                    endif
479                endif 
480            enddo !k
481        enddo !j
482        deallocate(frac_ls,prec_ls,frac_cv,prec_cv)
483       
484        if (gbx%use_precipitation_fluxes) then
485            ! convert precipitation flux into mixing ratio
486            call pf_to_mr(Npoints,Nlevels,Ncolumns,gbx%rain_ls,gbx%snow_ls,gbx%grpl_ls, &
487                        gbx%rain_cv,gbx%snow_cv,sgx%prec_frac,gbx%p,gbx%T, &
488                        sghydro%mr_hydro(:,:,:,I_LSRAIN),sghydro%mr_hydro(:,:,:,I_LSSNOW),sghydro%mr_hydro(:,:,:,I_LSGRPL), &
489                        sghydro%mr_hydro(:,:,:,I_CVRAIN),sghydro%mr_hydro(:,:,:,I_CVSNOW))
490        endif
491   !++++++++++ CRM mode ++++++++++
492   else
493      sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
494      sghydro%Reff(:,1,:,:) = gbx%Reff
495      !--------- Clouds -------
496      where ((gbx%dtau_s > 0.0))
497             sgx%frac_out(:,1,:) = 1  ! Subgrid cloud array. Dimensions (Npoints,Ncolumns,Nlevels)
498      endwhere
499   endif ! Ncolumns > 1
500 
501   
502   !++++++++++ Simulator ++++++++++
503    call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,stradar,stlidar)
504
505    ! Deallocate subgrid arrays
506    call free_cosp_sghydro(sghydro)
507END SUBROUTINE COSP_ITER
508
509END MODULE MOD_COSP
Note: See TracBrowser for help on using the repository browser.