source: LMDZ6/trunk/libf/phylmd/cosp2/quickbeam.F90

Last change on this file 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! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50module 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
97contains
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_wp*g_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_wp*g_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._wp*log10(z_ray(1:nprof,1:ngate)), 1._wp*R_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   
33289  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   
37099  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
385end module quickbeam
386
387
Note: See TracBrowser for help on using the repository browser.