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 a0c35ec7763081afce3ccaaa91ade2576f2f8193..0000000000000000000000000000000000000000
--- 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 215743cfb0e67a32125a895a48f7114e5099eba5..0000000000000000000000000000000000000000
--- 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 8d078cca1e9be027612fc25162d4c1b1ff36ea82..0000000000000000000000000000000000000000
--- 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 7a1348ad92996ce8ca3a900ed37b7914994ac00f..0000000000000000000000000000000000000000
--- 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 e0cdd612165c9401974af01bd5de8a1092087509..0000000000000000000000000000000000000000
--- 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 045d7a395ce45b36020e410d3542098bb676afe3..0000000000000000000000000000000000000000
--- 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 a26a5740d51a439791235810147b420531d1d4a7..0000000000000000000000000000000000000000
--- 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 f8758823086cffa742e46871378dbd83542dba5d..0000000000000000000000000000000000000000
--- 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 e2b5bca94223aad523ff94b3fc99d66334a36dfe..0000000000000000000000000000000000000000
--- 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 693f0ea061c906961b09fd259623803621febb0c..0000000000000000000000000000000000000000
--- 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 c2295d630cb84690228df926cbe42accfbee4916..0000000000000000000000000000000000000000
--- 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
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/mpl_abort_mod.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/mpl_abort_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..9c5f7e1a9e21332c1abd325fcc243e11945054f6
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/mpl_abort_mod.F90
@@ -0,0 +1,68 @@
+! (C) Copyright 2014- ECMWF.
+!
+! This software is licensed under the terms of the Apache Licence Version 2.0
+! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation
+! nor does it submit to any jurisdiction.
+
+MODULE MPL_ABORT_MOD
+
+USE MPL_DATA_MODULE
+USE MPL_MPIF
+USE OML_MOD
+USE YOMABRT, ONLY : MAB_CNT
+USE SDL_MOD, ONLY : SDL_TRACEBACK, SDL_DISABORT
+#ifdef NAG
+USE F90_UNIX_IO, ONLY: FLUSH
+#endif
+USE PARKIND1 ,ONLY : JPIM ,JPRB
+USE YOMHOOK , ONLY : LHOOK
+
+PRIVATE
+PUBLIC MPL_ABORT
+
+CONTAINS
+
+SUBROUTINE MPL_ABORT(CDMESSAGE)
+
+IMPLICIT NONE
+
+CHARACTER*(*),INTENT(IN),OPTIONAL :: CDMESSAGE
+INTEGER(KIND=JPIM) :: IRETURN_CODE,IERROR,ITID,INUMTH
+
+CHARACTER(LEN=80) :: CLTRBK
+
+ITID=OML_MY_THREAD()
+INUMTH=OML_MAX_THREADS()
+
+IF (MPL_UNIT > 0) CALL FLUSH(MPL_UNIT)
+!------Traceback from only one thread
+!$OMP CRITICAL (CRIT_MPL_ABORT)
+!$OMP FLUSH(MAB_CNT)
+IF (MAB_CNT == 0) THEN
+ WRITE(MPL_ERRUNIT,'(A,I6,A,I6)') 'MPL_ABORT: CALLED FROM PROCESSOR ',MPL_RANK,' THRD',ITID
+ IF(PRESENT(CDMESSAGE)) THEN
+ WRITE(MPL_ERRUNIT,*) 'MPL_ABORT: THRD',ITID," ",CDMESSAGE
+ ENDIF
+ MAB_CNT=1
+!$OMP FLUSH(MAB_CNT)
+
+#if defined(__INTEL_COMPILER)
+ CALL GET_ENVIRONMENT_VARIABLE("EC_LINUX_TRBK",CLTRBK)
+#else
+ CLTRBK='1'
+#endif
+ IF (LHOOK .AND. CLTRBK=='1') THEN
+ !CALL TABORT() ! should not hang and calls DrHook's error traceback processing (more robust nowadays)
+ ELSE
+ !CALL SDL_TRACEBACK(ITID) ! this will no longer hang with Intel compiler because intel tracebackqq is called, not linux traceback
+ ENDIF
+ENDIF
+!$OMP END CRITICAL (CRIT_MPL_ABORT)
+CALL SDL_DISABORT(MPL_COMM_OML(ITID))
+
+END SUBROUTINE MPL_ABORT
+
+END MODULE MPL_ABORT_MOD
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/sdl_mod.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/sdl_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..0b3fb317b8e696a092080564a2d8157300ae9bfd
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/drhook/drhook/sdl_mod.F90
@@ -0,0 +1,172 @@
+! (C) Copyright 2014- ECMWF.
+!
+! This software is licensed under the terms of the Apache Licence Version 2.0
+! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation
+! nor does it submit to any jurisdiction.
+
+MODULE SDL_MOD
+
+! Interface between user applications and system-dependent intrinsic
+! routines, provided by the computer vendors.
+
+! All routines which wish to call these routines must contain:
+! USE SDL_MOD
+
+! Author :
+! ------
+! 11-Apr-2005 R. El Khatib *METEO-FRANCE*
+! 26-Apr-2006 S.T.Saarinen Dr.Hook trace, calls to EC_RAISE, Intel/ifort traceback
+
+USE PARKIND1 ,ONLY : JPIM ,JPRB
+!USE YOMHOOK ,ONLY : LHOOK ,DR_HOOK
+USE YOMHOOK ,ONLY : LHOOK
+USE OML_MOD ,ONLY : OML_MY_THREAD
+USE MPL_MPIF ,ONLY : MPI_COMM_WORLD
+IMPLICIT NONE
+
+SAVE
+
+PRIVATE
+
+INTEGER, PARAMETER :: SIGABRT = 6 ! Hardcoded
+
+PUBLIC :: SDL_SRLABORT, SDL_DISABORT, SDL_TRACEBACK
+
+CONTAINS
+
+!-----------------------------------------------------------------------------
+SUBROUTINE SDL_TRACEBACK(KTID)
+
+! Purpose :
+! -------
+! Traceback
+
+! KTID : thread
+
+INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL :: KTID
+INTEGER(KIND=JPIM) ITID, IPRINT_OPTION, ILEVEL
+CHARACTER(LEN=80) :: CLTRBK
+#ifdef NECSX
+CHARACTER(LEN=*), PARAMETER :: CLNECMSG = '*** Calling NEC traceback ***'
+#endif
+IF (PRESENT(KTID)) THEN
+ ITID = KTID
+ELSE
+ ITID = OML_MY_THREAD()
+ENDIF
+IF (LHOOK) THEN
+ IPRINT_OPTION = 2
+ ILEVEL = 0
+ !CALL C_DRHOOK_PRINT(0, ITID, IPRINT_OPTION, ILEVEL) ! from drhook.c
+ENDIF
+#if defined(VPP)
+ CALL ERRTRA
+ IF (PRESENT(KTID)) CALL SLEEP(28)
+#elif defined(RS6K)
+ WRITE(0,*)'SDL_TRACEBACK: Calling XL_TRBK, THRD = ',ITID
+ CALL XL__TRBK()
+ WRITE(0,*)'SDL_TRACEBACK: Done XL_TRBK, THRD = ',ITID
+#elif defined(__INTEL_COMPILER)
+
+ CALL GET_ENVIRONMENT_VARIABLE("EC_LINUX_TRBK",CLTRBK)
+ IF (CLTRBK=='1') THEN
+ WRITE(0,*)'SDL_TRACEBACK: Calling LINUX_TRBK, THRD = ',ITID
+ !CALL LINUX_TRBK() ! See ifsaux/utilities/linuxtrbk.c
+ WRITE(0,*)'SDL_TRACEBACK: Done LINUX_TRBK, THRD = ',ITID
+ ELSE
+ WRITE(0,*)'SDL_TRACEBACK: Calling INTEL_TRBK, THRD = ',ITID
+ CALL INTEL_TRBK() ! See ifsaux/utilities/gentrbk.F90
+ WRITE(0,*)'SDL_TRACEBACK: Done INTEL_TRBK, THRD = ',ITID
+ ENDIF
+#elif defined(__NEC__)
+ ! A traceback using gdb-debugger, if available AND
+ ! activated via 'export GDBDEBUGGER=1'
+ WRITE(0,*)'SDL_TRACEBACK: Calling GDB_TRBK, THRD = ',ITID
+ CALL GDB_TRBK() ! See ifsaux/utilities/linuxtrbk.c
+ WRITE(0,*)'SDL_TRACEBACK: Done GDB_TRBK, THRD = ',ITID
+#elif defined(LINUX) || defined(SUN4)
+ WRITE(0,*)'SDL_TRACEBACK: Calling LINUX_TRBK, THRD = ',ITID
+ !CALL LINUX_TRBK() ! See ifsaux/utilities/linuxtrbk.c
+ WRITE(0,*)'SDL_TRACEBACK: Done LINUX_TRBK, THRD = ',ITID
+#elif defined(NECSX)
+! MESPUT writes out onto unit 6
+ WRITE(6,*)'SDL_TRACEBACK: Calling NEC/MESPUT, THRD = ',ITID
+ CALL NECSX_TRBK(CLNECMSG)
+ CALL FLUSH(6)
+ WRITE(6,*)'SDL_TRACEBACK: Done NEC/MESPUT, THRD = ',ITID
+#else
+ WRITE(0,*)'SDL_TRACEBACK: No proper traceback implemented.'
+ ! A traceback using dbx-debugger, if available AND
+ ! activated via 'export DBXDEBUGGER=1'
+ WRITE(0,*)'SDL_TRACEBACK: Calling DBX_TRBK, THRD = ',ITID
+ CALL DBX_TRBK() ! See ifsaux/utilities/linuxtrbk.c
+ WRITE(0,*)'SDL_TRACEBACK: Done DBX_TRBK, THRD = ',ITID
+ ! A traceback using gdb-debugger, if available AND
+ ! activated via 'export GDBDEBUGGER=1'
+ WRITE(0,*)'SDL_TRACEBACK: Calling GDB_TRBK, THRD = ',ITID
+ CALL GDB_TRBK() ! See ifsaux/utilities/linuxtrbk.c
+ WRITE(0,*)'SDL_TRACEBACK: Done GDB_TRBK, THRD = ',ITID
+#endif
+
+END SUBROUTINE SDL_TRACEBACK
+!-----------------------------------------------------------------------------
+SUBROUTINE SDL_SRLABORT
+
+! Purpose :
+! -------
+! To abort in serial environment
+
+!CALL EC_RAISE(SIGABRT)
+STOP 'SDL_SRLABORT'
+
+END SUBROUTINE SDL_SRLABORT
+!-----------------------------------------------------------------------------
+SUBROUTINE SDL_DISABORT(KCOMM)
+
+! Purpose :
+! -------
+! To abort in distributed environment
+
+! KCOMM : communicator
+
+INTEGER(KIND=JPIM), INTENT(IN) :: KCOMM
+
+INTEGER(KIND=JPIM) :: IRETURN_CODE,IERROR
+CHARACTER(LEN=80) :: CLJOBID
+CHARACTER(LEN=80) :: CLTRBK
+
+#ifdef VPP
+
+CALL VPP_ABORT()
+
+#else
+#if defined(__INTEL_COMPILER)
+! Intel compiler seems to hang in MPI_ABORT -- on all but the failing task(s)
+! ... when linux trbk is used. REK
+CALL GET_ENVIRONMENT_VARIABLE("EC_LINUX_TRBK",CLTRBK)
+IF (CLTRBK=='1') THEN
+IF (LHOOK) THEN
+ CALL GET_ENVIRONMENT_VARIABLE("SLURM_JOBID",CLJOBID)
+ IF (CLJOBID /= ' ') THEN
+ CALL SYSTEM("set -x; sleep 10; scancel --signal=TERM "//trim(CLJOBID)//" &")
+ ENDIF
+ENDIF
+ENDIF
+#endif
+
+IRETURN_CODE=SIGABRT
+!CALL MPI_ABORT(KCOMM,IRETURN_CODE,IERROR)
+CALL MPI_ABORT(MPI_COMM_WORLD,IRETURN_CODE,IERROR) ! Tracked by the supervisor/process-damager (manager) -- KCOMM /= MPI_COMM_WORLD may hang as sub-communicator
+
+#endif
+
+!CALL EC_RAISE(SIGABRT) ! In case ever ends up here
+STOP 'SDL_DISABORT'
+
+END SUBROUTINE SDL_DISABORT
+!-----------------------------------------------------------------------------
+
+END MODULE SDL_MOD
diff --git a/src/MNH/forcing.f90 b/src/MNH/forcing.f90
index 14a65ad60ebc6668b4a38851dc32ab1b06259041..522ac749e2f8ac6960198d55ea56ff653a6e76e3 100644
--- a/src/MNH/forcing.f90
+++ b/src/MNH/forcing.f90
@@ -163,6 +163,7 @@ use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, lbudg
USE MODD_CONF
USE MODD_CST
USE MODD_DYN
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_FRC
USE MODD_LUNIT
USE MODD_PARAMETERS
@@ -737,11 +738,13 @@ END IF
!
!* 4.2.1 integration of the tendency forcing for uv
!
-IF ( LTEND_UV_FRC ) THEN
- PRUS(:,:,:) = PRUS(:,:,:) + MXM(PRHODJ) * ZTENDUF(:,:,:)
- PRVS(:,:,:) = PRVS(:,:,:) + MYM(PRHODJ) * ZTENDVF(:,:,:)
+IF (.NOT. LOCEAN) THEN
+ IF ( LTEND_UV_FRC ) THEN
+ PRUS(:,:,:) = PRUS(:,:,:) + MXM(PRHODJ) * ZTENDUF(:,:,:)
+ PRVS(:,:,:) = PRVS(:,:,:) + MYM(PRHODJ) * ZTENDVF(:,:,:)
+ END IF
END IF
-!
+ !
!* 4.3 integration of the thermal and geostrophic wind
!
IF( LCORIO ) THEN
diff --git a/src/MNH/ibm_detect.f90 b/src/MNH/ibm_detect.f90
index b88a80445aad736413f45e64431eae8f3a3c3db1..3919727cc4297cc1ad3f7ebabbbe1b349333f47d 100644
--- a/src/MNH/ibm_detect.f90
+++ b/src/MNH/ibm_detect.f90
@@ -145,6 +145,7 @@ SUBROUTINE IBM_DETECT(PPHI)
ALLOCATE(I_INDE_GHOST(IIB:IIE,IJB:IJE,IKB:IKE,4))
ALLOCATE(ZIBM_TESTING(IIU,IJU,IKU,4))
ALLOCATE(Z_PHI(8),ZPROD(6),ZVECT(3),Z_IMG(8,3),Z_GHO(8,3),I_INDE_TEMPO(3),I_INDE_TEMPO2(3))
+ ALLOCATE(Z_NORM_GHOST(IIU,IJU,IKU,3),Z_NORM_TEMPO(IIU,IJU,IKU,3))
!
!------------------------------------------------------------------------------
!
@@ -464,7 +465,7 @@ SUBROUTINE IBM_DETECT(PPHI)
ALLOCATE(XIBM_IMAGE_P(I_DIME_GHOSTP,I_NUMB_LAYERP,1,3,3))
XIBM_IMAGE_P(:,:,:,:,:) = 0.
!
- ALLOCATE(Z_NORM_GHOST(IIU,IJU,IKU,3),Z_NORM_TEMPO(IIU,IJU,IKU,3),Z_NORM_TEMP1(IIU,IJU,IKU,4),Z_NORM_TEMP2(IIU,IJU,IKU), &
+ ALLOCATE(Z_NORM_TEMP1(IIU,IJU,IKU,4),Z_NORM_TEMP2(IIU,IJU,IKU), &
Z_NORM_TEMP3(IIU,IJU,IKU))
ALLOCATE(ZPHI(IIU,IJU,IKU,4))
ZPHI = 0.
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 2d81ad2b5d816dc3b1c07af42069ece9976fad73..1f51553ec73a01ae39a89c179a4d6df80acad088 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -13,8 +13,8 @@ INTERFACE
!
USE MODD_IO, ONLY: TFILEDATA
!
-INTEGER, INTENT(IN) :: KMI ! Model Index
-TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE ! Initial file
+INTEGER, INTENT(IN) :: KMI ! Model Index
+TYPE(TFILEDATA), POINTER, INTENT(IN) :: TPINIFILE ! Initial file
!
END SUBROUTINE INI_MODEL_n
!
@@ -494,8 +494,8 @@ IMPLICIT NONE
!* 0.1 declarations of arguments
!
!
-INTEGER, INTENT(IN) :: KMI ! Model Index
-TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE ! Initial file
+INTEGER, INTENT(IN) :: KMI ! Model Index
+TYPE(TFILEDATA), POINTER, INTENT(IN) :: TPINIFILE ! Initial file
!
!* 0.2 declarations of local variables
!
@@ -604,18 +604,19 @@ ILUOUT = TLUOUT%NLU
! --------------
!* 2.1 Read number of forcing fields
!
-IF (LFORCING) THEN ! Retrieve the number of time-dependent forcings.
+IF (.NOT. LOCEAN) THEN
+ IF (LFORCING) THEN ! Retrieve the number of time-dependent forcings.
CALL IO_Field_read(TPINIFILE,'FRC',NFRC,IRESP)
IF ( (IRESP /= 0) .OR. (NFRC <=0) ) THEN
WRITE(ILUOUT,'(A/A)') &
"INI_MODEL_n ERROR: you want to read forcing variables from FMfile", &
" but no fields have been found by IO_Field_read"
-!callabortstop
+ !callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','INI_MODEL_n','')
END IF
-END IF
-!
-! Modif PP for time evolving adv forcing
+ END IF
+ !
+ ! Modif PP for time evolving adv forcing
IF ( L2D_ADV_FRC ) THEN ! Retrieve the number of time-dependent forcings.
WRITE(ILUOUT,FMT=*) "INI_MODEL_n ENTER ADV_FORCING"
CALL IO_Field_read(TPINIFILE,'NADVFRC1',NADVFRC,IRESP)
@@ -627,9 +628,9 @@ END IF
CALL PRINT_MSG(NVERB_FATAL,'GEN','INI_MODEL_n','')
END IF
WRITE(ILUOUT,*) 'NADVFRC = ', NADVFRC
-END IF
-!
-IF ( L2D_REL_FRC ) THEN ! Retrieve the number of time-dependent forcings.
+ END IF
+ !
+ IF ( L2D_REL_FRC ) THEN ! Retrieve the number of time-dependent forcings.
WRITE(ILUOUT,FMT=*) "INI_MODEL_n ENTER REL_FORCING"
CALL IO_Field_read(TPINIFILE,'NRELFRC1',NRELFRC,IRESP)
IF ( (IRESP /= 0) .OR. (NRELFRC <=0) ) THEN
@@ -640,6 +641,7 @@ IF ( L2D_REL_FRC ) THEN ! Retrieve the number of time-dependent forcings.
CALL PRINT_MSG(NVERB_FATAL,'GEN','INI_MODEL_n','')
END IF
WRITE(ILUOUT,*) 'NRELFRC = ', NRELFRC
+ END IF
END IF
!* 2.2 Checks the position of vertical absorbing layer
!
@@ -2647,7 +2649,7 @@ IF (CSURF=='EXTE' .AND. (CPROGRAM=='MESONH' .OR. CPROGRAM=='DIAG ')) THEN
ENDIF
ELSE
! case after a spawning
- CINIFILEPGD = TPINIFILE%CNAME
+ TINIFILEPGD => TPINIFILE
END IF
!
CALL GOTO_SURFEX(KMI)
@@ -2978,16 +2980,20 @@ END IF
!
!* 33. Auto-coupling Atmos-Ocean LES NH
!
-IF (LCOUPLES) THEN
- ALLOCATE(XSSUFL_C(IIU,IJU,1)); XSSUFL_C=0.0
- ALLOCATE(XSSVFL_C(IIU,IJU,1)); XSSVFL_C=0.0
- ALLOCATE(XSSTFL_C(IIU,IJU,1)); XSSTFL_C=0.0
- ALLOCATE(XSSRFL_C(IIU,IJU,1)); XSSRFL_C=0.
-ELSE
- ALLOCATE(XSSUFL_C(0,0,0))
- ALLOCATE(XSSVFL_C(0,0,0))
- ALLOCATE(XSSTFL_C(0,0,0))
- ALLOCATE(XSSRFL_C(0,0,0))
+! Atmos Flux at interface
+! Allocate to a non-zero size only if LCOUPLES=T. It must be allocated only once (=>IF KMI==1)
+IF ( KMI == 1 ) THEN
+ IF ( LCOUPLES ) THEN
+ ALLOCATE( XSSUFL(IIU,IJU) ); XSSUFL=0.0
+ ALLOCATE( XSSVFL(IIU,IJU) ); XSSVFL=0.0
+ ALLOCATE( XSSTFL(IIU,IJU) ); XSSTFL=0.0
+ ALLOCATE( XSSRFL(IIU,IJU) ); XSSRFL=0.0
+ ELSE
+ ALLOCATE( XSSVFL(0,0) )
+ ALLOCATE( XSSTFL(0,0) )
+ ALLOCATE( XSSRFL(0,0) )
+ ALLOCATE( XSSUFL(0,0) )
+ END IF
END IF
!
END SUBROUTINE INI_MODEL_n
diff --git a/src/MNH/lesn.f90 b/src/MNH/lesn.f90
index 86b86a3e6e3a0b572e2686f371365d227c0713fb..7f1de616e99774a64fa58d65f0eca19adb8811d6 100644
--- a/src/MNH/lesn.f90
+++ b/src/MNH/lesn.f90
@@ -3277,7 +3277,7 @@ ELSE IF (CBL_HEIGHT_DEF=='FRI') THEN
+( XLES_SUBGRID_WV (:,NLES_CURRENT_TCOUNT,1) &
+XLES_RESOLVED_WV(:,NLES_CURRENT_TCOUNT,1))**2 )
ZFRIC_SURF = XLES_USTAR(NLES_CURRENT_TCOUNT)**2
- CALL BL_DEPTH_DIAG(YLDIMPHYEX,ZFRIC_SURF, XLES_ZS, &
+ CALL BL_DEPTH_DIAG(YLDIMPHYEX,SIZE(ZFRIC_LES),ZFRIC_SURF, XLES_ZS, &
ZFRIC_LES, XLES_Z, &
XFTOP_O_FSURF,XLES_BL_HEIGHT(NLES_CURRENT_TCOUNT))
END IF
diff --git a/src/MNH/mean_field.f90 b/src/MNH/mean_field.f90
index 573046c23188f3d2f7edd9804864543ceebb933b..a8a5cd98c508dc08f14cfb3f405ae686c3d5a230 100644
--- a/src/MNH/mean_field.f90
+++ b/src/MNH/mean_field.f90
@@ -67,6 +67,8 @@ USE MODD_CST
USE MODD_PASPOL
USE MODE_THERMO
USE MODI_SHUMAN
+! To play with MPI
+USE MODD_ARGSLIST_ll, ONLY: LIST_ll
!
! * EOL
USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
@@ -77,6 +79,8 @@ USE MODD_EOL_ADR, ONLY: XFAERO_BLEQ_RA_GLB, XFAERO_BLEQ_RA_SUM
USE MODD_EOL_ADR, ONLY: XAOA_BLEQ_GLB, XAOA_BLEQ_SUM
USE MODD_EOL_ALM, ONLY: XFAERO_RE_SUM, XFAERO_RE_GLB
!
+USE MODI_SHUMAN, ONLY: MXF, MYF, MZF
+!
USE MODE_MODELN_HANDLER
USE MODI_UPDATE_WELFORD
!
@@ -97,7 +101,6 @@ INTEGER :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
INTEGER :: JI,JJ,JK ! Loop indexes
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZUMEAN_OLD
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWMEAN_OLD
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHMEAN_OLD
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQSAT_W
@@ -115,6 +118,9 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZRH_P_MAXCOL
!
INTEGER :: IMI !Current model index
!
+! For communication
+TYPE(LIST_ll), POINTER :: TZFIELDS_Cov_ll ! Field list of Wind for exchange beofre covariances
+INTEGER :: IINFO ! Info integer
!
!-----------------------------------------------------------------------
!
@@ -228,10 +234,8 @@ END IF
IF (LCOV_FIELD) THEN
ALLOCATE( ZUMEAN_OLD (IIU, IJU, IKU) )
ALLOCATE( ZWMEAN_OLD (IIU, IJU, IKU) )
- ALLOCATE( ZTHMEAN_OLD(IIU, IJU, IKU) )
ZUMEAN_OLD (:,:,:) = XUM_MEAN (:,:,:)
ZWMEAN_OLD (:,:,:) = XWM_MEAN (:,:,:)
- ZTHMEAN_OLD(:,:,:) = XTHM_MEAN(:,:,:)
END IF
! Welford method for variables whom we compute variances
!
@@ -244,12 +248,24 @@ END IF
!
! Welford method for covariances
!
- IF (LCOV_FIELD) THEN
- XUV_MEAN = XUV_MEAN + (PUT-ZUMEAN_OLD)*(PVT-XVM_MEAN)
- XUW_MEAN = XUW_MEAN + (PUT-ZUMEAN_OLD)*(PWT-XWM_MEAN)
- XVW_MEAN = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PVT-XVM_MEAN)
- XWTH_MEAN = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PTHT-XTHM_MEAN)
- DEALLOCATE( ZUMEAN_OLD, ZWMEAN_OLD, ZTHMEAN_OLD )
+ IF (LCOV_FIELD) THEN
+!
+ NULLIFY(TZFIELDS_Cov_ll)
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PUT, 'MEAN_FIELD::PUT')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PWT, 'MEAN_FIELD::PWT')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PVT, 'MEAN_FIELD::PVT')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,ZUMEAN_OLD, 'MEAN_FIELD::ZUMEAN_OLD')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,ZWMEAN_OLD, 'MEAN_FIELD::ZWMEAN_OLD')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,XVM_MEAN, 'MEAN_FIELD::XVM_MEAN')
+ CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,XWM_MEAN, 'MEAN_FIELD::XWM_MEAN')
+ CALL UPDATE_HALO_ll(TZFIELDS_Cov_ll,IINFO)
+ CALL CLEANLIST_ll( TZFIELDS_Cov_ll)
+!
+ XUV_MEAN = XUV_MEAN + (MXF(PUT)-MXF(ZUMEAN_OLD))*(MYF(PVT)-MYF(XVM_MEAN))
+ XUW_MEAN = XUW_MEAN + (MXF(PUT)-MXF(ZUMEAN_OLD))*(MZF(PWT)-MZF(XWM_MEAN))
+ XVW_MEAN = XVW_MEAN + (MZF(PWT)-MZF(ZWMEAN_OLD))*(MYF(PVT)-MYF(XVM_MEAN))
+ XWTH_MEAN = XWTH_MEAN + (MZF(PWT)-MZF(ZWMEAN_OLD))*(PTHT-XTHM_MEAN)
+ DEALLOCATE( ZUMEAN_OLD, ZWMEAN_OLD )
END IF
!
!-----------------------------------------------------------------------
diff --git a/src/MNH/modd_oceanh.f90 b/src/MNH/modd_oceanh.f90
index e9936173ec35b1d57c9884cd19cd98774c9d0334..0716105b6ca32e75aad10eefc2768504d0e9dea2 100644
--- a/src/MNH/modd_oceanh.f90
+++ b/src/MNH/modd_oceanh.f90
@@ -40,8 +40,7 @@ IMPLICIT NONE
INTEGER, SAVE :: NFRCLT ! number of sea surface forcings PLUS 1
INTEGER, SAVE :: NINFRT ! Interval in second between forcings
TYPE (DATE_TIME), SAVE, DIMENSION(:), ALLOCATABLE :: TFRCLT ! date/time of sea surface forcings
-REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XSSUFL,XSSVFL,XSSTFL,XSSOLA ! Time evol Flux U V T Solar_Rad at sea surface
-REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XSSUFL_XY,XSSVFL_XY,XSSTFL_XY! XY flux shape
+REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XSSUFL,XSSVFL,XSSTFL,XSSRFL,XSSOLA ! Time evol Flux U,V,T,SW at sea surface
REAL, SAVE, DIMENSION(:), ALLOCATABLE :: XSSUFL_T,XSSVFL_T,XSSTFL_T,XSSOLA_T ! given time forcing fluxes
!
END MODULE MODD_OCEANH
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 209b07bd4db2a30db80b55888af34ac5f6b2f035..e988f0afef109a750eee64c3aefa107d2fad0e12 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -2355,7 +2355,7 @@ IF (OEXIT) THEN
END IF
!
CALL IO_File_close(TINIFILE)
- IF (CSURF=="EXTE") CALL IO_File_close(TINIFILEPGD)
+ IF ( CSURF == "EXTE" .AND. ASSOCIATED( TINIFILEPGD ) ) CALL IO_File_close( TINIFILEPGD )
!
!* 28.1 print statistics!
!
diff --git a/src/MNH/nhoa_coupln.f90 b/src/MNH/nhoa_coupln.f90
index 072dfbeb3ca7aa64c133a8d0b2855e4a3d1fdea8..b156d2a4ccbf8f63a16220bf99a8ec96ec48cbbc 100644
--- a/src/MNH/nhoa_coupln.f90
+++ b/src/MNH/nhoa_coupln.f90
@@ -60,13 +60,15 @@ USE MODE_MODELN_HANDLER
USE MODD_PARAMETERS
USE MODD_NESTING
USE MODD_CST
+USE MODD_REF, ONLY: XRHODREFZ
USE MODD_REF_n ! modules relative to the outer model $n
USE MODD_FIELD_n
+USE MODD_METRICS_n
USE MODD_CONF
USE MODD_PARAM_n
USE MODD_TURB_n
USE MODD_DYN_n, ONLY : LOCEAN
-USE MODD_REF, ONLY: LCOUPLES
+USE MODD_OCEANH
!
IMPLICIT NONE
!
@@ -80,7 +82,7 @@ INTEGER, INTENT(IN) :: KKU !
!
!* 0.2 declarations of local variables
!
-INTEGER :: IIB,IIE,IJB,IJE,IIU,IJU
+INTEGER :: IIU,IJU
INTEGER :: IKE
INTEGER :: ILBX,ILBY,ILBX2,ILBY2
INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop COUNTer
@@ -93,7 +95,8 @@ INTEGER :: IINFO_ll, IDIMX, IDIMY
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPUA,ZCOUPVA,ZCOUPTA
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPUO,ZCOUPVO,ZCOUPTO
!surf flux local work space
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPTFL,ZCOUPUFL,ZCOUPVFL
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPTFL,ZCOUPUFL,ZCOUPVFL,ZCOUPQFL
+REAL :: ZCDN ! CONSTANT DRAG COEFF
CHARACTER(LEN=4) :: ZINIT_TYPE
!
!---Coupled OA MesoNH----------------------------------------------------------------------------
@@ -101,11 +104,7 @@ CHARACTER(LEN=4) :: ZINIT_TYPE
! --------------
IIU = SIZE(XRHODJ,1)
IJU = SIZE(XRHODJ,2)
-
-! allocate flux local array
-ALLOCATE(ZCOUPTFL(IIU,IJU))
-ALLOCATE(ZCOUPUFL(IIU,IJU))
-ALLOCATE(ZCOUPVFL(IIU,IJU))
+ZCDN = 1.2E-3
! allocate sfc variable local array
ALLOCATE(ZCOUPUA(IIU,IJU))
ALLOCATE(ZCOUPVA(IIU,IJU))
@@ -113,7 +112,11 @@ ALLOCATE(ZCOUPTA(IIU,IJU))
ALLOCATE(ZCOUPUO(IIU,IJU))
ALLOCATE(ZCOUPVO(IIU,IJU))
ALLOCATE(ZCOUPTO(IIU,IJU))
-! values in ocean sfc
+ALLOCATE(ZCOUPTFL(IIU,IJU))
+ALLOCATE(ZCOUPUFL(IIU,IJU))
+ALLOCATE(ZCOUPVFL(IIU,IJU))
+ALLOCATE(ZCOUPQFL(IIU,IJU))
+! loading values in ocean sfc
IKE=KKU-JPVEXT
ZCOUPUO(:,:)= XUT(:,:,IKE)
ZCOUPVO(:,:)= XVT(:,:,IKE)
@@ -121,37 +124,31 @@ ZCOUPTO(:,:)= XTHT(:,:,IKE)
!
! we are going to the atmos model i.e. Model 1
CALL GOTO_MODEL(KDAD)
-IIB=1
-IIE=IIU
-IJB=1
-IJE=IJU
-!
! compute gradient between ocean & atmosphere
ZCOUPUA(:,:)= XUT(:,:,2)-ZCOUPUO(:,:)
ZCOUPVA(:,:)= XVT(:,:,2)-ZCOUPVO(:,:)
ZCOUPTA(:,:)= XTHT(:,:,2)-ZCOUPTO(:,:)
-!
-! sfc flux computation * RHO AIR !!!!
-! flux vu atmosp
-!
-ZCOUPTFL(:,:) = -1.2*1.E-3* SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPTA(:,:)
-ZCOUPUFL(:,:) = -1.2*1.E-3* SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPUA(:,:)
-ZCOUPVFL(:,:) = -1.2*1.E-3* SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPVA(:,:)
-!
-XSSUFL_C(:,:,1)= ZCOUPUFL(:,:)
-XSSVFL_C(:,:,1)= ZCOUPVFL(:,:)
-XSSTFL_C(:,:,1)= ZCOUPTFL(:,:)
-!
-!
-! We are going back in the ocean model
-! same sign & unit at the top of ocean model and the bottom of atmospheric model
+!!!!!!!!!!!!
+! Simple neutral bulk flux with Z0 fixed
+ZCOUPTFL(:,:) = -ZCDN*SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPTA(:,:)&
+ *XCPD*XRHODREFZ(2) ! in W/M2 (as WASP)
+ZCOUPUFL(:,:) = -ZCDN*XRHODREFZ(2)* SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPUA(:,:)
+ZCOUPVFL(:,:) = -ZCDN*XRHODREFZ(2)* SQRT(ZCOUPUA(:,:)**2 +ZCOUPVA(:,:)**2) * ZCOUPVA(:,:)
+! EVAP-PR TO BE ADDED HERE
+ZCOUPQFL(:,:)= 0.
+! keep in W/M2, N/m2, ..
+XSSUFL= ZCOUPUFL
+XSSVFL= ZCOUPVFL
+XSSTFL= ZCOUPTFL
+XSSRFL= ZCOUPQFL
+!
+! back in the ocean model
+!oce flux are directly computed in phys_param routine
+! same sign & unit fluxes at top of ocean model & bottom of atmospheric model
! rho_atmos * (w'u')_atmos = rho_ocean * (u'w')_ocean
-! rho_atmos *Cp_atmos* (u'w')_atmos = rho_ocean * CP_ocean * (u'w')_ocean
+! rho_atmos *Cp_atmos* (T'w')_atmos = rho_ocean * CP_ocean * (T'w')_ocean
!
CALL GOTO_MODEL(KMI)
-XSSUFL_C(:,:,1)= ZCOUPUFL(:,:)/XRH00OCEAN
-XSSVFL_C(:,:,1)= ZCOUPVFL(:,:)/XRH00OCEAN
-XSSTFL_C(:,:,1)= ZCOUPTFL(:,:)*1004./(3900.*XRH00OCEAN)
DEALLOCATE(ZCOUPUA,ZCOUPVA,ZCOUPUO,ZCOUPVO,ZCOUPTA,ZCOUPTO)
DEALLOCATE(ZCOUPTFL)
!
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index edd731e427f3ff0b3b502f8ef82bfb74052ade40..ec037233cbe8673a3ba5815ece7f3839400a3e0b 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -317,7 +317,7 @@ USE MODD_PRECIP_n
use modd_precision, only: MNHTIME
USE MODD_RADIATIONS_n
USE MODD_RAIN_ICE_DESCR_n, ONLY: XRTMIN
-USE MODD_REF, ONLY: LCOUPLES
+USE MODD_REF, ONLY: LCOUPLES,XRHODREFZ
USE MODD_REF_n
USE MODD_SALT
USE MODD_SHADOWS_n
@@ -865,24 +865,32 @@ END IF
! Sfc turbulent fluxes & Radiative tendency due to SW penetrating ocean
!
IF (LCOUPLES) THEN
-ZSFU(:,:)= XSSUFL_C(:,:,1)
-ZSFV(:,:)= XSSVFL_C(:,:,1)
-ZSFTH(:,:)= XSSTFL_C(:,:,1)
-ZSFRV(:,:)=XSSRFL_C(:,:,1)
-ELSE
-IF (LOCEAN) THEN
-!
+! Flux computed in NHOA/WASP and are given in W/M2
+! convert in SI as in MESO-NH & for Ocean & Atmos model
+ IF (LOCEAN) THEN
+ ZSFU= XSSUFL/CST%XRH00OCEAN
+ ZSFV= XSSVFL/CST%XRH00OCEAN
+ ZSFTH= XSSTFL/(CST%XCL*CST%XRH00OCEAN)
+!to do: to add factor s/(1-s)
+ ZSFRV=-XSSRFL/CST%XRH00OCEAN
+ ELSE
+ ZSFU= XSSUFL/XRHODREFZ(2)
+ ZSFV= XSSVFL/XRHODREFZ(2)
+ ZSFTH=XSSTFL/(XCPD*XRHODREFZ(2))
+ ZSFRV=XSSRFL/(CST%XLVTT*XRHODREFZ(2))
+ ENDIF
+ELSE ! no coupled
+ IF (LOCEAN) THEN
ALLOCATE( ZIZOCE(IKU)); ZIZOCE(:)=0.
ALLOCATE( ZPROSOL1(IKU))
ALLOCATE( ZPROSOL2(IKU))
- ALLOCATE(XSSOLA(IIU,IJU))
! Time interpolation
JSW = INT(TDTCUR%xtime/REAL(NINFRT))
ZSWA = TDTCUR%xtime/REAL(NINFRT)-REAL(JSW)
ZSFRV = 0.
ZSFTH = (XSSTFL_T(JSW+1)*(1.-ZSWA)+XSSTFL_T(JSW+2)*ZSWA)
- ZSFU = (XSSUFL_T(JSW+1)*(1.-ZSWA)+XSSUFL_T(JSW+2)*ZSWA)
- ZSFV = (XSSVFL_T(JSW+1)*(1.-ZSWA)+XSSVFL_T(JSW+2)*ZSWA)
+ ZSFU = (XSSUFL_T(JSW+1)*(1.-ZSWA)+XSSUFL_T(JSW+2)*ZSWA)
+ ZSFV = (XSSVFL_T(JSW+1)*(1.-ZSWA)+XSSVFL_T(JSW+2)*ZSWA)
!
ZIZOCE(IKU) = XSSOLA_T(JSW+1)*(1.-ZSWA)+XSSOLA_T(JSW+2)*ZSWA
ZPROSOL1(IKU) = CST%XROC*ZIZOCE(IKU)
@@ -896,12 +904,11 @@ IF (LOCEAN) THEN
XRTHS(:,:,JKM) = XRTHS(:,:,JKM) + XRHODJ(:,:,JKM)*ZIZOCE(JKM)
END DO
if ( TBUCONF%LBUDGET_th ) call Budget_store_end ( TBUDGETS(NBUDGET_TH), 'OCEAN', xrths(:, :, :) )
- DEALLOCATE (XSSOLA)
DEALLOCATE( ZIZOCE)
DEALLOCATE (ZPROSOL1)
DEALLOCATE (ZPROSOL2)
-END IF! LOCEAN NO LCOUPLES
-END IF!NO LCOUPLES
+END IF! TEST LOCEAN NO LCOUPLES
+END IF! TEST LCOUPLES
!
!
CALL SECOND_MNH2(ZTIME2)
@@ -1324,7 +1331,7 @@ ELSE ! case no SURFEX (CSURF logical)
ZCD_ROOF = 0.
ZSFRV_WALL = 0.
ZSFRV_ROOF = 0.
- IF (.NOT.LOCEAN) THEN
+ IF (.NOT.(LOCEAN.OR.LCOUPLES))THEN
ZSFTH = 0.
ZSFRV = 0.
ZSFSV = 0.
diff --git a/src/MNH/read_field.f90 b/src/MNH/read_field.f90
index 00593e5415708626fd57a918fba957f323943fd7..bd4a70fcf039c3bb9d9016621b09ff65afecb5c0 100644
--- a/src/MNH/read_field.f90
+++ b/src/MNH/read_field.f90
@@ -1324,7 +1324,7 @@ IF (LOCEAN .AND. (.NOT.LCOUPLES) .AND. (KOCEMI==1)) THEN
CLONGNAME = 'SSOLA', &
CUNITS = 'kg m3 K m s-1', &
CDIR = '--', &
- CCOMMENT = 'sfc solar flux to force ocean LES', &
+ CCOMMENT = 'sfc solar flux at sfc to force ocean LES', &
NGRID = 0, &
NTYPE = TYPEREAL, &
NDIMS = 1, &
diff --git a/src/MNH/set_rsou.f90 b/src/MNH/set_rsou.f90
index 4526251833db633bf65d0808d8f86cf6ad973943..e9a4dc78fe7ce1bdc111db69603d1dce3b4f90e7 100644
--- a/src/MNH/set_rsou.f90
+++ b/src/MNH/set_rsou.f90
@@ -272,6 +272,7 @@ USE MODD_NETCDF
USE MODD_OCEANH
USE MODD_PARAMETERS, ONLY: JPHEXT
USE MODD_TYPE_DATE
+USE MODD_TIME_n, ONLY : TDTCUR
!
USE MODE_ll
USE MODE_MSG
@@ -400,7 +401,8 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZOC_DEPTH
REAL, DIMENSION(:), ALLOCATABLE :: ZOC_LE,ZOC_H
REAL, DIMENSION(:), ALLOCATABLE :: ZOC_SW_DOWN,ZOC_SW_UP,ZOC_LW_DOWN,ZOC_LW_UP
REAL, DIMENSION(:), ALLOCATABLE :: ZOC_TAUX,ZOC_TAUY
-
+!
+REAL :: ZJZTIME ! TIME(HOUR) READ in PRE_IDEA1.NAM
!--------------------------------------------------------------------------------
!
!* 1. PROLOGUE : INITIALIZE SOME CONSTANTS, RETRIEVE LOGICAL
@@ -565,9 +567,10 @@ SELECT CASE(YKIND)
TFRCLT(JKT)= ZFRCLT(JKT)
XSSUFL_T(JKT)=ZSSUFL_T(JKT)/XRH00OCEAN
XSSVFL_T(JKT)=ZSSVFL_T(JKT)/XRH00OCEAN
- ! working in SI
- XSSTFL_T(JKT)=ZSSTFL_T(JKT) /(3900.*XRH00OCEAN)
- XSSOLA_T(JKT)=ZSSOLA_T(JKT) /(3900.*XRH00OCEAN)
+
+ ! working in SI
+ XSSTFL_T(JKT)=ZSSTFL_T(JKT) /(XCL*XRH00OCEAN)
+ XSSOLA_T(JKT)=ZSSOLA_T(JKT) /(XCL*XRH00OCEAN)
END DO
DEALLOCATE(ZFRCLT)
DEALLOCATE(ZSSUFL_T)
@@ -579,11 +582,12 @@ SELECT CASE(YKIND)
! 2.0.2 Ocean standard initialize from netcdf files
! U,V,T,S at Z levels + Forcings at model TOP (sea surface)
!--------------------------------------------------------------------------------
-!
+! WRITTEN TO READ NETCDF file for a DWL case of DYNAMO-CINDY
CASE ('STANDOCE')
-!
XP00=XP00OCEAN
- READ(ILUPRE,*) ZPTOP ! P_atmosphere at sfc =P top domain
+! HOUR OF FIRST FORCING on 14/11 TO USE AS READ in PRE_IDEA1.NAM
+ ZJZTIME= INT(TDTCUR%xtime)
+ READ(ILUPRE,*) ZPTOP ! P_atmosphere at sfc =P top ocean domain
READ(ILUPRE,*) YINFILE, YINFISF
WRITE(ILUOUT,FMT=*) 'Netcdf files to read:', YINFILE, YINFISF
! Open file containing initial profiles
@@ -670,6 +674,8 @@ SELECT CASE(YKIND)
! ZHEIGHT used only in set_ rsou, defined as such ZHEIGHT(ILEVELM)=H_model
! Z oriented in same time to have a model domain axis going
! from 0m (ocean bottom/model bottom) towards H (ocean sfc/model top)
+ WRITE(ILUOUT,FMT=*) 'End gridmodel comput in trans domain: JKM U V ZHEIGHTU', &
+ JKM,ZU(JKM),ZV(JKM),ZHEIGHTU(JKM)
END DO
!
DEALLOCATE(ZOC_TEMPERATURE)
@@ -728,30 +734,36 @@ SELECT CASE(YKIND)
WRITE(ILUOUT,FMT=*) JKM, ZOC_LE(JKM), ZOC_H(JKM),ZOC_SW_DOWN(JKM),ZOC_SW_UP(JKM),&
ZOC_LW_DOWN(JKM),ZOC_LW_UP(JKM),ZOC_TAUX(JKM),ZOC_TAUY(JKM)
ENDDO
- ! IFRCLT FORCINGS at sea surface
- IFRCLT=idimlen
+ ! IFRCLT FORCINGS at sea surface are givean each 10 min (6/hour)
+! start at ZJZTIME read in PRE_IDEA1.nam; skip forcings given from 13/11 22LT
+ZJZTIME=(ZJZTIME+2)*6
+ IFRCLT=idimlen-ZJZTIME
+WRITE(ILUOUT,FMT=*) 'FORCINGS GIVEN FROM 13/11/11 22LT; FORCING Number USED ZJZTIME=', ZJZTIME
+WRITE(ILUOUT,FMT=*) idimlen,'NB of given forcing', 'nb of first used forcing= ',IFRCLT
ALLOCATE(ZFRCLT(IFRCLT))
ALLOCATE(ZSSUFL_T(IFRCLT)); ZSSUFL_T = 0.0
ALLOCATE(ZSSVFL_T(IFRCLT)); ZSSVFL_T = 0.0
ALLOCATE(ZSSTFL_T(IFRCLT)); ZSSTFL_T = 0.0
ALLOCATE(ZSSOLA_T(IFRCLT)); ZSSOLA_T = 0.0
- DO JKT=1,IFRCLT
+ DO JKM=1,IFRCLT
+ JKT=JKM +ZJZTIME
! Initial file for CINDY-DYNAMO: all fluxes correspond to the absolute value (>0)
! modele ocean: axe z dirigé du bas vers la sfc de l'océan
! => flux dirigé vers le haut (positif ocean vers l'atmopshere i.e. bas vers le haut)
- ZSSOLA_T(JKT)=ZOC_SW_DOWN(JKT)-ZOC_SW_UP(JKT)
- ZSSTFL_T(JKT)=(ZOC_LW_DOWN(JKT)-ZOC_LW_UP(JKT)-ZOC_LE(JKT)-ZOC_H(JKT))
+ ZSSOLA_T(JKM)=ZOC_SW_DOWN(JKT)-ZOC_SW_UP(JKT)
+ ZSSTFL_T(JKM)=-(ZOC_LW_DOWN(JKT)-ZOC_LW_UP(JKT)-ZOC_LE(JKT)-ZOC_H(JKT))
! assume that Tau given on file is along Ox
! rho_air UW_air = rho_ocean UW_ocean= N/m2
! uw_ocean
- ZSSUFL_T(JKT)=ZOC_TAUX(JKT)
- ZSSVFL_T(JKT)=ZOC_TAUY(JKT)
- WRITE(ILUOUT,FMT=*) 'Forcing Nb Sol NSol UW_oc VW',&
- JKT,ZSSOLA_T(JKT),ZSSTFL_T(JKT),ZSSUFL_T(JKT),ZSSVFL_T(JKT)
+ ZSSUFL_T(JKM)=ZOC_TAUX(JKT)
+ ZSSVFL_T(JKM)=ZOC_TAUY(JKT)
+ WRITE(ILUOUT,FMT=*) 'Forcing Nb Sol NSol UW_oc VW',&
+ JKM,ZSSOLA_T(JKM),ZSSTFL_T(JKM),ZSSUFL_T(JKM),ZSSVFL_T(JKM)
ENDDO
! Allocate and Writing the corresponding variables in module MODD_OCEAN_FRC
NFRCLT=IFRCLT
! value to read later on file ?
+! interval in s between 2 given forcings
NINFRT=600
ALLOCATE(TFRCLT(NFRCLT))
ALLOCATE(XSSUFL_T(NFRCLT));XSSUFL_T(:)=0.
@@ -759,14 +771,13 @@ SELECT CASE(YKIND)
ALLOCATE(XSSTFL_T(NFRCLT));XSSTFL_T(:)=0.
ALLOCATE(XSSOLA_T(NFRCLT));XSSOLA_T(:)=0.
! on passe en unités SI, signe, etc pour le modele ocean
- ! W/m2 => SI : /(CP_mer * rho_mer)
- ! a revoir dans tt le code pour mettre de svaleurs plus exactes
+ ! W/m2 => SI : m/S *K /(CP_mer * rho_mer)
DO JKT=1,NFRCLT
TFRCLT(JKT)= ZFRCLT(JKT)
XSSUFL_T(JKT)=ZSSUFL_T(JKT)/XRH00OCEAN
XSSVFL_T(JKT)=ZSSVFL_T(JKT)/XRH00OCEAN
- XSSTFL_T(JKT)=ZSSTFL_T(JKT) /(3900.*XRH00OCEAN)
- XSSOLA_T(JKT)=ZSSOLA_T(JKT) /(3900.*XRH00OCEAN)
+ XSSTFL_T(JKT)=ZSSTFL_T(JKT) /(XCL *XRH00OCEAN)
+ XSSOLA_T(JKT)=ZSSOLA_T(JKT) /(XCL *XRH00OCEAN)
END DO
DEALLOCATE(ZFRCLT)
DEALLOCATE(ZSSUFL_T)
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index 5d6517151f1749475da0f4ed554eafcf22efb772..f7f185e05649dd30a4d50bf95a6b7d0734fc227d 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -759,38 +759,6 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CUNITS = 'm s-1'
TZFIELD%CCOMMENT = 'X_Y_Z_U component of max wind'
CALL IO_Field_write(TPFILE,TZFIELD,XUM_MAX)
-!
- IF (LCOV_FIELD) THEN
- !
- TZFIELD%CMNHNAME = 'UVME'
- TZFIELD%CLONGNAME = 'UVME'
- TZFIELD%CUNITS = 'm2 s-2'
- TZFIELD%CCOMMENT = 'X_Y_Z_UV component of mean wind variance'
- ZWORK3D = XUV_MEAN/MEAN_COUNT
- CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
- !
- TZFIELD%CMNHNAME = 'UWME'
- TZFIELD%CLONGNAME = 'UWME'
- TZFIELD%CUNITS = 'm2 s-2'
- TZFIELD%CCOMMENT = 'X_Y_Z_UW component of mean wind variance'
- ZWORK3D = XUW_MEAN/MEAN_COUNT
- CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
- !
- TZFIELD%CMNHNAME = 'VWME'
- TZFIELD%CLONGNAME = 'VWME'
- TZFIELD%CUNITS = 'm2 s-2'
- TZFIELD%CCOMMENT = 'X_Y_Z_VW component of mean wind variance'
- ZWORK3D = XVW_MEAN/MEAN_COUNT
- CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
- !
- TZFIELD%CMNHNAME = 'WTHME'
- TZFIELD%CLONGNAME = 'WTHME'
- TZFIELD%CUNITS = 'm2 s-2'
- TZFIELD%CCOMMENT = 'X_Y_Z_WTH component of mean wind variance'
- ZWORK3D = XWTH_MEAN/MEAN_COUNT
- CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
- !
- END IF
!
TZFIELD = TFIELDMETADATA( &
CMNHNAME = 'generic for mean_count variables', & !Temporary name to ease identification
@@ -1030,6 +998,38 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CCOMMENT = 'X_Y_Z_max kinetic energy'
CALL IO_Field_write(TPFILE,TZFIELD,XTKEM_MAX)
END IF
+!
+ IF (LCOV_FIELD) THEN
+ !
+ TZFIELD%CMNHNAME = 'UVME'
+ TZFIELD%CLONGNAME = 'UVME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UV component of mean wind variance'
+ ZWORK3D = XUV_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'UWME'
+ TZFIELD%CLONGNAME = 'UWME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UW component of mean wind variance'
+ ZWORK3D = XUW_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'VWME'
+ TZFIELD%CLONGNAME = 'VWME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_VW component of mean wind variance'
+ ZWORK3D = XVW_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'WTHME'
+ TZFIELD%CLONGNAME = 'WTHME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_WTH component of mean wind variance'
+ ZWORK3D = XWTH_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ END IF
!
END IF
!
diff --git a/src/Makefile.MESONH.mk b/src/Makefile.MESONH.mk
index c5f5aeadda9e915c1b2790d524cc47f1267ed8b1..caebcb7e1ad3de5cd15794506423becc93991b50 100644
--- a/src/Makefile.MESONH.mk
+++ b/src/Makefile.MESONH.mk
@@ -145,7 +145,7 @@ CPPFLAGS_RAD = -DMNH_ECRAD -DVER_ECRAD=$(VER_ECRAD)
INC_RAD += -I$(B)LIB/RAD/ecrad-$(VERSION_ECRAD)/include
ifeq "$(VER_ECRAD)" "140"
INC_RAD += -I$(B)LIB/RAD/ecrad-$(VERSION_ECRAD)/drhook/include
-#IGNORE_DEP_MASTER += yomhook.D
+IGNORE_DEP_MASTER += yomhook_dummy.D
endif
ifneq "$(VER_ECRAD)" "140"
IGNORE_DEP_MASTER += read_albedo_data.D read_emiss_data.D
diff --git a/src/PHYEX/turb/mode_bl_depth_diag.f90 b/src/PHYEX/turb/mode_bl_depth_diag.f90
index 2cae3a3fe2cf48cfd1f0a9a4d55c52b57abde3cf..daafb106e2d2ebc995bb4add8c7aef8d6a9738bd 100644
--- a/src/PHYEX/turb/mode_bl_depth_diag.f90
+++ b/src/PHYEX/turb/mode_bl_depth_diag.f90
@@ -12,7 +12,7 @@ END INTERFACE
!
CONTAINS
!
-SUBROUTINE BL_DEPTH_DIAG_3D(D,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG3D)
+SUBROUTINE BL_DEPTH_DIAG_3D(D,KDIM, PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG3D)
USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
!
!
@@ -56,10 +56,11 @@ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
IMPLICIT NONE
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KDIM ! Vertical dimensions (LES grid vs full grid)
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSURF ! surface flux
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PZS ! orography
-REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PFLUX ! flux
-REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZZ ! altitude of flux points
+REAL, DIMENSION(D%NIJT,KDIM), INTENT(IN) :: PFLUX ! flux
+REAL, DIMENSION(D%NIJT,KDIM), INTENT(IN) :: PZZ ! altitude of flux points
REAL, INTENT(IN) :: PFTOP_O_FSURF! Flux at BL top / Surface flux
REAL, DIMENSION(D%NIJT), INTENT(OUT) :: BL_DEPTH_DIAG3D
!
@@ -106,7 +107,7 @@ BL_DEPTH_DIAG3D(:) = BL_DEPTH_DIAG3D(:) / (1. - PFTOP_O_FSURF)
IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_3D',1,ZHOOK_HANDLE)
END SUBROUTINE BL_DEPTH_DIAG_3D
!
-SUBROUTINE BL_DEPTH_DIAG_1D(D,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG1D)
+SUBROUTINE BL_DEPTH_DIAG_1D(D,KDIM, PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG1D)
USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
!
USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
@@ -114,31 +115,37 @@ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
IMPLICIT NONE
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KDIM ! Vertical dimensions (LES grid vs full grid)
REAL, INTENT(IN) :: PSURF ! surface flux
REAL, INTENT(IN) :: PZS ! orography
-REAL, DIMENSION(D%NKT), INTENT(IN) :: PFLUX ! flux
-REAL, DIMENSION(D%NKT), INTENT(IN) :: PZZ ! altitude of flux points
+REAL, DIMENSION(KDIM), INTENT(IN) :: PFLUX ! flux
+REAL, DIMENSION(KDIM), INTENT(IN) :: PZZ ! altitude of flux points
REAL, INTENT(IN) :: PFTOP_O_FSURF! Flux at BL top / Surface flux
REAL, INTENT(OUT) :: BL_DEPTH_DIAG1D
+INTEGER :: JK ! loop counters
+INTEGER :: IKB,IIJB,IIJE,IKL
+REAL :: ZFLX ! flux at top of BL
!
-REAL, DIMENSION(1,1) :: ZSURF
-REAL, DIMENSION(1,1) :: ZZS
-REAL, DIMENSION(1,1,D%NKT) :: ZFLUX
-REAL, DIMENSION(1,1,D%NKT) :: ZZZ
-REAL, DIMENSION(1,1) :: ZBL_DEPTH_DIAG
-!
-INTEGER :: IKT
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',0,ZHOOK_HANDLE)
-IKT=D%NKT
-ZSURF = PSURF
-ZZS = PZS
-ZFLUX(1,1,:) = PFLUX(:)
-ZZZ (1,1,:) = PZZ (:)
!
-CALL BL_DEPTH_DIAG_3D(D,ZSURF,ZZS,ZFLUX,ZZZ,PFTOP_O_FSURF,ZBL_DEPTH_DIAG)
+IKB=D%NKTB
+IKL=D%NKL
!
-BL_DEPTH_DIAG1D = ZBL_DEPTH_DIAG(1,1)
+BL_DEPTH_DIAG1D = 0.
+IF (PSURF/=0.) THEN
+ DO JK=IKB,KDIM-1,IKL
+ IF (PZZ(JK-IKL)>PZS) THEN
+ ZFLX = PSURF * PFTOP_O_FSURF
+ IF ( (PFLUX(JK)-ZFLX)*(PFLUX(JK-IKL)-ZFLX) <= 0. ) THEN
+ BL_DEPTH_DIAG1D = (PZZ (JK-IKL) - PZS) &
+ + (PZZ (JK) - PZZ (JK-IKL)) &
+ * (ZFLX - PFLUX(JK-IKL) ) &
+ / (PFLUX(JK) - PFLUX(JK-IKL) )
+ END IF
+ END IF
+ END DO
+END IF
!
!-------------------------------------------------------------------------------
!
diff --git a/src/PHYEX/turb/mode_sbl_depth.f90 b/src/PHYEX/turb/mode_sbl_depth.f90
index 8da65ac5643dc573d4e5706c1a4b76cbb8456464..948948c07295d73c570f98965799ee77928bea01 100644
--- a/src/PHYEX/turb/mode_sbl_depth.f90
+++ b/src/PHYEX/turb/mode_sbl_depth.f90
@@ -109,7 +109,7 @@ ZUSTAR2(:) = SQRT(ZWU(:)**2+ZWV(:)**2)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZWIND(:,:)=SQRT(PFLXU(:,:)**2+PFLXV(:,:)**2)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
-CALL BL_DEPTH_DIAG(D,ZUSTAR2,PZZ(:,IKB),ZWIND,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_DYN)
+CALL BL_DEPTH_DIAG(D,D%NKT,ZUSTAR2,PZZ(:,IKB),ZWIND,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_DYN)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZSBL_DYN(:) = CSTURB%XSBL_O_BL * ZSBL_DYN(:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
@@ -118,7 +118,7 @@ ZSBL_DYN(:) = CSTURB%XSBL_O_BL * ZSBL_DYN(:)
!
!* BL and SBL diagnosed with buoyancy flux criteria
!
-CALL BL_DEPTH_DIAG(D,ZQ0,PZZ(:,IKB),PWTHV,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_THER)
+CALL BL_DEPTH_DIAG(D,D%NKT,ZQ0,PZZ(:,IKB),PWTHV,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_THER)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZSBL_THER(:)= CSTURB%XSBL_O_BL * ZSBL_THER(:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
diff --git a/src/SURFEX/teb_spartacus.F90 b/src/SURFEX/teb_spartacus.F90
index df4a2b1c96caf396d5fbd06d5183d67808dca19a..b0d9f6a7ff750a491f3c2d02558a2864975db36d 100644
--- a/src/SURFEX/teb_spartacus.F90
+++ b/src/SURFEX/teb_spartacus.F90
@@ -183,9 +183,9 @@ SUBROUTINE TEB_SPARTACUS (TOP, SPAOP, T, B, DMT, GDP, SB, TPN, PDIR_SW, PSCA_SW,
REAL :: ZLAD_MEAN
REAL, DIMENSION(SIZE(SB%XZ,1)) :: ZINTERPOL
!
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_DIRECT_SW
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_DIFFUSE_SW
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_LW
+ REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_DIRECT_SW
+ REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_DIFFUSE_SW
+ REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE :: ZTOP_FLUX_DN_LW
!
REAL :: ZDZ_SP
REAL :: ZH_COUNT