source: LMDZ4/branches/LMDZ4-dev/libf/cosp/cosp.F90 @ 1269

Last change on this file since 1269 was 1262, checked in by idelkadi, 15 years ago

-Rajout des fichiers cosp_input_nl.txt et cosp_output_nl.txt pour controler les parametres entrees sorties COSP
-Rajout du repertoire cosp contenant les sources de COSP

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