Philippe 08/09/2023: sensors: add dry air density Rhod for flyers and stations...

Philippe 08/09/2023: sensors: add dry air density Rhod for flyers and stations (was already available for profilers)

src/MNH/aircraft_balloon_evol.f90

@@ -849,7 +849,6 @@ REAL                           :: ZV_BAL   ! horizontal wind speed at balloon lo
 REAL, DIMENSION(SIZE(PZ,3))    :: ZZ       ! altitude of model levels at station location
 REAL, DIMENSION(SIZE(PR,1),SIZE(PR,2),SIZE(PR,3))    :: ZR
 
-
 TPFLYER%NMODELHIST(ISTORE) = TPFLYER%NMODEL
 
 TPFLYER%XX(ISTORE) = TPFLYER%XX_CUR
@@ -889,6 +888,8 @@ END DO
 
 TPFLYER%XFFZ  (:,ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT( SQRT(PU**2+PV**2) )
 
+TPFLYER%XRHOD (:,ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT( PRHODREF )
+
 IF (CCLOUD=="LIMA") THEN
   TPFLYER%XCIZ  (:,ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT( PSV(:,:,:,NSV_LIMA_NI) )
   TPFLYER%XCCZ  (:,ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT( PSV(:,:,:,NSV_LIMA_NC) )
@@ -906,6 +907,10 @@ CALL Sensor_rare_compute( TPFLYER, ISTORE, PR, PSV, PRHODREF, PCIT, ZTH_EXN, ZZ,
 
 ! vertical wind
 TPFLYER%XWZ  (:,ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT( ZWM(:,:,:) )
+
+! Dry air density at flyer position
+TPFLYER%XRHOD_SENSOR(ISTORE) = TPFLYER%INTERP_FROM_MASSPOINT( PRHODREF )
+
 IF (SIZE(PTKE)>0) TPFLYER%XTKE  (1,ISTORE) = TPFLYER%INTERP_FROM_MASSPOINT( PTKE )
 IF ( CRAD /= 'NONE' ) TPFLYER%XTSRAD(ISTORE) = TPFLYER%INTERP_HOR_FROM_MASSPOINT(PTS )
 TPFLYER%XTKE_DISS(ISTORE) = TPFLYER%INTERP_FROM_MASSPOINT( XCURRENT_TKE_DISS )

src/MNH/modd_sensor.f90

@@ -101,7 +101,9 @@ MODULE MODD_SENSOR
       REAL, DIMENSION(:,:,:), ALLOCATABLE :: XR      ! r*(n)
       REAL, DIMENSION(:,:,:), ALLOCATABLE :: XSV     ! Sv*(n)
       REAL, DIMENSION(:),     ALLOCATABLE :: XTSRAD  ! surface temperature Ts(n)
+      REAL, DIMENSION(:),     ALLOCATABLE :: XRHOD_SENSOR ! density of dry air at sensor position
 
+      REAL, DIMENSION(:,:),   ALLOCATABLE :: XRHOD      ! density of dry air
       REAL, DIMENSION(:,:),   ALLOCATABLE :: XCIZ       ! Ice concentration
       REAL, DIMENSION(:,:),   ALLOCATABLE :: XCCZ       ! Cloud concentration (LIMA)
       REAL, DIMENSION(:,:),   ALLOCATABLE :: XCRZ       ! Rain concentration (LIMA)
@@ -214,8 +216,11 @@ MODULE MODD_SENSOR
       ELSE
         ALLOCATE( TPSENSOR%XTSRAD(0) )
       END IF
+      ALLOCATE( TPSENSOR%XRHOD_SENSOR(KSTORE) ) ; IVARSIZE = IVARSIZE + 1
 
       IF ( OVERTPROF ) THEN
+        ALLOCATE( TPSENSOR%XRHOD(IKU, KSTORE) )      ; IVARSIZE = IVARSIZE + IKU
+
         IF ( CCLOUD == 'LIMA') THEN
           ALLOCATE( TPSENSOR%XCIZ    (IKU, KSTORE) ) ; IVARSIZE = IVARSIZE + IKU
           ALLOCATE( TPSENSOR%XCCZ    (IKU, KSTORE) ) ; IVARSIZE = IVARSIZE + IKU
@@ -235,6 +240,7 @@ MODULE MODD_SENSOR
         ALLOCATE( TPSENSOR%XCRARE    (IKU, KSTORE) ) ; IVARSIZE = IVARSIZE + IKU
         ALLOCATE( TPSENSOR%XCRARE_ATT(IKU, KSTORE) ) ; IVARSIZE = IVARSIZE + IKU
       ELSE
+        ALLOCATE( TPSENSOR%XRHOD   (0, 0) )
         ALLOCATE( TPSENSOR%XCIZ    (0, 0) )
         ALLOCATE( TPSENSOR%XCCZ    (0, 0) )
         ALLOCATE( TPSENSOR%XCRZ    (0, 0) )
@@ -255,6 +261,8 @@ MODULE MODD_SENSOR
       TPSENSOR%XR        (:,:,:) = XUNDEF
       TPSENSOR%XSV       (:,:,:) = XUNDEF
       TPSENSOR%XTSRAD    (:)     = XUNDEF_SFX
+      TPSENSOR%XRHOD_SENSOR(:)   = XNEGUNDEF
+      TPSENSOR%XRHOD     (:,:)   = XNEGUNDEF
       TPSENSOR%XCIZ      (:,:)   = XUNDEF
       TPSENSOR%XCCZ      (:,:)   = XUNDEF
       TPSENSOR%XCRZ      (:,:)   = XUNDEF
@@ -288,6 +296,8 @@ MODULE MODD_SENSOR
       DEALLOCATE( TPSENSOR%XR         )
       DEALLOCATE( TPSENSOR%XSV        )
       DEALLOCATE( TPSENSOR%XTSRAD     )
+      DEALLOCATE( TPSENSOR%XRHOD_SENSOR )
+      DEALLOCATE( TPSENSOR%XRHOD      )
       DEALLOCATE( TPSENSOR%XCIZ       )
       DEALLOCATE( TPSENSOR%XCCZ       )
       DEALLOCATE( TPSENSOR%XCRZ       )
@@ -1002,8 +1012,10 @@ MODULE MODD_SENSOR
       IF ( CRAD /= 'NONE' ) THEN
         PBUFFER(KPOS:KPOS+ISTORES-1) = TPSENSOR%XTSRAD(1:ISTORES) ; KPOS = KPOS + ISTORES
       END IF
+      PBUFFER(KPOS:KPOS+ISTORES-1) = TPSENSOR%XRHOD_SENSOR(1:ISTORES) ; KPOS = KPOS + ISTORES
 
-      IF ( SIZE( TPSENSOR%XIWCZ ) > 0 ) THEN
+      IF ( SIZE( TPSENSOR%XRHOD ) > 0 ) THEN
+        PBUFFER(KPOS:KPOS+IKU*ISTORES-1) = RESHAPE( TPSENSOR%XRHOD     (:,1:ISTORES), [IKU*ISTORES] ) ; KPOS = KPOS + IKU * ISTORES
         IF ( CCLOUD(1:3) == 'ICE') THEN
           PBUFFER(KPOS:KPOS+IKU*ISTORES-1) = RESHAPE( TPSENSOR%XCIZ    (:,1:ISTORES), [IKU*ISTORES] ) ; KPOS = KPOS + IKU * ISTORES
         END IF
@@ -1127,8 +1139,10 @@ MODULE MODD_SENSOR
       IF ( CRAD /= 'NONE' ) THEN
         TPSENSOR%XTSRAD(1:KSTORE) = PBUFFER(KPOS:KPOS+KSTORE-1) ; KPOS = KPOS + KSTORE
       END IF
+      TPSENSOR%XRHOD_SENSOR(1:KSTORE) = PBUFFER(KPOS:KPOS+KSTORE-1) ; KPOS = KPOS + KSTORE
 
-      IF ( SIZE( TPSENSOR%XIWCZ ) > 0 ) THEN
+      IF ( SIZE( TPSENSOR%XRHOD ) > 0 ) THEN
+        TPSENSOR%XRHOD     (:,1:KSTORE) = RESHAPE( PBUFFER(KPOS:KPOS+IKU*KSTORE-1), [ IKU, KSTORE ] ) ; KPOS = KPOS + IKU * KSTORE
         IF ( CCLOUD(1:3) == 'ICE' ) THEN
           TPSENSOR%XCIZ    (:,1:KSTORE) = RESHAPE( PBUFFER(KPOS:KPOS+IKU*KSTORE-1), [ IKU, KSTORE ] ) ; KPOS = KPOS + IKU * KSTORE
         END IF

src/MNH/modd_type_statprof.f90

@@ -110,7 +110,6 @@ TYPE, EXTENDS( TSTATPROFDATA ) ::  TPROFILERDATA
   REAL, DIMENSION(:,:),   ALLOCATABLE :: XTHV       ! thv(n)
   REAL, DIMENSION(:,:),   ALLOCATABLE :: XVISIGUL   ! VISI GULTEPE(n)
   REAL, DIMENSION(:,:),   ALLOCATABLE :: XVISIKUN   ! VISI KUNKEL(n)
-  REAL, DIMENSION(:,:),   ALLOCATABLE :: XRHOD      ! density of dry air/moist air
   REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAER       ! AER*(n) aerosol extinction
 
   REAL, DIMENSION(:), ALLOCATABLE :: XIWV ! integrated water vpour(n)
@@ -322,7 +321,6 @@ SUBROUTINE DATA_ARRAYS_ALLOCATE_PROFILER( TPSENSOR, KSTORE )
     ELSE
       ALLOCATE( TPSENSOR%XVISIKUN  (0, 0) )
     END IF
-    ALLOCATE( TPSENSOR%XRHOD     (IKU, KSTORE) )        ; IVARSIZE = IVARSIZE + IKU
     ALLOCATE( TPSENSOR%XAER      (IKU, KSTORE, NAER ) ) ; IVARSIZE = IVARSIZE + IKU * NAER
 
     ALLOCATE( TPSENSOR%XIWV(KSTORE) ) ; IVARSIZE = IVARSIZE + 1
@@ -341,7 +339,6 @@ SUBROUTINE DATA_ARRAYS_ALLOCATE_PROFILER( TPSENSOR, KSTORE )
     TPSENSOR%XTHV      (:,:) = XUNDEF
     IF ( CCLOUD == 'C2R2' .OR. CCLOUD == 'KHKO' )  TPSENSOR%XVISIGUL(:,:) = XUNDEF
     IF ( CCLOUD /= 'NONE' .AND. CCLOUD /= 'REVE' ) TPSENSOR%XVISIKUN(:,:) = XUNDEF
-    TPSENSOR%XRHOD     (:,:) = XUNDEF
     TPSENSOR%XAER      (:,:,:) = XUNDEF
 
     TPSENSOR%XIWV(:) = XUNDEF
@@ -370,7 +367,6 @@ SUBROUTINE DATA_ARRAYS_ALLOCATE_PROFILER( TPSENSOR, KSTORE )
     DEALLOCATE( TPSENSOR%XTHV       )
     DEALLOCATE( TPSENSOR%XVISIGUL   )
     DEALLOCATE( TPSENSOR%XVISIKUN   )
-    DEALLOCATE( TPSENSOR%XRHOD      )
     DEALLOCATE( TPSENSOR%XAER       )
 
     DEALLOCATE( TPSENSOR%XIWV )
@@ -524,7 +520,6 @@ SUBROUTINE DATA_ARRAYS_ALLOCATE_PROFILER( TPSENSOR, KSTORE )
     IF ( CCLOUD /= 'NONE' .AND. CCLOUD /= 'REVE' ) THEN
       PBUFFER(KPOS:KPOS+IKU*ISTORES-1) = RESHAPE( TPSENSOR%XVISIKUN(:,1:ISTORES), [IKU*ISTORES] ) ; KPOS = KPOS + IKU * ISTORES
     END IF
-    PBUFFER(KPOS:KPOS+IKU*ISTORES-1) = RESHAPE( TPSENSOR%XRHOD(:,1:ISTORES), [IKU*ISTORES] ) ; KPOS = KPOS + IKU * ISTORES
     PBUFFER(KPOS:KPOS+IKU*ISTORES*NAER-1) = RESHAPE( TPSENSOR%XAER(:,1:ISTORES,:), [IKU*ISTORES*NAER] )
     KPOS = KPOS + IKU * ISTORES * NAER
 
@@ -650,7 +645,6 @@ SUBROUTINE DATA_ARRAYS_ALLOCATE_PROFILER( TPSENSOR, KSTORE )
     IF ( CCLOUD /= 'NONE' .AND. CCLOUD /= 'REVE' ) THEN
       TPSENSOR%XVISIKUN(:,1:KSTORE) = RESHAPE( PBUFFER(KPOS:KPOS+IKU*KSTORE-1), [ IKU, KSTORE ] ) ; KPOS = KPOS + IKU * KSTORE
     END IF
-    TPSENSOR%XRHOD(:,1:KSTORE) = RESHAPE( PBUFFER(KPOS:KPOS+IKU*KSTORE-1), [ IKU, KSTORE ] ) ; KPOS = KPOS + IKU * KSTORE
     TPSENSOR%XAER(:,1:KSTORE,:) = RESHAPE( PBUFFER(KPOS:KPOS+IKU*KSTORE*NAER-1), [ IKU, KSTORE, NAER ] )
     KPOS = KPOS + IKU * KSTORE * NAER
 

src/MNH/modeln.f90

@@ -2190,7 +2190,7 @@ END IF
 !               --------------------------------
 !
 IF ( LSTATION ) &
-  CALL STATION_n( XZZ, XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST )
+  CALL STATION_n( XZZ, XRHODREF, XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST )
 !
 !---------------------------------------------------------
 !

src/MNH/profilern.f90

@@ -366,6 +366,12 @@ PROFILER: DO JP = 1, NUMBPROFILER_LOC
     TPROFILERS(JP)%XAER(:,IN,JSV) = TPROFILERS(JP)%INTERP_HOR_FROM_MASSPOINT( ZWORK2(:,:,:,JSV) )
   END DO
   IF (SIZE(PTKE)>0) TPROFILERS(JP)%XTKE  (:,IN) = TPROFILERS(JP)%INTERP_HOR_FROM_MASSPOINT( PTKE )
+
+  ! XRHOD_SENSOR is not computed for profilers because not very useful
+  ! If needed, the interpolation must also be done vertically
+  ! (and therefore the vertical interpolation coefficients have to be computed)
+  ! TPROFILERS(JP)%XRHOD_SENSOR(IN) = ...
+
   IF ( CRAD /= 'NONE' ) TPROFILERS(JP)%XTSRAD(IN) = TPROFILERS(JP)%INTERP_HOR_FROM_MASSPOINT( PTS )
   !
   IF ( LDIAG_SURFRAD_PROF ) CALL STATPROF_DIAG_SURFRAD(TPROFILERS(JP), IN )

src/MNH/stationn.f90

@@ -9,11 +9,12 @@ MODULE MODI_STATION_n
 !
 INTERFACE
 !
-      SUBROUTINE STATION_n( PZ,                             &
+      SUBROUTINE STATION_n( PZ, PRHODREF,                   &
                             PU, PV, PW, PTH, PR, PSV, PTKE, &
                             PTS, PP )
 !
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
+REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF ! dry air density of the reference state
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PW     ! vertical wind
@@ -33,7 +34,7 @@ END INTERFACE
 END MODULE MODI_STATION_n
 !
 !     #######################################################
-      SUBROUTINE STATION_n( PZ,                             &
+      SUBROUTINE STATION_n( PZ, PRHODREF,                   &
                             PU, PV, PW, PTH, PR, PSV, PTKE, &
                             PTS, PP )
 !     #######################################################
@@ -101,6 +102,7 @@ IMPLICIT NONE
 !
 !
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
+REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF ! dry air density of the reference state
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PW     ! vertical wind
@@ -165,6 +167,8 @@ STATION: DO JS = 1, NUMBSTAT_LOC
     TSTATIONS(JS)%XSV(1,IN,JSV) = TSTATIONS(JS)%INTERP_HOR_FROM_MASSPOINT( PSV(:,:,JK,JSV) )
   END DO
 
+  TSTATIONS(JS)%XRHOD_SENSOR(IN) = TSTATIONS(JS)%INTERP_HOR_FROM_MASSPOINT( PRHODREF(:,:,JK) )
+
   IF (SIZE(PTKE)>0) TSTATIONS(JS)%XTKE(1,IN) = TSTATIONS(JS)%INTERP_HOR_FROM_MASSPOINT( PTKE(:,:,JK) )
   IF ( CRAD /= 'NONE' ) TSTATIONS(JS)%XTSRAD(IN) = TSTATIONS(JS)%INTERP_HOR_FROM_MASSPOINT( PTS )
   TSTATIONS(JS)%XZS  = TSTATIONS(JS)%INTERP_HOR_FROM_MASSPOINT( PZ(:,:,1+JPVEXT))

src/MNH/write_aircraft_balloon.f90

@@ -205,7 +205,7 @@ USE MODD_SALT,             ONLY: LSALT, NMODE_SLT
 
 use mode_aircraft_balloon, only: Aircraft_balloon_longtype_get
 USE MODE_MODELN_HANDLER,   ONLY: GET_CURRENT_MODEL_INDEX
-use mode_sensor,           only: Add_dust_data, Add_fixpoint, Add_orilam_data, Add_point, Add_profile, Add_salt_data, &
+use mode_sensor,           only: Add_dust_data, Add_orilam_data, Add_point, Add_profile, Add_salt_data, &
                                  Sensor_current_processes_number_get, &
                                  ccomment, ctitle, cunit, xwork6, &
                                  Sensor_write_workarrays_allocate, Sensor_write_workarrays_deallocate
@@ -246,7 +246,7 @@ IF ( IMI /= TPFLYER%NMODEL ) RETURN
 !
 IKU = SIZE(TPFLYER%XRTZ,1) !number of vertical levels
 !
-IPROC = 9 + IRR + SIZE(TPFLYER%XSV,3) + 2 + SIZE(TPFLYER%XSVW_FLUX,2)
+IPROC = 10 + IRR + SIZE(TPFLYER%XSV,3) + 2 + SIZE(TPFLYER%XSVW_FLUX,2)
 IF ( IRR > 1 ) IPROC = IPROC + 1
 IF ( SIZE( TPFLYER%XTKE ) > 0 ) IPROC = IPROC + 1
 IPROC = IPROC + 1 ! TKE_DISS
@@ -318,7 +318,9 @@ IF ( IRR > 1 ) THEN !cloud water is present
   call Add_point( 'LWC', 'cloud liquid water content', 'g m-3', ZLWC(:) )
   DEALLOCATE( ZLWC, ZRHO )
 END IF
-!
+
+call Add_point( 'Rhod', 'Density of dry air', 'kg m-3', tpflyer%xrhod_sensor )
+
 IF (SIZE(TPFLYER%XTKE)>0) call Add_point( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tpflyer%xtke(1,:) )
 !
 call Add_point( 'H_FLUX',  'sensible flux', 'W m-2', tpflyer%xthw_flux(:) )
@@ -420,7 +422,7 @@ call Sensor_write_workarrays_deallocate( )
 !----------------------------------------------------------------------------
 !Treat vertical profiles
 
-IPROCZ = 8 + IRR
+IPROCZ = 9 + IRR
 IF ( CCLOUD == 'LIMA' )     IPROCZ = IPROCZ + 3
 IF ( CCLOUD(1:3) == 'ICE' ) IPROCZ = IPROCZ + 1
 
@@ -441,6 +443,7 @@ call Add_profile( 'FF', 'Horizontal wind', 'm s-1', tpflyer%xffz(:,:) )
 call Add_profile( 'IWC', 'Ice water content',    'kg m-3', tpflyer%xiwcz(:,:) )
 call Add_profile( 'LWC', 'Liquid water content', 'kg m-3', tpflyer%xlwcz(:,:) )
 
+call Add_profile( 'Rhod', 'Density of dry air',  'kg m-3', tpflyer%xrhod(:,:) )
 IF ( CCLOUD == 'LIMA' ) THEN
   call Add_profile( 'CCLOUDT', 'liquid cloud concentration', 'kg-1', tpflyer%xccz(:,:) )
   call Add_profile( 'CRAINT',  'Rain concentration',         'kg-1', tpflyer%xcrz(:,:) )

src/MNH/write_profilern.f90

@@ -145,7 +145,7 @@ if ( nrr >= 5 ) call Add_profile( 'Rs', 'Snow mixing ratio',               'kg k
 if ( nrr >= 6 ) call Add_profile( 'Rg', 'Graupel mixing ratio',            'kg kg-1', tpprofiler%xr(:,:,6) )
 if ( nrr >= 7 ) call Add_profile( 'Rh', 'Hail mixing ratio',               'kg kg-1', tpprofiler%xr(:,:,7) )
 
-call Add_profile( 'Rhod', 'Density of dry air in moist', 'kg m-3', tpprofiler%xrhod )
+call Add_profile( 'Rhod', 'Density of dry air', 'kg m-3', tpprofiler%xrhod )
 if ( cturb == 'TKEL') &
   call Add_profile( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tpprofiler%xtke )
 
@@ -265,6 +265,10 @@ IF( CRAD /= 'NONE' )  IPROC = IPROC + 1 !Tsrad term
 
 call Sensor_write_workarrays_allocate( 1, istore, iproc )
 
+! xrhod_sensor NOT computed => not written
+! if needed, please add its computation in STATION_n
+! call Add_point( 'Rhod', 'Density of dry air', 'kg m-3', tpstation%xrhod_sensor )
+
 if ( ldiag_surfrad_prof ) call Add_diag_surfrad_data( tpprofiler )
 
 call Add_point( 'IWV', 'Integrated Water Vapour',   'kg m-2', tpprofiler%xiwv )

src/MNH/write_stationn.f90

@@ -70,7 +70,7 @@ type(tfieldmetadata_base), dimension(:), allocatable :: tzfields
 !
 !----------------------------------------------------------------------------
 !
-IPROC = 5 + SIZE(TPSTATION%XR,3) + SIZE(TPSTATION%XSV,3)
+IPROC = 6 + SIZE(TPSTATION%XR,3) + SIZE(TPSTATION%XSV,3)
 
 IF ( CTURB == 'TKEL' ) IPROC = IPROC + 1
 IF (LDIAG_SURFRAD_STAT) THEN
@@ -123,6 +123,8 @@ do jrr = 1, SIZE( tpstation%xr, 3 )
   end select
 end do
 
+call Add_point( 'Rhod', 'Density of dry air', 'kg m-3', tpstation%xrhod_sensor )
+
 if ( cturb == 'TKEL' ) call Add_point( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tpstation%xtke(1,:) )
 
 if ( nsv > 0 ) then