Newer
Older
WAUTELET Philippe
committed
!MNH_LIC Copyright 2004-2019 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
WAUTELET Philippe
committed
!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,PH_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 concentration interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),OPTIONAL,INTENT(IN) :: PH_RAY ! hail 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,PH_RAY)
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
84
! ##############################
!
!!**** *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, ONLY: TFILEDATA
USE MODD_PARAM_ICE, ONLY: LSNOW_T_I=>LSNOW_T
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,XNS_I=>XNS,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,&
XALPHAH_I=>XALPHAH,XNUH_I=>XNUH,XDH_I=>XDH,XLBEXH_I=>XLBEXH,&
XLBH_I=>XLBH,XCCH_I=>XCCH,XAH_I=>XAH,XBH_I=>XBH,XCXH_I=>XCXH,XCH_I=>XCH,&
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, &
XLBDAS_MAX_I=>XLBDAS_MAX,XLBDAS_MIN_I=>XLBDAS_MIN,XTRANS_MP_GAMMAS_I=>XTRANS_MP_GAMMAS
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,XNS_L=>XNS,XAS_L=>XAS,XBS_L=>XBS,&
XCXS_L=>XCXS,XCS_L=>XCS,&
XLBDAS_MAX_L=>XLBDAS_MAX,XLBDAS_MIN_L=>XLBDAS_MIN,XTRANS_MP_GAMMAS_L=>XTRANS_MP_GAMMAS
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,&
VIE Benoît
committed
XRTMIN_L=>XRTMIN, LSNOW_T_L=>LSNOW_T
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_TMAT
!
USE MODE_ARF
USE MODE_FSCATTER
USE MODE_READTMAT
USE MODI_GAMMA, ONLY:GAMMA
!
USE MODD_LUNIT
WAUTELET Philippe
committed
USE MODE_IO_FILE, ONLY: IO_File_close, IO_File_open
USE MODE_IO_MANAGE_STRUCT, ONLY: IO_File_add2list
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 concentration interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),OPTIONAL,INTENT(IN) :: PH_RAY ! hail 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,ZAEHINT ! total attenuation horizontal
REAL :: ZAVRINT,ZAVSINT,ZAVGINT,ZAVHINT ! 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
VIE Benoît
committed
REAL :: ZN ! Number concentration
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_H,ZS22_CARRE_H,ZRE_S22S11_H,ZIM_S22S11_H
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,ZCH ! coefficients to take into account fall speeds when simulating Doppler winds
REAL, DIMENSION(:,:,:,:),ALLOCATABLE :: ZCONC_BIN
LOGICAL :: LPART_MASK ! indicates a partial mask along the beam
!
INTEGER :: IZER,IZEI,IZES,IZEG
INTEGER :: IVDOP,IRHV,IPDP,IDHV
INTEGER :: IAER,IAEI,IAES,IAEG
INTEGER :: IAVR,IAVI,IAVS,IAVG
INTEGER :: IATR,IATI,IATS,IATG
INTEGER :: IRHR, IRHS, IRHG, IZDA, IZDS, IZDG, IKDR, IKDS, IKDG
INTEGER :: IZEH, IRHH,IKDH,IZDH ! hail
INTEGER :: IAEH,IAVH,IATH
!
!for ZSNR threshold
REAL ::ZDISTRAD,ZSNR,ZSNR_R,ZSNR_S,ZSNR_I,ZSNR_G,ZSNR_H,ZZHH,ZZE_R,ZZE_I,ZZE_S,ZZE_G,ZZE_H
LOGICAL :: GTHRESHOLD_V, GTHRESHOLD_Z,GTHRESHOLD_ZR,GTHRESHOLD_ZI,GTHRESHOLD_ZS,GTHRESHOLD_ZG,GTHRESHOLD_ZH
!--------- TO READ T-MATRIX TABLE --------
CHARACTER(LEN=6) :: YBAND
CHARACTER(LEN=1) ::YTYPE
CHARACTER(LEN=1),DIMENSION(5) :: YTAB_TYPE
CHARACTER(LEN=25),DIMENSION(5) :: YFILE_COEFINT
REAL,DIMENSION(5) :: 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(5) :: 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
REAL, DIMENSION(:),ALLOCATABLE :: ZTC_T_H ,ZELEV_T_H ,ZFW_T_H,ZM_T_H,ZS11_CARRE_T_H,ZS22_CARRE_T_H,ZRE_S22S11_T_H
REAL, DIMENSION(:),ALLOCATABLE :: ZIM_S22S11_T_H,ZIM_S22FT_T_H,ZIM_S11FT_T_H,ZRE_S22FMS11FT_T_H
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,GHAIL
REAL,DIMENSION(5) :: ZCC_MIN,ZCC_MAX, ZCC_STEP
INTEGER,DIMENSION(5):: 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,ZNS
REAL :: ZCCG,ZLBG,ZLBEXG,ZDG,ZALPHAG,ZNUG,ZAG,ZBG,ZCXG
REAL :: ZCCH,ZLBH,ZLBEXH,ZDH,ZALPHAH,ZNUH,ZAH,ZBH,ZCXH
REAL :: ZLBI,ZLBEXI,ZDI,ZALPHAI,ZNUI,ZAI,ZBI
REAL,DIMENSION(:),ALLOCATABLE :: ZRTMIN
WAUTELET Philippe
committed
TYPE(TFILEDATA),POINTER :: TZFILE
!--------------
WAUTELET Philippe
committed
TZFILE => NULL()
!
IF (PRESENT(PCR_RAY)) THEN
GLIMA=.TRUE.
ELSE
GLIMA=.FALSE.
ENDIF
IF (PRESENT(PH_RAY)) THEN
GHAIL=.TRUE.
ELSE
GHAIL=.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
ZS11_CARRE_H=0
ZS22_CARRE_H=0
ZRE_S22S11_H=0
ZIM_S22S11_H=0
! Initialisation varibales microphysiques
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
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 (GHAIL) THEN
ZCCH=XCCH_I
ZCXH=XCXH_I
ZLBH=XLBH_I
ZLBEXH=XLBEXH_I
ZDH=XDH_I
ZALPHAH=XALPHAH_I
ZNUH=XNUH_I
ZAH=XAH_I
ZBH=XBH_I
ENDIF
!
! initialisation of refractivity indices
! 1 : ZHH
! 2 : ZDR
! 3 : KDP
! 4 : CSR
IZER=5 ! ZER
IZEI=IZER+1 ! ZEI
IZES=IZEI+1 ! ZES
IZEG=IZES+1 ! ZEG
IF (GHAIL) THEN
IZEH=IZEG+1 !ZEH
IVDOP=IZEH+1 !VRU
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
IF (LATT) THEN
IF (GHAIL) THEN
IAER=IVDOP+1
IAEI=IAER+1
IAES=IAEI+1
IAEG=IAES+1
IAEH=IAEG+1
IAVR=IAEH+1
IAVI=IAVR+1
IAVS=IAVI+1
IAVG=IAVS+1
IAVH=IAVG+1
IATR=IAVH+1
IATI=IATR+1
IATS=IATI+1
IATG=IATS+1
IATH=IATG+1
IRHV=IATH+1
ELSE
IAER=IVDOP+1
IAEI=IAER+1
IAES=IAEI+1
IAEG=IAES+1
IAVR=IAEG+1
IAVI=IAVR+1
IAVS=IAVI+1
IAVG=IAVS+1
IATR=IAVG+1
IATI=IATR+1
IATS=IATI+1
IATG=IATS+1
IRHV=IATG+1
ENDIF
ELSE
IRHV=IVDOP+1
ENDIF
IPDP=IRHV+1
IDHV=IPDP+1
IRHR=IDHV+1
IRHS=IRHR+1
IRHG=IRHS+1
IF (GHAIL) THEN
IRHH=IRHG+1
IZDA=IRHH+1
ELSE
IZDA=IRHG+1
ENDIF
IZDS=IZDA+1
IZDG=IZDS+1
IF (GHAIL) THEN
IZDH=IZDG+1
IKDR=IZDH+1
ELSE
IKDR=IZDG+1
ENDIF
IKDS=IKDR+1
IKDG=IKDS+1
IF (GHAIL) THEN
IKDH=IKDG+1
ENDIF
!
!
!
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)
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.
IF (GHAIL) ZREFL(:,:,:,:,:,:,IATH)=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 (GHAIL) THEN
WRITE(ILUOUT0,*) "XCCH:",ZCCH
WRITE(ILUOUT0,*) "XLBH:",ZLBH
WRITE(ILUOUT0,*) "XLBEXH:",ZLBEXH
ENDIF
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'
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
! 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.
! hail
ZELEV_MIN(5)=0.0
ZELEV_STEP(5)=4.0
ZELEV_MAX(5)=12.0
ZTC_MIN(5)=-20.0
ZTC_STEP(5)=1.0
ZTC_MAX(5)=30.0
ZFW_MIN(5)=0.
ZFW_STEP(5)=0.1
ZFW_MAX(5)=0.0
ZCC_MIN(5)=1.
ZCC_STEP(5)=1.
ZCC_MAX(5)=1.
DO JT=1,5
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
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.0532 or XLAM_RAD=0.0319"
WRITE(ILUOUT0,*) "change XLAM_RAD in namelist or compute new tmatrice table"
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING','')
ENDIF
!************ fichiers Min Max Pas et Coef Tmat ***********
DO JT=1,5 !types (r, s, g, w, h)
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
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)))
IF (GHAIL) THEN
!lookup tables for hail
ALLOCATE (ZTC_T_H(INB_LINE(5)),ZELEV_T_H(INB_LINE(5)),ZFW_T_H(INB_LINE(5)),ZM_T_H(INB_LINE(5)),&
ZS11_CARRE_T_H(INB_LINE(5)),ZS22_CARRE_T_H(INB_LINE(5)), ZRE_S22S11_T_H(INB_LINE(5)),ZIM_S22S11_T_H(INB_LINE(5)),&
ZRE_S22FMS11FT_T_H(INB_LINE(5)),ZIM_S22FT_T_H(INB_LINE(5)),ZIM_S11FT_T_H(INB_LINE(5)))
ENDIF
!===== 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(TZFILE,KRESP=IRESP)
WAUTELET Philippe
committed
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(TZFILE)
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(TZFILE,KRESP=IRESP)
WAUTELET Philippe
committed
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(TZFILE)
WAUTELET Philippe
committed
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(TZFILE,KRESP=IRESP)
WAUTELET Philippe
committed
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(TZFILE)
WAUTELET Philippe
committed
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(TZFILE,KRESP=IRESP)
WAUTELET Philippe
committed
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(TZFILE)
WAUTELET Philippe
committed
TZFILE => NULL()
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)
!hail
IF (GHAIL) THEN
WAUTELET Philippe
committed
CALL IO_File_add2list(TZFILE,YFILE_COEFINT(5),'TXT','READ')
CALL IO_File_open(TZFILE,KRESP=IRESP)
IUNIT = TZFILE%NLU
IF ( IRESP /= 0 ) THEN
WRITE(YMSG,*) "problem opening file ",TRIM(YFILE_COEFINT(5))
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING',YMSG)
ENDIF
ILINE=1
DO WHILE (ILINE .LE. INB_LINE(5))
READ( UNIT=IUNIT, FMT=6003,IOSTAT=IRESP ) ZTC_T_H(ILINE),ZELEV_T_H(ILINE),&
ZFW_T_H(ILINE),ZM_T_H(ILINE),ZS11_CARRE_T_H(ILINE),ZS22_CARRE_T_H(ILINE),ZRE_S22S11_T_H(ILINE),&
ZIM_S22S11_T_H(ILINE),ZRE_S22FMS11FT_T_H(ILINE),ZIM_S22FT_T_H(ILINE),ZIM_S11FT_T_H(ILINE)
ILINE=ILINE+1
ENDDO
WAUTELET Philippe
committed
CALL IO_File_close(TZFILE)
TZFILE => NULL()
WRITE(ILUOUT0,*) "NLIGNE hail",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
WRITE(ILUOUT0,*) "ZTC_T_H(ILINE),ZELEV_T_H(ILINE)", ZTC_T_H(ILINE),ZELEV_T_H(ILINE)
WRITE(ILUOUT0,*) "ZM_T_H(ILINE),ZS11_CARRE_T_H(ILINE),ZS22_CARRE_T_H(ILINE),ZRE_S22S11_T_H(ILINE)",&
ZM_T_W(ILINE),ZS11_CARRE_T_H(ILINE),ZS22_CARRE_T_H(ILINE),ZRE_S22S11_T_H(ILINE)
WRITE(ILUOUT0,*) "ZIM_S22S11_T_H(ILINE),ZRE_S22FMS11FT_T_H(ILINE),ZIM_S22FT_T_H(ILINE),ZIM_S11FT_T_H(ILINE)",&
ZIM_S22S11_T_H(ILINE),ZRE_S22FMS11FT_T_H(ILINE),ZIM_S22FT_T_H(ILINE),ZIM_S11FT_T_H(ILINE)
ENDIF
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
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.
ZAEHINT=1.
ZAVHINT=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)))
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