diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index f78bb8fca8794c83669ebcac571fe7215e111e18..88f9822e93a37490ac63d4aa850c863631727cb3 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -10,7 +10,8 @@ CONTAINS
& PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, &
& PTSTEP,PIMPL,PEXPL, &
& HTURBLEN,HTURBDIM, &
- & TPFILE,OTURB_DIAG,OLES_CALL,ODIAG_IN_RUN, &
+ & TPFILE,OTURB_DIAG,OLES_CALL,ODIAG_IN_RUN,OOCEAN, &
+ & PSFUM,PSFVM, &
& PTP,PRTKES,PRTHLS,PCOEF_DISS,PTDIFF,PTDISS,PRTKEMS,&
& TBUDGETS, KBUDGETS, &
& PEDR, PTR,PDISS, PCURRENT_TKE_DISS )
@@ -183,6 +184,7 @@ TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
LOGICAL, INTENT(IN) :: OLES_CALL !
LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some
LOGICAL, INTENT(IN) :: ODIAG_IN_RUN ! switch to activate online diagnostics (mesonh)
+LOGICAL, INTENT(IN) :: OOCEAN ! switch for Ocean model version
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT):: PDP ! Dyn. prod. of TKE
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTRH
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTP ! Ther. prod. of TKE
@@ -199,6 +201,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT), OPTIONAL :: PTR ! Transp
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT), OPTIONAL :: PDISS ! Dissipation of TKE
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT), OPTIONAL :: PEDR ! EDR
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT), OPTIONAL :: PCURRENT_TKE_DISS ! if ODIAG_IN_RUN in mesonh
+REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFUM,PSFVM ! momentum sfc flux
!
!
!
@@ -266,11 +269,17 @@ END IF
!* 2.2 Explicit TKE sources except horizontal turbulent transport
!
! extrapolate the dynamic production with a 1/Z law from its value at the
-! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
-!
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-PDP(IIJB:IIJE,IKB) = PDP(IIJB:IIJE,IKB) * (1. + PDZZ(IIJB:IIJE,IKB+D%NKL)/PDZZ(IIJB:IIJE,IKB))
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
+IF (OOCEAN) THEN
+ ! W(IKE) value stored in PDP(IKE) to the mass localization of tke(IKE)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ PDP(IIJB:IIJE,IKE) = PDP(IIJB:IIJE,IKE) * (1. + PDZZ(IIJB:IIJE,IKE)/PDZZ(IIJB:IIJE,IKE+1))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ELSE
+ ! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ PDP(IIJB:IIJE,IKB) = PDP(IIJB:IIJE,IKB) * (1. + PDZZ(IIJB:IIJE,IKB+D%NKL)/PDZZ(IIJB:IIJE,IKB))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+END IF
!
! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
! + (Dynamical Production) + (Thermal Production) - (dissipation)
@@ -288,6 +297,14 @@ ZSOURCE(IIJB:IIJE,1:D%NKT) = ( PRTKES(IIJB:IIJE,1:D%NKT) + PRTKEMS(IIJB:IIJE,1:
!* 2.2 implicit vertical TKE transport
!
!
+! To add here in ZSOURCE surface flux of TKE
+!(assumed to be 0 for ATM,
+IF (OOCEAN) THEN
+ !for ocean:wave breaking simple/very rough param wE = 100 Ustar**3 where ustar is the Tau_atmi/rhocea
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE (IIJB:IIJE,IKE)=ZSOURCE(IIJB:IIJE,IKE)-1.E2*((PSFUM(IIJB:IIJE)**2 + PSFVM(IIJB:IIJE)**2)**1.5) /PDZZ(IIJB:IIJE,IKE)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+END IF
! Compute the vector giving the elements just under the diagonal for the
! matrix inverted in TRIDIAG
!
diff --git a/src/common/turb/mode_turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
index 76c03b560d84f4c4ecf19bd465929367a5a75629..5e34881e45b528f43b234bbc16ee55fdb516c93a 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -14,7 +14,7 @@ SUBROUTINE TURB_VER(D,CST,CSTURB,TURBN,KRR,KRRL,KRRI, &
PTSTEP, TPFILE, &
PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
PCOSSLOPE,PSINSLOPE, &
- PRHODJ,PTHVREF, &
+ PRHODJ,PTHVREF,PSFUM,PSFVM, &
PSFTHM,PSFRM,PSFSVM,PSFTHP,PSFRP,PSFSVP, &
PCDUEFF,PTAU11M,PTAU12M,PTAU33M, &
PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
@@ -288,6 +288,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: MFMOIST ! moist mass flux
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHVREF ! ref. state Virtual
! Potential Temperature
!
+REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFUM,PSFVM ! surface fluxes
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
REAL, DIMENSION(D%NIJT,KSV), INTENT(IN) :: PSFSVM ! t - deltat
!
@@ -560,7 +561,7 @@ CALL TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KSV,O2D,OFLAT, &
HTURBDIM,PIMPL,PEXPL,PTSTEP,TPFILE, &
PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
PCOSSLOPE,PSINSLOPE, &
- PRHODJ, &
+ PRHODJ,PSFUM,PSFVM, &
PCDUEFF,PTAU11M,PTAU12M,PTAU33M, &
PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
PTKEM,ZLM,MFMOIST,ZWU,ZWV, &
diff --git a/src/common/turb/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index fe103fa6b97870fb2c670de56cd746bce9fbd57c..c78e26457f7a0b6cf6ded482546ac5b410dc7d7b 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -6,13 +6,14 @@ MODULE MODE_TURB_VER_DYN_FLUX
IMPLICIT NONE
CONTAINS
SUBROUTINE TURB_VER_DYN_FLUX(D,CST,CSTURB,TURBN,KSV,O2D,OFLAT,&
- OTURB_FLX,KRR, OOCEAN,OHARAT,OCOUPLES,OLES_CALL,&
+ OTURB_FLX,KRR,OOCEAN,OHARAT,OCOUPLES,OLES_CALL,&
HTURBDIM,PIMPL,PEXPL, &
PTSTEP, &
TPFILE, &
PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
PCOSSLOPE,PSINSLOPE, &
PRHODJ, &
+ PSFUM,PSFVM, &
PCDUEFF,PTAU11M,PTAU12M,PTAU33M, &
PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
PTKEM,PLM,MFMOIST,PWU,PWV, &
@@ -271,6 +272,7 @@ REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODJ ! dry density * gri
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
!
+REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFUM,PSFVM ! normal momentum sfc flux
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PCDUEFF ! Cd * || u || at time t
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PTAU11M ! <uu> in the axes linked
! to the maximum slope direction and the surface normal and the binormal
@@ -331,7 +333,8 @@ INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV,JIJ,JK ! scalar loop counter
REAL, DIMENSION(D%NIJT) :: ZCOEFFLXU, &
- ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM
+ ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM, &
+ ZFLUXSFCU,ZFLUXSFCV
! coefficients for the surface flux
! evaluation and copy of PUSLOPEM and
! PVSLOPEM in local 3D arrays
@@ -382,6 +385,8 @@ ENDIF
!
ZUSLOPEM(IIJB:IIJE)=PUSLOPEM(IIJB:IIJE)
ZVSLOPEM(IIJB:IIJE)=PVSLOPEM(IIJB:IIJE)
+ZFLUXSFCU(IIJB:IIJE)=PSFUM(IIJB:IIJE)
+ZFLUXSFCV(IIJB:IIJE)=PSFVM(IIJB:IIJE)
!
!----------------------------------------------------------------------------
!
@@ -411,7 +416,7 @@ ZA(IIJB:IIJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIJB:IIJE,1:D%NKT)* ZWORK4(IIJB
ZCOEFFLXU(IIJB:IIJE) = PCDUEFF(IIJB:IIJE) * (PDIRCOSZW(IIJB:IIJE)**2 - ZDIRSINZW(IIJB:IIJE)**2) &
* PCOSSLOPE(IIJB:IIJE)
ZCOEFFLXV(IIJB:IIJE) = PCDUEFF(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE)
-
+!
! prepare the implicit scheme coefficients for the surface flux
ZCOEFS(IIJB:IIJE)= ZCOEFFLXU(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) &
+ZCOEFFLXV(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE)
@@ -422,60 +427,53 @@ ZWORK11D(IIJB:IIJE)=ZCOEFS(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB)
CALL MXM2D_PHY(D,ZWORK11D,ZCOEFS)
!
!
-! ZSOURCE= FLUX /DZ
+ZSOURCE(IIJB:IIJE,IKTB+1:IKTE-1) = 0.
+! ZSOURCE= sfc FLUX /DZ
+! Sfx flux assumed to be in SI & at vorticity point
CALL MXM_PHY(D,PRHODJ,ZWORK1)
-IF (OOCEAN) THEN ! OCEAN MODEL ONLY
- ! Sfx flux assumed to be in SI & at vorticity point
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- IF (OCOUPLES) THEN
- ZSOURCE(IIJB:IIJE,IKE) = PSSUFL_C(IIJB:IIJE)/PDZZ(IIJB:IIJE,IKE) &
- *0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE))
- ELSE
- ZSOURCE(IIJB:IIJE,IKE) = PSSUFL(IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKE) = ZSOURCE(IIJB:IIJE,IKE) /PDZZ(IIJB:IIJE,IKE) &
- *0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE) )
- ENDIF
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- !No flux at the ocean domain bottom
- ZSOURCE(IIJB:IIJE,IKB) = 0.
- ZSOURCE(IIJB:IIJE,IKTB+1:IKTE-1) = 0
!
-ELSE !ATMOS MODEL ONLY
- IF (OCOUPLES) THEN
+IF (OOCEAN) THEN ! Ocean model
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKB) = PSSUFL_C(IIJB:IIJE)/PDZZ(IIJB:IIJE,IKB) &
- * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
+ ZWORK21D(IIJB:IIJE) = ZFLUXSFCU(IIJB:IIJE)/PDZZ(IIJB:IIJE,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ELSE
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ! compute the explicit tangential flux at the W point
- ZSOURCE(IIJB:IIJE,IKB) = &
- PTAU11M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
- -PTAU12M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
- -PTAU33M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE)
+ CALL MXM2D_PHY(D,ZWORK21D,ZWORK31D)
+ !
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE(IIJB:IIJE,IKE) = ZWORK31D(IIJB:IIJE) &
+ *0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ !
+ ! Zero flux at the ocean domain bottom
+ ZSOURCE(IIJB:IIJE,IKB) = 0.
+ !
+ELSE ! Atmosphere
+ ! Compute the explicit tangential flux at the W point
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE(IIJB:IIJE,IKB) = &
+ PTAU11M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
+ -PTAU12M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
+ -PTAU33M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE)
!
- ! add the vertical part or the surface flux at the U,W vorticity point
+ ! add the vertical part or the surface flux at the U,W vorticity point
!
- ZWORK31D(IIJB:IIJE) = ZSOURCE(IIJB:IIJE,IKB)/PDZZ(IIJB:IIJE,IKB)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- CALL MXM2D_PHY(D,ZWORK31D,ZWORK41D)
- ZWORK51D(IIJB:IIJE)= ZCOEFFLXU(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
- *ZUSLOPEM(IIJB:IIJE) &
- -ZCOEFFLXV(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
- *ZVSLOPEM(IIJB:IIJE)
- CALL MXM2D_PHY(D,ZWORK51D,ZWORK61D)
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKB) = &
- ( ZWORK41D(IIJB:IIJE) &
- + ZWORK61D(IIJB:IIJE) &
- - ZCOEFS(IIJB:IIJE) * PUM(IIJB:IIJE,IKB) * PIMPL &
- ) * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ENDIF
+ ZWORK31D(IIJB:IIJE) = ZSOURCE(IIJB:IIJE,IKB)/PDZZ(IIJB:IIJE,IKB)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MXM2D_PHY(D,ZWORK31D,ZWORK41D)
+ ZWORK51D(IIJB:IIJE)= ZCOEFFLXU(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
+ *ZUSLOPEM(IIJB:IIJE) &
+ -ZCOEFFLXV(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
+ *ZVSLOPEM(IIJB:IIJE)
+ CALL MXM2D_PHY(D,ZWORK51D,ZWORK61D)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE(IIJB:IIJE,IKB) = &
+ ( ZWORK41D(IIJB:IIJE) &
+ + ZWORK61D(IIJB:IIJE) &
+ - ZCOEFS(IIJB:IIJE) * PUM(IIJB:IIJE,IKB) * PIMPL &
+ ) * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
- ZSOURCE(IIJB:IIJE,IKTB+1:IKTE-1) = 0.
ZSOURCE(IIJB:IIJE,IKE) = 0.
-ENDIF !end ocean or atmosphere cases
+ENDIF
!
! Obtention of the split U at t+ deltat
!
@@ -495,8 +493,7 @@ CALL MXM_PHY(D,PDZZ,ZWORK5)
PRUS(IIJB:IIJE,1:D%NKT)= PRUS(IIJB:IIJE,1:D%NKT)+ZWORK1(IIJB:IIJE,1:D%NKT)*(ZRES(IIJB:IIJE,1:D%NKT) &
- PUM(IIJB:IIJE,1:D%NKT))/PTSTEP
!
-!
-!* 5.2 Partial Dynamic Production
+!* 5.2 Partial TKE Dynamic Production
!
! vertical flux of the U wind component
!
@@ -504,26 +501,20 @@ ZFLXZ(IIJB:IIJE,1:D%NKT) = -ZCMFS * ZWORK2(IIJB:IIJE,1:D%NKT) * ZWORK4(IIJB:
/ ZWORK5(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
-! surface flux
-CALL MXM_PHY(D,PDZZ,ZWORK1)
-CALL MXM_PHY(D,PRHODJ,ZWORK2)
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-ZFLXZ(IIJB:IIJE,IKB) = ZWORK1(IIJB:IIJE,IKB) * &
- ( ZSOURCE(IIJB:IIJE,IKB) &
- +ZCOEFS(IIJB:IIJE) * ZRES(IIJB:IIJE,IKB) * PIMPL &
- ) / 0.5 / ( 1. + ZWORK2(IIJB:IIJE,D%NKA)/ ZWORK2(IIJB:IIJE,IKB) )
-!
-ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
-!
-IF (OOCEAN) THEN !ocean model at phys sfc (ocean domain top)
-!
-!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZFLXZ(IIJB:IIJE,IKE) = ZWORK1(IIJB:IIJE,IKE) * &
- ZSOURCE(IIJB:IIJE,IKE) &
- / 0.5 / ( 1. + ZWORK2(IIJB:IIJE,D%NKU)/ ZWORK2(IIJB:IIJE,IKE) )
- ZFLXZ(IIJB:IIJE,D%NKU) = ZFLXZ(IIJB:IIJE,IKE)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
+IF (OOCEAN) THEN
+ ZFLXZ(IIJB:IIJE,IKE+1) = ZFLUXSFCU(IIJB:IIJE)
+ELSE
+ ! surface flux
+ CALL MXM_PHY(D,PDZZ,ZWORK1)
+ CALL MXM_PHY(D,PRHODJ,ZWORK2)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZFLXZ(IIJB:IIJE,IKB) = ZWORK1(IIJB:IIJE,IKB) * &
+ ( ZSOURCE(IIJB:IIJE,IKB) &
+ +ZCOEFS(IIJB:IIJE) * ZRES(IIJB:IIJE,IKB) * PIMPL &
+ ) / 0.5 / ( 1. + ZWORK2(IIJB:IIJE,D%NKA)/ ZWORK2(IIJB:IIJE,IKB) )
+ !
+ ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -545,8 +536,8 @@ END IF
!
PWU(IIJB:IIJE,1:D%NKT) = ZFLXZ(IIJB:IIJE,1:D%NKT)
!
-! Contribution to the dynamic production of TKE
-! compute the dynamic production at the mass point
+! Contribution to the TKE dynamic production of TKE
+! (computed at mass point)
!
CALL GZ_U_UW_PHY(D,PUM,PDZZ,ZWORK1)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
@@ -558,23 +549,30 @@ CALL MZF_PHY(D,ZWORK3,ZWORK4)
PDP(IIJB:IIJE,1:D%NKT) = -ZWORK4(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
-! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
+! Special cases near surface
CALL MXM_PHY(D,PDZZ,ZWORK1)
-ZWORK2(IIJB:IIJE,IKB) = ZFLXZ(IIJB:IIJE,IKB+D%NKL) * (PUM(IIJB:IIJE,IKB+D%NKL)-PUM(IIJB:IIJE,IKB)) &
- / ZWORK1(IIJB:IIJE,IKB+D%NKL)
-CALL MXF_PHY(D,ZWORK2,ZWORK3)
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-PDP(IIJB:IIJE,IKB) = -ZWORK3(IIJB:IIJE,IKB)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
-!
IF (OOCEAN) THEN
- ZWORK2(IIJB:IIJE,IKE) = ZFLXZ(IIJB:IIJE,IKE-D%NKL) * (PUM(IIJB:IIJE,IKE)-PUM(IIJB:IIJE,IKE-D%NKL)) &
+ ! evaluate the dynamic production at w(IKE) and store in PDP(IKE)
+ ! before to be extrapolated in tke_eps routine
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK2(IIJB:IIJE,IKE) = ZFLXZ(IIJB:IIJE,IKE) * (PUM(IIJB:IIJE,IKE)-PUM(IIJB:IIJE,IKE-D%NKL)) &
/ ZWORK1(IIJB:IIJE,IKE-D%NKL)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL MXF_PHY(D,ZWORK2,ZWORK3)
- ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
- !$mnh_expand_array(JIJ=IIJB:IIJE)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
PDP(IIJB:IIJE,IKE) = -ZWORK3(IIJB:IIJE,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ELSE ! Atmosphere
+ ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK2(IIJB:IIJE,IKB) = ZFLXZ(IIJB:IIJE,IKB+D%NKL) * (PUM(IIJB:IIJE,IKB+D%NKL)-PUM(IIJB:IIJE,IKB)) &
+ / ZWORK1(IIJB:IIJE,IKB+D%NKL)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MXF_PHY(D,ZWORK2,ZWORK3)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ PDP(IIJB:IIJE,IKB) = -ZWORK3(IIJB:IIJE,IKB)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+!
END IF
!
! Storage in the LES configuration
@@ -654,7 +652,7 @@ IF(HTURBDIM=='3DIM') THEN
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
END IF
!
- ! Complete the Dynamical production with the W wind component
+ ! Complete the TKE dynamical production with the W wind contribution
!
CALL GX_W_UW_PHY(D,OFLAT,PWM,PDXX,PDZZ,PDZX, ZWORK1)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
@@ -666,48 +664,39 @@ IF(HTURBDIM=='3DIM') THEN
ZA(IIJB:IIJE,1:D%NKT) = -ZWORK3(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
+ ! Special cases near surface
+ CALL DXM_PHY(D,PWM,ZWORK1)
+ IF (OOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE) in PDP(IKE)
+ !
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKE) * ZWORK1(IIJB:IIJE,IKE) &
+ / (0.5*(PDXX(IIJB:IIJE,IKE-D%NKL)+PDXX(IIJB:IIJE,IKE)))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MXF2D_PHY(D,ZWORK31D,ZWORK41D)
+ ZA(IIJB:IIJE,IKE) = ZWORK41D(IIJB:IIJE)
!
- ! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
-
- CALL DXM_PHY(D,PWM,ZWORK1)
- !
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKB+2*D%NKL)-PWM(IIJB:IIJE,IKB+D%NKL)) &
- / (PDZZ(IIJB:IIJE,IKB+2*D%NKL)+PDZZ(IIJB:IIJE,IKB+D%NKL)) &
- + (PWM(IIJB:IIJE,IKB+D%NKL)-PWM(IIJB:IIJE,IKB)) &
- / (PDZZ(IIJB:IIJE,IKB+D%NKL)+PDZZ(IIJB:IIJE,IKB))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- !
- CALL MXM2D_PHY(D,ZWORK21D,ZWORK51D)
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKB+D%NKL) &
- * ( ZWORK1(IIJB:IIJE,IKB+D%NKL) - ZWORK51D(IIJB:IIJE) &
- * PDZX(IIJB:IIJE,IKB+D%NKL) ) &
- / (0.5*(PDXX(IIJB:IIJE,IKB+D%NKL)+PDXX(IIJB:IIJE,IKB)))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- CALL MXF2D_PHY(D,ZWORK31D,ZWORK41D)
- ZA(IIJB:IIJE,IKB) = ZWORK41D(IIJB:IIJE)
- !
-IF (OOCEAN) THEN
- ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
- !
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKE-2*D%NKL)-PWM(IIJB:IIJE,IKE-D%NKL)) &
- / (PDZZ(IIJB:IIJE,IKE-2*D%NKL)+PDZZ(IIJB:IIJE,IKE-D%NKL)) &
- + (PWM(IIJB:IIJE,IKE-D%NKL)-PWM(IIJB:IIJE,IKE )) &
- / (PDZZ(IIJB:IIJE,IKE-D%NKL)+PDZZ(IIJB:IIJE,IKE ))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- !
- CALL MXM2D_PHY(D,ZWORK21D,ZWORK51D)
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKE-D%NKL) &
- * ( ZWORK1(IIJB:IIJE,IKE-D%NKL) - ZWORK51D(IIJB:IIJE) &
- * PDZX(IIJB:IIJE,IKE-D%NKL) ) &
- / (0.5*(PDXX(IIJB:IIJE,IKE-D%NKL)+PDXX(IIJB:IIJE,IKE)))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- CALL MXF2D_PHY(D,ZWORK31D,ZWORK41D)
- ZA(IIJB:IIJE,IKE) = ZWORK41D(IIJB:IIJE)
-END IF
+ ELSE !Atmosphere
+ ! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
+ !
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKB+2*D%NKL)-PWM(IIJB:IIJE,IKB+D%NKL)) &
+ / (PDZZ(IIJB:IIJE,IKB+2*D%NKL)+PDZZ(IIJB:IIJE,IKB+D%NKL)) &
+ + (PWM(IIJB:IIJE,IKB+D%NKL)-PWM(IIJB:IIJE,IKB)) &
+ / (PDZZ(IIJB:IIJE,IKB+D%NKL)+PDZZ(IIJB:IIJE,IKB))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ !
+ CALL MXM2D_PHY(D,ZWORK21D,ZWORK51D)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKB+D%NKL) &
+ * ( ZWORK1(IIJB:IIJE,IKB+D%NKL) - ZWORK51D(IIJB:IIJE) &
+ * PDZX(IIJB:IIJE,IKB+D%NKL) ) &
+ / (0.5*(PDXX(IIJB:IIJE,IKB+D%NKL)+PDXX(IIJB:IIJE,IKB)))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MXF2D_PHY(D,ZWORK31D,ZWORK41D)
+ ZA(IIJB:IIJE,IKB) = ZWORK41D(IIJB:IIJE)
+ !
+ END IF
!
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
PDP(IIJB:IIJE,1:D%NKT)=PDP(IIJB:IIJE,1:D%NKT)+ZA(IIJB:IIJE,1:D%NKT)
@@ -760,7 +749,7 @@ END IF
!----------------------------------------------------------------------------
!
!
-!* 6. SOURCES OF V,W WIND COMPONENTS AND COMPLETE 1D DYNAMIC PRODUCTION
+!* 6. SOURCES OF V,W WIND COMPONENTS AND COMPLETE 1D TKE DYNAMIC PRODUCTION
! -----------------------------------------------------------------
!
!* 6.1 Source of V wind component
@@ -797,63 +786,63 @@ ZWORK11D(IIJB:IIJE)=ZCOEFS(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL MYM2D_PHY(D,ZWORK11D,ZCOEFS)
!
-CALL MYM_PHY(D,PRHODJ,ZWORK1) ! it was MXM(PRHODJ) : bug corrected now ? REPRO55 OCEAN
+! No flux in SOURCE TERM NULL OUTSIDE BC
+ZSOURCE(IIJB:IIJE,IKB+1:IKE-1) = 0.
+! Surface case
+CALL MYM_PHY(D,PRHODJ,ZWORK1)
IF (OOCEAN) THEN ! Ocean case
- IF (OCOUPLES) THEN
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKE) = PSSVFL_C(IIJB:IIJE)/PDZZ(IIJB:IIJE,IKE) &
- *0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ELSE
+ ZCOEFFLXU(IIJB:IIJE) = PCDUEFF(IIJB:IIJE)
+ ZCOEFFLXV(IIJB:IIJE) = PCDUEFF(IIJB:IIJE)
+ ZCOEFS(IIJB:IIJE)=ZCOEFFLXU(IIJB:IIJE)
+ ! average this flux to be located at the U,W vorticity point
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK11D(IIJB:IIJE) = ZCOEFS(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKE)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MYM2D_PHY(D,ZWORK11D,ZCOEFS)
+ !
+ ZWORK11D(IIJB:IIJE) = ZFLUXSFCV(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKE)
+ CALL MYM2D_PHY(D,ZWORK11D,ZWORK21D)
+ !
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKE) = PSSVFL(IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKE) = ZSOURCE(IIJB:IIJE,IKE)/PDZZ(IIJB:IIJE,IKE) &
+ ZSOURCE(IIJB:IIJE,IKE) = ZWORK21D(IIJB:IIJE) &
*0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE))
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
- END IF
!No flux at the ocean domain bottom
ZSOURCE(IIJB:IIJE,IKB) = 0.
+!
ELSE ! Atmos case
+!
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK31D(IIJB:IIJE) = ZCOEFFLXU(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
- *ZUSLOPEM(IIJB:IIJE) &
- +ZCOEFFLXV(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
+ ZWORK31D(IIJB:IIJE) = ZCOEFFLXU(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
+ *ZUSLOPEM(IIJB:IIJE) &
+ +ZCOEFFLXV(IIJB:IIJE) / PDZZ(IIJB:IIJE,IKB) &
*ZVSLOPEM(IIJB:IIJE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL MYM2D_PHY(D,ZWORK31D,ZWORK61D)
-!
- IF (.NOT.OCOUPLES) THEN ! only atmosp without coupling
+ !
! compute the explicit tangential flux at the W point
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKB) = &
- PTAU11M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
- +PTAU12M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
- -PTAU33M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE)
- !
- ZWORK31D(IIJB:IIJE) = ZSOURCE(IIJB:IIJE,IKB)/PDZZ(IIJB:IIJE,IKB)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- CALL MYM2D_PHY(D,ZWORK31D,ZWORK51D)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE(IIJB:IIJE,IKB) = &
+ PTAU11M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
+ +PTAU12M(IIJB:IIJE) * PCOSSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) &
+ -PTAU33M(IIJB:IIJE) * PSINSLOPE(IIJB:IIJE) * ZDIRSINZW(IIJB:IIJE) * PDIRCOSZW(IIJB:IIJE)
+ !
+ ZWORK31D(IIJB:IIJE) = ZSOURCE(IIJB:IIJE,IKB)/PDZZ(IIJB:IIJE,IKB)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MYM2D_PHY(D,ZWORK31D,ZWORK51D)
!
! add the vertical part or the surface flux at the V,W vorticity point
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKB) = &
- ( ZWORK51D(IIJB:IIJE) &
- + ZWORK61D(IIJB:IIJE) &
- - ZCOEFS(IIJB:IIJE) * PVM(IIJB:IIJE,IKB) * PIMPL &
- ) * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZSOURCE(IIJB:IIJE,IKB) = &
+ ( ZWORK51D(IIJB:IIJE) &
+ + ZWORK61D(IIJB:IIJE) &
+ - ZCOEFS(IIJB:IIJE) * PVM(IIJB:IIJE,IKB) * PIMPL &
+ ) * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
- ELSE !atmosphere when coupling
- ! input flux assumed to be in SI and at vorticity point
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZSOURCE(IIJB:IIJE,IKB) = -PSSVFL_C(IIJB:IIJE)/(1.*PDZZ(IIJB:IIJE,IKB)) &
- * 0.5 * ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ENDIF
!No flux at the atmosphere top
ZSOURCE(IIJB:IIJE,IKE) = 0.
ENDIF ! End of Ocean or Atmospher Cases
-ZSOURCE(IIJB:IIJE,IKTB+1:IKTE-1) = 0.
!
! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(D,PVM,ZA,ZCOEFS,PTSTEP,PEXPL,PIMPL, &
@@ -881,22 +870,16 @@ ZFLXZ(IIJB:IIJE,1:D%NKT) = -ZCMFS * ZWORK2(IIJB:IIJE,1:D%NKT) * ZWORK4(IIJB:II
/ ZWORK5(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-ZFLXZ(IIJB:IIJE,IKB) = ZWORK5(IIJB:IIJE,IKB) * &
- ( ZSOURCE(IIJB:IIJE,IKB) &
- +ZCOEFS(IIJB:IIJE) * ZRES(IIJB:IIJE,IKB) * PIMPL &
- ) / 0.5 / ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
-!
-!
-ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
-!
IF (OOCEAN) THEN
+ ZFLXZ(IIJB:IIJE,IKE+1) = ZFLUXSFCV(IIJB:IIJE)
+ELSE
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZFLXZ(IIJB:IIJE,IKE) = ZWORK5(IIJB:IIJE,IKE) * &
- ZSOURCE(IIJB:IIJE,IKE) &
- / 0.5 / ( 1. + ZWORK1(IIJB:IIJE,D%NKU) / ZWORK1(IIJB:IIJE,IKE) )
- ZFLXZ(IIJB:IIJE,D%NKU) = ZFLXZ(IIJB:IIJE,IKE)
+ ZFLXZ(IIJB:IIJE,IKB) = ZWORK5(IIJB:IIJE,IKB) * &
+ ( ZSOURCE(IIJB:IIJE,IKB) &
+ +ZCOEFS(IIJB:IIJE) * ZRES(IIJB:IIJE,IKB) * PIMPL &
+ ) / 0.5 / ( 1. + ZWORK1(IIJB:IIJE,D%NKA) / ZWORK1(IIJB:IIJE,IKB) )
+ !
+ ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
@@ -919,8 +902,8 @@ END IF
!
PWV(IIJB:IIJE,1:D%NKT) = ZFLXZ(IIJB:IIJE,1:D%NKT)
!
-! Contribution to the dynamic production of TKE
-! compute the dynamic production contribution at the mass point
+! Contribution to the TKE dynamical production
+! computed at mass point
!
CALL GZ_V_VW_PHY(D,PVM,PDZZ,ZWORK1)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
@@ -932,27 +915,30 @@ CALL MZF_PHY(D,ZWORK3,ZWORK4)
ZA(IIJB:IIJE,1:D%NKT) = -ZWORK4(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
-! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
+! Special cases at surface
CALL MYM_PHY(D,PDZZ,ZWORK1)
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-ZWORK2(IIJB:IIJE,IKB) = ZFLXZ(IIJB:IIJE,IKB+D%NKL) * (PVM(IIJB:IIJE,IKB+D%NKL)-PVM(IIJB:IIJE,IKB)) &
- / ZWORK1(IIJB:IIJE,IKB+D%NKL)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
-CALL MYF_PHY(D,ZWORK2,ZWORK3)
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-ZA(IIJB:IIJE,IKB) = -ZWORK3(IIJB:IIJE,IKB)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
-!
IF (OOCEAN) THEN
- ! evaluate the dynamic production at w(IKE-D%NKL) in PDP(IKE)
+ ! evaluate the dynamic production at w(IKE) in PDP(IKE)
+ ! before extrapolation done in routine tke_eps_source
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK2(IIJB:IIJE,IKE) = ZFLXZ(IIJB:IIJE,IKE-D%NKL) * (PVM(IIJB:IIJE,IKE)-PVM(IIJB:IIJE,IKE-D%NKL)) &
+ ZWORK2(IIJB:IIJE,IKE) = ZFLXZ(IIJB:IIJE,IKE) * (PVM(IIJB:IIJE,IKE)-PVM(IIJB:IIJE,IKE-D%NKL)) &
/ ZWORK1(IIJB:IIJE,IKE-D%NKL)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL MYF_PHY(D,ZWORK2,ZWORK3)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZA(IIJB:IIJE,IKE) = -ZWORK3(IIJB:IIJE,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
+!
+ELSE ! Atmosphere
+ ! evaluate the dynamic production at w(IKB+D%NKL) in PDP(IKB)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK2(IIJB:IIJE,IKB) = ZFLXZ(IIJB:IIJE,IKB+D%NKL) * (PVM(IIJB:IIJE,IKB+D%NKL)-PVM(IIJB:IIJE,IKB)) &
+ / ZWORK1(IIJB:IIJE,IKB+D%NKL)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MYF_PHY(D,ZWORK2,ZWORK3)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZA(IIJB:IIJE,IKB) = -ZWORK3(IIJB:IIJE,IKB)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
@@ -985,12 +971,13 @@ END IF
!
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZFLXZ(IIJB:IIJE,D%NKA) = 2 * ZFLXZ(IIJB:IIJE,IKB) - ZFLXZ(IIJB:IIJE,IKB+D%NKL) ! extrapolation
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
IF (OOCEAN) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZFLXZ(IIJB:IIJE,D%NKU) = 2 * ZFLXZ(IIJB:IIJE,IKE) - ZFLXZ(IIJB:IIJE,IKE-D%NKL) ! extrapolation
+ ZFLXZ(IIJB:IIJE,IKE+D%NKL) = 2 * ZFLXZ(IIJB:IIJE,IKE) - ZFLXZ(IIJB:IIJE,IKE-D%NKL) ! extrapolation
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ ELSE
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZFLXZ(IIJB:IIJE,D%NKA) = 2 * ZFLXZ(IIJB:IIJE,IKB) - ZFLXZ(IIJB:IIJE,IKB+D%NKL) ! extrapolation
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
@@ -1045,43 +1032,36 @@ IF(HTURBDIM=='3DIM') THEN
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
ZA(IIJB:IIJE,1:D%NKT) = -ZWORK3(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
- !
- CALL DYM_PHY(D,PWM,ZWORK1)
- !
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKB+2*D%NKL )-PWM(IIJB:IIJE,IKB+D%NKL)) &
- / (PDZZ(IIJB:IIJE,IKB+2*D%NKL)+PDZZ(IIJB:IIJE,IKB+D%NKL)) &
- + (PWM(IIJB:IIJE,IKB+D%NKL)-PWM(IIJB:IIJE,IKB)) &
- / (PDZZ(IIJB:IIJE,IKB+D%NKL)+PDZZ(IIJB:IIJE,IKB))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- !
- CALL MYM2D_PHY(D,ZWORK21D,ZWORK51D)
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK31D(IIJB:IIJE ) = - ZFLXZ(IIJB:IIJE,IKB+D%NKL) &
- * ( ZWORK1(IIJB:IIJE,IKB+D%NKL) - ZWORK51D(IIJB:IIJE ) &
- * PDZY(IIJB:IIJE,IKB+D%NKL) ) &
- / (0.5*(PDYY(IIJB:IIJE,IKB+D%NKL)+PDYY(IIJB:IIJE,IKB)))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- CALL MYF2D_PHY(D,ZWORK31D,ZWORK41D)
- ZA(IIJB:IIJE,IKB) = ZWORK41D(IIJB:IIJE)
- !
+ !
+ CALL DYM_PHY(D,PWM,ZWORK1)
+ ! Special case near surface
IF (OOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE) and stored in PDP(IKE)
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKE) * ZWORK1(IIJB:IIJE,IKE) &
+ / (0.5*(PDYY(IIJB:IIJE,IKE-D%NKL)+PDYY(IIJB:IIJE,IKE)))
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ CALL MYF2D_PHY(D,ZWORK31D,ZWORK41D)
+ ZA(IIJB:IIJE,IKE) = ZWORK41D(IIJB:IIJE)
+ ELSE ! Atmosphere
+ ! evaluate the dynamic production at w(IKB+KKL) and stored in PDP(IKB)
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKE-2*D%NKL)-PWM(IIJB:IIJE,IKE-D%NKL)) &
- / (PDZZ(IIJB:IIJE,IKE-2*D%NKL)+PDZZ(IIJB:IIJE,IKE-D%NKL)) &
- + (PWM(IIJB:IIJE,IKE-D%NKL)-PWM(IIJB:IIJE,IKE )) &
- / (PDZZ(IIJB:IIJE,IKE-D%NKL)+PDZZ(IIJB:IIJE,IKE ))
+ ZWORK21D(IIJB:IIJE) = (PWM(IIJB:IIJE,IKB+2*D%NKL )-PWM(IIJB:IIJE,IKB+D%NKL)) &
+ / (PDZZ(IIJB:IIJE,IKB+2*D%NKL)+PDZZ(IIJB:IIJE,IKB+D%NKL)) &
+ + (PWM(IIJB:IIJE,IKB+D%NKL)-PWM(IIJB:IIJE,IKB)) &
+ / (PDZZ(IIJB:IIJE,IKB+D%NKL)+PDZZ(IIJB:IIJE,IKB))
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
CALL MYM2D_PHY(D,ZWORK21D,ZWORK51D)
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZWORK31D(IIJB:IIJE) = - ZFLXZ(IIJB:IIJE,IKE-D%NKL) &
- * ( ZWORK1(IIJB:IIJE,IKE-D%NKL) - ZWORK51D(IIJB:IIJE ) &
- * PDZY(IIJB:IIJE,IKE-D%NKL) ) &
- / (0.5*(PDYY(IIJB:IIJE,IKE-D%NKL)+PDYY(IIJB:IIJE,IKE)))
+ ZWORK31D(IIJB:IIJE ) = - ZFLXZ(IIJB:IIJE,IKB+D%NKL) &
+ * ( ZWORK1(IIJB:IIJE,IKB+D%NKL) - ZWORK51D(IIJB:IIJE ) &
+ * PDZY(IIJB:IIJE,IKB+D%NKL) ) &
+ / (0.5*(PDYY(IIJB:IIJE,IKB+D%NKL)+PDYY(IIJB:IIJE,IKB)))
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL MYF2D_PHY(D,ZWORK31D,ZWORK41D)
- ZA(IIJB:IIJE,IKE) = ZWORK41D(IIJB:IIJE)
+ ZA(IIJB:IIJE,IKB) = ZWORK41D(IIJB:IIJE)
+ !
END IF
!
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index c15331683b07680354725c1f51bf21a39a3a9d40..9398968be719efaa3f7fd4f1c81a75047c393138 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -573,22 +573,13 @@ IF (GFTHR) THEN
* ZWORK2(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
END IF
-! compute interface flux
-IF (OCOUPLES) THEN ! Autocoupling O-A LES
- IF (OOCEAN) THEN ! ocean model in coupled case
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZF(IIJB:IIJE,IKE) = (PSSTFL_C(IIJB:IIJE)+PSSRFL_C(IIJB:IIJE)) &
- *0.5* ( 1. + PRHODJ(IIJB:IIJE,D%NKU)/PRHODJ(IIJB:IIJE,IKE) )
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ELSE ! atmosph model in coupled case
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZF(IIJB:IIJE,IKB) = PSSTFL_C(IIJB:IIJE) &
- *0.5* ( 1. + PRHODJ(IIJB:IIJE,D%NKA)/PRHODJ(IIJB:IIJE,IKB) )
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ENDIF
-!
-ELSE ! No coupling O and A cases
- ! atmosp bottom
+! specialcase for surface
+IF (OOCEAN) THEN ! ocean model in coupled case
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZF(IIJB:IIJE,IKE+1) = PSFTHM(IIJB:IIJE) &
+ *0.5* ( 1. + PRHODJ(IIJB:IIJE,D%NKU)/PRHODJ(IIJB:IIJE,IKE) )
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ELSE ! 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
@@ -606,18 +597,14 @@ ELSE ! No coupling O and A cases
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
- IF (OOCEAN) THEN
- !$mnh_expand_array(JIJ=IIJB:IIJE)
- ZF(IIJB:IIJE,IKE) = PSSTFL(IIJB:IIJE) *0.5*(1. + PRHODJ(IIJB:IIJE,D%NKU) / PRHODJ(IIJB:IIJE,IKE))
- !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- ELSE !end ocean case (in nocoupled case)
! atmos top
#ifdef REPRO48
#else
ZF(IIJB:IIJE,IKE)=0.
+ !TODO merge : the following solution must be kept :
+ !ZF((IIJB:IIJE,IKE+1)=0.
#endif
- END IF
-END IF !end no coupled cases
+END IF
!
! Compute the split conservative potential temperature at t+deltat
CALL TRIDIAG_THERMO(D,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
@@ -674,12 +661,13 @@ IF (TURBN%LHGRAD) THEN
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
END IF
!
-!$mnh_expand_array(JIJ=IIJB:IIJE)
-ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
IF (OOCEAN) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE)
- ZFLXZ(IIJB:IIJE,D%NKU) = ZFLXZ(IIJB:IIJE,IKE)
+ ZFLXZ(IIJB:IIJE,IKE+1) = ZFLXZ(IIJB:IIJE,IKE)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ELSE
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!
@@ -696,8 +684,8 @@ PWTH(IIJB:IIJE,IKB)=0.5*(ZFLXZ(IIJB:IIJE,IKB)+ZFLXZ(IIJB:IIJE,IKB+D%NKL))
IF (OOCEAN) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE)
PWTH(IIJB:IIJE,IKE)=0.5*(ZFLXZ(IIJB:IIJE,IKE)+ZFLXZ(IIJB:IIJE,IKE+D%NKL))
- PWTH(IIJB:IIJE,D%NKA)=0.
- PWTH(IIJB:IIJE,D%NKU)=ZFLXZ(IIJB:IIJE,D%NKU)
+ PWTH(IIJB:IIJE,D%NKA)=0.
+ PWTH(IIJB:IIJE,D%NKU)=PWTH(IIJB:IIJE,IKE)! not used
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
ELSE
!$mnh_expand_array(JIJ=IIJB:IIJE)
@@ -979,17 +967,12 @@ IF (KRR /= 0) THEN
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
END IF
!
- ! compute interface flux
- IF (OCOUPLES) THEN ! coupling NH O-A
- IF (OOCEAN) THEN ! ocean model in coupled case
- ! evap effect on salinity to be added later !!!
- ZF(IIJB:IIJE,IKE) = 0.
- ELSE ! atmosph model in coupled case
- ZF(IIJB:IIJE,IKB) = 0.
- ! AJOUTER FLUX EVAP SUR MODELE ATMOS
- ENDIF
- !
- ELSE ! No coupling NH OA case
+ !special case at sfc
+ IF (OOCEAN) THEN
+ ! General ocean case
+ ! salinity/evap effect to be added later !!!!!
+ ZF(IIJB:IIJE,IKE) = 0.
+ ELSE ! atmosp case
! atmosp bottom
!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
! (in presence of slopes)
@@ -1009,19 +992,12 @@ IF (KRR /= 0) THEN
* 0.5 * (1. + PRHODJ(IIJB:IIJE,D%NKA) / PRHODJ(IIJB:IIJE,IKB))
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
- !
- IF (OOCEAN) THEN
- ! General ocean case
- ! salinity/evap effect to be added later !!!!!
- ZF(IIJB:IIJE,IKE) = 0.
- ELSE !end ocean case (in nocoupled case)
! atmos top
#ifdef REPRO48
#else
ZF(IIJB:IIJE,IKE)=0.
#endif
END IF
- END IF!end no coupled cases
! Compute the split conservative potential temperature at t+deltat
CALL TRIDIAG_THERMO(D,PRM(:,:,1),ZF,ZDFDDRDZ,PTSTEP,PIMPL,&
PDZZ,PRHODJ,PRP)
@@ -1035,11 +1011,11 @@ IF (KRR /= 0) THEN
!
! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
IF (TURBN%LHGRAD) THEN
- DO JK=1,D%NKU
+! DO JK=1,D%NKU
! !$mnh_expand_array(JIJ=IIJB:IIJE)
! ZALT(IIJB:IIJE,JK) = PZZ(IIJB:IIJE,JK)-XZS(IIJB:IIJE) !TODO TO BE ADDED AS ARGS OF TURB
! !$mnh_end_expand_array(JIJ=IIJB:IIJE)
- END DO
+! END DO
CALL MZM_PHY(D,PRHODJ,ZWORK1)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
ZWORK2(IIJB:IIJE,1:D%NKT) = ZWORK1(IIJB:IIJE,1:D%NKT)*ZF_LEONARD(IIJB:IIJE,1:D%NKT)
@@ -1080,6 +1056,10 @@ IF (KRR /= 0) THEN
ZFLXZ(IIJB:IIJE,D%NKA) = ZFLXZ(IIJB:IIJE,IKB)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
+ IF (OOCEAN) THEN
+ ZFLXZ(IIJB:IIJE,D%NKU) = ZFLXZ(IIJB:IIJE,IKE)
+ END IF
+ !
DO JK=IKTB+1,IKTE-1
!$mnh_expand_array(JIJ=IIJB:IIJE)
PWRC(IIJB:IIJE,JK)=0.5*(ZFLXZ(IIJB:IIJE,JK)+ZFLXZ(IIJB:IIJE,JK+D%NKL))
@@ -1087,10 +1067,20 @@ IF (KRR /= 0) THEN
END DO
!$mnh_expand_array(JIJ=IIJB:IIJE)
PWRC(IIJB:IIJE,IKB)=0.5*(ZFLXZ(IIJB:IIJE,IKB)+ZFLXZ(IIJB:IIJE,IKB+D%NKL))
- PWRC(IIJB:IIJE,D%NKA)=0.5*(ZFLXZ(IIJB:IIJE,D%NKA)+ZFLXZ(IIJB:IIJE,D%NKA+D%NKL))
- PWRC(IIJB:IIJE,IKE)=PWRC(IIJB:IIJE,IKE-D%NKL)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
+ IF (OOCEAN) THEN
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ PWRC(IIJB:IIJE,IKE)=0.5*(ZFLXZ(IIJB:IIJE,IKE)+ZFLXZ(IIJB:IIJE,IKE+D%NKL))
+ PWRC(IIJB:IIJE,D%NKA)=0.
+ PWRC(IIJB:IIJE,IKE+1)=ZFLXZ(IIJB:IIJE,IKE+1)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ ELSE
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ PWRC(IIJB:IIJE,D%NKA)=0.5*(ZFLXZ(IIJB:IIJE,D%NKA)+ZFLXZ(IIJB:IIJE,D%NKA+D%NKL))
+ PWRC(IIJB:IIJE,IKE)=PWRC(IIJB:IIJE,IKE-D%NKL)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ ENDIF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
! stores the conservative mixing ratio vertical flux
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 12d7e552f3ab9a388135c433337fb1c405790c7e..c457964847f5afa93c0e0434017b316ba4043bc6 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -879,8 +879,12 @@ IF (HPROGRAM/='AROME ') THEN
!
CALL UPDATE_ROTATE_WIND(D,ZUSLOPE,ZVSLOPE,HLBCX,HLBCY)
ELSE
- ZUSLOPE=PUT(IIJB:IIJE,D%NKA)
- ZVSLOPE=PVT(IIJB:IIJE,D%NKA)
+ ZUSLOPE(IIJB:IIJE)=PUT(IIJB:IIJE,D%NKA)
+ ZVSLOPE(IIJB:IIJE)=PVT(IIJB:IIJE,D%NKA)
+END IF
+IF (OOCEAN) THEN
+ ZUSLOPE(IIJB:IIJE)=PUT(IIJB:IIJE,D%NKU-1)
+ ZVSLOPE(IIJB:IIJE)=PVT(IIJB:IIJE,D%NKU-1)
END IF
!
!
@@ -893,15 +897,21 @@ ZCDUEFF(IIJB:IIJE) =-SQRT ( (PSFU(IIJB:IIJE)**2 + PSFV(IIJB:IIJE)**2) /
#else
(CST%XMNH_TINY + ZUSLOPE(IIJB:IIJE)**2 + ZVSLOPE(IIJB:IIJE)**2 ) )
#endif
+!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
!* 4.6 compute the surface tangential fluxes
!
-ZTAU11M(IIJB:IIJE) =2./3.*( (1.+ (PZZ(IIJB:IIJE,IKB+D%NKL)-PZZ(IIJB:IIJE,IKB)) &
+IF (OOCEAN) THEN
+ ZTAU11M(IIJB:IIJE)=0.
+ELSE
+ !$mnh_expand_array(JIJ=IIJB:IIJE)
+ ZTAU11M(IIJB:IIJE) =2./3.*( (1.+ (PZZ(IIJB:IIJE,IKB+D%NKL)-PZZ(IIJB:IIJE,IKB)) &
/(PDZZ(IIJB:IIJE,IKB+D%NKL)+PDZZ(IIJB:IIJE,IKB)) &
) *PTKET(IIJB:IIJE,IKB) &
-0.5 *PTKET(IIJB:IIJE,IKB+D%NKL) &
)
-!$mnh_end_expand_array(JIJ=IIJB:IIJE)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE)
+END IF
ZTAU12M(IIJB:IIJE) =0.0
ZTAU22M(IIJB:IIJE) =ZTAU11M(IIJB:IIJE)
ZTAU33M(IIJB:IIJE) =ZTAU11M(IIJB:IIJE)
@@ -989,7 +999,7 @@ CALL TURB_VER(D, CST,CSTURB,TURBN,KRR, KRRL, KRRI, &
PTSTEP,TPFILE, &
PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
PCOSSLOPE,PSINSLOPE, &
- PRHODJ,PTHVREF, &
+ PRHODJ,PTHVREF,PSFU,PSFV, &
PSFTH,PSFRV,PSFSV,PSFTH,PSFRV,PSFSV, &
ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU33M, &
PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT, &
@@ -1176,7 +1186,8 @@ CALL TKE_EPS_SOURCES(D,CST,CSTURB,BUCONF,HPROGRAM, &
& PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, &
& PTSTEP,PIMPL,ZEXPL, &
& HTURBLEN,HTURBDIM, &
- & TPFILE,OTURB_DIAG,OLES_CALL,ODIAG_IN_RUN, &
+ & TPFILE,OTURB_DIAG,OLES_CALL,ODIAG_IN_RUN,OOCEAN, &
+ & PSFU,PSFV, &
& PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,ZRTKEMS,&
& TBUDGETS,KBUDGETS, PEDR=PEDR, PTR=PTR,PDISS=PDISS, &
& PCURRENT_TKE_DISS=PCURRENT_TKE_DISS )
diff --git a/src/mesonh/ext/phys_paramn.f90 b/src/mesonh/ext/phys_paramn.f90
index 410592b19777f6cb30ad4b96c256896ea9e75a74..4e66a3cf118f4f887ad8e988824a7c7d29c3df52 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -839,28 +839,34 @@ END IF
!* 1.6 Ocean case:
! Sfc turbulent fluxes & Radiative tendency due to SW penetrating ocean
!
-IF (LOCEAN .AND. (.NOT.LCOUPLES)) THEN
+IF (LCOUPLES) THEN
+ZSFU(:,:)= XSSUFL_C(:,:,1)
+ZSFV(:,:)= XSSVFL_C(:,:,1)
+ZSFTH(:,:)= XSSTFL_C(:,:,1)
+ZSFRV(:,:)=XSSRFL_C(:,:,1)
+ELSE
+IF (LOCEAN) 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)
+ ZSFRV = 0.
+ ZSFTH = (XSSTFL_T(JSW+1)*(1.-ZSWA)+XSSTFL_T(JSW+2)*ZSWA)
+ ZSFU = (XSSUFL_T(JSW+1)*(1.-ZSWA)+XSSUFL_T(JSW+2)*ZSWA)
+ ZSFV = (XSSVFL_T(JSW+1)*(1.-ZSWA)+XSSVFL_T(JSW+2)*ZSWA)
!
ZIZOCE(IKU) = XSSOLA_T(JSW+1)*(1.-ZSWA)+XSSOLA_T(JSW+2)*ZSWA
ZPROSOL1(IKU) = CST%XROC*ZIZOCE(IKU)
ZPROSOL2(IKU) = (1.-CST%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)
+! WRITE(ILUOUT,*)'ZSWA JSW TDTCUR XTSTEP FT FU FV SolarR(IKU)', NINFRT, ZSWA,JSW,&
+! TDTCUR%xtime, XTSTEP, ZSFTH(2,2), ZSFU(2,2),ZSFV(2,2),ZIZOCE(IKU)
+ WRITE(ILUOUT,*)'XSSTP1,XSSTP,NINFRT,ZSWA,JSW,TDTCUR%xtime,ZSFT', &
+ XSSTFL_T(JSW+1),XSSTFL_T(JSW),NINFRT,ZSWA,JSW, TDTCUR%xtime,ZSFTH(2,2)
END IF
if ( TBUCONF%LBUDGET_th ) call Budget_store_init( TBUDGETS(NBUDGET_TH), 'OCEAN', xrths(:, :, :) )
DO JKM=IKU-1,2,-1
@@ -871,10 +877,12 @@ IF (LOCEAN .AND. (.NOT.LCOUPLES)) THEN
XRTHS(:,:,JKM) = XRTHS(:,:,JKM) + XRHODJ(:,:,JKM)*ZIZOCE(JKM)
END DO
if ( TBUCONF%LBUDGET_th ) call Budget_store_end ( TBUDGETS(NBUDGET_TH), 'OCEAN', xrths(:, :, :) )
+ DEALLOCATE (XSSOLA)
DEALLOCATE( ZIZOCE)
DEALLOCATE (ZPROSOL1)
DEALLOCATE (ZPROSOL2)
-END IF
+END IF! LOCEAN NO LCOUPLES
+END IF!NO LCOUPLES
!
!
CALL SECOND_MNH2(ZTIME2)
@@ -1277,14 +1285,18 @@ IF (CSURF=='EXTE') THEN
XVT(:,:,1+JPVEXT) = XVT(:,:,1+JPVEXT) - XVTRANS
END IF
!
-ELSE
- ZSFTH = 0.
- ZSFRV = 0.
+ELSE ! case no SURFEX (CSURF logical)
ZSFSV = 0.
ZSFCO2 = 0.
- ZSFU = 0.
- ZSFV = 0.
-END IF
+ IF (.NOT.LOCEAN) THEN
+ ZSFTH = 0.
+ ZSFRV = 0.
+ ZSFSV = 0.
+ ZSFCO2 = 0.
+ ZSFU = 0.
+ ZSFV = 0.
+ END IF
+END IF !CSURF
!
CALL SECOND_MNH2(ZTIME2)
!
@@ -1528,7 +1540,7 @@ IF (LOCEAN .AND. LDEEPOC) THEN
END DO
DO JJ=IJB,IJE
DO JI=IIB,IIE
- IF ( ZDIST(JI,JJ) > 1.) XSSTFL(JI,JJ)=0.
+ IF ( ZDIST(JI,JJ) > 1.) ZSFTH(JI,JJ)=0.
END DO
END DO
END IF !END DEEP OCEAN CONV CASE
@@ -1623,7 +1635,6 @@ CALL SECOND_MNH2(ZTIME4)
XUT(:,:,:) = XUT(:,:,:) - XUTRANS
XVT(:,:,:) = XVT(:,:,:) - XVTRANS
END IF
-
IF (CMF_CLOUD == 'STAT') THEN
XSIGS =SQRT( XSIGS**2 + ZSIGMF**2 )
ENDIF
@@ -1645,14 +1656,6 @@ 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