Changeset 12991 for NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-05-29T12:58:31+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90
r12702 r12991 52 52 USE prtctl ! Print control 53 53 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 54 USE sbcwave ! Surface boundary waves 54 55 55 56 IMPLICIT NONE … … 62 63 ! !!** Namelist namzdf_tke ** 63 64 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 64 66 INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) 65 67 REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] … … 74 76 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 75 77 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 76 80 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 77 81 REAL(wp) :: rn_eice ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4 … … 201 205 REAL(wp) :: zus , zwlc , zind ! - - 202 206 REAL(wp) :: zzd_up, zzd_lw ! - - 207 REAL(wp) :: ztaui, ztauj, z1_norm 203 208 INTEGER , DIMENSION(jpi,jpj) :: imlc 204 REAL(wp), DIMENSION(jpi,jpj) :: zhlc, zfr_i 209 REAL(wp), DIMENSION(jpi,jpj) :: zhlc, zfr_i, zWlc2 205 210 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 206 211 !!-------------------------------------------------------------------- … … 210 215 zfact2 = 1.5_wp * rn_Dt * rn_ediss 211 216 zfact3 = 0.5_wp * rn_ediss 217 ! 218 zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 212 219 ! 213 220 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 253 260 ! 254 261 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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) 256 263 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 257 264 ! 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 259 303 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 260 304 DO jk = 2, jpk 261 305 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 262 306 END DO 263 ! !* finite Langmuir Circulation depth264 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )307 ! 308 ! !- compare LHS to RHS of Eq.47 265 309 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 266 310 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 269 312 END_3D 270 313 ! ! finite LC depth … … 272 315 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 273 316 END_2D 317 ! 274 318 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 275 319 DO_2D_00_00 276 zus = zcof * SQRT( taum(ji,jj) )! Stokes drift320 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 277 321 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 278 322 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. … … 327 371 & ) * wmask(ji,jj,jk) 328 372 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 ! 329 421 ! !* 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 331 424 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 332 425 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 337 432 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) 338 433 END_3D … … 340 435 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 341 436 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 343 438 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 344 439 END_3D 345 DO_3D_00_00( 2, jpkm1 ) 440 DO_3D_00_00( 2, jpkm1 ) ! set the minimum value of tke 346 441 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 347 442 END_3D … … 436 531 zmxld(:,:,:) = rmxl_min 437 532 ! 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 445 547 ENDIF 446 548 ! … … 550 652 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 551 653 & 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 553 656 !!---------------------------------------------------------------------- 554 657 !
Note: See TracChangeset
for help on using the changeset viewer.
