Ignore:
Timestamp:
2020-05-29T12:58:31+02:00 (5 months ago)
Author:
emanuelaclementi
Message:

Included wave-current processes following Couvelard et al 2019 and Gurvan code -ticket #2155

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90

    r12702 r12991  
    5252   USE prtctl         ! Print control 
    5353   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     54   USE sbcwave        ! Surface boundary waves 
    5455 
    5556   IMPLICIT NONE 
     
    6263   !                      !!** Namelist  namzdf_tke  ** 
    6364   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not 
     65   LOGICAL  ::   ln_mxhsw  ! mixing length scale surface value as a fonction of wave height 
    6466   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    6567   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m] 
     
    7476   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3) 
    7577   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
     78   INTEGER  ::   nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 
     79   INTEGER  ::   nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 
    7680   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    7781   REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
     
    201205      REAL(wp) ::   zus   , zwlc  , zind       !   -         - 
    202206      REAL(wp) ::   zzd_up, zzd_lw             !   -         - 
     207      REAL(wp) ::   ztaui, ztauj, z1_norm 
    203208      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    204       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     209      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i, zWlc2 
    205210      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    206211      !!-------------------------------------------------------------------- 
     
    210215      zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    211216      zfact3 = 0.5_wp       * rn_ediss 
     217      ! 
     218      zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 
    212219      ! 
    213220      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    253260      ! 
    254261      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    255       IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002) 
     262      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke (Axell JGR 2002) 
    256263         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    257264         ! 
    258          !                        !* total energy produce by LC : cumulative sum over jk 
     265         !                       !* Langmuir velocity scale 
     266         ! 
     267         IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available 
     268            !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2 
     269            !                                ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s 
     270            !                                ! more precisely, it is the dot product that must be used : 
     271            !                                !     1/2  (W_lc)^2 = MAX( u* u_s + v* v_s , 0 )   only the positive part 
     272!!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 
     273!!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect ! 
     274            DO_2D_00_00 
     275!!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  ) 
     276                  zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) 
     277            END_2D 
     278! 
     279!  Projection of Stokes drift in the wind stress direction 
     280! 
     281            DO_2D_00_00 
     282                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 
     283                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) 
     284                  z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1) 
     285                  zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 
     286            END_2D 
     287         CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. ) 
     288! 
     289         ELSE                          ! Surface Stokes drift deduced from surface stress 
     290            !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44) 
     291            !                                ! using |tau| = rho_air Cd |U_10m|^2 , it comes: 
     292            !                                ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2   and thus: 
     293            !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 
     294            zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag 
     295            DO_2D_11_11 
     296               zWlc2(ji,jj) = zcof * taum(ji,jj) 
     297            END_2D 
     298            ! 
     299         ENDIF 
     300         ! 
     301         !                       !* Depth of the LC circulation  (Axell 2002, Eq.47) 
     302         !                             !- LHS of Eq.47 
    259303         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 
    260304         DO jk = 2, jpk 
    261305            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 
    262306         END DO 
    263          !                        !* finite Langmuir Circulation depth 
    264          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
     307         ! 
     308         !                             !- compare LHS to RHS of Eq.47 
    265309         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    266310         DO_3DS_11_11( jpkm1, 2, -1 ) 
    267             zus  = zcof * taum(ji,jj) 
    268             IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
     311            IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk 
    269312         END_3D 
    270313         !                               ! finite LC depth 
     
    272315            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 
    273316         END_2D 
     317         ! 
    274318         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    275319         DO_2D_00_00 
    276             zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     320            zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift 
    277321            zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    278322            IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     
    327371            &                                ) * wmask(ji,jj,jk) 
    328372      END_3D 
     373      ! 
     374      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     375      !                            ! choose to keep coherence with previous estimation of 
     376      !                            !  Surface boundary condition on tke 
     377      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     378      ! [EC] Should we keep this?? 
     379      IF ( ln_isfcav ) THEN 
     380         DO_2D_11_11                 ! en(mikt(ji,jj))   = rn_emin 
     381            en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 
     382         END_2D 
     383      END IF 
     384 
     385      IF ( cpl_phioc .and. ln_phioc )  THEN 
     386         SELECT CASE (nn_bc_surf) !! Dirichlet Boundary Condition using surface TKE flux from waves 
     387 
     388         CASE ( 0 ) 
     389 
     390            DO_2D_00_00            ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0) 
     391               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp 
     392               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1) 
     393               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of 
     394               zd_lw(ji,jj,1) = 1._wp              ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) 
     395               zd_up(ji,jj,1) = 0._wp 
     396            END_2D 
     397 
     398         CASE ( 1 ) 
     399            DO_2D_00_00 
     400               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp 
     401               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) 
     402               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) ) 
     403               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) 
     404               zd_lw(ji,jj,2) = 0._wp 
     405               zd_up(ji,jj,1) = 0._wp 
     406            END_2D 
     407! 
     408         END SELECT 
     409 
     410      ELSE  ! TKE Dirichlet boundary condition (without wave coupling) 
     411 
     412         DO_2D_00_00            ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0) 
     413            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     414            zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of 
     415            zd_lw(ji,jj,1) = 1._wp             ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) 
     416            zd_up(ji,jj,1) = 0._wp 
     417         END_2D 
     418 
     419      ENDIF 
     420      ! 
    329421      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    330       DO_3D_00_00( 3, jpkm1 ) 
     422!      DO_3D_00_00( 3, jpkm1 ) 
     423      DO_3D_00_00( 2, jpkm1 )         ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    331424         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    332425      END_3D 
    333       DO_2D_00_00 
    334          zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    335       END_2D 
    336       DO_3D_00_00( 3, jpkm1 ) 
     426!XC : commented to allow for neumann boundary condition 
     427!      DO_2D_00_00 
     428!         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     429!      END_2D 
     430!      DO_3D_00_00( 3, jpkm1 ) 
     431      DO_3D_00_00( 2, jpkm1 )         ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    337432         zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    338433      END_3D 
     
    340435         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    341436      END_2D 
    342       DO_3DS_00_00( jpk-2, 2, -1 ) 
     437      DO_3DS_00_00( jpk-2, 2, -1 )    ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    343438         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    344439      END_3D 
    345       DO_3D_00_00( 2, jpkm1 ) 
     440      DO_3D_00_00( 2, jpkm1 )        ! set the minimum value of tke 
    346441         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    347442      END_3D 
     
    436531      zmxld(:,:,:)  = rmxl_min 
    437532      ! 
    438       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
    439          zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 
    440          DO_2D_00_00 
    441             zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    442          END_2D 
    443       ELSE  
    444          zmxlm(:,:,1) = rn_mxl0 
     533      IF(ln_sdw .AND. ln_mxhsw) THEN 
     534         zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height) 
     535         ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6 
     536         zcoef       = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff 
     537         zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height) 
     538      ELSE 
     539         IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
     540            zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 
     541            DO_2D_00_00 
     542               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
     543            END_2D 
     544         ELSE  
     545            zmxlm(:,:,1) = rn_mxl0 
     546         ENDIF 
    445547      ENDIF 
    446548      ! 
     
    550652         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,          & 
    551653         &                 rn_mxl0 , nn_pdl   , ln_drg , ln_lc    , rn_lc,   & 
    552          &                 nn_etau , nn_htau  , rn_efr , rn_eice   
     654         &                 nn_etau , nn_htau  , rn_efr , rn_eice  ,          & 
     655         &                 nn_bc_surf, nn_bc_bot, ln_mxhsw 
    553656      !!---------------------------------------------------------------------- 
    554657      ! 
Note: See TracChangeset for help on using the changeset viewer.