Context Navigation

quickbeam.F90 @ 5072

Last change on this file since 5072 was 3358, checked in by idelkadi, 6 years ago
Implementation de la nouvelle version COSPv2 dans LMDZ. Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 17.6 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: DBZE_BINS,DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH, &
53	R_UNDEF,cloudsat_histRef,use_vgrid,vgrid_zl,vgrid_zu
54	USE MOD_COSP_STATS, ONLY: COSP_LIDAR_ONLY_CLOUD,hist1D,COSP_CHANGE_VERTICAL_GRID
55	implicit none
56
57	integer,parameter :: &
58	maxhclass = 20, & ! Qucikbeam maximum number of hydrometeor classes.
59	nRe_types = 550, & ! Quickbeam maximum number or Re size bins allowed in N and Z_scaled look up table.
60	nd = 85, & ! Qucikbeam number of discrete particles used in construction DSDs.
61	mt_ntt = 39, & ! Quickbeam number of temperatures in mie LUT.
62	Re_BIN_LENGTH = 10, & ! Quickbeam minimum Re interval in scale LUTs
63	Re_MAX_BIN = 250 ! Quickbeam maximum Re interval in scale LUTs
64	real(wp),parameter :: &
65	dmin = 0.1, & ! Quickbeam minimum size of discrete particle
66	dmax = 10000. ! Quickbeam maximum size of discrete particle
67
68	!djs logical :: radar_at_layer_one ! If true radar is assume to be at the edge
69	! of the first layer, if the first layer is the
70	! surface than a ground-based radar. If the
71	! first layer is the top-of-atmosphere, then
72	! a space borne radar.
73
74	! ##############################################################################################
75	type radar_cfg
76	! Radar properties
77	real(wp) :: freq,k2
78	integer :: nhclass ! Number of hydrometeor classes in use
79	integer :: use_gas_abs, do_ray
80	logical :: radar_at_layer_one ! If true radar is assume to be at the edge
81	! of the first layer, if the first layer is the
82	! surface than a ground-based radar. If the
83	! first layer is the top-of-atmosphere, then
84	! a space borne radar.
85
86	! Variables used to store Z scale factors
87	character(len=240) :: scale_LUT_file_name
88	logical :: load_scale_LUTs, update_scale_LUTs
89	logical, dimension(maxhclass,nRe_types) :: N_scale_flag
90	logical, dimension(maxhclass,mt_ntt,nRe_types) :: Z_scale_flag,Z_scale_added_flag
91	real(wp),dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled
92	real(wp),dimension(maxhclass,nd,nRe_types) :: fc, rho_eff
93	real(wp),dimension(Re_MAX_BIN) :: base_list,step_list
94
95	end type radar_cfg
96
97	contains
98	! ######################################################################################
99	! SUBROUTINE quickbeam_subcolumn
100	! ######################################################################################
101	!subroutine quickbeam_subcolumn(rcfg,nprof,ngate,hgt_matrix,z_vol,kr_vol,g_vol,&
102	! a_to_vol,g_to_vol,dBZe,Ze_non,Ze_ray)
103	subroutine quickbeam_subcolumn(rcfg,nprof,ngate,hgt_matrix,z_vol,kr_vol,g_vol,dBZe)
104
105	! INPUTS
106	type(radar_cfg),intent(inout) :: &
107	rcfg ! Derived type for radar simulator setup
108	integer,intent(in) :: &
109	nprof, & ! Number of hydrometeor profiles
110	ngate ! Number of vertical layers
111	real(wp),intent(in),dimension(nprof,ngate) :: &
112	hgt_matrix, & ! Height of hydrometeors (km)
113	z_vol, & ! Effective reflectivity factor (mm^6/m^3)
114	kr_vol, & ! Attenuation coefficient hydro (dB/km)
115	g_vol ! Attenuation coefficient gases (dB/km)
116
117	! OUTPUTS
118	real(wp), intent(out),dimension(nprof,ngate) :: &
119	! Ze_non, & ! Radar reflectivity without attenuation (dBZ)
120	! Ze_ray, & ! Rayleigh reflectivity (dBZ)
121	! g_to_vol, & ! Gaseous atteunation, radar to vol (dB)
122	! a_to_vol, & ! Hydromets attenuation, radar to vol (dB)
123	dBZe ! Effective radar reflectivity factor (dBZ)
124
125	! LOCAL VARIABLES
126	integer :: k,pr,start_gate,end_gate,d_gate
127	real(wp),dimension(nprof,ngate) :: &
128	Ze_non, & ! Radar reflectivity without attenuation (dBZ)
129	Ze_ray, & ! Rayleigh reflectivity (dBZ)
130	g_to_vol, & ! Gaseous atteunation, radar to vol (dB)
131	a_to_vol, & ! Hydromets attenuation, radar to vol (dB)
132	z_ray ! Reflectivity factor, Rayleigh only (mm^6/m^3)
133
134	! Load scaling matricies from disk -- but only the first time this subroutine is called
135	if(rcfg%load_scale_LUTs) then
136	call load_scale_LUTs(rcfg)
137	rcfg%load_scale_LUTs=.false.
138	rcfg%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
139	endif
140
141	! Initialization
142	g_to_vol = 0._wp
143	a_to_vol = 0._wp
144
145	! Loop over each range gate (ngate) ... starting with layer closest to the radar !
146	if(rcfg%radar_at_layer_one) then
147	start_gate = 1
148	end_gate = ngate
149	d_gate = 1
150	else
151	start_gate = ngate
152	end_gate = 1
153	d_gate = -1
154	endif
155	do k=start_gate,end_gate,d_gate
156	! Loop over each profile (nprof)
157	do pr=1,nprof
158	! Attenuation due to hydrometeors between radar and volume
159
160	! NOTE old scheme integrates attenuation only for the layers ABOVE
161	! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
162	! which may be a problem. ROJ
163	! in the new scheme I assign half the attenuation to the current layer
164	if(d_gate==1) then
165	! dheight calcuations assumes hgt_matrix points are the cell mid-points.
166	if (k>2) then
167	! add to previous value to half of above layer + half of current layer
168	a_to_vol(pr,k)= a_to_vol(pr,k-1) + &
169	(kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
170	else
171	a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
172	endif
173	else ! d_gate==-1
174	if(k<ngate) then
175	! Add to previous value half of above layer + half of current layer
176	a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
177	(kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
178	else
179	a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
180	endif
181	endif
182
183	! Attenuation due to gaseous absorption between radar and volume
184	if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr .eq. 1)) then
185	if (d_gate==1) then
186	if (k>1) then
187	! Add to previous value to half of above layer + half of current layer
188	g_to_vol(pr,k) = g_to_vol(pr,k-1) + &
189	0.5(g_vol(pr,k-1)+g_vol(pr,k))(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
190	else
191	g_to_vol(pr,k)= 0.5_wpg_vol(pr,k)(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
192	endif
193	else ! d_gate==-1
194	if (k<ngate) then
195	! Add to previous value to half of above layer + half of current layer
196	g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
197	0.5_wp(g_vol(pr,k+1)+g_vol(pr,k))(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
198	else
199	g_to_vol(pr,k)= 0.5_wpg_vol(pr,k)(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
200	endif
201	endif
202	elseif(rcfg%use_gas_abs == 2) then
203	! Using value calculated for the first column
204	g_to_vol(pr,k) = g_to_vol(1,k)
205	elseif (rcfg%use_gas_abs == 0) then
206	g_to_vol(pr,k) = 0._wp
207	endif
208	enddo ! End loop over pr (profile)
209	enddo ! End loop of k (range gate)
210
211	! Compute Rayleigh reflectivity, and full, attenuated reflectivity
212	if(rcfg%do_ray == 1) then
213	where(z_ray(1:nprof,1:ngate) > 0._wp)
214	Ze_ray(1:nprof,1:ngate) = 10._wp*log10(z_ray(1:nprof,1:ngate))
215	elsewhere
216	Ze_Ray(1:nprof,1:ngate) = 0._wp
217	endwhere
218	!djs Ze_ray(1:nprof,1:ngate) = merge(10._wplog10(z_ray(1:nprof,1:ngate)), 1._wpR_UNDEF, z_ray(1:nprof,1:ngate) > 0._wp)
219	else
220	Ze_ray(1:nprof,1:ngate) = R_UNDEF
221	end if
222
223	where(z_vol(1:nprof,1:ngate) > 0._wp)
224	Ze_non(1:nprof,1:ngate) = 10._wp*log10(z_vol(1:nprof,1:ngate))
225	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)
226	elsewhere
227	dBZe(1:nprof,1:ngate) = R_UNDEF
228	Ze_non(1:nprof,1:ngate) = R_UNDEF
229	end where
230
231	! Save any updates made
232	if (rcfg%update_scale_LUTs) call save_scale_LUTs(rcfg)
233
234	end subroutine quickbeam_subcolumn
235	! ######################################################################################
236	! SUBROUTINE quickbeam_column
237	! ######################################################################################
238	subroutine quickbeam_column(npoints,ncolumns,nlevels,llm,Ze_tot,zlev,zlev_half,cfad_ze)
239	! Inputs
240	integer,intent(in) :: &
241	npoints, & ! Number of horizontal grid points
242	ncolumns, & ! Number of subcolumns
243	nlevels, & ! Number of vertical layers in OLD grid
244	llm ! NUmber of vertical layers in NEW grid
245	real(wp),intent(in),dimension(npoints,ncolumns,Nlevels) :: &
246	Ze_tot !
247	real(wp),intent(in),dimension(npoints,Nlevels) :: &
248	zlev ! Model full levels
249	real(wp),intent(in),dimension(npoints,Nlevels+1) :: &
250	zlev_half ! Model half levels
251
252	! Outputs
253	real(wp),intent(inout),dimension(npoints,DBZE_BINS,llm) :: &
254	cfad_ze !
255
256	! Local variables
257	integer :: i,j
258	real(wp),dimension(npoints,ncolumns,llm) :: ze_totFlip
259
260	if (use_vgrid) then
261	! Regrid in the vertical
262	call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),&
263	zlev_half(:,nlevels:1:-1),Ze_tot(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),&
264	vgrid_zu(llm:1:-1),Ze_totFlip(:,:,llm:1:-1),log_units=.true.)
265
266	! Effective reflectivity histogram
267	do i=1,Npoints
268	do j=1,llm
269	cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_totFlip(i,:,j),DBZE_BINS,cloudsat_histRef)
270	enddo
271	enddo
272	where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
273
274	else
275	! Effective reflectivity histogram
276	do i=1,Npoints
277	do j=1,llm
278	cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_tot(i,:,j),DBZE_BINS,cloudsat_histRef)
279	enddo
280	enddo
281	where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
282	endif
283
284	end subroutine quickbeam_column
285	! ##############################################################################################
286	! ##############################################################################################
287
288
289	! ##############################################################################################
290	! ##############################################################################################
291	subroutine load_scale_LUTs(rcfg)
292
293	type(radar_cfg), intent(inout) :: rcfg
294	logical :: LUT_file_exists
295	integer :: i,j,k,ind
296
297	! Load scale LUT from file
298	inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
299	exist=LUT_file_exists)
300
301	if(.not.LUT_file_exists) then
302	write(,) '*************************************************'
303	write(,) 'Warning: Could NOT FIND radar LUT file: ', &
304	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
305	write(,) 'Will calculated LUT values as needed'
306	write(,) '*************************************************'
307	return
308	else
309	OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
310	form='unformatted', &
311	err= 89, &
312	access='DIRECT',&
313	recl=28)
314	write(,) 'Loading radar LUT file: ', &
315	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
316
317	do i=1,maxhclass
318	do j=1,mt_ntt
319	do k=1,nRe_types
320	ind = i+(j-1)maxhclass+(k-1)(nRe_types*mt_ntt)
321	read(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
322	rcfg%Ze_scaled(i,j,k), &
323	rcfg%Zr_scaled(i,j,k), &
324	rcfg%kr_scaled(i,j,k)
325	enddo
326	enddo
327	enddo
328	close(unit=12)
329	return
330	endif
331
332	89 write(,) 'Error: Found but could NOT READ radar LUT file: ', &
333	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
334
335	end subroutine load_scale_LUTs
336
337	! ##############################################################################################
338	! ##############################################################################################
339	subroutine save_scale_LUTs(rcfg)
340	type(radar_cfg), intent(inout) :: rcfg
341	logical :: LUT_file_exists
342	integer :: i,j,k,ind
343
344	inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
345	exist=LUT_file_exists)
346
347	OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
348	form='unformatted',err= 99,access='DIRECT',recl=28)
349
350	write(,) 'Creating or Updating radar LUT file: ', &
351	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
352
353	do i=1,maxhclass
354	do j=1,mt_ntt
355	do k=1,nRe_types
356	ind = i+(j-1)maxhclass+(k-1)(nRe_types*mt_ntt)
357	if(.not.LUT_file_exists .or. rcfg%Z_scale_added_flag(i,j,k)) then
358	rcfg%Z_scale_added_flag(i,j,k)=.false.
359	write(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
360	rcfg%Ze_scaled(i,j,k), &
361	rcfg%Zr_scaled(i,j,k), &
362	rcfg%kr_scaled(i,j,k)
363	endif
364	enddo
365	enddo
366	enddo
367	close(unit=12)
368	return
369
370	99 write(,) 'Error: Unable to create/update radar LUT file: ', &
371	trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
372	return
373
374	end subroutine save_scale_LUTs
375	! ##############################################################################################
376	! ##############################################################################################
377	subroutine quickbeam_init()
378
379
380	end subroutine quickBeam_init
381	! ##############################################################################################
382	! ##############################################################################################
383
384
385	end module quickbeam
386
387

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

Download in other formats:

Original Format