Changeset 4644 for branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

Timestamp:

2014-05-15T15:56:53+02:00 (10 years ago)

Author:

acc

Message:

Branch 2014/dev_r4642_WavesWG #1323. Import of surface wave components from the 2013/dev_ECMWF_waves branch + a few compatability changes and some mislaid documentation

File:

• branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (9 diffs)

branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -26 +26 @@
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!                 !  2012-12  (oyvind.breivik@ecwmf.int) Changed to wave breaking source term
    !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
@@ -42 +43 @@
    USE domvvl         ! domain: variable volume layer
    USE sbc_oce        ! surface boundary condition: ocean
+   USE sbcwave        ! wave parameters
+   USE sbcwave_ecmwf  ! ECMWF wave parameters
    USE zdf_oce        ! vertical physics: ocean variables
    USE zdfmxl         ! vertical physics: mixed layer
@@ -78 +81 @@
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+   REAL(wp) ::   rn_alpavg = 100.0_wp      ! Average alpha (normalized, water-side energy flux, only used in
+                                           ! case we use the source term formulation with no wave model information
+                                           ! NB! The option of running with source term and no wave model has
+                                           ! not been properly tested yet but is included for consistency (Oyvind Breivik, 2013-07-16)
+   REAL(wp) ::   rn_deptmaxtkew = 50.0_wp   ! maximum depth [m] to be affected by TKE from breaking waves 
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
@@ -85 +93 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   e0             !: top level turbulent kinetic energy   [m2/s2] from waves
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dfac           !: wave depth scaling factor   
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
@@ -109 +119 @@
       !!                ***  FUNCTION zdf_tke_alloc  ***
       !!----------------------------------------------------------------------
+      INTEGER :: iexalloc
       ALLOCATE(                                                                    &
 #if defined key_c1d
@@ -118 +129 @@
          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
          !
+      IF ( ln_wavetke ) THEN
+         ALLOCATE( e0(jpi,jpj) , dfac(jpi,jpj), STAT= iexalloc )
+         ! Setting arrays to zero initially.
+         e0(:,:) = 0.0_wp
+         dfac(:,:) = 0.0_wp
+         !
+         zdf_tke_alloc = zdf_tke_alloc + iexalloc
+      ENDIF
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
       IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
@@ -219 +238 @@
       REAL(wp) ::   zus   , zwlc  , zind            !    -         -
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
+      REAL(wp) ::   lambda !    -         -
+      REAL(wp) ::   zphio_flux              !    -         -
+      REAL(wp) ::   sbr, z0, fb                     !    -         -
+
 !!bfr      REAL(wp) ::   zebot                           !    -         -
       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
@@ -239 +262 @@
       !                     !  Surface boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      
+      ! Wave parameters read?
+      IF ( ln_wavetke ) THEN
+         ! Breaking-wave TKE injection as flux
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Roughness length proportional to wave height, capped at Hs = 0.1 m
+               z0 = 0.5_wp*MAX(swh_wavepar(ji,jj), 0.1_wp)
+               lambda = 2.38_wp/z0
+               ! The energy flux from the wave model in physical units
+               zphio_flux = MAX(phioc_wavepar(ji,jj), 0.0_wp)
+               ! TKE at surface (Mellor and Blumberg, 2004)
+               e0(ji,jj) = 0.5_wp*( 15.8_wp*zphio_flux/rau0 )**0.67_wp
+               ! Reducing the surface TKE by scaling with T-depth of first vertical level
+               dfac(ji,jj) = MAX( (2.0_wp*lambda*fsdept(ji,jj,1) )/3.0_wp, 0.00001_wp )
+               ! TKE surface boundary condition weighted by the thickness of the first level
+               en(ji,jj,1) = MAX( rn_emin0, e0(ji,jj) * ( 1.0_wp - EXP(-dfac(ji,jj)) ) / dfac(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+         ! If wave parameters are not read, use rn_ebb, corresponding to alpha = 100.0
+      ELSE
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Estimate average energy flux from std value found in rn_ebb (corresponds to alpha=100.0)
+               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+      ENDIF
+ 
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -443 +488 @@
       ENDIF
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      IF ( ln_wavetke ) THEN
+         CALL lbc_lnk( e0,     'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+         CALL lbc_lnk( dfac,   'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+      ENDIF
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer