diff --git a/src/arome/aux/modd_frc.f90 b/src/arome/aux/modd_frc.F90
similarity index 100%
rename from src/arome/aux/modd_frc.f90
rename to src/arome/aux/modd_frc.F90
diff --git a/src/arome/aux/modd_oceanh.f90 b/src/arome/aux/modd_oceanh.F90
similarity index 100%
rename from src/arome/aux/modd_oceanh.f90
rename to src/arome/aux/modd_oceanh.F90
diff --git a/src/arome/aux/mode_gather_ll.F90 b/src/arome/aux/mode_gather_ll.F90
index dd922f659edf6f298379fa0a7acb29e7e3060acb..88c37bbd9acd428168479ea952092f5799ebbf36 100644
--- a/src/arome/aux/mode_gather_ll.F90
+++ b/src/arome/aux/mode_gather_ll.F90
@@ -1,12 +1,28 @@
MODULE MODE_GATHER_ll
IMPLICIT NONE
+
+INTERFACE GATHERALL_FIELD_ll
+ MODULE PROCEDURE &
+ GATHERALL_X1, GATHERALL_X3
+END INTERFACE
+
CONTAINS
-SUBROUTINE GATHERALL_FIELD_ll(HDIR,PSEND,PRECV,KRESP)
+SUBROUTINE GATHERALL_X3(HDIR,PSEND,PRECV,KRESP)
CHARACTER(LEN=*), INTENT(IN) :: HDIR
REAL,DIMENSION(:,:,:), INTENT(IN) :: PSEND
REAL,DIMENSION(:,:,:), INTENT(INOUT):: PRECV
INTEGER, INTENT(INOUT):: KRESP
CALL ABORT
-END SUBROUTINE GATHERALL_FIELD_ll
+END SUBROUTINE GATHERALL_X3
+!
+SUBROUTINE GATHERALL_X1(HDIR,PSEND,PRECV,KRESP)
+CHARACTER(LEN=*), INTENT(IN) :: HDIR
+REAL,DIMENSION(:), INTENT(IN) :: PSEND
+REAL,DIMENSION(:), INTENT(INOUT):: PRECV
+INTEGER, INTENT(INOUT):: KRESP
+
+CALL ABORT
+END SUBROUTINE GATHERALL_X1
+!
END MODULE MODE_GATHER_ll
diff --git a/src/common/aux/modd_blowsnow.f90 b/src/common/aux/modd_blowsnow.F90
similarity index 100%
rename from src/common/aux/modd_blowsnow.f90
rename to src/common/aux/modd_blowsnow.F90
diff --git a/src/common/aux/modd_dimn.f90 b/src/common/aux/modd_dimn.F90
similarity index 100%
rename from src/common/aux/modd_dimn.f90
rename to src/common/aux/modd_dimn.F90
diff --git a/src/common/aux/modd_gridn.f90 b/src/common/aux/modd_gridn.F90
similarity index 100%
rename from src/common/aux/modd_gridn.f90
rename to src/common/aux/modd_gridn.F90
diff --git a/src/common/aux/modd_metricsn.f90 b/src/common/aux/modd_metricsn.F90
similarity index 100%
rename from src/common/aux/modd_metricsn.f90
rename to src/common/aux/modd_metricsn.F90
diff --git a/src/common/aux/modd_ref.f90 b/src/common/aux/modd_ref.F90
similarity index 100%
rename from src/common/aux/modd_ref.f90
rename to src/common/aux/modd_ref.F90
diff --git a/src/common/aux/modd_turbn.f90 b/src/common/aux/modd_turbn.F90
similarity index 61%
rename from src/common/aux/modd_turbn.f90
rename to src/common/aux/modd_turbn.F90
index 8c35fd9d4be7bd61167ee2e6aba4a4fb40e521a5..03cb317eb33be9eab48ea00e6d47556336501b3b 100644
--- a/src/common/aux/modd_turbn.f90
+++ b/src/common/aux/modd_turbn.F90
@@ -49,8 +49,6 @@
!
USE MODD_PARAMETERS, ONLY: JPMODELMAX
IMPLICIT NONE
-
-TYPE TURB_t
!
!
REAL :: XIMPL ! implicitness degree for the vertical terms of
@@ -107,105 +105,4 @@ TYPE TURB_t
! negative value : applied everywhere
! 0.000001 applied only inside the clouds ri+rc > 10**-6 kg/kg
!
-END TYPE TURB_t
-
-TYPE(TURB_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: TURB_MODEL
-
-REAL, POINTER :: XIMPL=>NULL()
-REAL, POINTER :: XKEMIN=>NULL()
-REAL, POINTER :: XCEDIS=>NULL()
-REAL, POINTER :: XCADAP=>NULL()
-CHARACTER (LEN=4), POINTER :: CTURBLEN=>NULL()
-CHARACTER (LEN=4), POINTER :: CTURBDIM=>NULL()
-LOGICAL, POINTER :: LTURB_FLX=>NULL()
-LOGICAL, POINTER :: LTURB_DIAG=>NULL()
-LOGICAL, POINTER :: LSUBG_COND=>NULL()
-LOGICAL, POINTER :: LSIGMAS=>NULL()
-LOGICAL, POINTER :: LSIG_CONV=>NULL()
-LOGICAL, POINTER :: LRMC01=>NULL()
-CHARACTER(LEN=4),POINTER :: CTOM=>NULL()
-CHARACTER(LEN=4),POINTER :: CSUBG_AUCV=>NULL()
-CHARACTER(LEN=80),POINTER :: CSUBG_AUCV_RI=>NULL()
-CHARACTER(LEN=80),POINTER :: CCONDENS=>NULL()
-CHARACTER(LEN=4),POINTER :: CLAMBDA3=>NULL()
-CHARACTER(LEN=80),POINTER :: CSUBG_MF_PDF=>NULL()
-REAL, DIMENSION(:,:), POINTER :: XBL_DEPTH=>NULL()
-REAL, DIMENSION(:,:), POINTER :: XSBL_DEPTH=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XWTHVMF=>NULL()
-REAL, POINTER :: VSIGQSAT=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XDYP=>NULL()
-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()
-REAL, POINTER :: XALTHGRAD=>NULL()
-REAL, POINTER :: XCLDTHOLD=>NULL()
-
-CONTAINS
-
-SUBROUTINE TURB_GOTO_MODEL(KFROM, KTO)
-INTEGER, INTENT(IN) :: KFROM, KTO
-!
-! Save current state for allocated arrays
-!
-!TURB_MODEL(KFROM)%XBL_DEPTH=>XBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
-!TURB_MODEL(KFROM)%XSBL_DEPTH=>XSBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
-!TURB_MODEL(KFROM)%XWTHVMF=>XWTHVMF !Done in FIELDLIST_GOTO_MODEL
-TURB_MODEL(KFROM)%XDYP=>XDYP
-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
-XKEMIN=>TURB_MODEL(KTO)%XKEMIN
-XCEDIS=>TURB_MODEL(KTO)%XCEDIS
-XCADAP=>TURB_MODEL(KTO)%XCADAP
-CTURBLEN=>TURB_MODEL(KTO)%CTURBLEN
-CTURBDIM=>TURB_MODEL(KTO)%CTURBDIM
-LTURB_FLX=>TURB_MODEL(KTO)%LTURB_FLX
-LTURB_DIAG=>TURB_MODEL(KTO)%LTURB_DIAG
-LSUBG_COND=>TURB_MODEL(KTO)%LSUBG_COND
-LSIGMAS=>TURB_MODEL(KTO)%LSIGMAS
-LSIG_CONV=>TURB_MODEL(KTO)%LSIG_CONV
-LRMC01=>TURB_MODEL(KTO)%LRMC01
-CTOM=>TURB_MODEL(KTO)%CTOM
-CSUBG_AUCV=>TURB_MODEL(KTO)%CSUBG_AUCV
-CSUBG_AUCV_RI=>TURB_MODEL(KTO)%CSUBG_AUCV_RI
-CCONDENS=>TURB_MODEL(KTO)%CCONDENS
-CLAMBDA3=>TURB_MODEL(KTO)%CLAMBDA3
-CSUBG_MF_PDF=>TURB_MODEL(KTO)%CSUBG_MF_PDF
-!XBL_DEPTH=>TURB_MODEL(KTO)%XBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
-!XSBL_DEPTH=>TURB_MODEL(KTO)%XSBL_DEPTH !Done in FIELDLIST_GOTO_MODEL
-!XWTHVMF=>TURB_MODEL(KTO)%XWTHVMF !Done in FIELDLIST_GOTO_MODEL
-VSIGQSAT=>TURB_MODEL(KTO)%VSIGQSAT
-XDYP=>TURB_MODEL(KTO)%XDYP
-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
-XALTHGRAD=>TURB_MODEL(KTO)%XALTHGRAD
-XCLDTHOLD=>TURB_MODEL(KTO)%XCLDTHOLD
-
-END SUBROUTINE TURB_GOTO_MODEL
-
END MODULE MODD_TURB_n
diff --git a/src/common/aux/tools.f90 b/src/common/aux/tools.F90
similarity index 100%
rename from src/common/aux/tools.f90
rename to src/common/aux/tools.F90
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index 2004edfcdf2108f3d7bc02c9ca25643aea7d1d5e..6a84eb9b3cb90b711360c9ffa7f765e65ca905b8 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -962,12 +962,10 @@ D_M3_WTH_W2TH_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_W2TH_O_DDTDZ(:,:,IKE)
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2TH_O_DDTDZ',1,ZHOOK_HANDLE)
END FUNCTION D_M3_WTH_W2TH_O_DDTDZ
!----------------------------------------------------------------------------
-FUNCTION M3_WTH_W2R(KKA,KKU,KKL,PREDTH1,PREDR1,PD,PKEFF,PTKE,PBLL_O_E,PEMOIST,PDTDZ)
+FUNCTION M3_WTH_W2R(KKA,KKU,KKL,PD,PKEFF,PTKE,PBLL_O_E,PEMOIST,PDTDZ)
INTEGER, INTENT(IN) :: KKA
INTEGER, INTENT(IN) :: KKU
INTEGER, INTENT(IN) :: KKL
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1
REAL, DIMENSION(:,:,:), INTENT(IN) :: PD
REAL, DIMENSION(:,:,:), INTENT(IN) :: PKEFF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
@@ -1018,12 +1016,10 @@ D_M3_WTH_W2R_O_DDTDZ(:,:,IKE+1)=D_M3_WTH_W2R_O_DDTDZ(:,:,IKE)
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WTH_W2R_O_DDTDZ',1,ZHOOK_HANDLE)
END FUNCTION D_M3_WTH_W2R_O_DDTDZ
!----------------------------------------------------------------------------
-FUNCTION M3_WTH_WR2(KKA,KKU,KKL,PREDTH1,PREDR1,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTDZ)
+FUNCTION M3_WTH_WR2(KKA,KKU,KKL,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTDZ)
INTEGER, INTENT(IN) :: KKA
INTEGER, INTENT(IN) :: KKU
INTEGER, INTENT(IN) :: KKL
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1
REAL, DIMENSION(:,:,:), INTENT(IN) :: PD
REAL, DIMENSION(:,:,:), INTENT(IN) :: PKEFF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
@@ -1798,12 +1794,10 @@ D_M3_WR_W2R_O_DDRDZ = D_M3_WTH_W2TH_O_DDTDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PBLL_O
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WR_W2R_O_DDRDZ',1,ZHOOK_HANDLE)
END FUNCTION D_M3_WR_W2R_O_DDRDZ
!----------------------------------------------------------------------------
-FUNCTION M3_WR_W2TH(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PBLL_O_E,PETHETA,PDRDZ)
+FUNCTION M3_WR_W2TH(KKA,KKU,KKL,PD,PKEFF,PTKE,PBLL_O_E,PETHETA,PDRDZ)
INTEGER, INTENT(IN) :: KKA
INTEGER, INTENT(IN) :: KKU
INTEGER, INTENT(IN) :: KKL
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1
REAL, DIMENSION(:,:,:), INTENT(IN) :: PD
REAL, DIMENSION(:,:,:), INTENT(IN) :: PKEFF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
@@ -1814,7 +1808,7 @@ FUNCTION M3_WR_W2TH(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PBLL_O_E,PETHETA,PD
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WR_W2TH',0,ZHOOK_HANDLE)
-M3_WR_W2TH = M3_WTH_W2R(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PBLL_O_E,PETHETA,PDRDZ)
+M3_WR_W2TH = M3_WTH_W2R(KKA,KKU,KKL,PD,PKEFF,PTKE,PBLL_O_E,PETHETA,PDRDZ)
!
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WR_W2TH',1,ZHOOK_HANDLE)
END FUNCTION M3_WR_W2TH
@@ -1839,12 +1833,10 @@ D_M3_WR_W2TH_O_DDRDZ = D_M3_WTH_W2R_O_DDTDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:D_M3_WR_W2TH_O_DDRDZ',1,ZHOOK_HANDLE)
END FUNCTION D_M3_WR_W2TH_O_DDRDZ
!----------------------------------------------------------------------------
-FUNCTION M3_WR_WTH2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDRDZ)
+FUNCTION M3_WR_WTH2(KKA,KKU,KKL,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDRDZ)
INTEGER, INTENT(IN) :: KKA
INTEGER, INTENT(IN) :: KKU
INTEGER, INTENT(IN) :: KKL
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1
REAL, DIMENSION(:,:,:), INTENT(IN) :: PD
REAL, DIMENSION(:,:,:), INTENT(IN) :: PKEFF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
@@ -1858,7 +1850,7 @@ FUNCTION M3_WR_WTH2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WR_WTH2',0,ZHOOK_HANDLE)
-M3_WR_WTH2 = M3_WTH_WR2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDRDZ)
+M3_WR_WTH2 = M3_WTH_WR2(KKA,KKU,KKL,PD,PKEFF,PTKE,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDRDZ)
!
IF (LHOOK) CALL DR_HOOK('MODE_PRANDTL:M3_WR_WTH2',1,ZHOOK_HANDLE)
END FUNCTION M3_WR_WTH2
diff --git a/src/common/turb/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index 4473628de977d93638286ef291b179b7dc68e47a..b4e25e766dc3d99a601f39c39bd2d1612e709bfc 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -316,17 +316,12 @@ REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: &
ZFLXZ, & ! vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
-INTEGER :: ILENCH ! Length of comment string in LFIFM file
INTEGER :: IIB,IIE, & ! I index values for the Beginning and End
IJB,IJE, & ! mass points of the domain in the 3 direct.
IKB,IKE !
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV ! scalar loop counter
-CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
-CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM
! coefficients for the surface flux
@@ -340,15 +335,17 @@ TYPE(TFIELDDATA) :: TZFIELD
!
!* 1. PRELIMINARIES
! -------------
+ZA=XUNDEF
+PDP=XUNDEF
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
+!
IIU=SIZE(PUM,1)
-IIE=IIU-JPHEXT
-IIB=1+JPHEXT
IJU=SIZE(PUM,2)
-IJE=IJU-JPHEXT
-IJB=1+JPHEXT
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
IKT=SIZE(PUM,3)
@@ -360,9 +357,7 @@ IKTE=IKT-JPVEXT_TURB
ZSOURCE = 0.
ZFLXZ = 0.
ZCMFS = XCMFS
-IF (LHARAT)THEN
- ZCMFS=1.
-ENDIF
+IF (LHARAT) ZCMFS=1.
!
ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
! compute the coefficients for the uncentred gradient computation near the
@@ -371,9 +366,9 @@ ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
! With LHARATU length scale and TKE are at half levels so remove MZM
!
IF (LHARAT) THEN
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
@@ -394,7 +389,6 @@ ZA(:,:,:) = -PTSTEP * ZCMFS * &
MXM( ZKEFF ) * MXM(MZM(PRHODJ, KKA, KKU, KKL)) / &
MXM( PDZZ )**2
!
-IF (CPROGRAM/='AROME ') ZA(1,:,:)=ZA(IIE,:,:)
!
! Compute the source of U wind component
!
@@ -411,27 +405,50 @@ 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)) )
+ -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)) )
+ 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 splitted U at t+ deltat
+! Obtention of the split U at t+ deltat
!
CALL TRIDIAG_WIND(KKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MXM(PRHODJ),ZSOURCE,ZRES)
@@ -456,6 +473,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
@@ -485,7 +508,15 @@ PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
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
!
@@ -504,8 +535,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 &
@@ -535,6 +570,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
@@ -573,7 +623,6 @@ ZA(:,:,:) = - PTSTEP * ZCMFS * &
MYM( PDZZ )**2
!
!
-IF(CPROGRAM/='AROME ') ZA(:,1,:)=ZA(:,IJE,:)
!
! Compute the source of V wind component
! compute the coefficient between the vertical flux and the 2 components of the
@@ -589,26 +638,46 @@ 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 without coupling
+ ! 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 when coupling
+ ! input flux assumed to be in SI and at vorticity point
+ 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 splitted V at t+ deltat
+! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(KKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MYM(PRHODJ),ZSOURCE,ZRES)
!
@@ -632,6 +701,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'
@@ -663,6 +739,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
@@ -682,6 +766,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
@@ -709,9 +796,23 @@ IF(HTURBDIM=='3DIM') THEN
/(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB )) &
) &
* PDZY(:,:,IKB+KKL:IKB+KKL) &
- ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB))) &
+ ) / (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
@@ -737,15 +838,6 @@ IF(HTURBDIM=='3DIM') THEN
!
END IF
!
-! complete the dynamic production at the marginal points
-IF (CPROGRAM/='AROME ') THEN
- PDP(:,:,KKA)= -999.
- PDP(:,:,KKU)= -999.
- PDP(:,1,:)= PDP(:,IJE,:)
- PDP(:,IJE+1,:)= PDP(:,IJB,:)
- PDP(1,:,:)= PDP(IIE,:,:)
- PDP(IIE+1,:,:)= PDP(IIB,:,:)
-END IF
!
!----------------------------------------------------------------------------
!
diff --git a/src/common/turb/mode_turb_ver_sv_corr.F90 b/src/common/turb/mode_turb_ver_sv_corr.F90
index f140c0792e47cc3315f874e56bd9f2241598d69e..7de442316cf8bdb5210ef16ce366845d4da4eb8a 100644
--- a/src/common/turb/mode_turb_ver_sv_corr.F90
+++ b/src/common/turb/mode_turb_ver_sv_corr.F90
@@ -130,6 +130,13 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',0,ZHOOK_HANDLE)
CALL SECOND_MNH(ZTIME1)
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
+!
DO JSV=1,NSV
!
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
@@ -139,7 +146,7 @@ DO JSV=1,NSV
IF (LLES_CALL) THEN
! approximation: diagnosed explicitely (without implicit term)
ZFLXZ(:,:,:) = PPSI_SV(:,:,:,JSV)*GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)**2
- ZFLXZ(:,:,:) = XCHF / ZCSVD * PLM * PLEPS * MZF(ZFLXZ(:,:,:), KKA, KKU, KKL)
+ ZFLXZ(:,:,:) = ZCSV / ZCSVD * PLM * PLEPS * MZF(ZFLXZ(:,:,:), KKA, KKU, KKL)
CALL LES_MEAN_SUBGRID(-2.*ZCSVD*SQRT(PTKEM)*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV) )
CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) )
END IF
@@ -149,7 +156,7 @@ DO JSV=1,NSV
IF (LLES_CALL) THEN
! approximation: diagnosed explicitely (without implicit term)
ZA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM)
- ZFLXZ(:,:,:)= ( XCSHF * PPHI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ ZFLXZ(:,:,:)= ( XCSHF * PPHI3 + ZCSV * PPSI_SV(:,:,:,JSV) ) &
* GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL) &
* GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCTSVD) * MZF(ZFLXZ, KKA, KKU, KKL)
@@ -158,7 +165,7 @@ DO JSV=1,NSV
!
IF (KRR>=1) THEN
ZA(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM)
- ZFLXZ(:,:,:)= ( XCHF * PPSI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ ZFLXZ(:,:,:)= ( ZCSV * PPSI3 + ZCSV * PPSI_SV(:,:,:,JSV) ) &
* GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL) &
* GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCQSVD) * MZF(ZFLXZ, KKA, KKU, KKL)
diff --git a/src/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index 9ef220cf8e7ca6883f9e689bcfeaa8cdfd3041bf..41516708838d42b6a2bfd2851c00d66f6f17ebb7 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -290,9 +290,6 @@ REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
ZFLXZ, & ! vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
-INTEGER :: ILENCH ! Length of comment string in LFIFM file
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
@@ -300,8 +297,6 @@ INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV ! loop counters
INTEGER :: JK ! loop
INTEGER :: ISV ! number of scalar var.
-CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
-CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
!
REAL :: ZTIME1, ZTIME2
@@ -326,12 +321,17 @@ IKTB =1+JPVEXT_TURB
ISV=SIZE(PSVM,4)
!
IF (LHARAT) THEN
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
-
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
!----------------------------------------------------------------------------
!
!* 8. SOURCES OF PASSIVE SCALAR VARIABLES
@@ -342,15 +342,16 @@ DO JSV=1,ISV
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
!
! Preparation of the arguments for TRIDIAG
-IF (LHARAT) THEN
- ZA(:,:,:) = -PTSTEP* &
- ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
- PDZZ**2
-ELSE
- ZA(:,:,:) = -PTSTEP*XCHF*PPSI_SV(:,:,:,JSV) * &
- ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
- PDZZ**2
-ENDIF
+ IF (LHARAT) THEN
+ ZA(:,:,:) = -PTSTEP* &
+ ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
+ PDZZ**2
+ ELSE
+ ZA(:,:,:) = -PTSTEP*ZCSV*PPSI_SV(:,:,:,JSV) * &
+ ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
+ PDZZ**2
+ ENDIF
+ ZSOURCE(:,:,:) = 0.
!
! Compute the sources for the JSVth scalar variable
@@ -381,7 +382,7 @@ ENDIF
IF ( (OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL ) THEN
! Diagnostic of the cartesian vertical flux
!
- ZFLXZ(:,:,:) = -XCHF * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), KKA, KKU, KKL) / PDZZ * &
+ ZFLXZ(:,:,:) = -ZCSV * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), KKA, KKU, KKL) / PDZZ * &
DZM(PIMPL*ZRES(:,:,:) + PEXPL*PSVM(:,:,:,JSV), KKA, KKU, KKL)
! surface flux
!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index 7a79162a45544c8ea8ebb7eab1e52629cad69006..7b7b6a4d42973c39fa4168a2f241554811658d93 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -228,6 +228,7 @@ USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!
USE MODD_CST
+USE MODD_CONF, ONLY: CPROGRAM
USE MODD_CTURB
USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
USE MODD_GRID_n, ONLY: XZS, XXHAT, XYHAT
@@ -248,7 +249,7 @@ USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
USE MODI_GRADIENT_W
USE MODI_GRADIENT_M
-USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM
+USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM, MYF, MXF
USE MODI_TRIDIAG
USE MODI_LES_MEAN_SUBGRID
USE MODI_TRIDIAG_THERMO
@@ -364,19 +365,43 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
ZF, & ! Flux in dTh/dt =-dF/dz (evaluated at t-1)(or rt instead of Th)
ZDFDDTDZ, & ! dF/d(dTh/dz)
ZDFDDRDZ, & ! dF/d(dr/dz)
- Z3RDMOMENT ! 3 order term in flux or variance equation
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
-INTEGER :: ILENCH ! Length of comment string in LFIFM file
+ Z3RDMOMENT,& ! 3 order term in flux or variance equation
+ ZF_LEONARD,& ! Leonard terms
+ ZRWTHL, &
+ ZRWRNP, &
+ ZCLD_THOLD
+!
+REAL,DIMENSION(SIZE(XZS,1),SIZE(XZS,2),KKU) :: ZALT
+!
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
-CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
-CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+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)
!
-REAL :: ZTIME1, ZTIME2
!
+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
@@ -392,6 +417,38 @@ TYPE(TFIELDDATA) :: TZFIELD
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',0,ZHOOK_HANDLE)
+!
+! 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,IIU,IJU)
+ 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
IKTB =1+JPVEXT_TURB
@@ -405,11 +462,22 @@ GUSERV = (KRR/=0)
!
IF (LHARAT) THEN
! LHARAT so TKE and length scales at half levels!
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
+! Define a cloud mask with ri and rc (used after with a threshold) for Leonard terms
+!
+IF(LHGRAD) THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2) + PRM(:,:,:,4)
+ ELSE
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2)
+ END IF
+ END IF
+END IF
!
! Flags for 3rd order quantities
!
@@ -437,12 +505,22 @@ END IF
! Compute the turbulent flux F and F' at time t-dt.
!
IF (LHARAT) THEN
-ZF (:,:,:) = -ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
-ZDFDDTDZ(:,:,:) = -ZKEFF
+ ZF (:,:,:) = -ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ ZDFDDTDZ(:,:,:) = -ZKEFF
ELSE
-ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
-ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
-ENDIF
+ ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
+END IF
+!
+IF (LHGRAD) THEN
+ ! Compute the Leonard terms for thl
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_LEONARD (:,:,:)= XCOEFHGRADTHL*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX,XDZZ,XDZX,KKA,KKU,KKL))&
+ *MZM(GX_M_M(PTHLM(:,:,:),XDXX,XDZZ,XDZX,KKA, KKU, KKL), KKA, KKU, KKL) &
+ + MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,KKA,KKU,KKL)) &
+ *MZM(GY_M_M(PTHLM(:,:,:),XDYY,XDZZ,XDZY,KKA, KKU, KKL), KKA, KKU, KKL) )
+END IF
!
! Effect of 3rd order terms in temperature flux (at flux point)
!
@@ -466,7 +544,7 @@ END IF
!
! d(w'2r')/dz
IF (GFWR) THEN
- ZF = ZF + M3_WTH_W2R(KKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,&
+ ZF = ZF + M3_WTH_W2R(KKA,KKU,KKL,PD,ZKEFF,&
& PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * PFWR
ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2R_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
& PD,ZKEFF,PTKEM,PBLL_O_E,PEMOIST) * PFWR
@@ -474,7 +552,7 @@ END IF
!
! d(w'r'2)/dz
IF (GFR2) THEN
- ZF = ZF + M3_WTH_WR2(KKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,PTKEM,&
+ ZF = ZF + M3_WTH_WR2(KKA,KKU,KKL,PD,ZKEFF,PTKEM,&
& PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTH_DZ) * MZM(PFR2, KKA, KKU, KKL)
ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WR2_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,PD,&
& ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST) * MZM(PFR2, KKA, KKU, KKL)
@@ -489,46 +567,95 @@ IF (GFTHR) THEN
ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTHR_O_DDTDZ(Z3RDMOMENT,PREDTH1,&
& PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR, KKA, KKU, KKL)
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
+#ifdef REPRO48
+#else
+ ZF(:,:,IKE)=0.
+#endif
+ END IF
+END IF !end no coupled cases
!
-! Compute the splitted conservative potential temperature at t+deltat
+! Compute the split conservative potential temperature at t+deltat
CALL TRIDIAG_THERMO(KKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
PRHODJ,PTHLP)
!
! Compute the equivalent tendency for the conservative potential temperature
-PRTHLS(:,:,:)= PRTHLS(:,:,:) + &
- PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+!
+ZRWTHL(:,:,:)= PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZRWTHL(:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:), KKA, KKU, KKL)*ZF_LEONARD(:,:,:),XDZZ,&
+ KKA, KKU, KKL)
+ END WHERE
+END IF
+!
+PRTHLS(:,:,:)= PRTHLS(:,:,:) + ZRWTHL(:,:,:)
!
!* 2.2 Partial Thermal Production
!
! Conservative potential temperature flux :
!
ZFLXZ(:,:,:) = ZF &
- + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ
+ + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ
+! replace the flux by the Leonard terms
+IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZFLXZ(:,:,:) = ZF_LEONARD(:,:,:)
+ END WHERE
+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))
+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(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
-
+END IF
+!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
! stores the conservative potential temperature vertical flux
TZFIELD%CMNHNAME = 'THW_FLX'
@@ -545,34 +672,44 @@ 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, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
- PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+IF (LOCEAN) THEN
+ PTP(:,:,:)= XG*XALPHAOC * MZF(ZFLXZ,KKA, KKU, KKL )
ELSE
- PTP(:,:,:)= PBETA * MZF(ZFLXZ, KKA, KKU, KKL)
-END IF
+ IF (KRR /= 0) THEN
+ PTP(:,:,:) = PBETA * MZF( MZM(PETHETA,KKA, KKU, KKL) * ZFLXZ,KKA, KKU, KKL )
+ PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ ELSE
+ PTP(:,:,:)= PBETA * MZF( ZFLXZ,KKA, KKU, KKL )
+ END IF
+END IF
!
! Buoyancy flux at flux points
!
PWTHV = MZM(PETHETA, KKA, KKU, KKL) * 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
! Correction for qc and qi negative in AROME
-!IF ( KRRL >= 1 ) THEN
-! IF ( KRRI >= 1 ) THEN
-! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
-! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
-! *(1.0-PFRAC_ICE(:,:,:))
-! PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
-! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
-! *PFRAC_ICE(:,:,:)
-! ELSE
-! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
-! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL)
-! END IF
-!END IF
+IF(CPROGRAM/='AROME ') THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
+ *(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
+ *PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL)
+ END IF
+ END IF
+END IF
!
!* 2.4 Storage in LES configuration
!
@@ -626,6 +763,16 @@ IF (KRR /= 0) THEN
ZDFDDRDZ(:,:,:) = -XCSHF*ZKEFF*D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
ENDIF
!
+ ! Compute Leonard Terms for Cloud mixing ratio
+ IF (LHGRAD) THEN
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_LEONARD (:,:,:)= XCOEFHGRADRM*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX,XDZZ,XDZX,KKA,KKU,KKL)) &
+ *MZM(GX_M_M(PRM(:,:,:,1),XDXX,XDZZ,XDZX,KKA,KKU,KKL),KKA,KKU,KKL) &
+ +MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY,KKA,KKU,KKL)) &
+ *MZM(GY_M_M(PRM(:,:,:,1),XDYY,XDZZ,XDZY,KKA,KKU,KKL),KKA,KKU,KKL) )
+ END IF
+ !
! Effect of 3rd order terms in temperature flux (at flux point)
!
! d(w'2r')/dz
@@ -648,7 +795,7 @@ IF (KRR /= 0) THEN
!
! d(w'2th')/dz
IF (GFWTH) THEN
- ZF = ZF + M3_WR_W2TH(KKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,&
+ ZF = ZF + M3_WR_W2TH(KKA,KKU,KKL,PD,ZKEFF,&
& PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * PFWTH
ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
& PD,ZKEFF,PTKEM,PBLL_O_E,PETHETA) * PFWTH
@@ -656,7 +803,7 @@ IF (KRR /= 0) THEN
!
! d(w'th'2)/dz
IF (GFTH2) THEN
- ZF = ZF + M3_WR_WTH2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,PTKEM,&
+ ZF = ZF + M3_WR_WTH2(KKA,KKU,KKL,PD,ZKEFF,PTKEM,&
& PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDR_DZ) * MZM(PFTH2, KKA, KKU, KKL)
ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTH2_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
&ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA) * MZM(PFTH2, KKA, KKU, KKL)
@@ -672,28 +819,64 @@ IF (KRR /= 0) THEN
& PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, KKA, KKU, KKL)
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
!
- ! Compute the splitted conservative potential temperature at t+deltat
+ 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
+#ifdef REPRO48
+#else
+ ZF(:,:,IKE)=0.
+#endif
+ 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)
!
! Compute the equivalent tendency for the conservative mixing ratio
- PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRHODJ(:,:,:) * &
- (PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ZRWRNP (:,:,:) = PRHODJ(:,:,:)*(PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZRWRNP (:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:),KKA,KKU,KKL)*ZF_LEONARD(:,:,:),XDZZ,KKA,KKU,KKL)
+ END WHERE
+ END IF
+ !
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + ZRWRNP (:,:,:)
!
!* 3.2 Complete thermal production
!
@@ -702,6 +885,13 @@ IF (KRR /= 0) THEN
ZFLXZ(:,:,:) = ZF &
+ PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ
!
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZFLXZ(:,:,:) = ZF_LEONARD(:,:,:)
+ END WHERE
+ END IF
+ !
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
DO JK=IKTB+1,IKTE-1
@@ -728,31 +918,40 @@ IF (KRR /= 0) THEN
END IF
!
! Contribution of the conservative water flux to the Buoyancy flux
- ZA(:,:,:) = PBETA * MZF(MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
- ZA(:,:,IKB) = PBETA(:,:,IKB) * PEMOIST(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
- PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
+ IF (LOCEAN) THEN
+ ZA(:,:,:)= -XG*XBETAOC * MZF(ZFLXZ, KKA, KKU, KKL )
+ ELSE
+ ZA(:,:,:) = PBETA * MZF( MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ,KKA,KKU,KKL )
+ 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, KKA, KKU, KKL) * 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
! Correction of qc and qi negative for AROME
-! IF ( KRRL >= 1 ) THEN
-! IF ( KRRI >= 1 ) THEN
-! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
-! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
-! *(1.0-PFRAC_ICE(:,:,:))
-! PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
-! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
-! *PFRAC_ICE(:,:,:)
-! ELSE
-! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
-! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL)
-! END IF
-! END IF
+IF(CPROGRAM/='AROME ') THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL ) &
+ *(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL ) &
+ *PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL )
+ END IF
+ END IF
+END IF
!
!* 3.4 Storage in LES configuration
!
@@ -826,6 +1025,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
!
!----------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',1,ZHOOK_HANDLE)
diff --git a/src/mesonh/turb/mode_turb_ver_dyn_flux.f90 b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
index f103db1a7fe2ade1b33834179c41d1a0fe048c83..f698d06a1d3775d9443cc8bc669825ccf16a9037 100644
--- a/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
@@ -316,17 +316,12 @@ REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: &
ZFLXZ, & ! vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
-INTEGER :: ILENCH ! Length of comment string in LFIFM file
INTEGER :: IIB,IIE, & ! I index values for the Beginning and End
IJB,IJE, & ! mass points of the domain in the 3 direct.
IKB,IKE !
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV ! scalar loop counter
-CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
-CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM
! coefficients for the surface flux
@@ -340,15 +335,17 @@ TYPE(TFIELDDATA) :: TZFIELD
!
!* 1. PRELIMINARIES
! -------------
+ZA=XUNDEF
+PDP=XUNDEF
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
+!
IIU=SIZE(PUM,1)
-IIE=IIU-JPHEXT
-IIB=1+JPHEXT
IJU=SIZE(PUM,2)
-IJE=IJU-JPHEXT
-IJB=1+JPHEXT
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
IKT=SIZE(PUM,3)
@@ -360,9 +357,7 @@ IKTE=IKT-JPVEXT_TURB
ZSOURCE = 0.
ZFLXZ = 0.
ZCMFS = XCMFS
-IF (LHARAT)THEN
- ZCMFS=1.
-ENDIF
+IF (LHARAT) ZCMFS=1.
!
ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
! compute the coefficients for the uncentred gradient computation near the
@@ -371,9 +366,9 @@ ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
! With LHARATU length scale and TKE are at half levels so remove MZM
!
IF (LHARAT) THEN
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
@@ -394,7 +389,6 @@ ZA(:,:,:) = -PTSTEP * ZCMFS * &
MXM( ZKEFF ) * MXM(MZM(PRHODJ, KKA, KKU, KKL)) / &
MXM( PDZZ )**2
!
-IF (CPROGRAM/='AROME ') ZA(1,:,:)=ZA(IIE,:,:)
!
! Compute the source of U wind component
!
@@ -411,27 +405,50 @@ 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)) )
+ -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)) )
+ 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 splitted U at t+ deltat
+! Obtention of the split U at t+ deltat
!
CALL TRIDIAG_WIND(KKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MXM(PRHODJ),ZSOURCE,ZRES)
@@ -456,6 +473,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
@@ -485,7 +508,15 @@ PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
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
!
@@ -504,8 +535,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 &
@@ -535,6 +570,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
@@ -573,7 +623,6 @@ ZA(:,:,:) = - PTSTEP * ZCMFS * &
MYM( PDZZ )**2
!
!
-IF(CPROGRAM/='AROME ') ZA(:,1,:)=ZA(:,IJE,:)
!
! Compute the source of V wind component
! compute the coefficient between the vertical flux and the 2 components of the
@@ -589,26 +638,46 @@ 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 without coupling
+ ! 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 when coupling
+ ! input flux assumed to be in SI and at vorticity point
+ 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 splitted V at t+ deltat
+! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(KKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MYM(PRHODJ),ZSOURCE,ZRES)
!
@@ -632,6 +701,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'
@@ -663,6 +739,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
@@ -682,6 +766,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
@@ -709,9 +796,23 @@ IF(HTURBDIM=='3DIM') THEN
/(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB )) &
) &
* PDZY(:,:,IKB+KKL:IKB+KKL) &
- ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB))) &
+ ) / (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
@@ -737,15 +838,6 @@ IF(HTURBDIM=='3DIM') THEN
!
END IF
!
-! complete the dynamic production at the marginal points
-IF (CPROGRAM/='AROME ') THEN
- PDP(:,:,KKA)= -999.
- PDP(:,:,KKU)= -999.
- PDP(:,1,:)= PDP(:,IJE,:)
- PDP(:,IJE+1,:)= PDP(:,IJB,:)
- PDP(1,:,:)= PDP(IIE,:,:)
- PDP(IIE+1,:,:)= PDP(IIB,:,:)
-END IF
!
!----------------------------------------------------------------------------
!
diff --git a/src/mesonh/turb/mode_turb_ver_thermo_flux.f90 b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
index f96471652ea613d4522acbd4053c0d4c3a0ad748..78a188588c7593ee21fadffc2690fb9c34a41c36 100644
--- a/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
@@ -228,6 +228,7 @@ USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!
USE MODD_CST
+USE MODD_CONF, ONLY: CPROGRAM
USE MODD_CTURB
USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
USE MODD_GRID_n, ONLY: XZS, XXHAT, XYHAT
@@ -248,7 +249,7 @@ USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
USE MODI_GRADIENT_W
USE MODI_GRADIENT_M
-USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM
+USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM, MYF, MXF
USE MODI_TRIDIAG
USE MODI_LES_MEAN_SUBGRID
USE MODI_TRIDIAG_THERMO
@@ -364,19 +365,43 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
ZF, & ! Flux in dTh/dt =-dF/dz (evaluated at t-1)(or rt instead of Th)
ZDFDDTDZ, & ! dF/d(dTh/dz)
ZDFDDRDZ, & ! dF/d(dr/dz)
- Z3RDMOMENT ! 3 order term in flux or variance equation
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
-INTEGER :: ILENCH ! Length of comment string in LFIFM file
+ Z3RDMOMENT,& ! 3 order term in flux or variance equation
+ ZF_LEONARD,& ! Leonard terms
+ ZRWTHL, &
+ ZRWRNP, &
+ ZCLD_THOLD
+!
+REAL,DIMENSION(SIZE(XZS,1),SIZE(XZS,2),KKU) :: ZALT
+!
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
-CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
-CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+INTEGER :: JI, JJ ! loop indexes
!
-REAL :: ZTIME1, ZTIME2
+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)
+!
+!
+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
@@ -392,6 +417,38 @@ TYPE(TFIELDDATA) :: TZFIELD
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',0,ZHOOK_HANDLE)
+!
+! 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,IIU,IJU)
+ 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
IKTB =1+JPVEXT_TURB
@@ -405,11 +462,22 @@ GUSERV = (KRR/=0)
!
IF (LHARAT) THEN
! LHARAT so TKE and length scales at half levels!
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
+! Define a cloud mask with ri and rc (used after with a threshold) for Leonard terms
+!
+IF(LHGRAD) THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2) + PRM(:,:,:,4)
+ ELSE
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2)
+ END IF
+ END IF
+END IF
!
! Flags for 3rd order quantities
!
@@ -437,12 +505,22 @@ END IF
! Compute the turbulent flux F and F' at time t-dt.
!
IF (LHARAT) THEN
-ZF (:,:,:) = -ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
-ZDFDDTDZ(:,:,:) = -ZKEFF
+ ZF (:,:,:) = -ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ ZDFDDTDZ(:,:,:) = -ZKEFF
ELSE
-ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
-ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
-ENDIF
+ ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
+END IF
+!
+IF (LHGRAD) THEN
+ ! Compute the Leonard terms for thl
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_LEONARD (:,:,:)= XCOEFHGRADTHL*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX, XDZZ, XDZX))&
+ *MZM(GX_M_M(PTHLM(:,:,:),XDXX,XDZZ,XDZX)) &
+ + MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY)) &
+ *MZM(GY_M_M(PTHLM(:,:,:),XDYY,XDZZ,XDZY)) )
+END IF
!
! Effect of 3rd order terms in temperature flux (at flux point)
!
@@ -489,49 +567,94 @@ IF (GFTHR) THEN
ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTHR_O_DDTDZ(Z3RDMOMENT,PREDTH1,&
& PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR, KKA, KKU, KKL)
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
-!
-#ifdef REPRO55
-ZF(:,:,IKE)=0.
+ IF (LOCEAN) THEN
+ ZF(:,:,IKE) = XSSTFL(:,:) *0.5*(1. + PRHODJ(:,:,KKU) / PRHODJ(:,:,IKE))
+ ELSE !end ocean case (in nocoupled case)
+ ! atmos top
+#ifdef REPRO48
+#else
+ ZF(:,:,IKE)=0.
#endif
-! Compute the splitted conservative potential temperature at t+deltat
+ 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,&
PRHODJ,PTHLP)
!
! Compute the equivalent tendency for the conservative potential temperature
-PRTHLS(:,:,:)= PRTHLS(:,:,:) + &
- PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+!
+ZRWTHL(:,:,:)= PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZRWTHL(:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:))*ZF_LEONARD(:,:,:),XDZZ)
+ END WHERE
+END IF
+!
+PRTHLS(:,:,:)= PRTHLS(:,:,:) + ZRWTHL(:,:,:)
!
!* 2.2 Partial Thermal Production
!
! Conservative potential temperature flux :
!
ZFLXZ(:,:,:) = ZF &
- + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ
+ + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM) / PDZZ
+! replace the flux by the Leonard terms
+IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZFLXZ(:,:,:) = ZF_LEONARD(:,:,:)
+ END WHERE
+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))
+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(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
-
+END IF
+!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
! stores the conservative potential temperature vertical flux
TZFIELD%CMNHNAME = 'THW_FLX'
@@ -548,36 +671,44 @@ 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, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
- PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+IF (LOCEAN) THEN
+ PTP(:,:,:)= XG*XALPHAOC * MZF(ZFLXZ,KKA, KKU, KKL )
ELSE
- PTP(:,:,:)= PBETA * MZF(ZFLXZ, KKA, KKU, KKL)
-END IF
+ IF (KRR /= 0) THEN
+ PTP(:,:,:) = PBETA * MZF( MZM(PETHETA,KKA, KKU, KKL) * ZFLXZ,KKA, KKU, KKL )
+ PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ ELSE
+ PTP(:,:,:)= PBETA * MZF( ZFLXZ,KKA, KKU, KKL )
+ END IF
+END IF
!
! Buoyancy flux at flux points
!
PWTHV = MZM(PETHETA, KKA, KKU, KKL) * 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
-!
-#ifdef REPRO55
-IF ( KRRL >= 1 ) THEN
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
- PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
- *(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
- PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
- *PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
- PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL)
- END IF
+! Correction for qc and qi negative in AROME
+IF(CPROGRAM/='AROME ') THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
+ *(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL) &
+ *PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PATHETA*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL)
+ END IF
+ END IF
END IF
-#endif
!
!* 2.4 Storage in LES configuration
!
@@ -631,6 +762,16 @@ IF (KRR /= 0) THEN
ZDFDDRDZ(:,:,:) = -XCSHF*ZKEFF*D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
ENDIF
!
+ ! Compute Leonard Terms for Cloud mixing ratio
+ IF (LHGRAD) THEN
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_LEONARD (:,:,:)= XCOEFHGRADRM*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX, XDZZ, XDZX)) &
+ *MZM(GX_M_M(PRM(:,:,:,1),XDXX,XDZZ,XDZX)) &
+ +MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY)) &
+ *MZM(GY_M_M(PRM(:,:,:,1),XDYY,XDZZ,XDZY)) )
+ END IF
+ !
! Effect of 3rd order terms in temperature flux (at flux point)
!
! d(w'2r')/dz
@@ -677,31 +818,64 @@ IF (KRR /= 0) THEN
& PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, KKA, KKU, KKL)
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
!
-#ifdef REPRO55
- ZF(:,:,IKE)=0.
+ 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
+#ifdef REPRO48
+#else
+ ZF(:,:,IKE)=0.
#endif
- ! Compute the splitted conservative potential temperature at t+deltat
+ 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)
!
! Compute the equivalent tendency for the conservative mixing ratio
- PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRHODJ(:,:,:) * &
- (PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ZRWRNP (:,:,:) = PRHODJ(:,:,:)*(PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZRWRNP (:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:),KKA,KKU,KKL)*ZF_LEONARD(:,:,:),XDZZ,KKA,KKU,KKL)
+ END WHERE
+ END IF
+ !
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + ZRWRNP (:,:,:)
!
!* 3.2 Complete thermal production
!
@@ -710,6 +884,13 @@ IF (KRR /= 0) THEN
ZFLXZ(:,:,:) = ZF &
+ PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ
!
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZFLXZ(:,:,:) = ZF_LEONARD(:,:,:)
+ END WHERE
+ END IF
+ !
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
DO JK=IKTB+1,IKTE-1
@@ -736,33 +917,40 @@ IF (KRR /= 0) THEN
END IF
!
! Contribution of the conservative water flux to the Buoyancy flux
- ZA(:,:,:) = PBETA * MZF(MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
- ZA(:,:,IKB) = PBETA(:,:,IKB) * PEMOIST(:,:,IKB) * &
- 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
- PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
+ IF (LOCEAN) THEN
+ ZA(:,:,:)= -XG*XBETAOC * MZF(ZFLXZ, KKA, KKU, KKL )
+ ELSE
+ ZA(:,:,:) = PBETA * MZF( MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ,KKA,KKU,KKL )
+ 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, KKA, KKU, KKL) * 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
-!
-#ifdef REPRO55
- IF ( KRRL >= 1 ) THEN
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
- PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL ) &
- *(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
- PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL ) &
- *PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
- PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ, KKA, KKU, KKL )
- END IF
- END IF
-#endif
+! Correction of qc and qi negative for AROME
+IF(CPROGRAM/='AROME ') THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL ) &
+ *(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL ) &
+ *PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+ PRHODJ*PAMOIST*2.*PSRCM*DZF(ZFLXZ/PDZZ,KKA,KKU,KKL )
+ END IF
+ END IF
+END IF
!
!* 3.4 Storage in LES configuration
!
@@ -836,6 +1024,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
!
!----------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',1,ZHOOK_HANDLE)