Context Navigation

cosp.F90 @ 5048

Last change on this file since 5048 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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

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

Download in other formats:

Original Format