Skip to content
Snippets Groups Projects
Commit 43fb460f authored by Gaelle DELAUTIER's avatar Gaelle DELAUTIER
Browse files

Gaelle 4/5/18 : add radar simulator for ICE4

parent e7724a7b
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/MNH/radar_scattering.f90
+ 419
− 71
View file @ 43fb460f
!MNH_LIC Copyright 2004-2018 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC Copyright 1994-2014 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.
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$ $Date$
!-----------------------------------------------------------------
! ######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)
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
......@@ -27,14 +31,15 @@ REAL,DIMENSION(:,:,:,:,:), INTENT(INOUT) :: PZE ! 5D matrix (iradar, ielev, i
!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
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)
PS_RAY,PG_RAY,PVDOP_RAY,PELEV,PX_H,PX_V,PW_H,PW_V,PZE,PBU_MASK_RAY,PCR_RAY,PH_RAY)
! ##############################
!
!!**** *RADAR_SCATTERING* - computes radar reflectivities.
......@@ -104,6 +109,17 @@ USE MODD_CST
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_LUNIT
USE MODD_PARAMETERS
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,&
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
!!LIMA
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,&
......@@ -116,26 +132,20 @@ USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XN
!!LIMA
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
USE MODE_FSCATTER
USE MODE_READTMAT
USE MODE_FGAU , ONLY:GAULAG
USE MODI_GAMMA, ONLY:GAMMA
!
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll,IO_FILE_OPEN_ll
USE MODE_FSCATTER
USE MODE_IO_ll
USE MODD_LUNIT
USE MODE_IO_MANAGE_STRUCT, ONLY: IO_FILE_ADD2LIST
USE MODE_MSG
USE MODE_READTMAT
!
USE MODI_GAMMA, ONLY:GAMMA
!
IMPLICIT NONE
!
......@@ -159,7 +169,8 @@ REAL,DIMENSION(:,:,:,:,:), INTENT(INOUT) :: PZE ! gate equivalent reflectivity
! 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
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 :
!
......@@ -173,8 +184,8 @@ REAL, DIMENSION(:,:,:,:,:,:,:),ALLOCATABLE :: ZREFL
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 :: 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
!
......@@ -199,6 +210,7 @@ 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
......@@ -209,30 +221,35 @@ 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 :: ZCR,ZCI,ZCS,ZCG,ZCH ! 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
INTEGER :: IMAX
LOGICAL :: LPART_MASK ! indicates a partial mask along the beam
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
!
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,ZZHH,ZZE_R,ZZE_I,ZZE_S,ZZE_G
LOGICAL :: GTHRESHOLD_V, GTHRESHOLD_Z,GTHRESHOLD_ZR,GTHRESHOLD_ZI,GTHRESHOLD_ZS,GTHRESHOLD_ZG
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(4) :: YTAB_TYPE
CHARACTER(LEN=25),DIMENSION(4) :: YFILE_COEFINT
CHARACTER(LEN=1),DIMENSION(5) :: YTAB_TYPE
CHARACTER(LEN=25),DIMENSION(5) :: YFILE_COEFINT
REAL,DIMENSION(4) :: ZELEV_MIN,ZELEV_MAX,ZELEV_STEP,&
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(4) :: INB_ELEV,INB_TC,INB_FW,INB_LINE
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
......@@ -245,7 +262,8 @@ REAL, DIMENSION(:),ALLOCATABLE :: ZIM_S22S11_T_R, ZIM_S22S11_T_S, ZIM_S22S11_T_G
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
......@@ -257,9 +275,9 @@ INTEGER :: ILUOUT0,IUNIT
!
! MODIF GAELLE POUR LIMA
!
LOGICAL :: GLIMA
REAL,DIMENSION(4) :: ZCC_MIN,ZCC_MAX, ZCC_STEP
INTEGER,DIMENSION(4):: INB_CC
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
......@@ -270,6 +288,7 @@ 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 :: ZCCH,ZLBH,ZLBEXH,ZDH,ZALPHAH,ZNUH,ZAH,ZBH,ZCXH
REAL :: ZLBI,ZLBEXI,ZDI,ZALPHAI,ZNUI,ZAI,ZBI
REAL,DIMENSION(:),ALLOCATABLE :: ZRTMIN
CHARACTER(LEN=100) :: YMSG
......@@ -285,6 +304,11 @@ IF (PRESENT(PCR_RAY)) THEN
ELSE
GLIMA=.FALSE.
ENDIF
IF (PRESENT(PH_RAY)) THEN
GHAIL=.TRUE.
ELSE
GHAIL=.FALSE.
ENDIF
!
!
!
......@@ -304,6 +328,10 @@ ENDIF
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
......@@ -374,26 +402,97 @@ ELSE ! ICE3
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
ENDIF
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"
!
! 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
ELSE
IRHOHV=10
IVDOP=IZEG+1 !VRU
END IF
IPHIDP=IRHOHV+1
IDELTAHV=IPHIDP+1
IRHR=IDELTAHV+1
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
IZDA=IRHG+1
IF (GHAIL) THEN
IRHH=IRHG+1
IZDA=IRHH+1
ELSE
IZDA=IRHG+1
ENDIF
IZDS=IZDA+1
IZDG=IZDS+1
IKDR=IZDG+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)
......@@ -413,12 +512,14 @@ IF(LFALL) THEN
ZCI=XC_I_I
ZCS=XCS_I
ZCG=XCG_I
IF (GHAIL) ZCH=XCH_I
ENDIF
ELSE
ZCR=0.
ZCI=0.
ZCS=0.
ZCG=0.
IF (GHAIL) ZCH=0.
END IF
! Calculation of nodes and weights for the Gauss-Laguerre quadrature
......@@ -428,7 +529,7 @@ IF(NDIFF/=0) THEN
CALL GAULAG(NPTS_GAULAG,ZX,ZW)
END IF
!
IVDOP=9 !index of Doppler Velocity (VRU) in ZREFL
!
IMAX=SIZE(PZE,5)
WRITE(ILUOUT0,*) "-----------------"
WRITE(ILUOUT0,*) "Radar scattering"
......@@ -441,6 +542,7 @@ 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
......@@ -468,6 +570,11 @@ 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
......@@ -493,7 +600,7 @@ IF(NDIFF==7) THEN
YTAB_TYPE(2)='s'
YTAB_TYPE(3)='g'
YTAB_TYPE(4)='w'
!
YTAB_TYPE(5)='h'
! definition des paramètres de lecture de la table T-matrice
! all mixing ratio
ZEXPM_MIN=-7.
......@@ -545,8 +652,20 @@ IF(NDIFF==7) THEN
ZCC_MIN(4)=1.
ZCC_STEP(4)=1.
ZCC_MAX(4)=1.
DO JT=1,4
! 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
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
......@@ -579,13 +698,13 @@ DO JI=1,INBRAD
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,*) "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','')
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SCATTERING','')
ENDIF
!************ fichiers Min Max Pas et Coef Tmat ***********
DO JT=1,4 !types (r, s, g, w)
DO JT=1,5 !types (r, s, g, w, h)
YTYPE=YTAB_TYPE(JT)
IF (JT .EQ. 1) THEN
......@@ -619,7 +738,12 @@ DO JI=1,INBRAD
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,&
......@@ -641,7 +765,7 @@ DO JI=1,INBRAD
ILINE=ILINE+1
ENDDO
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
TZFILE => NULL()
WRITE(ILUOUT0,*) "NLIGNE rain",ILINE
ILINE=2
WRITE(ILUOUT0,*) "ILINE=",ILINE
......@@ -732,6 +856,33 @@ DO JI=1,INBRAD
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
CALL IO_FILE_ADD2LIST(TZFILE,YFILE_COEFINT(5),'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(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
CALL IO_FILE_CLOSE_ll(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
ENDIF !END IF T-MATRIX => END OF LOOKUP TABLE READING
!============== loop over elevations =================
......@@ -760,6 +911,8 @@ DO JI=1,INBRAD
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.
......@@ -1465,6 +1618,178 @@ DO JI=1,INBRAD
END IF ! **************** end GRAUPEL (end IF SIZE(PG_RAY,1) > 0)
!-----------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
!**********************************
!**********************************
!**********************************
!**********************************
!---------------------------------------------------------------------------------------------------
!* 6. HAIL
! -------
!
!
IF (GHAIL) THEN
ZM=PRHODREF_RAY(JI,JEL,JAZ,JL,JH,JV)*PH_RAY(JI,JEL,JAZ,JL,JH,JV) !graupel content
IF(ZM > ZM_MIN) THEN
YTYPE='h'
ZQMI=SQRT(QEPSI(MIN(PT_RAY(JI,JEL,JAZ,JL,JH,JV),XTT),XLIGHTSPEED/XLAM_RAD(JI)))
ZQMW=SQRT(QEPSW(MAX(PT_RAY(JI,JEL,JAZ,JL,JH,JV),XTT),XLIGHTSPEED/XLAM_RAD(JI)))
!ini_radar.f90 : ZCXG = -0.5 XBG = 2.8 ( Xj et bj tab 2.1 p 24)
!ini_rain_ice.f90 : XLBEXG = 1.0/(XCXG-XBG) XAG = 19.6 (aj tab 2.1 p 24)
!XLBG = ( XAG*XCCG*MOMG(XALPHAG,XNUG,XBG) )**(-XLBEXG) (eq 2.6 p 23)
ZFW=0 !????????
ZLBDA=ZLBH*(PRHODREF_RAY(JI,JEL,JAZ,JL,JH,JV)*PH_RAY(JI,JEL,JAZ,JL,JH,JV))**ZLBEXH
!XTT : température du point triple de l'eau (273.16 K <=> 0.1 °C)
IF(PT_RAY(JI,JEL,JAZ,JL,JH,JV) > XTT) THEN ! mixture of ice and water
ZFRAC_ICE = .85 !(see p 68)
ELSE ! only ice
ZFRAC_ICE=1.
END IF
! from eq 3.77 p 68
!XRHOLW=1000 (initialized in ini_cst.f90)
ZDMELT_FACT=6.*ZAG/(XPI*XRHOLW*((1.-ZFRAC_ICE)+ZFRAC_ICE*0.92))
ZEXP=2.*ZBH
!Calculation of the refractive index from Bohren and Battan (3.72 p66)
ZQB=2.*ZQMW**2*(2.*ZQMI**2*LOG(ZQMI/ZQMW)/(ZQMI**2-ZQMW**2)-1.)/(ZQMI**2-ZQMW**2) !Beta (3.73 p66)
ZQM=SQRT(((1.-ZFRAC_ICE)*ZQMW**2+ZFRAC_ICE*ZQB*ZQMI**2)/(1.-ZFRAC_ICE+ZFRAC_ICE*ZQB)) ! Bohren & Battan (1982) 3.72 p66
ZQK=(ZQM**2-1.)/(ZQM**2+2.)
!Rayleigh, Rayleigh for ellipsoides or Rayleigh 6th order
IF(NDIFF==0.OR.NDIFF==3.OR.NDIFF==4) THEN
ZREFLOC(1:2)=ABS(ZQK)**2/.93*ZDMELT_FACT**2*1.E18*ZCCH*ZLBDA**(ZCXH-ZEXP)*MOMG(ZALPHAH,ZNUH,ZEXP)
ZREFLOC(3)=0.
IF(LWREFL) THEN ! weighting by reflectivities
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
-ZCH*SIN(PELEV(JI,JEL,JL,JV))*ABS(ZQK)**2/.93*ZDMELT_FACT**2&
*1.E18*ZCCH*ZLBDA**(ZCXH-ZEXP-ZDH)*MOMG(ZALPHAH,ZNUH,ZEXP+ZDH)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)+ZCCH*ZLBDA**ZCXH
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
-ZCH*SIN(PELEV(JI,JEL,JL,JV))&
*ZCCH*ZLBDA**(ZCXH-ZDH)*MOMG(ZALPHAH,ZNUH,ZDH)
END IF !end IF(LWREFL)
IF(LATT) THEN
IF(NDIFF==0.OR.NDIFF==3) THEN
ZAETMP(:)=ZCCH*ZLBDA**ZCXH*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK) &
*MOMG(ZALPHAH,ZNUH,ZBH)/ZLBDA**ZBH)
ELSE
ZAETMP(:)=ZCCH*ZLBDA**ZCXH*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK) &
*MOMG(ZALPHAH,ZNUH,ZBH)/ZLBDA**ZBH&
+ZDMELT_FACT**(5./3.)*XPI**4/15./XLAM_RAD(JI)**3 &
*AIMAG(ZQK**2*(ZQM**4+27.*ZQM**2+38.) &
/(2.*ZQM**2+3.))*MOMG(ZALPHAH,ZNUH,5.*ZBH/3.)/ZLBDA**(5.*ZBH/3.)&
+ZDMELT_FACT**2 *2.*XPI**5/3. /XLAM_RAD(JI)**4*REAL(ZQK**2) &
*MOMG(ZALPHAH,ZNUH,2.*ZBH) /ZLBDA**(2.*ZBH))
END IF ! end IF(NDIFF==0.OR.NDIFF==3)
END IF ! end IF(LATT)
ZRE_S22S11_H=0
ZIM_S22S11_H=0
ZS22_CARRE_H=0
ZS11_CARRE_H=0
!******************************* NDIFF==7 TmatInt ************************************
ELSE IF(NDIFF==7) THEN
ZREFLOC(:)=0
IF(LATT) ZAETMP(:)=0
CALL CALC_KTMAT(PELEV(JI,JEL,JL,JV), PT_RAY(JI,JEL,JAZ,JL,JH,JV),&
ZFW,ZM,&
ZELEV_MIN(3),ZELEV_MAX(3),ZELEV_STEP(3),&
ZTC_MIN(3),ZTC_MAX(3),ZTC_STEP(3),&
ZFW_MIN(3),ZFW_MAX(3),ZFW_STEP(3),&
ZEXPM_MIN,ZEXPM_MAX,ZEXPM_STEP,&
ITMAT,ZELEV_RED,ZTC_RED,ZFW_RED,ZM_RED)
IF (ITMAT(1) .NE. -NUNDEF) THEN
DO JIND=1,SIZE(KMAT_COEF,2),1
KMAT_COEF(1,JIND)=ZS11_CARRE_T_H(ITMAT(JIND))
KMAT_COEF(2,JIND)=ZS22_CARRE_T_H(ITMAT(JIND))
KMAT_COEF(3,JIND)=ZRE_S22S11_T_H(ITMAT(JIND))
KMAT_COEF(4,JIND)=ZIM_S22S11_T_H(ITMAT(JIND))
KMAT_COEF(5,JIND)=ZRE_S22FMS11FT_T_H(ITMAT(JIND))
KMAT_COEF(6,JIND)=ZIM_S22FT_T_H(ITMAT(JIND))
KMAT_COEF(7,JIND)=ZIM_S11FT_T_H(ITMAT(JIND))
ENDDO
CALL INTERPOL(ZELEV_RED,ZTC_RED,ZFW_RED,ZM_RED,KMAT_COEF,ZS11_CARRE_H,ZS22_CARRE_H,&
ZRE_S22S11_H,ZIM_S22S11_H,ZRE_S22FMS11F,ZIM_S22FT,ZIM_S11FT)
ELSE
ZS11_CARRE_H=0
ZS22_CARRE_H=0
ZRE_S22S11_H=0
ZIM_S22S11_H=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_H
ZREFLOC(2)=1.E18*(XLAM_RAD(JI))**4/(XPI**5*.93)*4*XPI*ZS11_CARRE_H
ZREFLOC(3)=180.E3/XPI*XLAM_RAD(JI)*ZRE_S22FMS11F
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
-ZCH*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(1) &
*1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCCH/4./ZLBDA**(3+ZDH)
IF(LATT) THEN
ZAETMP(1)=ZIM_S22FT*XLAM_RAD(JI)*2
ZAETMP(2)=ZIM_S11FT*XLAM_RAD(JI)*2
END IF
ELSE ! Mie (NDIFF=1)
ZREFLOC(:)=0.
IF(LATT) ZAETMP(:)=0.
DO JJ=1,NPTS_GAULAG ! ****** Gauss-Laguerre quadrature
ZD=ZX(JJ)**(1./ZALPHAH)/ZLBDA
ZDE=ZDMELT_FACT**(1./3.)*ZD**(ZBH/3.)
CALL BHMIE(XPI/XLAM_RAD(JI)*ZDE,ZQM,ZQEXT(1),ZQSCA,ZQBACK(1))
ZQBACK(2)=ZQBACK(1)
ZQEXT(2)=ZQEXT(1) ! modif Clotilde 23/04/2012
ZQBACK(3)=0.
ZREFLOC(1:3)=ZREFLOC(1:3)+ZQBACK(1:3)*ZX(JJ)**(ZNUH-1.+2.*ZBH/3./ZALPHAH)*ZW(JJ)
ZREFLOC(4)=ZREFLOC(4)+ZQBACK(1)*ZX(JJ)**(ZNUH-1.+2.*ZBH/3./ZALPHAH+ZDH/ZALPHAH)*ZW(JJ)
IF(LATT) ZAETMP(:)=ZAETMP(:)+ZQEXT(:)*ZX(JJ)**(ZNUH-1.+2.*ZBH/3./ZALPHAH)*ZW(JJ)
END DO ! ****** end loop on diameter (Gauss-Laguerre)
ZREFLOC(1:2)=ZREFLOC(1:2)*1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCH &
*ZLBDA**(ZCXH-2.*ZBH/3.)/(4.*GAMMA(ZNUH)*.93)*ZDMELT_FACT**(2./3.)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP) &
+PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
-ZCH*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(4) &
*1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCH &
*ZLBDA**(ZCXH-2.*ZBH/3.-ZDH)/(4.*GAMMA(ZNUH)*.93)*ZDMELT_FACT**(2./3.)
IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZCCH*ZLBDA**(ZCXH-2.*ZBH/3.)/(4.*GAMMA(ZNUH)) &
*ZDMELT_FACT**(2./3.)
ZRE_S22S11_H=0
ZIM_S22S11_H=0
ZS22_CARRE_H=0
ZS11_CARRE_H=0 !0 in case of Mie
END IF !**************** end loop for each type of diffusion : IF(NDIFF==0.OR.NDIFF==3.OR.NDIFF==4)
ZREFL(JI,JEL,JAZ,JL,JH,JV,1:3)=ZREFL(JI,JEL,JAZ,JL,JH,JV,1:3)+ZREFLOC(1:3)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZEH)=ZREFLOC(1) ! z_e due to graupel
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDH)=ZREFLOC(2) !Zvv for ZDR due to graupel
ZREFL(JI,JEL,JAZ,JL,JH,JV,IKDH)=ZREFLOC(3) !Zvv for ZDR due to graupel
IF (ZS22_CARRE_H*ZS11_CARRE_H .GT. 0) THEN
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHH)=SQRT(ZRE_S22S11_H**2+ZIM_S22S11_H**2)/SQRT(ZS22_CARRE_H*ZS11_CARRE_H)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHH)=1
END IF
IF(LATT) THEN
ZAELOC(JI,JEL,JAZ,JL,JH,JV,:)=ZAELOC(JI,JEL,JAZ,JL,JH,JV,:)+ZAETMP(:)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IAEH)=ZAETMP(1)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IAVH)=ZAETMP(2)
IF(JL>1) THEN
ZAEHINT=ZAEHINT*EXP(-2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IAEH)*XSTEP_RAD)
ZAVHINT=ZAVHINT*EXP(-2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IAVH)*XSTEP_RAD)
END IF
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZEH)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IZEH)*ZAEHINT ! Z_g attenuated
ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDH)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IZDH)*ZAVHINT ! Z_g attenuated
END IF !end IF(LATT)
END IF !**************** IF(PH_RAY(JI,JEL,JAZ,JL,JH,JV) > XRTMIN(6))
! Total attenuation even if no hydrometeors
IF(LATT.AND.JL>1) ZREFL(JI,JEL,JAZ,JL,JH,JV,IATH)=ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IATH) &
*EXP(-2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IAEH)*XSTEP_RAD)
END IF ! **************** end HAIL (end IF SIZE(PH_RAY,1) > 0)
!-----------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
!**********************************
!**********************************
!**********************************
!**********************************
IF(LWREFL) THEN ! weighting by reflectivities
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
......@@ -1476,26 +1801,26 @@ DO JI=1,INBRAD
ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)/ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)&
+PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)
END IF
!Calculation of Phidp (ZREFL(JI,JEL,JAZ,JL,JH,JV,IPHIDP) is initialized to 0 before the loop
IF (JL>1) ZREFL(JI,JEL,JAZ,JL,JH,JV,IPHIDP)=ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IPHIDP)+ &
!Calculation of Phidp (ZREFL(JI,JEL,JAZ,JL,JH,JV,IPDP) is initialized to 0 before the loop
IF (JL>1) ZREFL(JI,JEL,JAZ,JL,JH,JV,IPDP)=ZREFL(JI,JEL,JAZ,JL-1,JH,JV,IPDP)+ &
2.*ZREFL(JI,JEL,JAZ,JL-1,JH,JV,3)*XSTEP_RAD*1D-3
!Calculation of RhoHV and DeltaHV
ZRE_S22S11_T=ZRE_S22S11_R+ZRE_S22S11_I+ZRE_S22S11_S+ZRE_S22S11_G
ZIM_S22S11_T=ZIM_S22S11_R+ZIM_S22S11_I+ZIM_S22S11_S+ZIM_S22S11_G
ZS22_CARRE_T=ZS22_CARRE_R+ZS22_CARRE_I+ZS22_CARRE_S+ZS22_CARRE_G
ZS11_CARRE_T=ZS11_CARRE_R+ZS11_CARRE_I+ZS11_CARRE_S+ZS11_CARRE_G
ZRE_S22S11_T=ZRE_S22S11_R+ZRE_S22S11_I+ZRE_S22S11_S+ZRE_S22S11_G+ZRE_S22S11_H
ZIM_S22S11_T=ZIM_S22S11_R+ZIM_S22S11_I+ZIM_S22S11_S+ZIM_S22S11_G+ZIM_S22S11_H
ZS22_CARRE_T=ZS22_CARRE_R+ZS22_CARRE_I+ZS22_CARRE_S+ZS22_CARRE_G+ZS22_CARRE_H
ZS11_CARRE_T=ZS11_CARRE_R+ZS11_CARRE_I+ZS11_CARRE_S+ZS11_CARRE_G+ZS11_CARRE_H
!RhoHV
IF ((ZS22_CARRE_T*ZS11_CARRE_T)>0.) THEN
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHOHV)=SQRT(ZRE_S22S11_T**2+ZIM_S22S11_T**2)/SQRT(ZS22_CARRE_T*ZS11_CARRE_T)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHV)=SQRT(ZRE_S22S11_T**2+ZIM_S22S11_T**2)/SQRT(ZS22_CARRE_T*ZS11_CARRE_T)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHOHV)=-XUNDEF
ZREFL(JI,JEL,JAZ,JL,JH,JV,IRHV)=-XUNDEF
END IF
!DeltaHV
IF (ZRE_S22S11_T/=0) THEN
ZREFL(JI,JEL,JAZ,JL,JH,JV,IDELTAHV)=180/XPI*ATAN(ZIM_S22S11_T/ZRE_S22S11_T)
ZREFL(JI,JEL,JAZ,JL,JH,JV,IDHV)=180/XPI*ATAN(ZIM_S22S11_T/ZRE_S22S11_T)
ELSE
ZREFL(JI,JEL,JAZ,JL,JH,JV,IDELTAHV)=0
ZREFL(JI,JEL,JAZ,JL,JH,JV,IDHV)=0
END IF
ELSE !if temperature is not defined
ZREFL(JI,JEL,JAZ,JL,JH,JV,1:2)=XVALGROUND
......@@ -1521,11 +1846,16 @@ DO JI=1,INBRAD
!lookup tables for wet graupel
DEALLOCATE (ZTC_T_W,ZELEV_T_W,ZM_T_W,ZS11_CARRE_T_W,ZS22_CARRE_T_W,&
ZRE_S22S11_T_W,ZIM_S22S11_T_W,ZRE_S22FMS11FT_T_W,ZIM_S22FT_T_W,ZIM_S11FT_T_W)
IF (GHAIL) THEN
!lookup tables for hail
DEALLOCATE (ZTC_T_H,ZELEV_T_H,ZM_T_H,ZS11_CARRE_T_H,ZS22_CARRE_T_H,&
ZRE_S22S11_T_H,ZIM_S22S11_T_H,ZRE_S22FMS11FT_T_H,ZIM_S22FT_T_H,ZIM_S11FT_T_H)
ENDIF
END DO !JI
!
! attenuation in dB/km
IF(LATT) ZREFL(:,:,:,:,:,:,IAER:IAEG)=4343.*2.*ZREFL(:,:,:,:,:,:,IAER:IAEG) ! horizontal specific attenuation
IF(LATT) ZREFL(:,:,:,:,:,:,IAVR:IAVG)=4343.*2.*ZREFL(:,:,:,:,:,:,IAVR:IAVG) ! vertical specific attenuation
IF(LATT) ZREFL(:,:,:,:,:,:,IAER:IAEH)=4343.*2.*ZREFL(:,:,:,:,:,:,IAER:IAEH) ! horizontal specific attenuation
IF(LATT) ZREFL(:,:,:,:,:,:,IAVR:IAVH)=4343.*2.*ZREFL(:,:,:,:,:,:,IAVR:IAVH) ! vertical specific attenuation
! convective/stratiform
ZREFL(:,:,:,:,:,:,4)=PBU_MASK_RAY(:,:,:,:,:,:) ! CSR
! /convective/stratiform
......@@ -1546,7 +1876,7 @@ DO JI=1,INBRAD
DO JEL=1,IEL
DO JAZ=1,INBAZIM
IF (LATT) ZAETOT(:,:,1:2)=1.
PZE(JI,JEL,JAZ,1,IPHIDP)=0
PZE(JI,JEL,JAZ,1,IPDP)=0
DO JL=1,INBSTEPMAX
! if no undef point in gate JL and at least one point where T is defined
IF(COUNT(ZREFL(JI,JEL,JAZ,JL,:,:,1)==-XUNDEF)==0.AND. &
......@@ -1604,11 +1934,13 @@ DO JI=1,INBRAD
ZSNR_I=-XUNDEF
ZSNR_S=-XUNDEF
ZSNR_G=-XUNDEF
ZSNR_H=-XUNDEF
ZZHH=PZE(JI,JEL,JAZ,JL,1)
ZZE_R=PZE(JI,JEL,JAZ,JL,IZER)
ZZE_I=PZE(JI,JEL,JAZ,JL,IZEI)
ZZE_S=PZE(JI,JEL,JAZ,JL,IZES)
ZZE_G=PZE(JI,JEL,JAZ,JL,IZEG)
IF (GHAIL) ZZE_H=PZE(JI,JEL,JAZ,JL,IZEH)
ZDISTRAD=JL*XSTEP_RAD !radar distance in meters
IF (LSNRT) THEN
IF (ZZHH/=XVALGROUND .AND. ZZHH/=-XUNDEF.AND.ZZHH/=0) THEN
......@@ -1626,12 +1958,18 @@ DO JI=1,INBRAD
IF (ZZE_G/=XVALGROUND .AND. ZZE_G/=-XUNDEF.AND.ZZE_G/=0) THEN
ZSNR_G=10*LOG10(ZZE_G)-20*LOG10(ZDISTRAD/(100*10**3))
END IF
IF (GHAIL) THEN
IF (ZZE_H/=XVALGROUND .AND. ZZE_H/=-XUNDEF.AND.ZZE_H/=0) THEN
ZSNR_H=10*LOG10(ZZE_H)-20*LOG10(ZDISTRAD/(100*10**3))
END IF
END IF
GTHRESHOLD_V=(ZSNR>=XSNRMIN)
GTHRESHOLD_Z=GTHRESHOLD_V
GTHRESHOLD_ZR=(ZSNR_R>=XSNRMIN)
GTHRESHOLD_ZI=(ZSNR_I>=XSNRMIN)
GTHRESHOLD_ZS=(ZSNR_S>=XSNRMIN)
GTHRESHOLD_ZG=(ZSNR_G>=XSNRMIN)
IF (GHAIL) GTHRESHOLD_ZH=(ZSNR_H>=XSNRMIN)
ELSE
GTHRESHOLD_V=(ZZHH>=10**(XREFLVDOPMIN/10.))
GTHRESHOLD_Z=(ZZHH>=10**(XREFLMIN/10.))
......@@ -1639,6 +1977,7 @@ DO JI=1,INBRAD
GTHRESHOLD_ZI=(ZZE_I>=10**(XREFLMIN/10.))
GTHRESHOLD_ZS=(ZZE_S>=10**(XREFLMIN/10.))
GTHRESHOLD_ZG=(ZZE_G>=10**(XREFLMIN/10.))
IF (GHAIL) GTHRESHOLD_ZH=(ZZE_H>=10**(XREFLMIN/10.))
END IF
!--- Doppler velocities
IF(GTHRESHOLD_V) THEN
......@@ -1661,7 +2000,7 @@ DO JI=1,INBRAD
!--- Zhh, Zvv et variables globales
IF(GTHRESHOLD_Z .EQV. .FALSE.) THEN
PZE(JI,JEL,JAZ,JL,1:4)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IRHOHV:IDELTAHV)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IRHV:IDHV)=-XUNDEF
END IF
!--- ZER, ZDA, KDR, RHR
IF(GTHRESHOLD_ZR .EQV. .FALSE.) THEN
......@@ -1685,6 +2024,15 @@ DO JI=1,INBRAD
PZE(JI,JEL,JAZ,JL,IKDG)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IRHG)=-XUNDEF
END IF
!--- ZEH, ZDH, KDH, RHH
IF (GHAIL) THEN
IF(GTHRESHOLD_ZH .EQV. .FALSE.) THEN
PZE(JI,JEL,JAZ,JL,IZEH)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IZDH)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IKDH)=-XUNDEF
PZE(JI,JEL,JAZ,JL,IRHH)=-XUNDEF
END IF
END IF
!--- ZEI
IF(GTHRESHOLD_ZI .EQV. .FALSE.) THEN
PZE(JI,JEL,JAZ,JL,IZEI)=-XUNDEF
......
src/MNH/radar_simulator.f90
+ 296
− 50
View file @ 43fb460f
!MNH_LIC Copyright 2004-2018 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC Copyright 1994-2014 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.
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/radar_simulator.f90,v $ $Revision: 1.1.2.3.2.1.12.2.2.2 $ $Date: 2015/09/16 14:31:20 $
!-----------------------------------------------------------------
! ###########################
MODULE MODI_RADAR_SIMULATOR
! ###########################
......@@ -214,6 +218,7 @@ REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZR_RAY ! rain mi
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZI_RAY ! pristine ice mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZS_RAY ! snow/aggregates mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZG_RAY ! graupel mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZH_RAY ! hail mixing ratio interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZCIT_RAY ! pristine ice concentration interpolated along the rays
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZRHODREF_RAY
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE :: ZVDOP_RAY
......@@ -239,15 +244,20 @@ REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZBU_MASK_RAY
REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE :: ZN_RAY,ZDNDZ_RAY ! refractivity and its vertical gradient in radar coordinates
REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZN,ZDNDZ ! index of ZN_RAY,ZDNDZ_RAY in ZZE
INTEGER,PARAMETER :: IZER=5,IZEI=6,IZES=7,IZEG=8
INTEGER,PARAMETER :: IVDOP=9
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
INTEGER :: IZER,IZEI,IZES,IZEG
INTEGER :: IVDOP
INTEGER :: IAER,IAEI,IAES,IAEG
INTEGER :: IAVR,IAVI,IAVS,IAVG
INTEGER :: IATR,IATI,IATS,IATG
INTEGER :: IRHV,IPDP,IDHV
INTEGER :: IRHR, IRHS, IRHG, IZDA, IZDS, IZDG, IKDR, IKDS, IKDG
INTEGER :: IHAS,IRFR,IDNZ
REAL :: ZZSTEP
INTEGER :: IZEH, IRHH,IKDH,IZDH ! hail
INTEGER :: IAEH,IAVH,IATH
INTEGER :: IMR,IMI,IMS,IMG,IMH,ICIT,ITEM,ICRT
LOGICAL :: GHAIL
INTEGER,DIMENSION(:,:,:,:),ALLOCATABLE ::IREFL_CART_NBH
!
!Modif pour LIMA
!
......@@ -291,6 +301,12 @@ IF (PRESENT(PCRT)) THEN
ELSE
GLIMA=.FALSE.
ENDIF
IF (SIZE(PRT,4)== 7 ) THEN
GHAIL=.TRUE.
ELSE
GHAIL=.FALSE.
ENDIF
!
!* 1.2 Some constants and parameters
!
......@@ -372,47 +388,110 @@ DO JH=1,NPTS_H
END DO
ZELEV(:,:,:,:)=ZELEV(:,:,:,:)*ZRDSDG ! in radian
! initialisation of refractivity indices
IRFR=1 ! this is used down there in the interpolation part
IDNZ=1 ! this is used down there in the interpolation part
!IHAS=10 !number of calculated radar variables (on the radar projection)
!IHAS=13 !add of RhoHV, PhiDP, DeltaHV in ZZE
IHAS=22 !add of RHR-RHG, ZDA-ZDG, KDR-KDG
IF(LREFR) THEN
IRFR=IHAS+7 ! add of TEM:
!"HAS","M_R","M_I","M_S","M_G","CIT","TEM"
! puis "RFR"
IF (GLIMA) IRFR=IRFR+1 ! add "CRT"
ENDIF
IF(LDNDZ) THEN
IF(LREFR) THEN
IDNZ=IHAS+8 !+7 !17 ! refractivity vertical gradient
ELSE
IDNZ=IHAS+7 !+6 !16 ! refractivity vertical gradient
END IF
IF (GLIMA) IDNZ=IDNZ+1
END IF
IF(LATT) THEN
IRFR=IRFR+12 !add of AVR-AVG (vertical specific attenuation)
IDNZ=IDNZ+12
IHAS=IHAS+12
! 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
ELSE
IVDOP=IZEG+1 !VRU
END IF
IF (LATT) THEN
IRHV=22 !"ZHH","ZDR","KDP","CSR","ZER","ZEI","ZES","ZEG","VRU"
!"AER","AEI","AES","AEG","AVR","AVI","AVS","AVG","ATR","ATI","ATS","ATG"
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=10
END IF
IRHV=IVDOP+1
ENDIF
IPDP=IRHV+1
IDHV=IPDP+1
IRHR=IDHV+1
IRHS=IRHR+1
IRHG=IRHS+1
IZDA=IRHG+1
IF (GHAIL) THEN
IRHH=IRHG+1
IZDA=IRHH+1
ELSE
IZDA=IRHG+1
ENDIF
IZDS=IZDA+1
IZDG=IZDS+1
IKDR=IZDG+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
IHAS=IKDH+1
ELSE
IHAS=IKDG+1
ENDIF
IMR=IHAS+1
IMI=IMR+1
IMS=IMI+1
IMG=IMS+1
IF (GHAIL) THEN
IMH=IMG+1
ICIT=IMH+1
ELSE
ICIT=IMG+1
ENDIF
ITEM=ICIT+1
IF (GLIMA) THEN
ICRT=ITEM+1
IF(LREFR) THEN
IRFR=ICRT+1
IF(LDNDZ) IDNZ=IRFR+1
ELSE
IF(LDNDZ) IDNZ=ICRT+1
ENDIF
ELSE
IF(LREFR) THEN
IRFR=ITEM+1
IF(LDNDZ) IDNZ=IRFR+1
ELSE
IF(LDNDZ) IDNZ=ITEM+1
ENDIF
ENDIF
!
!----------------------------------------------------------------------------------------
!* 2. RAY TRACING DEFINITION
......@@ -528,6 +607,8 @@ DO JI=1,NBRAD
IF(XRPK<0.) ZLAT(JH,JV)=-ZLAT(JH,JV) ! projection from north pole
IF(ABS(ZRPK-1.)>1.E-10 .AND. ABS(COS(ZRDSDG*ZLAT(JH,JV)))<1.E-10) THEN
WRITE(ILUOUT0,*) 'Error in projection : '
WRITE(ILUOUT0,*) 'pole in the domain, but not with stereopolar projection'
!callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','RADAR_SIMULATOR',&
'Error in projection: pole in the domain, but not with stereopolar projection')
......@@ -696,6 +777,13 @@ IF (SIZE(PRT,4)>5) THEN
ALLOCATE(ZG_RAY(NBRAD,IIELV,NBAZIM,NBSTEPMAX+1,NPTS_H,NPTS_V))
TVARRAD(INVAR)%P=>ZG_RAY(:,:,:,:,:,:)
END IF
! HAIL
IF (SIZE(PRT,4)>6) THEN
INVAR=INVAR+1
TVARMOD(INVAR)%P=>PRT(:,:,:,7)
ALLOCATE(ZH_RAY(NBRAD,IIELV,NBAZIM,NBSTEPMAX+1,NPTS_H,NPTS_V))
TVARRAD(INVAR)%P=>ZH_RAY(:,:,:,:,:,:)
END IF
! convective/stratiform
TVARMOD(INVAR+1)%P=>ZBU_MASK(:,:,:)
ALLOCATE(ZBU_MASK_RAY(NBRAD,IIELV,NBAZIM,NBSTEPMAX+1,NPTS_H,NPTS_V))
......@@ -767,11 +855,21 @@ ALLOCATE(ZZE(NBRAD,IIELV,NBAZIM,NBSTEPMAX+1,SIZE(PREFL_CART(:,:,:,:,:),5)))
!-----------------------------------------------------------------------
IF (GLIMA) THEN
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
ZELEV,ZX_H,ZX_V,ZW_H,ZW_V,ZZE(:,:,:,:,1:IHAS-1),ZBU_MASK_RAY,ZCRT_RAY)
IF (GHAIL) THEN
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
ZELEV,ZX_H,ZX_V,ZW_H,ZW_V,ZZE(:,:,:,:,1:IHAS-1),ZBU_MASK_RAY,PCR_RAY=ZCRT_RAY,PH_RAY=ZH_RAY)
ELSE
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
ZELEV,ZX_H,ZX_V,ZW_H,ZW_V,ZZE(:,:,:,:,1:IHAS-1),ZBU_MASK_RAY,PCR_RAY=ZCRT_RAY)
ENDIF
ELSE
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
IF (GHAIL) THEN
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
ZELEV,ZX_H,ZX_V,ZW_H,ZW_V,ZZE(:,:,:,:,1:IHAS-1),ZBU_MASK_RAY,PH_RAY=ZH_RAY)
ELSE
CALL RADAR_SCATTERING(ZT_RAY,ZRHODREF_RAY,ZR_RAY,ZI_RAY,ZCIT_RAY,ZS_RAY,ZG_RAY,ZVDOP_RAY, &
ZELEV,ZX_H,ZX_V,ZW_H,ZW_V,ZZE(:,:,:,:,1:IHAS-1),ZBU_MASK_RAY)
ENDIF
ENDIF
DEALLOCATE(ZVDOP_RAY)
! convective/stratiform
......@@ -792,7 +890,7 @@ END WHERE
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZR_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+1)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,IMR)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
! M_i
WHERE(ZRHODREF_RAY/=-XUNDEF .AND. ZI_RAY/=-XUNDEF)
ZWORK=ZRHODREF_RAY*ZI_RAY
......@@ -802,7 +900,7 @@ END WHERE
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZI_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+2)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,IMI)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
! M_s
WHERE(ZRHODREF_RAY/=-XUNDEF .AND. ZS_RAY/=-XUNDEF)
ZWORK=ZRHODREF_RAY*ZS_RAY
......@@ -812,7 +910,7 @@ END WHERE
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZS_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+3)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,IMS)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
! M_g
WHERE(ZRHODREF_RAY/=-XUNDEF .AND. ZG_RAY/=-XUNDEF)
ZWORK=ZRHODREF_RAY*ZG_RAY
......@@ -822,8 +920,24 @@ END WHERE
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZG_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+4)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,IMG)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
DEALLOCATE(ZWORK)
IF (GHAIL) THEN
ALLOCATE(ZWORK(SIZE(ZRHODREF_RAY,1),SIZE(ZRHODREF_RAY,2),&
SIZE(ZRHODREF_RAY,3),SIZE(ZRHODREF_RAY,4),&
SIZE(ZRHODREF_RAY,5),SIZE(ZRHODREF_RAY,6)))
! M_h
WHERE(ZRHODREF_RAY/=-XUNDEF .AND. ZH_RAY/=-XUNDEF)
ZWORK=ZRHODREF_RAY*ZH_RAY
ELSEWHERE
ZWORK=-XUNDEF
END WHERE
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZH_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IMH)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
DEALLOCATE(ZWORK)
ENDIF
! CIT
IF (GLIMA)THEN
ALLOCATE(ZWORK(SIZE(ZRHODREF_RAY,1),SIZE(ZRHODREF_RAY,2),&
......@@ -838,12 +952,14 @@ IF (GLIMA)THEN
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+5)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,ICIT)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
DEALLOCATE(ZWORK)
ELSE
ZZE(:,:,:,:,IHAS+5)=ZCIT_RAY(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,ICIT)=ZCIT_RAY(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ENDIF
ZZE(:,:,:,:,IHAS+6)=ZT_RAY(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2) ! temperature
! temperature
ZZE(:,:,:,:,ITEM)=ZT_RAY(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
!CRT
IF (GLIMA)THEN
ALLOCATE(ZWORK(SIZE(ZRHODREF_RAY,1),SIZE(ZRHODREF_RAY,2),&
SIZE(ZRHODREF_RAY,3),SIZE(ZRHODREF_RAY,4),&
......@@ -856,7 +972,7 @@ IF (GLIMA)THEN
WHERE(ZRHODREF_RAY==XVALGROUND .OR. ZCRT_RAY==XVALGROUND)
ZWORK=XVALGROUND
END WHERE
ZZE(:,:,:,:,IHAS+7)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
ZZE(:,:,:,:,ICRT)=ZWORK(:,:,:,:,(NPTS_H+1)/2,(NPTS_V+1)/2)
DEALLOCATE(ZWORK)
ENDIF
IF(LREFR) THEN
......@@ -888,12 +1004,14 @@ IF (LCART_RAD) THEN !if cartesian interpolation
ALLOCATE(IREFL_CART_NBI(NBRAD,IIELV,2*NMAX,2*NMAX))
ALLOCATE(IREFL_CART_NBS(NBRAD,IIELV,2*NMAX,2*NMAX))
ALLOCATE(IREFL_CART_NBG(NBRAD,IIELV,2*NMAX,2*NMAX))
IF (GHAIL) ALLOCATE(IREFL_CART_NBH(NBRAD,IIELV,2*NMAX,2*NMAX))
PREFL_CART(:,:,:,:,:)=0.
IREFL_CART_NB(:,:,:,:)=0
IREFL_CART_NBR(:,:,:,:)=0
IREFL_CART_NBS(:,:,:,:)=0
IREFL_CART_NBI(:,:,:,:)=0
IREFL_CART_NBG(:,:,:,:)=0
IF (GHAIL) IREFL_CART_NBH(:,:,:,:)=0
IVDOP_CART_NB(:,:,:,:)=0
!
!* 5.1 reflectivity on a cartesian grid (this is the way DSO/CMR creates BUFRs)
......@@ -903,8 +1021,8 @@ IF (LCART_RAD) THEN !if cartesian interpolation
DO JEL=1,IEL
DO JAZ=1,NBAZIM
DO JL=1,NBSTEPMAX+1
IF ((ZZE(JI,JEL,JAZ,JL,IHAS+6)/=XVALGROUND).AND.(ZZE(JI,JEL,JAZ,JL,IHAS+6)/=-XUNDEF)) THEN !conversion en °C
ZZE(JI,JEL,JAZ,JL,IHAS+6)=ZZE(JI,JEL,JAZ,JL,IHAS+6)-273.15
IF ((ZZE(JI,JEL,JAZ,JL,ITEM)/=XVALGROUND).AND.(ZZE(JI,JEL,JAZ,JL,ITEM)/=-XUNDEF)) THEN !conversion en °C
ZZE(JI,JEL,JAZ,JL,ITEM)=ZZE(JI,JEL,JAZ,JL,ITEM)-273.15
ENDIF
IXGRID=CEILING(NMAX+((JL-1)*XSTEP_RAD*SIN(ZAZIM_BASE(JAZ))/XGRID))
IYGRID=CEILING(NMAX+((JL-1)*XSTEP_RAD*COS(ZAZIM_BASE(JAZ))/XGRID))
......@@ -1196,7 +1314,63 @@ IF (LCART_RAD) THEN !if cartesian interpolation
END IF
IREFL_CART_NBG(JI,JEL,IXGRID,IYGRID)=IREFL_CART_NBG(JI,JEL,IXGRID,IYGRID)+1
END IF
!**********************************************
!**** HAIL VARIABLES ****
!**********************************************
IF (GHAIL) THEN
IF(ZZE(JI,JEL,JAZ,JL,IZEH)==XVALGROUND.OR.PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH)==XVALGROUND &
.OR.(LREFR.AND.ZZE(JI,JEL,JAZ,JL,IRFR)==XVALGROUND) & ! case for refractivity at boundaries
.OR.(LDNDZ.AND.ZZE(JI,JEL,JAZ,JL,IDNZ)==XVALGROUND) & ! case for refractivity gradient at origin
) THEN ! if any XVALGROUND in the pixel for ZES -> pixel set to XVALGROUND for all snow variables
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH)=XVALGROUND
PREFL_CART(JI,JEL,IXGRID,IYGRID,IRHH)=XVALGROUND
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZDH)=XVALGROUND
PREFL_CART(JI,JEL,IXGRID,IYGRID,IKDH)=XVALGROUND
IF (LATT) THEN
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAEH)=XVALGROUND
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAVH)=XVALGROUND
PREFL_CART(JI,JEL,IXGRID,IYGRID,IATH)=XVALGROUND
END IF
IREFL_CART_NBG(JI,JEL,IXGRID,IYGRID)=1
!-xundef for ZEG
ELSE IF(ZZE(JI,JEL,JAZ,JL,IZEH)==-XUNDEF.OR.PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH)==-XUNDEF &
.OR.(LREFR.AND.ZZE(JI,JEL,JAZ,JL,IRFR)==-XUNDEF) & ! case for refractivity at boundaries
.OR.(LDNDZ.AND.ZZE(JI,JEL,JAZ,JL,IDNZ)==-XUNDEF) & ! case for refractivity gradient at origin
) THEN ! if any -XUNDEF in the pixel for ZHH-> pixel set to -XUNDEF for all general variables
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH)=-XUNDEF
PREFL_CART(JI,JEL,IXGRID,IYGRID,IRHH)=-XUNDEF
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZDH)=-XUNDEF
PREFL_CART(JI,JEL,IXGRID,IYGRID,IKDH)=-XUNDEF
IF (LATT) THEN
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAEH)=-XUNDEF
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAVH)=-XUNDEF
PREFL_CART(JI,JEL,IXGRID,IYGRID,IATH)=-XUNDEF
END IF
IREFL_CART_NBG(JI,JEL,IXGRID,IYGRID)=1
!if no xvalground and no -xundef for REFL: incrementation of all polar pixels inside the cartesian pixel
ELSE
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IZEH) &
+ZZE(JI,JEL,JAZ,JL,IZEH)
PREFL_CART(JI,JEL,IXGRID,IYGRID,IRHH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IRHH) &
+ZZE(JI,JEL,JAZ,JL,IRHH)
PREFL_CART(JI,JEL,IXGRID,IYGRID,IZDH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IZDH) &
+ZZE(JI,JEL,JAZ,JL,IZDH)
PREFL_CART(JI,JEL,IXGRID,IYGRID,IKDH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IKDH) &
+ZZE(JI,JEL,JAZ,JL,IKDH)
IF (LATT) THEN
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAEH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IAEH) &
+ZZE(JI,JEL,JAZ,JL,IAEH)
PREFL_CART(JI,JEL,IXGRID,IYGRID,IAVH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IAVH) &
+ZZE(JI,JEL,JAZ,JL,IAVH)
PREFL_CART(JI,JEL,IXGRID,IYGRID,IATH)=PREFL_CART(JI,JEL,IXGRID,IYGRID,IATH) &
+ZZE(JI,JEL,JAZ,JL,IATH)
END IF
IREFL_CART_NBH(JI,JEL,IXGRID,IYGRID)=IREFL_CART_NBH(JI,JEL,IXGRID,IYGRID)+1
END IF
END IF
END DO !JL
END DO !JAZ
END DO !JEL
......@@ -1458,7 +1632,61 @@ DEALLOCATE(ZZE)
ENDIF
ENDIF
END IF !******** END GRAUPEL
!**** HAIL VARIABLES
IF (GHAIL) THEN
IF(IREFL_CART_NBH(JI,JEL,JH,JV) == 0) THEN ! no values inside the cartesian pixel
IF((JH+SIGN(.5,JH-.5-NMAX)-.5-NMAX)**2+(JV+SIGN(.5,JV-.5-NMAX)-.5-NMAX)**2>=NMAX**2) THEN
! out of range
PREFL_CART(JI,JEL,JH,JV,IZEH)=XVALGROUND
PREFL_CART(JI,JEL,JH,JV,IRHH)=XVALGROUND
PREFL_CART(JI,JEL,JH,JV,IZDH)=XVALGROUND
PREFL_CART(JI,JEL,JH,JV,IKDH)=XVALGROUND
IF (LATT) THEN
PREFL_CART(JI,JEL,JH,JV,IAEH)=XVALGROUND
PREFL_CART(JI,JEL,JH,JV,IAVH)=XVALGROUND
PREFL_CART(JI,JEL,JH,JV,IATH)=XVALGROUND
END IF
ELSE
PREFL_CART(JI,JEL,JH,JV,IZEH)=0
PREFL_CART(JI,JEL,JH,JV,IRHH)=0
PREFL_CART(JI,JEL,JH,JV,IZDH)=0
PREFL_CART(JI,JEL,JH,JV,IKDH)=0
IF (LATT) THEN
PREFL_CART(JI,JEL,JH,JV,IAEH)=0
PREFL_CART(JI,JEL,JH,JV,IAVH)=0
PREFL_CART(JI,JEL,JH,JV,IATH)=0
ENDIF
WRITE(ILUOUT0,*) "Warning: some pixels have no reflectivity; increase XGRID or decrease XSTEP_RAD"
END IF
ELSE
PREFL_CART(JI,JEL,JH,JV,IZEH)=PREFL_CART(JI,JEL,JH,JV,IZEH)/IREFL_CART_NBH(JI,JEL,JH,JV)
PREFL_CART(JI,JEL,JH,JV,IRHH)=PREFL_CART(JI,JEL,JH,JV,IRHH)/IREFL_CART_NBH(JI,JEL,JH,JV)
PREFL_CART(JI,JEL,JH,JV,IZDH)=PREFL_CART(JI,JEL,JH,JV,IZDH)/IREFL_CART_NBH(JI,JEL,JH,JV)
PREFL_CART(JI,JEL,JH,JV,IKDH)=PREFL_CART(JI,JEL,JH,JV,IKDH)/IREFL_CART_NBH(JI,JEL,JH,JV)
IF (LATT) THEN
PREFL_CART(JI,JEL,JH,JV,IAEH)=PREFL_CART(JI,JEL,JH,JV,IAEH)/IREFL_CART_NBH(JI,JEL,JH,JV)
PREFL_CART(JI,JEL,JH,JV,IAVH)=PREFL_CART(JI,JEL,JH,JV,IAVH)/IREFL_CART_NBH(JI,JEL,JH,JV)
PREFL_CART(JI,JEL,JH,JV,IATH)=PREFL_CART(JI,JEL,JH,JV,IATH)/IREFL_CART_NBH(JI,JEL,JH,JV)
END IF
! --------- Converting to dB -----------
IF(PREFL_CART(JI,JEL,JH,JV,IZEH)> 0) THEN
PREFL_CART(JI,JEL,JH,JV,IZEH)=10.*LOG10(PREFL_CART(JI,JEL,JH,JV,IZEH)) ! Z_equiv in dBZ
IF(PREFL_CART(JI,JEL,JH,JV,IZDH) > 0.) THEN
PREFL_CART(JI,JEL,JH,JV,IZDH)=PREFL_CART(JI,JEL,JH,JV,IZEH) &
-10.*LOG10(PREFL_CART(JI,JEL,JH,JV,IZDH)) ! Zdr=Z_HH-Z_VV
ENDIF
ELSE IF (PREFL_CART(JI,JEL,JH,JV,IZEH)== 0) THEN
PREFL_CART(JI,JEL,JH,JV,IZES)=-XUNDEF
END IF
IF(LATT) THEN
IF(PREFL_CART(JI,JEL,JH,JV,IATH)>0.) THEN
PREFL_CART(JI,JEL,JH,JV,IATH)=10.*LOG10(PREFL_CART(JI,JEL,JH,JV,IATH))
ENDIF
ENDIF
END IF !******** END HAIL
END IF
END DO
END DO
END DO
......@@ -1506,8 +1734,8 @@ ELSE ! if polar output
DO JAZ=1,NBAZIM
DO JL=1,NBSTEPMAX+1
!conversion en deg celcius
IF (PREFL_CART(JI,JEL,JAZ,JL,IHAS+6)/=-XUNDEF .AND. PREFL_CART(JI,JEL,JAZ,JL,IHAS+6)/=XVALGROUND) &
PREFL_CART(JI,JEL,JAZ,JL,IHAS+6)=PREFL_CART(JI,JEL,JAZ,JL,IHAS+6)-273.15
IF (PREFL_CART(JI,JEL,JAZ,JL,ITEM)/=-XUNDEF .AND. PREFL_CART(JI,JEL,JAZ,JL,ITEM)/=XVALGROUND) &
PREFL_CART(JI,JEL,JAZ,JL,ITEM)=PREFL_CART(JI,JEL,JAZ,JL,ITEM)-273.15
!------ ZHH and ZDR
IF(PREFL_CART(JI,JEL,JAZ,JL,1)> 0) THEN
......@@ -1583,6 +1811,24 @@ ELSE ! if polar output
PREFL_CART(JI,JEL,JAZ,JL,IATG)=10.*LOG10(PREFL_CART(JI,JEL,JAZ,JL,IATG))
END IF
END IF !------ END GRAUPEL
! --------- HAIL -----------
IF (GHAIL) THEN
IF(PREFL_CART(JI,JEL,JAZ,JL,IZEH)> 0) THEN
PREFL_CART(JI,JEL,JAZ,JL,IZEH)=10.*LOG10(PREFL_CART(JI,JEL,JAZ,JL,IZEH)) ! Z_equiv in dBZ
IF(PREFL_CART(JI,JEL,JAZ,JL,IZDH) > 0.) THEN
PREFL_CART(JI,JEL,JAZ,JL,IZDH)=PREFL_CART(JI,JEL,JAZ,JL,IZEH) &
-10.*LOG10(PREFL_CART(JI,JEL,JAZ,JL,IZDH)) ! Zdr=Z_HH-Z_VV
ENDIF
ELSE IF (PREFL_CART(JI,JEL,JAZ,JL,IZEH)== 0) THEN
PREFL_CART(JI,JEL,JAZ,JL,IZES)=-XUNDEF
END IF
IF(LATT) THEN
IF(PREFL_CART(JI,JEL,JAZ,JL,IATH)>0.) THEN
PREFL_CART(JI,JEL,JAZ,JL,IATH)=10.*LOG10(PREFL_CART(JI,JEL,JAZ,JL,IATH))
END IF
END IF !------ END HAIL
END IF
END DO
END DO
END DO
......
src/MNH/write_lfifm1_for_diag.f90
+ 48
− 9
View file @ 43fb460f
......@@ -3608,7 +3608,13 @@ IF(LRADAR .AND. LUSERR) THEN
IF (CCLOUD=='LIMA') INBOUT=INBOUT+1 ! rain concentration CRT
IF(LREFR) INBOUT=INBOUT+1 !+refractivity
IF(LDNDZ) INBOUT=INBOUT+1 !+refractivity vertical gradient
IF(LATT) INBOUT=INBOUT+12 !+AER-AEG AVR-AVG (vertical specific attenuation) and ATR-ATG
IF(LATT) INBOUT=INBOUT+12 !+AER-AEG AVR-AVG (vertical specific attenuation) and ATR-ATG
IF ( CCLOUD=='ICE4' ) THEN
INBOUT=INBOUT+5 ! HAIL ZEH RHH ZDH KDH M_H
IF (LATT) THEN
INBOUT=INBOUT+3 ! AEH AVH ATH
ENDIF
END IF
WRITE(ILUOUT0,*) "Nombre de variables dans ZWORK42 en sortie de radar_simulator:",INBOUT
IF (LCART_RAD) THEN
......@@ -3625,18 +3631,51 @@ IF(LRADAR .AND. LUSERR) THEN
CALL RADAR_SIMULATOR(XUT,XVT,XWT,XRT,XCIT,XRHODREF,ZTEMP,XPABSM,ZWORK42,ZWORK43)
ENDIF
ALLOCATE(YRAD(INBOUT))
YRAD(1:9)=(/"ZHH","ZDR","KDP","CSR","ZER","ZEI","ZES","ZEG","VRU"/)
ICURR=10
YRAD(1:8)=(/"ZHH","ZDR","KDP","CSR","ZER","ZEI","ZES","ZEG"/)
ICURR=9
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR)="ZEH"
ICURR=ICURR+1
END IF
YRAD(ICURR)="VRU"
ICURR=ICURR+1
IF(LATT) THEN
YRAD(ICURR:ICURR+11)=(/"AER","AEI","AES","AEG","AVR","AVI","AVS","AVG","ATR","ATI","ATS","ATG"/)
ICURR=ICURR+12
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR:ICURR+14)=(/"AER","AEI","AES","AEG","AEH","AVR","AVI","AVS","AVG","AVH","ATR","ATI","ATS","ATG","ATH"/)
ICURR=ICURR+15
ELSE
YRAD(ICURR:ICURR+11)=(/"AER","AEI","AES","AEG","AVR","AVI","AVS","AVG","ATR","ATI","ATS","ATG"/)
ICURR=ICURR+12
END IF
END IF
YRAD(ICURR:ICURR+2)=(/"RHV","PDP","DHV"/)
ICURR=ICURR+3
YRAD(ICURR:ICURR+8)=(/"RHR","RHS","RHG","ZDA","ZDS","ZDG","KDR","KDS","KDG"/)
ICURR=ICURR+9
YRAD(ICURR:ICURR+6)=(/"HAS","M_R","M_I","M_S","M_G","CIT","TEM"/)
ICURR=ICURR+7
YRAD(ICURR:ICURR+2)=(/"RHR","RHS","RHG"/)
ICURR=ICURR+3
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR)="RHH"
ICURR=ICURR+1
END IF
YRAD(ICURR:ICURR+2)=(/"ZDA","ZDS","ZDG"/)
ICURR=ICURR+3
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR)="ZDH"
ICURR=ICURR+1
END IF
YRAD(ICURR:ICURR+2)=(/"KDR","KDS","KDG"/)
ICURR=ICURR+3
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR)="KDH"
ICURR=ICURR+1
END IF
YRAD(ICURR:ICURR+4)=(/"HAS","M_R","M_I","M_S","M_G"/)
ICURR=ICURR+5
IF (CCLOUD=='ICE4') THEN
YRAD(ICURR)="M_H"
ICURR=ICURR+1
END IF
YRAD(ICURR:ICURR+1)=(/"CIT","TEM"/)
ICURR=ICURR+2
IF (CCLOUD=='LIMA') THEN
YRAD(ICURR)="CRT"
ICURR=ICURR+1
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme