From a85b2491e646d52201f9ff617c9344fe2e85c655 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Mon, 17 Jun 2024 14:47:22 +0200
Subject: [PATCH] Quentin 17/06/2024: Cleaning old version of ECRAD

---
 .../ifs/cloud_overlap_decorr_len.F90          | 140 ----
 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cos_sza.F90   | 345 ----------
 .../ifs/ice_effective_radius.F90              | 196 ------
 .../ifs/liquid_effective_radius.F90           | 300 ---------
 .../ecrad-1.0.1_mnh/ifs/radiation_scheme.F90  | 625 ------------------
 .../ecrad-1.0.1_mnh/ifs/radiation_setup.F90   | 368 -----------
 src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomcst.F90 |  29 -
 src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomlun.F90 |  15 -
 .../RAD/ecrad-1.0.1_mnh/ifsrrtm/rrtm_kgb1.F90 | 355 ----------
 .../ifsrrtm/srtm_gas_optical_depth.F90        | 339 ----------
 .../ecrad-1.0.1_mnh/ifsrrtm/srtm_kgb16.F90    | 184 ------
 11 files changed, 2896 deletions(-)
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cloud_overlap_decorr_len.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cos_sza.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/ice_effective_radius.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/liquid_effective_radius.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_scheme.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_setup.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomcst.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomlun.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/rrtm_kgb1.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_gas_optical_depth.F90
 delete mode 100644 src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_kgb16.F90

diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cloud_overlap_decorr_len.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cloud_overlap_decorr_len.F90
deleted file mode 100644
index a0c35ec77..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cloud_overlap_decorr_len.F90
+++ /dev/null
@@ -1,140 +0,0 @@
-!      #################################
-       MODULE MODI_CLOUD_OVERLAP_DECORR_LEN
-!      #################################
-INTERFACE
-
-SUBROUTINE CLOUD_OVERLAP_DECORR_LEN &
-     & (KIDIA, KFDIA, KLON, PGEMU, NDECOLAT, &
-     &  PDECORR_LEN_EDGES_KM, PDECORR_LEN_WATER_KM, PDECORR_LEN_RATIO)
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON) ! Sine of latitude
-INTEGER(KIND=JPIM),INTENT(IN) :: NDECOLAT
-REAL(KIND=JPRB), INTENT(OUT)           :: PDECORR_LEN_EDGES_KM(KLON)
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_WATER_KM(KLON)
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_RATIO
-
-END SUBROUTINE CLOUD_OVERLAP_DECORR_LEN
-END INTERFACE
-END MODULE MODI_CLOUD_OVERLAP_DECORR_LEN
-
-SUBROUTINE CLOUD_OVERLAP_DECORR_LEN &
-     & (KIDIA, KFDIA, KLON, PGEMU, NDECOLAT, &
-     &  PDECORR_LEN_EDGES_KM, PDECORR_LEN_WATER_KM, PDECORR_LEN_RATIO)
-
-! CLOUD_OVERLAP_DECORR_LEN
-!
-! PURPOSE
-! -------
-!   Calculate the cloud overlap decorrelation length as a function of
-!   latitude for use in the radiation scheme
-!
-! INTERFACE
-! ---------
-!   CLOUD_OVERLAP_DECORR_LEN is called from RADLSWR and RADIATION_SCHEME
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF (using code extracted from radlswr.F90)
-!   Original: 2016-02-16
-!
-! MODIFICATIONS
-! -------------
-!
-! -------------------------------------------------------------------
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-USE YOMCST   , ONLY : RPI
-USE MODD_PARAM_ECRAD_N   , ONLY : XDECORR_CF,XDECORR_CW
-
-! -------------------------------------------------------------------
-
-IMPLICIT NONE
-
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-
-! *** Single-level variables 
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON) ! Sine of latitude
-
-! *** Configuration variable controlling the overlap scheme
-INTEGER(KIND=JPIM),INTENT(IN) :: NDECOLAT
-
-! OUTPUT ARGUMENTS
-
-! *** Decorrelation lengths for cloud edges and cloud water content,
-! *** in km
-REAL(KIND=JPRB), INTENT(OUT)           :: PDECORR_LEN_EDGES_KM(KLON)
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_WATER_KM(KLON)
-  
-! Ratio of water-content to cloud-edge decorrelation lengths
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_RATIO
-
-! LOCAL VARIABLES
-
-REAL(KIND=JPRB) :: ZRADIANS_TO_DEGREES, ZABS_LAT_DEG, ZCOS_LAT
-
-INTEGER(KIND=JPIM) :: JL
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',0,ZHOOK_HANDLE)
-  
-! -------------------------------------------------------------------
-IF (NDECOLAT == 0) THEN
-
-  ! Decorrelation lengths are constant values
-  PDECORR_LEN_EDGES_KM(KIDIA:KFDIA) = XDECORR_CF
-  IF (PRESENT(PDECORR_LEN_WATER_KM)) THEN
-    PDECORR_LEN_WATER_KM(KIDIA:KFDIA) = XDECORR_CW
-  ENDIF
-  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
-    PDECORR_LEN_RATIO = XDECORR_CW / XDECORR_CF
-  ENDIF
-
-ELSE
-
-  ZRADIANS_TO_DEGREES = 180.0_JPRB / RPI
-
-  IF (NDECOLAT == 1) THEN
-    ! Shonk et al. (2010) Eq. 13 formula
-    DO JL = KIDIA,KFDIA
-      ZABS_LAT_DEG = ABS(ASIN(PGEMU(JL)) * ZRADIANS_TO_DEGREES)
-      PDECORR_LEN_EDGES_KM(JL) = 2.899_JPRB - 0.02759_JPRB * ZABS_LAT_DEG
-    ENDDO
-  ELSE ! NDECOLAT == 2
-    DO JL = KIDIA,KFDIA
-      ! Shonk et al. (2010) but smoothed over the equator
-      ZCOS_LAT = COS(ASIN(PGEMU(JL)))
-      PDECORR_LEN_EDGES_KM(JL) = 0.75_JPRB + 2.149_JPRB * ZCOS_LAT*ZCOS_LAT
-    ENDDO
-  ENDIF
-
-  ! Both NDECOLAT = 1 and 2 assume that the decorrelation length for
-  ! cloud water content is half that for cloud edges
-  IF (PRESENT(PDECORR_LEN_WATER_KM)) THEN
-    PDECORR_LEN_WATER_KM(KIDIA:KFDIA) = PDECORR_LEN_EDGES_KM(KIDIA:KFDIA) * 0.5_JPRB
-  ENDIF
-
-  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
-    PDECORR_LEN_RATIO = 0.5_JPRB
-  ENDIF
-
-ENDIF
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',1,ZHOOK_HANDLE)
-
-END SUBROUTINE CLOUD_OVERLAP_DECORR_LEN
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cos_sza.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cos_sza.F90
deleted file mode 100644
index 215743cfb..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/cos_sza.F90
+++ /dev/null
@@ -1,345 +0,0 @@
-SUBROUTINE COS_SZA(KSTART,KEND,KCOL,PGEMU,PGELAM,LDRADIATIONTIMESTEP,PMU0)
-
-!**** *COS_SZA*   
-
-!     Purpose.
-!     --------
-!        Compute the cosine of the solar zenith angle.  Note that this
-!        is needed for three different things: (1) as input to the
-!        radiation scheme in which it is used to compute the path
-!        length of the direct solar beam through the atmosphere, (2)
-!        every timestep to scale the solar fluxes by the incoming
-!        solar radiation at top-of-atmosphere, and (3) to compute the
-!        albedo of the ocean.  For (1) we ideally want an average
-!        value for the duration of a radiation timestep while for (2)
-!        we want an average value for the duration of a model
-!        timestep.
-
-!**   Interface.
-!     ----------
-!        *CALL* *COS_SZA(...)
-
-!        Explicit arguments : 
-!        ------------------
-!            PGEMU - Sine of latitude
-!            PGELAM - Geographic longitude in radians
-!            LDRadiationTimestep - Is this for a radiation timestep?
-!            PMU0 - Output cosine of solar zenith angle
-
-!        Implicit arguments :
-!        --------------------
-!            YRRIP%RWSOVR, RWSOVRM - Solar time for model/radiation timesteps
-!            RCODECM, RSIDECM - Sine/cosine of solar declination
-!            YRERAD%LAverageSZA - Average solar zenith angle in time interval?
-!            YRRIP%TSTEP - Model timestep in seconds
-!            YRERAD%NRADFR - Radiation frequency in timesteps
-
-!     Method.
-!     -------
-!        Compute cosine of the solar zenith angle, mu0, from lat, lon
-!        and solar time using standard formula.  If
-!        YRERAD%LAverageSZA=FALSE then this is done at a single time,
-!        which is assumed to be the mid-point of either the model or
-!        the radiation timestep.  If YRERAD%LAverageSZA=TRUE then we
-!        compute the average over the model timestep exactly by first
-!        computing sunrise/sunset times. For radiation timesteps, mu0
-!        is to be used to compute the path length of the direct solar
-!        beam through the atmosphere, and the fluxes are subsequently
-!        weighted by mu0.  Therefore night-time values are not used,
-!        so we average mu0 only when the sun is above the horizon.
-
-!     Externals.
-!     ----------
-
-!     Reference.
-!     ----------
-!        ECMWF Research Department documentation of the IFS
-!
-!        See also: Zhou, L., M. Zhang, Q. Bao, and Y. Liu (2015), On
-!        the incident solar radiation in CMIP5
-!        models. Geophys. Res. Lett., 42, 1930–1935. doi:
-!        10.1002/2015GL063239.
-
-!     Author.
-!     -------
-!      Robin Hogan, ECMWF, May 2015
-
-!     Modifications:
-!     --------------
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-USE YOMCST   , ONLY : RPI, RDAY
-! USE YOMRIP   , ONLY : YRRIP
-USE MODD_RADIATIONS_N , ONLY : XCOSDEL, XSINDEL, XZENITH ! no difference between yoerip and yomrip here
-! USE YOERIP   , ONLY : YRERIP
-USE MODD_PARAM_ECRAD_N  , ONLY : NRADFR, LCENTREDTIMESZA, LAVERAGESZA
-USE YOMLUN   , ONLY : NULOUT
-
-USE MODD_PARAM_RAD_n,  ONLY : XDTRAD
-USE MODD_RADIATIONS_n, ONLY : XSINDEL, XCOSDEL, XTSIDER
-USE MODD_TIME_n,       ONLY : TDTRAD_FULL
-USE MODD_DYN_n,       ONLY : XTSTEP
-
-!     ------------------------------------------------------------------
-
-IMPLICIT NONE
-
-INTEGER(KIND=JPIM),INTENT(IN) :: KSTART      ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KEND        ! Last column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KCOL        ! Number of columns in arrays
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KCOL) ! Sine of latitude
-REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KCOL)! Longitude in radians
-LOGICAL, INTENT(IN) :: LDRADIATIONTIMESTEP   ! Is this for a radiation timestep?
-REAL(KIND=JPRB),  INTENT(OUT) :: PMU0(KCOL)  ! Cosine of solar zenith angle
-
-! Solar time at the start and end of the time interval
-REAL(KIND=JPRB) :: ZSOLARTIMESTART, ZSOLARTIMEEND
-
-! The time of half a model/radiation timestep, in radians
-REAL(KIND=JPRB) :: ZHALFTIMESTEP
-
-! For efficiency we precompute sin(solar declination)*sin(latitude)
-REAL(KIND=JPRB) :: ZSINDECSINLAT(KSTART:KEND)
-!...and cos(solar declination)*cos(latitude)
-REAL(KIND=JPRB) :: ZCOSDECCOSLAT(KSTART:KEND)
-! ...and cosine of latitude
-REAL(KIND=JPRB) :: ZCOSLAT(KSTART:KEND)
-
-! Tangent of solar declination
-REAL(KIND=JPRB) :: ZTANDEC
-
-! Hour angles (=local solar time in radians plus pi)
-REAL(KIND=JPRB) :: ZHOURANGLESTART, ZHOURANGLEEND
-REAL(KIND=JPRB) :: ZHOURANGLESUNSET, ZCOSHOURANGLESUNSET
-
-INTEGER(KIND=JPIM) :: JCOL        ! Column index
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-REAL(KIND=JPRB) :: ZTIME,ZUT
-REAL(KIND=JPRB) :: ZTUT,ZSOLANG
-
-!     ------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('COS_SZA',0,ZHOOK_HANDLE)
-
-! An average solar zenith angle can only be computed if the solar time
-! is centred on the time interval
-IF (LAVERAGESZA .AND. .NOT. LCENTREDTIMESZA) THEN
-  WRITE(NULOUT,*) 'ERROR IN COS_SZA: LAverageSZA=TRUE but LCentredTimeSZA=FALSE'
-  CALL ABOR1('COS_SZA: ABOR1 CALLED')
-ENDIF
-
-DO JCOL = KSTART,KEND
-  ZCOSLAT(JCOL) = SQRT(1.0_JPRB - PGEMU(JCOL)**2)
-ENDDO
-
-! Computation of solar hour angle from sunposn
-ZTIME     = TDTRAD_FULL%XTIME + 0.5*XDTRAD
-ZUT       = MOD( 24.0+MOD(ZTIME/3600.,24.0),24.0 )
-ZTUT = ZUT - XTSIDER
-ZSOLANG = (ZTUT-12.0)*15.0*(RPI/180.)          ! hour angle in radians
-
-IF (LDRADIATIONTIMESTEP) THEN
-  ! Compute the effective cosine of solar zenith angle for a radiation
-  ! timestep
-
-  ! Precompute quantities that may be used more than once
-  DO JCOL = KSTART,KEND
-    ZSINDECSINLAT(JCOL) = XSINDEL * PGEMU(JCOL)  ! this is the current declination, not that of radiation (averaged over time step)
-    ZCOSDECCOSLAT(JCOL) = XCOSDEL * ZCOSLAT(JCOL)
-  ENDDO
-
-  IF (.NOT. LAVERAGESZA) THEN
-    ! Original method: compute the value at the centre of the
-    ! radiation timestep (assuming that LCentredTimeSZA=TRUE - see
-    ! updtim.F90)
-    DO JCOL = KSTART,KEND
-      ! It would be more efficient to do it like this...
-      ! PMU0(JCOL)=MAX(0.0_JPRB, ZSinDecSinLat(JCOL) &
-      !      & - ZCosDecCosLat(JCOL) * COS(YRERIP%RWSOVRM + PGELAM(JCOL)))
-      ! ...but for bit reproducibility with previous cycle we do it
-      ! like this:
-      PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) &
-           & - XCOSDEL*COS(ZSOLANG)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) &
-           & + XCOSDEL*SIN(ZSOLANG)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL)))
-    ENDDO
-
-  ELSE
-    ! Compute the average MU0 for the period of the radiation
-    ! timestep, excluding times when the sun is below the horizon
-
-    ! First compute the sine and cosine of the times of the start and
-    ! end of the radiation timestep
-    ZHALFTIMESTEP = XTSTEP * REAL(NRADFR) * RPI / RDAY
-    ZSOLARTIMESTART = ZSOLANG - ZHALFTIMESTEP
-    ZSOLARTIMEEND   = ZSOLANG + ZHALFTIMESTEP
-
-    ! Compute tangent of solar declination, with check in case someone
-    ! simulates a planet completely tipped over
-    ZTANDEC = XSINDEL / MAX(XCOSDEL, 1.0E-12)
-
-    DO JCOL = KSTART,KEND
-      ! Sunrise equation: cos(hour angle at sunset) =
-      ! -tan(declination)*tan(latitude)
-      ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) &
-           &              / MAX(ZCOSLAT(JCOL), 1.0E-12)
-      IF (ZCOSHOURANGLESUNSET > 1.0) THEN
-        ! Perpetual darkness
-        PMU0(JCOL) = 0.0_JPRB
-      ELSE
-        ! Compute hour angle at start and end of time interval,
-        ! ensuring that the hour angle of the centre of the time
-        ! window is in the range -PI to +PI (equivalent to ensuring
-        ! that local solar time = solar time + longitude is in the
-        ! range 0 to 2PI)
-        IF (ZSOLANG + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN
-          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI
-          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - RPI 
-        ELSE
-          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI
-          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - 3.0_JPRB*RPI
-        ENDIF
-
-        IF (ZCOSHOURANGLESUNSET >= -1.0) THEN
-          ! Not perpetual daylight or perpetual darkness, so we need
-          ! to check for sunrise or sunset lying within the time
-          ! interval
-          ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET)
-          IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET &
-               & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN
-            ! The time interval is either completely before sunrise or
-            ! completely after sunset
-            PMU0(JCOL) = 0.0_JPRB
-            CYCLE
-          ENDIF
-
-          ! Bound the start and end hour angles by sunrise and sunset
-          ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, &
-               &                MIN(ZHOURANGLESTART, ZHOURANGLESUNSET))
-          ZHOURANGLEEND   = MAX(-ZHOURANGLESUNSET, &
-               &                MIN(ZHOURANGLEEND,   ZHOURANGLESUNSET))
-        ENDIF
-
-        IF (ZHOURANGLEEND - ZHOURANGLESTART > 1.0E-8) THEN
-          ! Compute average MU0 in the interval ZHourAngleStart to
-          ! ZHourAngleEnd
-          PMU0(JCOL) = ZSINDECSINLAT(JCOL) &
-               & + (ZCOSDECCOSLAT(JCOL) &
-               &    * (SIN(ZHOURANGLEEND) - SIN(ZHOURANGLESTART))) &
-               & / (ZHOURANGLEEND - ZHOURANGLESTART)
-
-          ! Just in case...
-          IF (PMU0(JCOL) < 0.0_JPRB) THEN
-            PMU0(JCOL) = 0.0_JPRB
-          ENDIF
-        ELSE
-          ! Too close to sunrise/sunset for a reliable calculation
-          PMU0(JCOL) = 0.0_JPRB
-        ENDIF
-
-      ENDIF
-    ENDDO
-  ENDIF
-
-ELSE
-  ! Compute the cosine of solar zenith angle for a model timestep
-
-  ! Precompute quantities that may be used more than once
-  DO JCOL = KSTART,KEND
-    ZSINDECSINLAT(JCOL) = XSINDEL * PGEMU(JCOL)
-    ZCOSDECCOSLAT(JCOL) = XCOSDEL * ZCOSLAT(JCOL)
-  ENDDO
-
-  IF (.NOT. LAVERAGESZA) THEN
-    ! Original method: compute the value at the centre of the
-    ! model timestep
-    DO JCOL = KSTART,KEND
-      ! It would be more efficient to do it like this...
-      ! PMU0(JCOL) = MAX(0.0_JPRB, ZSinDecSinLat(JCOL)        &
-      !      & - ZCosDecCosLat(JCOL)*COS(YRRIP%RWSOVR + PGELAM(JCOL)))
-      ! ...but for bit reproducibility with previous cycle we do it
-      ! like this:
-      PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) &
-           & - XCOSDEL*COS(ZSOLANG)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) &
-           & + XCOSDEL*SIN(ZSOLANG)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL)))
-    ENDDO
-
-  ELSE
-    ! Compute the average MU0 for the period of the model timestep
-
-    ! First compute the sine and cosine of the times of the start and
-    ! end of the model timestep
-    ZHALFTIMESTEP   = XTSTEP * RPI / RDAY
-    ZSOLARTIMESTART = ZSOLANG - ZHALFTIMESTEP
-    ZSOLARTIMEEND   = ZSOLANG + ZHALFTIMESTEP
-
-    ! Compute tangent of solar declination, with check in case someone
-    ! simulates a planet completely tipped over
-    ZTANDEC = XSINDEL / MAX(XCOSDEL, 1.0E-12)
-
-    DO JCOL = KSTART,KEND
-      ! Sunrise equation: cos(hour angle at sunset) =
-      ! -tan(declination)*tan(latitude)
-      ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) &
-           &              / MAX(ZCOSLAT(JCOL), 1.0E-12)
-      IF (ZCOSHOURANGLESUNSET > 1.0) THEN
-        ! Perpetual darkness
-        PMU0(JCOL) = 0.0_JPRB
-      ELSE
-        ! Compute hour angle at start and end of time interval,
-        ! ensuring that the hour angle of the centre of the time
-        ! window is in the range -PI to +PI (equivalent to ensuring
-        ! that local solar time = solar time + longitude is in the
-        ! range 0 to 2PI)
-        IF (ZSOLANG + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN
-          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI
-          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - RPI 
-        ELSE
-          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI
-          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - 3.0_JPRB*RPI
-        ENDIF
-
-        IF (ZCOSHOURANGLESUNSET >= -1.0) THEN
-          ! Not perpetual daylight or perpetual darkness, so we need
-          ! to check for sunrise or sunset lying within the time
-          ! interval
-          ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET)
-          IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET &
-               & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN
-            ! The time interval is either completely before sunrise or
-            ! completely after sunset
-            PMU0(JCOL) = 0.0_JPRB
-            CYCLE
-          ENDIF
-
-          ! Bound the start and end hour angles by sunrise and sunset
-          ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, &
-               &                MIN(ZHOURANGLESTART, ZHOURANGLESUNSET))
-          ZHOURANGLEEND   = MAX(-ZHOURANGLESUNSET, &
-               &                MIN(ZHOURANGLEEND,   ZHOURANGLESUNSET))
-        ENDIF
-
-        ! Compute average MU0 in the model timestep, although the
-        ! numerator considers only the time from ZHourAngleStart to
-        ! ZHourAngleEnd that the sun is above the horizon
-        PMU0(JCOL) = (ZSINDECSINLAT(JCOL) * (ZHOURANGLEEND-ZHOURANGLESTART)   &
-           & + ZCOSDECCOSLAT(JCOL)*(SIN(ZHOURANGLEEND)-SIN(ZHOURANGLESTART))) &
-           & / (2.0_JPRB * ZHALFTIMESTEP)
-
-        ! This shouldn't ever result in negative values, but just in
-        ! case
-        IF (PMU0(JCOL) < 0.0_JPRB) THEN
-          PMU0(JCOL) = 0.0_JPRB
-        ENDIF
-
-      ENDIF
-    ENDDO
-  ENDIF
-
-ENDIF
-
-
-!     ------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('COS_SZA',1,ZHOOK_HANDLE)
-END SUBROUTINE COS_SZA
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/ice_effective_radius.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/ice_effective_radius.F90
deleted file mode 100644
index 8d078cca1..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/ice_effective_radius.F90
+++ /dev/null
@@ -1,196 +0,0 @@
-!      #################################
-       MODULE MODI_ICE_EFFECTIVE_RADIUS
-!      #################################
-INTERFACE
-
-SUBROUTINE ICE_EFFECTIVE_RADIUS &
-     & (KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
-     &  PRE_UM)
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
-
-! *** Variables on model levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
-REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)  ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)       ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)      ! (kg/kg)
-
-! *** Single level variable
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON) ! Sine of latitude
-
-! OUTPUT ARGUMENT
-! Effective radius
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_UM(KLON,KLEV) ! (microns)
-
-END SUBROUTINE ICE_EFFECTIVE_RADIUS
-END INTERFACE
-END MODULE MODI_ICE_EFFECTIVE_RADIUS
-
-SUBROUTINE ICE_EFFECTIVE_RADIUS &
-     & (KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
-     &  PRE_UM)
-     
-! ICE_EFFECTIVE_RADIUS
-!
-! PURPOSE
-! -------
-!   Calculate effective radius of ice clouds
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF (using code extracted from radlswr.F90)
-!   Original: 2016-02-24
-!
-! MODIFICATIONS
-! -------------
-!
-!
-! -------------------------------------------------------------------
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-! USE YOERAD   , ONLY : YRERAD
-USE YOM_YGFL , ONLY : YGFL
-! USE YOECLDP  , ONLY : YRECLDP
-USE YOERDU   , ONLY : REPLOG, REPSCW
-USE YOMLUN   , ONLY : NULERR
-USE YOMCST   , ONLY : RD, RPI, RTT
-USE MODD_PARAM_ECRAD_n , ONLY : NRADIP
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
-
-! *** Variables on model levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
-REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)  ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)       ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)      ! (kg/kg)
-
-! *** Single level variable
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON) ! Sine of latitude
-
-! OUTPUT ARGUMENT
-! Effective radius
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_UM(KLON,KLEV) ! (microns)
-
-! LOCAL VARIABLES
-
-REAL(KIND=JPRB) :: ZIWC_INCLOUD_GM3 ! In-cloud ice+snow water content in g m-3
-REAL(KIND=JPRB) :: ZAIR_DENSITY_GM3 ! Air density in g m-3
-
-REAL(KIND=JPRB) :: ZTEMPERATURE_C   ! Temperature, degrees Celcius
-REAL(KIND=JPRB) :: ZTEMP_FACTOR     ! Temperature, Kelvin minus 83.15
-REAL(KIND=JPRB) :: ZAIWC, ZBIWC     ! Factors in empirical relationship
-REAL(KIND=JPRB) :: ZDEFAULT_RE_UM   ! Default effective radius in microns 
-REAL(KIND=JPRB) :: ZDIAMETER_UM     ! Effective diameter in microns
-
-! Min effective diameter in microns; may vary with latitude
-REAL(KIND=JPRB) :: ZMIN_DIAMETER_UM(KLON)
-
-INTEGER :: JL, JK
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-! -------------------------------------------------------------------
-
-#include "abor1.intfb.h"
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('ICE_EFFECTIVE_RADIUS',0,ZHOOK_HANDLE)
-
-! -------------------------------------------------------------------
-SELECT CASE(NRADIP)
-CASE(0)
-  ! Ice effective radius fixed at 40 microns
-  PRE_UM(KIDIA:KFDIA,:) = 40.0_JPRB  
-
-CASE(1,2)
-  ! Ice effective radius from Liou and Ou (1994)
-  DO JK = 1,KLEV
-    DO JL = KIDIA,KFDIA
-      ! Convert Kelvin to Celcius, preventing positive numbers
-      ZTEMPERATURE_C = MIN(PTEMPERATURE(JL,JK) - RTT, -0.1)
-      ! Liou and Ou's empirical formula
-      PRE_UM(JL,JK) = 326.3_JPRB + ZTEMPERATURE_C * (12.42_JPRB &
-           &  + ZTEMPERATURE_C * (0.197_JPRB + ZTEMPERATURE_C * 0.0012_JPRB))
-      IF (NRADIP == 1) THEN
-        ! Original Liou and Ou (1994) bounds of 40-130 microns
-        PRE_UM(JL,JK) = MAX(PRE_UM(JL,JK), 40.0_JPRB)
-        PRE_UM(JL,JK) = MIN(PRE_UM(JL,JK),130.0_JPRB)
-      ELSE
-        ! Formulation following Jakob, Klein modifications to ice
-        ! content
-        PRE_UM(JL,JK) = MAX(PRE_UM(JL,JK), 30.0_JPRB)
-        PRE_UM(JL,JK) = MIN(PRE_UM(JL,JK), 60.0_JPRB)
-      ENDIF
-    ENDDO
-  ENDDO
-
-CASE(3)
-  ! Ice effective radius = f(T,IWC) from Sun and Rikus (1999), revised
-  ! by Sun (2001)
-
-  ! Default effective radius is computed from an effective diameter of
-  ! 80 microns; note that multiplying by re2de actually converts from
-  ! effective diameter to effective radius.
-  ZDEFAULT_RE_UM = 80.0_JPRB * XRE2DE
-
-  ! Minimum effective diameter may vary with latitude
-  IF (NMINICE == 0) THEN
-    ! Constant effective diameter
-    ZMIN_DIAMETER_UM(KIDIA:KFDIA) = XRMINICE
-  ELSE
-    ! Ice effective radius varies with latitude, smaller at poles
-    DO JL = KIDIA,KFDIA
-      ZMIN_DIAMETER_UM(JL) = 20.0_JPRB + (XRMINICE - 20.0_JPRB) &
-           &                          * COS(ASIN(PGEMU(JL)))
-    ENDDO
-  ENDIF
-
-  DO JK = 1,KLEV
-    DO JL = KIDIA,KFDIA
-      IF (PCLOUD_FRAC(JL,JK) > 0.001_JPRB &
-           &  .AND. (PQ_ICE(JL,JK)+PQ_SNOW(JL,JK)) > 0.0_JPRB) THEN
-        ZAIR_DENSITY_GM3 = 1000.0_JPRB * PPRESSURE(JL,JK) / (RD*PTEMPERATURE(JL,JK))
-        ZIWC_INCLOUD_GM3 = ZAIR_DENSITY_GM3 * (PQ_ICE(JL,JK) + PQ_SNOW(JL,JK)) &
-             &           / PCLOUD_FRAC(JL,JK)
-        ZTEMPERATURE_C = PTEMPERATURE(JL,JK) - RTT
-        ! Sun, 2001 (corrected from Sun & Rikus, 1999)
-        ZAIWC = 45.8966_JPRB * ZIWC_INCLOUD_GM3**0.2214_JPRB
-        ZBIWC = 0.7957_JPRB  * ZIWC_INCLOUD_GM3**0.2535_JPRB
-        ZDIAMETER_UM = (1.2351_JPRB + 0.0105_JPRB * ZTEMPERATURE_C) &
-             & * (ZAIWC + ZBIWC*(PTEMPERATURE(JL,JK) - 83.15_JPRB))
-        ZDIAMETER_UM = MIN ( MAX( ZDIAMETER_UM, ZMIN_DIAMETER_UM(JL)), 155.0_JPRB)
-        PRE_UM(JL,JK) = ZDIAMETER_UM * XRE2DE
-      ELSE
-        PRE_UM(JL,JK) = ZDEFAULT_RE_UM
-      ENDIF
-    ENDDO
-  ENDDO
-  
-CASE DEFAULT
-  WRITE(NULERR,'(A,I0,A)') 'ICE EFFECTIVE RADIUS OPTION NRADLP=',NRADIP,' NOT AVAILABLE'
-  CALL ABOR1('ERROR IN ICE_EFFECTIVE_RADIUS')
-
-END SELECT
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('ICE_EFFECTIVE_RADIUS',1,ZHOOK_HANDLE)
-  
-END SUBROUTINE ICE_EFFECTIVE_RADIUS
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/liquid_effective_radius.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/liquid_effective_radius.F90
deleted file mode 100644
index 7a1348ad9..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/liquid_effective_radius.F90
+++ /dev/null
@@ -1,300 +0,0 @@
-MODULE MODI_LIQUID_EFFECTIVE_RADIUS
-
-INTERFACE
-
-SUBROUTINE LIQUID_EFFECTIVE_RADIUS &
-     & (KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQ, PQ_RAIN, &
-     &  PLAND_FRAC, PCCN_LAND, PCCN_SEA, &
-     &  PRE_UM)
-
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-
-
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
-
-! *** Variables on model levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
-REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQ(KLON,KLEV)       ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)      ! (kg/kg)
-
-! *** Single-level variables 
-REAL(KIND=JPRB),   INTENT(IN) :: PLAND_FRAC(KLON)        ! 1=land, 0=sea
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
-
-! OUTPUT ARGUMENT
-
-! Effective radius
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_UM(KLON,KLEV) ! (microns)
-
-END SUBROUTINE LIQUID_EFFECTIVE_RADIUS
-END INTERFACE
-END MODULE MODI_LIQUID_EFFECTIVE_RADIUS
-
-SUBROUTINE LIQUID_EFFECTIVE_RADIUS &
-     & (KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQ, PQ_RAIN, &
-     &  PLAND_FRAC, PCCN_LAND, PCCN_SEA, &
-     &  PRE_UM)
-
-! LIQUID_EFFECTIVE_RADIUS
-!
-! PURPOSE
-! -------
-!   Calculate effective radius of liquid clouds
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF (using code extracted from radlswr.F90)
-!   Original: 2015-09-24
-!
-! MODIFICATIONS
-! -------------
-!
-!
-! -------------------------------------------------------------------
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-!USE YOERAD   , ONLY : YRERAD
-USE MODD_PARAM_ECRAD_n , ONLY : NRADLP, NAERMACC, NMCVAR, XCCNSEA, XCCNLND, NAERCLD, &
-                            &   NACTAERO, LCCNO, LCCNL
-
-USE YOM_YGFL , ONLY : YGFL
-! USE YOECLDP  , ONLY : YRECLDP
-USE YOERDU   , ONLY : REPLOG, REPSCW
-USE YOMLUN   , ONLY : NULERR
-USE YOMCST   , ONLY : RD, RPI
-
-! -------------------------------------------------------------------
-
-IMPLICIT NONE
-
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
-
-! *** Variables on model levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
-REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQ(KLON,KLEV)       ! (kg/kg)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)      ! (kg/kg)
-
-! *** Single-level variables 
-REAL(KIND=JPRB),   INTENT(IN) :: PLAND_FRAC(KLON)        ! 1=land, 0=sea
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
-
-! OUTPUT ARGUMENT
-
-! Effective radius
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_UM(KLON,KLEV) ! (microns)
-
-! PARAMETERS
-
-! Minimum and maximum effective radius, in microns
-REAL(KIND=JPRB), PARAMETER :: PP_MIN_RE_UM =  4.0_JPRB
-REAL(KIND=JPRB), PARAMETER :: PP_MAX_RE_UM = 30.0_JPRB
-
-! LOCAL VARIABLES
-INTEGER :: IRADLP ! ID of effective radius scheme to use
-INTEGER :: NACTIVE_AEROSOL ! Number of active aerosol
-REAL(KIND=JPRB) :: ZCCN    ! CCN concentration (units?)
-
-REAL(KIND=JPRB) :: ZSPECTRAL_DISPERSION
-REAL(KIND=JPRB) :: ZNTOT_CM3 ! Number conc in cm-3
-REAL(KIND=JPRB) :: ZRE_CUBED
-REAL(KIND=JPRB) :: ZLWC_GM3, ZRWC_GM3 ! In-cloud liquid, rain content in g m-3
-REAL(KIND=JPRB) :: ZAIR_DENSITY_GM3   ! Air density in g m-3
-REAL(KIND=JPRB) :: ZRAIN_RATIO        ! Ratio of rain to liquid water content
-REAL(KIND=JPRB) :: ZWOOD_FACTOR, ZRATIO
-REAL(KIND=JPRB) :: ZNUM, ZDEN, ZNTOT,ZD
-
-INTEGER :: JL, JK
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-! -------------------------------------------------------------------
-
-#include "abor1.intfb.h"
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('LIQUID_EFFECTIVE_RADIUS',0,ZHOOK_HANDLE)
-
-! -------------------------------------------------------------------
-
-! Reproduce logic from RADLSWR
-NACTIVE_AEROSOL = NACTAERO
-IF (NACTAERO == 0 .AND. NAERMACC == 1) NACTIVE_AEROSOL = NMCVAR
-IRADLP = NRADLP
-IF (NACTIVE_AEROSOL >= 12 .AND. NAERCLD > 0 ) IRADLP=3 
-
-SELECT CASE(IRADLP)
-CASE(0)
-  ! Very old parameterization as a function of pressure, used in ERA-15
-  PRE_UM(KIDIA:KFDIA,:) = 10.0_JPRB &
-       &  + (100000.0_JPRB-PPRESSURE(KIDIA:KFDIA,:))*3.5_JPRB
-  
-CASE(1)
-  ! Simple distinction between land (10um) and ocean (13um) by Zhang
-  ! and Rossow
-  DO JL = KIDIA,KFDIA
-    IF (PLAND_FRAC(JL) < 0.5_JPRB) THEN
-      PRE_UM(JL,:) = 13.0_JPRB
-    ELSE
-      PRE_UM(JL,:) = 10.0_JPRB
-    ENDIF
-  ENDDO
-  
-CASE(2)
-  ! Martin et al. (JAS 1994)
-  ! Old ECMWF version
-  DO JL = KIDIA,KFDIA 
-    IF (PLAND_FRAC(JL) < 0.5_JPRB) THEN
-        ZCCN=150.
-        ZD=0.33
-        ZNTOT=-1.15E-03*ZCCN*ZCCN+0.963*ZCCN+5.30
-    ELSE
-        ZCCN=900.
-        ZD=0.43
-        ZNTOT=-2.10E-04*ZCCN*ZCCN+0.568*ZCCN-27.9
-    ENDIF
-    
-   ! Only consider cloudy regions
-    DO JK = 1,KLEV
-        IF (PCLOUD_FRAC(JL,JK) >= 0.001_JPRB &
-            &  .AND. (PQ_LIQ(JL,JK)+PQ_RAIN(JL,JK)) > 0.0_JPRB) THEN
-
-            ! Compute liquid and rain water contents
-            ZAIR_DENSITY_GM3 = 1000.0_JPRB * PPRESSURE(JL,JK) &
-                    &           / (RD*PTEMPERATURE(JL,JK))
-            ! In-cloud mean water contents found by dividing by cloud
-            ! fraction
-            ZLWC_GM3 = ZAIR_DENSITY_GM3 * PQ_LIQ(JL,JK)  / PCLOUD_FRAC(JL,JK)
-            ZRWC_GM3 = ZAIR_DENSITY_GM3 * PQ_RAIN(JL,JK) / PCLOUD_FRAC(JL,JK)
-    
-            ZNUM=3.*(ZLWC_GM3+ZRWC_GM3)*(1.+3.*ZD*ZD)**2 ! conversion mixing ratio into LWC
-            ZDEN=4.*RPI*ZNTOT*(1.+ZD*ZD)**3
-            PRE_UM(JL,JK) =100.*(ZNUM/ZDEN)**0.333
-            PRE_UM(JL,JK)=MAX(PRE_UM(JL,JK), 4.)
-            PRE_UM(JL,JK)=MIN(PRE_UM(JL,JK),16.)
-        ENDIF
-    ENDDO  
-    
-  ENDDO 
-! print*, "PRE_UM"
-! print*, PRE_UM 
-! pause
-! 
-! print*, "PQ_LIQ+PQ_RAIN"
-! print*, PQ_LIQ+PQ_RAIN
-! pause 
-
- 
-CASE(3) 
- ! Improved version in ECRAD
-  DO JL = KIDIA,KFDIA
-    ! First compute the cloud droplet concentration
-    IF (PLAND_FRAC(JL) < 0.5_JPRB) THEN
-      ! Sea case
-      IF (LCCNO) THEN
-        ZCCN = PCCN_SEA(JL)
-      ELSE
-        ZCCN = XCCNSEA
-      ENDIF
-      ZSPECTRAL_DISPERSION = 0.77_JPRB
-      ! Cloud droplet concentration in cm-3 (activated CCN) over
-      ! ocean
-      ZNTOT_CM3 = -1.15E-03_JPRB*ZCCN*ZCCN + 0.963_JPRB*ZCCN + 5.30_JPRB
-    ELSE
-      ! Land case
-      IF (LCCNL) THEN 
-        ZCCN=PCCN_LAND(JL)
-      ELSE  
-        ZCCN=XCCNLND
-      ENDIF
-      ZSPECTRAL_DISPERSION = 0.69_JPRB
-      ! Cloud droplet concentration in cm-3 (activated CCN) over
-      ! land
-      ZNTOT_CM3 = -2.10E-04_JPRB*ZCCN*ZCCN + 0.568_JPRB*ZCCN - 27.9_JPRB
-    ENDIF
-    
-    ZRATIO = (0.222_JPRB/ZSPECTRAL_DISPERSION)**0.333_JPRB
-    
-    DO JK = 1,KLEV
-
-      ! Only consider cloudy regions
-      IF (PCLOUD_FRAC(JL,JK) >= 0.001_JPRB &
-           &  .AND. (PQ_LIQ(JL,JK)+PQ_RAIN(JL,JK)) > 0.0_JPRB) THEN
-
-        ! Compute liquid and rain water contents
-        ZAIR_DENSITY_GM3 = 1000.0_JPRB * PPRESSURE(JL,JK) &
-             &           / (RD*PTEMPERATURE(JL,JK))
-        ! In-cloud mean water contents found by dividing by cloud
-        ! fraction
-        ZLWC_GM3 = ZAIR_DENSITY_GM3 * PQ_LIQ(JL,JK)  / PCLOUD_FRAC(JL,JK)
-        ZRWC_GM3 = ZAIR_DENSITY_GM3 * PQ_RAIN(JL,JK) / PCLOUD_FRAC(JL,JK)
-      
-        ! Wood's (2000, eq. 19) adjustment to Martin et al's
-        ! parameterization
-        IF (ZLWC_GM3 > REPSCW) THEN
-          ZRAIN_RATIO = ZRWC_GM3 / ZLWC_GM3
-          ZWOOD_FACTOR = ((1.0_JPRB + ZRAIN_RATIO)**0.666_JPRB) &
-               &     / (1.0_JPRB + 0.2_JPRB * ZRATIO*ZRAIN_RATIO)
-        ELSE
-          ZWOOD_FACTOR = 1.0_JPRB
-        ENDIF
-      
-        ! g m-3 and cm-3 units cancel out with density of water
-        ! 10^6/(1000*1000); need a factor of 10^6 to convert to
-        ! microns and cubed root is factor of 100 which appears in
-        ! equation below
-        ZRE_CUBED = (3.0_JPRB * (ZLWC_GM3 + ZRWC_GM3)) &
-             &    / (4.0_JPRB*RPI*ZNTOT_CM3*ZSPECTRAL_DISPERSION)
-        IF (ZRE_CUBED > REPLOG) THEN
-          PRE_UM(JL,JK) = ZWOOD_FACTOR*100.0_JPRB*EXP(0.333_JPRB*LOG(ZRE_CUBED))
-          ! Make sure effective radius is bounded in range 4-30 microns
-          PRE_UM(JL,JK) = MAX(PP_MIN_RE_UM, MIN(PRE_UM(JL,JK), PP_MAX_RE_UM))
-        ELSE
-          PRE_UM(JL,JK) = PP_MIN_RE_UM
-        ENDIF
-
-      ELSE
-        ! Cloud fraction or liquid+rain water content too low to
-        ! consider this a cloud
-        PRE_UM(JL,JK) = PP_MIN_RE_UM
-
-      ENDIF
-
-    ENDDO
-    
-  ENDDO
-  
-
-  
-CASE DEFAULT
-  WRITE(NULERR,'(A,I0,A)') 'LIQUID EFFECTIVE RADIUS OPTION IRADLP=',IRADLP,' NOT AVAILABLE'
-  CALL ABOR1('ERROR IN LIQUID_EFFECTIVE_RADIUS')
-END SELECT
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('LIQUID_EFFECTIVE_RADIUS',1,ZHOOK_HANDLE)
-  
-END SUBROUTINE LIQUID_EFFECTIVE_RADIUS
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_scheme.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_scheme.F90
deleted file mode 100644
index e0cdd6121..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_scheme.F90
+++ /dev/null
@@ -1,625 +0,0 @@
-!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/ecmwf_radiation_vers2.f90,v $ $Revision: 1.3.2.4.2.2.2.1 $
-! masdev4_7 BUG1 2007/06/15 17:47:17
-!-----------------------------------------------------------------
-!      #################################
-       MODULE MODI_RADIATION_SCHEME
-!      #################################
-
-CONTAINS
-
-SUBROUTINE RADIATION_SCHEME &
-     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, &
-     &  PSOLAR_IRRADIANCE, &
-     &  PMU0, PTEMPERATURE_SKIN, PALBEDO_DIF, PALBEDO_DIR, &
-     &  PEMIS, PEMIS_WINDOW, &
-     &  PCCN_LAND, PCCN_SEA, &
-     &  PGELAM, PGEMU, PLAND_SEA_MASK, &
-     &  PPRESSURE, PTEMPERATURE, &
-     &  PPRESSURE_H, PTEMPERATURE_H, &
-     &  PQ, PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, PCCL4, PO3_DP, &
-     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_RAIN, PQ_SNOW, &
-     &  PAEROSOL_OLD, PAEROSOL, &
-     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
-     &  PFLUX_SW_SURF, PFLUX_LW_SURF, PFLUX_SW_SURF_CLEAR, PFLUX_LW_SURF_CLEAR, &
-     &  PFLUX_DIR_SURF, PFLUX_DIR_SURF_CLEAR, PFLUX_DIR_SURF_INTO_SUN, &
-     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
-     &  PFLUX_SW_DN_TOA,PFLUX_SW_UP_TOA,PFLUX_LW_UP_TOA, &
-     &  PFLUX_SW_UP_TOA_CLEAR,PFLUX_LW_UP_TOA_CLEAR, &
-     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_UP, PFLUX_LW_UP, &
-     &  PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
-     &  PRE_LIQUID_UM, PRE_ICE_UM, &
-     &  PEMIS_OUT, PLWDERIVATIVE, &
-     &  PSWDIFFUSEBAND, PSWDIRECTBAND)
-
-! RADIATION_SCHEME - Interface to modular radiation scheme
-!
-! PURPOSE
-! -------
-!   The modular radiation scheme is contained in a separate
-!   library. This routine puts the the IFS arrays into appropriate
-!   objects, computing the additional data that is required, and sends
-!   it to the radiation scheme.  It returns net fluxes and surface
-!   flux components needed by the rest of the model. 
-!
-!   Lower case is used for variables and types taken from the
-!   radiation library
-!
-! INTERFACE
-! ---------
-!    RADIATION_SCHEME is called from RADLSWR. The
-!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
-!    should have been run first.
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF
-!   Original: 2015-09-16
-!
-! MODIFICATIONS
-! -------------
-!
-! TO DO
-! -----
-!
-!-----------------------------------------------------------------------
-
-! Modules from ifs or ifsaux libraries
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-! USE YOERAD   , ONLY : YRERAD ! does not exist in MNH
-
-! MNH
-USE MODD_PARAM_ECRAD_n, ONLY : NAERMACC, NDECOLAT, XCLOUD_FRAC_STD,  &  ! replace YRERAD to get attributes
-                            &  LAPPROXLWUPDATE, LAPPROXSWUPDATE
-USE MODD_RADIATIONS_n , ONLY : NSWB_MNH,NSWB_OLD                         
-USE MODE_THERMO  , ONLY : QSAT
-USE MODD_DYN_n , ONLY : XTSTEP, NSTOP
-USE MODD_TIME , ONLY : TDTEXP
-USE MODD_TIME_n , ONLY : TDTMOD,TDTCUR
-USE MODI_ICE_EFFECTIVE_RADIUS
-USE MODI_LIQUID_EFFECTIVE_RADIUS
-USE MODI_CLOUD_OVERLAP_DECORR_LEN
-USE MODD_LUNIT_n , ONLY : TLUOUT
-! MNH
-USE YOMLUN    ,ONLY : NULOUT
-USE RADIATION_SETUP, ONLY : rad_config, &
-     &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
-     &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
-     &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
-     &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT
-     
-! USE YOMRIP0  , ONLY : NINDAT ! does not exist in MNH
-! USE YOMCT3   , ONLY : NSTEP  ! does not exist in MNH
-! USE YOMRIP   , ONLY : YRRIP  ! does not exist in MNH
-USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
-
-
-! Modules from radiation library
-USE radiation_single_level,   ONLY : single_level_type
-USE radiation_thermodynamics, ONLY : thermodynamics_type
-USE radiation_gas
-USE radiation_cloud,          ONLY : cloud_type
-USE radiation_aerosol,        ONLY : aerosol_type
-USE radiation_flux,           ONLY : flux_type
-USE radiation_interface,      ONLY : radiation, set_gas_units
-USE radiation_save,           ONLY : save_inputs
-
-IMPLICIT NONE
-
-! INPUT ARGUMENTS
-
-! *** Array dimensions and ranges
-INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
-INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
-INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL ! Number of aerosol types
-
-! *** Single-level fields
-REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
-REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
-! Diffuse and direct components of surface shortwave albedo
-REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSWB_OLD)
-REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSWB_OLD)
-! Longwave emissivity outside and inside the window region
-REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
-REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
-! Longitude (radians), sine of latitude
-REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
-REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
-! Land-sea mask
-REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON) 
-
-! *** Variables on full levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
-! *** Variables on half levels
-REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
-REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
-
-! *** Gas mass mixing ratios on full levels
-REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCO2(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCH4(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PN2O(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PNO2(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCFC11(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCFC12(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCCL4(KLON,KLEV) 
-REAL(KIND=JPRB),   INTENT(IN) :: PO3_DP(KLON,KLEV) ! (Pa*kg/kg) !
-
-! *** Cloud fraction and hydrometeor mass mixing ratios
-REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
-REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
-
-! *** Aerosol mass mixing ratios
-REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV) ! Tegen in optical thickness at 550 nm
-REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
-
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON) 
-REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON) 
-
-! OUTPUT ARGUMENTS
-
-! *** Net fluxes on half-levels (W m-2)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1) 
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1) 
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1) 
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1) 
-
-! *** Surface flux components (W m-2)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_SURF(KLON,NSWB_MNH) 
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_SURF(KLON) 
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_SURF_CLEAR(KLON,NSWB_MNH)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_SURF_CLEAR(KLON)
-! Direct component of surface flux into horizontal plane
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_SURF(KLON,NSWB_MNH)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_SURF_CLEAR(KLON,NSWB_MNH)
-! As PFLUX_DIR but into a plane perpendicular to the sun
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_SURF_INTO_SUN(KLON,NSWB_MNH)
-
-! *** Ultraviolet and photosynthetically active radiation (W m-2)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
-
-! *** Other single-level diagnostics
-! Top-of-atmosphere fluxes (W m-2)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_TOA(KLON)
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_TOA(KLON) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_TOA(KLON) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_TOA_CLEAR(KLON) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_TOA_CLEAR(KLON) ! QL
-
-! Total fluxes - QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1) ! QL
-REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1) ! QL
-
-! Cloud effective radii in microns
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_LIQUID_UM(KLON,KLEV)
-REAL(KIND=JPRB),  INTENT(OUT) :: PRE_ICE_UM(KLON,KLEV)
-
-! Diagnosed longwave surface emissivity across the whole spectrum
-REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   
-
-! Partial derivative of total-sky longwave upward flux at each level
-! with respect to upward flux at surface, used to correct heating
-! rates at gridpoints/timesteps between calls to the full radiation
-! scheme.  Note that this version uses the convention of level index
-! increasing downwards, unlike the local variable ZLwDerivative that
-! is returned from the LW radiation scheme.
-REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
-
-! Surface diffuse and direct downwelling shortwave flux in each
-! shortwave albedo band, used in RADINTG to update the surface fluxes
-! accounting for high-resolution albedo information
-REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSWB_MNH)
-REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSWB_MNH)
-
-! LOCAL VARIABLES
-TYPE(single_level_type)   :: single_level
-TYPE(thermodynamics_type) :: thermodynamics
-TYPE(gas_type)            :: gas
-TYPE(cloud_type)          :: cloud
-TYPE(aerosol_type)        :: aerosol
-TYPE(flux_type)           :: flux
-
-! Mass mixing ratio of ozone (kg/kg)
-REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)
-
-! Cloud overlap decorrelation length for cloud boundaries in km
-REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)
-
-! Ratio of cloud overlap decorrelation length for cloud water
-! inhomogeneities to that for cloud boundaries (typically 0.5)
-REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO
-
-! The surface net longwave flux if the surface was a black body, used
-! to compute the effective broadband surface emissivity
-REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
-
-! Layer mass in kg m-2
-REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
-
-! Time integers
-INTEGER :: ITIM, IDAY
-
-! Loop indices
-INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-
-
-! Import time functions for iseed calculation
-!#include "../include/fcttim.func.h"
-! #include "liquid_effective_radius.intfb.h"
-! #include "ice_effective_radius.intfb.h"
-! #include "cloud_overlap_decorr_len.intfb.h"
-
-! PRINT*,"RADIATION_SCHEME"
-
-IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
-
-! Redirect OUTPUT to MNH OUTPUT_LISTING
-NULOUT = TLUOUT%NLU
-
-! Allocate memory in radiation objects
-CALL single_level%allocate(KLON, NSWB_OLD, 2, &                ! ncol, nalbedobands, nemisbands
-     &                     use_sw_albedo_direct=.TRUE.)        ! mu0, tskin, albedo, emiss, seed for McICA 
-CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)   ! p, t, qsat; hlevel, nlev+1
-CALL gas%allocate(KLON, KLEV)                                  ! mixing ratios
-CALL cloud%allocate(KLON, KLEV)                                ! ql, qi, radl, radi, clfr, cw-sig, covlp
-IF (NAERMACC > 0) THEN
-  CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL) ! MACC climatology
-ELSE
-  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology  ! mixing ratios per class
-ENDIF
-CALL flux%allocate(rad_config, 1, KLON, KLEV)                  ! output fluxes depend on config; cloud cover
-                                                               ! do_lw and do_sw set to .true. ?
-                                                               
-! Set thermodynamic profiles: simply copy over the half-level
-! pressure and temperature
-thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
-thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
-
-! IFS currently sets the half-level temperature at the surface to be
-! equal to the skin temperature. The radiation scheme takes as input
-! only the half-level temperatures and assumes the Planck function to
-! vary linearly in optical depth between half levels. In the lowest
-! atmospheric layer, where the atmospheric temperature can be much
-! cooler than the skin temperature, this can lead to significant
-! differences between the effective temperature of this lowest layer
-! and the true value in the model.
-!
-! We may approximate the temperature profile in the lowest model level
-! as piecewise linear between the top of the layer T[k-1/2], the
-! centre of the layer T[k] and the base of the layer Tskin.  The mean
-! temperature of the layer is then 0.25*T[k-1/2] + 0.5*T[k] +
-! 0.25*Tskin, which can be achieved by setting the atmospheric
-! temperature at the half-level corresponding to the surface as
-! follows:
-thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
-     &  = PTEMPERATURE(KIDIA:KFDIA,KLEV) &
-     &  + 0.5_JPRB * (PTEMPERATURE_H(KIDIA:KFDIA,KLEV+1) &
-     &               -PTEMPERATURE_H(KIDIA:KFDIA,KLEV))
-
-! Alternatively we respect the model's atmospheric temperature in the
-! lowest model level by setting the temperature at the lowest
-! half-level such that the mean temperature of the layer is correct:
-!thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
-!     &  = 2.0_JPRB * PTEMPERATURE(KIDIA:KFDIA,KLEV) &
-!     &             - PTEMPERATURE_H(KIDIA:KFDIA,KLEV)
-
-! Compute saturation specific humidity, used to hydrate aerosols. The
-! "2" for the last argument indicates that the routine is not being
-! called from within the convection scheme.
-!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
-!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2)  
-     
-!MNH     
-thermodynamics%h2o_sat_liq(:,:) = QSAT(REAL(PPRESSURE), REAL(PTEMPERATURE))
-thermodynamics%h2o_sat_liq(:,:) = thermodynamics%h2o_sat_liq(:,:) &
-                                & / (1.+thermodynamics%h2o_sat_liq(:,:)) ! mixing ratio => spec humid
-!MNH
-
-! Alternative approximate version using temperature and pressure from
-! the thermodynamics structure
-!CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
-
-! Set single-level fileds
-single_level%solar_irradiance              = PSOLAR_IRRADIANCE
-single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
-single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
-single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
-single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
-
-! pause
-
-! Longwave emissivity is in two bands
-single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
-single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)
-
-! Create the relevant seed from date and time get the starting day
-! and number of minutes since start
-! IDAY = NDD(NINDAT)  ! NINDAT is AAAAMMDD initial date NDD extract DD as INTEGER
-IDAY = TDTEXP%DATE%NDAY ! MNH
-ITIM = NINT((TDTMOD%xtime-TDTCUR%xtime) / 60.0_JPRB) ! YRRIP contains timestep infos ; number of minutes since beginning
-DO JLON = KIDIA, KFDIA
-  ! This method gives a unique value for roughly every 1-km square
-  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
-  ! latitude in degrees, which we multiply by 100 to give a unique
-  ! value for roughly every km. PGELAM*60*100 gives a unique number
-  ! for roughly every km of longitude around the equator, which we
-  ! multiply by 180*100 so there is no overlap with the latitude
-  ! values.  The result can be contained in a 32-byte integer (but
-  ! since random numbers are generated with the help of integer
-  ! overflow, it should not matter if the number did overflow).
-  single_level%iseed(JLON) = ITIM + IDAY & 
-       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
-       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
-ENDDO
-
-! Set cloud fields
-cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
-cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
-cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
-
-! Compute effective radii and convert to metres - Could take it from LIMA
-CALL LIQUID_EFFECTIVE_RADIUS(KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQUID, PQ_RAIN, &
-     &  PLAND_SEA_MASK, PCCN_LAND, PCCN_SEA, &
-     &  PRE_LIQUID_UM) 
-    
-cloud%re_liq(KIDIA:KFDIA,:) = PRE_LIQUID_UM(KIDIA:KFDIA,:) * 1.0e-6_JPRB ! conversion into meters
-
-CALL ICE_EFFECTIVE_RADIUS(KIDIA, KFDIA, KLON, KLEV, &
-     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
-     &  PRE_ICE_UM)
-     
-cloud%re_ice(KIDIA:KFDIA,:) = PRE_ICE_UM(KIDIA:KFDIA,:) * 1.0e-6_JPRB
-
-! Get the cloud overlap decorrelation length (for cloud boundaries),
-! in km, according to the parameterization specified by NDECOLAT,
-! and insert into the "cloud" object. Also get the ratio of
-! decorrelation lengths for cloud water content inhomogeneities and
-! cloud boundaries, and set it in the "rad_config" object.
-CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, NDECOLAT, &
-     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
-     
-rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO
-
-DO JLON = KIDIA,KFDIA
-  CALL cloud%set_overlap_param(thermodynamics, &
-       &                       ZDECORR_LEN_KM(JLON)*1000.0_JPRB, &
-       &                       istartcol=JLON, iendcol=JLON)
-ENDDO
-
-! Cloud water content fractional standard deviation is configurable
-! from namelist NAERAD but must be globally constant. Before it was
-! hard coded at 1.0.
-CALL cloud%create_fractional_std(KLON, KLEV, XCLOUD_FRAC_STD)
-
-! By default mid and high cloud effective size is 10 km
-CALL cloud%create_inv_cloud_effective_size(KLON,KLEV,1.0_JPRB/10000.0_JPRB)
-! But for boundary clouds (eta > 0.8) we set it to 1 km
-DO JLEV = 1,KLEV
-  DO JLON = KIDIA,KFDIA
-    IF (PPRESSURE(JLON,JLEV) > 0.8_JPRB * PPRESSURE_H(JLON,KLEV+1)) THEN
-      cloud%inv_cloud_effective_size(JLON,JLEV) = 1.0e-3_JPRB
-    ENDIF
-  ENDDO
-ENDDO
-
-
-! Compute the dry mass of each layer neglecting humidity effects, in
-! kg m-2, needed to scale some of the aerosol inputs
-CALL thermodynamics%get_layer_mass(ZLAYER_MASS)
-
-! Copy over aerosol mass mixing ratio
-IF (NAERMACC > 0) THEN
-
-  ! MACC aerosol climatology - this is already in mass mixing ratio
-  ! units with the required array orientation so we can copy it over
-  ! directly
-  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)
-
-  ! Add the tropospheric and stratospheric backgrounds contained in the
-  ! old Tegen arrays - this is very ugly!
-  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
-    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
-         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
-         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
-         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
-  ENDIF
-  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
-    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
-         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
-         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
-         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
-  ENDIF
-
-ELSE
-
-  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
-  ! 550-nm optical depth in each layer. The optics data file
-  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
-  ! coefficient, but a scaling factor that the 550-nm optical depth
-  ! should be multiplied by to obtain the optical depth in each
-  ! spectral band.  Therefore, in order for the units to work out, we
-  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
-  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
-  ! need to permute the array.
-  DO JLEV = 1,KLEV
-    DO JAER = 1,6
-      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
-         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
-         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
-    ENDDO
-  ENDDO
-
-ENDIF
-
-! Convert ozone Pa*kg/kg to kg/kg
-DO JLEV = 1,KLEV
-  DO JLON = KIDIA,KFDIA
-    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
-         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
-  ENDDO
-ENDDO
-
-! Insert gas mixing ratios
-CALL gas%put(IH2O,    IMassMixingRatio, PQ)
-CALL gas%put(ICO2,    IMassMixingRatio, PCO2)
-CALL gas%put(ICH4,    IMassMixingRatio, PCH4)
-CALL gas%put(IN2O,    IMassMixingRatio, PN2O)
-CALL gas%put(ICFC11,  IMassMixingRatio, PCFC11)
-CALL gas%put(ICFC12,  IMassMixingRatio, PCFC12)
-CALL gas%put(IHCFC22, IMassMixingRatio, PHCFC22)
-CALL gas%put(ICCL4,   IMassMixingRatio, PCCL4)
-CALL gas%put(IO3,     IMassMixingRatio, ZO3)
-CALL gas%put_well_mixed(IO2, IVolumeMixingRatio, 0.20944_JPRB) ! fixed volume mixing ratio for oxygen
-
-! Ensure the units of the gas mixing ratios are what is required by
-! the gas absorption model
-call set_gas_units(rad_config, gas)
-! 
-! Call radiation scheme
-CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
-     &  single_level, thermodynamics, gas, cloud, aerosol, flux)
-
-! Compute required output fluxes
-! First the net fluxes
-PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
-PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
-PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
-     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
-PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
-     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)
-     
-    
-! print*,"flux%sw_dn_surf_band(:,:)",flux%sw_dn_surf_band(:,:)
-!Now the surface fluxes
-PFLUX_SW_SURF      (KIDIA:KFDIA,:) = TRANSPOSE(flux%sw_dn_surf_band        (:,KIDIA:KFDIA)) ! modif Quentin
-PFLUX_LW_SURF      (KIDIA:KFDIA) = flux%lw_dn             (KIDIA:KFDIA,KLEV+1)
-PFLUX_SW_SURF_CLEAR(KIDIA:KFDIA,:) = TRANSPOSE(flux%sw_dn_surf_clear_band  (:,KIDIA:KFDIA)) ! modif Quentin
-PFLUX_LW_SURF_CLEAR(KIDIA:KFDIA) = flux%lw_dn_clear       (KIDIA:KFDIA,KLEV+1)
-PFLUX_DIR_SURF     (KIDIA:KFDIA,:) = TRANSPOSE(flux%sw_dn_direct_surf_band      (:,KIDIA:KFDIA)) ! modif Quentin
-PFLUX_DIR_SURF_CLEAR  (KIDIA:KFDIA,:) = TRANSPOSE(flux%sw_dn_direct_surf_clear_band (:,KIDIA:KFDIA)) ! modif Quentin
-PFLUX_DIR_SURF_INTO_SUN(KIDIA:KFDIA,:) = 0.0_JPRB
-DO JBAND = 1,NSWB_MNH
-    WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))    
-       PFLUX_DIR_SURF_INTO_SUN(KIDIA:KFDIA,JBAND) = PFLUX_DIR_SURF(KIDIA:KFDIA,JBAND) / PMU0(KIDIA:KFDIA)     
-    END WHERE
-END DO 
-
-! Top-of-atmosphere downwelling flux
-PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
-
-! Top-of-atmosphere upwelling fluxes - Q.L.
-PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
-PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
-PFLUX_SW_UP_TOA_CLEAR(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,1)
-PFLUX_LW_UP_TOA_CLEAR(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,1)
-
-
-
-! Total fluxes - QL
-! print*,"flux%sw_dn(KIDIA:KFDIA,:)",flux%sw_dn(KIDIA:KFDIA,:)
-
-PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
-PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
-PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
-PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
-PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
-PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
-PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
-PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
-
-
-! print*, "flux%sw_dn_band(:,1,1:10)"
-! print*,shape(flux%sw_dn_band)
-! print*, flux%sw_dn_band(1,1,1:5)
-! print*, flux%sw_dn_band(2,1,1:5)
-! print*, flux%sw_dn_band(3,1,1:5)
-! print*, flux%sw_dn_band(4,1,1:5)
-! print*, flux%sw_dn_band(5,1,1:5)
-! pause
-
-! Compute UV fluxes as weighted sum of appropriate shortwave bands
-PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
-DO JBAND = 1,NWEIGHT_UV
-  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
-       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
-ENDDO
-
-! Compute photosynthetically active radiation similarly
-PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
-PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
-DO JBAND = 1,NWEIGHT_PAR
-  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
-       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
-  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
-       &  + WEIGHT_PAR(JBAND) &
-       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
-ENDDO
-
-! Compute effective broadband emissivity
-ZBLACK_BODY_NET_LW = PFLUX_LW_SURF(KIDIA:KFDIA) &
-     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
-PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
-WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5) 
-  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
-END WHERE
-! Copy longwave derivatives
-IF (LAPPROXLWUPDATE) THEN
-  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
-END IF
-
-! Store the shortwave downwelling fluxes in each albedo band
-
-IF (LAPPROXSWUPDATE) THEN
-  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
-  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
-  DO JBAND = 1,rad_config%n_bands_sw
-    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
-    DO JLON = KIDIA,KFDIA
-      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
-           &  + flux%sw_dn_surf_band(JBAND,JLON) &
-           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
-      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
-           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
-    ENDDO
-  ENDDO
-ENDIF
-
-CALL single_level%deallocate
-CALL thermodynamics%deallocate
-CALL gas%deallocate
-CALL cloud%deallocate
-CALL aerosol%deallocate
-CALL flux%deallocate
-
-IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
-
-END SUBROUTINE RADIATION_SCHEME
-
-END MODULE MODI_RADIATION_SCHEME
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_setup.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_setup.F90
deleted file mode 100644
index 045d7a395..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifs/radiation_setup.F90
+++ /dev/null
@@ -1,368 +0,0 @@
-MODULE RADIATION_SETUP
-
-! RADIATION_SETUP - Setting up modular radiation scheme
-!
-! PURPOSE
-! -------
-!   The modular radiation scheme is contained in a separate
-!   library. SETUP_RADIATION_SCHEME in this module sets up a small
-!   number of global variables needed to store the information for it.
-!
-!   Lower case is used for variables and types taken from the
-!   radiation library
-!
-! INTERFACE
-! ---------
-!   SETUP_RADIATION_SCHEME is called from SUECRAD.  The radiation
-!   scheme is actually run using the RADIATION_SCHEME routine (not in
-!   this module).
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF
-!   Original: 2015-09-16
-!
-! MODIFICATIONS
-! -------------
-!
-!-----------------------------------------------------------------------
-
-  USE PARKIND1,         ONLY : JPRB
-  USE radiation_config, ONLY : config_type, &
-       &                       ISolverMcICA, ISolverSpartacus, &
-       &                       ILiquidModelSlingo, ILiquidModelSOCRATES, &
-       &                       IIceModelFu, IIceModelBaran, &
-       &                       IOverlapExponential,IOverlapExponentialRandom
-
-  ! Store configuration information for the radiation scheme in a
-  ! global variable
-  !type(config_type) :: rad_config
-  USE MODD_PARAM_ECRAD_n , ONLY : rad_config
-
- IMPLICIT NONE
-
-  ! Ultraviolet weightings
-  INTEGER         :: NWEIGHT_UV
-  INTEGER         :: IBAND_UV(100)
-  REAL(KIND=JPRB) :: WEIGHT_UV(100)
-  ! Photosynthetically active radiation weightings
-  INTEGER         :: NWEIGHT_PAR
-  INTEGER         :: IBAND_PAR(100)
-  REAL(KIND=JPRB) :: WEIGHT_PAR(100)
-
-  ! Background aerosol is specified in an ugly way: using the old
-  ! Tegen fields that are in terms of optical depth, and converted to
-  ! mass mixing ratio via the relevant mass-extinction coefficient
-  INTEGER, PARAMETER :: ITYPE_TROP_BG_AER = 8 ! hydrophobic organic
-  INTEGER, PARAMETER :: ITYPE_STRAT_BG_AER=12 ! non-absorbing sulphate
-  REAL(KIND=JPRB)    :: TROP_BG_AER_MASS_EXT
-  REAL(KIND=JPRB)    :: STRAT_BG_AER_MASS_EXT
-
-CONTAINS
-
-  ! This routine copies information between the IFS radiation
-  ! configuration (stored in global variables) and the radiation
-  ! configuration of the modular radiation scheme (stored in
-  ! rad_config).  The optional input logical LOUTPUT controls whether
-  ! to print lots of information during the setup stage (default is
-  ! no).
-  SUBROUTINE SETUP_RADIATION_SCHEME(LOUTPUT)
-
-    USE YOMHOOK,  ONLY : LHOOK, DR_HOOK
-    USE YOMLUN,   ONLY : NULNAM, NULOUT, NULERR
-    USE YOESRTWN, ONLY : NMPSRTM
-!    USE YOERAD,   ONLY : YRERAD
-    USE MODD_PARAM_ECRAD_n , ONLY : LAPPROXLWUPDATE, NAERMACC, NLIQOPT, NICEOPT,  &
-                                 &  NLWSOLVER, NSWSOLVER, NLWSCATTERING, NOVLP,CDATADIR 
-   
-    USE radiation_interface,      ONLY : setup_radiation
-    USE radiation_aerosol_optics, ONLY : dry_aerosol_sw_mass_extinction
-
-!#include "posname.intfb.h"
-
-    ! Whether or not to provide information on the radiation scheme
-    ! configuration
-    LOGICAL, INTENT(IN), OPTIONAL :: LOUTPUT
-
-    ! Verbosity of configuration information 0=none, 1=warning,
-    ! 2=info, 3=progress, 4=detailed, 5=debug
-    INTEGER :: IVERBOSESETUP
-    INTEGER :: ISTAT
-
-    REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',0,ZHOOK_HANDLE)
-
-    ! *** GENERAL SETUP ***
-
-    ! Configure verbosity of setup of radiation scheme
-    IVERBOSESETUP = 4 ! Provide plenty of information
-    IF (PRESENT(LOUTPUT)) THEN
-      IF (.NOT. LOUTPUT) THEN
-        IVERBOSESETUP = 1 ! Warnings and errors only
-      ENDIF
-    ENDIF
-    rad_config%iverbosesetup = IVERBOSESETUP
-
-    IF (IVERBOSESETUP > 1) THEN
-      WRITE(NULOUT,'(a)') '-------------------------------------------------------------------------------'
-      WRITE(NULOUT,'(a)') 'RADIATION_SETUP'
-    ENDIF
-
-    ! Normal operation of the radiation scheme displays only errors
-    ! and warnings
-    rad_config%iverbose = 1
-
-    ! For the time being, ensure a valid default directory name
-    rad_config%directory_name = CDATADIR
-
-    ! Do we do Hogan and Bozzo (2014) approximate longwave updates?
-    rad_config%do_lw_derivatives = LAPPROXLWUPDATE
-
-    ! Surface spectral fluxes are needed for spectral shortwave albedo
-    ! calculation
-    rad_config%do_surface_sw_spectral_flux = .TRUE.
-
-
-    ! *** SETUP GAS OPTICS ***
-
-    ! Assume IFS has already set-up RRTM, so the setup_radiation
-    ! routine below does not have to
-    rad_config%do_setup_ifsrrtm = .FALSE.
-!     rad_config%do_setup_ifsrrtm = .TRUE.
-
-
-    ! *** SETUP CLOUD OPTICS ***
-
-    ! Setup liquid optics
-    IF (NLIQOPT == 2) THEN
-      rad_config%i_liq_model = ILiquidModelSlingo
-    ELSEIF (NLIQOPT == 3) THEN
-      rad_config%i_liq_model = ILiquidModelSOCRATES
-    ELSE
-      WRITE(NULERR,'(a,i0)') 'Unavailable liquid optics model in modular radiation scheme: NLIQOPT=', &
-           &  NLIQOPT
-      CALL ABOR1('RADIATION_SETUP: error interpreting NLIQOPT')   
-    ENDIF
-
-    ! Setup ice optics
-    IF (NICEOPT == 3) THEN
-      rad_config%i_ice_model = IIceModelFu
-    ELSEIF (NICEOPT == 4) THEN
-      rad_config%i_ice_model = IIceModelBaran
-    ELSE
-      WRITE(NULERR,'(a,i0)') 'Unavailable ice optics model in modular radiation scheme: NICEOPT=', &
-           &  NICEOPT
-      CALL ABOR1('RADIATION_SETUP: error interpreting NICEOPT')   
-    ENDIF
-
-    ! For consistency with earlier versions of the IFS radiation
-    ! scheme, we perform shortwave delta-Eddington scaling *after* the
-    ! merge of the cloud, aerosol and gas optical properties.  Set
-    ! this to "false" to do the scaling on the cloud and aerosol
-    ! properties separately before merging with gases. Note that this
-    ! is not compatible with the SPARTACUS solver.
-    rad_config%do_sw_delta_scaling_with_gases = .TRUE.
-
-    ! Use Exponential-Exponential cloud overlap to match original IFS
-    ! implementation of Raisanen cloud generator
-    ! rad_config%i_overlap_scheme = IOverlapExponential
-    rad_config%i_overlap_scheme = NOVLP ! QL
-
-    ! *** SETUP AEROSOLS ***
-
-    rad_config%use_aerosols = .TRUE.
-
-    IF (NAERMACC > 0) THEN
-      ! Using MACC climatology - in this case the aerosol optics file
-      ! will be chosen automatically
-
-      ! 12 IFS aerosol classes: 1-3 Sea salt, 4-6 Boucher desert dust,
-      ! 7 hydrophilic organics, 8 hydrophobic organics, 9&10
-      ! hydrophobic black carbon, 11 ammonium sulphate, 12 inactive
-      ! SO2
-      rad_config%n_aerosol_types = 12
-
-      ! Indices to the aerosol optical properties in
-      ! aerosol_ifs_rrtm_*.nc, for each class, where negative numbers
-      ! index hydrophilic aerosol types and positive numbers index
-      ! hydrophobic aerosol types
-      rad_config%i_aerosol_type_map = 0 ! There can be up to 256 types
-      rad_config%i_aerosol_type_map(1:12) = (/ &
-           &  -1, &  ! Sea salt, size bin 1 (OPAC)
-           &  -2, &  ! Sea salt, size bin 2 (OPAC)
-           &  -3, &  ! Sea salt, size bin 3 (OPAC)
-           &   7, &  ! Desert dust, size bin 1 (Woodward 2001)
-           &   8, &  ! Desert dust, size bin 2 (Woodward 2001)
-           &   9, &  ! Desert dust, size bin 3 (Woodward 2001)
-           &  -4, &  ! Hydrophilic organic matter (OPAC)
-           &  10, &  ! Hydrophobic organic matter (OPAC)
-           &  11, &  ! Black carbon (Boucher)
-           &  11, &  ! Black carbon (Boucher)
-           &  -5, &  ! Ammonium sulphate (OPAC)
-           &  14 /)  ! Stratospheric sulphate (hand edited from OPAC)
-
-      ! Background aerosol mass-extinction coefficients are obtained
-      ! after the configuration files have been read - see later in
-      ! this routine.
-
-    ELSE
-      ! Using Tegen climatology
-      rad_config%n_aerosol_types = 6
-      rad_config%i_aerosol_type_map = 0 ! There can be up to 256 types
-      rad_config%i_aerosol_type_map(1:6) = (/ &
-           &  1, &  ! Continental background
-           &  2, &  ! Maritime
-           &  3, &  ! Desert
-           &  4, &  ! Urban
-           &  5, &  ! Volcanic active
-           &  6 /)  ! Stratospheric background
-
-      ! Manually set the aerosol optics file name (the directory will
-      ! be added automatically)
-      rad_config%aerosol_optics_override_file_name = 'aerosol_ifs_rrtm_tegen.nc'
-    ENDIF
-
-    ! *** SETUP SOLVER ***
-
-    ! 3D effects are off by default
-    rad_config%do_3d_effects = .FALSE.
-
-    ! Select longwave solver
-    SELECT CASE (NLWSOLVER)
-    CASE(0)
-      rad_config%i_solver_lw = ISolverMcICA
-    CASE(1)
-      rad_config%i_solver_lw = ISolverSpartacus
-    CASE(2)
-      rad_config%i_solver_lw = ISolverSpartacus
-      rad_config%do_3d_effects = .TRUE.
-    CASE DEFAULT
-      WRITE(NULERR,'(a,i0)') 'Unknown value for NLWSOLVER: ', NLWSOLVER
-      CALL ABOR1('RADIATION_SETUP: error interpreting NLWSOLVER')
-    END SELECT
-
-    ! Select shortwave solver
-
-    SELECT CASE (NSWSOLVER)
-    CASE(0)
-      rad_config%i_solver_sw = ISolverMcICA
-    CASE(1)
-      rad_config%i_solver_sw = ISolverSpartacus
-      rad_config%do_3d_effects = .FALSE.
-      IF (NLWSOLVER == 2) THEN
-        CALL ABOR1('RADIATION_SETUP: cannot represent 3D effects in LW but not SW')
-      ENDIF
-    CASE(2)
-      rad_config%i_solver_sw = ISolverSpartacus
-      rad_config%do_3d_effects = .TRUE.
-      IF (NLWSOLVER == 1) THEN
-        CALL ABOR1('RADIATION_SETUP: cannot represent 3D effects in SW but not LW')
-      ENDIF
-    CASE DEFAULT
-      WRITE(NULERR,'(a,i0)') 'Unknown value for NSWSOLVER: ', NSWSOLVER
-      CALL ABOR1('RADIATION_SETUP: error interpreting NSWSOLVER')
-    END SELECT
-
-    ! SPARTACUS solver requires delta scaling to be done separately
-    ! for clouds & aerosols
-    IF (rad_config%i_solver_sw == ISolverSpartacus) THEN
-      rad_config%do_sw_delta_scaling_with_gases = .FALSE.
-    ENDIF
-
-    ! Do we represent longwave scattering?
-    rad_config%do_lw_cloud_scattering = .FALSE.
-    rad_config%do_lw_aerosol_scattering = .FALSE.
-    SELECT CASE (NLWSCATTERING)
-    CASE(1)
-      rad_config%do_lw_cloud_scattering = .TRUE.
-    CASE(2)
-      rad_config%do_lw_cloud_scattering = .TRUE.
-      IF (NAERMACC > 0) THEN
-        ! Tegen climatology omits data required to do longwave
-        ! scattering by aerosols, so only turn this on with a more
-        ! recent scattering database
-        rad_config%do_lw_aerosol_scattering = .TRUE.
-      ENDIF
-    END SELECT
-
-
-    ! *** IMPLEMENT SETTINGS ***
-
-    ! For advanced configuration, the configuration data for the
-    ! "radiation" project can specified directly in the namelist.
-    ! However, the variable naming convention is not consistent with
-    ! the rest of the IFS.  For basic configuration there are specific
-    ! variables in the NAERAD namelist available in the YRERAD
-    ! structure.
-    ! CALL POSNAME(NULNAM, 'RADIATION', ISTAT)
-    ISTAT = 1 ! QL no namelist used, all in NAM_PARAM_ECRAD
-    SELECT CASE (ISTAT)
-      CASE(0)
-        CALL rad_config%read(unit=NULNAM)
-      CASE(1)
-        WRITE(NULOUT,'(a)') 'Namelist RADIATION not found, using settings from NAERAD only'
-      CASE DEFAULT
-        CALL ABOR1('RADIATION_SETUP: error reading RADIATION section of namelist file')
-    END SELECT
-
-    ! Print configuration
-    IF (IVERBOSESETUP > 1) THEN
-      WRITE(NULOUT,'(a)') 'Radiation scheme settings:'
-      CALL rad_config%print(IVERBOSE=IVERBOSESETUP)
-    ENDIF
-
-    ! Use configuration data to set-up radiation scheme, including
-    ! reading scattering datafiles
-    CALL setup_radiation(rad_config)
-
-    ! Populate the mapping between the 14 RRTM shortwave bands and the
-    ! 6 albedo inputs. The mapping according to the stated wavelength
-    ! ranges of the 6-band model does not match the hard-wired mapping
-    ! in NMPSRTM, but only the hard-wired values produce sensible
-    ! results...
-    rad_config%i_albedo_from_band_sw = NMPSRTM
-    !    call rad_config%define_sw_albedo_intervals(6, &
-    !         &  (/ 0.25e-6_jprb, 0.44e-6_jprb, 1.19e-6_jprb, &
-    !         &     2.38e-6_jprb, 4.00e-6_jprb /),  (/ 1,2,3,4,5,6 /))
-    
-    ! Likewise between the 16 RRTM longwave bands and the 2 emissivity
-    ! inputs (info taken from rrtm_ecrt_140gp_mcica.F90) representing
-    ! outside and inside the window region of the spectrum
-    ! rad_config%i_emiss_from_band_lw = (/ 1,1,1,1,1,2,2,2,1,1,1,1,1,1,1,1 /)
-    call rad_config%define_lw_emiss_intervals(3, &
-         &  (/ 8.0e-6_jprb,13.0e-6_jprb /),  (/ 1,2,1 /))
-
-    ! Get spectral weightings for UV and PAR
-    call rad_config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb, &
-         &  NWEIGHT_UV, IBAND_UV, WEIGHT_UV, 'ultraviolet')
-    call rad_config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb, &
-         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
-         &  'photosynthetically active radiation, PAR')
-
-    IF (NAERMACC > 0) THEN
-      ! With the MACC aerosol climatology we need to add in the
-      ! background aerosol afterwards using the Tegen arrays.  In this
-      ! case we first configure the background aerosol mass-extinction
-      ! coefficient at 550 nm, which corresponds to the 10th RRTMG
-      ! shortwave band.
-      TROP_BG_AER_MASS_EXT  = dry_aerosol_sw_mass_extinction(rad_config, &
-           &                                   ITYPE_TROP_BG_AER, 10)
-      STRAT_BG_AER_MASS_EXT = dry_aerosol_sw_mass_extinction(rad_config, &
-           &                                   ITYPE_STRAT_BG_AER, 10)
-      
-      WRITE(NULOUT,'(a,i0)') 'Tropospheric bacground uses aerosol type ', &
-           &                 ITYPE_TROP_BG_AER
-      WRITE(NULOUT,'(a,i0)') 'Stratospheric bacground uses aerosol type ', &
-           &                 ITYPE_STRAT_BG_AER
-    ENDIF
-      
-    IF (IVERBOSESETUP > 1) THEN
-      WRITE(NULOUT,'(a)') '-------------------------------------------------------------------------------'
-    ENDIF
-
-    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',1,ZHOOK_HANDLE)
-
-  END SUBROUTINE SETUP_RADIATION_SCHEME
-
-END MODULE RADIATION_SETUP
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomcst.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomcst.F90
deleted file mode 100644
index a26a5740d..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomcst.F90
+++ /dev/null
@@ -1,29 +0,0 @@
-MODULE YOMCST
-
-USE PARKIND1  ,ONLY : JPRB
-
-IMPLICIT NONE
-
-SAVE
-
-! * RPI          : number Pi
-REAL(KIND=JPRB), PARAMETER :: RPI = 3.14159265358979323846_JPRB
-! * RSIGMA       : Stefan-Bolzman constant
-REAL(KIND=JPRB), PARAMETER :: RSIGMA = 5.67037321e-8_JPRB ! W m-2 K-4
-! * RG           : gravity constant
-REAL(KIND=JPRB), PARAMETER :: RG = 9.80665_JPRB ! m s-2
-! * RD           : R_dry (dry air constant)
-REAL(KIND=JPRB), PARAMETER :: RD = 287.058_JPRB! J kg-1 K-1
-! * RMD          : dry air molar mass
-REAL(KIND=JPRB), PARAMETER :: RMD = 28.9644_JPRB
-! * RMV          : vapour water molar mass
-REAL(KIND=JPRB), PARAMETER :: RMV = 18.0153_JPRB
-! * RMO3         : ozone molar mass
-REAL(KIND=JPRB), PARAMETER :: RMO3 = 47.9942_JPRB
-! * RI0          : solar constant
-REAL(KIND=JPRB), PARAMETER :: RI0 = 1366.0_JPRB
-! * RDAY          : day duration in s
-REAL(KIND=JPRB), PARAMETER :: RDAY = 86400_JPRB
-! * RTT          : freezing temperature
-REAL(KIND=JPRB), PARAMETER :: RTT=273.16_JPRB
-END MODULE YOMCST
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomlun.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomlun.F90
deleted file mode 100644
index f87588230..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsaux/yomlun.F90
+++ /dev/null
@@ -1,15 +0,0 @@
-MODULE YOMLUN
-
-USE PARKIND1,      ONLY : JPIM
-USE YOMLUN_IFSAUX, ONLY : NULOUT, NULERR
-
-IMPLICIT NONE
-
-SAVE
-
-INTEGER(KIND=JPIM) :: NULRAD = 145
-
-INTEGER(KIND=JPIM) :: NULNAM  =  4
-
-!     ------------------------------------------------------------------
-END MODULE YOMLUN
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/rrtm_kgb1.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/rrtm_kgb1.F90
deleted file mode 100644
index e2b5bca94..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/rrtm_kgb1.F90
+++ /dev/null
@@ -1,355 +0,0 @@
-SUBROUTINE RRTM_KGB1(DIRECTORY)
-
-!     Originally by Eli J. Mlawer, Atmospheric & Environmental Research.
-!     BAND 1:  10-250 cm-1 (low - H2O; high - H2O)
-!     Reformatted for F90 by JJMorcrette, ECMWF
-!     R. Elkhatib 12-10-2005 Split for faster and more robust compilation.
-!     G.Mozdzynski March 2011 read constants from files
-!     ABozzo May 2013 update to RRTMG v4.85
-!     band 1:  10-350 cm-1
-!     T. Wilhelmsson and K. Yessad (Oct 2013) Geometry and setup refactoring.
-!     ------------------------------------------------------------------
-
-USE PARKIND1  ,ONLY : JPRB
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE YOMLUN    ,ONLY : NULRAD , NULOUT
-USE MPL_MODULE,ONLY : MPL_BROADCAST
-USE YOMTAG    ,ONLY : MTAGRAD
-USE YOMMP0    , ONLY : NPROC, MYPROC
-
-USE YOERRTO1 , ONLY : KAO     ,KBO     ,SELFREFO   ,FRACREFAO ,&
- & FRACREFBO  ,FORREFO, KAO_MN2, KBO_MN2
-
-!     ------------------------------------------------------------------
-
-IMPLICIT NONE
-
-CHARACTER(LEN=*), INTENT(IN) :: DIRECTORY
-
-!CHARACTER(LEN = 80) :: CLZZZ
-CHARACTER(LEN = 255) :: CLF1
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-#include "abor1.intfb.h"
-
-IF (LHOOK) CALL DR_HOOK('RRTM_KGB1',0,ZHOOK_HANDLE)
-
-IF( MYPROC==1 )THEN
-  !CALL GETENV("DATA",CLZZZ)
-  !IF(CLZZZ /= " ") THEN
-  !  CLF1=TRIM(CLZZZ) // "/RADRRTM"
-  CLF1 = DIRECTORY // "/RADRRTM"
-  WRITE(NULOUT,'(A,A)') 'Reading ',TRIM(CLF1)
-  ! RRTM and SRTM files from ecrad are in big-endian format.
-  ! Here they are read as big-endian at opening because otherwise MNH assumes littel-endian
-  ! No need for complation option export GFORTRAN_CONVERT_UNIT="little_endian;big_endian:145"
-  OPEN(NULRAD,FILE=TRIM(CLF1),FORM="UNFORMATTED",ACTION="READ",access='sequential',ERR=1000,CONVERT="BIG_ENDIAN")
-  !ELSE
-  !  OPEN(NULRAD,FILE='RADRRTM',FORM="UNFORMATTED",ACTION="READ",ERR=1000)
-  !ENDIF
-  READ(NULRAD,ERR=1001) KAO,KBO
-ENDIF
-IF( NPROC>1 )THEN
-  CALL MPL_BROADCAST (KAO,MTAGRAD,1,CDSTRING='RRTM_KGB1:')
-  CALL MPL_BROADCAST (KBO,MTAGRAD,1,CDSTRING='RRTM_KGB1:')
-ENDIF
-
-! Planck fraction mapping level: P = 212.7250 mbar, T = 223.06 K
-FRACREFAO(:) = (/ &
- & 2.1227E-01_JPRB,1.8897E-01_JPRB,1.3934E-01_JPRB,1.1557E-01_JPRB,9.5282E-02_JPRB,8.3359E-02_JPRB, &
- & 6.5333E-02_JPRB,5.2016E-02_JPRB,3.4272E-02_JPRB,4.0257E-03_JPRB,3.1857E-03_JPRB,2.6014E-03_JPRB, &
- & 1.9141E-03_JPRB,1.2612E-03_JPRB,5.3169E-04_JPRB,7.6476E-05_JPRB/)
-
-! Planck fraction mapping level: P = 212.7250 mbar, T = 223.06 K
-! These Planck fractions were calculated using lower atmosphere
-! parameters.
-FRACREFBO(:) = (/ &
- & 2.1227E-01_JPRB,1.8897E-01_JPRB,1.3934E-01_JPRB,1.1557E-01_JPRB,9.5282E-02_JPRB,8.3359E-02_JPRB, &
- & 6.5333E-02_JPRB,5.2016E-02_JPRB,3.4272E-02_JPRB,4.0257E-03_JPRB,3.1857E-03_JPRB,2.6014E-03_JPRB, &
- & 1.9141E-03_JPRB,1.2612E-03_JPRB,5.3169E-04_JPRB,7.6476E-05_JPRB/)
-
-!     The array FORREFO contains the coefficient of the water vapor
-!     foreign-continuum (including the energy term).  The first 
-!     index refers to reference temperature (296,260,224,260) and 
-!     pressure (970,475,219,3 mbar) levels.  The second index 
-!     runs over the g-channel (1 to 16).
-
-      FORREFO(1,:) = (/ &
-     & 3.6742e-02_JPRB,1.0664e-01_JPRB,2.6132e-01_JPRB,2.7906e-01_JPRB,2.8151e-01_JPRB,2.7465e-01_JPRB, &
-     & 2.8530e-01_JPRB,2.9123e-01_JPRB,3.0697e-01_JPRB,3.1801e-01_JPRB,3.2444e-01_JPRB,2.7746e-01_JPRB, &
-     & 3.1994e-01_JPRB,2.9750e-01_JPRB,2.1226e-01_JPRB,1.2847e-01_JPRB/)
-      FORREFO(2,:) = (/ &
-     & 4.0450e-02_JPRB,1.1085e-01_JPRB,2.9205e-01_JPRB,3.1934e-01_JPRB,3.1739e-01_JPRB,3.1450e-01_JPRB, &
-     & 3.2797e-01_JPRB,3.2223e-01_JPRB,3.3099e-01_JPRB,3.4800e-01_JPRB,3.4046e-01_JPRB,3.5700e-01_JPRB, &
-     & 3.8264e-01_JPRB,3.6679e-01_JPRB,3.3481e-01_JPRB,3.2113e-01_JPRB/)
-      FORREFO(3,:) = (/ &
-     & 4.6952e-02_JPRB,1.1999e-01_JPRB,3.1473e-01_JPRB,3.7015e-01_JPRB,3.6913e-01_JPRB,3.6352e-01_JPRB, &
-     & 3.7754e-01_JPRB,3.7402e-01_JPRB,3.7113e-01_JPRB,3.7720e-01_JPRB,3.8365e-01_JPRB,4.0876e-01_JPRB, &
-     & 4.2968e-01_JPRB,4.4186e-01_JPRB,4.3468e-01_JPRB,4.7083e-01_JPRB/)
-      FORREFO(4,:) = (/ &
-     & 7.0645e-02_JPRB,1.6618e-01_JPRB,2.8516e-01_JPRB,3.1819e-01_JPRB,3.0131e-01_JPRB,2.9552e-01_JPRB, &
-     & 2.8972e-01_JPRB,2.9348e-01_JPRB,2.8668e-01_JPRB,2.8483e-01_JPRB,2.8130e-01_JPRB,2.7757e-01_JPRB, &
-     & 2.9735e-01_JPRB,3.1684e-01_JPRB,3.0681e-01_JPRB,3.6778e-01_JPRB/)
-
-
-!     ------------------------------------------------------------------
-
-!     The array KAO contains absorption coefs at the 16 chosen g-values 
-!     for a range of pressure levels > ~100mb and temperatures.  The first
-!     index in the array, JT, which runs from 1 to 5, corresponds to 
-!     different temperatures.  More specifically, JT = 3 means that the 
-!     data are for the corresponding TREF for this  pressure level, 
-!     JT = 2 refers to the temperatureTREF-15, JT = 1 is for TREF-30, 
-!     JT = 4 is for TREF+15, and JT = 5 is for TREF+30.  The second 
-!     index, JP, runs from 1 to 13 and refers to the corresponding 
-!     pressure level in PREF (e.g. JP = 1 is for a pressure of 1053.63 mb).  
-!     The third index, IG, goes from 1 to 16, and tells us which 
-!     g-interval the absorption coefficients are for.
-
-
-
-!     The array KBO contains absorption coefs at the 16 chosen g-values 
-!     for a range of pressure levels < ~100mb and temperatures. The first 
-!     index in the array, JT, which runs from 1 to 5, corresponds to 
-!     different temperatures.  More specifically, JT = 3 means that the 
-!     data are for the reference temperature TREF for this pressure 
-!     level, JT = 2 refers to the temperature TREF-15, JT = 1 is for
-!     TREF-30, JT = 4 is for TREF+15, and JT = 5 is for TREF+30.  
-!     The second index, JP, runs from 13 to 59 and refers to the JPth
-!     reference pressure level (see taumol.f for the value of these
-!     pressure levels in mb).  The third index, IG, goes from 1 to 16,
-!     and tells us which g-interval the absorption coefficients are for.
-
-
-
-     KAO_MN2(:, 1) = (/ &
-     & 5.12042E-08_JPRB, 5.51239E-08_JPRB, 5.93436E-08_JPRB, 6.38863E-08_JPRB, 6.87767E-08_JPRB, &
-     & 7.40415E-08_JPRB, 7.97093E-08_JPRB, 8.58110E-08_JPRB, 9.23797E-08_JPRB, 9.94513E-08_JPRB, &
-     & 1.07064E-07_JPRB, 1.15260E-07_JPRB, 1.24083E-07_JPRB, 1.33581E-07_JPRB, 1.43807E-07_JPRB, &
-     & 1.54815E-07_JPRB, 1.66666E-07_JPRB, 1.79424E-07_JPRB, 1.93159E-07_JPRB/)
-      KAO_MN2(:, 2) = (/ &
-     & 2.30938E-07_JPRB, 2.41696E-07_JPRB, 2.52955E-07_JPRB, 2.64738E-07_JPRB, 2.77071E-07_JPRB, &
-     & 2.89978E-07_JPRB, 3.03486E-07_JPRB, 3.17623E-07_JPRB, 3.32419E-07_JPRB, 3.47904E-07_JPRB, &
-     & 3.64111E-07_JPRB, 3.81072E-07_JPRB, 3.98824E-07_JPRB, 4.17402E-07_JPRB, 4.36846E-07_JPRB, &
-     & 4.57196E-07_JPRB, 4.78494E-07_JPRB, 5.00784E-07_JPRB, 5.24112E-07_JPRB/)
-      KAO_MN2(:, 3) = (/ &
-     & 6.70458E-07_JPRB, 7.04274E-07_JPRB, 7.39795E-07_JPRB, 7.77109E-07_JPRB, 8.16304E-07_JPRB, &
-     & 8.57476E-07_JPRB, 9.00724E-07_JPRB, 9.46154E-07_JPRB, 9.93876E-07_JPRB, 1.04400E-06_JPRB, &
-     & 1.09666E-06_JPRB, 1.15197E-06_JPRB, 1.21008E-06_JPRB, 1.27111E-06_JPRB, 1.33522E-06_JPRB, &
-     & 1.40256E-06_JPRB, 1.47331E-06_JPRB, 1.54761E-06_JPRB, 1.62567E-06_JPRB/)
-      KAO_MN2(:, 4) = (/ &
-     & 1.84182E-06_JPRB, 1.89203E-06_JPRB, 1.94360E-06_JPRB, 1.99658E-06_JPRB, 2.05101E-06_JPRB, &
-     & 2.10692E-06_JPRB, 2.16435E-06_JPRB, 2.22335E-06_JPRB, 2.28396E-06_JPRB, 2.34622E-06_JPRB, &
-     & 2.41017E-06_JPRB, 2.47587E-06_JPRB, 2.54337E-06_JPRB, 2.61270E-06_JPRB, 2.68392E-06_JPRB, &
-     & 2.75708E-06_JPRB, 2.83224E-06_JPRB, 2.90944E-06_JPRB, 2.98875E-06_JPRB/)
-      KAO_MN2(:, 5) = (/ &
-     & 3.41996E-06_JPRB, 3.32758E-06_JPRB, 3.23770E-06_JPRB, 3.15024E-06_JPRB, 3.06515E-06_JPRB, &
-     & 2.98235E-06_JPRB, 2.90180E-06_JPRB, 2.82341E-06_JPRB, 2.74715E-06_JPRB, 2.67294E-06_JPRB, &
-     & 2.60074E-06_JPRB, 2.53049E-06_JPRB, 2.46214E-06_JPRB, 2.39563E-06_JPRB, 2.33092E-06_JPRB, &
-     & 2.26796E-06_JPRB, 2.20670E-06_JPRB, 2.14709E-06_JPRB, 2.08910E-06_JPRB/)
-      KAO_MN2(:, 6) = (/ &
-     & 3.38746E-06_JPRB, 3.25966E-06_JPRB, 3.13669E-06_JPRB, 3.01836E-06_JPRB, 2.90449E-06_JPRB, &
-     & 2.79491E-06_JPRB, 2.68947E-06_JPRB, 2.58801E-06_JPRB, 2.49037E-06_JPRB, 2.39642E-06_JPRB, &
-     & 2.30601E-06_JPRB, 2.21902E-06_JPRB, 2.13530E-06_JPRB, 2.05475E-06_JPRB, 1.97723E-06_JPRB, &
-     & 1.90264E-06_JPRB, 1.83086E-06_JPRB, 1.76179E-06_JPRB, 1.69532E-06_JPRB/)
-      KAO_MN2(:, 7) = (/ &
-     & 3.17530E-06_JPRB, 3.07196E-06_JPRB, 2.97199E-06_JPRB, 2.87527E-06_JPRB, 2.78170E-06_JPRB, &
-     & 2.69118E-06_JPRB, 2.60360E-06_JPRB, 2.51887E-06_JPRB, 2.43690E-06_JPRB, 2.35759E-06_JPRB, &
-     & 2.28087E-06_JPRB, 2.20664E-06_JPRB, 2.13483E-06_JPRB, 2.06536E-06_JPRB, 1.99814E-06_JPRB, &
-     & 1.93312E-06_JPRB, 1.87021E-06_JPRB, 1.80934E-06_JPRB, 1.75046E-06_JPRB/)
-      KAO_MN2(:, 8) = (/ &
-     & 2.84701E-06_JPRB, 2.77007E-06_JPRB, 2.69521E-06_JPRB, 2.62237E-06_JPRB, 2.55150E-06_JPRB, &
-     & 2.48254E-06_JPRB, 2.41545E-06_JPRB, 2.35017E-06_JPRB, 2.28666E-06_JPRB, 2.22486E-06_JPRB, &
-     & 2.16473E-06_JPRB, 2.10623E-06_JPRB, 2.04930E-06_JPRB, 1.99392E-06_JPRB, 1.94003E-06_JPRB, &
-     & 1.88760E-06_JPRB, 1.83659E-06_JPRB, 1.78695E-06_JPRB, 1.73866E-06_JPRB/)
-      KAO_MN2(:, 9) = (/ &
-     & 2.79917E-06_JPRB, 2.73207E-06_JPRB, 2.66658E-06_JPRB, 2.60266E-06_JPRB, 2.54027E-06_JPRB, &
-     & 2.47937E-06_JPRB, 2.41994E-06_JPRB, 2.36192E-06_JPRB, 2.30530E-06_JPRB, 2.25004E-06_JPRB, &
-     & 2.19610E-06_JPRB, 2.14346E-06_JPRB, 2.09208E-06_JPRB, 2.04193E-06_JPRB, 1.99298E-06_JPRB, &
-     & 1.94520E-06_JPRB, 1.89857E-06_JPRB, 1.85306E-06_JPRB, 1.80864E-06_JPRB/)
-      KAO_MN2(:,10) = (/ &
-     & 2.74910E-06_JPRB, 2.64462E-06_JPRB, 2.54412E-06_JPRB, 2.44743E-06_JPRB, 2.35442E-06_JPRB, &
-     & 2.26495E-06_JPRB, 2.17887E-06_JPRB, 2.09606E-06_JPRB, 2.01641E-06_JPRB, 1.93978E-06_JPRB, &
-     & 1.86606E-06_JPRB, 1.79514E-06_JPRB, 1.72692E-06_JPRB, 1.66129E-06_JPRB, 1.59815E-06_JPRB, &
-     & 1.53742E-06_JPRB, 1.47899E-06_JPRB, 1.42278E-06_JPRB, 1.36871E-06_JPRB/)
-      KAO_MN2(:,11) = (/ &
-     & 2.63952E-06_JPRB, 2.60263E-06_JPRB, 2.56626E-06_JPRB, 2.53039E-06_JPRB, 2.49503E-06_JPRB, &
-     & 2.46016E-06_JPRB, 2.42578E-06_JPRB, 2.39188E-06_JPRB, 2.35845E-06_JPRB, 2.32549E-06_JPRB, &
-     & 2.29299E-06_JPRB, 2.26094E-06_JPRB, 2.22934E-06_JPRB, 2.19819E-06_JPRB, 2.16747E-06_JPRB, &
-     & 2.13717E-06_JPRB, 2.10731E-06_JPRB, 2.07786E-06_JPRB, 2.04882E-06_JPRB/)
-      KAO_MN2(:,12) = (/ &
-     & 2.94106E-06_JPRB, 2.82819E-06_JPRB, 2.71966E-06_JPRB, 2.61528E-06_JPRB, 2.51492E-06_JPRB, &
-     & 2.41841E-06_JPRB, 2.32560E-06_JPRB, 2.23635E-06_JPRB, 2.15053E-06_JPRB, 2.06800E-06_JPRB, &
-     & 1.98863E-06_JPRB, 1.91232E-06_JPRB, 1.83893E-06_JPRB, 1.76836E-06_JPRB, 1.70049E-06_JPRB, &
-     & 1.63524E-06_JPRB, 1.57248E-06_JPRB, 1.51214E-06_JPRB, 1.45411E-06_JPRB/)
-      KAO_MN2(:,13) = (/ &
-     & 2.94607E-06_JPRB, 2.87369E-06_JPRB, 2.80309E-06_JPRB, 2.73422E-06_JPRB, 2.66705E-06_JPRB, &
-     & 2.60152E-06_JPRB, 2.53760E-06_JPRB, 2.47526E-06_JPRB, 2.41445E-06_JPRB, 2.35513E-06_JPRB, &
-     & 2.29726E-06_JPRB, 2.24082E-06_JPRB, 2.18577E-06_JPRB, 2.13207E-06_JPRB, 2.07969E-06_JPRB, &
-     & 2.02859E-06_JPRB, 1.97875E-06_JPRB, 1.93014E-06_JPRB, 1.88272E-06_JPRB/)
-      KAO_MN2(:,14) = (/ &
-     & 2.58051E-06_JPRB, 2.48749E-06_JPRB, 2.39782E-06_JPRB, 2.31139E-06_JPRB, 2.22807E-06_JPRB, &
-     & 2.14775E-06_JPRB, 2.07033E-06_JPRB, 1.99570E-06_JPRB, 1.92376E-06_JPRB, 1.85441E-06_JPRB, &
-     & 1.78756E-06_JPRB, 1.72313E-06_JPRB, 1.66101E-06_JPRB, 1.60114E-06_JPRB, 1.54342E-06_JPRB, &
-     & 1.48778E-06_JPRB, 1.43415E-06_JPRB, 1.38245E-06_JPRB, 1.33262E-06_JPRB/)
-      KAO_MN2(:,15) = (/ &
-     & 3.03447E-06_JPRB, 2.88559E-06_JPRB, 2.74401E-06_JPRB, 2.60938E-06_JPRB, 2.48135E-06_JPRB, &
-     & 2.35961E-06_JPRB, 2.24384E-06_JPRB, 2.13375E-06_JPRB, 2.02906E-06_JPRB, 1.92951E-06_JPRB, &
-     & 1.83484E-06_JPRB, 1.74481E-06_JPRB, 1.65921E-06_JPRB, 1.57780E-06_JPRB, 1.50039E-06_JPRB, &
-     & 1.42677E-06_JPRB, 1.35677E-06_JPRB, 1.29020E-06_JPRB, 1.22690E-06_JPRB/)
-      KAO_MN2(:,16) = (/ &
-     & 1.48655E-06_JPRB, 1.48283E-06_JPRB, 1.47913E-06_JPRB, 1.47543E-06_JPRB, 1.47174E-06_JPRB, &
-     & 1.46806E-06_JPRB, 1.46439E-06_JPRB, 1.46072E-06_JPRB, 1.45707E-06_JPRB, 1.45343E-06_JPRB, &
-     & 1.44979E-06_JPRB, 1.44617E-06_JPRB, 1.44255E-06_JPRB, 1.43894E-06_JPRB, 1.43534E-06_JPRB, &
-     & 1.43176E-06_JPRB, 1.42817E-06_JPRB, 1.42460E-06_JPRB, 1.42104E-06_JPRB/)
-      KBO_MN2(:, 1) = (/ &
-     & 5.12042E-08_JPRB, 5.51239E-08_JPRB, 5.93436E-08_JPRB, 6.38863E-08_JPRB, 6.87767E-08_JPRB, &
-     & 7.40415E-08_JPRB, 7.97093E-08_JPRB, 8.58110E-08_JPRB, 9.23797E-08_JPRB, 9.94513E-08_JPRB, &
-     & 1.07064E-07_JPRB, 1.15260E-07_JPRB, 1.24083E-07_JPRB, 1.33581E-07_JPRB, 1.43807E-07_JPRB, &
-     & 1.54815E-07_JPRB, 1.66666E-07_JPRB, 1.79424E-07_JPRB, 1.93159E-07_JPRB/)
-      KBO_MN2(:, 2) = (/ &
-     & 2.30938E-07_JPRB, 2.41696E-07_JPRB, 2.52955E-07_JPRB, 2.64738E-07_JPRB, 2.77071E-07_JPRB, &
-     & 2.89978E-07_JPRB, 3.03486E-07_JPRB, 3.17623E-07_JPRB, 3.32419E-07_JPRB, 3.47904E-07_JPRB, &
-     & 3.64111E-07_JPRB, 3.81072E-07_JPRB, 3.98824E-07_JPRB, 4.17402E-07_JPRB, 4.36846E-07_JPRB, &
-     & 4.57196E-07_JPRB, 4.78494E-07_JPRB, 5.00784E-07_JPRB, 5.24112E-07_JPRB/)
-      KBO_MN2(:, 3) = (/ &
-     & 6.70458E-07_JPRB, 7.04274E-07_JPRB, 7.39795E-07_JPRB, 7.77109E-07_JPRB, 8.16304E-07_JPRB, &
-     & 8.57476E-07_JPRB, 9.00724E-07_JPRB, 9.46154E-07_JPRB, 9.93876E-07_JPRB, 1.04400E-06_JPRB, &
-     & 1.09666E-06_JPRB, 1.15197E-06_JPRB, 1.21008E-06_JPRB, 1.27111E-06_JPRB, 1.33522E-06_JPRB, &
-     & 1.40256E-06_JPRB, 1.47331E-06_JPRB, 1.54761E-06_JPRB, 1.62567E-06_JPRB/)
-      KBO_MN2(:, 4) = (/ &
-     & 1.84182E-06_JPRB, 1.89203E-06_JPRB, 1.94360E-06_JPRB, 1.99658E-06_JPRB, 2.05101E-06_JPRB, &
-     & 2.10692E-06_JPRB, 2.16435E-06_JPRB, 2.22335E-06_JPRB, 2.28396E-06_JPRB, 2.34622E-06_JPRB, &
-     & 2.41017E-06_JPRB, 2.47587E-06_JPRB, 2.54337E-06_JPRB, 2.61270E-06_JPRB, 2.68392E-06_JPRB, &
-     & 2.75708E-06_JPRB, 2.83224E-06_JPRB, 2.90944E-06_JPRB, 2.98875E-06_JPRB/)
-      KBO_MN2(:, 5) = (/ &
-     & 3.41996E-06_JPRB, 3.32758E-06_JPRB, 3.23770E-06_JPRB, 3.15024E-06_JPRB, 3.06515E-06_JPRB, &
-     & 2.98235E-06_JPRB, 2.90180E-06_JPRB, 2.82341E-06_JPRB, 2.74715E-06_JPRB, 2.67294E-06_JPRB, &
-     & 2.60074E-06_JPRB, 2.53049E-06_JPRB, 2.46214E-06_JPRB, 2.39563E-06_JPRB, 2.33092E-06_JPRB, &
-     & 2.26796E-06_JPRB, 2.20670E-06_JPRB, 2.14709E-06_JPRB, 2.08910E-06_JPRB/)
-      KBO_MN2(:, 6) = (/ &
-     & 3.38746E-06_JPRB, 3.25966E-06_JPRB, 3.13669E-06_JPRB, 3.01836E-06_JPRB, 2.90449E-06_JPRB, &
-     & 2.79491E-06_JPRB, 2.68947E-06_JPRB, 2.58801E-06_JPRB, 2.49037E-06_JPRB, 2.39642E-06_JPRB, &
-     & 2.30601E-06_JPRB, 2.21902E-06_JPRB, 2.13530E-06_JPRB, 2.05475E-06_JPRB, 1.97723E-06_JPRB, &
-     & 1.90264E-06_JPRB, 1.83086E-06_JPRB, 1.76179E-06_JPRB, 1.69532E-06_JPRB/)
-      KBO_MN2(:, 7) = (/ &
-     & 3.17530E-06_JPRB, 3.07196E-06_JPRB, 2.97199E-06_JPRB, 2.87527E-06_JPRB, 2.78170E-06_JPRB, &
-     & 2.69118E-06_JPRB, 2.60360E-06_JPRB, 2.51887E-06_JPRB, 2.43690E-06_JPRB, 2.35759E-06_JPRB, &
-     & 2.28087E-06_JPRB, 2.20664E-06_JPRB, 2.13483E-06_JPRB, 2.06536E-06_JPRB, 1.99814E-06_JPRB, &
-     & 1.93312E-06_JPRB, 1.87021E-06_JPRB, 1.80934E-06_JPRB, 1.75046E-06_JPRB/)
-      KBO_MN2(:, 8) = (/ &
-     & 2.84701E-06_JPRB, 2.77007E-06_JPRB, 2.69521E-06_JPRB, 2.62237E-06_JPRB, 2.55150E-06_JPRB, &
-     & 2.48254E-06_JPRB, 2.41545E-06_JPRB, 2.35017E-06_JPRB, 2.28666E-06_JPRB, 2.22486E-06_JPRB, &
-     & 2.16473E-06_JPRB, 2.10623E-06_JPRB, 2.04930E-06_JPRB, 1.99392E-06_JPRB, 1.94003E-06_JPRB, &
-     & 1.88760E-06_JPRB, 1.83659E-06_JPRB, 1.78695E-06_JPRB, 1.73866E-06_JPRB/)
-      KBO_MN2(:, 9) = (/ &
-     & 2.79917E-06_JPRB, 2.73207E-06_JPRB, 2.66658E-06_JPRB, 2.60266E-06_JPRB, 2.54027E-06_JPRB, &
-     & 2.47937E-06_JPRB, 2.41994E-06_JPRB, 2.36192E-06_JPRB, 2.30530E-06_JPRB, 2.25004E-06_JPRB, &
-     & 2.19610E-06_JPRB, 2.14346E-06_JPRB, 2.09208E-06_JPRB, 2.04193E-06_JPRB, 1.99298E-06_JPRB, &
-     & 1.94520E-06_JPRB, 1.89857E-06_JPRB, 1.85306E-06_JPRB, 1.80864E-06_JPRB/)
-      KBO_MN2(:,10) = (/ &
-     & 2.74910E-06_JPRB, 2.64462E-06_JPRB, 2.54412E-06_JPRB, 2.44743E-06_JPRB, 2.35442E-06_JPRB, &
-     & 2.26495E-06_JPRB, 2.17887E-06_JPRB, 2.09606E-06_JPRB, 2.01641E-06_JPRB, 1.93978E-06_JPRB, &
-     & 1.86606E-06_JPRB, 1.79514E-06_JPRB, 1.72692E-06_JPRB, 1.66129E-06_JPRB, 1.59815E-06_JPRB, &
-     & 1.53742E-06_JPRB, 1.47899E-06_JPRB, 1.42278E-06_JPRB, 1.36871E-06_JPRB/)
-      KBO_MN2(:,11) = (/ &
-     & 2.63952E-06_JPRB, 2.60263E-06_JPRB, 2.56626E-06_JPRB, 2.53039E-06_JPRB, 2.49503E-06_JPRB, &
-     & 2.46016E-06_JPRB, 2.42578E-06_JPRB, 2.39188E-06_JPRB, 2.35845E-06_JPRB, 2.32549E-06_JPRB, &
-     & 2.29299E-06_JPRB, 2.26094E-06_JPRB, 2.22934E-06_JPRB, 2.19819E-06_JPRB, 2.16747E-06_JPRB, &
-     & 2.13717E-06_JPRB, 2.10731E-06_JPRB, 2.07786E-06_JPRB, 2.04882E-06_JPRB/)
-      KBO_MN2(:,12) = (/ &
-     & 2.94106E-06_JPRB, 2.82819E-06_JPRB, 2.71966E-06_JPRB, 2.61528E-06_JPRB, 2.51492E-06_JPRB, &
-     & 2.41841E-06_JPRB, 2.32560E-06_JPRB, 2.23635E-06_JPRB, 2.15053E-06_JPRB, 2.06800E-06_JPRB, &
-     & 1.98863E-06_JPRB, 1.91232E-06_JPRB, 1.83893E-06_JPRB, 1.76836E-06_JPRB, 1.70049E-06_JPRB, &
-     & 1.63524E-06_JPRB, 1.57248E-06_JPRB, 1.51214E-06_JPRB, 1.45411E-06_JPRB/)
-      KBO_MN2(:,13) = (/ &
-     & 2.94607E-06_JPRB, 2.87369E-06_JPRB, 2.80309E-06_JPRB, 2.73422E-06_JPRB, 2.66705E-06_JPRB, &
-     & 2.60152E-06_JPRB, 2.53760E-06_JPRB, 2.47526E-06_JPRB, 2.41445E-06_JPRB, 2.35513E-06_JPRB, &
-     & 2.29726E-06_JPRB, 2.24082E-06_JPRB, 2.18577E-06_JPRB, 2.13207E-06_JPRB, 2.07969E-06_JPRB, &
-     & 2.02859E-06_JPRB, 1.97875E-06_JPRB, 1.93014E-06_JPRB, 1.88272E-06_JPRB/)
-      KBO_MN2(:,14) = (/ &
-     & 2.58051E-06_JPRB, 2.48749E-06_JPRB, 2.39782E-06_JPRB, 2.31139E-06_JPRB, 2.22807E-06_JPRB, &
-     & 2.14775E-06_JPRB, 2.07033E-06_JPRB, 1.99570E-06_JPRB, 1.92376E-06_JPRB, 1.85441E-06_JPRB, &
-     & 1.78756E-06_JPRB, 1.72313E-06_JPRB, 1.66101E-06_JPRB, 1.60114E-06_JPRB, 1.54342E-06_JPRB, &
-     & 1.48778E-06_JPRB, 1.43415E-06_JPRB, 1.38245E-06_JPRB, 1.33262E-06_JPRB/)
-      KBO_MN2(:,15) = (/ &
-     & 3.03447E-06_JPRB, 2.88559E-06_JPRB, 2.74401E-06_JPRB, 2.60938E-06_JPRB, 2.48135E-06_JPRB, &
-     & 2.35961E-06_JPRB, 2.24384E-06_JPRB, 2.13375E-06_JPRB, 2.02906E-06_JPRB, 1.92951E-06_JPRB, &
-     & 1.83484E-06_JPRB, 1.74481E-06_JPRB, 1.65921E-06_JPRB, 1.57780E-06_JPRB, 1.50039E-06_JPRB, &
-     & 1.42677E-06_JPRB, 1.35677E-06_JPRB, 1.29020E-06_JPRB, 1.22690E-06_JPRB/)
-      KBO_MN2(:,16) = (/ &
-     & 1.48655E-06_JPRB, 1.48283E-06_JPRB, 1.47913E-06_JPRB, 1.47543E-06_JPRB, 1.47174E-06_JPRB, &
-     & 1.46806E-06_JPRB, 1.46439E-06_JPRB, 1.46072E-06_JPRB, 1.45707E-06_JPRB, 1.45343E-06_JPRB, &
-     & 1.44979E-06_JPRB, 1.44617E-06_JPRB, 1.44255E-06_JPRB, 1.43894E-06_JPRB, 1.43534E-06_JPRB, &
-     & 1.43176E-06_JPRB, 1.42817E-06_JPRB, 1.42460E-06_JPRB, 1.42104E-06_JPRB/)
-
-
-!     The array SELFREFO contains the coefficient of the water vapor
-!     self-continuum (including the energy term).  The first index
-!     refers to temperature in 7.2 degree increments.  For instance,
-!     JT = 1 refers to a temperature of 245.6, JT = 2 refers to 252.8,
-!     etc.  The second index runs over the g-channel (1 to 16).
-
-      SELFREFO(:, 1) = (/ &
-     & 2.16803e+00_JPRB, 1.98236e+00_JPRB, 1.81260e+00_JPRB, 1.65737e+00_JPRB, 1.51544e+00_JPRB, &
-     & 1.38567e+00_JPRB, 1.26700e+00_JPRB, 1.15850e+00_JPRB, 1.05929e+00_JPRB, 9.68576e-01_JPRB/)
-      SELFREFO(:, 2) = (/ &
-     & 3.70149e+00_JPRB, 3.43145e+00_JPRB, 3.18110e+00_JPRB, 2.94902e+00_JPRB, 2.73387e+00_JPRB, &
-     & 2.53441e+00_JPRB, 2.34951e+00_JPRB, 2.17810e+00_JPRB, 2.01919e+00_JPRB, 1.87188e+00_JPRB/)
-      SELFREFO(:, 3) = (/ &
-     & 6.17433e+00_JPRB, 5.62207e+00_JPRB, 5.11920e+00_JPRB, 4.66131e+00_JPRB, 4.24438e+00_JPRB, &
-     & 3.86474e+00_JPRB, 3.51906e+00_JPRB, 3.20430e+00_JPRB, 2.91769e+00_JPRB, 2.65672e+00_JPRB/)
-      SELFREFO(:, 4) = (/ &
-     & 6.56459e+00_JPRB, 5.94787e+00_JPRB, 5.38910e+00_JPRB, 4.88282e+00_JPRB, 4.42410e+00_JPRB, &
-     & 4.00848e+00_JPRB, 3.63190e+00_JPRB, 3.29070e+00_JPRB, 2.98155e+00_JPRB, 2.70145e+00_JPRB/)
-      SELFREFO(:, 5) = (/ &
-     & 6.49581e+00_JPRB, 5.91114e+00_JPRB, 5.37910e+00_JPRB, 4.89494e+00_JPRB, 4.45436e+00_JPRB, &
-     & 4.05344e+00_JPRB, 3.68860e+00_JPRB, 3.35660e+00_JPRB, 3.05448e+00_JPRB, 2.77956e+00_JPRB/)
-      SELFREFO(:, 6) = (/ &
-     & 6.50189e+00_JPRB, 5.89381e+00_JPRB, 5.34260e+00_JPRB, 4.84294e+00_JPRB, 4.39001e+00_JPRB, &
-     & 3.97944e+00_JPRB, 3.60727e+00_JPRB, 3.26990e+00_JPRB, 2.96409e+00_JPRB, 2.68687e+00_JPRB/)
-      SELFREFO(:, 7) = (/ &
-     & 6.64768e+00_JPRB, 6.01719e+00_JPRB, 5.44650e+00_JPRB, 4.92993e+00_JPRB, 4.46236e+00_JPRB, &
-     & 4.03914e+00_JPRB, 3.65605e+00_JPRB, 3.30930e+00_JPRB, 2.99543e+00_JPRB, 2.71134e+00_JPRB/)
-      SELFREFO(:, 8) = (/ &
-     & 6.43744e+00_JPRB, 5.87166e+00_JPRB, 5.35560e+00_JPRB, 4.88490e+00_JPRB, 4.45557e+00_JPRB, &
-     & 4.06397e+00_JPRB, 3.70679e+00_JPRB, 3.38100e+00_JPRB, 3.08384e+00_JPRB, 2.81281e+00_JPRB/)
-      SELFREFO(:, 9) = (/ &
-     & 6.55466e+00_JPRB, 5.99777e+00_JPRB, 5.48820e+00_JPRB, 5.02192e+00_JPRB, 4.59525e+00_JPRB, &
-     & 4.20484e+00_JPRB, 3.84759e+00_JPRB, 3.52070e+00_JPRB, 3.22158e+00_JPRB, 2.94787e+00_JPRB/)
-      SELFREFO(:,10) = (/ &
-     & 6.84510e+00_JPRB, 6.26933e+00_JPRB, 5.74200e+00_JPRB, 5.25902e+00_JPRB, 4.81667e+00_JPRB, &
-     & 4.41152e+00_JPRB, 4.04046e+00_JPRB, 3.70060e+00_JPRB, 3.38933e+00_JPRB, 3.10424e+00_JPRB/)
-      SELFREFO(:,11) = (/ &
-     & 6.83128e+00_JPRB, 6.25536e+00_JPRB, 5.72800e+00_JPRB, 5.24510e+00_JPRB, 4.80291e+00_JPRB, &
-     & 4.39799e+00_JPRB, 4.02722e+00_JPRB, 3.68770e+00_JPRB, 3.37681e+00_JPRB, 3.09212e+00_JPRB/)
-      SELFREFO(:,12) = (/ &
-     & 7.35969e+00_JPRB, 6.61719e+00_JPRB, 5.94960e+00_JPRB, 5.34936e+00_JPRB, 4.80968e+00_JPRB, &
-     & 4.32445e+00_JPRB, 3.88817e+00_JPRB, 3.49590e+00_JPRB, 3.14321e+00_JPRB, 2.82610e+00_JPRB/)
-      SELFREFO(:,13) = (/ &
-     & 7.50064e+00_JPRB, 6.80749e+00_JPRB, 6.17840e+00_JPRB, 5.60744e+00_JPRB, 5.08925e+00_JPRB, &
-     & 4.61894e+00_JPRB, 4.19210e+00_JPRB, 3.80470e+00_JPRB, 3.45310e+00_JPRB, 3.13399e+00_JPRB/)
-      SELFREFO(:,14) = (/ &
-     & 7.40801e+00_JPRB, 6.71328e+00_JPRB, 6.08370e+00_JPRB, 5.51316e+00_JPRB, 4.99613e+00_JPRB, &
-     & 4.52759e+00_JPRB, 4.10298e+00_JPRB, 3.71820e+00_JPRB, 3.36950e+00_JPRB, 3.05351e+00_JPRB/)
-      SELFREFO(:,15) = (/ &
-     & 7.51895e+00_JPRB, 6.68846e+00_JPRB, 5.94970e+00_JPRB, 5.29254e+00_JPRB, 4.70796e+00_JPRB, &
-     & 4.18795e+00_JPRB, 3.72538e+00_JPRB, 3.31390e+00_JPRB, 2.94787e+00_JPRB, 2.62227e+00_JPRB/)
-      SELFREFO(:,16) = (/ &
-     & 7.84774e+00_JPRB, 6.80673e+00_JPRB, 5.90380e+00_JPRB, 5.12065e+00_JPRB, 4.44138e+00_JPRB, &
-     & 3.85223e+00_JPRB, 3.34122e+00_JPRB, 2.89800e+00_JPRB, 2.51357e+00_JPRB, 2.18014e+00_JPRB/)
-
-
-
-
-
-IF (LHOOK) CALL DR_HOOK('RRTM_KGB1',1,ZHOOK_HANDLE)
-RETURN
-
-1000 CONTINUE
-CALL ABOR1("RRTM_KGB1:ERROR OPENING FILE RADRRTM")
-1001 CONTINUE
-CALL ABOR1("RRTM_KGB1:ERROR READING FILE RADRRTM")
-
-!     -----------------------------------------------------------------
-END SUBROUTINE RRTM_KGB1
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_gas_optical_depth.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_gas_optical_depth.F90
deleted file mode 100644
index 693f0ea06..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_gas_optical_depth.F90
+++ /dev/null
@@ -1,339 +0,0 @@
-#ifdef RS6K
-@PROCESS HOT(NOVECTOR) NOSTRICT
-#endif
-SUBROUTINE SRTM_GAS_OPTICAL_DEPTH &
- & ( KIDIA   , KFDIA   , KLEV    , PONEMINUS, &
- &   PRMU0, &
- &   KLAYTROP,&
- &   PCOLCH4  , PCOLCO2 , PCOLH2O , PCOLMOL  , PCOLO2 , PCOLO3 ,&
- &   PFORFAC , PFORFRAC , KINDFOR , PSELFFAC, PSELFFRAC, KINDSELF ,&
- &   PFAC00  , PFAC01   , PFAC10  , PFAC11 ,&
- &   KJP     , KJT      , KJT1 ,&
- !-- output arrays 
- &   POD, PSSA, PINCSOL)
-
-
-!**** *SRTM_GAS_OPTICAL_DEPTH* - SPECTRAL LOOP TO COMPUTE THE SHORTWAVE RADIATION FLUXES.
-
-!     PURPOSE.
-!     --------
-
-!          COMPUTE THE GAS OPTICAL DEPTH AT EACH SHORTWAVE G POINT
-
-!**   INTERFACE.
-!     ----------
-
-!          *SRTM_GAS_OPTICAL_DEPTH* IS CALLED FROM THE NEW RADIATION SCHEME
-
-!        IMPLICIT ARGUMENTS :
-!        --------------------
-
-!     ==== INPUTS ===
-!     ==== OUTPUTS ===
-
-!     METHOD.
-!     -------
-
-!     EXTERNALS.
-!     ----------
-
-!     REFERENCE.
-!     ----------
-
-!        SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
-!        DOCUMENTATION
-!     AUTHOR.
-!     -------
-!        ADAPTED FROM SRTM_SPCVRT_MCICA (BY JEAN-JACQUES MORCRETTE) BY
-!        ROBIN HOGAN
-!
-!     MODIFICATIONS.
-!     --------------
-!        ORIGINAL : 2015-07-16
-
-!     ------------------------------------------------------------------
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
-USE PARSRTM  , ONLY : JPB1, JPB2
-USE YOESRTM  , ONLY : JPGPT
-USE YOESRTWN , ONLY : NGC
-USE YOMDIMV  , ONLY : YRDIMV
-
-IMPLICIT NONE
-
-!     ------------------------------------------------------------------
-
-!*       0.1   ARGUMENTS
-!              ---------
-
-INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA, KFDIA
-INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PONEMINUS(KIDIA:KFDIA)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMU0(KIDIA:KFDIA)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KLAYTROP(KIDIA:KFDIA)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLCH4(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLCO2(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLH2O(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLMOL(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLO2(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLO3(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFORFAC(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFORFRAC(KIDIA:KFDIA,YRDIMV%NFLEVG)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KINDFOR(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PSELFFAC(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PSELFFRAC(KIDIA:KFDIA,YRDIMV%NFLEVG)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KINDSELF(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFAC00(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFAC01(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFAC10(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB)   ,INTENT(IN)    :: PFAC11(KIDIA:KFDIA,YRDIMV%NFLEVG)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KJP(KIDIA:KFDIA,YRDIMV%NFLEVG)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KJT(KIDIA:KFDIA,YRDIMV%NFLEVG)
-INTEGER(KIND=JPIM),INTENT(IN)    :: KJT1(KIDIA:KFDIA,YRDIMV%NFLEVG)
-
-REAL(KIND=JPRB)   ,INTENT(OUT)   :: POD(KIDIA:KFDIA,YRDIMV%NFLEVG,JPGPT) ! Optical depth
-REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSSA(KIDIA:KFDIA,YRDIMV%NFLEVG,JPGPT) ! Single scattering albedo
-REAL(KIND=JPRB)   ,INTENT(OUT)   :: PINCSOL(KIDIA:KFDIA,JPGPT) ! Incoming solar flux
-
-
-!     ------------------------------------------------------------------
-
-INTEGER(KIND=JPIM) :: IB1, IB2, IBM, IGT, IW(KIDIA:KFDIA), JB, JG, JK, JL, IC, ICOUNT
-
-INTEGER(KIND=JPIM) :: IND(KFDIA-KIDIA+1)
-
-
-!-- Output of SRTM_TAUMOLn routines
-REAL(KIND=JPRB) :: ZTAUG(KIDIA:KFDIA,YRDIMV%NFLEVG,16) ! Absorption optical depth
-REAL(KIND=JPRB) :: ZTAUR(KIDIA:KFDIA,YRDIMV%NFLEVG,16) ! Rayleigh optical depth
-REAL(KIND=JPRB) :: ZSFLXZEN(KIDIA:KFDIA,16) ! Incoming solar flux
-
-
-REAL(KIND=JPRB) :: ZTAU, ZPAO, ZPTO
-REAL(KIND=JPRB) :: ZPAOJ(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB) :: ZPTOJ(KIDIA:KFDIA,YRDIMV%NFLEVG)
-REAL(KIND=JPRB) :: ZRMU0D(KIDIA:KFDIA) 
- 
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-
-#include "srtm_taumol16.intfb.h"
-#include "srtm_taumol17.intfb.h"
-#include "srtm_taumol18.intfb.h"
-#include "srtm_taumol19.intfb.h"
-#include "srtm_taumol20.intfb.h"
-#include "srtm_taumol21.intfb.h"
-#include "srtm_taumol22.intfb.h"
-#include "srtm_taumol23.intfb.h"
-#include "srtm_taumol24.intfb.h"
-#include "srtm_taumol25.intfb.h"
-#include "srtm_taumol26.intfb.h"
-#include "srtm_taumol27.intfb.h"
-#include "srtm_taumol28.intfb.h"
-#include "srtm_taumol29.intfb.h"
-
-!     ------------------------------------------------------------------
-ASSOCIATE(NFLEVG=>YRDIMV%NFLEVG)
-IF (LHOOK) CALL DR_HOOK('SRTM_GAS_OPTICAL_DEPTH',0,ZHOOK_HANDLE)
-
-POD = 0.0	
-PSSA = 0.0
-PINCSOL = 0.0
-
-IB1=JPB1
-IB2=JPB2
-
-IC=0
-DO JL = KIDIA, KFDIA
-  IF (PRMU0(JL) > 0.0_JPRB) THEN
-    IC=IC+1
-    IND(IC)=JL
-    IW(JL)=0
-  ENDIF
-ENDDO
-ICOUNT=IC
-IF(ICOUNT==0)THEN
-  IF (LHOOK) CALL DR_HOOK('SRTM_SPCVRT_MCICA',1,ZHOOK_HANDLE)
-  RETURN
-ENDIF
-
-JB=IB1-1
-DO JB = IB1, IB2
-  DO IC=1,ICOUNT
-    JL=IND(IC)
-    IBM = JB-15
-    IGT = NGC(IBM)
-  ENDDO
-
-  !-- for each band, computes the gaseous and Rayleigh optical thickness 
-  !  for all g-points within the band
-
-  IF (JB == 16) THEN
-    CALL SRTM_TAUMOL16 &
-     & ( KIDIA   , KFDIA    , KLEV    ,&
-     &   PFAC00  , PFAC01   , PFAC10   , PFAC11   ,&
-     &   KJP     , KJT      , KJT1     , PONEMINUS,&
-     &   PCOLH2O , PCOLCH4  , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC , PSELFFRAC, KINDSELF, PFORFAC  , PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG    , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 17) THEN
-    CALL SRTM_TAUMOL17 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLCO2 , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 18) THEN
-    CALL SRTM_TAUMOL18 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLCH4 , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 19) THEN
-    CALL SRTM_TAUMOL19 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLCO2 , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 20) THEN
-    CALL SRTM_TAUMOL20 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     ,&
-     &   PCOLH2O , PCOLCH4 , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 21) THEN
-    CALL SRTM_TAUMOL21 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLCO2 , PCOLMOL  ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 22) THEN
-    CALL SRTM_TAUMOL22 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLMOL , PCOLO2   ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 23) THEN
-    CALL SRTM_TAUMOL23 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     ,&
-     &   PCOLH2O , PCOLMOL ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 24) THEN
-    CALL SRTM_TAUMOL24 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     , PONEMINUS ,&
-     &   PCOLH2O , PCOLMOL , PCOLO2   , PCOLO3 ,&
-     &   KLAYTROP, PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 25) THEN
-    !--- visible 16000-22650 cm-1   0.4415 - 0.6250 um
-    CALL SRTM_TAUMOL25 &
-     & ( KIDIA    , KFDIA   , KLEV     ,&
-     &   PFAC00   , PFAC01  , PFAC10 , PFAC11 ,&
-     &   KJP      , KJT     , KJT1   ,&
-     &   PCOLH2O  , PCOLMOL , PCOLO3 ,&
-     &   KLAYTROP ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR   , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 26) THEN
-    !--- UV-A 22650-29000 cm-1   0.3448 - 0.4415 um
-    CALL SRTM_TAUMOL26 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PCOLMOL ,KLAYTROP,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 27) THEN
-    !--- UV-B 29000-38000 cm-1   0.2632 - 0.3448 um
-    CALL SRTM_TAUMOL27 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP     , KJT     , KJT1     ,&
-     &   PCOLMOL , PCOLO3 ,&
-     &   KLAYTROP ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 28) THEN
-    !--- UV-C 38000-50000 cm-1   0.2000 - 0.2632 um
-    CALL SRTM_TAUMOL28 &
-     & ( KIDIA   , KFDIA   , KLEV    ,&
-     &   PFAC00  , PFAC01  , PFAC10 , PFAC11 ,&
-     &   KJP     , KJT     , KJT1   , PONEMINUS ,&
-     &   PCOLMOL , PCOLO2  , PCOLO3 ,&
-     &   KLAYTROP ,&
-     &   ZSFLXZEN, ZTAUG   , ZTAUR  , PRMU0     &
-     & )  
-
-  ELSEIF (JB == 29) THEN
-    CALL SRTM_TAUMOL29 &
-     & ( KIDIA    , KFDIA   , KLEV     ,&
-     &   PFAC00   , PFAC01  , PFAC10   , PFAC11 ,&
-     &   KJP      , KJT     , KJT1     ,&
-     &   PCOLH2O  , PCOLCO2 , PCOLMOL  ,&
-     &   KLAYTROP , PSELFFAC, PSELFFRAC, KINDSELF  , PFORFAC, PFORFRAC, KINDFOR ,&
-     &   ZSFLXZEN , ZTAUG   , ZTAUR    , PRMU0     &
-     & )  
-
-  ENDIF
-   
-  DO JG=1,IGT
-    DO IC=1,ICOUNT
-      JL=IND(IC)
-      IW(JL)=IW(JL)+1
-
-      ! Incoming solar flux into plane perp to incoming radiation
-      PINCSOL(JL,IW(JL)) = ZSFLXZEN(JL,JG)
-    ENDDO
-
-    DO JK=1,KLEV
-      DO IC=1,ICOUNT
-        JL=IND(IC)
-        POD (JL,JK,IW(JL)) = ZTAUR(JL,JK,JG) + ZTAUG(JL,JK,JG)
-        PSSA(JL,JK,IW(JL)) = ZTAUR(JL,JK,JG) / POD(JL,JK,IW(JL))
-      ENDDO
-    ENDDO
-
-  ENDDO   !-- end loop on JG (g point)
-
-ENDDO     !-- end loop on JB (band)
-
-!     ------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('SRTM_GAS_OPTICAL_DEPTH',1,ZHOOK_HANDLE)
-END ASSOCIATE
-END SUBROUTINE SRTM_GAS_OPTICAL_DEPTH
diff --git a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_kgb16.F90 b/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_kgb16.F90
deleted file mode 100644
index c2295d630..000000000
--- a/src/LIB/RAD/ecrad-1.0.1_mnh/ifsrrtm/srtm_kgb16.F90
+++ /dev/null
@@ -1,184 +0,0 @@
-SUBROUTINE SRTM_KGB16(DIRECTORY)
-
-!     Originally by J.Delamere, Atmospheric & Environmental Research.
-!     Revision: 2.4
-!     BAND 16:  2600-3000 cm-1 (low - H2O,CH4; high - nothing)
-!     Reformatted for F90 by JJMorcrette, ECMWF
-!     R. Elkhatib 12-10-2005 Split for faster and more robust compilation.
-!     G.Mozdzynski March 2011 read constants from files
-!     T. Wilhelmsson and K. Yessad (Oct 2013) Geometry and setup refactoring.
-!     ------------------------------------------------------------------
-
-USE PARKIND1  , ONLY : JPRB
-USE YOMHOOK   , ONLY : LHOOK, DR_HOOK
-USE YOMLUN    , ONLY : NULRAD , NULOUT
-USE YOMMP0    , ONLY : NPROC, MYPROC
-USE MPL_MODULE, ONLY : MPL_BROADCAST
-USE YOMTAG    , ONLY : MTAGRAD
-USE YOESRTA16 , ONLY : KA, KB, SELFREF, FORREF, SFLUXREF, RAYL, STRRAT1, LAYREFFR
-
-!     ------------------------------------------------------------------
-
-IMPLICIT NONE
-
-CHARACTER(LEN=*), INTENT(IN) :: DIRECTORY
-
-! KURUCZ
-!CHARACTER(LEN = 80) :: CLZZZ
-CHARACTER(LEN = 80) :: CLF1
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-#include "abor1.intfb.h"
-
-IF (LHOOK) CALL DR_HOOK('SRTM_KGB16',0,ZHOOK_HANDLE)
-
-IF( MYPROC==1 )THEN
-  !CALL GETENV("DATA",CLZZZ)
-  !IF(CLZZZ /= " ") THEN
-  !  CLF1=TRIM(CLZZZ)//"/RADSRTM"
-  CLF1 = DIRECTORY // "/RADSRTM"
-  WRITE(NULOUT,'(A,A)') 'Reading ',TRIM(CLF1)
-  ! RRTM and SRTM files from ecrad are in big-endian format.
-  ! Here they are covnerted into little-endian at opening
-  ! No need for compialtion option export GFORTRAN_CONVERT_UNIT="little_endian;big_endian:145"
-!   OPEN(NULRAD,FILE=TRIM(CLF1),FORM="UNFORMATTED",ACTION="READ",access='sequential',ERR=1000,CONVERT='swap')
-
-  OPEN(NULRAD,FILE=TRIM(CLF1),FORM="UNFORMATTED",ACTION="READ",access='sequential',ERR=1000,CONVERT='BIG_ENDIAN')
-  !ELSE
-  !  OPEN(NULRAD,FILE='RADSRTM',FORM="UNFORMATTED",ACTION="READ",ERR=1000)
-  !ENDIF
-  READ(NULRAD,ERR=1001) KA,KB
-ENDIF
-IF( NPROC>1 )THEN
-  CALL MPL_BROADCAST (KA,MTAGRAD,1,CDSTRING='SRTM_KGB16:')
-  CALL MPL_BROADCAST (KB,MTAGRAD,1,CDSTRING='SRTM_KGB16:')
-ENDIF
-
-SFLUXREF = (/ &
- & 1.92269_JPRB    , 1.72844_JPRB    , 1.64326_JPRB    , 1.58451_JPRB     &
- & , 1.44031_JPRB    , 1.25108_JPRB    , 1.02724_JPRB    , 0.776759_JPRB    &
- & , 0.534444_JPRB   , 5.87755E-02_JPRB, 4.86706E-02_JPRB, 3.87989E-02_JPRB &
- & , 2.84532E-02_JPRB, 1.82431E-02_JPRB, 6.92320E-03_JPRB, 9.70770E-04_JPRB /)  
-
-!     Rayleigh extinction coefficient at v = 2925 cm-1.
-RAYL = 2.91E-10_JPRB
-
-STRRAT1 = 252.131_JPRB
-
-LAYREFFR = 18
-
-!     ------------------------------------------------------------------
-
-!     The array KA contains absorption coefs at the 16 chosen g-values 
-!     for a range of pressure levels> ~100mb, temperatures, and binary
-!     species parameters (see taumol.f for definition).  The first 
-!     index in the array, JS, runs from 1 to 9, and corresponds to 
-!     different values of the binary species parameter.  For instance, 
-!     JS=1 refers to dry air, JS = 2 corresponds to the paramter value 1/8, 
-!     JS = 3 corresponds to the parameter value 2/8, etc.  The second index
-!     in the array, JT, which runs from 1 to 5, corresponds to different
-!     temperatures.  More specifically, JT = 3 means that the data are for
-!     the reference temperature TREF for this  pressure level, JT = 2 refers
-!     to TREF-15, JT = 1 is for TREF-30, JT = 4 is for TREF+15, and JT = 5
-!     is for TREF+30.  The third index, JP, runs from 1 to 13 and refers
-!     to the JPth reference pressure level (see taumol.f for these levels
-!     in mb).  The fourth index, IG, goes from 1 to 16, and indicates
-!     which g-interval the absorption coefficients are for.
-!     -----------------------------------------------------------------
-
-!     -----------------------------------------------------------------
-!     The array KB contains absorption coefs at the 16 chosen g-values 
-!     for a range of pressure levels < ~100mb and temperatures. The first 
-!     index in the array, JT, which runs from 1 to 5, corresponds to 
-!     different temperatures.  More specifically, JT = 3 means that the 
-!     data are for the reference temperature TREF for this pressure 
-!     level, JT = 2 refers to the temperature TREF-15, JT = 1 is for
-!     TREF-30, JT = 4 is for TREF+15, and JT = 5 is for TREF+30.  
-!     The second index, JP, runs from 13 to 59 and refers to the JPth
-!     reference pressure level (see taumol.f for the value of these
-!     pressure levels in mb).  The third index, IG, goes from 1 to 16,
-!     and tells us which g-interval the absorption coefficients are for.
-!     -----------------------------------------------------------------
-
-FORREF(:, 1) = (/ 0.525585E-05_JPRB, 0.527618E-05_JPRB, 0.746929E-04_JPRB /)
-FORREF(:, 2) = (/ 0.794660E-05_JPRB, 0.136902E-04_JPRB, 0.849878E-04_JPRB /)
-FORREF(:, 3) = (/ 0.197099E-04_JPRB, 0.733094E-04_JPRB, 0.121687E-03_JPRB /)
-FORREF(:, 4) = (/ 0.148274E-03_JPRB, 0.169776E-03_JPRB, 0.164848E-03_JPRB /)
-FORREF(:, 5) = (/ 0.230296E-03_JPRB, 0.210384E-03_JPRB, 0.182028E-03_JPRB /)
-FORREF(:, 6) = (/ 0.280575E-03_JPRB, 0.259217E-03_JPRB, 0.196080E-03_JPRB /)
-FORREF(:, 7) = (/ 0.329034E-03_JPRB, 0.291575E-03_JPRB, 0.207044E-03_JPRB /)
-FORREF(:, 8) = (/ 0.349989E-03_JPRB, 0.323471E-03_JPRB, 0.225712E-03_JPRB /)
-FORREF(:, 9) = (/ 0.366097E-03_JPRB, 0.321519E-03_JPRB, 0.253150E-03_JPRB /)
-FORREF(:,10) = (/ 0.383589E-03_JPRB, 0.355314E-03_JPRB, 0.262555E-03_JPRB /)
-FORREF(:,11) = (/ 0.375933E-03_JPRB, 0.372443E-03_JPRB, 0.261313E-03_JPRB /)
-FORREF(:,12) = (/ 0.370652E-03_JPRB, 0.382366E-03_JPRB, 0.250070E-03_JPRB /)
-FORREF(:,13) = (/ 0.375092E-03_JPRB, 0.379542E-03_JPRB, 0.265794E-03_JPRB /)
-FORREF(:,14) = (/ 0.389705E-03_JPRB, 0.384274E-03_JPRB, 0.322135E-03_JPRB /)
-FORREF(:,15) = (/ 0.372084E-03_JPRB, 0.390422E-03_JPRB, 0.370035E-03_JPRB /)
-FORREF(:,16) = (/ 0.437802E-03_JPRB, 0.373406E-03_JPRB, 0.373222E-03_JPRB /)
-
-!     -----------------------------------------------------------------
-!     The array SELFREF contains the coefficient of the water vapor
-!     self-continuum (including the energy term).  The first index
-!     refers to temperature in 7.2 degree increments.  For instance,
-!     JT = 1 refers to a temperature of 245.6, JT = 2 refers to 252.8,
-!     etc.  The second index runs over the g-channel (1 to 16).
-
-SELFREF(:, 1) = (/ &
- & 0.126758E-02_JPRB, 0.105253E-02_JPRB, 0.873963E-03_JPRB, 0.725690E-03_JPRB, 0.602573E-03_JPRB, &
- & 0.500344E-03_JPRB, 0.415458E-03_JPRB, 0.344973E-03_JPRB, 0.286447E-03_JPRB, 0.237849E-03_JPRB /)  
-SELFREF(:, 2) = (/ &
- & 0.144006E-02_JPRB, 0.118514E-02_JPRB, 0.975351E-03_JPRB, 0.802697E-03_JPRB, 0.660606E-03_JPRB, &
- & 0.543667E-03_JPRB, 0.447429E-03_JPRB, 0.368226E-03_JPRB, 0.303044E-03_JPRB, 0.249400E-03_JPRB /)  
-SELFREF(:, 3) = (/ &
- & 0.294018E-02_JPRB, 0.227428E-02_JPRB, 0.175920E-02_JPRB, 0.136077E-02_JPRB, 0.105258E-02_JPRB, &
- & 0.814189E-03_JPRB, 0.629789E-03_JPRB, 0.487153E-03_JPRB, 0.376821E-03_JPRB, 0.291478E-03_JPRB /)  
-SELFREF(:, 4) = (/ &
- & 0.395290E-02_JPRB, 0.348405E-02_JPRB, 0.307081E-02_JPRB, 0.270658E-02_JPRB, 0.238556E-02_JPRB, &
- & 0.210261E-02_JPRB, 0.185322E-02_JPRB, 0.163341E-02_JPRB, 0.143967E-02_JPRB, 0.126891E-02_JPRB /)  
-SELFREF(:, 5) = (/ &
- & 0.419122E-02_JPRB, 0.385638E-02_JPRB, 0.354829E-02_JPRB, 0.326481E-02_JPRB, 0.300398E-02_JPRB, &
- & 0.276399E-02_JPRB, 0.254317E-02_JPRB, 0.234000E-02_JPRB, 0.215305E-02_JPRB, 0.198104E-02_JPRB /)  
-SELFREF(:, 6) = (/ &
- & 0.495659E-02_JPRB, 0.456777E-02_JPRB, 0.420945E-02_JPRB, 0.387924E-02_JPRB, 0.357494E-02_JPRB, &
- & 0.329450E-02_JPRB, 0.303606E-02_JPRB, 0.279790E-02_JPRB, 0.257842E-02_JPRB, 0.237615E-02_JPRB /)  
-SELFREF(:, 7) = (/ &
- & 0.526981E-02_JPRB, 0.490687E-02_JPRB, 0.456893E-02_JPRB, 0.425426E-02_JPRB, 0.396126E-02_JPRB, &
- & 0.368844E-02_JPRB, 0.343441E-02_JPRB, 0.319788E-02_JPRB, 0.297764E-02_JPRB, 0.277256E-02_JPRB /)  
-SELFREF(:, 8) = (/ &
- & 0.575426E-02_JPRB, 0.531597E-02_JPRB, 0.491106E-02_JPRB, 0.453699E-02_JPRB, 0.419141E-02_JPRB, &
- & 0.387216E-02_JPRB, 0.357722E-02_JPRB, 0.330475E-02_JPRB, 0.305303E-02_JPRB, 0.282048E-02_JPRB /)  
-SELFREF(:, 9) = (/ &
- & 0.549881E-02_JPRB, 0.514328E-02_JPRB, 0.481074E-02_JPRB, 0.449970E-02_JPRB, 0.420877E-02_JPRB, &
- & 0.393665E-02_JPRB, 0.368213E-02_JPRB, 0.344406E-02_JPRB, 0.322138E-02_JPRB, 0.301310E-02_JPRB /)  
-SELFREF(:,10) = (/ &
- & 0.605357E-02_JPRB, 0.561246E-02_JPRB, 0.520349E-02_JPRB, 0.482432E-02_JPRB, 0.447278E-02_JPRB, &
- & 0.414686E-02_JPRB, 0.384469E-02_JPRB, 0.356453E-02_JPRB, 0.330479E-02_JPRB, 0.306398E-02_JPRB /)  
-SELFREF(:,11) = (/ &
- & 0.640504E-02_JPRB, 0.587858E-02_JPRB, 0.539540E-02_JPRB, 0.495194E-02_JPRB, 0.454492E-02_JPRB, &
- & 0.417136E-02_JPRB, 0.382850E-02_JPRB, 0.351382E-02_JPRB, 0.322501E-02_JPRB, 0.295993E-02_JPRB /)  
-SELFREF(:,12) = (/ &
- & 0.677803E-02_JPRB, 0.615625E-02_JPRB, 0.559152E-02_JPRB, 0.507859E-02_JPRB, 0.461271E-02_JPRB, &
- & 0.418957E-02_JPRB, 0.380524E-02_JPRB, 0.345617E-02_JPRB, 0.313913E-02_JPRB, 0.285116E-02_JPRB /)  
-SELFREF(:,13) = (/ &
- & 0.690347E-02_JPRB, 0.627003E-02_JPRB, 0.569472E-02_JPRB, 0.517219E-02_JPRB, 0.469761E-02_JPRB, &
- & 0.426658E-02_JPRB, 0.387509E-02_JPRB, 0.351953E-02_JPRB, 0.319659E-02_JPRB, 0.290328E-02_JPRB /)  
-SELFREF(:,14) = (/ &
- & 0.692680E-02_JPRB, 0.632795E-02_JPRB, 0.578087E-02_JPRB, 0.528109E-02_JPRB, 0.482452E-02_JPRB, &
- & 0.440742E-02_JPRB, 0.402638E-02_JPRB, 0.367828E-02_JPRB, 0.336028E-02_JPRB, 0.306977E-02_JPRB /)  
-SELFREF(:,15) = (/ &
- & 0.754894E-02_JPRB, 0.681481E-02_JPRB, 0.615207E-02_JPRB, 0.555378E-02_JPRB, 0.501367E-02_JPRB, &
- & 0.452609E-02_JPRB, 0.408593E-02_JPRB, 0.368857E-02_JPRB, 0.332986E-02_JPRB, 0.300603E-02_JPRB /)  
-SELFREF(:,16) = (/ &
- & 0.760689E-02_JPRB, 0.709755E-02_JPRB, 0.662232E-02_JPRB, 0.617891E-02_JPRB, 0.576519E-02_JPRB, &
- & 0.537917E-02_JPRB, 0.501899E-02_JPRB, 0.468293E-02_JPRB, 0.436938E-02_JPRB, 0.407682E-02_JPRB /)  
-
-IF (LHOOK) CALL DR_HOOK('SRTM_KGB16',1,ZHOOK_HANDLE)
-RETURN
-
-1000 CONTINUE
-CALL ABOR1("SRTM_KGB16:ERROR OPENING FILE RADSRTM")
-1001 CONTINUE
-CALL ABOR1("SRTM_KGB16:ERROR READING FILE RADSRTM")
-
-END SUBROUTINE SRTM_KGB16
-- 
GitLab