Context Navigation

cosp.F90 @ 1280

Last change on this file since 1280 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
25	MODULE MOD_COSP
26	USE MOD_COSP_TYPES
27	USE MOD_COSP_SIMULATOR
28	IMPLICIT NONE
29
30	CONTAINS
31
32
33	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34	!--------------------- SUBROUTINE COSP ---------------------------
35	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36	SUBROUTINE 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
278	END SUBROUTINE COSP
279
280	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281	!--------------------- SUBROUTINE COSP_ITER ----------------------
282	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283	SUBROUTINE 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 (NpointsNcolumnsNlevels < 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)
488	END SUBROUTINE COSP_ITER
489
490	END MODULE MOD_COSP

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format