Changeset 14007 for NEMO/trunk/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-12-02T15:30:51+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r13970 r14007 29 29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 30 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 70 74 INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 71 75 REAL(wp) :: rn_mxlice ! ice thickness value when scaling under sea-ice … … 81 85 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 82 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 83 89 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 84 90 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not … … 209 215 REAL(wp) :: zus , zwlc , zind ! - - 210 216 REAL(wp) :: zzd_up, zzd_lw ! - - 217 REAL(wp) :: ztaui, ztauj, z1_norm 211 218 INTEGER , DIMENSION(jpi,jpj) :: imlc 212 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3 219 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2 213 220 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 214 221 !!-------------------------------------------------------------------- … … 219 226 zfact2 = 1.5_wp * rn_Dt * rn_ediss 220 227 zfact3 = 0.5_wp * rn_ediss 228 ! 229 zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 221 230 ! 222 231 ! ice fraction considered for attenuation of langmuir & wave breaking … … 232 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 233 242 ! 234 DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 235 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 236 !! one way around would be to increase zbbirau 237 !! en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 238 !! & fr_i(ji,jj) * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 243 DO_2D( 0, 0, 0, 0 ) 239 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 240 248 END_2D 241 249 ! … … 274 282 ! 275 283 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 276 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) 277 285 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 278 286 ! 279 ! !* 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 280 325 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 281 326 DO jk = 2, jpk … … 283 328 & MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 284 329 END DO 285 ! !* finite Langmuir Circulation depth286 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )330 ! 331 ! !- compare LHS to RHS of Eq.47 287 332 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 288 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! Last w-level at which zpelc>=0.5*us*us 289 zus = zcof * taum(ji,jj) ! with us=0.016*wind(starting from jpk-1) 290 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 333 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 334 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 291 335 END_3D 292 336 ! ! finite LC depth … … 294 338 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 295 339 END_2D 340 ! 296 341 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 297 342 DO_2D( 0, 0, 0, 0 ) 298 zus = zcof * SQRT( taum(ji,jj) )! Stokes drift343 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 299 344 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 300 345 END_2D … … 351 396 & ) * wmask(ji,jj,jk) 352 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 ! 353 428 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 354 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1429 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 355 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) 356 431 END_3D 357 DO_2D( 0, 0, 0, 0 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 358 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 359 END_2D 360 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 ) 361 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) 362 438 END_3D … … 460 536 zmxlm(:,:,:) = rmxl_min 461 537 zmxld(:,:,:) = rmxl_min 538 ! 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 462 545 ! 463 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)464 ! 465 zraug = vkarmn * 2.e5_wp / ( rho0 * grav )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 ) 466 549 #if ! defined key_si3 && ! defined key_cice 467 DO_2D( 0, 0, 0, 0 ) ! No sea-ice468 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)469 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 470 553 #else 471 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 472 ! 473 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 ! 474 589 DO_2D( 0, 0, 0, 0 ) 475 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)590 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 476 591 END_2D 477 592 ! 478 CASE( 1 ) ! scaling with constant sea-ice thickness 479 DO_2D( 0, 0, 0, 0 ) 480 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 481 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 482 END_2D 483 ! 484 CASE( 2 ) ! scaling with mean sea-ice thickness 485 DO_2D( 0, 0, 0, 0 ) 486 #if defined key_si3 487 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 488 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 489 #elif defined key_cice 490 zmaxice = MAXVAL( h_i(ji,jj,:) ) 491 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 492 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 493 #endif 494 END_2D 495 ! 496 CASE( 3 ) ! scaling with max sea-ice thickness 497 DO_2D( 0, 0, 0, 0 ) 498 zmaxice = MAXVAL( h_i(ji,jj,:) ) 499 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 500 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 501 END_2D 502 ! 503 END SELECT 504 #endif 505 ! 506 DO_2D( 0, 0, 0, 0 ) 507 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 508 END_2D 509 ! 510 ELSE 511 zmxlm(:,:,1) = rn_mxl0 512 ENDIF 513 593 ELSE 594 zmxlm(:,:,1) = rn_mxl0 595 ENDIF 596 ENDIF 514 597 ! 515 598 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) … … 624 707 & rn_mxl0 , nn_mxlice, rn_mxlice, & 625 708 & nn_pdl , ln_lc , rn_lc , & 626 & nn_etau , nn_htau , rn_efr , nn_eice 709 & nn_etau , nn_htau , rn_efr , nn_eice , & 710 & nn_bc_surf, nn_bc_bot, ln_mxhsw 627 711 !!---------------------------------------------------------------------- 628 712 ! … … 666 750 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 667 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 668 758 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 669 759 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau
Note: See TracChangeset
for help on using the changeset viewer.
