Context Navigation

quickbeam.F90 @ 5209

Last change on this file since 5209 was 5185, checked in by abarral, 2 months ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 31.1 KB

Line
1	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2	! Copyright (c) 2015, Regents of the University of Colorado
3	! All rights reserved.
4
5	! Redistribution and use in source and binary forms, with or without modification, are
6	! permitted provided that the following conditions are met:
7
8	! 1. Redistributions of source code must retain the above copyright notice, this list of
9	! conditions and the following disclaimer.
10
11	! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12	! of conditions and the following disclaimer in the documentation and/or other
13	! materials provided with the distribution.
14
15	! 3. Neither the name of the copyright holder nor the names of its contributors may be
16	! used to endorse or promote products derived from this software without specific prior
17	! written permission.
18
19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20	! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21	! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22	! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24	! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26	! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28
29	! History
30	! 11/2005: John Haynes - Created
31	! 09/2006 placed into subroutine form (Roger Marchand,JMH)
32	! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
33	! 01/2008 'Do while' to determine if hydrometeor(s) present in volume
34	! changed for vectorization purposes (A. Bodas-Salcedo)
35
36	! 07/2010 V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT)
37	! ... All hydrometeor and radar simulator properties now included in hp structure
38	! ... hp structure should be initialized by call to radar_simulator_init prior
39	! ... to calling this subroutine.
40	! Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added
41	! ... Changes implement by Roj Marchand following work by Laura Fowler
42
43	! 10/2011 Modified ngate loop to go in either direction depending on flag
44	! hp%radar_at_layer_one. This affects the direction in which attenuation is summed.
45
46	! Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple
47	! summation. (Roger Marchand)
48	! May 2015 - D. Swales - Modified for COSPv2.0
49	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50	module quickbeam
51	USE COSP_KINDS, ONLY: wp
52	USE MOD_COSP_CONFIG, ONLY: R_UNDEF,cloudsat_histRef,use_vgrid,vgrid_zl,vgrid_zu,&
53	pClass_noPrecip, pClass_Rain1, pClass_Rain2, pClass_Rain3,&
54	pClass_Snow1, pClass_Snow2, pClass_Mixed1, pClass_Mixed2, &
55	pClass_Rain4, pClass_default, Zenonbinval, Zbinvallnd, &
56	N_HYDRO,nCloudsatPrecipClass !,cloudsat_preclvl !PREC_BUG
57	USE MOD_COSP_STATS, ONLY: COSP_LIDAR_ONLY_CLOUD,hist1D,COSP_CHANGE_VERTICAL_GRID
58	implicit none
59
60	integer,parameter :: &
61	maxhclass = 20, & ! Qucikbeam maximum number of hydrometeor classes.
62	nRe_types = 550, & ! Quickbeam maximum number or Re size bins allowed in N and Z_scaled look up table.
63	nd = 85, & ! Qucikbeam number of discrete particles used in construction DSDs.
64	mt_ntt = 39, & ! Quickbeam number of temperatures in mie LUT.
65	Re_BIN_LENGTH = 10, & ! Quickbeam minimum Re interval in scale LUTs
66	Re_MAX_BIN = 250 ! Quickbeam maximum Re interval in scale LUTs
67	real(wp),parameter :: &
68	dmin = 0.1, & ! Quickbeam minimum size of discrete particle
69	dmax = 10000. ! Quickbeam maximum size of discrete particle
70
71	!djs logical :: radar_at_layer_one ! If true radar is assume to be at the edge
72	! of the first layer, if the first layer is the
73	! surface than a ground-based radar. If the
74	! first layer is the top-of-atmosphere, then
75	! a space borne radar.
76
77	! ##############################################################################################
78	type radar_cfg
79	! Radar properties
80	real(wp) :: freq,k2
81	integer :: nhclass ! Number of hydrometeor classes in use
82	integer :: use_gas_abs, do_ray
83	logical :: radar_at_layer_one ! If true radar is assume to be at the edge
84	! of the first layer, if the first layer is the
85	! surface than a ground-based radar. If the
86	! first layer is the top-of-atmosphere, then
87	! a space borne radar.
88
89	! Variables used to store Z scale factors
90	character(len=240) :: scale_LUT_file_name
91	logical :: load_scale_LUTs, update_scale_LUTs
92	logical, dimension(maxhclass,nRe_types) :: N_scale_flag
93	logical, dimension(maxhclass,mt_ntt,nRe_types) :: Z_scale_flag,Z_scale_added_flag
94	real(wp),dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled
95	real(wp),dimension(maxhclass,nd,nRe_types) :: fc, rho_eff
96	real(wp),dimension(Re_MAX_BIN) :: base_list,step_list
97
98	end type radar_cfg
99
100	contains
101	! ######################################################################################
102	! SUBROUTINE quickbeam_subcolumn
103	! ######################################################################################
104	subroutine quickbeam_subcolumn(rcfg,nprof,ngate,hgt_matrix,z_vol,kr_vol,g_vol,dBZe,Ze_non)
105
106	! INPUTS
107	type(radar_cfg),intent(inout) :: &
108	rcfg ! Derived type for radar simulator setup
109	integer,intent(in) :: &
110	nprof, & ! Number of hydrometeor profiles
111	ngate ! Number of vertical layers
112	real(wp),intent(in),dimension(nprof,ngate) :: &
113	hgt_matrix, & ! Height of hydrometeors (km)
114	z_vol, & ! Effective reflectivity factor (mm^6/m^3)
115	kr_vol, & ! Attenuation coefficient hydro (dB/km)
116	g_vol ! Attenuation coefficient gases (dB/km)
117
118	! OUTPUTS
119	real(wp), intent(out),dimension(nprof,ngate) :: &
120	Ze_non, & ! Radar reflectivity without attenuation (dBZ)
121	dBZe ! Effective radar reflectivity factor (dBZ)
122
123	! LOCAL VARIABLES
124	integer :: k,pr,start_gate,end_gate,d_gate
125	real(wp),dimension(nprof,ngate) :: &
126	Ze_ray, & ! Rayleigh reflectivity (dBZ)
127	g_to_vol, & ! Gaseous atteunation, radar to vol (dB)
128	a_to_vol, & ! Hydromets attenuation, radar to vol (dB)
129	z_ray ! Reflectivity factor, Rayleigh only (mm^6/m^3)
130
131	! Load scaling matricies from disk -- but only the first time this subroutine is called
132	if(rcfg%load_scale_LUTs) then
133	call load_scale_LUTs(rcfg)
134	rcfg%load_scale_LUTs=.false.
135	rcfg%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
136	endif
137
138	! Initialization
139	g_to_vol = 0._wp
140	a_to_vol = 0._wp
141
142	! Loop over each range gate (ngate) ... starting with layer closest to the radar !
143	if(rcfg%radar_at_layer_one) then
144	start_gate = 1
145	end_gate = ngate
146	d_gate = 1
147	else
148	start_gate = ngate
149	end_gate = 1
150	d_gate = -1
151	endif
152	DO k=start_gate,end_gate,d_gate
153	! Loop over each profile (nprof)
154	DO pr=1,nprof
155	! Attenuation due to hydrometeors between radar and volume
156
157	! NOTE old scheme integrates attenuation only for the layers ABOVE
158	! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
159	! which may be a problem. ROJ
160	! in the new scheme I assign half the attenuation to the current layer
161	if(d_gate==1) then
162	! dheight calcuations assumes hgt_matrix points are the cell mid-points.
163	if (k>2) then
164	! add to previous value to half of above layer + half of current layer
165	a_to_vol(pr,k)= a_to_vol(pr,k-1) + &
166	(kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
167	else
168	a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
169	endif
170	else ! d_gate==-1
171	if(k<ngate) then
172	! Add to previous value half of above layer + half of current layer
173	a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
174	(kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
175	else
176	a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
177	endif
178	endif
179
180	! Attenuation due to gaseous absorption between radar and volume
181	if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .AND. pr .eq. 1)) then
182	if (d_gate==1) then
183	if (k>1) then
184	! Add to previous value to half of above layer + half of current layer
185	g_to_vol(pr,k) = g_to_vol(pr,k-1) + &
186	0.5(g_vol(pr,k-1)+g_vol(pr,k))(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
187	else
188	g_to_vol(pr,k)= 0.5_wpg_vol(pr,k)(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
189	endif
190	else ! d_gate==-1
191	if (k<ngate) then
192	! Add to previous value to half of above layer + half of current layer
193	g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
194	0.5_wp(g_vol(pr,k+1)+g_vol(pr,k))(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
195	else
196	g_to_vol(pr,k)= 0.5_wpg_vol(pr,k)(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
197	endif
198	endif
199	elseif(rcfg%use_gas_abs == 2) then
200	! Using value calculated for the first column
201	g_to_vol(pr,k) = g_to_vol(1,k)
202	elseif (rcfg%use_gas_abs == 0) then
203	g_to_vol(pr,k) = 0._wp
204	endif
205	enddo ! End loop over pr (profile)
206	enddo ! End loop of k (range gate)
207
208	! Compute Rayleigh reflectivity, and full, attenuated reflectivity
209	if(rcfg%do_ray == 1) then
210	where(z_ray(1:nprof,1:ngate) > 0._wp)
211	Ze_ray(1:nprof,1:ngate) = 10._wp*log10(z_ray(1:nprof,1:ngate))
212	elsewhere
213	Ze_Ray(1:nprof,1:ngate) = 0._wp
214	endwhere
215	else
216	Ze_ray(1:nprof,1:ngate) = R_UNDEF
217	end if
218
219	where(z_vol(1:nprof,1:ngate) > 0._wp)
220	Ze_non(1:nprof,1:ngate) = 10._wp*log10(z_vol(1:nprof,1:ngate))
221	dBZe(1:nprof,1:ngate) = Ze_non(1:nprof,1:ngate)-a_to_vol(1:nprof,1:ngate)-g_to_vol(1:nprof,1:ngate)
222	elsewhere
223	dBZe(1:nprof,1:ngate) = R_UNDEF
224	Ze_non(1:nprof,1:ngate) = R_UNDEF
225	end where
226
227	! Save any updates made
228	if (rcfg%update_scale_LUTs) call save_scale_LUTs(rcfg)
229
230	end subroutine quickbeam_subcolumn
231	! ######################################################################################
232	! SUBROUTINE quickbeam_column
233	! ######################################################################################
234	subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform, &
235	Ze_tot, Ze_tot_non, land, surfelev, t2m, fracPrecipIce, zlev, zlev_half, cfad_ze, & !PREC_BUG
236	cloudsat_precip_cover, cloudsat_pia)
237	! Inputs
238	integer,intent(in) :: &
239	npoints, & ! Number of horizontal grid points
240	ncolumns, & ! Number of subcolumns
241	nlevels, & ! Number of vertical layers in OLD grid
242	llm, & ! Number of vertical layers in NEW grid
243	DBZE_BINS ! Number of bins for cfad.
244	character(len=*),intent(in) :: &
245	platform ! Name of platform (e.g. cloudsat)
246	real(wp),dimension(Npoints),intent(in) :: &
247	land, & ! Land/Sea mask. (1/0)
248	surfelev, & ! Surface Elevation (m) !PREC_BUG
249	t2m ! Near-surface temperature
250	real(wp),dimension(Npoints,Ncolumns),intent(in) :: &
251	fracPrecipIce ! Fraction of precipitation which is frozen. (1)
252	real(wp),intent(in),dimension(npoints,ncolumns,Nlevels) :: &
253	Ze_tot, & ! Effective reflectivity factor (dBZ)
254	Ze_tot_non ! Effective reflectivity factor w/o attenuation (dBZ)
255	real(wp),intent(in),dimension(npoints,Nlevels) :: &
256	zlev ! Model full levels
257	real(wp),intent(in),dimension(npoints,Nlevels+1) :: &
258	zlev_half ! Model half levels
259
260	! Outputs
261	real(wp),intent(inout),dimension(npoints,DBZE_BINS,llm) :: &
262	cfad_ze !
263	real(wp),dimension(Npoints,nCloudsatPrecipClass),intent(out) :: &
264	cloudsat_precip_cover ! Model precip rate in by CloudSat precip flag
265	real(wp),dimension(Npoints),intent(out) :: &
266	cloudsat_pia ! Cloudsat path integrated attenuation
267
268	! Local variables
269	integer :: i,j
270	real(wp),dimension(npoints,ncolumns,llm) :: ze_toti,ze_noni
271	logical :: lcloudsat = .false.
272
273	! Which platforms to create diagnostics for?
274	if (platform .eq. 'cloudsat') lcloudsat=.true.
275
276	! Create Cloudsat diagnostics.
277	if (lcloudsat) then
278	if (use_vgrid) then
279	! Regrid in the vertical (NOTE This routine requires SFC-2-TOA ordering, so flip
280	! inputs and outputs to maintain TOA-2-SFC ordering convention in COSP2.)
281	call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),&
282	zlev_half(:,nlevels:1:-1),Ze_tot(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),&
283	vgrid_zu(llm:1:-1),Ze_toti(:,:,llm:1:-1),log_units=.true.)
284
285	! Effective reflectivity histogram
286	DO i=1,Npoints
287	DO j=1,llm
288	cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_toti(i,:,j),DBZE_BINS,cloudsat_histRef)
289	enddo
290	enddo
291	where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
292
293	! Compute cloudsat near-surface precipitation diagnostics
294	! First, regrid in the vertical Ze_tot_non.
295	call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),&
296	zlev_half(:,nlevels:1:-1),Ze_tot_non(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),&
297	vgrid_zu(llm:1:-1),Ze_noni(:,:,llm:1:-1),log_units=.true.)
298	! Not call routine to generate diagnostics.
299	call cloudsat_precipOccurence(Npoints, Ncolumns, llm, N_HYDRO, Ze_toti, Ze_noni, &
300	land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia) !PREC_BUG
301	else
302	! Effective reflectivity histogram
303	DO i=1,Npoints
304	DO j=1,llm
305	cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_tot(i,:,j),DBZE_BINS,cloudsat_histRef)
306	enddo
307	enddo
308	where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
309	endif
310	endif
311
312	end subroutine quickbeam_column
313	! ##############################################################################################
314	! ##############################################################################################
315	! ######################################################################################
316	! SUBROUTINE cloudsat_precipOccurence
317
318	! Notes from July 2016: Add precip flag also looped over subcolumns
319	! Modified by Tristan L'Ecuyer (TSL) to add precipitation flagging
320	! based on CloudSat's 2C-PRECIP-COLUMN algorithm (Haynes et al, JGR, 2009).
321	! To mimic the satellite algorithm, this code applies thresholds to non-attenuated
322	! reflectivities, Ze_non, consistent with those outlined in Haynes et al, JGR (2009).
323
324	! Procedures/Notes:
325
326	! (1) If the 2-way attenuation exceeds 40 dB, the pixel will be flagged as 'heavy rain'
327	! consistent with the multiple-scattering analysis of Battaglia et al, JGR (2008).
328	! (2) Rain, snow, and mixed precipitation scenes are separated according to the fraction
329	! of the total precipitation hydrometeor mixing ratio that exists as ice.
330	! (3) The entire analysis is applied to the range gate from 480-960 m to be consistent with
331	! CloudSat's 720 m ground-clutter.
332	! (4) Only a single flag is assigned to each model grid since there is no variation in
333	! hydrometeor contents across a single model level. Unlike CFADs, whose variation enters
334	! due to differening attenuation corrections from hydrometeors aloft, the non-attenuated
335	! reflectivities used in the computation of this flag cannot vary across sub-columns.
336
337	! radar_prec_flag = 1-Rain possible 2-Rain probable 3-Rain certain
338	! 4-Snow possible 5-Snow certain
339	! 6-Mixed possible 7-Mixed certain
340	! 8-Heavy Rain
341	! 9- default value
342
343	! Modified by Dustin Swales (University of Colorado) for use with COSP2.
344	! NOTE All inputs (Ze_out, Ze_non_out, fracPrecipIce) are at a single level from the
345	! statistical output grid used by Cloudsat. This level is controlled by the
346	! parameter cloudsat_preclvl, defined in src/cosp_config.F90
347	! ######################################################################################
348	subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_non_out, &
349	land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia) !PREC_BUG
350
351	! Inputs
352	integer,intent(in) :: &
353	Npoints, & ! Number of columns
354	Ncolumns, & ! Numner of subcolumns
355	Nhydro, & ! Number of hydrometeor types
356	llm ! Number of levels
357	real(wp),dimension(Npoints),intent(in) :: &
358	land, & ! Land/Sea mask. (1/0)
359	surfelev, & ! Surface Elevation (m) !PREC_BUG
360	t2m ! Near-surface temperature
361	real(wp),dimension(Npoints,Ncolumns,llm),intent(in) :: &
362	Ze_out, & ! Effective reflectivity factor (dBZ)
363	Ze_non_out ! Effective reflectivity factor, w/o attenuation (dBZ)
364	real(wp),dimension(Npoints,Ncolumns),intent(in) :: &
365	fracPrecipIce ! Fraction of precipitation which is frozen. (1)
366
367	! Outputs
368	real(wp),dimension(Npoints,nCloudsatPrecipClass),intent(out) :: &
369	cloudsat_precip_cover ! Model precip rate in by CloudSat precip flag
370	real(wp),dimension(Npoints),intent(out) :: &
371	cloudsat_pia ! Cloudsat path integrated attenuation
372
373	! Local variables
374	integer,dimension(Npoints,Ncolumns) :: &
375	cloudsat_pflag, & ! Subcolumn precipitation flag
376	cloudsat_precip_pia ! Subcolumn path integrated attenutation.
377	integer,dimension(Npoints) :: & !PREC_BUG
378	cloudsat_preclvl_index ! Altitude index for precip flags calculation !PREC_BUG
379	! in 40-level grid (2nd layer above surfelev) !PREC_BUG
380	integer :: pr,i,k,m,j,cloudsat_preclvl !PREC_BUG
381	real(wp) :: Zmax
382
383	! Initialize
384	cloudsat_pflag(:,:) = pClass_default
385	cloudsat_precip_pia(:,:) = 0._wp
386	cloudsat_precip_cover(:,:) = 0._wp
387	cloudsat_pia(:) = 0._wp
388	cloudsat_preclvl_index(:) = 0._wp !PREC_BUG
389
390
391	!!! Computing altitude index for precip flags calculation !PREC_BUG
392	cloudsat_preclvl_index(:) = 39 - floor( surfelev(:)/480. ) !PREC_BUG
393
394	! ######################################################################################
395	! SUBCOLUMN processing
396	! ######################################################################################
397	DO i=1, Npoints
398
399	cloudsat_preclvl = cloudsat_preclvl_index(i) !PREC_BUG
400	! print*, 'i, surfelev(i), cloudsat_preclvl : ', i, surfelev(i), cloudsat_preclvl !PREC_BUG
401
402	DO pr=1,Ncolumns
403	! 1) Compute the PIA in all profiles containing hydrometeors
404	if ( (Ze_non_out(i,pr,cloudsat_preclvl).gt.-100) .AND. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then
405	if ( (Ze_non_out(i,pr,cloudsat_preclvl).lt.100) .AND. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then
406	cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl)
407	endif
408	endif
409
410	! 2) Compute precipitation flag
411	! ################################################################################
412	! 2a) Oceanic points.
413	! ################################################################################
414	if (land(i) .eq. 0) then
415	! print*, 'aaa i, pr, fracPrecipIce(i,pr) : ', i, pr, fracPrecipIce(i,pr) !Artem
416	! Snow
417	if(fracPrecipIce(i,pr).gt.0.9) then
418	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
419	cloudsat_pflag(i,pr) = pClass_Snow2 ! TSL: Snow certain
420	endif
421	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).AND. &
422	Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
423	cloudsat_pflag(i,pr) = pClass_Snow1 ! TSL: Snow possible
424	endif
425	endif
426
427	! Mixed
428	if(fracPrecipIce(i,pr).gt.0.1.AND.fracPrecipIce(i,pr).le.0.9) then
429	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
430	cloudsat_pflag(i,pr) = pClass_Mixed2 ! TSL: Mixed certain
431	endif
432	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).AND. &
433	Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
434	cloudsat_pflag(i,pr) = pClass_Mixed1 ! TSL: Mixed possible
435	endif
436	endif
437
438	! Rain
439	if(fracPrecipIce(i,pr).le.0.1) then
440	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(1)) then
441	cloudsat_pflag(i,pr) = pClass_Rain3 ! TSL: Rain certain
442	endif
443	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(3).AND. &
444	Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(1)) then
445	cloudsat_pflag(i,pr) = pClass_Rain2 ! TSL: Rain probable
446	endif
447	if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).AND. &
448	Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(3)) then
449	cloudsat_pflag(i,pr) = pClass_Rain1 ! TSL: Rain possible
450	endif
451	if(cloudsat_precip_pia(i,pr).gt.40) then
452	cloudsat_pflag(i,pr) = pClass_Rain4 ! TSL: Heavy Rain
453	endif
454	endif
455
456	! No precipitation
457	if(Ze_non_out(i,pr,cloudsat_preclvl).le.-15) then
458	cloudsat_pflag(i,pr) = pClass_noPrecip ! TSL: Not Raining
459	endif
460	endif ! Ocean points
461
462	! ################################################################################
463	! 2b) Land points.
464	! ################################################################################
465	if (land(i) .eq. 1) then
466	! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out);
467	Zmax=maxval(Ze_out(i,pr,:))
468
469	! Snow (T<273)
470	if(t2m(i) .lt. 273._wp) then
471	if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(5)) then
472	cloudsat_pflag(i,pr) = pClass_Snow2 ! JEK: Snow certain
473	endif
474	if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .AND. &
475	Ze_out(i,pr,cloudsat_preclvl).le.Zbinvallnd(5)) then
476	cloudsat_pflag(i,pr) = pClass_Snow1 ! JEK: Snow possible
477	endif
478	endif
479
480	! Mized phase (273<T<275)
481	if(t2m(i) .ge. 273._wp .AND. t2m(i) .le. 275._wp) then
482	if ((Zmax .gt. Zbinvallnd(1) .AND. cloudsat_precip_pia(i,pr).gt.30) .or. &
483	(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(4))) then
484	cloudsat_pflag(i,pr) = pClass_Mixed2 ! JEK: Mixed certain
485	endif
486	if ((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .AND. &
487	Ze_out(i,pr,cloudsat_preclvl) .le. Zbinvallnd(4)) .AND. &
488	(Zmax .gt. Zbinvallnd(5)) ) then
489	cloudsat_pflag(i,pr) = pClass_Mixed1 ! JEK: Mixed possible
490	endif
491	endif
492
493	! Rain (T>275)
494	if(t2m(i) .gt. 275) then
495	if ((Zmax .gt. Zbinvallnd(1) .AND. cloudsat_precip_pia(i,pr).gt.30) .or. &
496	(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(2))) then
497	cloudsat_pflag(i,pr) = pClass_Rain3 ! JEK: Rain certain
498	endif
499	if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .AND. &
500	(Zmax .gt. Zbinvallnd(3))) then
501	cloudsat_pflag(i,pr) = pClass_Rain2 ! JEK: Rain probable
502	endif
503	if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .AND. &
504	(Zmax.lt.Zbinvallnd(3))) then
505	cloudsat_pflag(i,pr) = pClass_Rain1 ! JEK: Rain possible
506	endif
507	if(cloudsat_precip_pia(i,pr).gt.40) then
508	cloudsat_pflag(i,pr) = pClass_Rain4 ! JEK: Heavy Rain
509	endif
510	endif
511
512	! No precipitation
513	if(Ze_out(i,pr,cloudsat_preclvl).le.-15) then
514	cloudsat_pflag(i,pr) = pClass_noPrecip ! JEK: Not Precipitating
515	endif
516	endif ! Land points
517	enddo ! Sub-columns
518	enddo ! Gridpoints
519
520	! ######################################################################################
521	! COLUMN processing
522	! ######################################################################################
523
524	! Aggregate subcolumns
525	DO i=1,Npoints
526	! Gridmean precipitation fraction for each precipitation type
527	DO k=1,nCloudsatPrecipClass
528	if (any(cloudsat_pflag(i,:) .eq. k-1)) then
529	cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq. k-1)
530	endif
531	enddo
532
533	! Gridmean path integrated attenuation (pia)
534	cloudsat_pia(i)=sum(cloudsat_precip_pia(i,:))
535	enddo
536
537	! Normalize by number of subcolumns
538	where ((cloudsat_precip_cover /= R_UNDEF).AND.(cloudsat_precip_cover /= 0.0)) &
539	cloudsat_precip_cover = cloudsat_precip_cover / Ncolumns
540	where ((cloudsat_pia/= R_UNDEF).AND.(cloudsat_pia/= 0.0)) &
541	cloudsat_pia = cloudsat_pia / Ncolumns
542
543	end subroutine cloudsat_precipOccurence
544
545	! ##############################################################################################
546	! ##############################################################################################
547	subroutine load_scale_LUTs(rcfg)
548
549	type(radar_cfg), intent(inout) :: rcfg
550	logical :: LUT_file_exists
551	integer :: i,j,k,ind
552
553	! Load scale LUT from file
554	inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
555	exist=LUT_file_exists)
556
557	if(.not.LUT_file_exists) then
558	write(,) '*************************************************'
559	write(,) 'Warning: Could NOT FIND radar LUT file: ', &
560	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
561	write(,) 'Will calculated LUT values as needed'
562	write(,) '*************************************************'
563	return
564	else
565	OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
566	form='unformatted', &
567	err= 89, &
568	access='DIRECT',&
569	recl=28)
570	write(,) 'Loading radar LUT file: ', &
571	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
572
573	DO i=1,maxhclass
574	DO j=1,mt_ntt
575	DO k=1,nRe_types
576	ind = i+(j-1)maxhclass+(k-1)(nRe_types*mt_ntt)
577	read(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
578	rcfg%Ze_scaled(i,j,k), &
579	rcfg%Zr_scaled(i,j,k), &
580	rcfg%kr_scaled(i,j,k)
581	enddo
582	enddo
583	enddo
584	close(unit=12)
585	return
586	endif
587
588	89 write(,) 'Error: Found but could NOT READ radar LUT file: ', &
589	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
590
591	end subroutine load_scale_LUTs
592
593	! ##############################################################################################
594	! ##############################################################################################
595	subroutine save_scale_LUTs(rcfg)
596	type(radar_cfg), intent(inout) :: rcfg
597	logical :: LUT_file_exists
598	integer :: i,j,k,ind
599
600	inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
601	exist=LUT_file_exists)
602
603	OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
604	form='unformatted',err= 99,access='DIRECT',recl=28)
605
606	write(,) 'Creating or Updating radar LUT file: ', &
607	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
608
609	DO i=1,maxhclass
610	DO j=1,mt_ntt
611	DO k=1,nRe_types
612	ind = i+(j-1)maxhclass+(k-1)(nRe_types*mt_ntt)
613	if(.not.LUT_file_exists .or. rcfg%Z_scale_added_flag(i,j,k)) then
614	rcfg%Z_scale_added_flag(i,j,k)=.false.
615	write(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
616	rcfg%Ze_scaled(i,j,k), &
617	rcfg%Zr_scaled(i,j,k), &
618	rcfg%kr_scaled(i,j,k)
619	endif
620	enddo
621	enddo
622	enddo
623	close(unit=12)
624	return
625
626	99 write(,) 'Error: Unable to create/update radar LUT file: ', &
627	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
628	return
629
630	end subroutine save_scale_LUTs
631	! ##############################################################################################
632	! ##############################################################################################
633	subroutine quickbeam_init()
634
635
636	end subroutine quickBeam_init
637	! ##############################################################################################
638	! ##############################################################################################
639
640
641	end module quickbeam
642
643

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

Download in other formats:

Original Format