Newer
Older
!MNH_LIC Copyright 2004-2018 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! ######spl
MODULE MODI_RADAR_SCATTERING
! #############################
!
INTERFACE
SUBROUTINE RADAR_SCATTERING(PT_RAY,PRHODREF_RAY,PR_RAY,PI_RAY,PCIT_RAY,PS_RAY,PG_RAY,PVDOP_RAY, &
PELEV,PX_H,PX_V,PW_H,PW_V,PZE,PBU_MASK_RAY,PCR_RAY)
REAL, DIMENSION(:,:,:,:,:,:),INTENT(IN) :: PT_RAY ! temperature interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),INTENT(IN) :: PRHODREF_RAY !
REAL, DIMENSION(:,:,:,:,:,:),INTENT(IN) :: PR_RAY ! rainwater mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),INTENT(IN) :: PI_RAY ! pristine ice mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PCIT_RAY ! pristine ice concentration interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PS_RAY !aggregates mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PG_RAY ! graupel mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PVDOP_RAY !Doppler radial velocity interpolated along the rays
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PELEV ! elevation
REAL, DIMENSION(:), INTENT(IN) :: PX_H ! Gaussian horizontal nodes
REAL, DIMENSION(:), INTENT(IN) :: PX_V ! Gaussian vertical nodes
REAL, DIMENSION(:), INTENT(IN) :: PW_H ! Gaussian horizontal weights
REAL, DIMENSION(:), INTENT(IN) :: PW_V ! Gaussian vertical weights
REAL,DIMENSION(:,:,:,:,:), INTENT(INOUT) :: PZE ! 5D matrix (iradar, ielev, iaz, irangestep, ivar) containing the radar variables that will be calculated
!in polar or cartesian projection (same projection as the observation grid)
! convective/stratiform
REAL, DIMENSION(:,:,:,:,:,:),INTENT(INOUT) :: PBU_MASK_RAY
REAL, DIMENSION(:,:,:,:,:,:),OPTIONAL,INTENT(IN) :: PCR_RAY ! rainwater mixing ratio interpolated along the rays
END SUBROUTINE RADAR_SCATTERING
END INTERFACE
END MODULE MODI_RADAR_SCATTERING
!
! ######spl
SUBROUTINE RADAR_SCATTERING(PT_RAY,PRHODREF_RAY,PR_RAY,PI_RAY,PCIT_RAY, &
PS_RAY,PG_RAY,PVDOP_RAY,PELEV,PX_H,PX_V,PW_H,PW_V,PZE,PBU_MASK_RAY,PCR_RAY)
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
! ##############################
!
!!**** *RADAR_SCATTERING* - computes radar reflectivities.
!!
!! PURPOSE
!! -------
!! Compute equivalent reflectivities of a mixed phase cloud.
!!
!!** METHOD
!! ------
!! The reflectivities are computed using the n(D) * sigma(D) formula. The
!! equivalent reflectiviy is the sum of the reflectivity produced by the
!! the raindrops and the equivalent reflectivities of the ice crystals.
!! The latter are computed using the mass-equivalent diameter.
!! Four types of diffusion are possible : Rayleigh, Mie, T-matrix, and
!! Rayleigh-Gans (Kerker, 1969, Chap. 10; Battan, 1973, Sec. 5.4; van de
!! Hulst, 1981, Sec. 6.32; Doviak and Zrnic, 1993, p. 249; Bringi and
!! Chandrasekar, 2001, Chap. 2).
!! The integration over diameters for Mie and T-matrix methods is done by
!! using Gauss-Laguerre quadrature (Press et al. 1986). Attenuation is taken
!! into account by computing the extinction efficiency and correcting
!! reflectivities along the beam path.
!! Gaussian quadrature methods are used to model the beam broadening (Gauss-
!! Hermite or Gauss-Legendre, see Press et al. 1986).
!!
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST
!! XLIGHTSPEED
!! XPI
!! Module MODD_ARF
!!
!! REFERENCE
!! ---------
!! Press, W. H., B. P. Flannery, S. A. Teukolsky et W. T. Vetterling, 1986:
!! Numerical Recipes: The Art of Scientific Computing. Cambridge University
!! Press, 818 pp.
!! Probert-Jones, J. R., 1962 : The radar equation in meteorology. Quart.
!! J. Roy. Meteor. Soc., 88, 485-495.
!!
!! AUTHOR
!! ------
!! O. Caumont & V. Ducrocq * Meteo-France *
!!
!! MODIFICATIONS
!! -------------
!! Original 26/03/2004
!! O. Caumont 09/09/2009 minor changes to compute radial velocities when no
!! hydrometeors so as to emulate wind lidar
!! O. Caumont 21/12/2009 correction of bugs to compute KDP.
!! O. Caumont 11/02/2010 thresholding and conversion from linear to
!! log values after interpolation instead of before.
!! G.Tanguy 25/03/2010 Introduction of MODD_TMAT and ALLOCATE/DEALLOCATE
!! C.Augros 2014 New simulator for T matrice
!! G.Delautier 10/2014 : Mise a jour simulateur T-matrice pour LIMA
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
WAUTELET Philippe
committed
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_PARAM_LIMA_WARM, ONLY: XDR_L=>XDR,XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XCR_L=>XCR
USE MODD_PARAM_LIMA_COLD, ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
XDS_L=>XDS,XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS
USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
XRTMIN_L=>XRTMIN
USE MODD_RADAR, ONLY:XLAM_RAD,XSTEP_RAD,NBELEV,NDIFF,LATT,NPTS_GAULAG,LQUAD,XVALGROUND,NDGS, &
LFALL,LWBSCS,LWREFL,XREFLVDOPMIN,XREFLMIN,LSNRT,XSNRMIN
USE MODD_RAIN_ICE_DESCR, ONLY: XALPHAR_I=>XALPHAR,XNUR_I=>XNUR,XDR_I=>XDR,XLBEXR_I=>XLBEXR,&
XLBR_I=>XLBR,XCCR_I=>XCCR,XBR_I=>XBR,XCR_I=>XCR,&
XALPHAS_I=>XALPHAS,XNUS_I=>XNUS,XDS_I=>XDS,XLBEXS_I=>XLBEXS,&
XLBS_I=>XLBS,XCCS_I=>XCCS,XAS_I=>XAS,XBS_I=>XBS,XCXS_I=>XCXS,XCS_I=>XCS,&
XALPHAG_I=>XALPHAG,XNUG_I=>XNUG,XDG_I=>XDG,XLBEXG_I=>XLBEXG,&
XLBG_I=>XLBG,XCCG_I=>XCCG,XAG_I=>XAG,XBG_I=>XBG,XCXG_I=>XCXG,XCG_I=>XCG,&
XALPHAI_I=>XALPHAI,XNUI_I=>XNUI,XDI_I=>XDI,XLBEXI_I=>XLBEXI,&
XLBI_I=>XLBI,XAI_I=>XAI,XBI_I=>XBI,XC_I_I=>XC_I,&
XRTMIN_I=>XRTMIN
USE MODD_TMAT
!
USE MODE_ARF
WAUTELET Philippe
committed
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll,IO_FILE_OPEN_ll
WAUTELET Philippe
committed
USE MODE_IO_MANAGE_STRUCT, ONLY: IO_FILE_ADD2LIST
USE MODE_READTMAT
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
!
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PT_RAY ! temperature interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PRHODREF_RAY !
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PR_RAY ! rainwater mixing ratio interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PI_RAY ! pristine ice mixing ratio interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PCIT_RAY !pristine ice concentration interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PS_RAY !aggregates mixing ratio interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PG_RAY ! graupel mixing ratio interpolated along the rays
REAL,DIMENSION(:,:,:,:,:,:), INTENT(IN) :: PVDOP_RAY !Doppler radial velocity interpolated along the rays
REAL,DIMENSION(:,:,:,:), INTENT(IN) :: PELEV ! elevation
REAL,DIMENSION(:), INTENT(IN) :: PX_H ! Gaussian horizontal nodes
REAL,DIMENSION(:), INTENT(IN) :: PX_V ! Gaussian vertical nodes
REAL,DIMENSION(:), INTENT(IN) :: PW_H ! Gaussian horizontal weights
REAL,DIMENSION(:), INTENT(IN) :: PW_V ! Gaussian vertical weights
REAL,DIMENSION(:,:,:,:,:), INTENT(INOUT) :: PZE ! gate equivalent reflectivity factor (horizontal & vertical)
! convective/stratiform
REAL,DIMENSION(:,:,:,:,:,:),INTENT(INOUT) :: PBU_MASK_RAY
! /convective/stratiform
REAL, DIMENSION(:,:,:,:,:,:),OPTIONAL,INTENT(IN) :: PCR_RAY ! rainwater mixing ratio interpolated along the rays
!
!* 0.2 Declarations of local variables :
!
REAL, DIMENSION(:,:,:,:,:,:,:),ALLOCATABLE :: ZREFL
!1: ZHH (dBZ), 2: ZDR, 3: KDP, 4: CSR (0 pr air clair, 1 pour stratiforme, 2 pour convectif)
!5-8: ZER, ZEI, ZES,ZEG
!9 : VRU (vitesse radiale)
!10-13 : AER, AEI, AES, AEG
!14-17: ATR, ATI, ATS, ATG
!18-20: RhoHV, PhiDP, DeltaHV
REAL, DIMENSION(:,:,:,:,:,:,:),ALLOCATABLE :: ZAELOC ! local attenuation
REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZAETOT ! 1: total attenuation, 2: // vertical
REAL :: ZAERINT,ZAEIINT,ZAESINT,ZAEGINT ! total attenuation horizontal
REAL :: ZAVRINT,ZAVSINT,ZAVGINT ! total attenuation vertical
!
REAL,DIMENSION(:),ALLOCATABLE :: ZX,ZW ! Gauss-Laguerre points and weights
!
REAL,DIMENSION(4) :: ZREFLOC
REAL,DIMENSION(2) :: ZAETMP
REAL,DIMENSION(:),ALLOCATABLE :: ZVTEMP ! temp var for Gaussian quadrature 8 : r_r, 9 : r_i, 10 : r_s , 11 : r_g
REAL :: ZCXR=-1.0 ! for rain N ~ 1/N_0 (in Kessler parameterization)
REAL :: ZDMELT_FACT ! factor used to compute the equivalent melted diameter
REAL :: ZEQICE=0.224! factor used to convert the ice crystals reflectivity into an equivalent liquid water reflectivity (from Smith, JCAM 84)
REAL :: ZEXP ! anciliary parameter
REAL :: ZLBDA ! slope distribution parameter
REAL :: ZFRAC_ICE,ZD,ZDE ! auxiliary variables
REAL :: ZQSCA
REAL,DIMENSION(2) :: ZQEXT
REAL,DIMENSION(3) :: ZQBACK ! Q_b(HH),Q_b(VV) (backscattering efficiencies at horizontal and vertical polarizations, resp.)
!REAL :: P=DACOS(-1D0)
REAL :: ZRHOI ! pristine ice density (from m=a*D**b),
REAL :: ZRHOPI=916. !pure ice density (kg/m3)
COMPLEX :: ZNUM, ZDEN !for calculation of ice dielectri cconstant
COMPLEX :: ZQM,ZQMW,ZQMI,ZQK,ZQB, ZEPSI ! dielectric parameters
REAL :: ZS11_CARRE_R,ZS22_CARRE_R,ZRE_S22S11_R,ZIM_S22S11_R
REAL :: ZS11_CARRE_I,ZS22_CARRE_I,ZRE_S22S11_I,ZIM_S22S11_I
REAL :: ZS11_CARRE_S,ZS22_CARRE_S,ZRE_S22S11_S,ZIM_S22S11_S
REAL :: ZS11_CARRE_G,ZS22_CARRE_G,ZRE_S22S11_G,ZIM_S22S11_G
REAL :: ZS11_CARRE_T,ZS22_CARRE_T,ZRE_S22S11_T,ZIM_S22S11_T
REAL :: ZRE_S22FMS11F,ZIM_S22FT,ZIM_S11FT
REAL :: ZM
!
INTEGER :: INBRAD,IIELV,INBAZIM,INBSTEPMAX,INPTS_H,INPTS_V ! sizes of the arrays
INTEGER :: IEL
INTEGER :: JI,JL,JEL,JAZ,JH,JV,JJ,JT ! Loop variables of control
REAL :: ZLB ! depolarization factor along the spheroid symmetry axis
REAL :: ZCXI=0. ! should be defined with other parameters of microphysical scheme
REAL :: ZCR,ZCI,ZCS,ZCG ! coefficients to take into account fall speeds when simulating Doppler winds
REAL, DIMENSION(:,:,:,:),ALLOCATABLE :: ZCONC_BIN
INTEGER :: IVDOP,IMAX,IRHOHV,IPHIDP,IDELTAHV
INTEGER :: IRHR,IRHS,IRHG,IZDA,IZDS,IZDG,IKDR,IKDS,IKDG
LOGICAL :: LPART_MASK ! indicates a partial mask along the beam
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
INTEGER,PARAMETER :: IZER=5,IZEI=6,IZES=7,IZEG=8
INTEGER,PARAMETER :: IAER=10,IAEI=11,IAES=12,IAEG=13
INTEGER,PARAMETER :: IAVR=14,IAVI=15,IAVS=16,IAVG=17
INTEGER,PARAMETER :: IATR=18,IATI=19,IATS=20,IATG=21
!for ZSNR threshold
REAL ::ZDISTRAD,ZSNR,ZSNR_R,ZSNR_S,ZSNR_I,ZSNR_G,ZZHH,ZZE_R,ZZE_I,ZZE_S,ZZE_G
LOGICAL :: GTHRESHOLD_V, GTHRESHOLD_Z,GTHRESHOLD_ZR,GTHRESHOLD_ZI,GTHRESHOLD_ZS,GTHRESHOLD_ZG
!--------- TO READ T-MATRIX TABLE --------
CHARACTER(LEN=6) :: YBAND
CHARACTER(LEN=1) ::YTYPE
CHARACTER(LEN=1),DIMENSION(4) :: YTAB_TYPE
CHARACTER(LEN=25),DIMENSION(4) :: YFILE_COEFINT
REAL,DIMENSION(4) :: ZELEV_MIN,ZELEV_MAX,ZELEV_STEP,&
ZTC_MIN,ZTC_MAX,ZTC_STEP,ZFW_MIN,ZFW_MAX,ZFW_STEP
INTEGER :: IRESP,ILINE,INB_M
INTEGER,DIMENSION(4) :: INB_ELEV,INB_TC,INB_FW,INB_LINE
REAL, DIMENSION(:),ALLOCATABLE :: ZTC_T_R, ZTC_T_S, ZTC_T_G, ZTC_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZELEV_T_R, ZELEV_T_S, ZELEV_T_G, ZELEV_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZFW_T_S, ZFW_T_G, ZFW_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZM_T_R, ZM_T_S, ZM_T_G, ZM_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZS11_CARRE_T_R, ZS11_CARRE_T_S, ZS11_CARRE_T_G, ZS11_CARRE_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZS22_CARRE_T_R, ZS22_CARRE_T_S, ZS22_CARRE_T_G, ZS22_CARRE_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZRE_S22S11_T_R, ZRE_S22S11_T_S, ZRE_S22S11_T_G, ZRE_S22S11_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZIM_S22S11_T_R, ZIM_S22S11_T_S, ZIM_S22S11_T_G, ZIM_S22S11_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZIM_S22FT_T_R, ZIM_S22FT_T_S, ZIM_S22FT_T_G, ZIM_S22FT_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZIM_S11FT_T_R, ZIM_S11FT_T_S, ZIM_S11FT_T_G, ZIM_S11FT_T_W
REAL, DIMENSION(:),ALLOCATABLE :: ZRE_S22FMS11FT_T_R, ZRE_S22FMS11FT_T_S, ZRE_S22FMS11FT_T_G, ZRE_S22FMS11FT_T_W
INTEGER,DIMENSION(16):: ITMAT
REAL:: ZELEV_RED,ZTC_RED,ZM_RED,ZFW_RED
INTEGER :: JIND
REAL,DIMENSION(7,16) :: KMAT_COEF !matrice contenant tous les coef interpolés
!pour chaque val inf et sup de ELEV_t
REAL :: ZEXPM_MIN, ZEXPM_STEP, ZEXPM_MAX,ZM_MIN
REAL :: ZFW !water fraction inside melting graupel (ZFW=0 for rain, snow and dry graupel). used only with NDIFF=7: Tmatrix
INTEGER :: ILUOUT0,IUNIT
!
! MODIF GAELLE POUR LIMA
LOGICAL :: GLIMA
REAL,DIMENSION(4) :: ZCC_MIN,ZCC_MAX, ZCC_STEP
INTEGER,DIMENSION(4):: INB_CC
REAL, DIMENSION(:),ALLOCATABLE :: ZCC_T_R
REAL :: ZCC_RED
LOGICAL :: GCALC
REAL :: ZCC
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE :: ZM_6D,ZCC_6D
REAL :: ZC
!
REAL :: ZCCR,ZLBR,ZLBEXR,ZDR,ZALPHAR,ZNUR,ZBR
REAL :: ZCCS,ZLBS,ZLBEXS,ZDS,ZALPHAS,ZNUS,ZAS,ZBS,ZCXS
REAL :: ZCCG,ZLBG,ZLBEXG,ZDG,ZALPHAG,ZNUG,ZAG,ZBG,ZCXG
REAL :: ZLBI,ZLBEXI,ZDI,ZALPHAI,ZNUI,ZAI,ZBI
REAL,DIMENSION(:),ALLOCATABLE :: ZRTMIN
WAUTELET Philippe
committed
TYPE(TFILEDATA),POINTER :: TZFILE
!--------------
WAUTELET Philippe
committed
TZFILE => NULL()
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
!
IF (PRESENT(PCR_RAY)) THEN
GLIMA=.TRUE.
ELSE
GLIMA=.FALSE.
ENDIF
!
!
!
ZS11_CARRE_R=0
ZS22_CARRE_R=0
ZRE_S22S11_R=0
ZIM_S22S11_R=0
ZS11_CARRE_I=0
ZS22_CARRE_I=0
ZRE_S22S11_I=0
ZIM_S22S11_I=0
ZS11_CARRE_S=0
ZS22_CARRE_S=0
ZRE_S22S11_S=0
ZIM_S22S11_S=0
ZS11_CARRE_G=0
ZS22_CARRE_G=0
ZRE_S22S11_G=0
ZIM_S22S11_G=0
! Initialisation varibales microphysiques
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
IF (GLIMA) THEN ! LIMA
ZLBR=XLBR_L
ZLBEXR=XLBEXR_L
ZDR=XDR_L
ZALPHAR=XALPHAR_L
ZNUR=XNUR_L
ZBR=XBR_L
ZCCS=XCCS_L
ZCXS=XCXS_L
ZLBS=XLBS_L
ZLBEXS=XLBEXS_L
ZDS=XDS_L
ZALPHAS=XALPHAS_L
ZNUS=XNUS_L
ZAS=XAS_L
ZBS=XBS_L
ZCCG=XCCG_L
ZCXG=XCXG_L
ZLBG=XLBG_L
ZLBEXG=XLBEXG_L
ZDG=XDG_L
ZALPHAG=XALPHAG_L
ZNUG=XNUG_L
ZAG=XAG_L
ZBG=XBG_L
ZLBI=XLBI_L
ZLBEXI=XLBEXI_L
ZDI=XDI_L
ZALPHAI=XALPHAI_L
ZNUI=XNUI_L
ZAI=XAI_L
ZBI=XBI_L
ALLOCATE(ZRTMIN(SIZE(XRTMIN_L)))
ZRTMIN=XRTMIN_L
ELSE ! ICE3
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
ZCCR=XCCR_I
ZLBR=XLBR_I
ZLBEXR=XLBEXR_I
ZDR=XDR_I
ZALPHAR=XALPHAR_I
ZNUR=XNUR_I
ZBR=XBR_I
ZCCS=XCCS_I
ZCXS=XCXS_I
ZLBS=XLBS_I
ZLBEXS=XLBEXS_I
ZDS=XDS_I
ZALPHAS=XALPHAS_I
ZNUS=XNUS_I
ZAS=XAS_I
ZBS=XBS_I
ZCCG=XCCG_I
ZCXG=XCXG_I
ZLBG=XLBG_I
ZLBEXG=XLBEXG_I
ZDG=XDG_I
ZALPHAG=XALPHAG_I
ZNUG=XNUG_I
ZAG=XAG_I
ZBG=XBG_I
ZLBI=XLBI_I
ZLBEXI=XLBEXI_I
ZDI=XDI_I
ZALPHAI=XALPHAI_I
ZNUI=XNUI_I
ZAI=XAI_I
ZBI=XBI_I
ALLOCATE(ZRTMIN(SIZE(XRTMIN_I)))
ZRTMIN=XRTMIN_I
IF (LATT) THEN
IRHOHV=22 !au lieu de 18:
!"ZHH","ZDR","KDP","CSR","ZER","ZEI","ZES","ZEG","VRU"
!"AER","AEI","AES","AEG","AVR","AVI","AVS","AVG","ATR","ATI","ATS","ATG"
ELSE
IRHOHV=10
END IF
IPHIDP=IRHOHV+1
IDELTAHV=IPHIDP+1
IRHR=IDELTAHV+1
IRHS=IRHR+1
IRHG=IRHS+1
IZDA=IRHG+1
IZDS=IZDA+1
IZDG=IZDS+1
IKDR=IZDG+1
IKDS=IKDR+1
IKDG=IKDS+1
INBRAD=SIZE(PT_RAY,1)
IIELV=SIZE(PT_RAY,2)
INBAZIM=SIZE(PT_RAY,3)
INBSTEPMAX=SIZE(PT_RAY,4)
INPTS_H=SIZE(PT_RAY,5)
INPTS_V=SIZE(PT_RAY,6)
!
! Initialisation for radial winds
IF(LFALL) THEN
IF (GLIMA) THEN
ZCR=XCR_L
ZCI=XC_I_L
ZCS=XCS_L
ZCG=XCG_L
ELSE
ZCR=XCR_I
ZCI=XC_I_I
ZCS=XCS_I
ZCG=XCG_I
ELSE
ZCR=0.
ZCI=0.
ZCS=0.
ZCG=0.
! Calculation of nodes and weights for the Gauss-Laguerre quadrature
! for Mie and T-matrix and RG
ALLOCATE(ZX(NPTS_GAULAG),ZW(NPTS_GAULAG)) !NPTS_GAULAG : number of points for the quadrature
CALL GAULAG(NPTS_GAULAG,ZX,ZW)
IVDOP=9 !index of Doppler Velocity (VRU) in ZREFL
WRITE(ILUOUT0,*) "-----------------"
WRITE(ILUOUT0,*) "Radar scattering"
WRITE(ILUOUT0,*) "-----------------"
WRITE(ILUOUT0,*) 'Nombre de variables dans PZE: ',IMAX
ALLOCATE(ZREFL(INBRAD,IIELV,INBAZIM,INBSTEPMAX,INPTS_H,INPTS_V,IMAX))
ZREFL(:,:,:,:,:,:,:)=0.
IF(LATT) THEN
ZREFL(:,:,:,:,:,:,IATR:IATG)=1.
END IF
PZE(:,:,:,:,:)=0.
IF (LATT)THEN
ALLOCATE(ZAELOC(INBRAD,IIELV,INBAZIM,INBSTEPMAX,INPTS_H,INPTS_V,2))
ALLOCATE(ZAETOT(INPTS_H,INPTS_V,2))
ZAELOC(:,:,:,:,:,:,:)=0. ! initialization of attenuation stuff (alpha_e for first gate)
ZAETOT(:,:,:)=1. ! initialization of attenuation stuff (total attenuation)
WRITE(ILUOUT0,*) 'BEFORE LOOP DIFFUSION'
ALLOCATE(ZCONC_BIN(INBRAD,IIELV,INBAZIM,INBSTEPMAX))
ZCONC_BIN(:,:,:,:)=0.
WRITE(ILUOUT0,*) "XCCR:",ZCCR
WRITE(ILUOUT0,*) "XLBR:",ZLBR
WRITE(ILUOUT0,*) "XLBEXR:",ZLBEXR
WRITE(ILUOUT0,*) "XCCS:",ZCCS
WRITE(ILUOUT0,*) "XLBS:",ZLBS
WRITE(ILUOUT0,*) "XLBEXS:",ZLBEXS
WRITE(ILUOUT0,*) "XCCG:",ZCCG
WRITE(ILUOUT0,*) "XLBG:",ZLBG
WRITE(ILUOUT0,*) "XLBEXG:",ZLBEXG
IF (GLIMA .AND. NDIFF==7) THEN
IF (ZALPHAR/=1 .AND. ZNUR /=2.) THEN
WRITE(ILUOUT0,*) " ERROR : TMATRICE TABLE ARE MADE WITH XALPHAR=1 XNUR=2"
WRITE(ILUOUT0,*) " FOR CCLOUD=LIMA. PLEASE CHANGE THIS VALUES OR PROVIDE "
WRITE(ILUOUT0,*) " NEW TMATRICE TABLES "
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING','')
ENDIF
ELSE
IF (ZALPHAR/=1 .AND. ZNUR /=1.) THEN
WRITE(ILUOUT0,*) " ERROR : TMATRICE TABLE ARE MADE WITH XALPHAR=1 XNUR=1"
WRITE(ILUOUT0,*) " FOR CCLOUD=ICE3. PLEASE CHANGE THIS VALUEs OR PROVIDE "
WRITE(ILUOUT0,*) " NEW TMATRICE TABLES "
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING','')
ENDIF
ENDIF
!---------------------------------------------
! LOOP OVER EVERYTHING
!--------------------------------------------
IF(NDIFF==7) THEN
YTAB_TYPE(1)='r'
YTAB_TYPE(2)='s'
YTAB_TYPE(3)='g'
YTAB_TYPE(4)='w'
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
! definition des paramètres de lecture de la table T-matrice
! all mixing ratio
ZEXPM_MIN=-7.
ZEXPM_STEP=0.01
ZEXPM_MAX=-2.
ZM_MIN=10**ZEXPM_MIN
! rain
ZELEV_MIN(1)=0.0
ZELEV_STEP(1)=4.0
ZELEV_MAX(1)=12.0
ZTC_MIN(1)=-20.0
ZTC_STEP(1)=1.0
ZTC_MAX(1)=40.0
ZFW_MIN(1)=0.0
ZFW_STEP(1)=0.1
ZFW_MAX(1)=0.0
IF (GLIMA) THEN
ZCC_MIN(1)=1.8
ZCC_STEP(1)=0.02
ZCC_MAX(1)=6
ELSE
ZCC_MIN(1)=1.
ZCC_STEP(1)=1.
ZCC_MAX(1)=1.
ENDIF
! snow + graupel
ZELEV_MIN(2:3)=0.0
ZELEV_STEP(2:3)=12.0
ZELEV_MAX(2:3)=12.0
ZTC_MIN(2:3)=-70.0
ZTC_STEP(2:3)=1.0
ZTC_MAX(2:3)=10.0
ZFW_MIN(2:3)=0.0
ZFW_STEP(2:3)=0.1
ZFW_MAX(2:3)=0.0
ZCC_MIN(2:3)=1.
ZCC_STEP(2:3)=1.
ZCC_MAX(2:3)=1.
! wet graupel
ZELEV_MIN(4)=0.0
ZELEV_STEP(4)=4.0
ZELEV_MAX(4)=12.0
ZTC_MIN(4)=-10.0
ZTC_STEP(4)=1.0
ZTC_MAX(4)=10.0
ZFW_MIN(4)=0.0
ZFW_STEP(4)=0.1
ZFW_MAX(4)=1.0
ZCC_MIN(4)=1.
ZCC_STEP(4)=1.
ZCC_MAX(4)=1.
DO JT=1,4
INB_ELEV(JT)=NINT((ZELEV_MAX(JT)-ZELEV_MIN(JT))/ZELEV_STEP(JT))+1
INB_TC(JT)=NINT((ZTC_MAX(JT)-ZTC_MIN(JT))/ZTC_STEP(JT))+1
INB_FW(JT)=NINT((ZFW_MAX(JT)-ZFW_MIN(JT))/ZFW_STEP(JT))+1
INB_M=NINT((ZEXPM_MAX-ZEXPM_MIN)/ZEXPM_STEP)+1
INB_CC(JT)=NINT((ZCC_MAX(JT)-ZCC_MIN(JT))/ZCC_STEP(JT))+1
INB_LINE(JT)=INB_ELEV(JT)*INB_TC(JT)*INB_FW(JT)*INB_M*INB_CC(JT)
ENDDO
ENDIF
!---------------------------------------------
! LOOP OVER EVERYTHING
!--------------------------------------------
!============== loop over radars =================
WRITE(ILUOUT0,*) "INBRAD",INBRAD
DO JI=1,INBRAD
WRITE(ILUOUT0,*) "JI",JI
WRITE(ILUOUT0,*) "XLAM_RAD(JI):",XLAM_RAD(JI)
IF(NDIFF==7) THEN ! If T-MATRIX
!---------------------------------------------------------------------------------------------
! 0. LECTURE DES TABLES TMAT POUR PLUIE, NEIGE, GRAUPEL
! en fonction de la bande frequence
!---------------------------------------------------------------------------------------------
IF ( XLAM_RAD(JI)==0.1062) THEN
YBAND='S106.2'
ELSEIF (XLAM_RAD(JI) ==0.0532 ) THEN
YBAND='C053.2'
ELSEIF (XLAM_RAD(JI)==0.0319 ) THEN
YBAND='X031.9'
ELSE
WRITE(ILUOUT0,*) "ERROR RADAR_SCATTERING"
WRITE(ILUOUT0,*) "Tmatrice tables are only available for XLAM_RAD=0.1062"
WRITE(ILUOUT0,*) "or XLAM_RAD=0.053.2 or XLAM_RAD=0.031.8"
WRITE(ILUOUT0,*) "change XLAM_RAD in namelist or compute new tmatrice table"
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING','')
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
ENDIF
!************ fichiers Min Max Pas et Coef Tmat ***********
DO JT=1,4 !types (r, s, g, w)
YTYPE=YTAB_TYPE(JT)
IF (JT .EQ. 1) THEN
IF (GLIMA) THEN
YFILE_COEFINT(JT)='TmatCoefInt_LIMA_'//YBAND//YTYPE
ELSE
YFILE_COEFINT(JT)='TmatCoefInt_ICE3_'//YBAND//YTYPE
ENDIF
ELSE
YFILE_COEFINT(JT)='TmatCoefInt_'//YBAND//YTYPE
ENDIF
YFILE_COEFINT(JT)=TRIM(ADJUSTL(YFILE_COEFINT(JT)))
ENDDO
!lookup tables for rain
ALLOCATE (ZTC_T_R(INB_LINE(1)),ZELEV_T_R(INB_LINE(1)),ZCC_T_R(INB_LINE(1)),ZM_T_R(INB_LINE(1)),&
ZS11_CARRE_T_R(INB_LINE(1)),ZS22_CARRE_T_R(INB_LINE(1)), ZRE_S22S11_T_R(INB_LINE(1)),ZIM_S22S11_T_R(INB_LINE(1)),&
ZRE_S22FMS11FT_T_R(INB_LINE(1)),ZIM_S22FT_T_R(INB_LINE(1)),ZIM_S11FT_T_R(INB_LINE(1)))
!lookup tables for snow
ALLOCATE (ZTC_T_S(INB_LINE(2)),ZELEV_T_S(INB_LINE(2)),ZFW_T_S(INB_LINE(2)),ZM_T_S(INB_LINE(2)),&
ZS11_CARRE_T_S(INB_LINE(2)),ZS22_CARRE_T_S(INB_LINE(2)),ZRE_S22S11_T_S(INB_LINE(2)),ZIM_S22S11_T_S(INB_LINE(2)),&
ZRE_S22FMS11FT_T_S(INB_LINE(2)),ZIM_S22FT_T_S(INB_LINE(2)),ZIM_S11FT_T_S(INB_LINE(2)))
!lookup tables for graupel
ALLOCATE (ZTC_T_G(INB_LINE(3)),ZELEV_T_G(INB_LINE(3)),ZFW_T_G(INB_LINE(3)),ZM_T_G(INB_LINE(3)),&
ZS11_CARRE_T_G(INB_LINE(3)),ZS22_CARRE_T_G(INB_LINE(3)), ZRE_S22S11_T_G(INB_LINE(3)),ZIM_S22S11_T_G(INB_LINE(3)),&
ZRE_S22FMS11FT_T_G(INB_LINE(3)),ZIM_S22FT_T_G(INB_LINE(3)),ZIM_S11FT_T_G(INB_LINE(3)))
!lookup tables for wet graupel
ALLOCATE (ZTC_T_W(INB_LINE(4)),ZELEV_T_W(INB_LINE(4)),ZFW_T_W(INB_LINE(4)),ZM_T_W(INB_LINE(4)),&
ZS11_CARRE_T_W(INB_LINE(4)),ZS22_CARRE_T_W(INB_LINE(4)), ZRE_S22S11_T_W(INB_LINE(4)),ZIM_S22S11_T_W(INB_LINE(4)),&
ZRE_S22FMS11FT_T_W(INB_LINE(4)),ZIM_S22FT_T_W(INB_LINE(4)),ZIM_S11FT_T_W(INB_LINE(4)))
!===== Lecture des tables ===========
6003 FORMAT (E11.4,2X,E9.3,2X,E10.4,2X,E10.4,2X,E12.5,2X,E12.5,2X,&
E12.5,2X,E12.5,2X,E12.5,2X,E12.5,2X,E12.5)
!rain
WAUTELET Philippe
committed
CALL IO_FILE_ADD2LIST(TZFILE,YFILE_COEFINT(1),'TXT','READ')
CALL IO_FILE_OPEN_ll(TZFILE,KRESP=IRESP)
IUNIT = TZFILE%NLU
IF ( IRESP /= 0 ) THEN
WRITE(YMSG,*) "problem opening file ",TRIM(YFILE_COEFINT(1))
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING',YMSG)
ENDIF
ILINE=1
DO WHILE (ILINE .LE. INB_LINE(1))
READ( UNIT=IUNIT,FMT=6003, IOSTAT=IRESP ) ZTC_T_R(ILINE),ZELEV_T_R(ILINE),&
ZCC_T_R(ILINE),ZM_T_R(ILINE),ZS11_CARRE_T_R(ILINE),ZS22_CARRE_T_R(ILINE),ZRE_S22S11_T_R(ILINE),&
ZIM_S22S11_T_R(ILINE),ZRE_S22FMS11FT_T_R(ILINE),ZIM_S22FT_T_R(ILINE),ZIM_S11FT_T_R(ILINE)
ILINE=ILINE+1
ENDDO
WAUTELET Philippe
committed
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
WRITE(ILUOUT0,*) "NLIGNE rain",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
WRITE(ILUOUT0,*) "ZTC_T_R(ILINE),ZELEV_T_R(ILINE),ZCC_T_R(ILINE)",&
ZTC_T_R(ILINE),ZELEV_T_R(ILINE),ZCC_T_R(ILINE)
WRITE(ILUOUT0,*) "ZM_T_R(ILINE),ZS11_CARRE_T_R(ILINE),ZS22_CARRE_T_R(ILINE),ZRE_S22S11_T_R(ILINE)",&
ZM_T_R(ILINE),ZS11_CARRE_T_R(ILINE),ZS22_CARRE_T_R(ILINE),ZRE_S22S11_T_R(ILINE)
WRITE(ILUOUT0,*) "ZIM_S22S11_T_R(ILINE),ZRE_S22FMS11FT_T_R(ILINE),ZIM_S22FT_T_R(ILINE),ZIM_S11FT_T_R(ILINE)",&
ZIM_S22S11_T_R(ILINE),ZRE_S22FMS11FT_T_R(ILINE),ZIM_S22FT_T_R(ILINE),ZIM_S11FT_T_R(ILINE)
!snow
WAUTELET Philippe
committed
CALL IO_FILE_ADD2LIST(TZFILE,YFILE_COEFINT(2),'TXT','READ')
CALL IO_FILE_OPEN_ll(TZFILE,KRESP=IRESP)
IUNIT = TZFILE%NLU
IF ( IRESP /= 0 ) THEN
WRITE(YMSG,*) "problem opening file ",TRIM(YFILE_COEFINT(2))
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING',YMSG)
ENDIF
ILINE=1
DO WHILE (ILINE .LE. INB_LINE(2))
READ( UNIT=IUNIT,FMT=6003, IOSTAT=IRESP ) ZTC_T_S(ILINE),ZELEV_T_S(ILINE),&
ZFW_T_S(ILINE),ZM_T_S(ILINE),ZS11_CARRE_T_S(ILINE),ZS22_CARRE_T_S(ILINE),ZRE_S22S11_T_S(ILINE),&
ZIM_S22S11_T_S(ILINE),ZRE_S22FMS11FT_T_S(ILINE),ZIM_S22FT_T_S(ILINE),ZIM_S11FT_T_S(ILINE)
ILINE=ILINE+1
ENDDO
WAUTELET Philippe
committed
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
WRITE(ILUOUT0,*) "NLIGNE snow",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
WRITE(ILUOUT0,*) "ZTC_T_S(ILINE),ZELEV_T_S(ILINE),ZFW_T_S(ILINE)",&
ZTC_T_S(ILINE),ZELEV_T_S(ILINE),ZFW_T_S(ILINE)
WRITE(ILUOUT0,*) "ZM_T_S(ILINE),ZS11_CARRE_T_S(ILINE),ZS22_CARRE_T_S(ILINE),ZRE_S22S11_T_S(ILINE)",&
ZM_T_S(ILINE),ZS11_CARRE_T_S(ILINE),ZS22_CARRE_T_S(ILINE),ZRE_S22S11_T_S(ILINE)
WRITE(ILUOUT0,*) "ZIM_S22S11_T_S(ILINE),ZRE_S22FMS11FT_T_S(ILINE),ZIM_S22FT_T_S(ILINE),ZIM_S11FT_T_S(ILINE)",&
ZIM_S22S11_T_S(ILINE),ZRE_S22FMS11FT_T_S(ILINE),ZIM_S22FT_T_S(ILINE),ZIM_S11FT_T_S(ILINE)
!graupel
WAUTELET Philippe
committed
CALL IO_FILE_ADD2LIST(TZFILE,YFILE_COEFINT(3),'TXT','READ')
CALL IO_FILE_OPEN_ll(TZFILE,KRESP=IRESP)
IUNIT = TZFILE%NLU
IF ( IRESP /= 0 ) THEN
WRITE(YMSG,*) "problem opening file ",TRIM(YFILE_COEFINT(3))
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING',YMSG)
ENDIF
ILINE=1
DO WHILE (ILINE .LE. INB_LINE(3))
READ( UNIT=IUNIT, FMT=6003,IOSTAT=IRESP ) ZTC_T_G(ILINE),ZELEV_T_G(ILINE),&
ZFW_T_G(ILINE),ZM_T_G(ILINE),ZS11_CARRE_T_G(ILINE),ZS22_CARRE_T_G(ILINE),ZRE_S22S11_T_G(ILINE),&
ZIM_S22S11_T_G(ILINE),ZRE_S22FMS11FT_T_G(ILINE),ZIM_S22FT_T_G(ILINE),ZIM_S11FT_T_G(ILINE)
ILINE=ILINE+1
ENDDO
WAUTELET Philippe
committed
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
WRITE(ILUOUT0,*) "NLIGNE graupel",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
WRITE(ILUOUT0,*) "ZTC_T_G(ILINE),ZELEV_T_G(ILINE)",&
ZTC_T_G(ILINE),ZELEV_T_G(ILINE)
WRITE(ILUOUT0,*) "ZM_T_G(ILINE),ZS11_CARRE_T_G(ILINE),ZS22_CARRE_T_G(ILINE),ZRE_S22S11_T_G(ILINE)",&
ZM_T_G(ILINE),ZS11_CARRE_T_G(ILINE),ZS22_CARRE_T_G(ILINE),ZRE_S22S11_T_G(ILINE)
WRITE(ILUOUT0,*) "ZIM_S22S11_T_G(ILINE),ZRE_S22FMS11FT_T_G(ILINE),ZIM_S22FT_T_G(ILINE),ZIM_S11FT_T_G(ILINE)",&
ZIM_S22S11_T_G(ILINE),ZRE_S22FMS11FT_T_G(ILINE),ZIM_S22FT_T_G(ILINE),ZIM_S11FT_T_G(ILINE)
!wet graupel
WAUTELET Philippe
committed
CALL IO_FILE_ADD2LIST(TZFILE,YFILE_COEFINT(4),'TXT','READ')
CALL IO_FILE_OPEN_ll(TZFILE,KRESP=IRESP)
IUNIT = TZFILE%NLU
IF ( IRESP /= 0 ) THEN
WRITE(YMSG,*) "problem opening file ",TRIM(YFILE_COEFINT(4))
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING',YMSG)
ENDIF
ILINE=1
DO WHILE (ILINE .LE. INB_LINE(4))
READ( UNIT=IUNIT, FMT=6003,IOSTAT=IRESP ) ZTC_T_W(ILINE),ZELEV_T_W(ILINE),&
ZFW_T_W(ILINE),ZM_T_W(ILINE),ZS11_CARRE_T_W(ILINE),ZS22_CARRE_T_W(ILINE),ZRE_S22S11_T_W(ILINE),&
ZIM_S22S11_T_W(ILINE),ZRE_S22FMS11FT_T_W(ILINE),ZIM_S22FT_T_W(ILINE),ZIM_S11FT_T_W(ILINE)
ILINE=ILINE+1
ENDDO
WAUTELET Philippe
committed
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
WRITE(ILUOUT0,*) "NLIGNE wet graupel",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
WRITE(ILUOUT0,*) "ZTC_T_W(ILINE),ZELEV_T_W(ILINE)", ZTC_T_W(ILINE),ZELEV_T_W(ILINE)
WRITE(ILUOUT0,*) "ZM_T_W(ILINE),ZS11_CARRE_T_W(ILINE),ZS22_CARRE_T_W(ILINE),ZRE_S22S11_T_W(ILINE)",&
ZM_T_W(ILINE),ZS11_CARRE_T_W(ILINE),ZS22_CARRE_T_W(ILINE),ZRE_S22S11_T_W(ILINE)
WRITE(ILUOUT0,*) "ZIM_S22S11_T_W(ILINE),ZRE_S22FMS11FT_T_W(ILINE),ZIM_S22FT_T_W(ILINE),ZIM_S11FT_T_W(ILINE)",&
ZIM_S22S11_T_W(ILINE),ZRE_S22FMS11FT_T_W(ILINE),ZIM_S22FT_T_W(ILINE),ZIM_S11FT_T_W(ILINE)
ENDIF !END IF T-MATRIX => END OF LOOKUP TABLE READING
!============== loop over elevations =================
IEL=NBELEV(JI)
WRITE(ILUOUT0,*) "NBELEV(JI)",NBELEV(JI)
WRITE(ILUOUT0,*) "INPTS_V",INPTS_V
DO JEL=1,IEL
WRITE(ILUOUT0,*) "JEL",JEL
JL=1
JV=1
WRITE(ILUOUT0,*) "JL,JV",JL,JV
WRITE(ILUOUT0,*) "PELEV(JI,JEL,JL,JV)*180./XPI",PELEV(JI,JEL,JL,JV)*180./XPI
JL=INBSTEPMAX
JV=INPTS_V
WRITE(ILUOUT0,*) "JL,JV",JL,JV
WRITE(ILUOUT0,*) "PELEV(JI,JEL,JL,JV)*180./XPI",PELEV(JI,JEL,JL,JV)*180./XPI
!============== loop over azimuths =================
DO JAZ=1,INBAZIM
DO JH=1,INPTS_H !horizontal discretization of the beam
DO JV=1,INPTS_V ! vertical discretization (we go down to check partial masks)
IF(LATT) THEN
ZAERINT=1.
ZAVRINT=1.
ZAEIINT=1.
ZAESINT=1.
ZAVSINT=1.
ZAEGINT=1.
ZAVGINT=1.
END IF
!Loop over the ranges for one azimuth. If the range is masked, the reflectivity for all the consecutive ranges is set to 0
LPART_MASK=.FALSE.
LOOPJL: DO JL=1,INBSTEPMAX
IF(LPART_MASK) THEN ! THIS RAY IS MASKED
ZREFL(JI,JEL,JAZ,JL:INBSTEPMAX,JH,JV,1)=0.
EXIT LOOPJL
ELSE
! if not underground or outside of the MESO-NH domain (PT_RAY : temperature interpolated along the rays)
IF(PT_RAY(JI,JEL,JAZ,JL,JH,JV) /= -XUNDEF) THEN
!
!---------------------------------------------------------------------------------------------------
!* 2. RAINDROPS
! ---------
!
IF(SIZE(PR_RAY,1) > 0) THEN ! if PR_RAY is available for at least one radar
!contenu en hydrometeore
ZM=PRHODREF_RAY(JI,JEL,JAZ,JL,JH,JV)*PR_RAY(JI,JEL,JAZ,JL,JH,JV)
IF (GLIMA) ZCC=PRHODREF_RAY(JI,JEL,JAZ,JL,JH,JV)*PCR_RAY(JI,JEL,JAZ,JL,JH,JV)
!ZM_MIN : min value for rain content (10**-7 <=> Z=-26 dBZ)mixing ratio
IF (GLIMA) THEN
GCALC=((ZM > ZM_MIN).AND.(ZCC > 10**ZCC_MIN(1)))
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
GCALC=(ZM > ZM_MIN)
ENDIF
IF(GCALC ) THEN
!calculation of the dielectrique constant (permittitivité relative)
! for liquid water from function QEPSW
!(defined in mode_fscatter.f90 => equation 3.6 p 64)
YTYPE='r'
ZQMW=SQRT(QEPSW(PT_RAY(JI,JEL,JAZ,JL,JH,JV),XLIGHTSPEED/XLAM_RAD(JI)))
!ZLBDA : slope distribution parameter (equation 2.6 p 23)
IF (GLIMA) THEN
ZLBDA=( ZLBR*ZCC / ZM )**ZLBEXR
ELSE
ZLBDA=ZLBR*(ZM)**ZLBEXR
ENDIF
ZQK=(ZQMW**2-1.)/(ZQMW**2+2.) !dielectric factor (3.43 p 56)
ZFW=0 !Liquid water fraction (only for melting graupel => 0 for rain)
!compteur=compteur+1
!---------------------------------------------------
! ------------ DIFFUSION --------------
!---------------------------------------------------
!******************************* NDIFF=0 or 4 *********************************
IF(NDIFF==0.OR.NDIFF==4) THEN ! Rayleigh
!ZREFLOC(1:2) : Zh et Zv = int(sigma(D)*N(D)) (eq 1.6 p 16)
!with N(D) formulation (eq 2.2 p 23) and sigma Rayleigh (3.41 p 55)
!MOMG : gamma function defined in mong.f90
!XCCR = 1.E7; XLBEXR = -0.25! Marshall-Palmer law (radar_rain_ice.f90)
!ZCXR : -1 (Xi coeff in equation 2.3 p 23)
ZREFLOC(1:2)=1.E18*ZCCR*ZLBDA**(ZCXR-6.)*MOMG(ZALPHAR,ZNUR,6.)
IF(LWREFL) THEN ! weighting by reflectivities
!ZREFL(...,IVDOP)=radial velocity (IVDOP=9), weighted by reflectivity and
!taking into account raindrops fall velocity (ZCR = 842, XDR = 0.8 -> 2.8 p23 et 2.1 p24)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=-ZCR*SIN(PELEV(JI,JEL,JL,JV)) &
*1.E18*ZCCR*ZLBDA**(ZCXR-6.-ZDR)*MOMG(ZALPHAR,ZNUR,6.+ZDR)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZCCR*ZLBDA**ZCXR ! N0j of equation 2.3 p23 (density of particules)
!projection of fall velocity only
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=-ZCR*SIN(PELEV(JI,JEL,JL,JV)) &
*ZCCR*ZLBDA**(ZCXR-ZDR)*MOMG(ZALPHAR,ZNUR,ZDR)
END IF ! end weighting by reflectivities
IF(LATT) THEN ! Calculation of Extinction coefficient
IF(NDIFF==0) THEN ! Rayleigh 3rd order : calculation from equations
! 3.39 p55 : extinction coeff = int(extinction_section(D) * N(D))
! 2.2 and 2.3 p23: simplification of int(D**p * N(D)) and N0j
! 3.42 p57 : extinction_section(D)
ZAETMP(:)=ZCCR*ZLBDA**ZCXR*(XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
*MOMG(ZALPHAR,ZNUR,ZBR)/ZLBDA**ZBR)
ELSE ! Rayleigh 6th order ! eq 3.52 p 58 for extinction coefficient
ZAETMP(:)=ZCCR*ZLBDA**ZCXR*(XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
*MOMG(ZALPHAR,ZNUR,ZBR)/ZLBDA**ZBR &
+XPI**4/15./XLAM_RAD(JI)**3*AIMAG(ZQK**2*(ZQMW**4+27.*ZQMW**2+38.) &
/(2.*ZQMW**2+3.))*MOMG(ZALPHAR,ZNUR,5.*ZBR/3.)/ZLBDA**(5.*ZBR/3.)&
+2.*XPI**5/3. /XLAM_RAD(JI)**4*REAL(ZQK**2) &
*MOMG(ZALPHAR,ZNUR,2.*ZBR) /ZLBDA**(2.*ZBR))
END IF
END IF ! end IF(LATT)
ZRE_S22S11_R=0
ZIM_S22S11_R=0
ZS22_CARRE_R=0
ZS11_CARRE_R=0
!******************************* NDIFF==7 ************************************
ELSE IF(NDIFF==7) THEN !T-matrix
ZREFLOC(:)=0
IF(LATT) ZAETMP(:)=0
IF (GLIMA) THEN
CALL CALC_KTMAT_LIMA(PELEV(JI,JEL,JL,JV),&
PT_RAY(JI,JEL,JAZ,JL,JH,JV),ZCC,ZM,&
ZELEV_MIN(1),ZELEV_MAX(1),ZELEV_STEP(1),&
ZTC_MIN(1),ZTC_MAX(1),ZTC_STEP(1),&
ZCC_MIN(1),ZCC_MAX(1),ZCC_STEP(1),&
ZEXPM_MIN,ZEXPM_MAX,ZEXPM_STEP,&
ITMAT,ZELEV_RED,ZTC_RED,ZCC_RED,ZM_RED)
ELSE
CALL CALC_KTMAT(PELEV(JI,JEL,JL,JV),&
PT_RAY(JI,JEL,JAZ,JL,JH,JV),ZFW,ZM,&
ZELEV_MIN(1),ZELEV_MAX(1),ZELEV_STEP(1),&
ZTC_MIN(1),ZTC_MAX(1),ZTC_STEP(1),&
ZFW_MIN(1),ZFW_MAX(1),ZFW_STEP(1),&
ZEXPM_MIN,ZEXPM_MAX,ZEXPM_STEP,&
ITMAT,ZELEV_RED,ZTC_RED,ZFW_RED,ZM_RED)
ENDIF
IF (ITMAT(1) .NE. -NUNDEF) THEN
DO JIND=1,SIZE(KMAT_COEF,2),1
KMAT_COEF(1,JIND)=ZS11_CARRE_T_R(ITMAT(JIND))
KMAT_COEF(2,JIND)=ZS22_CARRE_T_R(ITMAT(JIND))
KMAT_COEF(3,JIND)=ZRE_S22S11_T_R(ITMAT(JIND))
KMAT_COEF(4,JIND)=ZIM_S22S11_T_R(ITMAT(JIND))
KMAT_COEF(5,JIND)=ZRE_S22FMS11FT_T_R(ITMAT(JIND))
KMAT_COEF(6,JIND)=ZIM_S22FT_T_R(ITMAT(JIND))
KMAT_COEF(7,JIND)=ZIM_S11FT_T_R(ITMAT(JIND))
ENDDO
IF (GLIMA) THEN
CALL INTERPOL(ZELEV_RED,ZTC_RED,ZCC_RED,ZM_RED,KMAT_COEF,ZS11_CARRE_R,ZS22_CARRE_R,&
ZRE_S22S11_R,ZIM_S22S11_R,ZRE_S22FMS11F,ZIM_S22FT,ZIM_S11FT)
ELSE
CALL INTERPOL(ZELEV_RED,ZTC_RED,ZFW_RED,ZM_RED,KMAT_COEF,ZS11_CARRE_R,ZS22_CARRE_R,&
ZRE_S22S11_R,ZIM_S22S11_R,ZRE_S22FMS11F,ZIM_S22FT,ZIM_S11FT)
ENDIF
ELSE
ZS11_CARRE_R=0
ZS22_CARRE_R=0
ZRE_S22S11_R=0
ZIM_S22S11_R=0
ZRE_S22FMS11F=0
ZIM_S22FT=0
ZIM_S11FT=0
END IF
ZREFLOC(1)=1.E18*(XLAM_RAD(JI))**4/(XPI**5*.93)*4*XPI*ZS22_CARRE_R
ZREFLOC(2)=1.E18*(XLAM_RAD(JI))**4/(XPI**5*.93)*4*XPI*ZS11_CARRE_R
ZREFLOC(3)=180.E3/XPI*XLAM_RAD(JI)*ZRE_S22FMS11F
IF (GLIMA) THEN
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
-ZCR*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(1) &
*1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCC/4./ZLBDA**(2+ZDR)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
-ZCR*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(1) &
*1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCCR/4./ZLBDA**(3+ZDR)
ENDIF
IF(LATT) THEN
ZAETMP(1)=ZIM_S22FT*XLAM_RAD(JI)*2
ZAETMP(2)=ZIM_S11FT*XLAM_RAD(JI)*2
END IF
!******************************* NDIFF=1 or 3 *********************************
! Gauss Laguerre integration
ELSE ! MIE OR T-MATRIX OR RAYLEIGH FOR ELLIPSOIDES
ZREFLOC(:)=0.
IF(LATT) ZAETMP(:)=0.
DO JJ=1,NPTS_GAULAG ! ****** Gauss-Laguerre quadrature
SELECT CASE(NDIFF)
CASE(1) ! *************** NDIFF=1 MIE *****************
! subroutine BHMIE defined in mode_fscatter.f90
! calculate extinction coefficient ZQEXT(1),scattering : ZQSCA
! and backscattering ZQBACK(1) on the horizontal plan (spheroid)
CALL BHMIE(XPI/XLAM_RAD(JI)*ZX(JJ)/ZLBDA,ZQMW,ZQEXT(1),ZQSCA,ZQBACK(1))
ZQBACK(2)=ZQBACK(1) !=> same because sphere
ZQEXT(2)=ZQEXT(1) ! modif Clotilde 23/04/2012
ZQBACK(3)=0. !=> 0 because sphere
CASE(3) !****************** NDIFF==3 RG RAYLEIGH FOR ELLIPSOIDES ***********************
IF(ARF(ZX(JJ)/ZLBDA)==1.) THEN
ZLB=1./3.
ELSE
ZLB=1./(ARF(ZX(JJ)/ZLBDA))**2-1. ! f**2
ZLB=(1.+ZLB)/ZLB*(1.-ATAN(SQRT(ZLB))/SQRT(ZLB)) ! lambda_b
IF(ZX(JJ)/ZLBDA>16.61E-3) PRINT*, 'Negative axis ratio; reduce NPTS_GAULAG.'
END IF
! equation 3.44 p 56 (ZX**4 instead of ZX**6 but ZQBACK is multiplied after by ZX**2)
ZQBACK(1)=4.*(XPI/XLAM_RAD(JI)*ZX(JJ)/ZLBDA)**4&
*ABS((ZQMW**2-1.)/3./(1.+.5*(1.-ZLB)*(ZQMW**2-1.)))**2
! equation 3.45 p 56
ZQBACK(2)=4.*(XPI/XLAM_RAD(JI)*ZX(JJ)/ZLBDA)**4*ABS((ZQMW**2-1.)/3.*&
(SIN(PELEV(JI,JEL,JL,JV))**2/(1.+.5*(1.-ZLB)*(ZQMW**2-1.))+& ! PELEV=PI+THETA_I
COS(PELEV(JI,JEL,JL,JV))**2/(1.+ZLB*(ZQMW**2-1.))) )**2 !
! KDP from equation 3.49
ZQBACK(3)=ZX(JJ)/ZLBDA**3*REAL((ZQMW**2-1.)**2*(3.*ZLB-1.)/(2.+(ZQMW**2-1.)*(ZLB+1.) &
+ZLB*(1.-ZLB)*(ZQMW**2-1.)**2))
IF(LATT) THEN
! equations 3.48 and 3.49 p57
ZQEXT(1)=4.*(XPI/XLAM_RAD(JI)*ZX(JJ)/ZLBDA)*AIMAG((ZQMW**2-1.)/3./(1.+.5*(1.-ZLB)*(ZQMW**2-1.)))
ZQEXT(2)=4.*(XPI/XLAM_RAD(JI)*ZX(JJ)/ZLBDA)*AIMAG((ZQMW**2-1.)/3.*&
(SIN(PELEV(JI,JEL,JL,JV))**2/(1.+.5*(1.-ZLB)*(ZQMW**2-1.))+& ! PELEV=PI+THETA_I
COS(PELEV(JI,JEL,JL,JV))**2/(1.+ZLB*(ZQMW**2-1.))))
END IF
END SELECT !end SELECT NDIFF
!incrementation of the reflectivity and Kdp(1,2,3,4 for Zh, Zv, )
!with the backscattering coefficients for each point of the GAULAG distribution
! or each diameter D
ZREFLOC(1:3)=ZREFLOC(1:3)+ZQBACK(1:3)*ZX(JJ)**2*ZW(JJ)
ZREFLOC(4)=ZREFLOC(4)+ZQBACK(1)*ZX(JJ)**(2+ZDR)*ZW(JJ)
!same for attenuation with extinction coefficient
IF(LATT) ZAETMP(:)=ZAETMP(:)+ZQEXT(:)*ZX(JJ)**2*ZW(JJ)
END DO ! ****** end loop Gauss-Laguerre quadrature
ZREFLOC(1:2)=1.E18*ZREFLOC(1:2)*(XLAM_RAD(JI)/XPI)**4/.93*ZCCR/4./ZLBDA**3
ZREFLOC(3)=ZREFLOC(3)*XPI**2/6./XLAM_RAD(JI)*ZCCR/ZLBDA &
*180.E3/XPI ! (in deg/km)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
-ZCR*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(4) &
*1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCCR/4./ZLBDA**(3+ZDR)
!********* for all cases with Gauss-Laguerre integration
ZRE_S22S11_R=0
ZIM_S22S11_R=0
ZS22_CARRE_R=0
ZS11_CARRE_R=0
IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZCCR*ZLBDA**(ZCXR-2.*ZBR/3.)/(4.*GAMMA(ZNUR))
END IF ! ****************** End if for each type of diffusion ************************
!incrementation of ZHH, ZDR and KDP
ZREFL(JI,JEL,JAZ,JL,JH,JV,1:3)=ZREFL(JI,JEL,JAZ,JL,JH,JV,1:3)+ZREFLOC(1:3)
! ZER (Z due to raindrops)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZER)=ZREFLOC(1)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDA)=ZREFLOC(2) !Zvv for ZDR due to rain
ZREFL(JI,JEL,JAZ,JL,JH,JV,IKDR)=ZREFLOC(3) !Zvv for ZDR due to rain
! RhoHV due to rain
IF (ZS22_CARRE_R*ZS11_CARRE_R .GT. 0) THEN
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHR)=SQRT(ZRE_S22S11_R**2+ZIM_S22S11_R**2)/SQRT(ZS22_CARRE_R*ZS11_CARRE_R)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHR)=1
END IF
IF(LATT) THEN
ZAELOC(JI,JEL,JAZ,JL,JH,JV,:)=ZAETMP(:) ! specific attenuation due to rain
ZREFL(JI,JEL,JAZ,JL,JH,JV,IAER)=ZAETMP(1)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IAVR)=ZAETMP(2)
! for ranges over 1, correction of attenuation on reflectivity due to rain
IF(JL>1) THEN
ZAERINT=ZAERINT*EXP(-2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IAER)*XSTEP_RAD)
ZAVRINT=ZAVRINT*EXP(-2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IAVR)*XSTEP_RAD)
END IF
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZER)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IZER)*ZAERINT ! Z_r attenuated
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDA)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDA)*ZAVRINT ! ZVr attenuated
END IF !end IF(LATT)
END IF
! mimimum rainwater mixing ratio
! Total attenuation even if no hydrometeors (equation 1.7 p 17)