Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

quickbeam.F90 @ 5018

Last change on this file since 5018 was 3491, checked in by idelkadi, 6 years ago

Integration of version 2 of the COSP simulator in LMDZ
This line, and those below, will be ignored--

M makegcm
M makelmdz
M makelmdz_fcm
M libf/phylmd/physiq_mod.F90
A libf/phylmd/cospv2
A libf/phylmd/cospv2/mo_rng.F90
A libf/phylmd/cospv2/quickbeam_optics.F90
A libf/phylmd/cospv2/cosp_cloudsat_interface.F90
A libf/phylmd/cospv2/cosp_config.F90
A libf/phylmd/cospv2/lidar_simulator.F90
A libf/phylmd/cospv2/prec_scops.F90
A libf/phylmd/cospv2/mrgrnk.F90
A libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
A libf/phylmd/cospv2/cosp_atlid_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
A libf/phylmd/cospv2/cosp_math_constants.F90
A libf/phylmd/cospv2/MISR_simulator.F90
A libf/phylmd/cospv2/modis_simulator.F90
A libf/phylmd/cospv2/math_lib.F90
A libf/phylmd/cospv2/cosp_grLidar532_interface.F90
A libf/phylmd/cospv2/cosp_errorHandling.F90
A libf/phylmd/cospv2/cosp_stats.F90
A libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
A libf/phylmd/cospv2/cosp_utils.F90
A libf/phylmd/cospv2/cosp_optics.F90
A libf/phylmd/cospv2/icarus.F90
A libf/phylmd/cospv2/scops.F90
A libf/phylmd/cospv2/optics_lib.F90
A libf/phylmd/cospv2/cosp_kinds.F90
A libf/phylmd/cospv2/cosp_calipso_interface.F90
A libf/phylmd/cospv2/quickbeam.F90
A libf/phylmd/cospv2/parasol.F90
A libf/phylmd/cospv2/cosp_phys_constants.F90
A libf/phylmd/cospv2/cosp.F90
A libf/phylmd/cospv2/array_lib.F90
A libf/phylmd/cospv2/cosp_isccp_interface.F90
A libf/phylmd/cospv2/cosp_parasol_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_interface.F90
A libf/phylmd/cospv2/cosp_misr_interface.F90
A libf/phylmd/cospv2/cosp_modis_interface.F90

File size: 31.2 KB

Rev	Line
[3491]	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