From 9a48a2d14aafa9f2dbd1d691d5d66d63a86d0056 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Fri, 26 Mar 2021 16:30:45 +0100
Subject: [PATCH] Quentin 26/03/2021: diag ISOALtitudes correction for cloud
and precipitation and add potential temperature + remove debug print
---
src/MNH/write_lfifm1_for_diag_supp.f90 | 84 ++++++++++++++++----------
src/SURFEX/prep_grib_grid.F90 | 1 -
2 files changed, 52 insertions(+), 33 deletions(-)
diff --git a/src/MNH/write_lfifm1_for_diag_supp.f90 b/src/MNH/write_lfifm1_for_diag_supp.f90
index 1f5bbee65..8bec988c7 100644
--- a/src/MNH/write_lfifm1_for_diag_supp.f90
+++ b/src/MNH/write_lfifm1_for_diag_supp.f90
@@ -1288,41 +1288,61 @@ IF (LISOAL .AND.XISOAL(1)/=0.) THEN
! *********************
! Cloud
! *********************
- IF (SIZE(XRT,4) >= 5) THEN
- ZWORK31(:,:,:) = (XRT(:,:,:,2)+XRT(:,:,:,4)+XRT(:,:,:,5))*1.E3
- CALL ZINTER(ZWORK31, XZZ, ZWAL, ZAL, IIU, IJU, IKU, IKB, IAL, XUNDEF)
- WHERE(ZWAL.EQ.XUNDEF) ZWAL=ZFILLVAL
- TZFIELD%CMNHNAME = 'ALT_CLOUD'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'ALT_CLOUD'
- TZFIELD%CUNITS = 'g kg-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_cloud ALT'
- TZFIELD%NGRID = 1
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZWAL)
- END IF
+ ZWORK31(:,:,:) = 0.
+ IF (SIZE(XRT,4) >= 2) ZWORK31(:,:,:) = XRT(:,:,:,2) ! Rc
+ IF (SIZE(XRT,4) >= 4) ZWORK31(:,:,:) = ZWORK31(:,:,:) + XRT(:,:,:,4) !Ri
+ ZWORK31(:,:,:) = ZWORK31(:,:,:)*1.E3
+ CALL ZINTER(ZWORK31, XZZ, ZWAL, ZAL, IIU, IJU, IKU, IKB, IAL, XUNDEF)
+ WHERE(ZWAL.EQ.XUNDEF) ZWAL=ZFILLVAL
+ TZFIELD%CMNHNAME = 'ALT_CLOUD'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ALT_CLOUD'
+ TZFIELD%CUNITS = 'g kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_cloud ALT'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWAL)
! *********************
! Precipitation
! *********************
- IF (SIZE(XRT,4) >= 6) THEN
- ZWORK31(:,:,:) = (XRT(:,:,:,2)+XRT(:,:,:,6))*1.E3
- CALL ZINTER(ZWORK31, XZZ, ZWAL, ZAL, IIU, IJU, IKU, IKB, IAL, XUNDEF)
- WHERE(ZWAL.EQ.XUNDEF) ZWAL=ZFILLVAL
- TZFIELD%CMNHNAME = 'ALT_PRECIP'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'ALT_PRECIP'
- TZFIELD%CUNITS = 'g kg-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_precipitation ALT'
- TZFIELD%NGRID = 1
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZWAL)
- END IF
+ ZWORK31(:,:,:) = 0.
+ IF (SIZE(XRT,4) >= 3) ZWORK31(:,:,:) = XRT(:,:,:,3) ! Rr
+ IF (SIZE(XRT,4) >= 5) ZWORK31(:,:,:) = ZWORK31(:,:,:) + XRT(:,:,:,5) !Rsnow
+ IF (SIZE(XRT,4) >= 6) ZWORK31(:,:,:) = ZWORK31(:,:,:) + XRT(:,:,:,6) !Rgraupel
+ IF (SIZE(XRT,4) >= 7) ZWORK31(:,:,:) = ZWORK31(:,:,:) + XRT(:,:,:,7) !Rhail
+ ZWORK31(:,:,:) = ZWORK31(:,:,:)*1.E3
+ CALL ZINTER(ZWORK31, XZZ, ZWAL, ZAL, IIU, IJU, IKU, IKB, IAL, XUNDEF)
+ WHERE(ZWAL.EQ.XUNDEF) ZWAL=ZFILLVAL
+ TZFIELD%CMNHNAME = 'ALT_PRECIP'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ALT_PRECIP'
+ TZFIELD%CUNITS = 'g kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_precipitation ALT'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWAL)
+! *********************
+! Potential temperature
+! *********************
+ CALL ZINTER(XTHT, XZZ, ZWAL, ZAL, IIU, IJU, IKU, IKB, IAL, XUNDEF)
+ WHERE(ZWAL.EQ.XUNDEF) ZWAL=ZFILLVAL
+ TZFIELD%CMNHNAME = 'ALT_THETA'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ALT_THETA'
+ TZFIELD%CUNITS = 'K'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_potential temperature ALT'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWAL)
! *********************
! Pressure
! *********************
diff --git a/src/SURFEX/prep_grib_grid.F90 b/src/SURFEX/prep_grib_grid.F90
index de1b37958..b8c5b9c93 100644
--- a/src/SURFEX/prep_grib_grid.F90
+++ b/src/SURFEX/prep_grib_grid.F90
@@ -494,7 +494,6 @@ SELECT CASE (HGRIDTYPE)
IF (IRET == 0 .OR. IMISSING/=1) THEN ! quasi-regular
CALL GRIB_GET(IGRIB,'pl',ININLO_GRIB)
XILO2=360.-360./(MAXVAL(ININLO_GRIB))
- print*,"XILO2=",XILO2
ENDIF
DEALLOCATE(ININLO_GRIB)
ENDIF
--
GitLab