diff --git a/src/LIB/SURCOUCHE/src/mode_field.f90 b/src/LIB/SURCOUCHE/src/mode_field.f90
index 041c3a85b1f926af62708756566e192c72a596ff..d3a1e491f775888a9ca30da15c2290be8e44fabe 100644
--- a/src/LIB/SURCOUCHE/src/mode_field.f90
+++ b/src/LIB/SURCOUCHE/src/mode_field.f90
@@ -14,6 +14,7 @@
! P. Wautelet 06/06/2019: bug correction in FIELDLIST_GOTO_MODEL (XLSTHM was overwritten if LUSERV=.FALSE. due to wrong IF block)
! P. Wautelet 19/06/2019: add Fieldlist_nmodel_resize subroutine + provide KMODEL to INI_FIELD_LIST when known
! P. Wautelet 23/01/2020: split in modd_field.f90 and mode_field.f90
+! JL Redelsperger 03/2021; add variables for Ocean LES and auto-coupled version
!-----------------------------------------------------------------
module mode_field
@@ -1047,6 +1048,32 @@ TFIELDLIST(IDX)%LTIMEDEP = .FALSE.
IDX = IDX+1
!
IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'LOCEAN'
+TFIELDLIST(IDX)%CSTDNAME = ''
+TFIELDLIST(IDX)%CLONGNAME = 'LOCEAN'
+TFIELDLIST(IDX)%CUNITS = ''
+TFIELDLIST(IDX)%CDIR = '--'
+TFIELDLIST(IDX)%CCOMMENT = 'Logical for Ocean MesoNH'
+TFIELDLIST(IDX)%NGRID = 0
+TFIELDLIST(IDX)%NTYPE = TYPELOG
+TFIELDLIST(IDX)%NDIMS = 0
+TFIELDLIST(IDX)%LTIMEDEP = .FALSE.
+IDX = IDX+1
+!
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'LCOUPLES'
+TFIELDLIST(IDX)%CSTDNAME = ''
+TFIELDLIST(IDX)%CLONGNAME = 'LCOUPLES'
+TFIELDLIST(IDX)%CUNITS = ''
+TFIELDLIST(IDX)%CDIR = '--'
+TFIELDLIST(IDX)%CCOMMENT = 'Logical for coupling O-A LES'
+TFIELDLIST(IDX)%NGRID = 0
+TFIELDLIST(IDX)%NTYPE = TYPELOG
+TFIELDLIST(IDX)%NDIMS = 0
+TFIELDLIST(IDX)%LTIMEDEP = .FALSE.
+IDX = IDX+1
+!
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
TFIELDLIST(IDX)%CMNHNAME = 'SURF'
TFIELDLIST(IDX)%CSTDNAME = ''
TFIELDLIST(IDX)%CLONGNAME = 'SURF'
@@ -1269,6 +1296,20 @@ ALLOCATE(TFIELDLIST(IDX)%TFIELD_X3D(IMODEL))
IDX = IDX+1
!
IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'PHIT'
+TFIELDLIST(IDX)%CSTDNAME = 'reduced pressure'
+TFIELDLIST(IDX)%CLONGNAME = 'PHIT'
+TFIELDLIST(IDX)%CUNITS = 'Pa'
+TFIELDLIST(IDX)%CDIR = 'XY'
+TFIELDLIST(IDX)%CCOMMENT = 'X_Y_Z_Reduced Pressure Oce/Shallow conv'
+TFIELDLIST(IDX)%NGRID = 1
+TFIELDLIST(IDX)%NTYPE = TYPEREAL
+TFIELDLIST(IDX)%NDIMS = 3
+TFIELDLIST(IDX)%LTIMEDEP = .TRUE.
+ALLOCATE(TFIELDLIST(IDX)%TFIELD_X3D(IMODEL))
+IDX = IDX+1
+!
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
TFIELDLIST(IDX)%CMNHNAME = 'RT'
TFIELDLIST(IDX)%CSTDNAME = ''
TFIELDLIST(IDX)%CLONGNAME = 'RT'
@@ -3449,6 +3490,31 @@ IDX = IDX+1
!END IF !LFILTERING
END IF !CPROGRAM==REAL .OR. LFICDF
!
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'NFRCLT'
+TFIELDLIST(IDX)%CSTDNAME = ''
+TFIELDLIST(IDX)%CLONGNAME = 'NFRCLT'
+TFIELDLIST(IDX)%CUNITS = 'nb entier de forc'
+TFIELDLIST(IDX)%CDIR = '--'
+TFIELDLIST(IDX)%CCOMMENT = 'nb de flux sfc forcant LES ocean'
+TFIELDLIST(IDX)%NGRID = 0
+TFIELDLIST(IDX)%NTYPE = TYPEINT
+TFIELDLIST(IDX)%NDIMS = 0
+TFIELDLIST(IDX)%LTIMEDEP = .FALSE.
+IDX = IDX+1
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'NINFRT'
+TFIELDLIST(IDX)%CSTDNAME = ''
+TFIELDLIST(IDX)%CLONGNAME = 'NINFRT'
+TFIELDLIST(IDX)%CUNITS = 'nb entier de forc'
+TFIELDLIST(IDX)%CDIR = '--'
+TFIELDLIST(IDX)%CCOMMENT = 'nb de flux sfc forcant LES ocean'
+TFIELDLIST(IDX)%NGRID = 0
+TFIELDLIST(IDX)%NTYPE = TYPEINT
+TFIELDLIST(IDX)%NDIMS = 0
+TFIELDLIST(IDX)%LTIMEDEP = .FALSE.
+IDX = IDX+1
+!
!
WRITE(YMSG,'("number of used fields=",I4," out of ",I4)') IDX-1,MAXFIELDS
CALL PRINT_MSG(NVERB_INFO,'GEN','INI_FIELD_LIST',TRIM(YMSG))
@@ -3724,6 +3790,7 @@ USE MODD_GRID_n
USE MODD_HURR_FIELD_n
USE MODD_LIMA_PRECIP_SCAVENGING_n
USE MODD_LSFIELD_n
+USE MODD_OCEANH
USE MODD_PARAM_n
USE MODD_PAST_FIELD_n
USE MODD_CH_PH_n
@@ -3792,6 +3859,7 @@ CALL FIND_FIELD_ID_FROM_MNHNAME('WT', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(K
CALL FIND_FIELD_ID_FROM_MNHNAME('THT', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XTHT
CALL FIND_FIELD_ID_FROM_MNHNAME('TKET', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XTKET
CALL FIND_FIELD_ID_FROM_MNHNAME('PABST',IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XPABST
+CALL FIND_FIELD_ID_FROM_MNHNAME('PHIT',IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XPHIT
CALL FIND_FIELD_ID_FROM_MNHNAME('RT', IID,IRESP); TFIELDLIST(IID)%TFIELD_X4D(KFROM)%DATA => XRT
!
IF (ASSOCIATED(XRT)) THEN
@@ -3926,7 +3994,7 @@ CALL FIND_FIELD_ID_FROM_MNHNAME('LSUM', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(K
CALL FIND_FIELD_ID_FROM_MNHNAME('LSVM', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XLSVM
CALL FIND_FIELD_ID_FROM_MNHNAME('LSWM', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XLSWM
CALL FIND_FIELD_ID_FROM_MNHNAME('LSTHM',IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XLSTHM
-IF ( LUSERV ) THEN
+IF (LUSERV) THEN
CALL FIND_FIELD_ID_FROM_MNHNAME('LSRVM',IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XLSRVM
END IF
CALL FIND_FIELD_ID_FROM_MNHNAME('LBXUM', IID,IRESP); TFIELDLIST(IID)%TFIELD_X3D(KFROM)%DATA => XLBXUM
@@ -4107,6 +4175,7 @@ CALL FIND_FIELD_ID_FROM_MNHNAME('WT', IID,IRESP); XWT => TFIELDLIST(IID)%TF
CALL FIND_FIELD_ID_FROM_MNHNAME('THT', IID,IRESP); XTHT => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
CALL FIND_FIELD_ID_FROM_MNHNAME('TKET', IID,IRESP); XTKET => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
CALL FIND_FIELD_ID_FROM_MNHNAME('PABST',IID,IRESP); XPABST => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
+CALL FIND_FIELD_ID_FROM_MNHNAME('PHIT', IID,IRESP); XPHIT => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
CALL FIND_FIELD_ID_FROM_MNHNAME('RT', IID,IRESP); XRT => TFIELDLIST(IID)%TFIELD_X4D(KTO)%DATA
!
IF (ASSOCIATED(XRT)) THEN
@@ -4253,7 +4322,7 @@ CALL FIND_FIELD_ID_FROM_MNHNAME('LSUM', IID,IRESP); XLSUM => TFIELDLIST(IID)%TF
CALL FIND_FIELD_ID_FROM_MNHNAME('LSVM', IID,IRESP); XLSVM => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
CALL FIND_FIELD_ID_FROM_MNHNAME('LSWM', IID,IRESP); XLSWM => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
CALL FIND_FIELD_ID_FROM_MNHNAME('LSTHM',IID,IRESP); XLSTHM => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
-IF ( LUSERV ) THEN
+IF (LUSERV) THEN
CALL FIND_FIELD_ID_FROM_MNHNAME('LSRVM',IID,IRESP); XLSRVM => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
END IF
CALL FIND_FIELD_ID_FROM_MNHNAME('LBXUM', IID,IRESP); XLBXUM => TFIELDLIST(IID)%TFIELD_X3D(KTO)%DATA
diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90
index 860afcf7985344f55f85f9b8efc2314ceb658774..31db7726781c590d4174cd1cf51af2cd2243fae2 100644
--- a/src/MNH/bl89.f90
+++ b/src/MNH/bl89.f90
@@ -72,6 +72,7 @@ END MODULE MODI_BL89
!! reversed vertical levels
!! Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
!! 01/2019 (Q. Rodier) support for RM17 mixing length
+!! 03/2021 (JL Redelsperger) Ocean model case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -81,6 +82,7 @@ USE MODD_CONF, ONLY: CPROGRAM
USE MODD_CST
USE MODD_CTURB
USE MODD_PARAMETERS
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_precision, only: MNHREAL
!
!
@@ -177,6 +179,7 @@ ELSE
ZSHEAR (:,JK) = RESHAPE(PSHEAR (:,:,JK),(/ IIU*IJU /) )
ZTKEM (:,JK) = RESHAPE(PTKEM (:,:,JK),(/ IIU*IJU /) )
ZG_O_THVREF(:,JK) = RESHAPE(XG/PTHVREF(:,:,JK),(/ IIU*IJU /) )
+ IF (LOCEAN) ZG_O_THVREF(:,JK) = XG * XALPHAOC
DO JRR=1,KRR
ZRM (:,JK,JRR) = RESHAPE(PRM (:,:,JK,JRR),(/ IIU*IJU /) )
END DO
diff --git a/src/MNH/compute_press_from_oceanbot.f90 b/src/MNH/compute_press_from_oceanbot.f90
new file mode 100644
index 0000000000000000000000000000000000000000..86fe83d26a4061f4247523cb208e2d7012248acb
--- /dev/null
+++ b/src/MNH/compute_press_from_oceanbot.f90
@@ -0,0 +1,232 @@
+!MNH_LIC Copyright 1994-2018 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.
+!-----------------------------------------------------------------
+! #####################################
+ MODULE MODI_COMPUTE_PRESS_FROM_OCEANBOT
+! #####################################
+INTERFACE COMPUTE_PRESS_FROM_OCEANBOT
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D(PRHO,PZFLUX,PSURF2D, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:), INTENT(IN) :: PSURF2D ! bottom pressure
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! press at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! press at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D
+!
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT1D(PRHO,PZFLUX,PSURF, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHO ! virtual potential temperature
+REAL, DIMENSION(:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, INTENT(IN) :: PSURF ! botttom press function
+REAL, DIMENSION(:), INTENT(OUT) :: PFLUX ! press at flux points
+REAL, DIMENSION(:), INTENT(OUT) :: PMASS ! press at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT1D
+!
+END INTERFACE
+END MODULE MODI_COMPUTE_PRESS_FROM_OCEANBOT
+! ######################################
+ MODULE MODI_COMPUTE_PRESS_FROM_OCEANBOT3D
+! ######################################
+INTERFACE
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D(PRHO,PZFLUX,PSURF2D, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:), INTENT(IN) :: PSURF2D! bot press function
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! bot press at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! bot press at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D
+!
+END INTERFACE
+END MODULE MODI_COMPUTE_PRESS_FROM_OCEANBOT3D
+! ########################################################################
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D(PRHO,PZFLUX,PSURF2D,PFLUX,PMASS)
+! ########################################################################
+!
+!!**** *COMPUTE_EXNER_FROM_OCEANBOT3D* - computation of hydrostatic
+!! pressure from ocean bottom
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!! 1 The local pressure is computed at flux points by integration of the hydrostatic
+!! relation from ground (PSURF2D) to top.
+!!
+!! dP= -g rho dz
+!!
+!! 2 The pressure at mass level is computed as follows and linearly
+!! extrapolated for the uppest non-physical level:
+!!
+!! ~ P(k+1)-P(k)
+!! P(k) = -----------------------
+!! lnP(k+1)-lnP(k)
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! function FMLOOK :to retrieve a logical unit number associated with a file
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! Module MODD_CONF : contains configuration variables for all models.
+!! NVERB : verbosity level for output-listing
+!! LTHINSHELL : logical for thinshell approximation
+!! Module MODD_LUNIT : contains logical unit names for all models
+!! CLUOUT0 : name of output-listing
+!! Module MODD_CST : contains physical constants
+ !! XG : gravity constant
+ !! Module MODD_PARAMETERS
+!! JPVEXT,JPHEXT
+!!
+!! REFERENCE
+!! ---------
+!!
+!! Book 2
+!!
+!! AUTHOR
+!! ------
+!!
+!! JLR
+!! MODIFICATIONS
+!! -------------
+!! Original Fev2021
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY : XG
+USE MODD_PARAMETERS, ONLY : JPVEXT
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of arguments
+! ------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:) , INTENT(IN) :: PSURF2D! bottom pressuren
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! pres at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! pres at mass points
+!
+!* 0.2 Declaration of local variables
+! ------------------------------
+INTEGER :: IKB,IKU,JK
+!-------------------------------------------------------------------------------
+!* 1. INITIALIZATIONS
+! ---------------
+IKB=JPVEXT+1
+IKU=SIZE(PZFLUX,3)
+!-------------------------------------------------------------------------------
+!* 2. COMPUTATION OF PRESSURE AT FLUX POINTS
+! ------------------------------------------------
+PFLUX(:,:,IKB)=PSURF2D(:,:)
+!
+ DO JK=IKB+1,IKU
+ PFLUX(:,:,JK)=PFLUX(:,:,JK-1) + XG*PRHO(:,:,JK-1)*(PZFLUX(:,:,JK-1)-PZFLUX(:,:,JK))
+ END DO
+ DO JK=IKB-1,1,-1
+ PFLUX(:,:,JK)=PFLUX(:,:,JK+1) + XG*PRHO(:,:,JK) *(PZFLUX(:,:,JK+1)-PZFLUX(:,:,JK))
+ END DO
+!-------------------------------------------------------------------------------
+!
+!* 3. COMPUTATION OF Pressure AT MASS POINTS
+! --------------------------------------------
+ PMASS(:,:,1:IKU-1)=(PFLUX(:,:,1:IKU-1)+PFLUX(:,:,2:IKU))*.5
+! (PFLUX(:,:,1:IKU-1)-PFLUX(:,:,2:IKU)) &
+! /(LOG(PFLUX(:,:,1:IKU-1))-LOG(PFLUX(:,:,2:IKU)))
+!
+!accute extrapolation not possible as level IKU is in atmosphere. Assume rho_air=1.2
+ PMASS(:,:,IKU)= PMASS(:,:,IKU-1) - XG*1.2 * ( PZFLUX(:,:,IKU)-PZFLUX(:,:,IKU-1) )
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT3D
+! ######################################################################
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT1D(PRHO,PZFLUX,PSURF,PFLUX,PMASS)
+! ######################################################################
+!
+!!**** *COMPUTE_PRESS_FROM_OCEANBOT1D* - computation of hydrostatic press eq
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! Book 2
+!!
+!! AUTHOR
+!! ------
+!!
+!! V.Masson Meteo-France
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 06/03/96
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_COMPUTE_PRESS_FROM_OCEANBOT3D
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of arguments
+! ------------------------
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHO ! virtual potential temperature
+REAL, DIMENSION(:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, INTENT(IN) :: PSURF ! ground Exner function
+REAL, DIMENSION(:), INTENT(OUT) :: PFLUX ! Exner function at flux points
+REAL, DIMENSION(:), INTENT(OUT) :: PMASS ! Exner function at mass points
+!
+!* 0.2 Declaration of local variables
+! ------------------------------
+!
+REAL, DIMENSION(1,1,SIZE(PZFLUX)) :: ZRHO ! virtual potential temperature
+REAL, DIMENSION(1,1,SIZE(PZFLUX)) :: ZZFLUX ! altitude of flux points
+REAL, DIMENSION(1,1) :: ZPSURF ! ground Exner function
+REAL, DIMENSION(1,1,SIZE(PZFLUX)) :: ZPFLUX ! Exner function at flux points
+REAL, DIMENSION(1,1,SIZE(PZFLUX)) :: ZPMASS ! Exner function at mass points
+!
+!-------------------------------------------------------------------------------
+!
+ZRHO(1,1,:)=PRHO(:)
+ZZFLUX(1,1,:)=PZFLUX(:)
+ZPSURF(1,1)=PSURF
+!
+CALL COMPUTE_PRESS_FROM_OCEANBOT3D(ZRHO,ZZFLUX,ZPSURF,ZPFLUX,ZPMASS)
+!
+PFLUX(:)=ZPFLUX(1,1,:)
+PMASS(:)=ZPMASS(1,1,:)
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANBOT1D
+
diff --git a/src/MNH/compute_press_from_oceansfc.f90 b/src/MNH/compute_press_from_oceansfc.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ba8f86c200a06f7b3b879acaaceafca30e9b700c
--- /dev/null
+++ b/src/MNH/compute_press_from_oceansfc.f90
@@ -0,0 +1,238 @@
+!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.
+!-----------------------------------------------------------------
+! ##################################
+ MODULE MODI_COMPUTE_PRESS_FROM_OCEANSFC
+! ##################################
+INTERFACE COMPUTE_PRESS_FROM_OCEANSFC
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D(PRHO,PZFLUX,PTOP2D, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:) , INTENT(IN) :: PTOP2D ! 2D Pressure at domain top
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! Pressure at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! Pressure at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D
+!
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC1D(PRHO,PZFLUX,PTOP, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, INTENT(IN) :: PTOP ! Pressure at domain top
+REAL, DIMENSION(:), INTENT(OUT) :: PFLUX ! Pressure at flux points
+REAL, DIMENSION(:), INTENT(OUT) :: PMASS ! Pressure at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC1D
+END INTERFACE
+END MODULE MODI_COMPUTE_PRESS_FROM_OCEANSFC
+! ####################################
+ MODULE MODI_COMPUTE_PRESS_FROM_OCEANSFC3D
+! ####################################
+INTERFACE
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D(PRHO,PZFLUX,PTOP2D, &
+ PFLUX,PMASS)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:) , INTENT(IN) :: PTOP2D ! Pressure at top domain
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! Pressure at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! Pressure at mass points
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D
+END INTERFACE
+END MODULE MODI_COMPUTE_PRESS_FROM_OCEANSFC3D
+! ####################################################################
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D(PRHO,PZFLUX,PTOP2D,PFLUX,PMASS)
+! ####################################################################
+!
+!!**** *COMPUTE_PRESS_FROM_OCEANSFC3D* - computation of hydrostatic
+!! from model top
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!! 1 The local Exner function in computed by integration of the hydrostatic
+!! relation from top (PEXNTOP2D) to bottom.
+!!
+!! dPI= -g/(Cpd thetav) dz
+!!
+!! 2 The Exner function at mass level is computed as follows and linearly
+!! extrapolated for the uppest non-physical level:
+!!
+!! ~ PI(k+1)-PI(k)
+!! PI(k) = -----------------------
+!! lnPI(k+1)-lnPI(k)
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! Module MODD_CONF : contains configuration variables for all models.
+!! NVERB : verbosity level for output-listing
+!! LTHINSHELL : logical for thinshell approximation
+!! Module MODD_LUNIT : contains logical unit names for all models
+!! CLUOUT0 : name of output-listing
+!! Module MODD_CST : contains physical constants
+!! XG : gravity constant
+!! XCPD: specific heat for dry air at constant pressure
+!! XRD : gas constant for dry air
+!! Module MODD_PARAMETERS
+!! JPVEXT,JPHEXT
+!!
+!! REFERENCE
+!! ---------
+!!
+!! Book 2
+!!
+!! AUTHOR
+!! ------
+!! JLR
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 01/01/21
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY : XG
+USE MODD_PARAMETERS, ONLY : JPVEXT
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of arguments
+! ------------------------
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, DIMENSION(:,:) , INTENT(IN) :: PTOP2D ! Pressure at domain
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFLUX ! Pressure at flux points
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMASS ! Pressure at mass points
+!
+!* 0.2 Declaration of local variables
+! ------------------------------
+!
+INTEGER :: IKB,IKE,IKU,JK
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+IKB=JPVEXT+1
+IKE=SIZE(PZFLUX,3)-JPVEXT
+IKU=SIZE(PZFLUX,3)
+!-------------------------------------------------------------------------------
+!* 2. COMPUTATION OF PRESSURE AT FLUX POINTS
+! ------------------------------------------------
+PFLUX(:,:,IKE+1)=PTOP2D(:,:)
+!
+ DO JK=IKE,1,-1
+ PFLUX(:,:,JK) = PFLUX(:,:,JK+1) + XG*PRHO(:,:,JK)*(PZFLUX(:,:,JK+1)-PZFLUX(:,:,JK))
+ END DO
+ DO JK=IKE+2,SIZE(PFLUX,3)
+ PFLUX(:,:,JK) = PFLUX(:,:,JK-1) + XG*PRHO(:,:,JK-1)*(PZFLUX(:,:,JK-1)-PZFLUX(:,:,JK))
+ END DO
+!-------------------------------------------------------------------------------
+!
+!* 3. COMPUTATION OF PRESSURE AT MASS POINTS
+! --------------------------------------------
+ PMASS(:,:,1:IKU-1)=(PFLUX(:,:,1:IKU-1)-PFLUX(:,:,2:IKU)) &
+ /(LOG(PFLUX(:,:,1:IKU-1))-LOG(PFLUX(:,:,2:IKU)))
+!accute extrapolation not possible as level IKU is in atmosphere. Assume rho_air=1.2
+ PMASS(:,:,IKU)= PMASS(:,:,IKU-1) - XG*1.2 * ( PZFLUX(:,:,IKU)-PZFLUX(:,:,IKU-1) )
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC3D
+! ##################################################################
+ SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC1D(PRHOD,PZFLUX,PTOP,PFLUX,PMASS)
+! ##################################################################
+!
+!!**** *COMPUTE_PRESS_FROM_OCEANSFC1D* - computation of hydrostatic
+!! Pressure from model top
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! Book 2
+!!
+!! AUTHOR
+!! ------
+!!
+!! JLR
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 01/01/21
+!-----------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_COMPUTE_PRESS_FROM_OCEANSFC3D
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of arguments
+! ------------------------
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHOD ! virtual potential temperature
+REAL, DIMENSION(:), INTENT(IN) :: PZFLUX ! altitude of flux points
+REAL, INTENT(IN) :: PTOP ! top pressure
+REAL, DIMENSION(:), INTENT(OUT) :: PFLUX ! Pressure at flux points
+REAL, DIMENSION(:), INTENT(OUT) :: PMASS ! Pressure at mass points
+!
+!* 0.2 Declaration of local variables
+! ------------------------------
+!
+!JUAN
+REAL, DIMENSION(3,1,SIZE(PZFLUX)) :: ZRHOD ! virtual potential temperature
+REAL, DIMENSION(3,1,SIZE(PZFLUX)) :: ZZFLUX ! altitude of flux points
+REAL, DIMENSION(3,1) :: ZPTOP ! top pressure
+REAL, DIMENSION(3,1,SIZE(PZFLUX)) :: ZPFLUX ! Pressure at flux points
+REAL, DIMENSION(3,1,SIZE(PZFLUX)) :: ZPMASS ! Pressure at mass points
+INTEGER :: JI ! loop index in I
+!JUAN
+!
+!-------------------------------------------------------------------------------
+!
+!JUAN
+DO JI=1,3
+ZRHOD(JI,1,:)=PRHOD(:)
+ZZFLUX(JI,1,:)=PZFLUX(:)
+ZPTOP(JI,1)=PTOP
+END DO
+!JUAN
+
+!
+CALL COMPUTE_PRESS_FROM_OCEANSFC3D(ZRHOD,ZZFLUX,ZPTOP,ZPFLUX,ZPMASS)
+!
+PFLUX(:)=ZPFLUX(2,1,:)
+PMASS(:)=ZPMASS(2,1,:)
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE COMPUTE_PRESS_FROM_OCEANSFC1D
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index f39f694235ad68d13968bc9bcef562f05e42aa95..2e9dd5730e72de76434299627dfbeabe0f4390c2 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -213,6 +213,7 @@ END MODULE MODI_DEFAULT_DESFM_n
! P-A Joulin 21/05/2021: add Wind turbines
! S. Riette 21/05/2021: add options to PDF subgrid scheme
! 05/2021 D. Ricard add the contribution of Leonard terms in the turbulence scheme
+!! JL Redelsperger 06/2021 add parameters allowing to active idealized oceanic convection
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -360,6 +361,7 @@ LUSERI = .FALSE.
LUSERS = .FALSE.
LUSERG = .FALSE.
LUSERH = .FALSE.
+LOCEAN = .FALSE.
!NSV = 0
!NSV_USER = 0
LUSECI = .FALSE.
@@ -798,6 +800,11 @@ IF (KMI == 1) THEN
XUTRANS = 0.0
XVTRANS = 0.0
LPGROUND_FRC = .FALSE.
+ LDEEPOC = .FALSE.
+ XCENTX_OC = 16000.
+ XCENTY_OC = 16000.
+ XRADX_OC = 8000.
+ XRADY_OC = 8000.
END IF
!
!-------------------------------------------------------------------------------
diff --git a/src/MNH/default_expre.f90 b/src/MNH/default_expre.f90
index 6c158a6c4e3d218dbae5664a8d1b8ebc00f8a0bf..f592e5205520b06147b8a117dab3c6f4b1d9bca7 100644
--- a/src/MNH/default_expre.f90
+++ b/src/MNH/default_expre.f90
@@ -95,6 +95,7 @@ END MODULE MODI_DEFAULT_EXPRE
!! add the uniform soil values 05/02/96 (J.Stein)
!! removes default values for ground variables 26/11/96 (V.Masson)
!! add default value for LBOUSS 11/07/13 (C.Lac)
+!! add default value LOCEAN LCOUPLES /03/21 (JL Redelsperger)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -103,6 +104,7 @@ USE MODD_CONF ! declarative modules
USE MODD_DIM_n
USE MODD_GRID
USE MODD_REF
+USE MODD_DYN_n, ONLY : LOCEAN
!
IMPLICIT NONE
!
@@ -140,7 +142,9 @@ XLATORI = 37.
!* 4. SET DEFAULT VALUES FOR MODD_REF :
! --------------------------------
!
-LBOUSS = .FALSE.
+LBOUSS = .FALSE.
+LOCEAN = .FALSE.
+LCOUPLES= .FALSE.
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/dyn_sources.f90 b/src/MNH/dyn_sources.f90
index efc0e518999897364ba4ef04a6a3461289d00189..e79ff37a97716739a3be8b48dddd54c033bc5c9a 100644
--- a/src/MNH/dyn_sources.f90
+++ b/src/MNH/dyn_sources.f90
@@ -158,7 +158,8 @@ use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, &
USE MODD_CONF
USE MODD_CST
USE MODD_DYN
-
+USE MODD_DYN_n, ONLY : LOCEAN
+!
use mode_budget, only: Budget_store_init, Budget_store_end
USE MODE_MPPDB
@@ -319,6 +320,7 @@ ELSE
ENDIF
!
IF( .NOT.L1D ) THEN
+ IF (.NOT. LOCEAN) THEN
!
IF (KRR > 0) THEN
if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
@@ -346,6 +348,7 @@ IF( .NOT.L1D ) THEN
if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
END IF
+ END IF
!
END IF
!
diff --git a/src/MNH/emoist.f90 b/src/MNH/emoist.f90
index f58a1b32d2f6d7901b6c25f61edad0a0ddde339f..99a5e1ce0b994bdc27aa85e86b93a75ad21b1ea0 100644
--- a/src/MNH/emoist.f90
+++ b/src/MNH/emoist.f90
@@ -80,12 +80,14 @@ FUNCTION EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM) RESULT(PEMOIST)
!! J. Stein Feb 28, 1996 optimization + Doctorization
!! J. Stein Spet 15, 1996 Amoist previously computed
!! J.-P. Pinty May 20, 2003 Improve EMOIST expression
+!! 03/2021 (JL Redelsperger) Ocean model case
!!
!! ----------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
!
IMPLICIT NONE
!
@@ -118,14 +120,21 @@ INTEGER :: JRR ! moist loop counter
!
!* 1. COMPUTE EMOIST
! --------------
-!
-!
-IF ( KRR == 0 ) THEN ! dry case
- PEMOIST(:,:,:) = 0.
-ELSE IF ( KRR == 1 ) THEN ! only vapor
+IF (LOCEAN) THEN
+ IF ( KRR == 0 ) THEN ! Unsalted
+ PEMOIST(:,:,:) = 0.
+ ELSE
+ PEMOIST(:,:,:) = 1. ! Salted case
+ END IF
+!
+ELSE
+!
+ IF ( KRR == 0 ) THEN ! dry case
+ PEMOIST(:,:,:) = 0.
+ ELSE IF ( KRR == 1 ) THEN ! only vapor
ZDELTA = (XRV/XRD) - 1.
PEMOIST(:,:,:) = ZDELTA*PTHLM(:,:,:)
-ELSE ! liquid water & ice present
+ ELSE ! liquid water & ice present
ZDELTA = (XRV/XRD) - 1.
ZRW(:,:,:) = PRM(:,:,:,1)
!
@@ -173,8 +182,9 @@ ELSE ! liquid water & ice present
/ (1. + ZRW(:,:,:)) &
) * PAMOIST(:,:,:) * 2. * PSRCM(:,:,:)
END IF
-END IF
+ END IF
!
+END IF
!---------------------------------------------------------------------------
!
END FUNCTION EMOIST
diff --git a/src/MNH/etheta.f90 b/src/MNH/etheta.f90
index 6c673cb978262ffe9cdb70251f4943f8b828cb0b..37f39873d08a4a5a37e032642d0ec36752909676 100644
--- a/src/MNH/etheta.f90
+++ b/src/MNH/etheta.f90
@@ -82,12 +82,13 @@ FUNCTION ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM) RESULT(PETHETA)
!! J. Stein Feb 28, 1996 optimization + Doctorization
!! J. Stein Sept 15, 1996 Atheta previously computed
!! J.-P. Pinty May 20, 2003 Improve ETHETA expression
-!!
+!! J.L Redelsperger 03, 2021 Ocean Model Case
!! ----------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
!
IMPLICIT NONE
!
@@ -125,12 +126,15 @@ INTEGER :: JRR ! moist loop counter
! --------------
!
!
-IF ( KRR == 0 ) THEN ! dry case
+IF (LOCEAN) THEN ! ocean case
+ PETHETA(:,:,:) = 1.
+ELSE
+ IF ( KRR == 0.) THEN ! dry case
PETHETA(:,:,:) = 1.
-ELSE IF ( KRR == 1 ) THEN ! only vapor
+ ELSE IF ( KRR == 1 ) THEN ! only vapor
ZDELTA = (XRV/XRD) - 1.
PETHETA(:,:,:) = 1. + ZDELTA*PRM(:,:,:,1)
-ELSE ! liquid water & ice present
+ ELSE ! liquid water & ice present
ZDELTA = (XRV/XRD) - 1.
ZRW(:,:,:) = PRM(:,:,:,1)
!
@@ -173,8 +177,9 @@ ELSE ! liquid water & ice present
/ (1. + ZRW(:,:,:)) &
) * PATHETA(:,:,:) * 2. * PSRCM(:,:,:)
END IF
-END IF
+ END IF
!
+END IF
!---------------------------------------------------------------------------
!
END FUNCTION ETHETA
diff --git a/src/MNH/gravity.f90 b/src/MNH/gravity.f90
index 72d0c1649e522a6b44541baa3b2ace50d7147061..27c8c3461514a02d0757bafaf6e5b7044bf1b447 100644
--- a/src/MNH/gravity.f90
+++ b/src/MNH/gravity.f90
@@ -100,7 +100,7 @@ END MODULE MODI_GRAVITY
!! C.Lac - March 2011 - Splitted from dyn_sources
!! Q.Rodier 06/15 correction on budget
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
-!!
+!! J.L. Redelsperger 03/2021 : Ocean model case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -108,6 +108,8 @@ END MODULE MODI_GRAVITY
!
USE MODD_CONF
USE MODD_CST
+USE MODD_REF
+USE MODD_DYN_n, ONLY : LOCEAN
!
USE MODI_SHUMAN
USE MODI_GET_HALO
@@ -144,6 +146,16 @@ REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: &
!
IF( .NOT.L1D ) THEN ! no buoyancy for 1D case
!
+ IF (LOCEAN) THEN !ocean case
+ CALL GET_HALO(PTHT)
+ IF(KRR > 0) THEN
+ CALL GET_HALO(PRT(:,:,:,1))
+ PRWS(:,:,:) = PRWS + XG * (XALPHAOC*MZM((PTHT - PTHVREF )*PRHODJ) &
+ - XBETAOC*MZM((PRT(:,:,:,1) - XSA00OCEAN)*PRHODJ) )
+ ELSE ! unsalted case
+ PRWS(:,:,:) = PRWS + XG * XALPHAOC*MZM((PTHT - PTHVREF )*PRHODJ )
+ END IF
+ELSE ! Atmospheric case
IF(KRR > 0) THEN
!
! compute the ratio : 1 + total water mass / dry air mass
@@ -156,7 +168,6 @@ IF( .NOT.L1D ) THEN ! no buoyancy for 1D case
END DO
!
! compute the virtual potential temperature when water is present in any form
-!
CALL GET_HALO(PTHT)
!
@@ -171,6 +182,7 @@ IF( .NOT.L1D ) THEN ! no buoyancy for 1D case
! compute the gravity term
!
PRWS(:,:,:) = PRWS + XG * MZM( ( (ZWORK2/PTHVREF) - 1. ) * PRHODJ )
+ END IF
!
! the extrapolation for the PTHT and the THVREF must be the same at the
! ground
diff --git a/src/MNH/ini_cst.f90 b/src/MNH/ini_cst.f90
index 8ce239471ca9e449b35764ceb4c08f93e2ff3f8e..6c266aac10677013c673f9ea6db4aaff4096c59e 100644
--- a/src/MNH/ini_cst.f90
+++ b/src/MNH/ini_cst.f90
@@ -112,6 +112,13 @@ XG = 9.80665
!* 4. REFERENCE PRESSURE
! -------------------
!
+! Ocean model cst same as in 1D/CMO SURFEX
+! values used in init_cst to overwrite XP00 and XTH00
+XRH00OCEAN =1024.
+XTH00OCEAN = 286.65
+XSA00OCEAN= 32.6
+XP00OCEAN = 201.E5
+!Atmospheric model
XP00 = 1.E5
XTH00 = 300.
!-------------------------------------------------------------------------------
@@ -151,7 +158,11 @@ XALPW = LOG(XESTT) + (XBETAW /XTT) + (XGAMW *LOG(XTT))
XGAMI = (XCI - XCPV) / XRV
XBETAI = (XLSTT/XRV) + (XGAMI * XTT)
XALPI = LOG(XESTT) + (XBETAI /XTT) + (XGAMI *LOG(XTT))
-
+! Values identical to ones used in CMO1D in SURFEX /could be modified
+! Coefficient of thermal expansion of water (K-1)
+XALPHAOC = 1.9E-4
+! Coeff of Haline contraction coeff (S-1)
+XBETAOC= 7.7475E-4
!
! Some machine precision value depending of real4/8 use
!
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index ef474931a1d41fac0bdb41bd9cc2506cce4976f3..e7dad9f4e8dfd044f03aef7edb721927b17abee6 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -292,6 +292,7 @@ END MODULE MODI_INI_MODEL_n
! S. Riette 04/2020: XHL* fields
! F. Auguste 02/2021: add IBM
! T.Nigel 02/2021: add turbulence recycling
+! J.L.Redelsperger 06/2011: OCEAN case
!---------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -370,6 +371,7 @@ USE MODD_NESTING, only: CDAD_NAME, NDAD, NDT_2_WAY, NDTRATIO, NDXRATIO
USE MODD_NSV
USE MODD_NSV
USE MODD_NUDGING_n, only: LNUDGING
+USE MODD_OCEANH
USE MODD_OUT_n
USE MODD_PARAMETERS
USE MODD_PARAM_KAFR_n
@@ -528,7 +530,7 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZIBM_LS ! LevelSet IBM
!
INTEGER, DIMENSION(:,:),ALLOCATABLE :: IINDEX ! indices of non-zero terms
INTEGER, DIMENSION(:),ALLOCATABLE :: IIND
-INTEGER :: JM
+INTEGER :: JM, JT
!
!------------------------------------------
! Dummy pointers needed to correct an ifort Bug
@@ -971,11 +973,19 @@ ALLOCATE(XDZZ(IIU,IJU,IKU))
!
!* 3.3 Modules MODD_REF and MODD_REF_n
!
-IF (KMI == 1) THEN
+! different reference state for O and A models
+! POur le moment meme etat de ref O et A
+!IF ((KMI == 1).OR.LCOUPLES) THEN
+IF (KMI==1) THEN
ALLOCATE(XRHODREFZ(IKU),XTHVREFZ(IKU))
+ELSE IF (LCOUPLES) THEN
+! in coupled O-A case, need different variables for ocean
+ ALLOCATE(XRHODREFZO(IKU),XTHVREFZO(IKU))
ELSE
!Do not allocate XRHODREFZ and XTHVREFZ because they are the same on all grids (not 'n' variables)
END IF
+!
+ALLOCATE(XPHIT(IIU,IJU,IKU))
ALLOCATE(XRHODREF(IIU,IJU,IKU))
ALLOCATE(XTHVREF(IIU,IJU,IKU))
ALLOCATE(XEXNREF(IIU,IJU,IKU))
@@ -1830,7 +1840,7 @@ IF (CCLOUD=='LIMA') CALL INIT_AEROSOL_PROPERTIES
! --------------------------------
!
CALL MPPDB_CHECK3D(XUT,"INI_MODEL_N-before read_field::XUT",PRECISION)
-CALL READ_FIELD(TPINIFILE,IIU,IJU,IKU, &
+CALL READ_FIELD(KMI,TPINIFILE,IIU,IJU,IKU, &
CGETTKET,CGETRVT,CGETRCT,CGETRRT,CGETRIT,CGETCIT,CGETZWS, &
CGETRST,CGETRGT,CGETRHT,CGETSVT,CGETSRCT,CGETSIGS,CGETCLDFR, &
CGETBL_DEPTH,CGETSBL_DEPTH,CGETPHC,CGETPHR,CUVW_ADV_SCHEME, &
@@ -2229,10 +2239,14 @@ ALLOCATE(ZSCA_ALB(IIU,IJU,NSWB_MNH))
ALLOCATE(ZEMIS (IIU,IJU,NLWB_MNH))
ALLOCATE(ZTSRAD (IIU,IJU))
!
-IF ((TPINIFILE%NMNHVERSION(1)==4 .AND. TPINIFILE%NMNHVERSION(2)>=6) .OR. TPINIFILE%NMNHVERSION(1)>4) THEN
- CALL IO_Field_read(TPINIFILE,'SURF',CSURF)
+IF (LCOUPLES.AND.(KMI>1))THEN
+ CSURF ="NONE"
ELSE
- CSURF = "EXTE"
+ IF ((TPINIFILE%NMNHVERSION(1)==4 .AND. TPINIFILE%NMNHVERSION(2)>=6) .OR. TPINIFILE%NMNHVERSION(1)>4) THEN
+ CALL IO_Field_read(TPINIFILE,'SURF',CSURF)
+ ELSE
+ CSURF = "EXTE"
+ END IF
END IF
!
!
@@ -2601,5 +2615,19 @@ IF (LMAIN_EOL .AND. KMI == NMODEL_EOL) THEN
END SELECT
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))
+END IF
+!
END SUBROUTINE INI_MODEL_n
diff --git a/src/MNH/ini_one_wayn.f90 b/src/MNH/ini_one_wayn.f90
index fffb5f8d6f5bd2a47ea253c49cd177eb765a4e8b..49932a2494dab21ec0048fb1f6f6c81090bfcef5 100644
--- a/src/MNH/ini_one_wayn.f90
+++ b/src/MNH/ini_one_wayn.f90
@@ -90,7 +90,7 @@ SUBROUTINE INI_ONE_WAY_n(KDAD,KMI, &
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
! P. Wautelet 03/05/2019: restructuration of one_wayn and ini_one_wayn
! P. Wautelet 04/06/2020: correct call to Set_conc_lima + initialize ZCONCM
-!
+! J-L Redelsperger 06/2021: add Ocean coupling
!------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -106,6 +106,7 @@ USE MODD_NSV, only: NSV_A, NSV_C1R3BEG_A, NSV_C1R3_A, NSV_C2R2BEG_A,
NSV_SLTBEG_A, NSV_SLTDEPBEG_A, NSV_SLTDEP_A, NSV_SLT_A, NSV_USER_A
USE MODD_PARAM_n, only: CCLOUD
+USE MODD_REF, ONLY: LCOUPLES
USE MODD_REF_n, only: XRHODJ, XRHODREF
!
use mode_bikhardt
@@ -183,6 +184,20 @@ REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCHEMM ! chemical concentrations
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCHEMMI ! chemical ice phase concentrations
!-------------------------------------------------------------------------------
!
+IF (LCOUPLES) THEN
+ PLBXUM=0.
+ PLBXVM=0.
+ PLBXWM=0.
+ PLBXTHM=0.
+ PLBYTHM=0.
+ PLBXTKEM=0.
+ PLBYTKEM =0.
+ PLBXRM =0.
+ PLBYRM=0.
+ PLBXSVM =0.
+ PLBYSVM=0.
+RETURN
+ENDIF
!* 0. INITIALISATION
!
CALL GOTO_MODEL(KDAD)
diff --git a/src/MNH/ini_segn.f90 b/src/MNH/ini_segn.f90
index e0819946a423f1d3026e954ffd933dadeafd362f..b852d82f0f227c7b6c6217967fd71af257636fb3 100644
--- a/src/MNH/ini_segn.f90
+++ b/src/MNH/ini_segn.f90
@@ -172,6 +172,7 @@ END MODULE MODI_INI_SEG_n
USE MODD_CONF
USE MODD_CONF_n, ONLY: CSTORAGE_TYPE
USE MODN_CONFZ
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_DYN
USE MODD_IO, ONLY: NVERB_FATAL, NVERB_WARNING, TFILE_OUTPUTLISTING, TFILEDATA
USE MODD_LUNIT
@@ -440,6 +441,8 @@ END IF
! routine which read related informations in the EXSEG descriptor in order to
! check coherence between both informations.
!
+CALL IO_Field_read(TPINIFILE,'LOCEAN',LOCEAN)
+!
CALL READ_EXSEG_n(KMI,TZFILE_DES,YCONF,GFLAT,GUSERV,GUSERC, &
GUSERR,GUSERI,GUSECI,GUSERS,GUSERG,GUSERH,GUSECHEM, &
GUSECHAQ,GUSECHIC,GCH_PH, &
diff --git a/src/MNH/ini_sizen.f90 b/src/MNH/ini_sizen.f90
index 5e21aeeaffa449c582206864c825cfdfaaa846be..caad447cacd9e7ecdeab6d7253e5b16c52dfcf29 100644
--- a/src/MNH/ini_sizen.f90
+++ b/src/MNH/ini_sizen.f90
@@ -112,6 +112,7 @@ USE MODD_LUNIT_n, ONLY: CINIFILE, CINIFILEPGD, TLUOUT
USE MODD_NESTING, ONLY: CMY_NAME, CDAD_NAME, NDAD, NDXRATIO_ALL, NDYRATIO_ALL, &
NXOR_ALL, NYOR_ALL, NXEND_ALL,NYEND_ALL
USE MODD_PARAMETERS, ONLY: JPMODELMAX, JPHEXT,JPVEXT
+USE MODD_REF, ONLY: LCOUPLES
!
USE MODE_IO, ONLY: IO_Pack_set
USE MODE_IO_FIELD_READ, only: IO_Field_read
@@ -164,6 +165,12 @@ IF (IRESP /= 0) THEN
END IF
!
IF ( KMI > 1 ) THEN
+ IF (LCOUPLES) THEN
+ IF (KMI==2) THEN
+ CMY_NAME(NDAD(KMI)) = CDAD_NAME(KMI)
+ WRITE(UNIT=ILUOUT,FMT=*) 'NDAD',NDAD(KMI),'changed in '//TRIM(CMY_NAME(NDAD(KMI)))//TRIM(CDAD_NAME(KMI)),KMI
+ END IF
+ END IF
IF ( TRIM(CDAD_NAME(KMI)) /= TRIM(CMY_NAME(NDAD(KMI))) ) THEN
WRITE(UNIT=ILUOUT,FMT=9005) NDAD(KMI)
WRITE(ILUOUT,FMT=*) ' THE INITIAL FM-File IS NOT CONSISTANT WITH THE ONE OF THE DAD MODEL!'
diff --git a/src/MNH/ini_tke_eps.f90 b/src/MNH/ini_tke_eps.f90
index c76c795b5e67772a14cd56384910110c2a55a010..5bdd3d5755ae688cf9d4de32d04a0d9f1844fd28 100644
--- a/src/MNH/ini_tke_eps.f90
+++ b/src/MNH/ini_tke_eps.f90
@@ -84,18 +84,20 @@ END MODULE MODI_INI_TKE_EPS
!! Aug 10, 1998 (N. Asencio) add parallel code
!! May 2006 Remove KEPS
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! March 2021 (JL Redelsperger) Add Ocean LES case)
!! -------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
-USE MODD_CTURB ! XLINI, XCED, XCMFS, XTKEMIN
-USE MODD_CST ! XG, XRD, XRV
-USE MODD_PARAMETERS ! JPVEXT
+USE MODD_CTURB, ONLY : XLINI, XCED, XCMFS, XTKEMIN, XCSHF
+USE MODD_CST, ONLY : XG, XRD, XRV, XALPHAOC
+USE MODD_PARAMETERS, ONLY : JPVEXT
+USE MODD_DYN_n, ONLY : LOCEAN
!
-USE MODI_SHUMAN ! DZF, MXF, MYF, MZM
+USE MODI_SHUMAN, ONLY : DZF, MXF, MYF, MZM
USE MODE_ll
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODD_ARGSLIST_ll,ONLY : LIST_ll
!
IMPLICIT NONE
!
@@ -148,11 +150,20 @@ IF (HGETTKET == 'INIT' ) THEN
PVT(:,:,IKE+1) = PVT(:,:,IKE)
!
! determines TKE
+ ! Equilibrium/Stationary/neutral 1D TKE equation
+IF (LOCEAN) THEN
+ PTKET(:,:,:)=(XLINI**2/XCED)*( &
+ XCMFS*( DZF(MXF(MZM(PUT)))**2 &
+ +DZF(MYF(MZM(PVT)))**2) / ZDELTZ &
+ -(XG*XALPHAOC)*XCSHF*DZF(MZM(PTHT)) &
+ ) / ZDELTZ
+ELSE
PTKET(:,:,:)=(XLINI**2/XCED)*( &
XCMFS*( DZF(MXF(MZM(PUT)))**2 &
+DZF(MYF(MZM(PVT)))**2) / ZDELTZ &
-(XG/PTHVREF)*XCSHF*DZF(MZM(PTHT)) &
) / ZDELTZ
+END IF
! positivity control
WHERE (PTKET < XTKEMIN) PTKET=XTKEMIN
!
diff --git a/src/MNH/init_mnh.f90 b/src/MNH/init_mnh.f90
index 7b7d49a3aa77d54df4e93992845fcadb816821c2..e275487920bccbf19c22ec93e4dc0e1b880f9279 100644
--- a/src/MNH/init_mnh.f90
+++ b/src/MNH/init_mnh.f90
@@ -76,13 +76,15 @@
!* 0. DECLARATIONS
! ------------
USE MODD_CONF
-USE MODD_DYN_n, ONLY: CPRESOPT,NITR ! only for spawning purpose
+USE MODD_CST, ONLY:XP00,XTH00,XP00OCEAN,XTH00OCEAN
+USE MODD_DYN_n, ONLY: CPRESOPT,NITR,LOCEAN ! only for spawning purpose
USE MODD_IO, ONLY: TFILE_OUTPUTLISTING, TPTR2FILE
USE MODD_LBC_n, ONLY: CLBCX,CLBCY ! only for spawning purpose
USE MODD_LUNIT
USE MODD_LUNIT_n
USE MODD_MNH_SURFEX_n
USE MODD_PARAMETERS
+USE MODD_REF
!
use mode_field, only: Alloc_field_scalars, Fieldlist_goto_model
USE MODE_IO_FILE, ONLY: IO_File_open
@@ -180,6 +182,13 @@ END IF
!
IF (CPROGRAM=='DIAG') CALL RESET_EXSEG()
!
+! UPDATE CONSTANTS FOR OCEAN MODEL
+DO JMI=1,JPMODELMAX
+ IF (LOCEAN) THEN
+ XP00=XP00OCEAN
+ XTH00=XTH00OCEAN
+ END IF
+END DO
!-------------------------------------------------------------------------------
!
!
diff --git a/src/MNH/modd_cst.f90 b/src/MNH/modd_cst.f90
index f6bc6c52f04bc35604a71222dbf354c17bda4d84..6832448e4008da354f7372b48bd03b6225af5a99 100644
--- a/src/MNH/modd_cst.f90
+++ b/src/MNH/modd_cst.f90
@@ -42,6 +42,7 @@
!! C. Mari 31/10/00 add NDAYSEC
!! V. Masson 01/03/03 add conductivity of ice
!! J.Escobar : 10/2017 : for real*4 , add XMNH_HUGE_12_LOG
+!! J.L. Redelsperger 03/2021 add constants for ocean penetrating solar
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -63,6 +64,8 @@ REAL,SAVE :: XRADIUS,XOMEGA ! Earth radius, earth rotation
REAL,SAVE :: XG ! Gravity constant
!
REAL,SAVE :: XP00 ! Reference pressure
+REAL,SAVE :: XP00OCEAN ! Reference pressurefor ocean model
+REAL,SAVE :: XRH00OCEAN ! Reference density for ocean model
!
REAL,SAVE :: XSTEFAN,XI0 ! Stefan-Boltzman constant, solar constant
!
@@ -83,8 +86,19 @@ REAL,SAVE :: XALPW,XBETAW,XGAMW ! Constants for saturation vapor
REAL,SAVE :: XALPI,XBETAI,XGAMI ! Constants for saturation vapor
! pressure function over solid ice
REAL,SAVE :: XCONDI ! thermal conductivity of ice (W m-1 K-1)
-REAL, SAVE :: XTH00 ! reference value for the potential
- ! temperature
+REAL,SAVE :: XALPHAOC ! thermal expansion coefficient for ocean (K-1)
+REAL,SAVE :: XBETAOC ! Haline contraction coeff for ocean (S-1)
+REAL,SAVE :: XTH00 ! reference value for the potential temperature
+REAL,SAVE :: XTH00OCEAN ! Ref value for pot temp in ocean model
+REAL,SAVE :: XSA00OCEAN ! Ref value for SAlinity in ocean model
+REAL,SAVE :: XROC=0.69! 3 coeffs for SW penetration in Ocean (Hoecker et al)
+REAL,SAVE :: XD1=1.1
+REAL,SAVE :: XD2=23.
+! Values used in SURFEX CMO
+!REAL,SAVE :: XROC=0.58
+!REAL,SAVE :: XD1=0.35
+!REAL,SAVE :: XD2=23.
+
REAL,SAVE :: XRHOLI ! Volumic mass of liquid water
!
INTEGER, SAVE :: NDAYSEC ! Number of seconds in a day
diff --git a/src/MNH/modd_dynn.f90 b/src/MNH/modd_dynn.f90
index c3b64335c44521991cc06f96a9825de59183ca7c..60d8cfc76f9bea08cd3737a559ae18ddd371f3d5 100644
--- a/src/MNH/modd_dynn.f90
+++ b/src/MNH/modd_dynn.f90
@@ -43,6 +43,7 @@
!! Modification 01/2016 (JP Pinty) Add LIMA
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! Modification 07/2017 (V. Vionnet) Add blowing snow variable
+!! Modification 03/2021 (JL Redelsperger) Add logical LOCEAN
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -55,6 +56,7 @@ TYPE DYN_t
!
INTEGER :: NSTOP ! Number of time step
REAL :: XTSTEP ! Time step
+ LOGICAL :: LOCEAN
!
!++++++++++++++++++++++++++++++++++
!PART USED BY THE PRESSURE SOLVER
@@ -194,6 +196,7 @@ CHARACTER(LEN=5), POINTER :: CPRESOPT=>NULL()
INTEGER, POINTER :: NITR=>NULL()
LOGICAL, POINTER :: LITRADJ=>NULL()
LOGICAL, POINTER :: LRES=>NULL()
+LOGICAL, POINTER :: LOCEAN=>NULL()
REAL, POINTER :: XRES=>NULL()
REAL, POINTER :: XRELAX=>NULL()
REAL, POINTER :: XDXHATM=>NULL()
@@ -297,6 +300,7 @@ CPRESOPT=>DYN_MODEL(KTO)%CPRESOPT
NITR=>DYN_MODEL(KTO)%NITR
LITRADJ=>DYN_MODEL(KTO)%LITRADJ
LRES=>DYN_MODEL(KTO)%LRES
+LOCEAN=>DYN_MODEL(KTO)%LOCEAN
XRES=>DYN_MODEL(KTO)%XRES
XRELAX=>DYN_MODEL(KTO)%XRELAX
XDXHATM=>DYN_MODEL(KTO)%XDXHATM
diff --git a/src/MNH/modd_fieldn.f90 b/src/MNH/modd_fieldn.f90
index e5a933bbae17aa91b3be17c35e5ebb969dca0272..3b75c92ecbe1330d283d24e1579d24c121415ab7 100644
--- a/src/MNH/modd_fieldn.f90
+++ b/src/MNH/modd_fieldn.f90
@@ -147,6 +147,7 @@ REAL, DIMENSION(:,:,:), POINTER :: XSSPRO=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTKET=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XRTKES=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XPABST=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XPHIT=>NULL()
REAL, DIMENSION(:,:,:,:), POINTER :: XRT=>NULL()
REAL, DIMENSION(:,:,:,:), POINTER :: XRRS=>NULL()
REAL, DIMENSION(:,:,:,:), POINTER :: XRRS_CLD=>NULL()
@@ -252,6 +253,7 @@ XRTHS=>FIELD_MODEL(KTO)%XRTHS
!XTKET=>FIELD_MODEL(KTO)%XTKET !Done in FIELDLIST_GOTO_MODEL
XRTKES=>FIELD_MODEL(KTO)%XRTKES
!XPABST=>FIELD_MODEL(KTO)%XPABST !Done in FIELDLIST_GOTO_MODEL
+!XPHIT=>FIELD_MODEL(KTO)%XPHIT
!XRT=>FIELD_MODEL(KTO)%XRT !Done in FIELDLIST_GOTO_MODEL
XRRS=>FIELD_MODEL(KTO)%XRRS
!XRRS_CLD=>FIELD_MODEL(KTO)%XRRS_CLD !Done in FIELDLIST_GOTO_MODEL
diff --git a/src/MNH/modd_frc.f90 b/src/MNH/modd_frc.f90
index a430ad6034fce90152a49a63052db13cb2f87e2b..1e4b6c7ff5b91c4440abf45ce1f5b9e2bfe29281 100644
--- a/src/MNH/modd_frc.f90
+++ b/src/MNH/modd_frc.f90
@@ -45,6 +45,7 @@
!! add SST and surface pressure forcing
!! 01/2004 V. Masson surface externalization: removes SST forcing
!! 09/2017 Q.Rodier add LTEND_UV_FRC
+!! 03/2021 JL Redelsperger Parameters defining sfc forcing shape for idealized ocean case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -83,7 +84,7 @@ REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: XTENDVFRC ! large scale V tendency
LOGICAL, SAVE :: LGEOST_UV_FRC ! enables geostrophic wind term
LOGICAL, SAVE :: LGEOST_TH_FRC ! enables thermal wind advection
LOGICAL, SAVE :: LTEND_THRV_FRC ! enables tendency forcing
-LOGICAL, SAVE :: LTEND_UV_FRC ! enables tendency forcing of the wind
+LOGICAL, SAVE :: LTEND_UV_FRC ! enables tendency forcing of the wind
LOGICAL, SAVE :: LVERT_MOTION_FRC ! enables prescribed a forced vertical
! transport for all prognostic variables
LOGICAL, SAVE :: LRELAX_THRV_FRC ! enables temp. and humidity relaxation
@@ -100,5 +101,10 @@ LOGICAL, SAVE :: LTRANS ! enables a Galilean translation of the
! domain of simulation
LOGICAL, SAVE :: LPGROUND_FRC ! enables surf. pressure forcing
!
+LOGICAL, SAVE :: LDEEPOC ! activates sfc forcing for ideal ocean deep conv
+REAL, SAVE :: XCENTX_OC ! center of sfc forc for ideal ocean
+REAL, SAVE :: XRADX_OC ! radius of sfc forc for ideal ocean
+REAL, SAVE :: XCENTY_OC ! center of sfc forc for ideal ocean
+REAL, SAVE :: XRADY_OC ! radius of sfc forc for ideal ocean
!
END MODULE MODD_FRC
diff --git a/src/MNH/modd_oceanh.f90 b/src/MNH/modd_oceanh.f90
new file mode 100644
index 0000000000000000000000000000000000000000..f53e2063ebba8c04c39d0d32e51137c3234839d9
--- /dev/null
+++ b/src/MNH/modd_oceanh.f90
@@ -0,0 +1,47 @@
+!MNH_LIC Copyright 1996-2021 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.
+!-----------------------------------------------------------------
+! #################
+ MODULE MODD_OCEANH
+! #################
+!
+!!**** *MODD_OCEAN* - declaration of variables used in ocean version
+!!
+!! PURPOSE
+!! -------
+! Declarative module for the variables
+!! at interface for OCEAN LES MESONH version including auto-coupling O-A LES
+!
+!!
+!!** IMPLICIT ARGUMENTS
+!! ------------------
+!! None
+!!
+!! AUTHOR
+!! ------
+!! JL Redelsperger LOPS
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 03/2021
+!
+!* 0. DECLARATIONS
+! ------------
+ !
+USE MODD_TYPE_DATE
+!
+IMPLICIT NONE
+!
+!* fields for Sea Sfc FORCINGs
+! ------------------
+!
+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_T,XSSVFL_T,XSSTFL_T,XSSOLA_T ! given time forcing fluxes
+!
+END MODULE MODD_OCEANH
diff --git a/src/MNH/modd_ref.f90 b/src/MNH/modd_ref.f90
index 8b2f932878b87b893353170624ea8473cf167f79..206898425df34d55fb0f8a0f1120fe7f0b327132 100644
--- a/src/MNH/modd_ref.f90
+++ b/src/MNH/modd_ref.f90
@@ -46,6 +46,13 @@ REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XRHODREFZ ! rhod(z) for referenc
REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XTHVREFZ ! Thetav(z) for reference
! state without orography
REAL,SAVE :: XEXNTOP ! Exner function at model top
+!
+! For coupled A-O case
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XRHODREFZO! rhod(z) for ocean ref state in coupled mode
+REAL,SAVE, DIMENSION(:), ALLOCATABLE, TARGET :: XTHVREFZO !Thetav(z) for ocean ref state in coupled mode
+REAL,SAVE :: XEXNTOPO ! Exner function at ocean model top in coupled mode
+!
LOGICAL, SAVE :: LBOUSS ! Boussinesq approximation
+LOGICAL, SAVE ::LCOUPLES ! AUTOCOUPLED ATMS-OCEAN LES VERSION
!
END MODULE MODD_REF
diff --git a/src/MNH/modd_turbn.f90 b/src/MNH/modd_turbn.f90
index ac0b49ae30cbb2ee5af3243352036ed35f538986..5b938a6bdec96d14a3135984b12337506f22c753 100644
--- a/src/MNH/modd_turbn.f90
+++ b/src/MNH/modd_turbn.f90
@@ -40,6 +40,7 @@
!! C.Lac Nov 2014 add terms of TKE production for LES diag
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! D. Ricard May 2021 add the switches for Leonard terms
+!! JL Redelsperger 03/2021 Add O-A flux for auto-coupled LES case
!!
!-------------------------------------------------------------------------------
!
@@ -94,6 +95,10 @@ TYPE TURB_t
REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL() ! Transport production of Kinetic energy
REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL() ! Dissipation of Kinetic energy
REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL() ! Mixing length
+ REAL, DIMENSION(:,:,:), POINTER :: XSSUFL_C=>NULL() ! O-A interface flux for u
+ REAL, DIMENSION(:,:,:), POINTER :: XSSVFL_C=>NULL() ! O-A interface flux for v
+ REAL, DIMENSION(:,:,:), POINTER :: XSSTFL_C=>NULL() ! O-A interface flux for theta
+ REAL, DIMENSION(:,:,:), POINTER :: XSSRFL_C=>NULL() ! O-A interface flux for vapor
LOGICAL :: LHGRAD ! logical switch for the computation of the Leornard Terms
REAL :: XCOEFHGRADTHL ! coeff applied to thl contribution
REAL :: XCOEFHGRADRM ! coeff applied to mixing ratio contribution
@@ -133,6 +138,10 @@ REAL, DIMENSION(:,:,:), POINTER :: XTHP=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSUFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSVFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSTFL_C=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XSSRFL_C=>NULL()
LOGICAL, POINTER :: LHGRAD=>NULL()
REAL, POINTER :: XCOEFHGRADTHL=>NULL()
REAL, POINTER :: XCOEFHGRADRM=>NULL()
@@ -154,6 +163,10 @@ TURB_MODEL(KFROM)%XTHP=>XTHP
TURB_MODEL(KFROM)%XTR=>XTR
TURB_MODEL(KFROM)%XDISS=>XDISS
TURB_MODEL(KFROM)%XLEM=>XLEM
+TURB_MODEL(KFROM)%XSSUFL_C=>XSSUFL_C
+TURB_MODEL(KFROM)%XSSVFL_C=>XSSVFL_C
+TURB_MODEL(KFROM)%XSSTFL_C=>XSSTFL_C
+TURB_MODEL(KFROM)%XSSRFL_C=>XSSRFL_C
!
! Current model is set to model KTO
XIMPL=>TURB_MODEL(KTO)%XIMPL
@@ -183,6 +196,10 @@ XTHP=>TURB_MODEL(KTO)%XTHP
XTR=>TURB_MODEL(KTO)%XTR
XDISS=>TURB_MODEL(KTO)%XDISS
XLEM=>TURB_MODEL(KTO)%XLEM
+XSSUFL_C=>TURB_MODEL(KTO)%XSSUFL_C
+XSSVFL_C=>TURB_MODEL(KTO)%XSSVFL_C
+XSSTFL_C=>TURB_MODEL(KTO)%XSSTFL_C
+XSSRFL_C=>TURB_MODEL(KTO)%XSSRFL_C
LHGRAD=>TURB_MODEL(KTO)%LHGRAD
XCOEFHGRADTHL=>TURB_MODEL(KTO)%XCOEFHGRADTHL
XCOEFHGRADRM=>TURB_MODEL(KTO)%XCOEFHGRADRM
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 0c5a12fef4a4452c662ba91efc49b836ede28e44..91b86154ebc81cc6916f2c245bb57367683a7ae0 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -271,6 +271,7 @@ END MODULE MODI_MODEL_n
! F. Auguste 01/02/2021: add IBM
! T. Nagel 01/02/2021: add turbulence recycling
! P. Wautelet 19/02/2021 add NEGA2 term for SV budgets
+! J.L. Redelsperger 03/2021, add Call NHOA_COUPLN (coupling O & A LES version)
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -349,6 +350,7 @@ USE MODD_PROFILER_n
USE MODD_RADIATIONS_n, ONLY: XTSRAD,XSCAFLASWD,XDIRFLASWD,XDIRSRFSWD, XAER, XDTHRAD
USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN
USE MODD_RECYCL_PARAM_n
+USE MODD_REF
USE MODD_REF_n
USE MODD_SALT, ONLY: LSALT
USE MODD_SERIES, ONLY: LSERIES
@@ -756,11 +758,14 @@ CALL SECOND_MNH2(ZTIME1)
!
ISYNCHRO = MODULO (KTCOUNT, NDTRATIO(IMI) ) ! test of synchronisation
!
-
-
-IF (IMI/=1 .AND. NDAD(IMI)/=IMI .AND. (ISYNCHRO==1 .OR. NDTRATIO(IMI) == 1) ) THEN
-!
- ! Use dummy pointers to correct an ifort BUG
+!
+IF (LCOUPLES.AND.LOCEAN) THEN
+ CALL NHOA_COUPL_n(NDAD(IMI),XTSTEP,IMI,KTCOUNT,IKU)
+END IF
+! No Gridnest in coupled OA LES for now
+IF (.NOT. LCOUPLES .AND. IMI/=1 .AND. NDAD(IMI)/=IMI .AND. (ISYNCHRO==1 .OR. NDTRATIO(IMI) == 1) ) THEN
+!
+! Use dummy pointers to correct an ifort BUG
DPTR_XBMX1=>XBMX1
DPTR_XBMX2=>XBMX2
DPTR_XBMX3=>XBMX3
diff --git a/src/MNH/modn_dynn.f90 b/src/MNH/modn_dynn.f90
index 0f6585f99f5adf5dd4cd7466946cd2b1cd1a2711..43d1f326be9763663dde2b5a225918b87fb757e3 100644
--- a/src/MNH/modn_dynn.f90
+++ b/src/MNH/modn_dynn.f90
@@ -55,6 +55,7 @@
!! Modifications 09/06/11 (Barthe) add LHORELAX_SVELEC in namelist
!! Modifications 15/06/11 (Lac) add LHORELAX for conditional sampling
!! Modifications 12/02/12 (Pialat/Tulet) add LHORELAX_SVFF for ForeFire scalar variables
+!! Modification 03/2021 (JL Redelsperger) add logical LOCEAN for ocean LES version
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -62,6 +63,7 @@
!
USE MODD_PARAMETERS, ONLY : JPSVMAX
USE MODD_DYN_n, ONLY : &
+ LOCEAN_n => LOCEAN, &
XTSTEP_n => XTSTEP, &
CPRESOPT_n => CPRESOPT, &
NITR_n => NITR, &
@@ -109,7 +111,8 @@ REAL ,SAVE :: XTSTEP
CHARACTER(LEN=5),SAVE :: CPRESOPT
INTEGER ,SAVE :: NITR
LOGICAL ,SAVE :: LITRADJ
-LOGICAL ,SAVE :: LRES
+LOGICAL ,SAVE :: LRES
+LOGICAL ,SAVE :: LOCEAN
REAL ,SAVE :: XRES
REAL ,SAVE :: XRELAX
LOGICAL, SAVE :: LHORELAX_UVWTH
@@ -156,12 +159,14 @@ NAMELIST/NAM_DYNn/XTSTEP,CPRESOPT,NITR,LITRADJ,LRES,XRES,XRELAX,LHORELAX_UVWTH,
#ifdef MNH_FOREFIRE
LHORELAX_SVFF, &
#endif
+ LOCEAN,&
NRIMX,NRIMY,XRIMKMAX,XT4DIFU, &
XT4DIFTH,XT4DIFSV
!
CONTAINS
!
SUBROUTINE INIT_NAM_DYNn
+ LOCEAN = LOCEAN_n
XTSTEP = XTSTEP_n
CPRESOPT = CPRESOPT_n
NITR = NITR_n
@@ -205,6 +210,7 @@ SUBROUTINE INIT_NAM_DYNn
END SUBROUTINE INIT_NAM_DYNn
SUBROUTINE UPDATE_NAM_DYNn
+ LOCEAN_n = LOCEAN
XTSTEP_n = XTSTEP
CPRESOPT_n = CPRESOPT
NITR_n = NITR
diff --git a/src/MNH/modn_frc.f90 b/src/MNH/modn_frc.f90
index 685ee4f243d5c526fed26540b16b1cabf6259430..50577a0e77e5344617a74c9923a86779f4ad434b 100644
--- a/src/MNH/modn_frc.f90
+++ b/src/MNH/modn_frc.f90
@@ -66,6 +66,11 @@ NAMELIST /NAM_FRC/ LGEOST_UV_FRC , &
LTRANS , &
XUTRANS , &
XVTRANS , &
- LPGROUND_FRC
+ LPGROUND_FRC ,&
+ LDEEPOC ,&
+ XCENTX_OC ,&
+ XCENTY_OC ,&
+ XRADX_OC ,&
+ XRADY_OC
!
END MODULE MODN_FRC
diff --git a/src/MNH/modn_nesting.f90 b/src/MNH/modn_nesting.f90
index 90a3f00212fc22a1a5b535a9353b1bb736c5571e..f7c12be2aa19db9db935d420c161d3746d1f0174 100644
--- a/src/MNH/modn_nesting.f90
+++ b/src/MNH/modn_nesting.f90
@@ -43,16 +43,18 @@
!!
!! MODIFICATIONS
!! -------------
-!! Original 16/08/95
+!! Original 16/08/95
+!! JL Redelsperger 03/2021 : Add Auto-coupled O-A LES case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_NESTING
+USE MODD_REF, ONLY: LCOUPLES
!
IMPLICIT NONE
!
-NAMELIST/NAM_NESTING/NDAD,NDTRATIO,XWAY
+NAMELIST/NAM_NESTING/NDAD,NDTRATIO,XWAY,LCOUPLES
!
END MODULE MODN_NESTING
diff --git a/src/MNH/nhoa_coupln.f90 b/src/MNH/nhoa_coupln.f90
new file mode 100644
index 0000000000000000000000000000000000000000..78810dea127ab5c4a7752bbbfcaf6ebecb719931
--- /dev/null
+++ b/src/MNH/nhoa_coupln.f90
@@ -0,0 +1,155 @@
+!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.
+!-----------------------------------------------------------------
+! ###################
+ MODULE MODI_NHOA_COUPL_n
+! ###################
+!
+INTERFACE
+!
+ SUBROUTINE NHOA_COUPL_n(KDAD,PTSTEP,KMI,KTCOUNT,KKU)
+!
+INTEGER, INTENT(IN) :: KDAD ! Number of the DAD model
+REAL, INTENT(IN) :: PTSTEP ! Time step
+INTEGER, INTENT(IN) :: KMI ! model number
+INTEGER, INTENT(IN) :: KKU
+INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop COUNTer
+ ! (=1 at the segment beginning)
+!
+END SUBROUTINE NHOA_COUPL_n
+END INTERFACE
+END MODULE MODI_NHOA_COUPL_n
+!
+! ####################################################################
+SUBROUTINE NHOA_COUPL_n(KDAD,PTSTEP,KMI,KTCOUNT,KKU)
+! ####################################################################
+!!
+!! PURPOSE
+!! -------
+! To compute the flux at the O-A interface in the auto-coupling O-A LES case
+!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: JPHEXT,JPVEXT
+!! Module MODD_FIELD$n : XUT,XVT,XWT,XRT,XTHT,XPABST
+!! Module MODD_REF$n : XRHODJ, XRVREF,XTHVREF, XRHODREF
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! JL Redelsperger 03/2021 Version 0
+!! MODIFICATIONS
+!! -------------
+!!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODE_ll
+USE MODE_MODELN_HANDLER
+!
+USE MODD_PARAMETERS
+USE MODD_NESTING
+USE MODD_CST
+USE MODD_REF_n ! modules relative to the outer model $n
+USE MODD_FIELD_n
+USE MODD_CONF
+USE MODD_PARAM_n
+USE MODD_TURB_n
+USE MODD_DYN_n, ONLY : LOCEAN
+USE MODD_REF, ONLY: LCOUPLES
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+INTEGER, INTENT(IN) :: KDAD ! Number of the DAD model
+REAL, INTENT(IN) :: PTSTEP ! Time step
+INTEGER, INTENT(IN) :: KMI ! model number
+INTEGER, INTENT(IN) :: KKU !
+!
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: IIB,IIE,IJB,IJE,IIU,IJU
+INTEGER :: IKE
+INTEGER :: ILBX,ILBY,ILBX2,ILBY2
+INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop COUNTer
+ ! (=1 at the segment beginning)
+!
+INTEGER :: JRR,JSV ! Loop index
+!
+INTEGER :: IINFO_ll, IDIMX, IDIMY
+! surface variables: wind, current, Temp
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPUA,ZCOUPVA,ZCOUPTA
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPUO,ZCOUPVO,ZCOUPTO
+!surf flux local work space
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOUPTFL,ZCOUPUFL,ZCOUPVFL
+CHARACTER(LEN=4) :: ZINIT_TYPE
+!
+!---Coupled OA MesoNH----------------------------------------------------------------------------
+!* 0. INITIALISATION
+! --------------
+! allocate flux local array
+ALLOCATE(ZCOUPTFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPUFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPVFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+! allocate sfc variable local array
+ALLOCATE(ZCOUPUA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPVA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPTA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPUO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPVO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPTO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+! values in ocean sfc
+IKE=KKU-JPVEXT
+ZCOUPUO(:,:)= XUT(:,:,IKE)
+ZCOUPVO(:,:)= XVT(:,:,IKE)
+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
+! rho_atmos * (w'u')_atmos = rho_ocean * (u'w')_ocean
+! rho_atmos *Cp_atmos* (u'w')_atmos = rho_ocean * CP_ocean * (u'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)
+!
+END SUBROUTINE NHOA_COUPL_n
+
diff --git a/src/MNH/one_wayn.f90 b/src/MNH/one_wayn.f90
index dae09f3c1fb37f07aab0a4caaa7996cf63b98b41..07ce225e37fe62abba178b9a61959628efac9a87 100644
--- a/src/MNH/one_wayn.f90
+++ b/src/MNH/one_wayn.f90
@@ -128,6 +128,7 @@ SUBROUTINE ONE_WAY_n(KDAD,PTSTEP,KMI,KTCOUNT, &
USE MODD_CH_MNHC_n, only: LUSECHAQ, LUSECHIC
USE MODD_CONF, only: CEQNSYS
USE MODD_CST, only: XCPD, XP00, XRD, XRV, XTH00
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_FIELD_n, only: XPABST, XRT, XSVT, XUT, XVT, XWT, XTHT, XTKET
USE MODD_NESTING, only: NXOR_ALL, NXEND_ALL, NYOR_ALL, NYEND_ALL
USE MODD_NSV, only: NSV_A, NSV_C1R3BEG_A, NSV_C1R3_A, NSV_C2R2BEG_A, NSV_C2R2_A, NSV_CHEMBEG_A, NSV_CHEMEND_A, &
@@ -139,6 +140,7 @@ USE MODD_NSV, only: NSV_A, NSV_C1R3BEG_A, NSV_C1R3_A, NSV_C2R2BEG_A,
USE MODD_PARAMETERS, only: JPHEXT, JPVEXT
USE MODD_PARAM_n, only: CCLOUD
+USE MODD_REF, ONLY : LCOUPLES
USE MODD_REF_n, only: XRHODJ, XRHODREF, XRVREF, XTHVREF
!
use mode_bikhardt
@@ -239,6 +241,21 @@ REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCHEMTI
!
integer :: igrid
!
+IF (LCOUPLES) THEN
+ PDRYMASST=0.
+ PDRYMASSS=0.
+ PLBXUS=0.
+ PLBXVS=0.
+ PLBXWS=0.
+ PLBXTHS=0.
+ PLBYTHS=0.
+ PLBXTKES=0.
+ PLBYTKES =0.
+ PLBXRS =0.
+ PLBYRS=0.
+ PLBXSVS =0.
+ PLBYSVS=0.
+ELSE
!-------------------------------------------------------------------------------
!
!* 0. INITIALISATION
@@ -790,6 +807,8 @@ DEALLOCATE(ZWORK)
DEALLOCATE(ZCOEFLIN_LBXM_RED,ZCOEFLIN_LBYM_RED,IKLIN_LBXM_RED,IKLIN_LBYM_RED)
!
!------------------------------------------------------------------------------
+ENDIF ! END LCOUPLES couplage
+!
CALL GOTO_MODEL(KMI)
!
END SUBROUTINE ONE_WAY_n
diff --git a/src/MNH/p_abs.f90 b/src/MNH/p_abs.f90
index 398caa3b9036055871a7e4afbd71f214097f3e95..c89380368cc5d84eb2e496438906c323c576648a 100644
--- a/src/MNH/p_abs.f90
+++ b/src/MNH/p_abs.f90
@@ -16,7 +16,7 @@ INTERFACE
!
SUBROUTINE P_ABS (KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
PTHT, PRT, PRHODJ, PRHODREF, PTHETAV, PTHVREF, &
- PRVREF, PEXNREF, PPHIT )
+ PRVREF, PEXNREF, PPHIT, PPHI0)
!
IMPLICIT NONE
!
@@ -44,7 +44,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF! Exner function of the
!
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHIT ! Perturbation of
! either the Exner function Pi or Pi * Cpd * THvref
-!
+REAL, INTENT(INOUT) :: PPHI0 ! Phi0 at time t !
!
END SUBROUTINE P_ABS
!
@@ -54,7 +54,7 @@ END MODULE MODI_P_ABS
! #######################################################################
SUBROUTINE P_ABS (KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
PTHT, PRT, PRHODJ, PRHODREF, PTHETAV, PTHVREF, &
- PRVREF, PEXNREF, PPHIT )
+ PRVREF, PEXNREF, PPHIT, PPHI0 )
! #######################################################################
!
!!**** *P_ABS * - routine to compute the absolute Exner pressure deviation PHI
@@ -108,6 +108,8 @@ END MODULE MODI_P_ABS
!! from Durran (1989), MAE and DUR respectively
!! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation
!! J. Colin 07/13 Add LBOUSS
+!! J.L Redelsperger 03/2021 Change of one step to pressure computation
+!! in order to perform Ocean runs (equivalent to LHE shallow convection)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -116,6 +118,7 @@ END MODULE MODI_P_ABS
USE MODD_CST
USE MODD_CONF
USE MODD_PARAMETERS
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_REF, ONLY : LBOUSS
!
USE MODE_ll
@@ -151,6 +154,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVREF ! vapor mixing ratio
REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF! Exner function of the
! reference state
!
+REAL, INTENT(INOUT) :: PPHI0 ! PHI0 at t
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHIT ! Perturbation of
! either the Exner function Pi or Pi * Cpd * THvref
!
@@ -357,7 +361,7 @@ ELSEIF( CEQNSYS == 'LHE' ) THEN
ENDWHERE
ENDIF
!
- ! compute the absolute pressure function
+ ! compute the absolute pressure function (LHE equation system case)
!
!
!
@@ -379,9 +383,18 @@ ELSEIF( CEQNSYS == 'LHE' ) THEN
ZMASSGUESS = SUM_DD_R2_ll(ZMASSGUESS_2D)
ZWATERMASST = SUM_DD_R2_ll(ZWATERMASST_2D)
!
- ZPHI0 = (PDRYMASST + ZWATERMASST - 2. * PREFMASS + ZMASSGUESS ) / PMASS_O_PHI0
- PPHIT(:,:,:) = PPHIT(:,:,:) + ZPHI0
- !
+!
+! case shallow bouss : to get the real pressure fluctuation
+! Eq 2.40 p15 : constant not resolved in poisson equation
+IF (.NOT. LOCEAN) THEN
+ PPHI0 = (PDRYMASST + ZWATERMASST - 2. * PREFMASS + ZMASSGUESS ) / PMASS_O_PHI0
+ELSE
+! PPHI0 = 0. => to be possibly modified for ocean LES case
+ PPHI0=0.
+END IF
+! following computation moved in PRESSURE routine (Eq 2.40 bis p15: Phi_total)
+! PPHIT(:,:,:) = PPHIT(:,:,:) + ZPHI0
+!
END IF
!
!-------------------------------------------------------------------------------
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index 5b0c9f50b45c371938296c04b5339bb9349bfc65..d99996d8bf042d61000603fc781d942bd9935b3a 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -236,6 +236,7 @@ END MODULE MODI_PHYS_PARAM_n
! P. Wautelet 21/11/2019: ZRG_HOUR and ZRAT_HOUR are now parameter arrays
!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!! 02/2021 F.Auguste: add IBM
+!! 03/2021: JL Redelsperger Add the SW flux penetration for Ocean model case
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -283,6 +284,7 @@ USE MODD_METRICS_n
USE MODD_MNH_SURFEX_n
USE MODD_NESTING, ONLY : XWAY,NDAD, NDXRATIO_ALL, NDYRATIO_ALL
USE MODD_NSV
+USE MODD_OCEANH
USE MODD_OUT_n
USE MODD_PARAM_C2R2, ONLY : LSEDC
USE MODD_PARAMETERS
@@ -298,6 +300,7 @@ USE MODD_PRECIP_n
use modd_precision, only: MNHTIME
USE MODD_RADIATIONS_n
USE MODD_RAIN_ICE_DESCR, ONLY : XRTMIN
+USE MODD_REF, ONLY : LCOUPLES
USE MODD_REF_n
USE MODD_SALT
USE MODD_SHADOWS_n
@@ -444,6 +447,11 @@ INTEGER :: IKIDM ! index loop
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSAVE_INPRR,ZSAVE_INPRS,ZSAVE_INPRG,ZSAVE_INPRH
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSAVE_INPRC,ZSAVE_PRCONV,ZSAVE_PRSCONV
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZSAVE_DIRFLASWD, ZSAVE_SCAFLASWD,ZSAVE_DIRSRFSWD
+! for ocean model
+INTEGER :: JKM , JSW ! vertical index loop
+REAL :: ZSWA,TINTSW ! index for SW interpolation and int time betwenn forcings (ocean model)
+REAL, DIMENSION(:), ALLOCATABLE :: ZIZOCE(:) ! Solar flux penetrating in ocean
+REAL, DIMENSION(:), ALLOCATABLE :: ZPROSOL1(:),ZPROSOL2(:) ! Funtions for penetrating solar flux
!
!-----------------------------------------------------------------------------
@@ -806,6 +814,48 @@ IF (CRAD /='NONE') THEN
if ( lbudget_th ) call Budget_store_end ( tbudgets(NBUDGET_TH), 'RAD', xrths(:, :, :) )
END IF
!
+!
+!* 1.6 Ocean case:
+! Sfc turbulent fluxes & Radiative tendency due to SW penetrating ocean
+!
+IF (LOCEAN .AND. (.NOT.LCOUPLES)) THEN
+!
+ ALLOCATE( ZIZOCE(IKU)); ZIZOCE(:)=0.
+ ALLOCATE( ZPROSOL1(IKU))
+ ALLOCATE( ZPROSOL2(IKU))
+ ALLOCATE(XSSUFL(IIU,IJU))
+ ALLOCATE(XSSVFL(IIU,IJU))
+ ALLOCATE(XSSTFL(IIU,IJU))
+ ALLOCATE(XSSOLA(IIU,IJU))
+! Time interpolation
+ JSW = INT(TDTCUR%xtime/REAL(NINFRT))
+ ZSWA = TDTCUR%xtime/REAL(NINFRT)-REAL(JSW)
+ XSSTFL = (XSSTFL_T(JSW+1)*(1.-ZSWA)+XSSTFL_T(JSW+2)*ZSWA)
+ XSSUFL = (XSSUFL_T(JSW+1)*(1.-ZSWA)+XSSUFL_T(JSW+2)*ZSWA)
+ XSSVFL = (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) = XROC*ZIZOCE(IKU)
+ ZPROSOL2(IKU) = (1.-XROC)*ZIZOCE(IKU)
+ IF(NVERB >= 5 ) THEN
+ WRITE(ILUOUT,*)'ZSWA JSW TDTCUR XTSTEP FT FU FV SolarR(IKU)', NINFRT, ZSWA,JSW,&
+ TDTCUR%xtime, XTSTEP, XSSTFL(2,2), XSSUFL(2,2),XSSVFL(2,2),ZIZOCE(IKU)
+ END IF
+ DO JKM=IKU-1,2,-1
+ ZPROSOL1(JKM) = ZPROSOL1(JKM+1)* exp(-XDZZ(2,2,JKM)/XD1)
+ ZPROSOL2(JKM) = ZPROSOL2(JKM+1)* exp(-XDZZ(2,2,JKM)/XD2)
+ ZIZOCE(JKM) = (ZPROSOL1(JKM+1)-ZPROSOL1(JKM) + ZPROSOL2(JKM+1)-ZPROSOL2(JKM))/XDZZ(2,2,JKM)
+! Adding to tep temp tendency, the solar radiation penetrating in ocean
+ if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'OCE', xrths(:, :, :) )
+ XRTHS(:,:,JKM) = XRTHS(:,:,JKM) + XRHODJ(:,:,JKM)*ZIZOCE(JKM)
+ if ( lbudget_th ) call Budget_store_end ( tbudgets(NBUDGET_TH), 'OCE', xrths(:, :, :) )
+ END DO
+ DEALLOCATE( ZIZOCE)
+ DEALLOCATE (ZPROSOL1)
+ DEALLOCATE (ZPROSOL2)
+END IF
+!
+!
CALL SECOND_MNH2(ZTIME2)
!
PRAD = PRAD + ZTIME2 - ZTIME1 &
@@ -1526,7 +1576,14 @@ PMAFL = PMAFL + ZTIME4 - ZTIME3 - ZTIME_LES_MF
PTIME_BU = PTIME_BU + XTIME_LES_BU_PROCESS + XTIME_BU_PROCESS
!
!
+!* deallocate sf flux array for ocean model (in grid nesting, dimensions can vary)
!
+IF (LOCEAN .AND. (.NOT. LCOUPLES)) THEN
+ DEALLOCATE(XSSUFL)
+ DEALLOCATE(XSSVFL)
+ DEALLOCATE(XSSTFL)
+ DEALLOCATE(XSSOLA)
+END IF
!-------------------------------------------------------------------------------
!
!* deallocation of variables used in more than one parameterization
diff --git a/src/MNH/prandtl.f90 b/src/MNH/prandtl.f90
index dc9e578612d8bac8a73a236086e6af7bc890c9a1..fa9aade6145431a2df5e5d5f8e366481b9c4575a 100644
--- a/src/MNH/prandtl.f90
+++ b/src/MNH/prandtl.f90
@@ -187,6 +187,7 @@ END MODULE MODI_PRANDTL
!! vertical levels
!! 2017-09 J.Escobar, use epsilon XMNH_TINY_12 for R*4
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! JL Redelsperger 03/2021 : adding Ocean case for temperature only
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -195,6 +196,7 @@ END MODULE MODI_PRANDTL
USE MODD_CST
USE MODD_CONF
USE MODD_CTURB
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_field, only: tfielddata, TYPEREAL
USE MODD_IO, ONLY: TFILEDATA
USE MODD_PARAMETERS
@@ -287,15 +289,22 @@ PEMOIST(:,:,KKA) = 2.*PEMOIST(:,:,IKB) - PEMOIST(:,:,IKB+KKL)
!
! 1.3 1D Redelsperger numbers
!
-PBLL_O_E(:,:,:) = MZM( XG / PTHVREF(:,:,:) * PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:) )
-IF (KRR /= 0) THEN ! moist case
- PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * PETHETA(:,:,:) * &
- & GZ_M_W(KKA,KKU,KKL,PTHLM,PDZZ)
- PREDR1(:,:,:) = XCTV*PBLL_O_E(:,:,:) * PEMOIST(:,:,:) * &
- & GZ_M_W(KKA,KKU,KKL,PRM(:,:,:,1),PDZZ)
-ELSE ! dry case
+IF (LOCEAN) THEN
+ PBLL_O_E(:,:,:) = MZM(XG *XALPHAOC* PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:) )
PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * GZ_M_W(KKA,KKU,KKL,PTHLM,PDZZ)
PREDR1(:,:,:) = 0.
+ELSE
+ PBLL_O_E(:,:,:) = MZM(XG / PTHVREF(:,:,:) * PLM(:,:,:) * PLEPS(:,:,:) / PTKEM(:,:,:) )
+ IF (KRR /= 0) THEN ! moist case
+ PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * PETHETA(:,:,:) * &
+ & GZ_M_W(KKA,KKU,KKL,PTHLM,PDZZ)
+ PREDR1(:,:,:) = XCTV*PBLL_O_E(:,:,:) * PEMOIST(:,:,:) * &
+ & GZ_M_W(KKA,KKU,KKL,PRM(:,:,:,1),PDZZ)
+ ELSE ! dry case
+ PREDTH1(:,:,:)= XCTV*PBLL_O_E(:,:,:) * GZ_M_W(KKA,KKU,KKL,PTHLM,PDZZ)
+ PREDR1(:,:,:) = 0.
+ END IF
+!
END IF
!
! 3. Limits on 1D Redelperger numbers
@@ -335,9 +344,11 @@ PREDR1 (:,:,:) = PREDR1 (:,:,:) * ZW1(:,:,:)
ZW2=SIGN(1.,PREDTH1(:,:,:))
PREDTH1(:,:,:)= ZW2(:,:,:) * MAX(XMNH_TINY_12, ZW2(:,:,:)*PREDTH1(:,:,:))
!
-IF (KRR /= 0) THEN ! dry case
- ZW2=SIGN(1.,PREDR1(:,:,:))
- PREDR1(:,:,:)= ZW2(:,:,:) * MAX(XMNH_TINY_12, ZW2(:,:,:)*PREDR1(:,:,:))
+IF (.NOT.LOCEAN) THEN
+ IF (KRR /= 0) THEN ! dry case
+ ZW2=SIGN(1.,PREDR1(:,:,:))
+ PREDR1(:,:,:)= ZW2(:,:,:) * MAX(XMNH_TINY_12, ZW2(:,:,:)*PREDR1(:,:,:))
+ END IF
END IF
!
!
@@ -448,57 +459,73 @@ IF(HTURBDIM=='1DIM') THEN
!
ELSE IF (L2D) THEN ! 3D case in a 2D model
!
- DO JSV=1,ISV
+ IF (LOCEAN) THEN
IF (KRR /= 0) THEN
- ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
+ ZW1 = MZM((XG *XALPHAOC * PLM * PLEPS / PTKEM)**2 ) *PETHETA
ELSE
- ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2)
- END IF
- PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV) + &
- ZW1* &
- MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
- GX_M_M(PTHLM,PDXX,PDZZ,PDZX) &
- )
-!
- IF (KRR /= 0) THEN
- PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV) + &
- ZW1 * PEMOIST * &
- MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
- GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
- )
- ELSE
- PRED2RS3(:,:,:,JSV) = 0.
- END IF
- ENDDO
+ ZW1 = MZM((XG *XALPHAOC * PLM * PLEPS / PTKEM)**2)
+ END IF
+ ELSE
+ DO JSV=1,ISV
+ IF (KRR /= 0) THEN
+ ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
+ ELSE
+ ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2)
+ END IF
+ PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV) + &
+ ZW1* &
+ MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
+ GX_M_M(PTHLM,PDXX,PDZZ,PDZX) &
+ )
+!
+ IF (KRR /= 0) THEN
+ PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV) + &
+ ZW1 * PEMOIST * &
+ MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
+ GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
+ )
+ ELSE
+ PRED2RS3(:,:,:,JSV) = 0.
+ END IF
+ ENDDO
+ END IF
!
ELSE ! 3D case in a 3D model
!
- DO JSV=1,ISV
- IF (KRR /= 0) THEN
- ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
- ELSE
- ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2)
- END IF
- PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV) + &
- ZW1* &
- MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
- GX_M_M(PTHLM,PDXX,PDZZ,PDZX) &
- +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)* &
- GY_M_M(PTHLM,PDYY,PDZZ,PDZY) &
- )
-!
- IF (KRR /= 0) THEN
- PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV) + &
- ZW1 * PEMOIST * &
- MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
- GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
- +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)* &
- GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY) &
- )
- ELSE
- PRED2RS3(:,:,:,JSV) = 0.
- END IF
- ENDDO
+ IF (LOCEAN) THEN
+ IF (KRR /= 0) THEN
+ ZW1 = MZM((XG *XALPHAOC * PLM * PLEPS / PTKEM)**2 ) *PETHETA
+ ELSE
+ ZW1 = MZM((XG *XALPHAOC * PLM * PLEPS / PTKEM)**2)
+ END IF
+ ELSE
+ DO JSV=1,ISV
+ IF (KRR /= 0) THEN
+ ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2 ) *PETHETA
+ ELSE
+ ZW1 = MZM( (XG / PTHVREF * PLM * PLEPS / PTKEM)**2)
+ END IF
+ PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV) + &
+ ZW1* &
+ MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
+ GX_M_M(PTHLM,PDXX,PDZZ,PDZX) &
+ +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)* &
+ GY_M_M(PTHLM,PDYY,PDZZ,PDZY) &
+ )
+!
+ IF (KRR /= 0) THEN
+ PRED2RS3(:,:,:,JSV) = PREDR1(:,:,:) * PREDS1(:,:,:,JSV) + &
+ ZW1 * PEMOIST * &
+ MZM(GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)* &
+ GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
+ +GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)* &
+ GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY) &
+ )
+ ELSE
+ PRED2RS3(:,:,:,JSV) = 0.
+ END IF
+ ENDDO
+ END IF
!
END IF ! end of HTURBDIM if-block
!
diff --git a/src/MNH/prep_ideal_case.f90 b/src/MNH/prep_ideal_case.f90
index 665d038d648d00b933f4ac38adfcd719c352c5b7..ffae04e797d279e43e196855fdbd4ea54cb72e46 100644
--- a/src/MNH/prep_ideal_case.f90
+++ b/src/MNH/prep_ideal_case.f90
@@ -320,6 +320,7 @@
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
! F. Auguste 02/2021 : add IBM
! P. Wautelet 09/03/2021: move some chemistry initializations to ini_nsv
+! Jean-Luc Redelsperger 03/2021 : : ocean LES case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -576,7 +577,7 @@ TYPE(TFILEDATA),POINTER :: TZEXPREFILE => NULL()
NAMELIST/NAM_CONF_PRE/ LTHINSHELL,LCARTESIAN, &! Declarations in MODD_CONF
LPACK, &!
NVERB,CIDEAL,CZS, &!+global variables initialized
- LBOUSS,LPERTURB, &! at their declarations
+ LBOUSS,LOCEAN,LPERTURB, &! at their declarations
LFORCING,CEQNSYS, &! at their declarations
LSHIFT,L2D_ADV_FRC,L2D_REL_FRC, &
NHALO , JPHEXT
@@ -976,6 +977,11 @@ ALLOCATE(XDZZ(NIU,NJU,NKU))
!
ALLOCATE(XRHODREFZ(NKU),XTHVREFZ(NKU))
XTHVREFZ(:)=0.0
+IF (LCOUPLES) THEN
+ ! Arrays for reference state different in ocean and atmosphere
+ ALLOCATE(XRHODREFZO(NKU),XTHVREFZO(NKU))
+ XTHVREFZO(:)=0.0
+END IF
IF(CEQNSYS == 'DUR') THEN
ALLOCATE(XRVREF(NIU,NJU,NKU))
ELSE
@@ -990,7 +996,11 @@ ALLOCATE(XLSUM(NIU,NJU,NKU))
ALLOCATE(XLSVM(NIU,NJU,NKU))
ALLOCATE(XLSWM(NIU,NJU,NKU))
ALLOCATE(XLSTHM(NIU,NJU,NKU))
-ALLOCATE(XLSRVM(NIU,NJU,NKU))
+IF ( NRR >= 1) THEN
+ ALLOCATE(XLSRVM(NIU,NJU,NKU))
+ELSE
+ ALLOCATE(XLSRVM(0,0,0))
+ENDIF
!
! allocate lateral boundary field used for coupling
!
@@ -1635,9 +1645,11 @@ IF(LPERTURB) CALL SET_PERTURB(TZEXPREFILE)
!
!* 5.9 Anelastic correction and pressure:
!
-CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
-IF ( .NOT. L1D ) CALL PRESSURE_IN_PREP(XDXX,XDYY,XDZX,XDZY,XDZZ)
-CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
+IF (.NOT.LOCEAN) THEN
+ CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
+ IF ( .NOT. L1D ) CALL PRESSURE_IN_PREP(XDXX,XDYY,XDZX,XDZY,XDZZ)
+ CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
+END IF
!
!
!* 5.10 Compute THETA, vapor and cloud mixing ratio
@@ -1654,19 +1666,28 @@ IF (CIDEAL == 'RSOU') THEN
ALLOCATE(ZRSATW(NIU,NJU,NKU))
ALLOCATE(ZRSATI(NIU,NJU,NKU))
ZRT=XRT(:,:,:,1)+XRT(:,:,:,2)+XRT(:,:,:,4)
+IF (LOCEAN) THEN
+ ZEXN(:,:,:)= 1.
+ ZT=XTHT
+ ZTHL=XTHT
+ ZCPH=XCPD+ XCPV * XRT(:,:,:,1)
+ ZLVOCPEXN = XLVTT
+ ZLSOCPEXN = XLSTT
+ELSE
ZEXN=(XPABST/XP00) ** (XRD/XCPD)
ZT=XTHT*(XPABST/XP00)**(XRD/XCPD)
ZCPH=XCPD+ XCPV * XRT(:,:,:,1)+ XCL *XRT(:,:,:,2) + XCI * XRT(:,:,:,4)
ZLVOCPEXN = (XLVTT + (XCPV-XCL) * (ZT-XTT))/(ZCPH*ZEXN)
ZLSOCPEXN = (XLSTT + (XCPV-XCI) * (ZT-XTT))/(ZCPH*ZEXN)
ZTHL=XTHT-ZLVOCPEXN*XRT(:,:,:,2)-ZLSOCPEXN*XRT(:,:,:,4)
+ CALL TH_R_FROM_THL_RT_3D('T',ZFRAC_ICE,XPABST,ZTHL,ZRT,XTHT,XRT(:,:,:,1), &
+ XRT(:,:,:,2),XRT(:,:,:,4),ZRSATW, ZRSATI)
+END IF
DEALLOCATE(ZEXN)
DEALLOCATE(ZT)
DEALLOCATE(ZCPH)
DEALLOCATE(ZLVOCPEXN)
DEALLOCATE(ZLSOCPEXN)
- CALL TH_R_FROM_THL_RT_3D('T',ZFRAC_ICE,XPABST,ZTHL,ZRT,XTHT,XRT(:,:,:,1), &
- XRT(:,:,:,2),XRT(:,:,:,4),ZRSATW, ZRSATI)
DEALLOCATE(ZTHL)
DEALLOCATE(ZRT)
! Coherence test
diff --git a/src/MNH/pressurez.f90 b/src/MNH/pressurez.f90
index 7cf850831db046798050853a1a2420b6bc15e302..23c74074db140ad39a6145cf95b1a7d7657ec246 100644
--- a/src/MNH/pressurez.f90
+++ b/src/MNH/pressurez.f90
@@ -220,6 +220,8 @@ END MODULE MODI_PRESSUREZ
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
! P. Wautelet 28/01/2020: use the new data structures and subroutines for budgets for U
+!! JL Redelsperger 03/2021 : Shallow convection case added in LHE case:
+!! working for both atmosphere and ocean cases
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -229,7 +231,8 @@ USE MODD_ARGSLIST_ll, ONLY: LIST_ll
use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets
USE MODD_CST
USE MODD_CONF
-USE MODD_DYN_n, ONLY: LRES, XRES
+USE MODD_DYN_n, ONLY: LRES, XRES,LOCEAN
+USE MODD_FIELD_n, ONLY: XPHIT
USE MODD_IBM_PARAM_n, ONLY : XIBM_LS,XIBM_SU,LIBM,NIBM_ITR,XIBM_EPSI
USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_MPIF
@@ -359,6 +362,7 @@ REAL, DIMENSION(SIZE(PPABST,1),SIZE(PPABST,2),SIZE(PPABST,3)) :: ZTHETAV, &
! MAE + DUR => Exner function perturbation
! LHE => Exner function perturbation * CPD * THVREF
!
+REAL :: ZPHI0
REAL :: ZRV_OV_RD ! XRV / XRD
REAL :: ZMAXVAL, ZMAXRES, ZMAX,ZMAX_ll ! for print
INTEGER, DIMENSION(3) :: IMAXLOC ! purpose
@@ -513,9 +517,15 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
ZPHIT(:,:,:)=(PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
!
ELSEIF(CEQNSYS=='LHE') THEN
- ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) &
- * XCPD * PTHVREF(:,:,:)
- !
+ IF ( .NOT. LOCEAN) THEN
+ ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) &
+ * XCPD * PTHVREF(:,:,:)
+ ELSE
+ ! Field at T- DT for LHE anelastic approx
+ ! not used in ocean case (flat LHE)
+ ZPHIT(:,:,:)=0.
+ END IF
+!
END IF
!
IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN .AND. .NOT. LIBM) THEN
@@ -781,12 +791,23 @@ IF ((ZMAX_ll > 1.E-12) .AND. KTCOUNT >0 ) THEN
!IF ( KTCOUNT >0 .AND. .NOT.LBOUSS ) THEN
CALL P_ABS ( KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
PTHT, PRT, PRHODJ, PRHODREF, ZTHETAV, PTHVREF, &
- PRVREF, PEXNREF, ZPHIT )
+ PRVREF, PEXNREF, ZPHIT, ZPHI0 )
!
IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
ELSEIF(CEQNSYS=='LHE') THEN
- PPABST(:,:,:)=XP00*(ZPHIT/(XCPD*PTHVREF)+PEXNREF)**(XCPD/XRD)
+ IF (.NOT. LOCEAN) THEN
+ ! Deep atmosphere case : computing of PI fluctuation ; ZPHI0 (computed in P_ABS routine) is added
+ XPHIT(:,:,:) = (ZPHIT+ZPHI0)/(XCPD*PTHVREF)
+ ! Absolute Pressure
+ PPABST(:,:,:)=XP00*(XPHIT(:,:,:)+PEXNREF)**(XCPD/XRD)
+ ! Computing press fluctuation
+ XPHIT(:,:,:) = PPABST(:,:,:) - XP00*PEXNREF**(XCPD/XRD)
+ ELSE
+! Shallow atmosphere ou ocean
+ XPHIT(:,:,:) = (ZPHIT+ZPHI0)*PRHODREF
+ PPABST(:,:,:)=XPHIT(:,:,:) + XP00*PEXNREF**(XCPD/XRD)
+ END IF
ENDIF
!
IF( HLBCX(1) == 'CYCL' ) THEN
diff --git a/src/MNH/read_field.f90 b/src/MNH/read_field.f90
index 4236b545058083091d043272a760ee8976c42953..e699bf02adbdbeb963970813fad79a9d30bdde82 100644
--- a/src/MNH/read_field.f90
+++ b/src/MNH/read_field.f90
@@ -9,7 +9,7 @@
!
INTERFACE
!
- SUBROUTINE READ_FIELD(TPINIFILE,KIU,KJU,KKU, &
+ SUBROUTINE READ_FIELD(KOCMI,TPINIFILE,KIU,KJU,KKU, &
HGETTKET,HGETRVT,HGETRCT,HGETRRT,HGETRIT,HGETCIT,HGETZWS, &
HGETRST,HGETRGT,HGETRHT,HGETSVT,HGETSRCT,HGETSIGS,HGETCLDFR, &
HGETBL_DEPTH,HGETSBL_DEPTH,HGETPHC,HGETPHR,HUVW_ADV_SCHEME, &
@@ -38,7 +38,7 @@ USE MODD_TIME ! for type DATE_TIME
!
!
TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE !Initial file
-INTEGER, INTENT(IN) :: KIU, KJU, KKU
+INTEGER, INTENT(IN) :: KIU, KJU, KKU,KOCMI
! array sizes in x, y and z directions
!
CHARACTER (LEN=*), INTENT(IN) :: HGETTKET, &
@@ -131,7 +131,7 @@ END INTERFACE
END MODULE MODI_READ_FIELD
!
! ########################################################################
- SUBROUTINE READ_FIELD(TPINIFILE,KIU,KJU,KKU, &
+ SUBROUTINE READ_FIELD(KOCEMI,TPINIFILE,KIU,KJU,KKU, &
HGETTKET,HGETRVT,HGETRCT,HGETRRT,HGETRIT,HGETCIT,HGETZWS, &
HGETRST,HGETRGT,HGETRHT,HGETSVT,HGETSRCT,HGETSIGS,HGETCLDFR, &
HGETBL_DEPTH,HGETSBL_DEPTH,HGETPHC,HGETPHR,HUVW_ADV_SCHEME, &
@@ -253,6 +253,7 @@ END MODULE MODI_READ_FIELD
! M. Leriche 10/06/2019: in restart case read all immersion modes for LIMA
!! F. Auguste 02/2021: add fields necessary for IBM
!! T. Nagel 02/2021: add fields necessary for turbulence recycling
+!! J.L. Redelsperger 03/2021: add necessary variables for Ocean LES case
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -269,6 +270,7 @@ USE MODD_CONF_n
USE MODD_CST
USE MODD_CTURB
USE MODD_DUST
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_ELEC_DESCR, ONLY: CELECNAMES
use modd_field, only: tfielddata, tfieldlist, TYPEDATE, TYPEREAL,TYPELOG,TYPEINT
USE MODD_FIELD_n, only: XZWS_DEFAULT
@@ -282,6 +284,7 @@ USE MODD_LATZ_EDFLX
USE MODD_LG, ONLY: CLGNAMES
USE MODD_LUNIT_N, ONLY: TLUOUT
USE MODD_NSV
+USE MODD_OCEANH
USE MODD_PARAM_C2R2, ONLY: LSUPSAT
!
USE MODD_PARAM_LIMA , ONLY: NMOD_CCN, LSCAV, LAERO_MASS, &
@@ -292,6 +295,7 @@ USE MODD_PARAM_n, ONLY: CSCONV
USE MODD_PASPOL
USE MODD_RAIN_C2R2_DESCR, ONLY: C2R2NAMES
USE MODD_RECYCL_PARAM_n
+USE MODD_REF, ONLY:LCOUPLES
USE MODD_SALT
USE MODD_TIME ! for type DATE_TIME
!
@@ -310,7 +314,7 @@ IMPLICIT NONE
!
!
TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE !Initial file
-INTEGER, INTENT(IN) :: KIU, KJU, KKU
+INTEGER, INTENT(IN) :: KIU, KJU, KKU,KOCEMI
! array sizes in x, y and z directions
!
CHARACTER (LEN=*), INTENT(IN) :: HGETTKET, &
@@ -1532,6 +1536,66 @@ END SELECT
!* 2.4 READ FORCING VARIABLES
! ----------------------
!
+! READ FIELD ONLY FOR MODEL1 (identical for all model in GN)
+IF (LOCEAN .AND. (.NOT.LCOUPLES) .AND. (KOCEMI==1)) THEN
+!
+ CALL IO_Field_read(TPINIFILE,'NFRCLT',NFRCLT)
+ CALL IO_Field_read(TPINIFILE,'NINFRT',NINFRT)
+!
+ TZFIELD%CMNHNAME = 'SSUFL_T'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'SSUFL'
+ TZFIELD%CUNITS = 'kg m-1 s-1'
+ TZFIELD%CDIR = '--'
+ TZFIELD%CCOMMENT = 'sfc stress along U to force ocean LES '
+ TZFIELD%NGRID = 0
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 1
+ TZFIELD%LTIMEDEP = .FALSE.
+ ALLOCATE(XSSUFL_T(NFRCLT))
+ CALL IO_Field_read(TPINIFILE,TZFIELD,XSSUFL_T(:))
+!
+ TZFIELD%CMNHNAME = 'SSVFL_T'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'SSVFL'
+ TZFIELD%CUNITS = 'kg m-1 s-1'
+ TZFIELD%CDIR = '--'
+ TZFIELD%CCOMMENT = 'sfc stress along V to force ocean LES '
+ TZFIELD%NGRID = 0
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 1
+ TZFIELD%LTIMEDEP = .FALSE.
+ALLOCATE(XSSVFL_T(NFRCLT))
+ CALL IO_Field_read(TPINIFILE,TZFIELD,XSSVFL_T(:))
+!
+ TZFIELD%CMNHNAME = 'SSTFL_T'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'SSTFL'
+ TZFIELD%CUNITS = 'kg m3 K m s-1'
+ TZFIELD%CDIR = '--'
+ TZFIELD%CCOMMENT = 'sfc total heat flux to force ocean LES '
+ TZFIELD%NGRID = 0
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 1
+ TZFIELD%LTIMEDEP = .FALSE.
+ ALLOCATE(XSSTFL_T(NFRCLT))
+ CALL IO_Field_read(TPINIFILE,TZFIELD,XSSTFL_T(:))
+!
+ TZFIELD%CMNHNAME = 'SSOLA_T'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'SSOLA'
+ TZFIELD%CUNITS = 'kg m3 K m s-1'
+ TZFIELD%CDIR = '--'
+ TZFIELD%CCOMMENT = 'sfc solar flux at sfc to force ocean LES '
+ TZFIELD%NGRID = 0
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 1
+ TZFIELD%LTIMEDEP = .FALSE.
+ ALLOCATE(XSSOLA_T(NFRCLT))
+ CALL IO_Field_read(TPINIFILE,TZFIELD,XSSOLA_T(:))
+!
+END IF ! ocean sfc forcing end
+
!
IF ( LFORCING ) THEN
DO JT=1,KFRC
@@ -1880,4 +1944,3 @@ END IF
!
!
END SUBROUTINE READ_FIELD
-
diff --git a/src/MNH/set_cstn.f90 b/src/MNH/set_cstn.f90
index 986526c82b5fd2812e3c20c9d54300592149a19e..22408fdc282d82461f1f57a711cdca15d1b17dac 100644
--- a/src/MNH/set_cstn.f90
+++ b/src/MNH/set_cstn.f90
@@ -170,6 +170,7 @@ END MODULE MODI_SET_CSTN
!
USE MODD_CONF
USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_GRID_n
USE MODD_IO, ONLY: TFILEDATA
USE MODD_LUNIT_n, ONLY: TLUOUT
@@ -393,8 +394,12 @@ END DO
!* 4.3 Compute Mixing ratio
!
! determines the pressure under the ground
-ZEXNGRDM= ( ZPGROUND / XP00) ** ZRDSCPD &
- - XG/XCPD / (0.5*(ZTHV(1)+ZTHVM(1))) * (ZZMASS_PROFILE(1) - ZHEIGHT(1))
+IF (LOCEAN) THEN
+ ZEXNGRDM= ( ZPGROUND / XP00) ** ZRDSCPD
+ELSE
+ ZEXNGRDM= ( ZPGROUND / XP00) ** ZRDSCPD &
+ - XG/XCPD / (0.5*(ZTHV(1)+ZTHVM(1))) * (ZZMASS_PROFILE(1) - ZHEIGHT(1))
+END IF
ZPGRDM = XP00 * ZEXNGRDM ** (1./ZRDSCPD)
ZPM(:) = PRESS_HEIGHT(ZZMASS_PROFILE(:),ZTHVM,ZPGRDM,ZTHVM(1),ZZMASS_PROFILE(1)) ! compute P
ZTVM(:) = ZTHVM(:) * (ZPM(:) / XP00) ** ZRDSCPD ! compute Tv
diff --git a/src/MNH/set_mass.f90 b/src/MNH/set_mass.f90
index a7b266aaa997e476fbc3fa3d32a820b47f56efbe..d38a9ab856a26ed7f9b1de607591743c09ea665a 100644
--- a/src/MNH/set_mass.f90
+++ b/src/MNH/set_mass.f90
@@ -121,6 +121,7 @@ SUBROUTINE SET_MASS(TPFILE,OPROFILE_IN_PROC, PZFLUX_PROFILE,
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! J-L Redelsperger 06/2021: Ocean case
!
!-------------------------------------------------------------------------------
!!
@@ -135,9 +136,11 @@ USE MODD_CST
USE MODD_REF
USE MODD_PARAMETERS
USE MODD_DIM_n
+USE MODD_DYN_n, ONLY : LOCEAN
!
USE MODE_GATHER_ll
USE MODE_ll
+USE MODE_THERMO
!
USE MODI_VER_INT_THERMO ! interface modules
USE MODI_WATER_SUM
@@ -147,6 +150,8 @@ USE MODI_VER_INT_DYN
USE MODI_SHUMAN
USE MODI_COMPUTE_EXNER_FROM_GROUND
USE MODI_COMPUTE_EXNER_FROM_TOP
+USE MODI_COMPUTE_PRESS_FROM_OCEANSFC
+USE MODI_COMPUTE_PRESS_FROM_OCEANBOT
USE MODI_SET_GEOSBAL
USE MODE_REPRO_SUM
USE MODE_MPPDB
@@ -191,7 +196,10 @@ REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZPMHP_MX ! pressu
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZRHOD_MX ! local rhod (mass level)
REAL,DIMENSION(SIZE(XZHAT)) :: ZRHOD_PROFILE ! local rhod (mass level) at initialization profile column
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZPMASS_MX ! pressure (mass level)
+REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZPFLUX_MX ! pressure (mass level)
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT)) :: ZEXNSURF2D_MX ! local Exner function at ground
+REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT)) :: ZPRESS2D_MX ! local pressure at ground
+REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT)) :: ZPRESSFC ! pressure at ocean sfc (ocen model case)
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZHEXNFLUX_MX ! local hyd. Exner function at flux points on the mixed grid
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZHEXNMASS_MX ! local hyd. Exner function at mass points on the mixed grid
!
@@ -224,6 +232,8 @@ REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZRHODJU ! horiz
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZRHODJV ! the MESONH Arakawa C grid
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZHEXNFLUX ! local hyd. Exner function at flux points (MNH grid)
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZHEXNMASS ! local hyd. Exner function at mass points (MNH grid)
+REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZPMASS ! local hyd. pres at mass points (MNH grid)
+REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZPFLUX ! local hyd. pres at flux points (MNH grid)
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZRHOD ! dry density on MESO-NH grid
!
!!$INTEGER :: IIBP,IIEP,IJBP,IJEP
@@ -279,17 +289,35 @@ ENDIF
!------------------------------
!* 2.2 compute exner function on mixed grid
!
-ZEXNSURF2D_MX(:,:)=(PPGROUND/XP00)**(XRD/XCPD)
-CALL COMPUTE_EXNER_FROM_GROUND(ZTHV3D_MX,PZFLUX_MX,&
- ZEXNSURF2D_MX,ZHEXNFLUX_MX,ZHEXNMASS_MX)
-ZEXNTOP2D(:,:)=ZHEXNFLUX_MX(:,:,IKE+1)
-ZPMASS_MX(:,:,:)=XP00*(ZHEXNMASS_MX(:,:,:))**(XCPD/XRD)
-ZRHOD_MX(:,:,:)=ZPMASS_MX(:,:,:)/(ZPMASS_MX(:,:,:)/XP00)**(XRD/XCPD) &
+ZEXNSURF2D_MX(:,:)=(PPGROUND/XP00)**(XRD/XCPD)
+!
+IF (LOCEAN) THEN
+ ZTHVREF3D(:,:,:) = ZTHV3D_MX(:,:,:)
+ ZRHOD_MX(:,:,:)= XRH00OCEAN*(1.-XALPHAOC*(ZTHV3D_MX(:,:,:)-XTH00OCEAN) &
+ +XBETAOC *(ZMR3D_MX(:,:,:,1)-XSA00OCEAN))
+ ZPRESS2D_MX(:,:)=PPGROUND
+ CALL COMPUTE_PRESS_FROM_OCEANBOT(ZRHOD_MX,PZFLUX_MX,ZPRESS2D_MX,ZPFLUX_MX,ZPMASS_MX)
+ ZHEXNFLUX_MX(:,:,:)=(ZPFLUX_MX(:,:,:)/XP00)**(XRD/XCPD)
+ ZHEXNMASS_MX(:,:,:)=(ZPMASS_MX(:,:,:)/XP00)**(XRD/XCPD)
+ ZEXNTOP2D(:,:)=ZHEXNFLUX_MX(:,:,IKE+1)
+ELSE
+ CALL COMPUTE_EXNER_FROM_GROUND(ZTHV3D_MX,PZFLUX_MX,&
+ ZEXNSURF2D_MX,ZHEXNFLUX_MX,ZHEXNMASS_MX)
+ ZEXNTOP2D(:,:)=ZHEXNFLUX_MX(:,:,IKE+1)
+ ZPMASS_MX(:,:,:)=XP00*(ZHEXNMASS_MX(:,:,:))**(XCPD/XRD)
+ENDIF
+!
+IF (LOCEAN) THEN
+ IF (LCOUPLES) THEN
+ XEXNTOPO=SUM_DD_R2_ll(ZHEXNFLUX_MX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
+ ELSE
+ XEXNTOP=SUM_DD_R2_ll(ZHEXNFLUX_MX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
+ END IF
+ELSE
+ ZRHOD_MX(:,:,:)=ZPMASS_MX(:,:,:)/(ZPMASS_MX(:,:,:)/XP00)**(XRD/XCPD) &
/(XRD*ZTHV3D_MX(:,:,:)*(1.+WATER_SUM(ZMR3D_MX(:,:,:,:))))
-
-XEXNTOP=SUM_DD_R2_ll(ZHEXNFLUX_MX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
-
-
+ XEXNTOP=SUM_DD_R2_ll(ZHEXNFLUX_MX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
+END IF
!------------------------------
!* 2.3 Rotate wind in model axis and take into account variations in x,y
! directions on the mixed grid
@@ -445,15 +473,17 @@ DEALLOCATE(ZNFLYZ_TOT,ZNFLYZ_TOT_ll)
!
IF (PRESENT(PCORIOZ)) THEN
- CALL SET_GEOSBAL(ZUW3D_FL,ZVW3D_FL,PTHVM,PMRM, &
+!To be modified later for ocean model case
+ CALL SET_GEOSBAL(ZUW3D_FL,ZVW3D_FL,PTHVM,PMRM, &
KILOC,KJLOC,OBOUSS,ZTHV3D,PCORIOZ)
CALL COMPUTE_EXNER_FROM_TOP(ZTHV3D,XZZ,ZEXNTOP2D,ZHEXNFLUX,ZHEXNMASS)
XPABSM(:,:,:)=XP00*ZHEXNMASS(:,:,:) ** (XCPD/XRD)
ELSE
-!
-! Interpolation of theta and r
!
- IF (SIZE(ZTHV3D_MX,3) > 3) THEN
+!No interpolation for ocean case (no bathimetry)
+! Interpolation of theta and r in atmos case
+!
+ IF (SIZE(ZTHV3D_MX,3) > 3) THEN
CALL VER_INT_THERMO(TPFILE,OSHIFT,ZTHV3D_MX,ZMR3D_MX,PZS_MX,PZS_MX,PZMASS_MX,&
PZFLUX_MX,ZPMHP_MX,ZEXNTOP2D, &
ZTHV3D,XRT,ZPMHP,ZDIAG)
@@ -462,7 +492,11 @@ ELSE
XRT = ZMR3D_MX
ZDIAG = 0.
END IF
- XTHT(:,:,:)=ZTHV3D(:,:,:)*(1.+WATER_SUM(XRT(:,:,:,:)))/(1.+XRV/XRD*XRT(:,:,:,1))
+ IF (LOCEAN) THEN
+ XTHT(:,:,:)=ZTHV3D(:,:,:)
+ ELSE
+ XTHT(:,:,:)=ZTHV3D(:,:,:)*(1.+WATER_SUM(XRT(:,:,:,:)))/(1.+XRV/XRD*XRT(:,:,:,1))
+ ENDIF
ZTHV3D(:,:,1)=ZTHV3D(:,:,2)
XTHT(:,:,1)=XTHT(:,:,2)
XRT(:,:,1,:)=XRT(:,:,2,:)
@@ -472,7 +506,6 @@ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHV3D, 'SET_MASS::ZTHV3D' )
CALL ADD3DFIELD_ll( TZFIELDS_ll, XRT(:,:,1,:), 'SET_MASS::XRT(:,:,1,:)' )
CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS_ll)
-
!
IF (NRR>=3) THEN
WHERE (XRT(:,:,:,3)<1.E-20)
@@ -489,17 +522,19 @@ CALL CLEANLIST_ll(TZFIELDS_ll)
CALL VER_INT_DYN(OSHIFT,ZRHODU_MX,ZRHODV_MX,PZFLUX_MX,PZMASS_MX,PZS_MX,ZRHODUA,ZRHODVA)
ZRHODJU(:,:,:)=MXM(ZRHODUA(:,:,:)*PJ(:,:,:))
ZRHODJV(:,:,:)=MYM(ZRHODVA(:,:,:)*PJ(:,:,:))
- CALL COMPUTE_EXNER_FROM_TOP(ZTHV3D,XZZ,ZEXNTOP2D,ZHEXNFLUX,ZHEXNMASS)
- XPABST(:,:,:)=ZPMHP(:,:,:) + XP00*ZHEXNMASS(:,:,:) ** (XCPD/XRD)
- ZRHOD(:,:,:)=XPABST(:,:,:)/(XPABST(:,:,:)/XP00)**(XRD/XCPD) &
- /(XRD*XTHT(:,:,:)*(1.+XRV/XRD*XRT(:,:,:,1)))
+ IF (.NOT.LOCEAN) THEN
+ CALL COMPUTE_EXNER_FROM_TOP(ZTHV3D,XZZ,ZEXNTOP2D,ZHEXNFLUX,ZHEXNMASS)
+ XPABST(:,:,:)=ZPMHP(:,:,:) + XP00*ZHEXNMASS(:,:,:) ** (XCPD/XRD)
+ ZRHOD(:,:,:)=XPABST(:,:,:)/(XPABST(:,:,:)/XP00)**(XRD/XCPD) /(XRD*XTHT(:,:,:)*(1.+XRV/XRD*XRT(:,:,:,1)))
+ ELSE
+ ZRHOD(:,:,:)=XRH00OCEAN*(1.-XALPHAOC*(XTHT(:,:,:)-XTH00OCEAN)+XBETAOC*(XRT(:,:,:,1)-XSA00OCEAN))
+ END IF
XUT(:,:,:)=ZRHODJU(:,:,:)/MXM(ZRHOD(:,:,:)*PJ(:,:,:))
XVT(:,:,:)=ZRHODJV(:,:,:)/MYM(ZRHOD(:,:,:)*PJ(:,:,:))
XWT(:,:,:)=0
CALL MPPDB_CHECK3DM("SET_MASS:XVT,ZRHODJV,PJ,ZRHODVA",PRECISION,&
& XVT,ZRHODJV,PJ,ZRHODVA )
ENDIF
-
!
!-------------------------------------------------------------------------------
!* 4. COMPUTE ANELASTIC REFERENCE (PV)
@@ -518,17 +553,32 @@ ELSE
DO JK = 1,IKU
CALL REDUCESUM_ll(XTHVREFZ(JK), IINFO_ll)
END DO
-
- XRHODREFZ(:) = XP00/ (XRD* XTHVREFZ(:))
- ZTHVREF3D(:,:,:)=XTHVREFZ(2)
- CALL COMPUTE_EXNER_FROM_GROUND(ZTHVREF3D,PZFLUX_MX,&
- ZEXNSURF2D_MX,ZHEXNFLUX,ZHEXNMASS)
-
- XEXNTOP=SUM_DD_R2_ll(ZHEXNFLUX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
-
- ZEXNTOP2D=ZHEXNFLUX(:,:,IKE+1)
- CALL COMPUTE_EXNER_FROM_TOP(ZTHVREF3D,XZZ,ZEXNTOP2D,ZHEXNFLUX,ZHEXNMASS)
- XPABST(:,:,:)=ZPMHP(:,:,:) + XP00*ZHEXNMASS(:,:,:) ** (XCPD/XRD)
+!
+ IF (LOCEAN) THEN
+! Ocean case boussinesq
+ IF (LCOUPLES) THEN
+ XRHODREFZO(:) = XRH00OCEAN
+ XTHVREFZ(:) = ZTHV3D(KILOC,KJLOC,IKU-3) ! XTHVREFZ is uniform
+ ZTHVREF3D(:,:,:)=XTHVREFZ(IKU-3)
+ ZPRESSFC(:,:)=XP00*XEXNTOPO**(XCPD/XRD)
+ ELSE
+ XRHODREFZ(:) = XRH00OCEAN
+ ZPRESSFC(:,:)=XP00*XEXNTOP**(XCPD/XRD)
+! on prend pour le moment la valeur de la couche mélangée
+ END IF
+ CALL COMPUTE_PRESS_FROM_OCEANSFC(ZRHOD,XZZ,ZPRESSFC,ZPFLUX,ZPMASS)
+ XPABST(:,:,:)= ZPMASS(:,:,:)
+!
+ ELSE
+! ATmos: rho = P/ (R Tv)
+ XRHODREFZ(:) = XP00/ (XRD* XTHVREFZ(:))
+ ZTHVREF3D(:,:,:)=XTHVREFZ(2)
+ XEXNTOP=SUM_DD_R2_ll(ZHEXNFLUX(IIB:IIE,IJB:IJE,IKE+1))/REAL(NIMAX_ll*NJMAX_ll)
+ ZEXNTOP2D=ZHEXNFLUX(:,:,IKE+1)
+ CALL COMPUTE_EXNER_FROM_TOP(ZTHVREF3D,XZZ,ZEXNTOP2D,ZHEXNFLUX,ZHEXNMASS)
+ XPABST(:,:,:)=ZPMHP(:,:,:) + XP00*ZHEXNMASS(:,:,:) ** (XCPD/XRD)
+ENDIF
+ ! end of bouss case
ENDIF
!---------------------------------------------------------------------------------
END SUBROUTINE SET_MASS
diff --git a/src/MNH/set_perturb.f90 b/src/MNH/set_perturb.f90
index 1f88a63cfff307d1dd499afd5f73d1d7d1f51d9e..ead3d16e74e5bfb12696e98a17ffebd7e60dfcd6 100644
--- a/src/MNH/set_perturb.f90
+++ b/src/MNH/set_perturb.f90
@@ -98,7 +98,8 @@ END MODULE MODI_SET_PERTURB
!! C.Lac, V.Masson 1/2018 : White noise in the LBC
!! Q.Rodier 10/2018 : move allocate(ZWHITE) for NKWH>2
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
-!
+!! J.L Redelsperger 03/2021 : : white noise in Ocean LES case at the top of domain(Sfc)
+!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -106,6 +107,7 @@ END MODULE MODI_SET_PERTURB
!
USE MODD_CST
USE MODD_CONF
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_DIM_n
USE MODD_FIELD_n
USE MODD_GRID_n
@@ -115,6 +117,7 @@ USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_LSFIELD_n
USE MODD_PARAMETERS
USE MODD_REF_n
+USE MODD_REF
USE MODD_VAR_ll , ONLY : NMNH_COMM_WORLD
!
USE MODE_GATHER_ll
@@ -196,6 +199,7 @@ INTEGER, DIMENSION(:), ALLOCATABLE :: i_seed
INTEGER :: ni_seed
!
INTEGER :: IXOR,IYOR,JI_ll,JJ_ll
+INTEGER :: INOISB,INOISE ! Loop indice for White noise
!
NAMELIST/NAM_PERT_PRE/CPERT_KIND,XAMPLITH, &! Perturbation parameters
XAMPLIRV,XCENTERZ,XRADX, &!
@@ -249,6 +253,13 @@ IJE_ll=IJU_ll-JPHEXT
!
CALL GET_OR_ll('B',IXOR,IYOR)
!
+IF (LOCEAN) THEN
+ INOISB=NKWH
+ INOISE=IKE
+ELSE
+ INOISB=2
+ INOISE=NKWH
+ENDIF
!-------------------------------------------------------------------------------
!
!* 2. COMPUTE THE PERTURBATION ON THETA :
@@ -375,8 +386,7 @@ SELECT CASE(CPERT_KIND)
! J.Escobar optim => need only identical random on all domain
!
ALLOCATE(ZWHITE(IIU,IJU))
-!
- DO JK = IKB, NKWH
+ DO JK = INOISB,INOISE
IKX = (NIMAX_ll+1)/2
ZX = 2*XPI/NIMAX_ll
ALLOCATE(ZCX_ll(IIU_ll,IKX))
diff --git a/src/MNH/set_ref.f90 b/src/MNH/set_ref.f90
index 925fd52efdf6c0189af56d44e594453beaf3bd85..5e355a10e401920849c8a6a738f82a8ff39a431e 100644
--- a/src/MNH/set_ref.f90
+++ b/src/MNH/set_ref.f90
@@ -150,6 +150,7 @@ END MODULE MODI_SET_REF
!! PRHODREF, PEXNREF, PTHVREF after computation
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! Jean-Luc Redelsperger 03/2021 : OCEAN LES case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -160,9 +161,11 @@ USE MODD_IO, ONLY: TFILEDATA
USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_PARAMETERS
USE MODD_REF
+USE MODD_DYN_n, ONLY : LOCEAN
!
USE MODE_IO_FIELD_READ, only: IO_Field_read
USE MODE_ll
+USE MODE_THERMO
USE MODE_MPPDB
USE MODE_REPRO_SUM
!
@@ -210,6 +213,7 @@ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZRHOREF
! Reference density
REAL, DIMENSION(SIZE(PZZ,3)) :: ZZHATM ! height of the mass levels
! in the transformed space (GCS transf.) or without orography
+REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZDENSOC,ZPFLUX,ZPMASS
!
INTEGER :: IIU ! Upper dimension in x direction
INTEGER :: IJU ! Upper dimension in y direction
@@ -263,7 +267,12 @@ IF (KMI == 1) THEN
LNEUTRAL=.FALSE.
IF (MAXVAL(XTHVREFZ(IKB:IKE))-MINVAL(XTHVREFZ(IKB:IKE)) < 1.E-10) LNEUTRAL=.TRUE.
END IF
-!
+!* Ref state diff for O & A in LES coupled mode
+IF (LCOUPLES .AND. LOCEAN) THEN
+ CALL IO_Field_read(TPINIFILE,'RHOREFZ',XRHODREFZO)
+ CALL IO_Field_read(TPINIFILE,'THVREFZ',XTHVREFZO)
+ CALL IO_Field_read(TPINIFILE,'EXNTOP', XEXNTOPO)
+END IF
!-------------------------------------------------------------------------------
!
!* 3. SET REFERENCE STATE WITH OROGRAPHY
@@ -285,7 +294,13 @@ CALL MPPDB_CHECK3D(ZZM,"SET_REF::ZZM",PRECISION)
!
!* 3.2 Interpolation
!
-DO JI = 1,SIZE(PZZ,1)
+IF (LCOUPLES .AND. LOCEAN) THEN
+ DO JK = 1,IKU
+ PTHVREF(:,:,JK) = XTHVREFZO(JK)
+ PRHODREF(:,:,JK)= XRHODREFZO(JK)
+ END DO
+ELSE
+ DO JI = 1,SIZE(PZZ,1)
DO JJ = 1,SIZE(PZZ,2)
!
DO JK = 1,IKU
@@ -316,22 +331,21 @@ DO JI = 1,SIZE(PZZ,1)
END IF
END DO
END DO
-END DO
+ END DO
+END IF
!
! change the extrapolation option for the thvref field to be consistent with
! the extrapolation option for the flottability at the ground and for rhodref
! to be consistent with the extrapolation to compute a divergence
PTHVREF(:,:,IKB-1) = PTHVREF(:,:,IKB)
PRHODREF(:,:,IKB-1) = PRHODREF(:,:,IKB)
-
CALL MPPDB_CHECK3D(PTHVREF,"SET_REF::PTHVREF",PRECISION)
CALL MPPDB_CHECK3D(PRHODREF,"SET_REF::PRHODREF",PRECISION)
!
!-------------------------------------------------------------------------------
!
-!* 4. COMPUTE EXNER FUNCTION
-! ----------------------
-!
+!* 4. COMPUTE EXNER FUNCTION AT MASS GRID POINT
+! ----------------------------------------
IF (LCARTESIAN .OR. LTHINSHELL) THEN
ZD1=0.
ELSE
@@ -340,24 +354,56 @@ ENDIF
!
ZGSCPD = XG/XCPD
!
-PEXNREF(:,:,IKE)=(XEXNTOP*(1.+ZD1*2./7.*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE))/ &
- (XRADIUS+(PZZ(:,:,IKE+1)+ZZM(:,:,IKE))/2.)) &
- + ZGSCPD/PTHVREF(:,:,IKE)*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE)))/ &
-(1.-ZD1*2./7.*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE))/(XRADIUS+(PZZ(:,:,IKE+1)+ZZM(:,:,IKE))/2.))
+IF (LOCEAN) THEN
+!--------------------------------
+! Pressure at domain top (Flux point !!!) saved in Press_mass above the ocen sfc
+ IF (LCOUPLES) THEN
+ ZPMASS(:,:,IKE+1)= XP00 *XEXNTOPO**(XCPD/XRD)
+ ELSE
+ ZPMASS(:,:,IKE+1)= XP00 *XEXNTOP**(XCPD/XRD)
+ ENDIF
+ ZPMASS(:,:,IKE) = ZPMASS(:,:,IKE+1) +XG*PRHODREF(:,:,IKE)*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE))
+ DO JK = IKE-1,1,-1
+ ZPMASS(:,:,JK) = ZPMASS(:,:,JK+1) + XG * &
+ .5*(PRHODREF(:,:,JK)+ PRHODREF(:,:,JK+1)) * (ZZM(:,:,JK+1) -ZZM(:,:,JK))
+ END DO
+!
+ IF (LCOUPLES) THEN
+ DO JK = IKE+1, IKU
+! Pressure above domain top (i.e. ocean sfc), i.e. in atmosphere (should be not used)
+ ZPMASS(:,:,JK) = XP00 *XEXNTOPO**(XCPD/XRD)
+ END DO
+ ELSE
+ DO JK = IKE+1, IKU
+! Pressure above domain top (i.e. ocean sfc), i.e. in atmosphere (should be not used)
+ ZPMASS(:,:,JK) = XP00 *XEXNTOP**(XCPD/XRD)
+ END DO
+ ENDIF
+ PEXNREF(:,:,:)= (ZPMASS(:,:,:)/XP00)**(XRD/XCPD)
+ ! OCEAN end
+ELSE
+ ! ATMOSPHERE
+ PEXNREF(:,:,IKE)=(XEXNTOP*(1.+ZD1*2./7.*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE))/ &
+ (XRADIUS+(PZZ(:,:,IKE+1)+ZZM(:,:,IKE))/2.)) &
+ + ZGSCPD/PTHVREF(:,:,IKE)*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE)))/ &
+ (1.-ZD1*2./7.*(PZZ(:,:,IKE+1)-ZZM(:,:,IKE))/(XRADIUS+(PZZ(:,:,IKE+1)+ZZM(:,:,IKE))/2.))
+!
+ DO JK = IKE-1, 1, -1
+ PEXNREF(:,:,JK)=(PEXNREF(:,:,JK+1)*(1.+ZD1*2./7.*(ZZM(:,:,JK+1) -ZZM(:,:,JK))/ &
+ (XRADIUS+PZZ(:,:,JK+1)))+ &
+ 2.*ZGSCPD/(PTHVREF(:,:,JK+1)+PTHVREF(:,:,JK))*(ZZM(:,:,JK+1) -ZZM(:,:,JK)))/&
+ (1.-ZD1*2./7.*(ZZM(:,:,JK+1) -ZZM(:,:,JK))/(XRADIUS+PZZ(:,:,JK+1)))
+ END DO
+!
+ DO JK = IKE+1, IKU
+ PEXNREF(:,:,JK)=(PEXNREF(:,:,JK-1)*(1.+ZD1*2./7.*(ZZM(:,:,JK-1) -ZZM(:,:,JK))/ &
+ (XRADIUS+PZZ(:,:,JK)))+ &
+ 2.*ZGSCPD/(PTHVREF(:,:,JK-1)+PTHVREF(:,:,JK))*(ZZM(:,:,JK-1) -ZZM(:,:,JK)))/&
+ (1.-ZD1*2./7.*(ZZM(:,:,JK-1) -ZZM(:,:,JK))/ (XRADIUS+PZZ(:,:,JK)))
+ END DO
!
-DO JK = IKE-1, 1, -1
- PEXNREF(:,:,JK)=(PEXNREF(:,:,JK+1)*(1.+ZD1*2./7.*(ZZM(:,:,JK+1) -ZZM(:,:,JK))/ &
- (XRADIUS+PZZ(:,:,JK+1)))+ &
- 2.*ZGSCPD/(PTHVREF(:,:,JK+1)+PTHVREF(:,:,JK))*(ZZM(:,:,JK+1) -ZZM(:,:,JK)))/&
- (1.-ZD1*2./7.*(ZZM(:,:,JK+1) -ZZM(:,:,JK))/(XRADIUS+PZZ(:,:,JK+1)))
-END DO
+END IF
!
-DO JK = IKE+1, IKU
- PEXNREF(:,:,JK)=(PEXNREF(:,:,JK-1)*(1.+ZD1*2./7.*(ZZM(:,:,JK-1) -ZZM(:,:,JK))/ &
- (XRADIUS+PZZ(:,:,JK)))+ &
- 2.*ZGSCPD/(PTHVREF(:,:,JK-1)+PTHVREF(:,:,JK))*(ZZM(:,:,JK-1) -ZZM(:,:,JK)))/&
- (1.-ZD1*2./7.*(ZZM(:,:,JK-1) -ZZM(:,:,JK))/ (XRADIUS+PZZ(:,:,JK)))
-END DO
!
CALL MPPDB_CHECK3D(PEXNREF,"SET_REF::PEXNREF",PRECISION)
!-------------------------------------------------------------------------------
@@ -372,8 +418,8 @@ IF (LBOUSS) THEN
ELSE
ZRHOREF(:,:,:) = PEXNREF(:,:,:) ** ZCVD_O_RD * XP00 / ( XRD * PTHVREF(:,:,:) )
ZRHOREF(:,:,1)=ZRHOREF(:,:,2) ! this avoids to obtain erroneous values for
+ ! rv at this last point
END IF
- ! rv at this last point
!
IF ( CEQNSYS == 'DUR' ) THEN
IF ( SIZE(PRVREF,1) == 0 ) THEN
@@ -402,13 +448,10 @@ CALL CLEANLIST_ll(TZFIELDS_ll)
CALL MPPDB_CHECK3D(ZRHOREF,"SET_REF::ZRHOREF",PRECISION)
IF ( SIZE(PRVREF,1) /= 0 ) CALL MPPDB_CHECK3D(PRVREF,"SET_REF::PRVREF",PRECISION)
CALL MPPDB_CHECK3D(PRHODJ,"SET_REF::PRHODJ",PRECISION)
-
-
!
!* 6. COMPUTES THE TOTAL MASS OF REFERENCE ATMOSPHERE
! -----------------------------------------------
!
-IF (CEQNSYS == "LHE" ) THEN
ZCVD_O_RDCPD = ZCVD_O_RD / XCPD
!
ALLOCATE(ZREFMASS_2D(IIB:IIE,IJB:IJE))
@@ -432,7 +475,6 @@ IF (CEQNSYS == "LHE" ) THEN
PMASS_O_PHI0 = SUM_DD_R2_ll(ZMASS_O_PHI0_2D)
!JUAN16
!
-END IF
!
!
!-------------------------------------------------------------------------------
@@ -480,7 +522,7 @@ IF ( HLBCY(1)=='OPEN' ) THEN
ENDDO
ENDIF
PLINMASS = PLINMASS + SUM_DD_R2_ll(ZLINMASS_S_2D)
-!
+ !
ALLOCATE( ZLINMASS_N_2D(IIB:IIE,IJE+1:IJE+1))
ZLINMASS_N_2D = 0.0
IF (LNORTH_ll(HSPLITTING='B')) THEN
diff --git a/src/MNH/set_rsou.f90 b/src/MNH/set_rsou.f90
index 353c6298bd02baa8db936c4df60565fc4235c013..8f42147d1318ebc51d7182ac1a486b19f5bac80d 100644
--- a/src/MNH/set_rsou.f90
+++ b/src/MNH/set_rsou.f90
@@ -102,7 +102,17 @@ END MODULE MODI_SET_RSOU
! (Pressure, U, V) ,
! (Pressure, T, Hu)
!
-!
+! For ocean-LES case the following kind of data is permitted
+!
+! YKIND = 'IDEALOCE' : ZGROUND (Water depth),PGROUND(Sfc Atmos Press),
+! TGROUND (SST), RGROUND (SSS)
+! (Depth , U, V) starting from sfc
+! (Depth, T, S)
+! (Time, LE, H, SW_d,SW_u,LW_d,LW_u,Stress_X,Stress_Y)
+!
+! YKIND = 'STANDOCE' : (Depth , Temp, Salinity, U, V) starting from sfc
+! (Time, LE, H, SW_d,SW_u,LW_d,LW_u,Stress_X,Stress_Y)
+!
!!** METHOD
!! ------
!! The radiosounding is first read, then data are converted in order to
@@ -242,6 +252,7 @@ END MODULE MODI_SET_RSOU
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 19/04/2019: removed unused dummy arguments and variables
+!! JL Redelsperger : 01/2021: Ocean LES cases added
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -250,11 +261,16 @@ END MODULE MODI_SET_RSOU
USE MODD_CONF
USE MODD_CONF_n
USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_FIELD_n
USE MODD_GRID
USE MODD_GRID_n
+USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_IO, ONLY: TFILEDATA
+USE MODD_NETCDF
+USE MODD_OCEANH
USE MODD_PARAMETERS, ONLY: JPHEXT
+USE MODD_TYPE_DATE
!
USE MODE_ll
USE MODE_MSG
@@ -269,6 +285,8 @@ USE MODI_THETAVPU_THETAVPM
USE MODI_TH_R_FROM_THL_RT_1D
USE MODI_VERT_COORD
!
+USE NETCDF ! for reading the NR files
+!
IMPLICIT NONE
!
!
@@ -291,10 +309,15 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PJ ! jacobien
!
!* 0.2 Declarations of local variables :
!
-INTEGER :: ILUPRE ! logical unit number
-!
-! variables read in EXPRE file at the RS levels
-!
+INTEGER :: ILUPRE ! logical unit number of the EXPRE return code
+INTEGER :: ILUOUT ! Logical unit number for output-listing
+! local variables for reading sea sfc flux forcing for ocean model
+INTEGER :: IFRCLT
+REAL, DIMENSION(:), ALLOCATABLE :: ZSSUFL_T,ZSSVFL_T,ZSSTFL_T,ZSSOLA_T !
+TYPE (DATE_TIME), DIMENSION(:), ALLOCATABLE :: ZFRCLT ! date/time of sea surface forcings
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! variables read in EXPRE file at the RS/CTD levels
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CHARACTER(LEN=8) :: YKIND ! Kind of variables in
! EXPRE FILE
INTEGER :: ILEVELU ! number of wind levels
@@ -332,7 +355,7 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZLSOCPEXN
REAL, DIMENSION(SIZE(XZHAT)) :: ZZFLUX_PROFILE ! altitude of flux points on the initialization columns
REAL, DIMENSION(SIZE(XZHAT)) :: ZZMASS_PROFILE ! altitude of mass points on the initialization columns
!
-! fieds on the grid of the model without orography
+! fields on the grid of the model without orography
!
REAL, DIMENSION(SIZE(XZHAT)) :: ZUW,ZVW ! Wind at w model grid levels
REAL, DIMENSION(SIZE(XZHAT)) :: ZMRM ! vapor mixing ratio at mass model
@@ -341,29 +364,42 @@ REAL, DIMENSION(SIZE(XZHAT)) :: ZMRCM,ZMRIM
REAL, DIMENSION(SIZE(XZHAT)) :: ZTHVM ! Temperature at mass model grid levels
REAL, DIMENSION(SIZE(XZHAT)) :: ZTHLM ! Thetal at mass model grid levels
REAL, DIMENSION(SIZE(XZHAT)) :: ZTHM ! Thetal at mass model grid levels
+REAL, DIMENSION(SIZE(XZHAT)) :: ZRHODM ! density at mass model grid level
REAL, DIMENSION(:), ALLOCATABLE :: ZMRT ! Total Vapor mixing ratio at mass levels on mixed grid
REAL, DIMENSION(:), ALLOCATABLE :: ZEXNMASS ! exner fonction at mass level
REAL, DIMENSION(:), ALLOCATABLE :: ZEXNFLUX ! exner fonction at flux level
REAL :: ZEXNSURF ! exner fonction at surface
+REAL, DIMENSION(:), ALLOCATABLE :: ZPREFLUX ! pressure at flux model grid level
REAL, DIMENSION(:), ALLOCATABLE :: ZFRAC_ICE ! ice fraction
REAL, DIMENSION(:), ALLOCATABLE :: ZRSATW, ZRSATI
REAL :: ZDZSDH,ZDZ1SDH,ZDZ2SDH ! interpolation
! working arrays
!
-INTEGER :: JK,JKLEV, JKU,JKM ! Loop indexes
+INTEGER :: JK,JKLEV,JKU,JKM,JKT,JJ,JI,JO,JLOOP ! Loop indexes
INTEGER :: IKU ! Upper bound in z direction
REAL :: ZRDSCPD,ZRADSDG, & ! Rd/Cpd, Pi/180.,
- ZRVSRD,ZRDSRV ! Rv/Rd, Rd/Rv
+ ZRVSRD,ZRDSRV, & ! Rv/Rd, Rd/Rv
+ ZPTOP ! Pressure at domain top
LOGICAL :: GUSERC ! use of input data cloud
-INTEGER :: IIB, IIE, IJB, IJE
+INTEGER :: IIB, IIE, IJB, IJE,IKB,IKE
INTEGER :: IXOR_ll, IYOR_ll
INTEGER :: IINFO_ll
LOGICAL :: GPROFILE_IN_PROC ! T : initialization profile is in current processor
+CHARACTER(LEN=100) :: YMSG
!
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT)) ::ZZS_LS
REAL,DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) ::ZZFLUX_MX,ZZMASS_MX ! mixed grid
-INTEGER :: JLOOP
!-------------------------------------------------------------------------------
+! For standard ocean version, reading external files
+CHARACTER(LEN=256) :: infile,infisf ! files to be read
+INTEGER :: NZ,NLATI,NLONGI,IDX,ITIME,ncid,varid,ndimp,ntcindy
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZOC_TEMPERATURE,ZOC_SALINITY,ZOC_U,ZOC_V
+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
+
+!--------------------------------------------------------------------------------
!
!* 1. PROLOGUE : INITIALIZE SOME CONSTANTS, RETRIEVE LOGICAL
! UNIT NUMBERS AND READ KIND OF DATA IN EXPRE FILE
@@ -380,31 +416,21 @@ ZRVSRD = XRV/XRD
ZRDSRV = XRD/XRV
!
!* 1.2 Retrieve logical unit numbers
-!
-!
+!
ILUPRE = TPEXPREFILE%NLU
+ILUOUT = TLUOUT%NLU
!
!* 1.3 Read data kind in EXPRE file
!
READ(ILUPRE,*) YKIND
-!
+WRITE(ILUOUT,*) 'YKIND read in set_rsou: ', YKIND
!
IF(LUSERC .AND. YKIND/='PUVTHDMR' .AND. YKIND/='ZUVTHDMR' .AND. YKIND/='ZUVTHLMR') THEN
CALL PRINT_MSG(NVERB_FATAL,'GEN','SET_RSOU','hydrometeors are not allowed for YKIND = '//trim(YKIND))
ENDIF
-! Demande Thierry Bergot Sept 2012
-!IF(LUSERC .AND.(YKIND == 'PUVTHDMR' .OR. YKIND == 'ZUVTHDMR').AND. .NOT. L1D) THEN
-! !callabortstop
-! CALL PRINT_MSG(NVERB_FATAL,'GEN','SET_RSOU','use of hydrometeors for YKIND=P(Z)UVTHDMR is only allowed in 1D case')
-!ENDIF
-!
-!IF(LUSERI .AND. YKIND=='ZUVTHLMR') THEN
-! !callabortstop
-! CALL PRINT_MSG(NVERB_FATAL,'GEN','SET_RSOU','use of ice for YKIND=ZUVTHLMR is not allowed')
-!ENDIF
!
IF(YKIND=='ZUVTHLMR' .AND. .NOT. LUSERC) THEN
- !callabortstop
+!callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','SET_RSOU','LUSERC=T is required for YKIND=ZUVTHLMR')
ENDIF
!
@@ -416,8 +442,345 @@ IF(LUSERC .AND. (YKIND == 'PUVTHDMR' .OR. YKIND == 'ZUVTHDMR')) GUSERC=.TRUE.
! --------------------------------------------------------
!
SELECT CASE(YKIND)
+!
+! 2.0.1 Ocean case 1
+!
+ CASE ('IDEALOCE')
+!
+ ! Read data in PRE_IDEA1.nam
+ ! Surface
+ WRITE(ILUOUT,FMT=*) 'Reading data for ideal ocean :IDEALOCE'
+ READ(ILUPRE,*) ZPTOP ! P_atmosphere at sfc =P top domain
+ READ(ILUPRE,*) ZTGROUND ! SST
+ READ(ILUPRE,*) ZMRGROUND ! SSS
+ WRITE(ILUOUT,FMT=*) 'Patm SST SSS', ZPTOP,ZTGROUND,ZMRGROUND
+ READ(ILUPRE,*) ILEVELU ! Read number of Current levels
+ ! Allocate required memory
+ ALLOCATE(ZHEIGHTU(ILEVELU),ZU(ILEVELU),ZV(ILEVELU))
+ ALLOCATE(ZOC_U(ILEVELU,1,1),ZOC_V(ILEVELU,1,1))
+ WRITE(ILUOUT,FMT=*) 'Level number for Current in data', ILEVELU
+ ! Read U and V at each wind level
+ DO JKU = 1,ILEVELU
+ READ(ILUPRE,*) ZHEIGHTU(JKU),ZOC_U(JKU,1,1),ZOC_V(JKU,1,1)
+ ! WRITE(ILUOUT,FMT=*) 'Leveldata D(m) under sfc: U_cur, V_cur', JKU, ZHEIGHTU(JKU),ZU(JKU),ZV(JKU)
+ END DO
+ DO JKU=1,ILEVELU
+ ! Z axis reoriented as in the model
+ IDX = ILEVELU-JKU+1
+ ZU(JKU) = ZOC_U(IDX,1,1)
+ ZV(JKU) = ZOC_V(IDX,1,1)
+ ! 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)
+ END DO
+ ! Read number of mass levels
+ READ(ILUPRE,*) ILEVELM
+ ! Allocate required memory
+ ALLOCATE(ZOC_DEPTH(ILEVELM))
+ ALLOCATE(ZHEIGHTM(ILEVELM))
+ ALLOCATE(ZTHL(ILEVELM),ZTH(ILEVELM),ZTHV(ILEVELM))
+ ALLOCATE(ZMR(ILEVELM),ZRT(ILEVELM))
+ ALLOCATE(ZOC_TEMPERATURE(ILEVELM,1,1),ZOC_SALINITY(ILEVELM,1,1))
+ ! Read T and S at each mass level
+ DO JKM= 2,ILEVELM
+ READ(ILUPRE,*) ZOC_DEPTH(JKM),ZOC_TEMPERATURE(JKM,1,1),ZOC_SALINITY(JKM,1,1)
+ END DO
+ ! Complete the mass arrays with the ground informations read in EXPRE file
+ ZOC_DEPTH(1) = 0.
+ ZOC_TEMPERATURE(1,1,1)= ZTGROUND
+ ZOC_SALINITY(1,1,1)= ZMRGROUND
+ !!!!!!!!!!!!!!!!!!!!!!!!Inversing Axis!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Going from the data (axis downward i.e inverse model) grid to the model grid (axis upward)
+ ! Uniform bathymetry; depth goes from ocean sfc downwards (data grid)
+ ! ZHEIGHT goes from the model domain bottom up to the sfc ocean (top of model domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ZZGROUND = 0.
+ ZTGROUND = ZOC_TEMPERATURE(ILEVELM,1,1)
+ ZMRGROUND = ZOC_SALINITY(ILEVELM,1,1)
+ DO JKM= 1,ILEVELM
+ ! Z upward axis (oriented as in the model), i.e.
+ ! going from 0m (ocean bottom/model bottom) upward to H (ocean sfc/model top)
+ ! ZHEIGHT used only in set_ rsou, defined as such ZHEIGHT(ILEVELM)=H_model
+ IDX = ILEVELM-JKM+1
+ ZTH(JKM) = ZOC_TEMPERATURE(IDX,1,1)
+ ZMR(JKM) = ZOC_SALINITY(IDX,1,1)
+ ZHEIGHTM(JKM)= ZOC_DEPTH(ILEVELM)- ZOC_DEPTH(IDX)
+ WRITE(ILUOUT,FMT=*) 'Model oriented initial data: JKM IDX depth T S ZHEIGHTM', &
+ JKM,IDX,ZOC_DEPTH(IDX),ZTH(JKM),ZMR(JKM),ZHEIGHTM(JKM)
+ END DO
+ ! mass levels of the RS
+ ZTHV = ZTH ! TV==THETA=TL
+ ZTHL = ZTH
+ ZRT = ZMR
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! READ Sea Surface Forcing !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Reading the forcings from prep_idea1.nam
+ READ(ILUPRE,*) IFRCLT ! Number of time-dependent forcing
+ IF (IFRCLT > 99*8) THEN
+ ! CAUTION: number of forcing times is limited by the WRITE format 99(8E10.3)
+ ! and also by the name of forcing variables (format I3.3)
+ ! You have to modify those if you need more forcing times
+ CALL PRINT_MSG(NVERB_FATAL,'IO','SET_RSOU','maximum forcing times NFRCLT is 99*8')
+ END IF
+!
+ WRITE(UNIT=ILUOUT,FMT='(" THERE ARE ",I2," SFC FLUX FORCINGs AT:")') 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
+ WRITE(ILUOUT,FMT='(A, I4)') "SET_RSOU/Reading Sea Surface forcing: Number=", JKT
+ READ(ILUPRE,*) ZFRCLT(JKT)%nyear, ZFRCLT(JKT)%nmonth, &
+ ZFRCLT(JKT)%nday, ZFRCLT(JKT)%xtime
+ READ(ILUPRE,*) ZSSUFL_T(JKT)
+ READ(ILUPRE,*) ZSSVFL_T(JKT)
+ READ(ILUPRE,*) ZSSTFL_T(JKT)
+ READ(ILUPRE,*) ZSSOLA_T(JKT)
+ END DO
!
-!* 2.1 STANDARD case : ZGROUND, PGROUND, TGROUND, TDGROUND
+ DO JKT = 1 , IFRCLT
+ WRITE(UNIT=ILUOUT,FMT='(F9.0, "s, date:", I3, "/", I3, "/", I5)') &
+ ZFRCLT(JKT)%xtime, ZFRCLT(JKT)%nday, &
+ ZFRCLT(JKT)%nmonth, ZFRCLT(JKT)%nyear
+ END DO
+ NINFRT= INT(ZFRCLT(2)%xtime)
+ WRITE(ILUOUT,FMT='(A)') &
+ "Number U-Stress, V-Stress, Heat turb Flux, Solar Flux Interval(s)",NINFRT
+ DO JKT = 1, IFRCLT
+ WRITE(ILUOUT,FMT='(I10,99(3F10.2))') JKT, ZSSUFL_T(JKT),ZSSVFL_T(JKT),ZSSTFL_T(JKT)
+ END DO
+ NFRCLT = IFRCLT
+ ALLOCATE(TFRCLT(NFRCLT))
+ ALLOCATE(XSSUFL_T(NFRCLT));XSSUFL_T(:)=0.
+ ALLOCATE(XSSVFL_T(NFRCLT));XSSVFL_T(:)=0.
+ ALLOCATE(XSSTFL_T(NFRCLT));XSSTFL_T(:)=0.
+ ALLOCATE(XSSOLA_T(NFRCLT));XSSOLA_T(:)=0.
+!
+ DO JKT=1,NFRCLT
+ 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)
+ END DO
+ DEALLOCATE(ZFRCLT)
+ DEALLOCATE(ZSSUFL_T)
+ DEALLOCATE(ZSSVFL_T)
+ DEALLOCATE(ZSSTFL_T)
+ DEALLOCATE(ZSSOLA_T)
+!
+!--------------------------------------------------------------------------------
+! 2.0.2 Ocean standard initialize from netcdf files
+! U,V,T,S at Z levels + Forcings at model TOP (sea surface)
+!--------------------------------------------------------------------------------
+!
+ CASE ('STANDOCE')
+!
+ READ(ILUPRE,*) ZPTOP ! P_atmosphere at sfc =P top domain
+ READ(ILUPRE,*) INFILE,INFISF
+ WRITE(ILUOUT,FMT=*) 'Netcdf files to read:',INFILE,INFISF
+ ! Open file containing initial profiles
+ CALL check(nf90_open(infile,NF90_NOWRITE,ncid), "opening NC file")
+ ! Reading dimensions and lengths
+ CALL check( nf90_inq_dimid(ncid, "depth",ndimp), "getting depth dimension id" )
+ CALL check( nf90_inquire_dimension(ncid, ndimp, len=NZ), "getting NZ " )
+ CALL check( nf90_inquire_dimension(ncid, 2, len=NLONGI), "getting NLONG " )
+ CALL check( nf90_inquire_dimension(ncid, 1, len=NLATI), "getting NLAT " )
+!
+ WRITE(ILUOUT,FMT=*) 'NB LEVLS READ NZ, NLONG NLAT ', NZ, NLONGI,NLATI
+ ALLOCATE(ZOC_TEMPERATURE(NLATI,NLONGI,NZ),ZOC_SALINITY(NLATI,NLONGI,NZ))
+ ALLOCATE(ZOC_U(NLATI,NLONGI,NZ),ZOC_V(NLATI,NLONGI,NZ))
+ ALLOCATE(ZOC_DEPTH(NZ))
+ WRITE(ILUOUT,FMT=*) 'NETCDF READING ==> Temp'
+ CALL check(nf90_inq_varid(ncid,"temperature",varid), "getting temp varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_TEMPERATURE), "reading temp")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> salinity'
+ CALL check(nf90_inq_varid(ncid,"salinity",varid), "getting salinity varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_SALINITY), "reading salinity")
+ WRITE(ILUOUT,FMT=*) 'Netcdf ==> Reading depth'
+ CALL check(nf90_inq_varid(ncid,"depth",varid), "getting depth varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_DEPTH), "reading depth")
+ WRITE(ILUOUT,FMT=*) 'depth: max min ', MAXVAL(ZOC_DEPTH),MINVAL(ZOC_DEPTH)
+ WRITE(ILUOUT,FMT=*) 'depth 1 nz: ', ZOC_DEPTH(1),ZOC_DEPTH(NZ)
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> Currents'
+ CALL check(nf90_inq_varid(ncid,"u",varid), "getting u varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_U), "reading u")
+ CALL check(nf90_inq_varid(ncid,"v",varid), "getting v varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_V), "reading v")
+ CALL check(nf90_close(ncid), "closing infile")
+ WRITE(ILUOUT,FMT=*) 'End of initial file reading'
+!
+ DO JKM=1,NZ
+ ZOC_TEMPERATURE(1,1,JKM)=ZOC_TEMPERATURE(1,1,JKM)+273.15
+ WRITE(ILUOUT,FMT=*) 'Z T(Kelvin) S(Sverdup) U V K',&
+ JKM,ZOC_DEPTH(JKM),ZOC_TEMPERATURE(1,1,JKM),ZOC_SALINITY(1,1,JKM),ZOC_U(1,1,JKM),ZOC_V(1,1,JKM), JKM
+ ENDDO
+ ! number of data levels
+ ILEVELM=NZ
+ ! Model bottom
+ ZTGROUND = ZOC_TEMPERATURE(1,1,ILEVELM)
+ ZMRGROUND = ZOC_SALINITY(1,1,ILEVELM)
+ ZZGROUND=0.
+ ! Allocate required memory
+ ALLOCATE(ZHEIGHTM(ILEVELM))
+ ALLOCATE(ZT(ILEVELM))
+ ALLOCATE(ZTV(ILEVELM))
+ ALLOCATE(ZMR(ILEVELM))
+ ALLOCATE(ZTHV(ILEVELM))
+ ALLOCATE(ZTHL(ILEVELM))
+ ALLOCATE(ZRT(ILEVELM))
+ ! Going from the inverse model grid (data) to the normal one
+ DO JKM= 1,ILEVELM
+ ! Z axis reoriented as in the model
+ IDX = ILEVELM-JKM+1
+ ZT(JKM) = ZOC_TEMPERATURE(1,1,IDX)
+ ZMR(JKM) = ZOC_SALINITY(1,1,IDX)
+ ! 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)
+ ! translation/inversion
+ ZHEIGHTM(JKM) = -ZOC_DEPTH(IDX) + ZOC_DEPTH(ILEVELM)
+ WRITE(ILUOUT,FMT=*) 'End gridmodel comput: JKM IDX depth T S ZHEIGHTM', &
+ JKM,IDX,ZOC_DEPTH(IDX),ZT(JKM),ZMR(JKM),ZHEIGHTM(JKM)
+ END DO
+ ! complete ther variables
+ ZTV = ZT
+ ZTHV = ZT
+ ZRT = ZMR
+ ZTHL = ZT
+ ZTH = ZT
+ ! INIT --- U V -----
+ ILEVELU=nz ! Same nb of levels for u,v,T,S
+ !Assume that current and temp are given at same level
+ ALLOCATE(ZHEIGHTU(ILEVELU))
+ ALLOCATE(ZU(ILEVELU),ZV(ILEVELU))
+ ZHEIGHTU=ZHEIGHTM
+ DO JKM= 1,ILEVELU
+ ! Z axis reoriented as in the model
+ IDX = ILEVELU-JKM+1
+ ZU(JKM) = ZOC_U(1,1,IDX)
+ ZV(JKM) = ZOC_V(1,1,IDX)
+ ! 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)
+ END DO
+!
+ DEALLOCATE(ZOC_TEMPERATURE)
+ DEALLOCATE(ZOC_SALINITY)
+ DEALLOCATE(ZOC_U)
+ DEALLOCATE(ZOC_V)
+ DEALLOCATE(ZOC_DEPTH)
+!
+ ! Reading/initializing surface forcings
+!
+ WRITE(ILUOUT,FMT=*) 'netcdf sfc forcings file to be read:',infisf
+ ! Open of sfc forcing file
+ CALL check(nf90_open(infisf,NF90_NOWRITE,ncid), "opening NC file")
+ ! Reading dimension and length
+ CALL check( nf90_inq_dimid(ncid,"t",ndimp), "getting time dimension id" )
+ CALL check( nf90_inquire_dimension(ncid, ndimp, len=ntcindy), "getting ntcindy " )
+!
+ WRITE(ILUOUT,FMT=*) 'nb sfc-forcing time ntcindy=',ntcindy
+ ALLOCATE(ZOC_LE(ntcindy))
+ ALLOCATE(ZOC_H(ntcindy))
+ ALLOCATE(ZOC_SW_DOWN(ntcindy))
+ ALLOCATE(ZOC_SW_UP(ntcindy))
+ ALLOCATE(ZOC_LW_DOWN(ntcindy))
+ ALLOCATE(ZOC_LW_UP(ntcindy))
+ ALLOCATE(ZOC_TAUX(ntcindy))
+ ALLOCATE(ZOC_TAUY(ntcindy))
+!
+ WRITE(ILUOUT,FMT=*)'Netcdf Reading ==> LE'
+ CALL check(nf90_inq_varid(ncid,"LE",varid), "getting LE varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_LE), "reading LE flux")
+ WRITE(ILUOUT,FMT=*)'Netcdf Reading ==> H'
+ CALL check(nf90_inq_varid(ncid,"H",varid), "getting H varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_H), "reading H flux")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> SW_DOWN'
+ CALL check(nf90_inq_varid(ncid,"SW_DOWN",varid), "getting SW_DOWN varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_SW_DOWN), "reading SW_DOWN")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> SW_UP'
+ CALL check(nf90_inq_varid(ncid,"SW_UP",varid), "getting SW_UP varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_SW_UP), "reading SW_UP")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> LW_DOWN'
+ CALL check(nf90_inq_varid(ncid,"LW_DOWN",varid), "getting LW_DOWN varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_LW_DOWN), "reading LW_DOWN")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> LW_UP'
+ CALL check(nf90_inq_varid(ncid,"LW_UP",varid), "getting LW_UP varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_LW_UP), "reading LW_UP")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> TAUX'
+ CALL check(nf90_inq_varid(ncid,"TAUX",varid), "getting TAUX varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_TAUX), "reading TAUX")
+ WRITE(ILUOUT,FMT=*) 'Netcdf Reading ==> TAUY'
+ CALL check(nf90_inq_varid(ncid,"TAUY",varid), "getting TAUY varid")
+ CALL check(nf90_get_var(ncid,varid,ZOC_TAUY), "reading TAUY")
+ CALL check(nf90_close(ncid), "closing infisfS")
+!
+ WRITE(ILUOUT,FMT=*) ' Forcing-Number LE H SW_down SW_up LW_down LW_up TauX TauY'
+ DO JKM=1,NTCINDY
+ 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=NTCINDY
+ 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
+ ! 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))
+ ! 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)
+ ENDDO
+ ! Allocate and Writing the corresponding variables in module MODD_OCEAN_FRC
+ NFRCLT=IFRCLT
+ ! value to read later on file ?
+ NINFRT=600
+ ALLOCATE(TFRCLT(NFRCLT))
+ ALLOCATE(XSSUFL_T(NFRCLT));XSSUFL_T(:)=0.
+ ALLOCATE(XSSVFL_T(NFRCLT));XSSVFL_T(:)=0.
+ 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
+ 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)
+ END DO
+ DEALLOCATE(ZFRCLT)
+ DEALLOCATE(ZSSUFL_T)
+ DEALLOCATE(ZSSVFL_T)
+ DEALLOCATE(ZSSTFL_T)
+ DEALLOCATE(ZSSOLA_T)
+ DEALLOCATE(ZOC_LE)
+ DEALLOCATE(ZOC_H)
+ DEALLOCATE(ZOC_SW_DOWN)
+ DEALLOCATE(ZOC_SW_UP)
+ DEALLOCATE(ZOC_LW_DOWN)
+ DEALLOCATE(ZOC_LW_UP)
+ DEALLOCATE(ZOC_TAUX)
+ DEALLOCATE(ZOC_TAUY)
+ ! END OCEAN STANDARD
+!
+!
+!* 2.1 ATMOSPHERIC STANDARD case : ZGROUND, PGROUND, TGROUND, TDGROUND
! (Pressure, dd, ff) ,
! (Pressure, T, Td)
!
@@ -448,8 +811,7 @@ SELECT CASE(YKIND)
ALLOCATE(ZHEIGHTM(ILEVELM)) ! Allocate memory for needed
ALLOCATE(ZTHV(ILEVELM)) ! arrays
ALLOCATE(ZMR(ILEVELM))
- ALLOCATE(ZTV(ILEVELM)) ! Allocate memory for intermediate
- ! arrays
+ ALLOCATE(ZTV(ILEVELM)) ! Allocate memory for intermediate arrays
ALLOCATE(ZTHL(ILEVELM))
ALLOCATE(ZRT(ILEVELM))
!
@@ -1203,20 +1565,36 @@ END DO
ALLOCATE(ZEXNFLUX(IKU))
ALLOCATE(ZEXNMASS(IKU))
ALLOCATE(ZPRESS(IKU))
+ALLOCATE(ZPREFLUX(IKU))
ALLOCATE(ZFRAC_ICE(IKU))
ALLOCATE(ZRSATW(IKU))
ALLOCATE(ZRSATI(IKU))
ALLOCATE(ZMRT(IKU))
ZMRT=ZMRM+ZMRCM+ZMRIM
-ZEXNSURF=(ZPGROUND/XP00)**(XRD/XCPD)
ZTHVM=ZTHLM
-DO JLOOP=1,20 ! loop for pression
- CALL COMPUTE_EXNER_FROM_GROUND(ZTHVM,ZZMASS_PROFILE(:),ZEXNSURF,ZEXNFLUX,ZEXNMASS)
- ZPRESS(:)=XP00*(ZEXNMASS(:))**(XCPD/XRD)
- CALL TH_R_FROM_THL_RT_1D('T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
- ZRSATW, ZRSATI)
- ZTHVM(:)=ZTHM(:)*(1.+XRV/XRD*ZMRM(:))/(1.+(ZMRM(:)+ZMRIM(:)+ZMRCM(:)))
-ENDDO
+!
+IF (LOCEAN) THEN
+ ZRHODM(:)=XRH00OCEAN*(1.-XALPHAOC*(ZTHLM(:) - XTH00OCEAN)&
+ +XBETAOC* (ZMRM(:) - XSA00OCEAN))
+ ZPREFLUX(IKU)=ZPTOP
+ DO JK=IKU-1,2,-1
+ ZPREFLUX(JK) = ZPREFLUX(JK+1) + XG*ZRHODM(JK)*(ZZFLUX_PROFILE(JK+1)-ZZFLUX_PROFILE(JK))
+ END DO
+ ZPGROUND=ZPREFLUX(2)
+ WRITE(ILUOUT,FMT=*)'ZPGROUND i.e. Pressure at ocean domain bottom',ZPGROUND
+ ZTHM=ZTHVM
+ELSE
+! Atmospheric case
+ ZEXNSURF=(ZPGROUND/XP00)**(XRD/XCPD)
+ DO JLOOP=1,20 ! loop for pression
+ CALL COMPUTE_EXNER_FROM_GROUND(ZTHVM,ZZMASS_PROFILE(:),ZEXNSURF,ZEXNFLUX,ZEXNMASS)
+ ZPRESS(:)=XP00*(ZEXNMASS(:))**(XCPD/XRD)
+ CALL TH_R_FROM_THL_RT_1D('T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
+ ZRSATW, ZRSATI)
+ ZTHVM(:)=ZTHM(:)*(1.+XRV/XRD*ZMRM(:))/(1.+(ZMRM(:)+ZMRIM(:)+ZMRCM(:)))
+ ENDDO
+ENDIF
+!
DEALLOCATE(ZEXNFLUX)
DEALLOCATE(ZEXNMASS)
DEALLOCATE(ZPRESS)
@@ -1224,7 +1602,6 @@ DEALLOCATE(ZFRAC_ICE)
DEALLOCATE(ZRSATW)
DEALLOCATE(ZRSATI)
DEALLOCATE(ZMRT)
-
!-------------------------------------------------------------------------------
!
!* 4. COMPUTE FIELDS ON THE MODEL GRID (WITH OROGRAPHY)
@@ -1234,6 +1611,20 @@ CALL SET_MASS(TPFILE,GPROFILE_IN_PROC, ZZFLUX_PROFILE, &
ZTHVM,ZMRM,ZUW,ZVW,OSHIFT,OBOUSS,PJ,HFUNU,HFUNV, &
PMRCM=ZMRCM,PMRIM=ZMRIM,PCORIOZ=PCORIOZ)
!
+DEALLOCATE(ZPREFLUX)
+DEALLOCATE(ZHEIGHTM)
+DEALLOCATE(ZTHV)
+DEALLOCATE(ZMR)
+DEALLOCATE(ZTHL)
!-------------------------------------------------------------------------------
+CONTAINS
+ SUBROUTINE CHECK(STATUS,LOC)
+ INTEGER, INTENT(IN) :: STATUS
+ CHARACTER(LEN=*), INTENT(IN) :: LOC
+ IF(STATUS /= NF90_NOERR) THEN
+ WRITE(ILUOUT,FMT=*) "Error at ", LOC
+ WRITE(ILUOUT,FMT=*) NF90_STRERROR(STATUS)
+ END IF
+ END SUBROUTINE check
!
END SUBROUTINE SET_RSOU
diff --git a/src/MNH/spawn_model2.f90 b/src/MNH/spawn_model2.f90
index 82c8f1ecac1eb09d844ef9c48dd3e30a7b0e3d5e..81baecd56f380ce62aca8e065907d4a6b5a6d9a7 100644
--- a/src/MNH/spawn_model2.f90
+++ b/src/MNH/spawn_model2.f90
@@ -1525,6 +1525,7 @@ IF (NVERB >= 5) THEN
WRITE(ILUOUT,*) 'SPAWN_MODEL2: NVERB=',NVERB
WRITE(ILUOUT,*) 'SPAWN_MODEL2: XLON0,XLAT0,XBETA=',XLON0,XLAT0,XBETA
WRITE(ILUOUT,*) 'SPAWN_MODEL2: LCARTESIAN=',LCARTESIAN
+ WRITE(ILUOUT,*) 'SPAWN_MODEL2: LOCEAN,LCOUPLES=',LOCEAN,LCOUPLES
IF(LCARTESIAN) THEN
WRITE(ILUOUT,*) 'SPAWN_MODEL2: No map projection used.'
ELSE
diff --git a/src/MNH/th_r_from_thl_rt_1d.f90 b/src/MNH/th_r_from_thl_rt_1d.f90
index fcec2372fd5ff3a140c04b04df94a31073645167..6149852c13b350153f94245a4a6bcb5b4df969d4 100644
--- a/src/MNH/th_r_from_thl_rt_1d.f90
+++ b/src/MNH/th_r_from_thl_rt_1d.f90
@@ -1,3 +1,7 @@
+!MNH_LIC Copyright 1994-2021 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.
! ######spl
MODULE MODI_TH_R_FROM_THL_RT_1D
! ###############################
@@ -67,7 +71,8 @@ REAL, DIMENSION(:), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturati
! ------------
!
USE MODI_COMPUTE_FRAC_ICE
-USE MODD_CST !, ONLY: XP00, XRD, XCPD, XCPV, XCL, XCI, XLVTT, XTT, XLSTT
+USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODE_THERMO
!
IMPLICIT NONE
@@ -131,8 +136,11 @@ PTH(:)=PTHL(:)+ZLVOCPEXN(:)*PRL(:)+ZLSOCPEXN(:)*PRI(:)
! ---------
DO II=1,JITER
- ZT(:)=PTH(:)*ZEXN(:)
-
+ IF (LOCEAN) THEN
+ ZT=PTH
+ ELSE
+ ZT(:)=PTH(:)*ZEXN(:)
+ END IF
!Computation of liquid/ice fractions
PFRAC_ICE(:) = 0.
WHERE(PRL(:)+PRI(:) > 1.E-20)
diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90
index cc09425a87ef4a81be15d4d2edb65a21a5492ed8..0074628147dd16a5e56f3fff71188cbffd4ee656 100644
--- a/src/MNH/turb.f90
+++ b/src/MNH/turb.f90
@@ -345,6 +345,7 @@ END MODULE MODI_TURB
! P. Wautelet + Benoit Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets
! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct
! R. Honnert/V. Masson 02/2021: new mixing length in the grey zone
+!! J.L. Redelsperger 03/2021: add Ocean LES case
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -358,12 +359,14 @@ use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, lbud
USE MODD_CONF
USE MODD_CST
USE MODD_CTURB
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_field, only: tfielddata, TYPEREAL
USE MODD_IO, ONLY: TFILEDATA
USE MODD_LES
USE MODD_NSV
USE MODD_PARAMETERS, ONLY: JPVEXT_TURB
USE MODD_PARAM_LIMA
+USE MODD_REF
USE MODD_TURB_n, ONLY: XCADAP
!
USE MODI_GRADIENT_M
@@ -643,7 +646,11 @@ END DO
!
!* 2.2 Exner function at t
!
-ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD)
+IF (LOCEAN) THEN
+ ZEXN(:,:,:) = 1.
+ELSE
+ ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD)
+END IF
!
!* 2.3 dissipative heating coeff a t
!
@@ -1512,15 +1519,26 @@ IF (.NOT. ORMC01) THEN
!
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
- DO JK=IKTB,IKTE
- ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
- -PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
- IF ( PLM(JI,JJ,JK)>ZD) THEN
- PLM(JI,JJ,JK)=ZD
- ELSE
- EXIT
- ENDIF
- END DO
+ IF (LOCEAN) THEN
+ DO JK=IKTE,IKTB,-1
+ ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
+ IF ( PLM(JI,JJ,JK)>ZD) THEN
+ PLM(JI,JJ,JK)=ZD
+ ELSE
+ EXIT
+ ENDIF
+ END DO
+ ELSE
+ DO JK=IKTB,IKTE
+ ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
+ -PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
+ IF ( PLM(JI,JJ,JK)>ZD) THEN
+ PLM(JI,JJ,JK)=ZD
+ ELSE
+ EXIT
+ ENDIF
+ END DO
+ ENDIF
END DO
END DO
END IF
@@ -1587,7 +1605,6 @@ END IF
ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
!
-! For dry simulations
IF (KRR>0) THEN
DO JK = IKTB+1,IKTE-1
DO JJ=1,SIZE(PUT,2)
@@ -1596,8 +1613,12 @@ IF (KRR>0) THEN
(PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK ))
ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK ,1))/PDZZ(JI,JJ,JK+KKL)+ &
(PRT(JI,JJ,JK ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK ))
- ZVAR=XG/PTHVREF(JI,JJ,JK)* &
+ IF (LOCEAN) THEN
+ ZVAR=XG*(XALPHAOC*ZDTHLDZ(JI,JJ,JK)-XBETAOC*ZDRTDZ(JI,JJ,JK))
+ ELSE
+ ZVAR=XG/PTHVREF(JI,JJ,JK)* &
(ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK))
+ END IF
!
IF (ZVAR>0.) THEN
PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
@@ -1606,14 +1627,18 @@ IF (KRR>0) THEN
END DO
END DO
END DO
-ELSE
+ELSE! For dry atmos or unsalted ocean runs
DO JK = IKTB+1,IKTE-1
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK ))/PDZZ(JI,JJ,JK+KKL)+ &
(PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK ))
- ZVAR=XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
- !
+ IF (LOCEAN) THEN
+ ZVAR= XG*XALPHAOC*ZDTHLDZ(JI,JJ,JK)
+ ELSE
+ ZVAR= XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
+ END IF
+!
IF (ZVAR>0.) THEN
PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
@@ -1631,8 +1656,12 @@ ELSE
ZDRTDZ(:,:,IKB)=0
ENDIF
!
-ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)* &
- (ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
+IF (LOCEAN) THEN
+ ZWORK2D(:,:)=XG*(XALPHAOC*ZDTHLDZ(:,:,IKB)-XBETAOC*ZDRTDZ(:,:,IKB))
+ELSE
+ ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)* &
+ (ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
+END IF
WHERE(ZWORK2D(:,:)>0.)
PLM(:,:,IKB)=MAX(XMNH_EPSILON,MIN( PLM(:,:,IKB), &
0.76* SQRT(PTKET(:,:,IKB)/ZWORK2D(:,:))))
@@ -1645,15 +1674,26 @@ IF (.NOT. ORMC01) THEN
!
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
- DO JK=IKTB,IKTE
- ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
- *PDIRCOSZW(JI,JJ)
- IF ( PLM(JI,JJ,JK)>ZD) THEN
- PLM(JI,JJ,JK)=ZD
- ELSE
- EXIT
- ENDIF
- END DO
+ IF (LOCEAN) THEN
+ DO JK=IKTE,IKTB,-1
+ ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
+ IF ( PLM(JI,JJ,JK)>ZD) THEN
+ PLM(JI,JJ,JK)=ZD
+ ELSE
+ EXIT
+ ENDIF
+ END DO
+ ELSE
+ DO JK=IKTB,IKTE
+ ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
+ *PDIRCOSZW(JI,JJ)
+ IF ( PLM(JI,JJ,JK)>ZD) THEN
+ PLM(JI,JJ,JK)=ZD
+ ELSE
+ EXIT
+ ENDIF
+ END DO
+ ENDIF
END DO
END DO
END IF
diff --git a/src/MNH/turb_ver.f90 b/src/MNH/turb_ver.f90
index fce78c562131d968b5d3d307a70f25d06c226716..ddd97f1470d3cccd5ce56dae4b2006b823a18b97 100644
--- a/src/MNH/turb_ver.f90
+++ b/src/MNH/turb_ver.f90
@@ -311,6 +311,7 @@ END MODULE MODI_TURB_VER
!! 10/2012 (J.Escobar) Bypass PGI bug , redefine some allocatable array inplace of automatic
!! 08/2014 (J.Escobar) Bypass PGI memory leak bug , replace IF statement with IF THEN ENDIF
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! JL Redelsperger 03/2021 : add Ocean LES case
!!--------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -318,6 +319,7 @@ END MODULE MODI_TURB_VER
!
USE MODD_CST
USE MODD_CTURB
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_field, only: tfielddata, TYPEREAL
USE MODD_IO, ONLY: TFILEDATA
USE MODD_PARAMETERS
@@ -463,9 +465,9 @@ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: &
!!$REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3),NSV) :: &
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: &
ZPSI_SV, & ! Prandtl number for scalars
- ZREDS1, & ! 1D Redeslperger number R_sv
- ZRED2THS, & ! 3D Redeslperger number R*2_thsv
- ZRED2RS ! 3D Redeslperger number R*2_rsv
+ ZREDS1, & ! 1D Redelsperger number R_sv
+ ZRED2THS, & ! 3D Redelsperger number R*2_thsv
+ ZRED2RS ! 3D Redelsperger number R*2_rsv
!
LOGICAL :: GUSERV ! flag to use water vapor
INTEGER :: IKB,IKE ! index value for the Beginning
@@ -509,7 +511,6 @@ ALLOCATE ( &
!
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
-
!
!
! 3D Redelsperger numbers
@@ -529,7 +530,11 @@ CALL PRANDTL(KKA,KKU,KKL,KRR,KRRI,OTURB_FLX, &
!
! Buoyancy coefficient
!
-ZBETA = XG/PTHVREF
+IF (LOCEAN) THEN
+ ZBETA = XG*XALPHAOC
+ELSE
+ ZBETA = XG/PTHVREF
+END IF
!
! Square root of Tke
!
diff --git a/src/MNH/turb_ver_dyn_flux.f90 b/src/MNH/turb_ver_dyn_flux.f90
index 10a40478b7dd7ba01fc2191d1646c78631448743..7d5983a0408da4bd324bb689a00a25db9ee6c2e7 100644
--- a/src/MNH/turb_ver_dyn_flux.f90
+++ b/src/MNH/turb_ver_dyn_flux.f90
@@ -280,6 +280,7 @@ END MODULE MODI_TURB_VER_DYN_FLUX
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! Q. Rodier 17/01/2019 : cleaning : remove cyclic conditions on DP and ZA
+!! JL Redelsperger 03/2021 : Add Ocean & O-A Autocoupling LES Cases
!!--------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -288,11 +289,15 @@ END MODULE MODI_TURB_VER_DYN_FLUX
USE MODD_CONF
USE MODD_CST
USE MODD_CTURB
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_field, only: tfielddata, TYPEREAL
USE MODD_IO, ONLY: TFILEDATA
USE MODD_LES
USE MODD_NSV
+USE MODD_OCEANH
USE MODD_PARAMETERS
+USE MODD_REF, ONLY : LCOUPLES
+USE MODD_TURB_n
!
!
USE MODI_GRADIENT_U
@@ -459,25 +464,48 @@ ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) &
! average this flux to be located at the U,W vorticity point
ZCOEFS(:,:,1:1)=MXM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
!
-! compute the explicit tangential flux at the W point
-ZSOURCE(:,:,IKB) = &
- PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
- -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
!
-! add the vertical part or the surface flux at the U,W vorticity point
-
-ZSOURCE(:,:,IKB:IKB) = &
- ( MXM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
- + MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+! ZSOURCE= FLUX /DZ
+IF (LOCEAN) THEN ! OCEAN MODEL ONLY
+ ! Sfx flux assumed to be in SI & at vorticity point
+ IF (LCOUPLES) THEN
+ ZSOURCE(:,:,IKE:IKE) = XSSUFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
+ *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)))
+ ELSE
+ ZSOURCE(:,:,IKE) = XSSUFL(:,:)
+ ZSOURCE(:,:,IKE:IKE) = ZSOURCE (:,:,IKE:IKE) /PDZZ(:,:,IKE:IKE) &
+ *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
+ ENDIF
+ !No flux at the ocean domain bottom
+ ZSOURCE(:,:,IKB) = 0.
+ ZSOURCE(:,:,IKTB+1:IKTE-1) = 0
+!
+ELSE !ATMOS MODEL ONLY
+ IF (LCOUPLES) THEN
+ ZSOURCE(:,:,IKB:IKB) = XSSUFL_C(:,:,1:1)/PDZZ(:,:,IKB:IKB) &
+ * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
+ ELSE
+ ! compute the explicit tangential flux at the W point
+ ZSOURCE(:,:,IKB) = &
+ PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+ -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
+ -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+!
+ ! add the vertical part or the surface flux at the U,W vorticity point
+!
+ ZSOURCE(:,:,IKB:IKB) = &
+ ( MXM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
+ + MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
*ZUSLOPEM(:,:,1:1) &
-ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
*ZVSLOPEM(:,:,1:1) ) &
- - ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL &
- ) * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
+ - ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL &
+ ) * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
+ ENDIF
!
-ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
-ZSOURCE(:,:,IKE) = 0.
+ ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
+ ZSOURCE(:,:,IKE) = 0.
+ENDIF !end ocean or atmosphere cases
!
! Obtention of the split U at t+ deltat
!
@@ -504,6 +532,12 @@ ZFLXZ(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB)) * &
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+IF (LOCEAN) THEN !ocean model at phys sfc (ocean domain top)
+ ZFLXZ(:,:,IKE:IKE) = MXM(PDZZ(:,:,IKE:IKE)) * &
+ ZSOURCE(:,:,IKE:IKE) &
+ / 0.5 / ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+END IF
!
IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
! stores the U wind component vertical flux
@@ -533,7 +567,15 @@ PDP(:,:,:) = - MZF( MXF ( ZFLXZ * GZ_U_UW(PUM,PDZZ) ) )
PDP(:,:,IKB:IKB) = - MXF ( &
ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PUM(:,:,IKB+KKL:IKB+KKL)-PUM(:,:,IKB:IKB)) &
/ MXM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
- )
+ )
+!
+IF (LOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+ PDP(:,:,IKE:IKE) = - MXF ( &
+ ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-KKL:IKE-KKL)) &
+ / MXM(PDZZ(:,:,IKE-KKL:IKE-KKL)) &
+ )
+END IF
!
! Storage in the LES configuration
!
@@ -552,8 +594,12 @@ END IF
!
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
- ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
! used to compute the W source at the ground
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ IF (LOCEAN) THEN
+ ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
+ END IF
+
!
IF (.NOT. LFLAT) THEN
PRWS(:,:,:)= PRWS &
@@ -583,6 +629,21 @@ IF(HTURBDIM=='3DIM') THEN
) / (0.5*(PDXX(:,:,IKB+KKL:IKB+KKL)+PDXX(:,:,IKB:IKB))) &
)
!
+IF (LOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+ ZA(:,:,IKE:IKE) = - MXF ( &
+ ZFLXZ(:,:,IKE-KKL:IKE-KKL) * &
+ ( DXM( PWM(:,:,IKE-KKL:IKE-KKL) ) &
+ -MXM( (PWM(:,:,IKE-2*KKL:IKE-2*KKL )-PWM(:,:,IKE-KKL:IKE-KKL)) &
+ /(PDZZ(:,:,IKE-2*KKL:IKE-2*KKL)+PDZZ(:,:,IKE-KKL:IKE-KKL)) &
+ +(PWM(:,:,IKE-KKL:IKE-KKL)-PWM(:,:,IKE:IKE )) &
+ /(PDZZ(:,:,IKE-KKL:IKE-KKL)+PDZZ(:,:,IKE:IKE )) &
+ ) &
+ * PDZX(:,:,IKE-KKL:IKE-KKL) &
+ ) / (0.5*(PDXX(:,:,IKE-KKL:IKE-KKL)+PDXX(:,:,IKE:IKE))) &
+ )
+END IF
+ !
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
! Storage in the LES configuration
@@ -636,24 +697,44 @@ ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) &
! average this flux to be located at the V,W vorticity point
ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
!
-! compute the explicit tangential flux at the W point
-ZSOURCE(:,:,IKB) = &
- PTAU11M(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
- +PTAU12M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
-!
-! add the vertical part or the surface flux at the V,W vorticity point
-ZSOURCE(:,:,IKB:IKB) = &
- ( MYM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
- + MYM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZUSLOPEM(:,:,1:1) &
- +ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
- *ZVSLOPEM(:,:,1:1) ) &
- - ZCOEFS(:,:,1:1) * PVM(:,:,IKB:IKB) * PIMPL &
- ) * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
-!
+IF (LOCEAN) THEN ! Ocean case
+ IF (LCOUPLES) THEN
+ ZSOURCE(:,:,IKE:IKE) = XSSVFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
+ *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+ ELSE
+ ZSOURCE(:,:,IKE) = XSSVFL(:,:)
+ ZSOURCE(:,:,IKE:IKE) = ZSOURCE(:,:,IKE:IKE)/PDZZ(:,:,IKE:IKE) &
+ *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+ END IF
+ !No flux at the ocean domain bottom
+ ZSOURCE(:,:,IKB) = 0.
+ELSE ! Atmos case
+ IF (.NOT.LCOUPLES) THEN ! only atmosp sans couplage
+ ! compute the explicit tangential flux at the W point
+ ZSOURCE(:,:,IKB) = &
+ PTAU11M(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+ +PTAU12M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) &
+ -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+!
+ ! add the vertical part or the surface flux at the V,W vorticity point
+ ZSOURCE(:,:,IKB:IKB) = &
+ ( MYM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
+ + MYM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZUSLOPEM(:,:,1:1) &
+ +ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZVSLOPEM(:,:,1:1) ) &
+ - ZCOEFS(:,:,1:1) * PVM(:,:,IKB:IKB) * PIMPL &
+ ) * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
+!
+ ELSE !atmosphere quand couplage
+ ! flux en input supposé etre en SI et en point de voticité
+ ZSOURCE(:,:,IKB:IKB) = -XSSVFL_C(:,:,1:1)/(1.*PDZZ(:,:,IKB:IKB)) &
+ * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
+ ENDIF
+ !No flux at the atmosphere top
+ ZSOURCE(:,:,IKE) = 0.
+ENDIF ! End of Ocean or Atmospher Cases
ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
-ZSOURCE(:,:,IKE) = 0.
!
! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(KKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
@@ -679,6 +760,13 @@ ZFLXZ(:,:,IKB:IKB) = MYM(PDZZ(:,:,IKB:IKB)) * &
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
+IF (LOCEAN) THEN
+ ZFLXZ(:,:,IKE:IKE) = MYM(PDZZ(:,:,IKE:IKE)) * &
+ ZSOURCE(:,:,IKE:IKE) &
+ / 0.5 / ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+END IF
+!
IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
! stores the V wind component vertical flux
TZFIELD%CMNHNAME = 'VW_VFLX'
@@ -710,6 +798,14 @@ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PVM(:,:,IKB+KKL:IKB+KKL)-PVM(:,:,IKB:IKB)) &
/ MYM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
)
!
+IF (LOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+ ZA(:,:,IKE:IKE) = - MYF ( &
+ ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-KKL:IKE-KKL)) &
+ / MYM(PDZZ(:,:,IKE-KKL:IKE-KKL)) &
+ )
+END IF
+!
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
! Storage in the LES configuration
@@ -729,6 +825,9 @@ END IF
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ IF (LOCEAN) THEN
+ ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
+ END IF
!
IF (.NOT. L2D) THEN
IF (.NOT. LFLAT) THEN
@@ -759,6 +858,20 @@ IF(HTURBDIM=='3DIM') THEN
) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB))) &
)
!
+ IF (LOCEAN) THEN
+ ZA(:,:,IKE:IKE) = - MYF ( &
+ ZFLXZ(:,:,IKE-KKL:IKE-KKL) * &
+ ( DYM( PWM(:,:,IKE-KKL:IKE-KKL) ) &
+ -MYM( (PWM(:,:,IKE-2*KKL:IKE-2*KKL)-PWM(:,:,IKE-KKL:IKE-KKL)) &
+ /(PDZZ(:,:,IKE-2*KKL:IKE-2*KKL)+PDZZ(:,:,IKE-KKL:IKE-KKL)) &
+ +(PWM(:,:,IKE-KKL:IKE-KKL)-PWM(:,:,IKE:IKE )) &
+ /(PDZZ(:,:,IKE-KKL:IKE-KKL)+PDZZ(:,:,IKE:IKE )) &
+ ) &
+ * PDZY(:,:,IKE-KKL:IKE-KKL) &
+ ) / (0.5*(PDYY(:,:,IKE-KKL:IKE-KKL)+PDYY(:,:,IKE:IKE))) &
+ )
+ END IF
+!
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
END IF
diff --git a/src/MNH/turb_ver_thermo_flux.f90 b/src/MNH/turb_ver_thermo_flux.f90
index 22fdeae7c5a21016639ef91217ece653545fba23..c9187c7190a8d8b72c779a14015db3a4b195c35b 100644
--- a/src/MNH/turb_ver_thermo_flux.f90
+++ b/src/MNH/turb_ver_thermo_flux.f90
@@ -323,11 +323,13 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
!! 2012-02 (Y. Seity) add possibility to run with reversed
!! vertical levels
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!! 2018 (D. Ricard) last version of HGRAD turbulence scheme
+!! 2021 (D. Ricard) last version of HGRAD turbulence scheme
!! Leronard terms instead of Reynolds terms
!! applied to vertical fluxes of r_np and Thl
!! for implicit version of turbulence scheme
!! corrections and cleaning
+!! JL Redelsperger : 03/2021: Ocean and Autocoupling O-A LES Cases
+!! Sfc flux shape for LDEEPOC Case
!!--------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -336,13 +338,19 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
USE MODD_CST
USE MODD_CTURB
use modd_field, only: tfielddata, TYPEREAL
-USE MODD_GRID_n, ONLY: XZS, XXHAT
+USE MODD_GRID_n, ONLY: XZS, XXHAT, XYHAT
USE MODD_IO, ONLY: TFILEDATA
USE MODD_METRICS_n
USE MODD_PARAMETERS
USE MODD_TURB_n, ONLY: LHGRAD, XCOEFHGRADTHL, XCOEFHGRADRM, XALTHGRAD, XCLDTHOLD
USE MODD_CONF
USE MODD_LES
+USE MODD_DIM_n
+USE MODD_DYN_n, ONLY : LOCEAN
+USE MODD_OCEANH
+USE MODD_REF, ONLY : LCOUPLES
+USE MODD_TURB_n
+USE MODD_FRC
!
USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
@@ -362,6 +370,8 @@ USE MODE_IO_FIELD_WRITE, only: IO_Field_write
USE MODE_PRANDTL
!
USE MODI_SECOND_MNH
+USE MODE_ll
+USE MODE_GATHER_ll
!
IMPLICIT NONE
!
@@ -475,10 +485,35 @@ INTEGER :: IKB,IKE ! I index values for the Beginning and End
! mass points of the domain in the 3 direct.
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: JI, JJ ! loop indexes
+!
+!
+INTEGER :: IIB,IJB ! Lower bounds of the physical
+ ! sub-domain in x and y directions
+INTEGER :: IIE,IJE ! Upper bounds of the physical
+ ! sub-domain in x and y directions
+!
+REAL, DIMENSION(:), ALLOCATABLE :: ZXHAT_ll ! Position x in the conformal
+ ! plane (array on the complete domain)
+REAL, DIMENSION(:), ALLOCATABLE :: ZYHAT_ll ! Position y in the conformal
+ ! plane (array on the complete domain)
+!
+!
+CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
+CHARACTER (LEN=LEN_HREC) :: YRECFM ! Name of the desired field in LFIFM file
!
REAL :: ZTIME1, ZTIME2
REAL :: ZDELTAX
-!
+REAL :: ZXBEG,ZXEND,ZYBEG,ZYEND ! Forcing size for ocean deep convection
+REAL, DIMENSION(SIZE(XXHAT),SIZE(XYHAT)) :: ZDIST ! distance
+ ! from the center of the cooling
+REAL :: ZFLPROV
+INTEGER :: JKM ! vertical index loop
+INTEGER :: JSW
+REAL :: ZSWA ! index for time flux interpolation
+!
+INTEGER :: IIU, IJU
+INTEGER :: IRESP
INTEGER :: JK
LOGICAL :: GUSERV ! flag to use water
LOGICAL :: GFTH2 ! flag to use w'th'2
@@ -491,6 +526,36 @@ TYPE(TFIELDDATA) :: TZFIELD
!
!* 1. PRELIMINARIES
! -------------
+! Size for a given proc & a given model
+IIU=SIZE(PTHLM,1)
+IJU=SIZE(PTHLM,2)
+!
+!! Compute Shape of sfc flux for Oceanic Deep Conv Case
+!
+IF (LOCEAN .AND. LDEEPOC) THEN
+ !* COMPUTES THE PHYSICAL SUBDOMAIN BOUNDS
+ ALLOCATE(ZXHAT_ll(NIMAX_ll+2*JPHEXT),ZYHAT_ll(NJMAX_ll+2*JPHEXT))
+ !compute ZXHAT_ll = position in the (0:Lx) domain 1 (Lx=Size of domain1 )
+ !compute XXHAT_ll = position in the (L0_subproc,Lx_subproc) domain for the current subproc
+ ! L0_subproc as referenced in the full domain 1
+ CALL GATHERALL_FIELD_ll('XX',XXHAT,ZXHAT_ll,IRESP)
+ CALL GATHERALL_FIELD_ll('YY',XYHAT,ZYHAT_ll,IRESP)
+ CALL GET_DIM_EXT_ll('B',IIU,IJU)
+ CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+ DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ ZDIST(JI,JJ) = SQRT( &
+ (( (XXHAT(JI)+XXHAT(JI+1))*0.5 - XCENTX_OC ) / XRADX_OC)**2 + &
+ (( (XYHAT(JJ)+XYHAT(JJ+1))*0.5 - XCENTY_OC ) / XRADY_OC)**2 &
+ )
+ END DO
+ END DO
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IF ( ZDIST(JI,JJ) > 1.) XSSTFL(JI,JJ)=0.
+ END DO
+ END DO
+END IF !END DEEP OCEAN CONV CASE
!
IKT =SIZE(PTHLM,3)
IKTE =IKT-JPVEXT_TURB
@@ -544,7 +609,6 @@ END IF
!
ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM)/PDZZ
ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
-
!
IF (LHGRAD) THEN
! Compute the Leonard terms for thl
@@ -601,21 +665,38 @@ IF (GFTHR) THEN
ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTHR_O_DDTDZ(Z3RDMOMENT,PREDTH1,&
& PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR)
END IF
+! compute interface flux
+IF (LCOUPLES) THEN ! Autocoupling O-A LES
+ IF (LOCEAN) THEN ! ocean model in coupled case
+ ZF(:,:,IKE) = (XSSTFL_C(:,:,1)+XSSRFL_C(:,:,1)) &
+ *0.5* ( 1. + PRHODJ(:,:,KKU)/PRHODJ(:,:,IKE) )
+ ELSE ! atmosph model in coupled case
+ ZF(:,:,IKB) = XSSTFL_C(:,:,1) &
+ *0.5* ( 1. + PRHODJ(:,:,KKA)/PRHODJ(:,:,IKB) )
+ ENDIF
+!
+ELSE ! No coupling O and A cases
+ ! atmosp bottom
+ !*In 3D, a part of the flux goes vertically,
+ ! and another goes horizontally (in presence of slopes)
+ !*In 1D, part of energy released in horizontal flux is taken into account in the vertical part
+ IF (HTURBDIM=='3DIM') THEN
+ ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
+ * PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ ELSE
+ ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
+ / PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ END IF
!
-!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
-! (in presence of slopes)
-!* in 1DIM case, the part of energy released in horizontal flux
-! is taken into account in the vertical part
-!
-IF (HTURBDIM=='3DIM') THEN
- ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
- * PDIRCOSZW(:,:) &
- * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
-ELSE
- ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
- / PDIRCOSZW(:,:) &
- * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
-END IF
+ IF (LOCEAN) THEN
+ ZF(:,:,IKE) = XSSTFL(:,:) *0.5*(1. + PRHODJ(:,:,KKU) / PRHODJ(:,:,IKE))
+ ELSE !end ocean case (in nocoupled case)
+ ! atmos top
+ ZF(:,:,IKE)=0.
+ END IF
+END IF !end no coupled cases
!
! Compute the split conservative potential temperature at t+deltat
CALL TRIDIAG_THERMO(KKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
@@ -650,15 +731,25 @@ IF (LHGRAD) THEN
END IF
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+IF (LOCEAN) THEN
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+END IF
!
- DO JK=IKTB+1,IKTE-1
- PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
- END DO
- PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
- PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
+DO JK=IKTB+1,IKTE-1
+ PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+END DO
+!
+PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
+!
+IF (LOCEAN) THEN
+ PWTH(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+KKL))
+ PWTH(:,:,KKA)=0.
+ PWTH(:,:,KKU)=ZFLXZ(:,:,KKU)
+ELSE
PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
-
-
+ PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
+END IF
+!
IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
! stores the conservative potential temperature vertical flux
TZFIELD%CMNHNAME = 'THW_FLX'
@@ -675,19 +766,27 @@ IF ( OTURB_FLX .AND. tpfile%lopened ) THEN
END IF
!
! Contribution of the conservative temperature flux to the buoyancy flux
-IF (KRR /= 0) THEN
- PTP(:,:,:) = PBETA * MZF( MZM(PETHETA) * ZFLXZ )
- PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+IF (LOCEAN) THEN
+ PTP(:,:,:)= XG*XALPHAOC * MZF(ZFLXZ )
ELSE
- PTP(:,:,:)= PBETA * MZF( ZFLXZ )
-END IF
+ IF (KRR /= 0) THEN
+ PTP(:,:,:) = PBETA * MZF( MZM(PETHETA) * ZFLXZ )
+ PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ ELSE
+ PTP(:,:,:)= PBETA * MZF( ZFLXZ )
+ END IF
+END IF
!
! Buoyancy flux at flux points
!
PWTHV = MZM(PETHETA) * ZFLXZ
PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
!
+IF (LOCEAN) THEN
+ ! temperature contribution to Buy flux
+ PWTHV(:,:,IKE) = PETHETA(:,:,IKE) * ZFLXZ(:,:,IKE)
+END IF
!* 2.3 Partial vertical divergence of the < Rc w > flux
!
IF ( KRRL >= 1 ) THEN
@@ -807,21 +906,42 @@ IF (KRR /= 0) THEN
& PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR)
END IF
!
- !* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
- ! (in presence of slopes)
- !* in 1DIM case, the part of energy released in horizontal flux
- ! is taken into account in the vertical part
- !
- IF (HTURBDIM=='3DIM') THEN
- ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
- * PDIRCOSZW(:,:) &
- * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
- ELSE
- ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
- / PDIRCOSZW(:,:) &
- * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
- END IF
+ ! compute interface flux
+ IF (LCOUPLES) THEN ! coupling NH O-A
+ IF (LOCEAN) THEN ! ocean model in coupled case
+ ! evap effect on salinity to be added later !!!
+ ZF(:,:,IKE) = 0.
+ ELSE ! atmosph model in coupled case
+ ZF(:,:,IKB) = 0.
+ ! AJOUTER FLUX EVAP SUR MODELE ATMOS
+ ENDIF
!
+ ELSE ! No coupling NH OA case
+ ! atmosp bottom
+ !* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
+ ! (in presence of slopes)
+ !* in 1DIM case, the part of energy released in horizontal flux
+ ! is taken into account in the vertical part
+ !
+ IF (HTURBDIM=='3DIM') THEN
+ ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
+ * PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ ELSE
+ ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
+ / PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ END IF
+ !
+ IF (LOCEAN) THEN
+ ! General ocean case
+ ! salinity/evap effect to be added later !!!!!
+ ZF(:,:,IKE) = 0.
+ ELSE !end ocean case (in nocoupled case)
+ ! atmos top
+ ZF(:,:,IKE)=0.
+ END IF
+ END IF!end no coupled cases
! Compute the split conservative potential temperature at t+deltat
CALL TRIDIAG_THERMO(KKA,KKU,KKL,PRM(:,:,:,1),ZF,ZDFDDRDZ,PTSTEP,PIMPL,&
PDZZ,PRHODJ,PRP)
@@ -882,15 +1002,22 @@ IF (KRR /= 0) THEN
END IF
!
! Contribution of the conservative water flux to the Buoyancy flux
- ZA(:,:,:) = PBETA * MZF( MZM(PEMOIST) * ZFLXZ )
- ZA(:,:,IKB) = PBETA(:,:,IKB) * PEMOIST(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
- PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
+ IF (LOCEAN) THEN
+ ZA(:,:,:)= -XG*XBETAOC * MZF(ZFLXZ )
+ ELSE
+ ZA(:,:,:) = PBETA * MZF( MZM(PEMOIST) * ZFLXZ )
+ ZA(:,:,IKB) = PBETA(:,:,IKB) * PEMOIST(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
+ END IF
!
! Buoyancy flux at flux points
!
PWTHV = PWTHV + MZM(PEMOIST) * ZFLXZ
PWTHV(:,:,IKB) = PWTHV(:,:,IKB) + PEMOIST(:,:,IKB) * ZFLXZ(:,:,IKB)
+ IF (LOCEAN) THEN
+ PWTHV(:,:,IKE) = PWTHV(:,:,IKE) + PEMOIST(:,:,IKE)* ZFLXZ(:,:,IKE)
+ END IF
!
!* 3.3 Complete vertical divergence of the < Rc w > flux
!
@@ -973,6 +1100,9 @@ IF ( ((OTURB_FLX .AND. tpfile%lopened) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
END IF
!
END IF !end of <w Rc>
+IF (LOCEAN.AND.LDEEPOC) THEN
+ DEALLOCATE(ZXHAT_ll,ZYHAT_ll)
+END IF
!
!----------------------------------------------------------------------------
END SUBROUTINE TURB_VER_THERMO_FLUX
diff --git a/src/MNH/ver_int_thermo.f90 b/src/MNH/ver_int_thermo.f90
index 1d3424a943e422cf37b2ec52dd4651470c25a108..7ec9f531dfe76f77385a65aded9d65e877339198 100644
--- a/src/MNH/ver_int_thermo.f90
+++ b/src/MNH/ver_int_thermo.f90
@@ -136,6 +136,7 @@ END MODULE MODI_VER_INT_THERMO
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 22/02/2019: replace Hollerith edit descriptor (deleted from Fortran 95 standard)
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! Jean-Luc Redelsperger 03/2021 OCEAN LES case
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -145,6 +146,7 @@ USE MODD_ARGSLIST_ll, ONLY: LIST_ll
USE MODD_CONF
USE MODD_CONF_n
USE MODD_CST
+USE MODD_DYN_n, ONLY : LOCEAN
USE MODD_GRID_n
USE MODD_IO, ONLY: TFILEDATA
USE MODD_LUNIT, ONLY: TLUOUT0
@@ -343,38 +345,42 @@ CALL MPPDB_CHECK3D(ZPMHPOHP_SH,"ver_int_thermo2a::ZPMHPOHP_SH",PRECISION)
!* 3.1 Computation of the shift profile
! --------------------------------
!
-ZTHVCLIMGR=3.5E-3 ! K/m
-CALL FREE_ATM_PROFILE(TPFILE,PTHV_MX,PZMASS_MX,PZS_LS,PZSMT_LS,ZTHVCLIMGR,ZTHV_FREE,ZZ_FREE)
-CALL MPPDB_CHECK3D(ZTHV_FREE,"VER_INT_THERMO:ZTHV_FREE",PRECISION)
+IF (LOCEAN) THEN
+ ZTHV_SH(:,:,:) = PTHV_MX(:,:,:)
+ELSE
+ ZTHVCLIMGR=3.5E-3 ! K/m
+ CALL FREE_ATM_PROFILE(TPFILE,PTHV_MX,PZMASS_MX,PZS_LS,PZSMT_LS,ZTHVCLIMGR,ZTHV_FREE,ZZ_FREE)
+ CALL MPPDB_CHECK3D(ZTHV_FREE,"VER_INT_THERMO:ZTHV_FREE",PRECISION)
!
!* 3.2 Computation of the value of thetav on the shifted grid
! ------------------------------------------------------
!
-CALL COEF_VER_INTERP_LIN(ZZ_FREE(:,:,:),PZMASS_MX(:,:,:))
-ZTHV_FREE_MX(:,:,:)=VER_INTERP_LIN(ZTHV_FREE(:,:,:),NKLIN(:,:,:),XCOEFLIN(:,:,:))
-CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"VER_INT_THERMO:ZTHV_FREE_MX",PRECISION)
+ CALL COEF_VER_INTERP_LIN(ZZ_FREE(:,:,:),PZMASS_MX(:,:,:))
+ ZTHV_FREE_MX(:,:,:)=VER_INTERP_LIN(ZTHV_FREE(:,:,:),NKLIN(:,:,:),XCOEFLIN(:,:,:))
+ CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"VER_INT_THERMO:ZTHV_FREE_MX",PRECISION)
!
-!20131113 add update_halo here
-CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"ver_int_thermo3a::ZTHV_FREE_MX",PRECISION)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHV_FREE_MX, 'VER_INT_THERMO::ZTHV_FREE_MX' )
- CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_ll)
-!20131112 check3d
-CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"ver_int_thermo2a::ZTHV_FREE_MX",PRECISION)
+ !20131113 add update_halo here
+ CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"ver_int_thermo3a::ZTHV_FREE_MX",PRECISION)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHV_FREE_MX, 'VER_INT_THERMO::ZTHV_FREE_MX' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !20131112 check3d
+ CALL MPPDB_CHECK3D(ZTHV_FREE_MX,"ver_int_thermo2a::ZTHV_FREE_MX",PRECISION)
!
-CALL COEF_VER_INTERP_LIN(ZZ_FREE(:,:,:),ZZMASS_SH(:,:,:))
-ZTHV_FREE_SH(:,:,:)=VER_INTERP_LIN(ZTHV_FREE(:,:,:),NKLIN(:,:,:),XCOEFLIN(:,:,:))
-CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"VER_INT_THERMO:ZTHV_FREE_SH",PRECISION)
+ CALL COEF_VER_INTERP_LIN(ZZ_FREE(:,:,:),ZZMASS_SH(:,:,:))
+ ZTHV_FREE_SH(:,:,:)=VER_INTERP_LIN(ZTHV_FREE(:,:,:),NKLIN(:,:,:),XCOEFLIN(:,:,:))
+ CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"VER_INT_THERMO:ZTHV_FREE_SH",PRECISION)
!
-!20131113 add update_halo here
-CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"ver_int_thermo3a::ZTHV_FREE_SH",PRECISION)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHV_FREE_SH, 'VER_INT_THERMO::ZTHV_FREE_SH' )
- CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_ll)
-!20131112 check3d
-CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"ver_int_thermo2a::ZTHV_FREE_SH",PRECISION)
-!
-ZTHV_SH(:,:,:) = PTHV_MX(:,:,:) - ZTHV_FREE_MX(:,:,:) + ZTHV_FREE_SH(:,:,:)
+ !20131113 add update_halo here
+ CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"ver_int_thermo3a::ZTHV_FREE_SH",PRECISION)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHV_FREE_SH, 'VER_INT_THERMO::ZTHV_FREE_SH' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !20131112 check3d
+ CALL MPPDB_CHECK3D(ZTHV_FREE_SH,"ver_int_thermo2a::ZTHV_FREE_SH",PRECISION)
+!
+ ZTHV_SH(:,:,:) = PTHV_MX(:,:,:) - ZTHV_FREE_MX(:,:,:) + ZTHV_FREE_SH(:,:,:)
+END IF
!
!20131113 add update_halo here
CALL MPPDB_CHECK3D(ZTHV_SH,"ver_int_thermo3a::ZTHV_SH",PRECISION)
@@ -393,29 +399,37 @@ CALL MPPDB_CHECK3D(ZTHV_SH,"ver_int_thermo2a::ZTHV_SH",PRECISION)
!* 4.1 Computation of relative humidity on the mixed grid
! --------------------------------------------------
!
-ZRV_MX(:,:,:)=MAX(PR_MX(:,:,:,1),1.E-10)
-ZTH_MX(:,:,:)=PTHV_MX(:,:,:)*(1.+WATER_SUM(PR_MX(:,:,:,:)))/(1.+XRV/XRD*ZRV_MX(:,:,:))
-!
-!20131113 add update_halo here
-CALL MPPDB_CHECK3D(ZTH_MX,"ver_int_thermo4a::ZTH_MX",PRECISION)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTH_MX, 'VER_INT_THERMO::ZTH_MX' )
- CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_ll)
-!20131112 check3d
-CALL MPPDB_CHECK3D(ZTH_MX,"ver_int_thermo4b::ZTH_MX",PRECISION)
+IF (LOCEAN) THEN
+ ZRV_MX(:,:,:) = PR_MX(:,:,:,1)
+ ZTH_MX(:,:,:) = PTHV_MX(:,:,:)
+ ZT_MX(:,:,:) = ZTH_MX(:,:,:)
+ ZES_MX(:,:,:) = SM_FOES(ZT_MX(:,:,:))
+ ZHU_MX(:,:,:) = 1.E-10
+ELSE
+ ZRV_MX(:,:,:)=MAX(PR_MX(:,:,:,1),1.E-10)
+ ZTH_MX(:,:,:)=PTHV_MX(:,:,:)*(1.+WATER_SUM(PR_MX(:,:,:,:)))/(1.+XRV/XRD*ZRV_MX(:,:,:))
!
-ZT_MX(:,:,:)=ZTH_MX(:,:,:)*ZEXNMASS_MX(:,:,:)
+ !20131113 add update_halo here
+ CALL MPPDB_CHECK3D(ZTH_MX,"ver_int_thermo4a::ZTH_MX",PRECISION)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTH_MX, 'VER_INT_THERMO::ZTH_MX' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !20131112 check3d
+ CALL MPPDB_CHECK3D(ZTH_MX,"ver_int_thermo4b::ZTH_MX",PRECISION)
!
-!20131113 add update_halo here
-CALL MPPDB_CHECK3D(ZT_MX,"ver_int_thermo4a::ZT_MX",PRECISION)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, ZT_MX, 'VER_INT_THERMO::ZT_MX' )
- CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_ll)
-!20131112 check3d
-CALL MPPDB_CHECK3D(ZT_MX,"ver_int_thermo4b::ZT_MX",PRECISION)
+ ZT_MX(:,:,:)=ZTH_MX(:,:,:)*ZEXNMASS_MX(:,:,:)
!
-ZES_MX(:,:,:)=SM_FOES(ZT_MX(:,:,:))
-ZHU_MX(:,:,:)=100.*ZP_MX(:,:,:)/(XRD/XRV/ZRV_MX(:,:,:)+1.)/ZES_MX(:,:,:)
+ !20131113 add update_halo here
+ CALL MPPDB_CHECK3D(ZT_MX,"ver_int_thermo4a::ZT_MX",PRECISION)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZT_MX, 'VER_INT_THERMO::ZT_MX' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !20131112 check3d
+ CALL MPPDB_CHECK3D(ZT_MX,"ver_int_thermo4b::ZT_MX",PRECISION)
+!
+ ZES_MX(:,:,:)=SM_FOES(ZT_MX(:,:,:))
+ ZHU_MX(:,:,:)=100.*ZP_MX(:,:,:)/(XRD/XRV/ZRV_MX(:,:,:)+1.)/ZES_MX(:,:,:)
+END IF
!
!20131113 add update_halo here
CALL MPPDB_CHECK3D(ZHU_MX,"ver_int_thermo4a::ZHU_MX",PRECISION)
@@ -549,17 +563,22 @@ END DO
CALL COMPUTE_EXNER_FROM_TOP(PTHV,XZZ,PEXNTOP2D,ZHEXN,ZHEXNMASS)
ZP(:,:,:) = PPMHP(:,:,:) + XP00 * ZHEXNMASS(:,:,:) ** (XCPD/XRD)
!
-!20131113 add update_halo here
-!CALL MPPDB_CHECK3D(ZP,"ver_int_thermo6a::ZP",PRECISION)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, ZP, 'VER_INT_THERMO::ZP' )
-CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
-CALL CLEANLIST_ll(TZFIELDS_ll)
-!20131112 check3d
-CALL MPPDB_CHECK3D(ZP,"ver_int_thermo6b::ZP",PRECISION)
-!
-PR(:,:,:,1)=SM_PMR_HU(ZP(:,:,:), &
- PTHV(:,:,:)*(ZP(:,:,:)/XP00)**(XRD/XCPD), &
- ZHU(:,:,:),PR(:,:,:,:),KITERMAX=100)
+IF (LOCEAN) THEN
+ !no interpolation for salinity
+ PR(:,:,:,1)=PR_MX(:,:,:,1)
+ELSE
+ !20131113 add update_halo here
+ !CALL MPPDB_CHECK3D(ZP,"ver_int_thermo6a::ZP",PRECISION)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZP, 'VER_INT_THERMO::ZP' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !20131112 check3d
+ CALL MPPDB_CHECK3D(ZP,"ver_int_thermo6b::ZP",PRECISION)
+!
+ PR(:,:,:,1)=SM_PMR_HU(ZP(:,:,:), &
+ PTHV(:,:,:)*(ZP(:,:,:)/XP00)**(XRD/XCPD), &
+ ZHU(:,:,:),PR(:,:,:,:),KITERMAX=100)
+END IF
CALL ADD3DFIELD_ll( TZFIELDS_ll, PR(:,:,:,1), 'VER_INT_THERMO::PR(:,:,:,1)' )
CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS_ll)
diff --git a/src/MNH/write_les_budgetn.f90 b/src/MNH/write_les_budgetn.f90
index da600ed6365bedd7889cd09f4a9bfebe804ddafd..a27560c73f1df8c904ef9c7a4d906a321c54f962 100644
--- a/src/MNH/write_les_budgetn.f90
+++ b/src/MNH/write_les_budgetn.f90
@@ -45,13 +45,15 @@ subroutine Write_les_budget_n( tpdiafile )
!! 06/11/02 (V. Masson) new LES budgets
! P. Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 15/10/2020: restructure Les_diachro calls to use tfield_metadata_base type
+! JL Redelsperger 03/21 modif buoyancy flix for OCEAN LES case
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
use modd_conf_n, only: luserv
-use modd_cst, only: xg
+use modd_cst, only: xg, XALPHAOC
+USE MODD_DYN_n, ONLY : LOCEAN
use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_TIME, &
NMNHDIM_BUDGET_TERM, NMNHDIM_UNUSED, &
tfield_metadata_base, TYPEREAL
@@ -763,6 +765,9 @@ ELSE
ZLES_BUDGET(:,:,ILES) = XG * XLES_SUBGRID_ThlThv(:,:,1) &
/ XLES_MEAN_Th (:,:,1)
END IF
+IF (LOCEAN) THEN
+ ZLES_BUDGET(:,:,ILES) = XG * XLES_SUBGRID_ThlThv(:,:,1) *XALPHAOC
+END IF
!
!* 3.6 residual of subgrid budget
! --------------------------
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 4448d579bf7d4e31da5e0cedbb4c5eb226440a45..50131fea21510726da1ce8d84a5489d6f0e5f07f 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -146,6 +146,7 @@ END MODULE MODI_WRITE_LFIFM1_FOR_DIAG
! S Bielli 02/2019: sea salt: significant sea wave height influences salt emission; 5 salt modes
! P. Wautelet 18/03/2020: remove ICE2 option
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
+! J.L Redelsperger 03/2021 Adding OCEAN LES Case and Autocoupled O-A LES
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -455,6 +456,8 @@ CALL IO_Field_write(TPFILE,'DTSEG',TDTSEG)
!
CALL IO_Field_write(TPFILE,'CARTESIAN',LCARTESIAN)
CALL IO_Field_write(TPFILE,'LBOUSS', LBOUSS)
+CALL IO_Field_write(TPFILE,'LOCEAN', LOCEAN)
+CALL IO_Field_write(TPFILE,'LCOUPLES', LCOUPLES)
!
IF (LCARTESIAN .AND. LWIND_ZM) THEN
LWIND_ZM=.FALSE.
@@ -462,9 +465,15 @@ IF (LCARTESIAN .AND. LWIND_ZM) THEN
END IF
!* 1.4 Reference state variables :
!
-CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZ)
-CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZ)
-CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOP)
+IF (LCOUPLES.AND.LOCEAN) THEN
+ CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZO)
+ CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZO)
+ CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOPO)
+ELSE
+ CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZ)
+ CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZ)
+ CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOP)
+END IF
!
CALL IO_Field_write(TPFILE,'RHODREF',XRHODREF)
CALL IO_Field_write(TPFILE,'THVREF', XTHVREF)
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index ee394f919dee89150d40158838f3e4f19adee935..03363082cf5405f48c941abbd35769fcc40d3ca7 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -178,6 +178,7 @@ END MODULE MODI_WRITE_LFIFM_n
! T.Nagel 02/2021 : Add turbulence recycling
! P. Wautelet 10/03/2021: use scalar variable names for dust and salt
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
+! J.L. Redelsperger 03/2021: Add OCEAN and auto-coupled O-A LES cases
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -229,6 +230,7 @@ USE MODD_HURR_FIELD_n
USE MODD_PREP_REAL, ONLY: CDUMMY_2D, XDUMMY_2D
USE MODD_DUST
USE MODD_SALT
+USE MODD_OCEANH
USE MODD_PASPOL
#ifdef MNH_FOREFIRE
USE MODD_FOREFIRE
@@ -435,6 +437,8 @@ CALL IO_Field_write(TPFILE,'L2D', L2D)
CALL IO_Field_write(TPFILE,'PACK', LPACK)
CALL IO_Field_write(TPFILE,'CARTESIAN',LCARTESIAN)
CALL IO_Field_write(TPFILE,'LBOUSS', LBOUSS)
+CALL IO_Field_write(TPFILE,'LOCEAN', LOCEAN)
+CALL IO_Field_write(TPFILE,'LCOUPLES', LCOUPLES)
!
CALL IO_Field_write(TPFILE,'SURF', CSURF)
CALL IO_Field_write(TPFILE,'CPL_AROME',LCPL_AROME)
@@ -1562,9 +1566,15 @@ END IF
!
!* 1.5 Reference state variables :
!
-CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZ)
-CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZ)
-CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOP)
+IF (LCOUPLES.AND.LOCEAN) THEN
+ CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZO)
+ CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZO)
+ CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOPO)
+ELSE
+ CALL IO_Field_write(TPFILE,'RHOREFZ',XRHODREFZ)
+ CALL IO_Field_write(TPFILE,'THVREFZ',XTHVREFZ)
+ CALL IO_Field_write(TPFILE,'EXNTOP', XEXNTOP)
+END IF
!
!
!* 1.6 Tendencies
@@ -1943,8 +1953,85 @@ IF(LBLOWSNOW) THEN
END IF
ENDIF
!
-!* 1.11 Forcing variables
+!* 1.11 Ocean LES variables
+!
+ IF ((.NOT.LCOUPLES).AND.LOCEAN) THEN
+WRITE (ILUOUT,*) 'LOCEAN NFRCLT', LOCEAN,NFRCLT,XSSUFL_T(1),XSSVFL_T(2),XSSTFL_T(1),XSSOLA_T(2)
+TZFIELD%CMNHNAME = 'NFRCLT'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'NFRCLT'
+TZFIELD%CUNITS = 'number of forc'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'nb de flux sfc forcant LES ocean'
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEINT
+TZFIELD%NDIMS = 0
+TZFIELD%LTIMEDEP = .FALSE.
+CALL IO_Field_write(TPFILE,TZFIELD,NFRCLT)
+!
+TZFIELD%CMNHNAME = 'NINFRT'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'NINFRT'
+TZFIELD%CUNITS = 'interv between forc'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'int between flux sfc forcant LES ocean'
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEINT
+TZFIELD%NDIMS = 0
+TZFIELD%LTIMEDEP = .FALSE.
+CALL IO_Field_write(TPFILE,TZFIELD,NINFRT)
+!
+TZFIELD%CMNHNAME = 'SSUFL_T'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'SSUFL'
+TZFIELD%CUNITS = 'kg m-1 s-1'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'sfc stress along U to force ocean LES '
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEREAL
+TZFIELD%NDIMS = 1
+TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_Field_write(TPFILE,TZFIELD,XSSUFL_T(:))
+!
+TZFIELD%CMNHNAME = 'SSVFL_T'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'SSVFL'
+TZFIELD%CUNITS = 'kg m-1 s-1'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'sfc stress along V to force ocean LES '
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEREAL
+TZFIELD%NDIMS = 1
+TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_Field_write(TPFILE,TZFIELD,XSSVFL_T(:))
+!
+TZFIELD%CMNHNAME = 'SSTFL_T'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'SSTFL'
+TZFIELD%CUNITS = 'kg m3 K m s-1'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'sfc total heat flux to force ocean LES '
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEREAL
+TZFIELD%NDIMS = 1
+TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_Field_write(TPFILE,TZFIELD,XSSTFL_T(:))
+!
+TZFIELD%CMNHNAME = 'SSOLA_T'
+TZFIELD%CSTDNAME = ''
+TZFIELD%CLONGNAME = 'SSOLA'
+TZFIELD%CUNITS = 'kg m3 K m s-1'
+TZFIELD%CDIR = '--'
+TZFIELD%CCOMMENT = 'sfc solar flux to force ocean LES '
+TZFIELD%NGRID = 0
+TZFIELD%NTYPE = TYPEREAL
+TZFIELD%NDIMS = 1
+TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_Field_write(TPFILE,TZFIELD,XSSOLA_T(:))
+!
+END IF ! ocean sfc forcing end
!
+!* 1.12 Forcing variables
!
IF (LFORCING) THEN
!
@@ -2230,7 +2317,7 @@ IF ( L2D_REL_FRC ) THEN
ENDDO
ENDIF
!
-!* 1.11bis Eddy Fluxes variables ! Modif PP
+!* 1.13 Eddy Fluxes variables ! Modif PP
!
IF ( LTH_FLX ) THEN
CALL IO_Field_write(TPFILE,'VT_FLX',XVTH_FLUX_M)
@@ -2239,13 +2326,13 @@ END IF
!
IF ( LUV_FLX) CALL IO_Field_write(TPFILE,'VU_FLX',XVU_FLUX_M)
!
-!* 1.12 Balloon variables
+!* 1.14 Balloon variables
!
!
IF (LFLYER) CALL WRITE_BALLOON_n(TPFILE)
!
!
-!* 1.13 Filtered variables for hurricane initialization
+!* 1.15 Filtered variables for hurricane initialization
!
!
IF ( CPROGRAM=='REAL ' ) THEN
@@ -2290,7 +2377,7 @@ IF ( CPROGRAM=='REAL ' ) THEN
!
END IF
!
-!* 1.14 Dummy variables in PREP_REAL_CASE
+!* 1.16 Dummy variables in PREP_REAL_CASE
!
IF (ALLOCATED(CDUMMY_2D)) THEN
TZFIELD%CSTDNAME = ''
@@ -2311,7 +2398,7 @@ IF ( CPROGRAM=='REAL ' ) THEN
!
END IF
!
-!* 1.15 Wind turbine variables
+!* 1.17 Wind turbine variables
!
! i) Main
!