From 669138c8d04a5840edbcd57a07ccd614150049c0 Mon Sep 17 00:00:00 2001
From: Juan Escobar <escj@aero.obs-mip.fr>
Date: Fri, 18 Jun 2021 17:10:08 +0200
Subject: [PATCH] Juan 18/06/2021:ECRAD-1.4.0/1.0.1 management of MNH files

---
 src/MNH/default_desfmn.f90       | 22 ++++++++++++
 src/MNH/ini_modeln.f90           | 33 +++++++++++++++--
 src/MNH/ini_radiations_ecrad.f90 | 13 ++++++-
 src/MNH/modd_param_ecradn.f90    | 35 ++++++++++++++++++
 src/MNH/modd_radiationsn.f90     | 15 ++++++++
 src/MNH/modn_param_ecradn.f90    | 62 ++++++++++++++++++++++++++++++--
 src/MNH/radiations.f90           | 14 +++++++-
 7 files changed, 186 insertions(+), 8 deletions(-)

diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 2e9dd5730..8d7e14740 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -239,6 +239,9 @@ USE MODD_LES
 USE MODD_PARAM_RAD_n
 #ifdef MNH_ECRAD
 USE MODD_PARAM_ECRAD_n
+#if ( VER_ECRAD == 140 ) 
+USE MODD_RADIATIONS_n , ONLY : NSWB_MNH, NLWB_MNH
+#endif
 #endif
 USE MODD_BLANK_n
 USE MODD_FRC
@@ -718,8 +721,27 @@ LFIX_DAT=.FALSE.
 !*      13bis.   SET DEFAULT VALUES FOR MODD_PARAM_ECRAD_n :
 !             ---------------------------------------
 !
+#if ( VER_ECRAD == 101 )
 NSWSOLVER = 0           ! 0: 'McICA 1: 'SPARTACUS' 2: 'SPARTACUS' + 3D effect                            
 NLWSOLVER = 0           ! 0: 'McICA 1: 'SPARTACUS' 2: 'SPARTACUS' + 3D effect 
+#endif
+#if ( VER_ECRAD == 140 )
+LSPEC_ALB = .FALSE.
+LSPEC_EMISS = .FALSE.
+
+
+!ALLOCATE(USER_ALB_DIFF(NSWB_MNH))
+!ALLOCATE(USER_ALB_DIR(NSWB_MNH))
+!ALLOCATE(USER_EMISS(NLWB_MNH))
+!PRINT*,USER_ALB_DIFF
+!USER_ALB_DIFF = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
+!USER_ALB_DIR = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
+!USER_EMISS = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
+SURF_TYPE="SNOW"
+
+NLWSOLVER = 1           ! 0: 'McICA 1: 'SPARTACUS' 2: 'SPARTACUS' + 3D effect 
+NSWSOLVER = 1          ! 0: 'McICA 1: 'SPARTACUS' 2: 'SPARTACUS' + 3D effect                            
+#endif
 ! LEFF3D         = .TRUE.
 ! LSIDEM         = .TRUE.
 NREG           = 3            ! Number of cloudy regions (3=TripleClouds)
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index e7dad9f4e..31e10ae5c 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -472,6 +472,11 @@ USE MODI_SUNPOS_n
 USE MODI_SURF_SOLAR_GEOM
 USE MODI_UPDATE_METRICS
 USE MODI_UPDATE_NSV
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+USE YOERDI   , ONLY :RCCO2
+#endif
+#endif
 !
 IMPLICIT NONE
 !
@@ -1402,15 +1407,31 @@ END IF
 ! Initialization of SW bands
 NSWB_OLD = 6 ! Number of bands in ECMWF original scheme (from Fouquart et Bonnel (1980))
              ! then modified through INI_RADIATIONS_ECMWF but remains equal to 6 practically
+
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+NLWB_OLD = 16 ! For XEMIS initialization (should be spectral in the future)
+#endif
+#endif
+
+NLWB_MNH = 16 ! For XEMIS initialization (should be spectral in the future)
+
 IF (CRAD == 'ECRA') THEN
     NSWB_MNH = 14
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+    NLWB_MNH = 16
+#endif
+#endif
 ELSE
     NSWB_MNH = NSWB_OLD
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+    NLWB_MNH = NLWB_OLD
+#endif
+#endif
 END IF
 
-NLWB_MNH = 16 ! For XEMIS initialization (should be spectral in the future)
-
-
 ALLOCATE(XSW_BANDS (NSWB_MNH)) 
 ALLOCATE(XLW_BANDS (NLWB_MNH)) 
 ALLOCATE(XZENITH   (IIU,IJU))
@@ -2223,8 +2244,14 @@ CALL INI_LW_SETUP (CRAD,NLWB_MNH,XLW_BANDS)
 !
 XCCO2 = 360.0E-06 * 44.0E-03 / XMD
 #ifdef MNH_ECRAD
+#if ( VER_ECRAD == 101 )
 RCCO2 = 360.0E-06 * 44.0E-03 / XMD
 #endif
+#if ( VER_ECRAD == 140 )
+XCCO2 = 0.001*360.0E-06 * 44.0E-03 / XMD
+RCCO2 = 0.001*360.0E-06 * 44.0E-03 / XMD
+#endif
+#endif
 !
 !
 !*      17.2   Externalized surface fields
diff --git a/src/MNH/ini_radiations_ecrad.f90 b/src/MNH/ini_radiations_ecrad.f90
index 6b3fdd2f2..be2571b92 100644
--- a/src/MNH/ini_radiations_ecrad.f90
+++ b/src/MNH/ini_radiations_ecrad.f90
@@ -165,9 +165,14 @@ LCCNO = .FALSE.       ! True if CCN over sea is diagnosed
 
 ! Constant cloud condensation nuclei over land and sea
 ! In ECMWF original code, those values were 900 and 150
+#if ( VER_ECRAD == 101 )
 XCCNLND = 900_JPRB     ! constant CCN over land in m-3 (needed for Martin et al., 1994 parameterization)
 XCCNSEA = 50_JPRB      ! constant CCN over sea in m-3
-
+#endif
+#if ( VER_ECRAD == 140 )
+XCCNLND = 900_JPRB     ! constant CCN over land in cm-3 (needed for Martin et al., 1994 parameterization)
+XCCNSEA = 50_JPRB      ! constant CCN over sea in cm-3 (IFS value, 150 originally in MNH)
+#endif
 ! NAERMACC is in the namelist
 ! NAERMACC = 0  -> Use of Tegen aerosol climatology
 ! NAERMACC = 1  -> Use of MACC aerosol classification
@@ -216,6 +221,12 @@ RCCFC12 = 484.E-12_JPRB
 RCCFC22 =   0.E-12_JPRB
 RCCCL4  =   0.E-12_JPRB
 
+#if ( VER_ECRAD == 140 )
+USER_ALB_DIFF(:) = 0.5
+USER_ALB_DIR(:) = 0 
+USER_EMISS(:) = 1
+#endif
+
 ! Radiation computed every NRADFR timesteps
 NRADFR = INT(XDTRAD/XTSTEP)
 
diff --git a/src/MNH/modd_param_ecradn.f90 b/src/MNH/modd_param_ecradn.f90
index a7bf0b344..a984d405b 100644
--- a/src/MNH/modd_param_ecradn.f90
+++ b/src/MNH/modd_param_ecradn.f90
@@ -42,6 +42,9 @@
 USE MODD_PARAMETERS, ONLY: JPMODELMAX
 USE PARKIND1 , ONLY : JPIM,JPRB
 #ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+USE MODD_RADIATIONS_n , ONLY : NSWB_MNH, NLWB_MNH
+#endif
 USE radiation_config, ONLY : config_type
 #endif
 IMPLICIT NONE
@@ -59,6 +62,18 @@ TYPE PARAM_ECRAD_t
   INTEGER(KIND=JPIM) :: NRADIP           ! 0: 40 mum, 1: Liou and Ou (1994), 2: Liou and Ou (1994) improved,
                                          ! 3: Sun and Rikus (1999)
   REAL(KIND=JPRB)    :: XCLOUD_FRAC_STD  ! Cloud water content horizontal fractional standard deviation in a gridbox  
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+  LOGICAL :: LSPEC_ALB
+  LOGICAL :: LSPEC_EMISS
+  REAL(KIND=JPRB), DIMENSION(14)    :: USER_ALB_DIFF
+  REAL(KIND=JPRB), DIMENSION(14)    :: USER_ALB_DIR
+  REAL(KIND=JPRB), DIMENSION(16)    :: USER_EMISS
+  
+  CHARACTER (LEN=4) :: SURF_TYPE
+#endif
+#endif  
+  
   INTEGER(KIND=JPIM) :: NLWSCATTERING    ! 0: No longwave scattering
                                          ! 1: Longwave scattering by clouds only
                                          ! 2: Longwave scattering by clouds and aerosols
@@ -192,6 +207,16 @@ REAL(KIND=JPRB), POINTER :: XRMINICE=>NULL()
 INTEGER(KIND=JPIM), POINTER :: NMINICE=>NULL() 
 INTEGER(KIND=JPIM), POINTER :: NDECOLAT=>NULL()
 REAL(KIND=JPRB), POINTER :: XCLOUD_FRAC_STD=>NULL()
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+LOGICAL, POINTER :: LSPEC_ALB=>NULL()
+LOGICAL, POINTER :: LSPEC_EMISS=>NULL()
+REAL(KIND=JPRB), DIMENSION(:), POINTER    :: USER_ALB_DIFF=>NULL()
+REAL(KIND=JPRB), DIMENSION(:), POINTER    :: USER_ALB_DIR=>NULL()
+REAL(KIND=JPRB), DIMENSION(:), POINTER    :: USER_EMISS=>NULL()
+CHARACTER(LEN=4), POINTER :: SURF_TYPE
+#endif
+#endif
 !INTEGER, POINTER :: NSW=>NULL() 
 !INTEGER, POINTER :: NSW_EC=>NULL() 
 REAL(KIND=JPRB), POINTER :: XCCH4=>NULL()
@@ -261,6 +286,16 @@ INTEGER, INTENT(IN) :: KFROM, KTO
     NMINICE=>PARAM_ECRAD_MODEL(KTO)%NMINICE
     NDECOLAT=>PARAM_ECRAD_MODEL(KTO)%NDECOLAT
     XCLOUD_FRAC_STD=>PARAM_ECRAD_MODEL(KTO)%XCLOUD_FRAC_STD
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+    LSPEC_ALB=>PARAM_ECRAD_MODEL(KTO)%LSPEC_ALB
+    LSPEC_EMISS=>PARAM_ECRAD_MODEL(KTO)%LSPEC_EMISS
+    USER_ALB_DIFF=>PARAM_ECRAD_MODEL(KTO)%USER_ALB_DIFF
+    USER_ALB_DIR=>PARAM_ECRAD_MODEL(KTO)%USER_ALB_DIR
+    USER_EMISS=>PARAM_ECRAD_MODEL(KTO)%USER_EMISS
+    SURF_TYPE=>PARAM_ECRAD_MODEL(KTO)%SURF_TYPE
+#endif
+#endif    
 !     NSW=>PARAM_ECRAD_MODEL(KTO)%NSW
 !     NSW_EC=>PARAM_ECRAD_MODEL(KTO)%NSW_EC
     XCCH4=>PARAM_ECRAD_MODEL(KTO)%XCCH4   
diff --git a/src/MNH/modd_radiationsn.f90 b/src/MNH/modd_radiationsn.f90
index ec51e9e32..c0953ce5e 100644
--- a/src/MNH/modd_radiationsn.f90
+++ b/src/MNH/modd_radiationsn.f90
@@ -59,6 +59,11 @@ TYPE RADIATIONS_t
   INTEGER :: NRAD    ! number of satellite radiances to synthesize
   INTEGER :: NAER    ! number od AERosol classes
   INTEGER :: NSWB_OLD    ! number of SW bands in ECMWF original code (usually 6)
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )  
+  INTEGER :: NLWB_OLD    ! number of LW bands for emissivity original code (usually 2)
+#endif
+#endif  
   INTEGER :: NSWB_MNH! number of SW bands practically used (14 if ECRAD, NSWB if original code) 
   INTEGER :: NLWB_MNH! number of LW bands practically used (16 if RRTM) 
   INTEGER :: NSTATM  ! index od the STAndard ATMosphere level just above
@@ -132,6 +137,11 @@ INTEGER, POINTER :: NFLUX=>NULL()
 INTEGER, POINTER :: NRAD=>NULL()
 INTEGER, POINTER :: NAER=>NULL()
 INTEGER, POINTER :: NSWB_OLD=>NULL()
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )  
+INTEGER, POINTER :: NLWB_OLD=>NULL()
+#endif
+#endif
 INTEGER, POINTER :: NSWB_MNH=>NULL()
 INTEGER, POINTER :: NLWB_MNH=>NULL()
 INTEGER, POINTER :: NSTATM=>NULL()
@@ -223,6 +233,11 @@ NFLUX=>RADIATIONS_MODEL(KTO)%NFLUX
 NRAD=>RADIATIONS_MODEL(KTO)%NRAD
 NAER=>RADIATIONS_MODEL(KTO)%NAER
 NSWB_OLD=>RADIATIONS_MODEL(KTO)%NSWB_OLD
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+NLWB_OLD=>RADIATIONS_MODEL(KTO)%NLWB_OLD
+#endif
+#endif
 NSWB_MNH=>RADIATIONS_MODEL(KTO)%NSWB_MNH
 NLWB_MNH=>RADIATIONS_MODEL(KTO)%NLWB_MNH
 NSTATM=>RADIATIONS_MODEL(KTO)%NSTATM
diff --git a/src/MNH/modn_param_ecradn.f90 b/src/MNH/modn_param_ecradn.f90
index 096bf9520..1cf2bfa47 100644
--- a/src/MNH/modn_param_ecradn.f90
+++ b/src/MNH/modn_param_ecradn.f90
@@ -14,7 +14,11 @@
 !
 !-------------------------------------------------------------------------------
 USE PARKIND1 , ONLY : JPIM,JPRB
-
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+USE MODD_RADIATIONS_n , ONLY : NSWB_MNH, NLWB_MNH
+#endif
+#endif
 !*       0.   DECLARATIONS
 !             ------------
 !
@@ -30,7 +34,18 @@ USE MODD_PARAM_ECRAD_n, ONLY: &
          NREG_n => NREG, &
          XCLOUD_FRAC_STD_n => XCLOUD_FRAC_STD, &
          NLWSCATTERING_n => NLWSCATTERING, &
-         NAERMACC_n => NAERMACC
+         NAERMACC_n => NAERMACC        
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+USE MODD_PARAM_ECRAD_n, ONLY: &
+         LSPEC_ALB_n => LSPEC_ALB, &
+         LSPEC_EMISS_n => LSPEC_EMISS, &
+       !  USER_ALB_DIFF_n => USER_ALB_DIFF, &
+       !  USER_ALB_DIR_n => USER_ALB_DIR, &
+       !  USER_EMISS_n => USER_EMISS
+         SURF_TYPE_n => SURF_TYPE
+#endif
+#endif
 !          EFF3D_n => EFF3D, &
 !          SIDEM_n => SIDEM, &
 !          LWCSCA_n => LWCSCA, &
@@ -54,6 +69,16 @@ INTEGER(KIND=JPIM), SAVE :: NREG
 REAL(KIND=JPRB), SAVE    :: XCLOUD_FRAC_STD
 INTEGER(KIND=JPIM), SAVE :: NLWSCATTERING
 INTEGER(KIND=JPIM), SAVE :: NAERMACC
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+LOGICAL, SAVE :: LSPEC_ALB
+LOGICAL, SAVE :: LSPEC_EMISS
+!REAL(KIND=JPRB), ALLOCATABLE, SAVE    :: USER_ALB_DIFF(:)
+!REAL(KIND=JPRB), ALLOCATABLE, SAVE    :: USER_ALB_DIR(:)
+!REAL(KIND=JPRB), ALLOCATABLE, SAVE    :: USER_EMISS(:)
+CHARACTER(LEN=4), SAVE   :: SURF_TYPE
+#endif
+#endif
 ! LOGICAL, SAVE :: EFF3D
 ! LOGICAL, SAVE :: SIDEM
 ! LOGICAL, SAVE :: LWCSCA
@@ -64,7 +89,18 @@ INTEGER(KIND=JPIM), SAVE :: NAERMACC
 !
 NAMELIST/NAM_PARAM_ECRADn/NSWSOLVER,NLWSOLVER,NRADLP,NRADIP,&
                           NLIQOPT,NICEOPT,NOVLP,NGAS,NREG,XCLOUD_FRAC_STD,&
-                          NLWSCATTERING, NAERMACC
+                          NLWSCATTERING, NAERMACC &
+#ifndef MNH_ECRAD
+                          & 
+#else
+#if ( VER_ECRAD == 140 )
+                          , LSPEC_ALB, LSPEC_EMISS, &
+!                          USER_ALB_DIFF, USER_ALB_DIR, USER_EMISS, &
+                          SURF_TYPE
+#else
+                          & 
+#endif
+#endif
 !
 CONTAINS
 !
@@ -81,6 +117,16 @@ SUBROUTINE INIT_NAM_PARAM_ECRADn
   XCLOUD_FRAC_STD = XCLOUD_FRAC_STD_n
   NLWSCATTERING = NLWSCATTERING_n
   NAERMACC = NAERMACC_n
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+  LSPEC_ALB = LSPEC_ALB_n
+  LSPEC_EMISS = LSPEC_EMISS_n
+!  USER_ALB_DIFF = USER_ALB_DIFF_n
+!  USER_ALB_DIR = USER_ALB_DIR_n
+!  USER_EMISS = USER_EMISS_n
+  SURF_TYPE = SURF_TYPE_n
+#endif
+#endif  
 END SUBROUTINE INIT_NAM_PARAM_ECRADn
 
 SUBROUTINE UPDATE_NAM_PARAM_ECRADn
@@ -96,6 +142,16 @@ SUBROUTINE UPDATE_NAM_PARAM_ECRADn
   XCLOUD_FRAC_STD_n = XCLOUD_FRAC_STD
   NLWSCATTERING_n = NLWSCATTERING
   NAERMACC_n = NAERMACC
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+  LSPEC_ALB_n = LSPEC_ALB
+  LSPEC_EMISS_n = LSPEC_EMISS
+!  USER_ALB_DIFF_n = USER_ALB_DIFF
+!  USER_ALB_DIR_n = USER_ALB_DIR
+!  USER_EMISS_n = USER_EMISS
+  SURF_TYPE_n = SURF_TYPE
+#endif
+#endif  
 END SUBROUTINE UPDATE_NAM_PARAM_ECRADn
 
 END MODULE MODN_PARAM_ECRAD_n
diff --git a/src/MNH/radiations.f90 b/src/MNH/radiations.f90
index cdb8b13cb..87c377fdd 100644
--- a/src/MNH/radiations.f90
+++ b/src/MNH/radiations.f90
@@ -711,6 +711,11 @@ DO JK=IKUP,KFLEV
   ZQVAVE(:,JK) = 0.5*( PSTATM(JK1,5)/PSTATM(JK1,4)+   &
                  PSTATM(JK2,5)/PSTATM(JK2,4)    )
 END DO
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+ZQVAVE(:,:) = 0.001*ZQVAVE(:,:)
+#endif
+#endif
 !
 !        2.1 pronostic water concentation fields (C2R2 coupling) 
 !
@@ -1100,6 +1105,11 @@ DO JJ=IJB,IJE
     ZO3AVE(IIJ,:)  = POZON (JI,JJ,:)           
   END DO
 END DO
+#ifdef MNH_ECRAD
+#if ( VER_ECRAD == 140 )
+POZON = POZON
+#endif
+#endif
 !
 !-------------------------------------------------------------------------------
 !
@@ -1151,6 +1161,8 @@ ELSE
      END DO
    ENDDO  
 END IF
+!PRINT*,"ZALBP",ZALBP
+!PRINT*,"ZALBD",ZALBD
 !
 !
 ! LW emissivity
@@ -1158,7 +1170,7 @@ ZEMIW(:,:)= ZEMIS(:,:)
 !
 !solar constant
 ZRII0= PCORSOL*XI0  ! solar constant multiplied by seasonal variations due to Earth-Sun distance
-!
+!PRINT*,XI0
 !
 !
 !*       5.2   ACCOUNTS FOR THE CLEAR-SKY APPROXIMATION
-- 
GitLab