Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/ZDF/zdftke.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/ZDF/zdftke.F90
r13295 r14037 28 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 29 29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 30 !! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition (ln_drg) 30 !! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition 31 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add wave coupling 32 ! ! following Couvelard et al., 2019 31 33 !!---------------------------------------------------------------------- 32 34 … … 58 60 USE prtctl ! Print control 59 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 62 USE sbcwave ! Surface boundary waves 60 63 61 64 IMPLICIT NONE … … 68 71 ! !!** Namelist namzdf_tke ** 69 72 LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not 73 LOGICAL :: ln_mxhsw ! mixing length scale surface value as a fonction of wave height 74 INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 75 REAL(wp) :: rn_mxlice ! ice thickness value when scaling under sea-ice 70 76 INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) 71 77 REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] 72 INTEGER :: nn_mxlice ! type of scaling under sea-ice73 REAL(wp) :: rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1)74 78 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 75 79 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) … … 79 83 REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2] 80 84 REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) 81 LOGICAL :: ln_drg ! top/bottom friction forcing flag82 85 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 83 86 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 87 INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 88 INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 84 89 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 85 REAL(wp) :: rn_eice ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/486 90 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 87 91 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 92 INTEGER :: nn_eice ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3) 88 93 89 94 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 200 205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 201 206 ! 202 INTEGER :: ji, jj, jk ! dummy loop arguments207 INTEGER :: ji, jj, jk ! dummy loop arguments 203 208 REAL(wp) :: zetop, zebot, zmsku, zmskv ! local scalars 204 209 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 205 210 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 206 REAL(wp) :: zbbrau, zri ! local scalars 207 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 208 REAL(wp) :: ztx2 , zty2 , zcof ! - - 209 REAL(wp) :: ztau , zdif ! - - 210 REAL(wp) :: zus , zwlc , zind ! - - 211 REAL(wp) :: zzd_up, zzd_lw ! - - 211 REAL(wp) :: zbbrau, zbbirau, zri ! local scalars 212 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 213 REAL(wp) :: ztx2 , zty2 , zcof ! - - 214 REAL(wp) :: ztau , zdif ! - - 215 REAL(wp) :: zus , zwlc , zind ! - - 216 REAL(wp) :: zzd_up, zzd_lw ! - - 217 REAL(wp) :: ztaui, ztauj, z1_norm 212 218 INTEGER , DIMENSION(jpi,jpj) :: imlc 213 REAL(wp), DIMENSION(jpi,jpj) :: z hlc, zfr_i219 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2 214 220 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 215 221 !!-------------------------------------------------------------------- 216 222 ! 217 zbbrau = rn_ebb / rho0 ! Local constant initialisation 218 zfact1 = -.5_wp * rn_Dt 219 zfact2 = 1.5_wp * rn_Dt * rn_ediss 220 zfact3 = 0.5_wp * rn_ediss 223 zbbrau = rn_ebb / rho0 ! Local constant initialisation 224 zbbirau = 3.75_wp / rho0 225 zfact1 = -.5_wp * rn_Dt 226 zfact2 = 1.5_wp * rn_Dt * rn_ediss 227 zfact3 = 0.5_wp * rn_ediss 228 ! 229 zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 230 ! 231 ! ice fraction considered for attenuation of langmuir & wave breaking 232 SELECT CASE ( nn_eice ) 233 CASE( 0 ) ; zice_fra(:,:) = 0._wp 234 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(:,:) * 10._wp ) 235 CASE( 2 ) ; zice_fra(:,:) = fr_i(:,:) 236 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 237 END SELECT 221 238 ! 222 239 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 223 240 ! ! Surface/top/bottom boundary condition on tke 224 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 225 ! 242 ! 226 243 DO_2D( 0, 0, 0, 0 ) 227 244 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 245 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) 246 zd_lw(ji,jj,1) = 1._wp 247 zd_up(ji,jj,1) = 0._wp 228 248 END_2D 229 249 ! … … 236 256 ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 237 257 ! 238 IF( ln_drg ) THEN!== friction used as top/bottom boundary condition on TKE239 ! 240 DO_2D( 0, 0, 0, 0 ) 258 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 259 ! 260 DO_2D( 0, 0, 0, 0 ) ! bottom friction 241 261 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 242 262 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) … … 246 266 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 247 267 END_2D 248 IF( ln_isfcav ) THEN ! top friction249 DO_2D( 0, 0, 0, 0 ) 268 IF( ln_isfcav ) THEN 269 DO_2D( 0, 0, 0, 0 ) ! top friction 250 270 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 251 271 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) … … 262 282 ! 263 283 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 264 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke !(Axell JGR 2002)284 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002) 265 285 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 266 286 ! 267 ! !* total energy produce by LC : cumulative sum over jk 287 ! !* Langmuir velocity scale 288 ! 289 IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available 290 ! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2 291 ! ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s 292 ! ! more precisely, it is the dot product that must be used : 293 ! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part 294 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 295 !!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect ! 296 DO_2D( 0, 0, 0, 0 ) 297 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) 298 zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) 299 END_2D 300 ! 301 ! Projection of Stokes drift in the wind stress direction 302 ! 303 DO_2D( 0, 0, 0, 0 ) 304 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 305 ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) 306 z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1) 307 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 308 END_2D 309 CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. ) 310 ! 311 ELSE ! Surface Stokes drift deduced from surface stress 312 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) 313 ! ! using |tau| = rho_air Cd |U_10m|^2 , it comes: 314 ! ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2 and thus: 315 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 316 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 317 DO_2D( 1, 1, 1, 1 ) 318 zWlc2(ji,jj) = zcof * taum(ji,jj) 319 END_2D 320 ! 321 ENDIF 322 ! 323 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 324 ! !- LHS of Eq.47 268 325 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 269 326 DO jk = 2, jpk … … 271 328 & MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 272 329 END DO 273 ! !* finite Langmuir Circulation depth274 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )330 ! 331 ! !- compare LHS to RHS of Eq.47 275 332 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 276 333 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 277 zus = zcof * taum(ji,jj) 278 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 334 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 279 335 END_3D 280 336 ! ! finite LC depth … … 282 338 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 283 339 END_2D 340 ! 284 341 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 285 342 DO_2D( 0, 0, 0, 0 ) 286 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 287 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 288 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 343 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 344 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 289 345 END_2D 290 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 291 IF ( zfr_i(ji,jj) /= 0. ) THEN 292 ! vertical velocity due to LC 346 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 347 IF ( zus3(ji,jj) /= 0._wp ) THEN 293 348 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 294 349 ! ! vertical velocity due to LC 295 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i350 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) 296 351 ! ! TKE Langmuir circulation source term 297 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * z fr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)352 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 298 353 ENDIF 299 354 ENDIF … … 309 364 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 310 365 ! 311 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri )366 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 312 367 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 313 368 ! ! local Richardson number … … 322 377 ENDIF 323 378 ! 324 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 379 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en 325 380 zcof = zfact1 * tmask(ji,jj,jk) 326 381 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 327 382 ! ! eddy coefficient (ensure numerical stability) 328 383 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 329 & / ( e3t(ji,jj,jk ,Kmm) & 330 & * e3w(ji,jj,jk ,Kmm) ) 384 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) ) 331 385 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 332 & / ( e3t(ji,jj,jk-1,Kmm) & 333 & * e3w(ji,jj,jk ,Kmm) ) 386 & / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) ) 334 387 ! 335 388 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 343 396 & ) * wmask(ji,jj,jk) 344 397 END_3D 398 ! 399 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 400 ! ! Surface boundary condition on tke if 401 ! ! coupling with waves 402 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 403 ! 404 IF ( cpl_phioc .and. ln_phioc ) THEN 405 SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 406 407 CASE ( 0 ) ! Dirichlet BC 408 DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 409 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 410 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) 411 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of 412 END_2D 413 414 CASE ( 1 ) ! Neumann BC 415 DO_2D( 0, 0, 0, 0 ) 416 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 417 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) 418 en(ji,jj,1) = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) ) 419 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) 420 zdiag(ji,jj,1) = 1._wp 421 zd_lw(ji,jj,2) = 0._wp 422 END_2D 423 424 END SELECT 425 426 ENDIF 427 ! 345 428 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 346 DO_3D( 0, 0, 0, 0, 3, jpkm1 )429 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 347 430 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 348 431 END_3D 349 DO_2D( 0, 0, 0, 0 ) 350 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 351 END_2D 352 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 432 !XC : commented to allow for neumann boundary condition 433 ! DO_2D( 0, 0, 0, 0 ) 434 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 435 ! END_2D 436 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 353 437 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) 354 438 END_3D 355 DO_2D( 0, 0, 0, 0 ) 439 DO_2D( 0, 0, 0, 0 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 356 440 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 357 441 END_2D … … 359 443 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 360 444 END_3D 361 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 445 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! set the minimum value of tke 362 446 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 363 447 END_3D … … 368 452 !!gm BUG : in the exp remove the depth of ssh !!! 369 453 !!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 370 371 454 ! 455 ! penetration is partly switched off below sea-ice if nn_eice/=0 456 ! 372 457 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 373 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 374 459 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 375 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) )* wmask(ji,jj,jk) * tmask(ji,jj,1)460 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 376 461 END_3D 377 462 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) … … 379 464 jk = nmln(ji,jj) 380 465 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 381 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) )* wmask(ji,jj,jk) * tmask(ji,jj,1)466 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 382 467 END_2D 383 468 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) … … 389 474 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 390 475 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 391 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)476 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 392 477 END_3D 393 478 ENDIF … … 452 537 zmxld(:,:,:) = rmxl_min 453 538 ! 454 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 455 ! 456 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 539 IF(ln_sdw .AND. ln_mxhsw) THEN 540 zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 541 ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6 542 zcoef = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff 543 zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 544 ELSE 545 ! 546 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 547 ! 548 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 457 549 #if ! defined key_si3 && ! defined key_cice 458 DO_2D( 0, 0, 0, 0 )459 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)460 END_2D550 DO_2D( 0, 0, 0, 0 ) ! No sea-ice 551 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 552 END_2D 461 553 #else 462 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 463 ! 464 CASE( 0 ) ! No scaling under sea-ice 554 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 555 ! 556 CASE( 0 ) ! No scaling under sea-ice 557 DO_2D( 0, 0, 0, 0 ) 558 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 559 END_2D 560 ! 561 CASE( 1 ) ! scaling with constant sea-ice thickness 562 DO_2D( 0, 0, 0, 0 ) 563 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 564 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 565 END_2D 566 ! 567 CASE( 2 ) ! scaling with mean sea-ice thickness 568 DO_2D( 0, 0, 0, 0 ) 569 #if defined key_si3 570 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 571 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 572 #elif defined key_cice 573 zmaxice = MAXVAL( h_i(ji,jj,:) ) 574 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 575 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 576 #endif 577 END_2D 578 ! 579 CASE( 3 ) ! scaling with max sea-ice thickness 580 DO_2D( 0, 0, 0, 0 ) 581 zmaxice = MAXVAL( h_i(ji,jj,:) ) 582 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 583 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 584 END_2D 585 ! 586 END SELECT 587 #endif 588 ! 465 589 DO_2D( 0, 0, 0, 0 ) 466 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)590 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 467 591 END_2D 468 592 ! 469 CASE( 1 ) ! scaling with constant sea-ice thickness 470 DO_2D( 0, 0, 0, 0 ) 471 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 472 END_2D 473 ! 474 CASE( 2 ) ! scaling with mean sea-ice thickness 475 DO_2D( 0, 0, 0, 0 ) 476 #if defined key_si3 477 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 478 #elif defined key_cice 479 zmaxice = MAXVAL( h_i(ji,jj,:) ) 480 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 481 #endif 482 END_2D 483 ! 484 CASE( 3 ) ! scaling with max sea-ice thickness 485 DO_2D( 0, 0, 0, 0 ) 486 zmaxice = MAXVAL( h_i(ji,jj,:) ) 487 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 488 END_2D 489 ! 490 END SELECT 491 #endif 492 ! 493 DO_2D( 0, 0, 0, 0 ) 494 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 495 END_2D 496 ! 497 ELSE 498 zmxlm(:,:,1) = rn_mxl0 499 ENDIF 500 593 ELSE 594 zmxlm(:,:,1) = rn_mxl0 595 ENDIF 596 ENDIF 501 597 ! 502 598 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) … … 533 629 ! 534 630 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 535 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 631 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! from the surface to the bottom : 536 632 zmxlm(ji,jj,jk) = & 537 633 & MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 538 634 END_3D 539 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 635 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface : 540 636 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 541 637 zmxlm(ji,jj,jk) = zemxl … … 544 640 ! 545 641 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 546 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 642 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! from the surface to the bottom : lup 547 643 zmxld(ji,jj,jk) = & 548 644 & MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 549 645 END_3D 550 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) 646 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown 551 647 zmxlm(ji,jj,jk) = & 552 648 & MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) … … 564 660 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 565 661 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 566 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 662 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points 567 663 zsqen = SQRT( en(ji,jj,jk) ) 568 664 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 573 669 ! 574 670 ! 575 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt671 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 576 672 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 577 673 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) … … 610 706 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 611 707 & rn_mxl0 , nn_mxlice, rn_mxlice, & 612 & nn_pdl , ln_drg , ln_lc , rn_lc, & 613 & nn_etau , nn_htau , rn_efr , rn_eice 708 & nn_pdl , ln_lc , rn_lc , & 709 & nn_etau , nn_htau , rn_efr , nn_eice , & 710 & nn_bc_surf, nn_bc_bot, ln_mxhsw 614 711 !!---------------------------------------------------------------------- 615 712 ! … … 637 734 WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl 638 735 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 736 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 639 737 IF( ln_mxl0 ) THEN 640 738 WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice 641 739 IF( nn_mxlice == 1 ) & 642 740 WRITE(numout,*) ' ice thickness when scaling under sea-ice rn_mxlice = ', rn_mxlice 643 ENDIF 644 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 645 WRITE(numout,*) ' top/bottom friction forcing flag ln_drg = ', ln_drg 741 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 742 CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice' 743 CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness' 744 CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness' 745 CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness' 746 CASE DEFAULT 747 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 748 END SELECT 749 ENDIF 646 750 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 647 751 WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc 752 IF ( cpl_phioc .and. ln_phioc ) THEN 753 SELECT CASE( nn_bc_surf) ! Type of scaling under sea-ice 754 CASE( 0 ) ; WRITE(numout,*) ' nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves' 755 CASE( 1 ) ; WRITE(numout,*) ' nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves' 756 END SELECT 757 ENDIF 648 758 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 649 759 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau 650 760 WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr 651 WRITE(numout,*) ' below sea-ice: =0 ON rn_eice = ', rn_eice 652 WRITE(numout,*) ' =4 OFF when ice fraction > 1/4 ' 653 IF( ln_drg ) THEN 654 WRITE(numout,*) 655 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:' 656 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top)= ', r_z0_top 657 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot)= ', r_z0_bot 658 ENDIF 761 WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice 762 SELECT CASE( nn_eice ) 763 CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on langmuir & surface wave breaking' 764 CASE( 1 ) ; WRITE(numout,*) ' ==>>> weigthed by 1-TANH( fr_i(:,:) * 10 )' 765 CASE( 2 ) ; WRITE(numout,*) ' ==>>> weighted by 1-fr_i(:,:)' 766 CASE( 3 ) ; WRITE(numout,*) ' ==>>> weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 767 CASE DEFAULT 768 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 769 END SELECT 659 770 WRITE(numout,*) 660 771 WRITE(numout,*) ' ==>>> critical Richardson nb with your parameters ri_cri = ', ri_cri … … 700 811 CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl) 701 812 ! 702 IF( lwxios ) THEN703 CALL iom_set_rstw_var_active('en')704 CALL iom_set_rstw_var_active('avt_k')705 CALL iom_set_rstw_var_active('avm_k')706 CALL iom_set_rstw_var_active('dissl')707 ENDIF708 813 END SUBROUTINE zdf_tke_init 709 814 … … 737 842 ! 738 843 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 739 CALL iom_get( numror, jpdom_auto, 'en' , en , ldxios = lrxios)740 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k , ldxios = lrxios)741 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k , ldxios = lrxios)742 CALL iom_get( numror, jpdom_auto, 'dissl', dissl , ldxios = lrxios)844 CALL iom_get( numror, jpdom_auto, 'en' , en ) 845 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k ) 846 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k ) 847 CALL iom_get( numror, jpdom_auto, 'dissl', dissl ) 743 848 ELSE ! start TKE from rest 744 849 IF(lwp) WRITE(numout,*) … … 759 864 ! ! ------------------- 760 865 IF(lwp) WRITE(numout,*) '---- tke_rst ----' 761 IF( lwxios ) CALL iom_swap( cwxios_context ) 762 CALL iom_rstput( kt, nitrst, numrow, 'en' , en , ldxios = lwxios ) 763 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios ) 764 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios ) 765 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl, ldxios = lwxios ) 766 IF( lwxios ) CALL iom_swap( cxios_context ) 866 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 867 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 868 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 869 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 767 870 ! 768 871 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
