source: trunk/libf/cosp/cosp.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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