diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index bb4ba8aeffac377ce44021b2f14a7edccb4a60ca..2b9faf5642cc6be5e6a697603e8712dba218d0e4 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -8,10 +8,10 @@ CONTAINS
SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP, &
& PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, &
& PTSTEP,PIMPL,PEXPL, &
- & HTURBLEN,HTURBDIM, &
- & TPFILE,OTURB_DIAG,ODIAG_IN_RUN, &
+ & HTURBLEN,HTURBDIM,OOCEAN, &
+ & TPFILE,OTURB_DIAG,ODIAG_IN_RUN,PSFUM,PSFVM, &
& PTP,PRTKES,PRTHLS,PCOEF_DISS,PTDIFF,PTDISS,PRTKEMS,&
- & TBUDGETS, KBUDGETS, &
+ & TBUDGETS, KBUDGETS, &
& PEDR, PTR,PDISS, PCURRENT_TKE_DISS )
! ##################################################################
!
@@ -178,8 +178,10 @@ CHARACTER(LEN=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
CHARACTER(LEN=4), INTENT(IN) :: HTURBLEN ! kind of mixing length
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some
-LOGICAL, INTENT(IN) :: ODIAG_IN_RUN ! switch to activate online diagnostics (mesonh)
- ! diagnostic fields in the syncronous FM-file
+ ! diagnostic fields in the syncronous FM-file
+LOGICAL, INTENT(IN) :: ODIAG_IN_RUN ! switch to activate online diagnostics (mesonh)
+LOGICAL, INTENT(IN) :: OOCEAN ! switch to activate LES OCEAN version
+
REAL, DIMENSION(:,:,:), INTENT(INOUT):: PDP ! Dyn. prod. of TKE
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTRH
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTP ! Ther. prod. of TKE
@@ -196,6 +198,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PTR ! Transport pro
REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PDISS ! Dissipation of TKE
REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PEDR ! EDR
REAL, DIMENSION(:,:,:), INTENT(INOUT), OPTIONAL :: PCURRENT_TKE_DISS ! if ODIAG_IN_RUN in mesonh
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFUM,PSFVM ! momentum sfc flux
!
!
!
@@ -263,8 +266,13 @@ 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)
-PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+KKL)/PDZZ(:,:,IKB))
+IF (LOCEAN) THEN
+ ! W(IKE) value stored in PDP(IKE) to the mass localization of tke(IKE)
+ PDP(:,:,IKE) = PDP(:,:,IKE) * (1. + PDZZ(:,:,IKE)/PDZZ(:,:,IKE+1))
+ELSE
+ ! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
+ PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+KKL)/PDZZ(:,:,IKB))
+END IF
!
! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
! + (Dynamical Production) + (Thermal Production) - (dissipation)
@@ -276,6 +284,12 @@ ZSOURCE(:,:,:) = ( PRTKES(:,:,:) + PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
!* 2.2 implicit vertical TKE transport
!
!
+! To add here in ZSOURCE surface flux of TKE
+!(assumed to be 0 for ATM,
+IF (LOCEAN) THEN
+ !for ocean:wave breaking simple/very rough param wE = 100 Ustar**3 where ustar is the Tau_atmi/rhocea
+ ZSOURCE (:,:,IKE)=ZSOURCE(:,:,IKE)-1.E2*((PSFUM(:,:)**2 + PSFVM(:,:)**2)**1.5) /PDZZ(:,:,IKE)
+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 49ed101bb374d335175ae8b26f446b5bdf427f44..393377797d4a4af6e802590439c6089404d96c7e 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -11,7 +11,7 @@ SUBROUTINE TURB_VER(KKA,KKU,KKL,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, &
@@ -275,6 +275,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual sc
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
! Potential Temperature
!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFUM,PSFVM ! surface fluxe
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! t - deltat
!
@@ -556,7 +557,7 @@ CALL TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
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 e4c1dff95ea98c9bd16bd7e63145bed84c94763e..9d018ff6e338b1a090d65ed8f0e0fbf17d94d871 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -13,6 +13,7 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
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, &
@@ -217,7 +218,6 @@ USE MODD_LES
USE MODD_NSV
USE MODD_OCEANH
USE MODD_PARAMETERS
-USE MODD_REF, ONLY : LCOUPLES
USE MODD_TURB_n
!
!
@@ -267,6 +267,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFUM,PSFVM ! normal momentum sfc flux
REAL, DIMENSION(:,:), INTENT(IN) :: PCDUEFF ! Cd * || u || at time t
REAL, DIMENSION(:,:), INTENT(IN) :: PTAU11M ! <uu> in the axes linked
! to the maximum slope direction and the surface normal and the binormal
@@ -325,6 +326,7 @@ REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
! coefficients for the surface flux
! evaluation and copy of PUSLOPEM and
! PVSLOPEM in local 3D arrays
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZFLUXSFCU,ZFLUXSFCV
INTEGER :: IIU,IJU ! size of array in x,y,z directions
!
REAL :: ZTIME1, ZTIME2, ZCMFS
@@ -372,6 +374,8 @@ ENDIF
!
ZUSLOPEM(:,:,1)=PUSLOPEM(:,:)
ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
+ZFLUXSFCU(:,:,1)=PSFUM(:,:)
+ZFLUXSFCV(:,:,1)=PSFVM(:,:)
!
!----------------------------------------------------------------------------
!
@@ -395,7 +399,7 @@ ZA(:,:,:) = -PTSTEP * ZCMFS * &
ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
* PCOSSLOPE(:,:)
ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PSINSLOPE(:,:)
-
+!
! prepare the implicit scheme coefficients for the surface flux
ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) &
+ZCOEFFLXV(:,:,1) * PSINSLOPE(:,:)
@@ -403,28 +407,16 @@ 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) )
!
-!
-! ZSOURCE= FLUX /DZ
-IF (OOCEAN) 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(:,:,IKTB+1:IKTE-1) = 0
+! ZSOURCE= sfc FLUX /DZ
+! Sfx flux assumed to be in SI & at vorticity point
+IF (OOCEAN) THEN ! Ocean model
+ ZSOURCE(:,:,IKE:IKE) =MXM( ZFLUXSFCU(:,:,1:1) / PDZZ(:,:,IKE:IKE) ) &
+ * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
+ ! Zero 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
+ELSE ! Atmosphere
+ ! Compute the explicit tangential flux at the W point
ZSOURCE(:,:,IKB) = &
PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
-PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
@@ -440,11 +432,9 @@ ELSE !ATMOS MODEL ONLY
*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.
+ ! Zero flux at upper BC for atmos case
ZSOURCE(:,:,IKE) = 0.
-ENDIF !end ocean or atmosphere cases
+ENDIF
!
! Obtention of the split U at t+ deltat
!
@@ -455,27 +445,21 @@ CALL TRIDIAG_WIND(KKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
!
PRUS(:,:,:)=PRUS(:,:,:)+MXM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
!
-!
-!* 5.2 Partial Dynamic Production
+!* 5.2 Partial TKE Dynamic Production
!
! vertical flux of the U wind component
!
ZFLXZ(:,:,:) = -ZCMFS * MXM(ZKEFF) * &
DZM(PIMPL*ZRES + PEXPL*PUM, KKA, KKU, KKL) / MXM(PDZZ)
-!
-! surface flux
-ZFLXZ(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB)) * &
+IF (OOCEAN) THEN
+ ZFLXZ(:,:,IKE+1) = ZFLUXSFCU(:,:,1)
+ELSE
+ ZFLXZ(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB)) * &
( ZSOURCE(:,:,IKB:IKB) &
+ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
) / 0.5 / ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
!
-ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
-
-IF (OOCEAN) 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)
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -497,23 +481,25 @@ END IF
!
PWU(:,:,:) = ZFLXZ(:,:,:)
!
-! 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)
!
PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
!
+! Special cases near surface
+IF (OOCEAN) THEN
+! evaluate the dynamic production at w(IKE) and store in PDP(IKE)
+! before to be extrapolated in tke_eps routine
+ PDP(:,:,IKE:IKE) = - MXF ( &
+ ZFLXZ(:,:,IKE:IKE) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-1:IKE-1)) &
+ / MXM(PDZZ(:,:,IKE-1:IKE-1)) &
+ )
+ELSE ! Atmosphere
! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
-PDP(:,:,IKB:IKB) = - MXF ( &
- ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PUM(:,:,IKB+KKL:IKB+KKL)-PUM(:,:,IKB:IKB)) &
+ PDP(:,:,IKB:IKB) = - MXF ( &
+ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PUM(:,:,IKB+KKL:IKB+KKL)-PUM(:,:,IKB:IKB)) &
/ MXM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
)
-!
-IF (OOCEAN) 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
@@ -535,10 +521,7 @@ IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
! used to compute the W source at the ground
ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
- IF (OOCEAN) THEN
- ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
- END IF
-
+ IF (OOCEAN) ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
!
IF (.NOT. LFLAT) THEN
PRWS(:,:,:)= PRWS &
@@ -550,11 +533,19 @@ IF(HTURBDIM=='3DIM') THEN
PRWS(:,:,:)= PRWS -DXF(MZM(MXM(PRHODJ) /PDXX, KKA, KKU, KKL) * ZFLXZ )
END IF
!
- ! Complete the Dynamical production with the W wind component
+ ! Complete the TKE dynamical production with the W wind contribution
!
ZA(:,:,:)=-MZF(MXF(ZFLXZ * GX_W_UW(PWM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)), KKA, KKU, KKL)
!
- !
+! Special cases near surface
+ IF (OOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE) and stored in PDP(IKE)
+ ZA(:,:,IKE:IKE) = - MXF ( &
+ ZFLXZ(:,:,IKE:IKE) * &
+ DXM( PWM(:,:,IKE:IKE) ) &
+ / (0.5*(PDXX(:,:,IKE-1:IKE-1)+PDXX(:,:,IKE:IKE))) &
+ )
+ ELSE !Atmosphere
! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
ZA(:,:,IKB:IKB) = - MXF ( &
ZFLXZ(:,:,IKB+KKL:IKB+KKL) * &
@@ -567,23 +558,9 @@ IF(HTURBDIM=='3DIM') THEN
* PDZX(:,:,IKB+KKL:IKB+KKL) &
) / (0.5*(PDXX(:,:,IKB+KKL:IKB+KKL)+PDXX(:,:,IKB:IKB))) &
)
+ END IF
!
-IF (OOCEAN) 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(:,:,:)
+ PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
! Storage in the LES configuration
!
@@ -609,7 +586,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
@@ -636,19 +613,23 @@ 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) )
!
+
+! No flux in SOURCE TERM NULL OUTSIDE BC
+ZSOURCE(:,:,IKB+1:IKE-1) = 0.
+! Surface case
IF (OOCEAN) 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.
+ ZCOEFFLXU(:,:,1) = PCDUEFF(:,:)
+ ZCOEFFLXV(:,:,1) = PCDUEFF(:,:)
+ ZCOEFS(:,:,1)=ZCOEFFLXU(:,:,1)
+! average this flux to be located at the U,W vorticity point
+ ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKE:IKE) )
+ ZSOURCE(:,:,IKE:IKE) = MYM( ZFLUXSFCV(:,:,1:1) / PDZZ(:,:,IKE:IKE) ) &
+ * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+!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(:,:) &
@@ -665,17 +646,11 @@ ELSE ! Atmos case
- 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
+!No flux at the atmosphere top
ZSOURCE(:,:,IKE) = 0.
-ENDIF ! End of Ocean or Atmospher Cases
-ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
-!
-! Obtention of the split V at t+ deltat
+ENDIF ! End of Ocean or Atmospheric Cases
+!
+! 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)
!
@@ -691,20 +666,17 @@ PRVS(:,:,:)=PRVS(:,:,:)+MYM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PVM(:,:,:))/PTSTEP
ZFLXZ(:,:,:) = -ZCMFS * MYM(ZKEFF) * &
DZM(PIMPL*ZRES + PEXPL*PVM, KKA, KKU, KKL) / MYM(PDZZ)
!
-ZFLXZ(:,:,IKB:IKB) = MYM(PDZZ(:,:,IKB:IKB)) * &
+!
+IF (OOCEAN) THEN
+ ZFLXZ(:,:,IKE+1) = ZFLUXSFCV(:,:,1)
+ELSE
+ ZFLXZ(:,:,IKB:IKB) = MYM(PDZZ(:,:,IKB:IKB)) * &
( ZSOURCE(:,:,IKB:IKB) &
+ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
) / 0.5 / ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
!
-!
-ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
-!
-IF (OOCEAN) 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
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
! stores the V wind component vertical flux
@@ -725,24 +697,25 @@ END IF
!
PWV(:,:,:) = ZFLXZ(:,:,:)
!
-! 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
!
ZA(:,:,:) = - MZF(MYF(ZFLXZ * GZ_V_VW(PVM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
!
-! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
-ZA(:,:,IKB:IKB) = &
- - MYF ( &
-ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PVM(:,:,IKB+KKL:IKB+KKL)-PVM(:,:,IKB:IKB)) &
- / MYM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
- )
-!
+! Special cases at surface
IF (OOCEAN) THEN
- ! evaluate the dynamic production at w(IKE-KKL) in PDP(IKE)
+! evaluate the dynamic production at w(IKE) and stored in PDP(IKE)
+ ! before extrapolation done in routine tke_eps_source
ZA(:,:,IKE:IKE) = - MYF ( &
- ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-KKL:IKE-KKL)) &
- / MYM(PDZZ(:,:,IKE-KKL:IKE-KKL)) &
+ ZFLXZ(:,:,IKE:IKE) * (PVM(:,:,IKE:IKE)-PVM(:,:,IKE-1:IKE-1)) &
+ / MYM(PDZZ(:,:,IKE-1:IKE-1)) &
)
+ELSE ! Atmosphere
+! evaluate the dynamic production at w(IKB+KL) and stored in PDP(IKB)
+ ZA(:,:,IKB:IKB) = - MYF ( &
+ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PVM(:,:,IKB+KKL:IKB+KKL)-PVM(:,:,IKB:IKB)) &
+ / MYM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
+ )
END IF
!
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
@@ -763,10 +736,11 @@ END IF
!
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
- ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
- IF (OOCEAN) THEN
- ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
- END IF
+ IF (OOCEAN) THEN
+ ZFLXZ(:,:,IKE+1) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-1) ! extrapolation
+ ELSE
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ END IF
!
IF (.NOT. L2D) THEN
IF (.NOT. LFLAT) THEN
@@ -780,11 +754,18 @@ IF(HTURBDIM=='3DIM') THEN
END IF
END IF
!
- ! Complete the Dynamical production with the W wind component
+ ! Complete the TKE dynamical production with the W wind component
IF (.NOT. L2D) THEN
ZA(:,:,:) = - MZF(MYF(ZFLXZ * GY_W_VW(PWM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)), KKA, KKU, KKL)
- !
- ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+ ! Special case near surface
+ IF (OOCEAN) THEN
+ ! evaluate the dynamic production at w(IKE) and stored in PDP(IKE)
+ ZA(:,:,IKE:IKE) = - MYF ( &
+ ZFLXZ(:,:,IKE:IKE) * DYM( PWM(:,:,IKE:IKE) ) &
+ / (0.5*(PDYY(:,:,IKE-1:IKE-1)+PDYY(:,:,IKE:IKE))) &
+ )
+ ELSE ! Atmosphere
+ ! evaluate the dynamic production at w(IKB+KKL) and stored in PDP(IKB)
ZA(:,:,IKB:IKB) = - MYF ( &
ZFLXZ(:,:,IKB+KKL:IKB+KKL) * &
( DYM( PWM(:,:,IKB+KKL:IKB+KKL) ) &
@@ -797,18 +778,6 @@ IF(HTURBDIM=='3DIM') THEN
) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB))) &
)
!
- IF (OOCEAN) 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(:,:,:)
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index 4aac5c81ae2b9ec635b0f59654fc4c69924a2596..ae6faa452abe03daa215f4bcceab3935dada304d 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -241,7 +241,6 @@ USE MODD_TURB_n, ONLY: LHGRAD, XCOEFHGRADTHL, XCOEFHGRADRM, XALTHGRAD, X
USE MODD_CONF
USE MODD_LES
USE MODD_OCEANH
-USE MODD_REF, ONLY: LCOUPLES
USE MODD_TURB_n
!
USE MODI_GRADIENT_U
@@ -528,18 +527,10 @@ 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 (OOCEAN) 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
+! specialcase for surface
+IF (OOCEAN) THEN ! ocean sfc (domain top)
+ ZF(:,:,IKE+1) = PSFTHM(:,:) *0.5*(1. + PRHODJ(:,:,KKU) / PRHODJ(:,:,IKE))
+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
@@ -552,17 +543,14 @@ ELSE ! No coupling O and A cases
/ PDIRCOSZW(:,:) &
* 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
END IF
-!
- IF (OOCEAN) 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.
+ !TODO merge : the following solution must be kept :
+ !ZF(:,:,IKE+1)=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,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
@@ -597,9 +585,10 @@ IF (LHGRAD) THEN
END WHERE
END IF
!
-ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
IF (OOCEAN) THEN
- ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+ ZFLXZ(:,:,IKE+1) = ZFLXZ(:,:,IKE)
+ELSE ! ATMOS
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
END IF
!
DO JK=IKTB+1,IKTE-1
@@ -610,8 +599,8 @@ PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
!
IF (OOCEAN) THEN
PWTH(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+KKL))
- PWTH(:,:,KKA)=0.
- PWTH(:,:,KKU)=ZFLXZ(:,:,KKU)
+ PWTH(:,:,KKA)=0.
+ PWTH(:,:,KKU)=PWTH(:,:,IKE)! not used
ELSE
PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
@@ -780,17 +769,12 @@ IF (KRR /= 0) THEN
& PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, KKA, KKU, KKL)
END IF
!
- ! compute interface flux
- IF (LCOUPLES) THEN ! coupling NH O-A
- IF (OOCEAN) THEN ! ocean model in coupled case
- ! evap effect on salinity to be added later !!!
- ZF(:,:,IKE) = 0.
- ELSE ! atmosph model in coupled case
- ZF(:,:,IKB) = 0.
- ! AJOUTER FLUX EVAP SUR MODELE ATMOS
- ENDIF
- !
- ELSE ! No coupling NH OA case
+ !special case at sfc
+ IF (OOCEAN) THEN
+ ! General ocean case
+ ! salinity/evap effect to be added later !!!!!
+ ZF(:,:,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)
@@ -806,19 +790,12 @@ IF (KRR /= 0) THEN
/ PDIRCOSZW(:,:) &
* 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
END IF
- !
- IF (OOCEAN) 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)
@@ -855,13 +832,23 @@ IF (KRR /= 0) THEN
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
+ IF (OOCEAN) THEN
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+ END IF
+ !
DO JK=IKTB+1,IKTE-1
PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
END DO
PWRC(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
- PWRC(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
- PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
!
+ IF (OOCEAN) THEN
+ PWRC(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+KKL))
+ PWRC(:,:,KKA)=0.
+ PWRC(:,:,IKE+1)=ZFLXZ(:,:,IKE+1)
+ ELSE
+ PWRC(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
+ PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
+ 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 63fd923053c57e7c79e2bcdd40e1ff242c83d09e..ace7ca393fe23c64b35fe4fea93090cee9c6cd75 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -788,6 +788,10 @@ ELSE
ZUSLOPE=PUT(:,:,KKA)
ZVSLOPE=PVT(:,:,KKA)
END IF
+IF (OOCEAN) THEN
+ ZUSLOPE=PUT(:,:,KKU-1)
+ ZVSLOPE=PVT(:,:,KKU-1)
+END IF
!
!
!* 4.2 compute the proportionality coefficient between wind and stress
@@ -801,11 +805,15 @@ ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / &
!
!* 4.6 compute the surface tangential fluxes
!
-ZTAU11M(:,:) =2./3.*( (1.+ (PZZ (:,:,IKB+KKL)-PZZ (:,:,IKB)) &
- /(PDZZ(:,:,IKB+KKL)+PDZZ(:,:,IKB)) &
- ) *PTKET(:,:,IKB) &
- -0.5 *PTKET(:,:,IKB+KKL) &
- )
+IF (OOCEAN) THEN
+ ZTAU11M(:,:)=0.
+ELSE
+ ZTAU11M(:,:) =2./3.*( (1.+ (PZZ (:,:,IKB+KKL)-PZZ (:,:,IKB)) &
+ /(PDZZ(:,:,IKB+KKL)+PDZZ(:,:,IKB)) &
+ ) *PTKET(:,:,IKB) &
+ -0.5 *PTKET(:,:,IKB+KKL) &
+ )
+END IF
ZTAU12M(:,:) =0.0
ZTAU22M(:,:) =ZTAU11M(:,:)
ZTAU33M(:,:) =ZTAU11M(:,:)
@@ -895,7 +903,7 @@ CALL TURB_VER(KKA,KKU,KKL,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, &
@@ -1070,8 +1078,8 @@ END IF
CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH, &
& PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, &
& PTSTEP,PIMPL,ZEXPL, &
- & HTURBLEN,HTURBDIM, &
- & TPFILE,OTURB_DIAG,ODIAG_IN_RUN, &
+ & HTURBLEN,HTURBDIM,OOCEAN, &
+ & TPFILE,OTURB_DIAG,ODIAG_IN_RUN,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 e77c51456e785de14018acab81cfbdb27fbb1052..9ad4e40a4ec0c238e5a823d914d64b8bb99de2c7 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -828,28 +828,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) = XROC*ZIZOCE(IKU)
ZPROSOL2(IKU) = (1.-XROC)*ZIZOCE(IKU)
IF(NVERB >= 5 ) THEN
- WRITE(ILUOUT,*)'ZSWA JSW TDTCUR XTSTEP FT FU FV SolarR(IKU)', NINFRT, ZSWA,JSW,&
- TDTCUR%xtime, XTSTEP, XSSTFL(2,2), XSSUFL(2,2),XSSVFL(2,2),ZIZOCE(IKU)
+! 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 ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'OCEAN', xrths(:, :, :) )
DO JKM=IKU-1,2,-1
@@ -860,10 +866,12 @@ IF (LOCEAN .AND. (.NOT.LCOUPLES)) THEN
XRTHS(:,:,JKM) = XRTHS(:,:,JKM) + XRHODJ(:,:,JKM)*ZIZOCE(JKM)
END DO
if ( 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)
@@ -1266,14 +1274,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)
!
@@ -1516,7 +1528,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
@@ -1603,7 +1615,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
@@ -1625,14 +1636,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