! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Copyright (c) 2015, Regents of the University of Colorado ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without modification, are ! permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of ! conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list ! of conditions and the following disclaimer in the documentation and/or other ! materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be ! used to endorse or promote products derived from this software without specific prior ! written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! History ! 11/2005: John Haynes - Created ! 09/2006 placed into subroutine form (Roger Marchand,JMH) ! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand) ! 01/2008 'Do while' to determine if hydrometeor(s) present in volume ! changed for vectorization purposes (A. Bodas-Salcedo) ! ! 07/2010 V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT) ! ... All hydrometeor and radar simulator properties now included in hp structure ! ... hp structure should be initialized by call to radar_simulator_init prior ! ... to calling this subroutine. ! Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added ! ... Changes implement by Roj Marchand following work by Laura Fowler ! ! 10/2011 Modified ngate loop to go in either direction depending on flag ! hp%radar_at_layer_one. This affects the direction in which attenuation is summed. ! ! Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple ! summation. (Roger Marchand) ! May 2015 - D. Swales - Modified for COSPv2.0 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% module quickbeam USE COSP_KINDS, ONLY: wp USE MOD_COSP_CONFIG, ONLY: R_UNDEF,cloudsat_histRef,use_vgrid,vgrid_zl,vgrid_zu,& pClass_noPrecip, pClass_Rain1, pClass_Rain2, pClass_Rain3,& pClass_Snow1, pClass_Snow2, pClass_Mixed1, pClass_Mixed2, & pClass_Rain4, pClass_default, Zenonbinval, Zbinvallnd, & N_HYDRO,nCloudsatPrecipClass !,cloudsat_preclvl !PREC_BUG USE MOD_COSP_STATS, ONLY: COSP_LIDAR_ONLY_CLOUD,hist1D,COSP_CHANGE_VERTICAL_GRID implicit none integer,parameter :: & maxhclass = 20, & ! Qucikbeam maximum number of hydrometeor classes. nRe_types = 550, & ! Quickbeam maximum number or Re size bins allowed in N and Z_scaled look up table. nd = 85, & ! Qucikbeam number of discrete particles used in construction DSDs. mt_ntt = 39, & ! Quickbeam number of temperatures in mie LUT. Re_BIN_LENGTH = 10, & ! Quickbeam minimum Re interval in scale LUTs Re_MAX_BIN = 250 ! Quickbeam maximum Re interval in scale LUTs real(wp),parameter :: & dmin = 0.1, & ! Quickbeam minimum size of discrete particle dmax = 10000. ! Quickbeam maximum size of discrete particle !djs logical :: radar_at_layer_one ! If true radar is assume to be at the edge ! of the first layer, if the first layer is the ! surface than a ground-based radar. If the ! first layer is the top-of-atmosphere, then ! a space borne radar. ! ############################################################################################## type radar_cfg ! Radar properties real(wp) :: freq,k2 integer :: nhclass ! Number of hydrometeor classes in use integer :: use_gas_abs, do_ray logical :: radar_at_layer_one ! If true radar is assume to be at the edge ! of the first layer, if the first layer is the ! surface than a ground-based radar. If the ! first layer is the top-of-atmosphere, then ! a space borne radar. ! Variables used to store Z scale factors character(len=240) :: scale_LUT_file_name logical :: load_scale_LUTs, update_scale_LUTs logical, dimension(maxhclass,nRe_types) :: N_scale_flag logical, dimension(maxhclass,mt_ntt,nRe_types) :: Z_scale_flag,Z_scale_added_flag real(wp),dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled real(wp),dimension(maxhclass,nd,nRe_types) :: fc, rho_eff real(wp),dimension(Re_MAX_BIN) :: base_list,step_list end type radar_cfg contains ! ###################################################################################### ! SUBROUTINE quickbeam_subcolumn ! ###################################################################################### subroutine quickbeam_subcolumn(rcfg,nprof,ngate,hgt_matrix,z_vol,kr_vol,g_vol,dBZe,Ze_non) ! INPUTS type(radar_cfg),intent(inout) :: & rcfg ! Derived type for radar simulator setup integer,intent(in) :: & nprof, & ! Number of hydrometeor profiles ngate ! Number of vertical layers real(wp),intent(in),dimension(nprof,ngate) :: & hgt_matrix, & ! Height of hydrometeors (km) z_vol, & ! Effective reflectivity factor (mm^6/m^3) kr_vol, & ! Attenuation coefficient hydro (dB/km) g_vol ! Attenuation coefficient gases (dB/km) ! OUTPUTS real(wp), intent(out),dimension(nprof,ngate) :: & Ze_non, & ! Radar reflectivity without attenuation (dBZ) dBZe ! Effective radar reflectivity factor (dBZ) ! LOCAL VARIABLES integer :: k,pr,start_gate,end_gate,d_gate real(wp),dimension(nprof,ngate) :: & Ze_ray, & ! Rayleigh reflectivity (dBZ) g_to_vol, & ! Gaseous atteunation, radar to vol (dB) a_to_vol, & ! Hydromets attenuation, radar to vol (dB) z_ray ! Reflectivity factor, Rayleigh only (mm^6/m^3) ! Load scaling matricies from disk -- but only the first time this subroutine is called if(rcfg%load_scale_LUTs) then call load_scale_LUTs(rcfg) rcfg%load_scale_LUTs=.false. rcfg%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run endif ! Initialization g_to_vol = 0._wp a_to_vol = 0._wp ! Loop over each range gate (ngate) ... starting with layer closest to the radar ! if(rcfg%radar_at_layer_one) then start_gate = 1 end_gate = ngate d_gate = 1 else start_gate = ngate end_gate = 1 d_gate = -1 endif do k=start_gate,end_gate,d_gate ! Loop over each profile (nprof) do pr=1,nprof ! Attenuation due to hydrometeors between radar and volume ! NOTE old scheme integrates attenuation only for the layers ABOVE ! the current layer ... i.e. 1 to k-1 rather than 1 to k ... ! which may be a problem. ROJ ! in the new scheme I assign half the attenuation to the current layer if(d_gate==1) then ! dheight calcuations assumes hgt_matrix points are the cell mid-points. if (k>2) then ! add to previous value to half of above layer + half of current layer a_to_vol(pr,k)= a_to_vol(pr,k-1) + & (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) else a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1)) endif else ! d_gate==-1 if(k<ngate) then ! Add to previous value half of above layer + half of current layer a_to_vol(pr,k) = a_to_vol(pr,k+1) + & (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) else a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)) endif endif ! Attenuation due to gaseous absorption between radar and volume if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr .eq. 1)) then if (d_gate==1) then if (k>1) then ! Add to previous value to half of above layer + half of current layer g_to_vol(pr,k) = g_to_vol(pr,k-1) + & 0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) else g_to_vol(pr,k)= 0.5_wp*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1)) endif else ! d_gate==-1 if (k<ngate) then ! Add to previous value to half of above layer + half of current layer g_to_vol(pr,k) = g_to_vol(pr,k+1) + & 0.5_wp*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) else g_to_vol(pr,k)= 0.5_wp*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)) endif endif elseif(rcfg%use_gas_abs == 2) then ! Using value calculated for the first column g_to_vol(pr,k) = g_to_vol(1,k) elseif (rcfg%use_gas_abs == 0) then g_to_vol(pr,k) = 0._wp endif enddo ! End loop over pr (profile) enddo ! End loop of k (range gate) ! Compute Rayleigh reflectivity, and full, attenuated reflectivity if(rcfg%do_ray == 1) then where(z_ray(1:nprof,1:ngate) > 0._wp) Ze_ray(1:nprof,1:ngate) = 10._wp*log10(z_ray(1:nprof,1:ngate)) elsewhere Ze_Ray(1:nprof,1:ngate) = 0._wp endwhere else Ze_ray(1:nprof,1:ngate) = R_UNDEF end if where(z_vol(1:nprof,1:ngate) > 0._wp) Ze_non(1:nprof,1:ngate) = 10._wp*log10(z_vol(1:nprof,1:ngate)) 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) elsewhere dBZe(1:nprof,1:ngate) = R_UNDEF Ze_non(1:nprof,1:ngate) = R_UNDEF end where ! Save any updates made if (rcfg%update_scale_LUTs) call save_scale_LUTs(rcfg) end subroutine quickbeam_subcolumn ! ###################################################################################### ! SUBROUTINE quickbeam_column ! ###################################################################################### subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform, & Ze_tot, Ze_tot_non, land, surfelev, t2m, fracPrecipIce, zlev, zlev_half, cfad_ze, & !PREC_BUG cloudsat_precip_cover, cloudsat_pia) ! Inputs integer,intent(in) :: & npoints, & ! Number of horizontal grid points ncolumns, & ! Number of subcolumns nlevels, & ! Number of vertical layers in OLD grid llm, & ! Number of vertical layers in NEW grid DBZE_BINS ! Number of bins for cfad. character(len=*),intent(in) :: & platform ! Name of platform (e.g. cloudsat) real(wp),dimension(Npoints),intent(in) :: & land, & ! Land/Sea mask. (1/0) surfelev, & ! Surface Elevation (m) !PREC_BUG t2m ! Near-surface temperature real(wp),dimension(Npoints,Ncolumns),intent(in) :: & fracPrecipIce ! Fraction of precipitation which is frozen. (1) real(wp),intent(in),dimension(npoints,ncolumns,Nlevels) :: & Ze_tot, & ! Effective reflectivity factor (dBZ) Ze_tot_non ! Effective reflectivity factor w/o attenuation (dBZ) real(wp),intent(in),dimension(npoints,Nlevels) :: & zlev ! Model full levels real(wp),intent(in),dimension(npoints,Nlevels+1) :: & zlev_half ! Model half levels ! Outputs real(wp),intent(inout),dimension(npoints,DBZE_BINS,llm) :: & cfad_ze ! real(wp),dimension(Npoints,nCloudsatPrecipClass),intent(out) :: & cloudsat_precip_cover ! Model precip rate in by CloudSat precip flag real(wp),dimension(Npoints),intent(out) :: & cloudsat_pia ! Cloudsat path integrated attenuation ! Local variables integer :: i,j real(wp),dimension(npoints,ncolumns,llm) :: ze_toti,ze_noni logical :: lcloudsat = .false. ! Which platforms to create diagnostics for? if (platform .eq. 'cloudsat') lcloudsat=.true. ! Create Cloudsat diagnostics. if (lcloudsat) then if (use_vgrid) then ! Regrid in the vertical (*NOTE* This routine requires SFC-2-TOA ordering, so flip ! inputs and outputs to maintain TOA-2-SFC ordering convention in COSP2.) call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),& zlev_half(:,nlevels:1:-1),Ze_tot(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),& vgrid_zu(llm:1:-1),Ze_toti(:,:,llm:1:-1),log_units=.true.) ! Effective reflectivity histogram do i=1,Npoints do j=1,llm cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_toti(i,:,j),DBZE_BINS,cloudsat_histRef) enddo enddo where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns ! Compute cloudsat near-surface precipitation diagnostics ! First, regrid in the vertical Ze_tot_non. call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),& zlev_half(:,nlevels:1:-1),Ze_tot_non(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),& vgrid_zu(llm:1:-1),Ze_noni(:,:,llm:1:-1),log_units=.true.) ! Not call routine to generate diagnostics. call cloudsat_precipOccurence(Npoints, Ncolumns, llm, N_HYDRO, Ze_toti, Ze_noni, & land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia) !PREC_BUG else ! Effective reflectivity histogram do i=1,Npoints do j=1,llm cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_tot(i,:,j),DBZE_BINS,cloudsat_histRef) enddo enddo where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns endif endif end subroutine quickbeam_column ! ############################################################################################## ! ############################################################################################## ! ###################################################################################### ! SUBROUTINE cloudsat_precipOccurence ! ! Notes from July 2016: Add precip flag also looped over subcolumns ! Modified by Tristan L'Ecuyer (TSL) to add precipitation flagging ! based on CloudSat's 2C-PRECIP-COLUMN algorithm (Haynes et al, JGR, 2009). ! To mimic the satellite algorithm, this code applies thresholds to non-attenuated ! reflectivities, Ze_non, consistent with those outlined in Haynes et al, JGR (2009). ! ! Procedures/Notes: ! ! (1) If the 2-way attenuation exceeds 40 dB, the pixel will be flagged as 'heavy rain' ! consistent with the multiple-scattering analysis of Battaglia et al, JGR (2008). ! (2) Rain, snow, and mixed precipitation scenes are separated according to the fraction ! of the total precipitation hydrometeor mixing ratio that exists as ice. ! (3) The entire analysis is applied to the range gate from 480-960 m to be consistent with ! CloudSat's 720 m ground-clutter. ! (4) Only a single flag is assigned to each model grid since there is no variation in ! hydrometeor contents across a single model level. Unlike CFADs, whose variation enters ! due to differening attenuation corrections from hydrometeors aloft, the non-attenuated ! reflectivities used in the computation of this flag cannot vary across sub-columns. ! ! radar_prec_flag = 1-Rain possible 2-Rain probable 3-Rain certain ! 4-Snow possible 5-Snow certain ! 6-Mixed possible 7-Mixed certain ! 8-Heavy Rain ! 9- default value ! ! Modified by Dustin Swales (University of Colorado) for use with COSP2. ! *NOTE* All inputs (Ze_out, Ze_non_out, fracPrecipIce) are at a single level from the ! statistical output grid used by Cloudsat. This level is controlled by the ! parameter cloudsat_preclvl, defined in src/cosp_config.F90 ! ###################################################################################### subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_non_out, & land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia) !PREC_BUG ! Inputs integer,intent(in) :: & Npoints, & ! Number of columns Ncolumns, & ! Numner of subcolumns Nhydro, & ! Number of hydrometeor types llm ! Number of levels real(wp),dimension(Npoints),intent(in) :: & land, & ! Land/Sea mask. (1/0) surfelev, & ! Surface Elevation (m) !PREC_BUG t2m ! Near-surface temperature real(wp),dimension(Npoints,Ncolumns,llm),intent(in) :: & Ze_out, & ! Effective reflectivity factor (dBZ) Ze_non_out ! Effective reflectivity factor, w/o attenuation (dBZ) real(wp),dimension(Npoints,Ncolumns),intent(in) :: & fracPrecipIce ! Fraction of precipitation which is frozen. (1) ! Outputs real(wp),dimension(Npoints,nCloudsatPrecipClass),intent(out) :: & cloudsat_precip_cover ! Model precip rate in by CloudSat precip flag real(wp),dimension(Npoints),intent(out) :: & cloudsat_pia ! Cloudsat path integrated attenuation ! Local variables integer,dimension(Npoints,Ncolumns) :: & cloudsat_pflag, & ! Subcolumn precipitation flag cloudsat_precip_pia ! Subcolumn path integrated attenutation. integer,dimension(Npoints) :: & !PREC_BUG cloudsat_preclvl_index ! Altitude index for precip flags calculation !PREC_BUG ! in 40-level grid (2nd layer above surfelev) !PREC_BUG integer :: pr,i,k,m,j,cloudsat_preclvl !PREC_BUG real(wp) :: Zmax ! Initialize cloudsat_pflag(:,:) = pClass_default cloudsat_precip_pia(:,:) = 0._wp cloudsat_precip_cover(:,:) = 0._wp cloudsat_pia(:) = 0._wp cloudsat_preclvl_index(:) = 0._wp !PREC_BUG !!! Computing altitude index for precip flags calculation !PREC_BUG cloudsat_preclvl_index(:) = 39 - floor( surfelev(:)/480. ) !PREC_BUG ! ###################################################################################### ! SUBCOLUMN processing ! ###################################################################################### do i=1, Npoints cloudsat_preclvl = cloudsat_preclvl_index(i) !PREC_BUG ! print*, 'i, surfelev(i), cloudsat_preclvl : ', i, surfelev(i), cloudsat_preclvl !PREC_BUG do pr=1,Ncolumns ! 1) Compute the PIA in all profiles containing hydrometeors if ( (Ze_non_out(i,pr,cloudsat_preclvl).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then if ( (Ze_non_out(i,pr,cloudsat_preclvl).lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl) endif endif ! 2) Compute precipitation flag ! ################################################################################ ! 2a) Oceanic points. ! ################################################################################ if (land(i) .eq. 0) then ! print*, 'aaa i, pr, fracPrecipIce(i,pr) : ', i, pr, fracPrecipIce(i,pr) !Artem ! Snow if(fracPrecipIce(i,pr).gt.0.9) then if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then cloudsat_pflag(i,pr) = pClass_Snow2 ! TSL: Snow certain endif if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. & Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then cloudsat_pflag(i,pr) = pClass_Snow1 ! TSL: Snow possible endif endif ! Mixed if(fracPrecipIce(i,pr).gt.0.1.and.fracPrecipIce(i,pr).le.0.9) then if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then cloudsat_pflag(i,pr) = pClass_Mixed2 ! TSL: Mixed certain endif if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. & Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then cloudsat_pflag(i,pr) = pClass_Mixed1 ! TSL: Mixed possible endif endif ! Rain if(fracPrecipIce(i,pr).le.0.1) then if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(1)) then cloudsat_pflag(i,pr) = pClass_Rain3 ! TSL: Rain certain endif if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(3).and. & Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(1)) then cloudsat_pflag(i,pr) = pClass_Rain2 ! TSL: Rain probable endif if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. & Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(3)) then cloudsat_pflag(i,pr) = pClass_Rain1 ! TSL: Rain possible endif if(cloudsat_precip_pia(i,pr).gt.40) then cloudsat_pflag(i,pr) = pClass_Rain4 ! TSL: Heavy Rain endif endif ! No precipitation if(Ze_non_out(i,pr,cloudsat_preclvl).le.-15) then cloudsat_pflag(i,pr) = pClass_noPrecip ! TSL: Not Raining endif endif ! Ocean points ! ################################################################################ ! 2b) Land points. ! ################################################################################ if (land(i) .eq. 1) then ! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out); Zmax=maxval(Ze_out(i,pr,:)) ! Snow (T<273) if(t2m(i) .lt. 273._wp) then if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(5)) then cloudsat_pflag(i,pr) = pClass_Snow2 ! JEK: Snow certain endif if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .and. & Ze_out(i,pr,cloudsat_preclvl).le.Zbinvallnd(5)) then cloudsat_pflag(i,pr) = pClass_Snow1 ! JEK: Snow possible endif endif ! Mized phase (273<T<275) if(t2m(i) .ge. 273._wp .and. t2m(i) .le. 275._wp) then if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. & (Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(4))) then cloudsat_pflag(i,pr) = pClass_Mixed2 ! JEK: Mixed certain endif if ((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .and. & Ze_out(i,pr,cloudsat_preclvl) .le. Zbinvallnd(4)) .and. & (Zmax .gt. Zbinvallnd(5)) ) then cloudsat_pflag(i,pr) = pClass_Mixed1 ! JEK: Mixed possible endif endif ! Rain (T>275) if(t2m(i) .gt. 275) then if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. & (Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(2))) then cloudsat_pflag(i,pr) = pClass_Rain3 ! JEK: Rain certain endif if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. & (Zmax .gt. Zbinvallnd(3))) then cloudsat_pflag(i,pr) = pClass_Rain2 ! JEK: Rain probable endif if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. & (Zmax.lt.Zbinvallnd(3))) then cloudsat_pflag(i,pr) = pClass_Rain1 ! JEK: Rain possible endif if(cloudsat_precip_pia(i,pr).gt.40) then cloudsat_pflag(i,pr) = pClass_Rain4 ! JEK: Heavy Rain endif endif ! No precipitation if(Ze_out(i,pr,cloudsat_preclvl).le.-15) then cloudsat_pflag(i,pr) = pClass_noPrecip ! JEK: Not Precipitating endif endif ! Land points enddo ! Sub-columns enddo ! Gridpoints ! ###################################################################################### ! COLUMN processing ! ###################################################################################### ! Aggregate subcolumns do i=1,Npoints ! Gridmean precipitation fraction for each precipitation type do k=1,nCloudsatPrecipClass if (any(cloudsat_pflag(i,:) .eq. k-1)) then cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq. k-1) endif enddo ! Gridmean path integrated attenuation (pia) cloudsat_pia(i)=sum(cloudsat_precip_pia(i,:)) enddo ! Normalize by number of subcolumns where ((cloudsat_precip_cover /= R_UNDEF).and.(cloudsat_precip_cover /= 0.0)) & cloudsat_precip_cover = cloudsat_precip_cover / Ncolumns where ((cloudsat_pia/= R_UNDEF).and.(cloudsat_pia/= 0.0)) & cloudsat_pia = cloudsat_pia / Ncolumns end subroutine cloudsat_precipOccurence ! ############################################################################################## ! ############################################################################################## subroutine load_scale_LUTs(rcfg) type(radar_cfg), intent(inout) :: rcfg logical :: LUT_file_exists integer :: i,j,k,ind ! Load scale LUT from file inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', & exist=LUT_file_exists) if(.not.LUT_file_exists) then write(*,*) '*************************************************' write(*,*) 'Warning: Could NOT FIND radar LUT file: ', & trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' write(*,*) 'Will calculated LUT values as needed' write(*,*) '*************************************************' return else OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',& form='unformatted', & err= 89, & access='DIRECT',& recl=28) write(*,*) 'Loading radar LUT file: ', & trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' do i=1,maxhclass do j=1,mt_ntt do k=1,nRe_types ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) read(12,rec=ind) rcfg%Z_scale_flag(i,j,k), & rcfg%Ze_scaled(i,j,k), & rcfg%Zr_scaled(i,j,k), & rcfg%kr_scaled(i,j,k) enddo enddo enddo close(unit=12) return endif 89 write(*,*) 'Error: Found but could NOT READ radar LUT file: ', & trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' end subroutine load_scale_LUTs ! ############################################################################################## ! ############################################################################################## subroutine save_scale_LUTs(rcfg) type(radar_cfg), intent(inout) :: rcfg logical :: LUT_file_exists integer :: i,j,k,ind inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', & exist=LUT_file_exists) OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',& form='unformatted',err= 99,access='DIRECT',recl=28) write(*,*) 'Creating or Updating radar LUT file: ', & trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' do i=1,maxhclass do j=1,mt_ntt do k=1,nRe_types ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) if(.not.LUT_file_exists .or. rcfg%Z_scale_added_flag(i,j,k)) then rcfg%Z_scale_added_flag(i,j,k)=.false. write(12,rec=ind) rcfg%Z_scale_flag(i,j,k), & rcfg%Ze_scaled(i,j,k), & rcfg%Zr_scaled(i,j,k), & rcfg%kr_scaled(i,j,k) endif enddo enddo enddo close(unit=12) return 99 write(*,*) 'Error: Unable to create/update radar LUT file: ', & trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' return end subroutine save_scale_LUTs ! ############################################################################################## ! ############################################################################################## subroutine quickbeam_init() end subroutine quickBeam_init ! ############################################################################################## ! ############################################################################################## end module quickbeam