- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r6225 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 28 29 !!---------------------------------------------------------------------- 29 #if defined key_zdftke || defined key_esopa30 #if defined key_zdftke 30 31 !!---------------------------------------------------------------------- 31 32 !! 'key_zdftke' TKE vertical physics … … 52 53 USE timing ! Timing 53 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 #if defined key_agrif 56 USE agrif_opa_interp 57 USE agrif_opa_update 58 #endif 54 59 55 60 IMPLICIT NONE … … 84 89 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 85 90 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2]87 91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 88 92 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k ! not enhanced Kz 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmu_k, avmv_k ! not enhanced Kz 93 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 91 94 #if defined key_c1d 92 95 ! !!** 1D cfg only ** ('key_c1d') … … 96 99 97 100 !! * Substitutions 98 # include "domzgr_substitute.h90"99 101 # include "vectopt_loop_substitute.h90" 100 102 !!---------------------------------------------------------------------- 101 !! NEMO/OPA 4.0 , NEMO Consortium (2011)103 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 102 104 !! $Id$ 103 105 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 114 116 & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & 115 117 #endif 116 & en (jpi,jpj,jpk) , htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 117 & avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk), & 118 & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc ) 118 & htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 119 & apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 119 120 ! 120 121 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) … … 172 173 !!---------------------------------------------------------------------- 173 174 ! 175 #if defined key_agrif 176 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 177 IF( .NOT.Agrif_Root() ) CALL Agrif_Tke 178 #endif 179 ! 174 180 IF( kt /= nit000 ) THEN ! restore before value to compute tke 175 181 avt (:,:,:) = avt_k (:,:,:) … … 188 194 avmv_k(:,:,:) = avmv(:,:,:) 189 195 ! 190 END SUBROUTINE zdf_tke 196 #if defined key_agrif 197 ! Update child grid f => parent grid 198 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 199 #endif 200 ! 201 END SUBROUTINE zdf_tke 191 202 192 203 … … 220 231 REAL(wp) :: zzd_up, zzd_lw ! - - 221 232 !!bfr REAL(wp) :: zebot ! - - 222 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 223 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 224 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 233 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 234 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 235 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 236 REAL(wp) :: zri ! local Richardson number 225 237 !!-------------------------------------------------------------------- 226 238 ! 227 239 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 228 240 ! 229 CALL wrk_alloc( jpi,jpj, imlc ) ! integer230 CALL wrk_alloc( jpi,jpj, zhlc )231 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw)241 CALL wrk_alloc( jpi,jpj, imlc ) ! integer 242 CALL wrk_alloc( jpi,jpj, zhlc ) 243 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 232 244 ! 233 245 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 236 248 zfact3 = 0.5_wp * rn_ediss 237 249 ! 250 ! 238 251 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 239 252 ! ! Surface boundary condition on tke 240 253 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 254 IF ( ln_isfcav ) THEN 255 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 256 DO ji = fs_2, fs_jpim1 ! vector opt. 257 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 258 END DO 259 END DO 260 END IF 241 261 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 262 DO ji = fs_2, fs_jpim1 ! vector opt. … … 256 276 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 257 277 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 258 !CDIR NOVERRCHK259 278 !! DO jj = 2, jpjm1 260 !CDIR NOVERRCHK261 279 !! DO ji = fs_2, fs_jpim1 ! vector opt. 262 280 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & … … 275 293 ! 276 294 ! !* total energy produce by LC : cumulative sum over jk 277 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)295 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1) 278 296 DO jk = 2, jpk 279 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)297 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk) 280 298 END DO 281 299 ! !* finite Langmuir Circulation depth … … 291 309 END DO 292 310 ! ! finite LC depth 293 # if defined key_vectopt_loop294 DO jj = 1, 1295 DO ji = 1, jpij ! vector opt. (forced unrolling)296 # else297 311 DO jj = 1, jpj 298 312 DO ji = 1, jpi 299 # endif 300 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 313 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj)) 301 314 END DO 302 315 END DO 303 316 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 304 !CDIR NOVERRCHK305 317 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 306 !CDIR NOVERRCHK 307 DO jj = 2, jpjm1 308 !CDIR NOVERRCHK 318 DO jj = 2, jpjm1 309 319 DO ji = fs_2, fs_jpim1 ! vector opt. 310 320 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 311 321 ! ! vertical velocity due to LC 312 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )313 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )322 zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) ) 323 zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) ) 314 324 ! ! TKE Langmuir circulation source term 315 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)325 en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 316 326 END DO 317 327 END DO … … 328 338 ! 329 339 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 330 DO jj = 1, jpj ! here avmu, avmv used as workspace 331 DO ji = 1, jpi 332 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 333 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 334 & / ( fse3uw_n(ji,jj,jk) & 335 & * fse3uw_b(ji,jj,jk) ) 336 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 337 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 338 & / ( fse3vw_n(ji,jj,jk) & 339 & * fse3vw_b(ji,jj,jk) ) 340 END DO 341 END DO 342 END DO 343 ! 340 DO jj = 1, jpjm1 341 DO ji = 1, fs_jpim1 ! vector opt. 342 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 343 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 344 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 345 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 346 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 347 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 348 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 349 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 350 END DO 351 END DO 352 END DO 353 ! 354 IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr 355 ! Note that zesh2 is also computed in the next loop. 356 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 357 DO jk = 2, jpkm1 358 DO jj = 2, jpjm1 359 DO ji = fs_2, fs_jpim1 ! vector opt. 360 ! ! shear prod. at w-point weightened by mask 361 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 362 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 363 ! ! local Richardson number 364 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 365 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 366 367 END DO 368 END DO 369 END DO 370 ! 371 ENDIF 372 ! 344 373 DO jk = 2, jpkm1 !* Matrix and right hand side in en 345 374 DO jj = 2, jpjm1 … … 347 376 zcof = zfact1 * tmask(ji,jj,jk) 348 377 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal 349 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) )378 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) ) 350 379 zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal 351 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )352 !! shear prod. at w-point weightened by mask353 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &354 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )355 380 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 381 ! ! shear prod. at w-point weightened by mask 382 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 383 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 384 ! 356 385 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 357 386 zd_lw(ji,jj,jk) = zzd_lw … … 360 389 ! ! right hand side in en 361 390 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 362 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) 391 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 392 & * wmask(ji,jj,jk) 363 393 END DO 364 394 END DO … … 373 403 END DO 374 404 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 375 DO ji = fs_2, fs_jpim1 405 DO ji = fs_2, fs_jpim1 ! vector opt. 376 406 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 377 407 END DO … … 385 415 END DO 386 416 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 387 DO ji = fs_2, fs_jpim1 417 DO ji = fs_2, fs_jpim1 ! vector opt. 388 418 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 389 419 END DO … … 399 429 DO jj = 2, jpjm1 400 430 DO ji = fs_2, fs_jpim1 ! vector opt. 401 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)431 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 402 432 END DO 403 433 END DO … … 407 437 ! ! TKE due to surface and internal wave breaking 408 438 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 439 !!gm BUG : in the exp remove the depth of ssh !!! 440 441 409 442 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 410 443 DO jk = 2, jpkm1 411 444 DO jj = 2, jpjm1 412 445 DO ji = fs_2, fs_jpim1 ! vector opt. 413 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - fsdepw(ji,jj,jk) / htau(ji,jj) ) &414 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)446 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) & 447 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 415 448 END DO 416 449 END DO … … 420 453 DO ji = fs_2, fs_jpim1 ! vector opt. 421 454 jk = nmln(ji,jj) 422 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - fsdepw(ji,jj,jk) / htau(ji,jj) ) &423 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)455 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) & 456 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 424 457 END DO 425 458 END DO 426 459 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 427 !CDIR NOVERRCHK428 460 DO jk = 2, jpkm1 429 !CDIR NOVERRCHK 430 DO jj = 2, jpjm1 431 !CDIR NOVERRCHK 461 DO jj = 2, jpjm1 432 462 DO ji = fs_2, fs_jpim1 ! vector opt. 433 463 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 434 464 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 435 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress465 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 436 466 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 437 467 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 438 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( - fsdepw(ji,jj,jk) / htau(ji,jj) ) &439 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)468 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) & 469 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 440 470 END DO 441 471 END DO … … 444 474 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 445 475 ! 446 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer447 CALL wrk_dealloc( jpi,jpj, zhlc )448 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw)476 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 477 CALL wrk_dealloc( jpi,jpj, zhlc ) 478 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 449 479 ! 450 480 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 490 520 INTEGER :: ji, jj, jk ! dummy loop indices 491 521 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 492 REAL(wp) :: zdku, z pdlr, zri, zsqen! - -522 REAL(wp) :: zdku, zri, zsqen ! - - 493 523 REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - - 494 524 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld … … 505 535 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 506 536 ! 537 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 538 zmxlm(:,:,:) = rmxl_min 539 zmxld(:,:,:) = rmxl_min 540 ! 507 541 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 508 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 509 zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) ) 510 ELSE ! surface set to the minimum value 542 DO jj = 2, jpjm1 543 DO ji = fs_2, fs_jpim1 544 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 545 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 546 END DO 547 END DO 548 ELSE 511 549 zmxlm(:,:,1) = rn_mxl0 512 550 ENDIF 513 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 514 ! 515 !CDIR NOVERRCHK 551 ! 516 552 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 517 !CDIR NOVERRCHK 518 DO jj = 2, jpjm1 519 !CDIR NOVERRCHK 553 DO jj = 2, jpjm1 520 554 DO ji = fs_2, fs_jpim1 ! vector opt. 521 555 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) … … 527 561 ! !* Physical limits for the mixing length 528 562 ! 529 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value563 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 530 564 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 531 565 ! 532 566 SELECT CASE ( nn_mxl ) 533 567 ! 568 !!gm Not sure of that coding for ISF.... 569 ! where wmask = 0 set zmxlm == e3w_n 534 570 CASE ( 0 ) ! bounded by the distance to surface and bottom 535 571 DO jk = 2, jpkm1 536 572 DO jj = 2, jpjm1 537 573 DO ji = fs_2, fs_jpim1 ! vector opt. 538 zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), & 539 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 540 zmxlm(ji,jj,jk) = zemxl 541 zmxld(ji,jj,jk) = zemxl 574 zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 575 & gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) ) 576 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 577 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 578 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 542 579 END DO 543 580 END DO … … 548 585 DO jj = 2, jpjm1 549 586 DO ji = fs_2, fs_jpim1 ! vector opt. 550 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )587 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) ) 551 588 zmxlm(ji,jj,jk) = zemxl 552 589 zmxld(ji,jj,jk) = zemxl … … 559 596 DO jj = 2, jpjm1 560 597 DO ji = fs_2, fs_jpim1 ! vector opt. 561 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )598 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 562 599 END DO 563 600 END DO … … 566 603 DO jj = 2, jpjm1 567 604 DO ji = fs_2, fs_jpim1 ! vector opt. 568 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )605 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 569 606 zmxlm(ji,jj,jk) = zemxl 570 607 zmxld(ji,jj,jk) = zemxl … … 577 614 DO jj = 2, jpjm1 578 615 DO ji = fs_2, fs_jpim1 ! vector opt. 579 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )616 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 580 617 END DO 581 618 END DO … … 584 621 DO jj = 2, jpjm1 585 622 DO ji = fs_2, fs_jpim1 ! vector opt. 586 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 587 END DO 588 END DO 589 END DO 590 !CDIR NOVERRCHK 623 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 624 END DO 625 END DO 626 END DO 591 627 DO jk = 2, jpkm1 592 !CDIR NOVERRCHK 593 DO jj = 2, jpjm1 594 !CDIR NOVERRCHK 628 DO jj = 2, jpjm1 595 629 DO ji = fs_2, fs_jpim1 ! vector opt. 596 630 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) … … 612 646 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 613 647 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 614 !CDIR NOVERRCHK615 648 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 616 !CDIR NOVERRCHK 617 DO jj = 2, jpjm1 618 !CDIR NOVERRCHK 649 DO jj = 2, jpjm1 619 650 DO ji = fs_2, fs_jpim1 ! vector opt. 620 651 zsqen = SQRT( en(ji,jj,jk) ) 621 652 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 622 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)623 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)653 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 654 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 624 655 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 625 656 END DO … … 628 659 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 629 660 ! 630 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- andv-points661 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 631 662 DO jj = 2, jpjm1 632 663 DO ji = fs_2, fs_jpim1 ! vector opt. 633 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk)634 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk)664 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 665 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 635 666 END DO 636 667 END DO … … 642 673 DO jj = 2, jpjm1 643 674 DO ji = fs_2, fs_jpim1 ! vector opt. 644 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 645 ! ! shear 646 zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) ) & 647 & + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) ) 648 zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) ) & 649 & + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) ) 650 ! ! local Richardson number 651 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear ) 652 zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) ) 653 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 654 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 655 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 675 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 656 676 # if defined key_c1d 657 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 658 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 677 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 678 !!gm bug NO zri here.... 679 !!gm remove the specific diag for c1d ! 680 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 659 681 # endif 660 682 END DO … … 740 762 ! 741 763 ! !* Check of some namelist values 742 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 743 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 744 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 745 #if ! key_coupled 746 IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 747 #endif 764 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 765 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 766 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 767 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 748 768 749 769 IF( ln_mxl0 ) THEN … … 765 785 ! !* set vertical eddy coef. to the background value 766 786 DO jk = 1, jpk 767 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)768 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)769 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)770 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)787 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 788 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 789 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 790 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 771 791 END DO 772 792 dissl(:,:,:) = 1.e-12_wp … … 819 839 en (:,:,:) = rn_emin * tmask(:,:,:) 820 840 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 841 ! 842 avt_k (:,:,:) = avt (:,:,:) 843 avm_k (:,:,:) = avm (:,:,:) 844 avmu_k(:,:,:) = avmu(:,:,:) 845 avmv_k(:,:,:) = avmv(:,:,:) 846 ! 821 847 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 822 848 ENDIF … … 824 850 en(:,:,:) = rn_emin * tmask(:,:,:) 825 851 DO jk = 1, jpk ! set the Kz to the background value 826 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)827 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)828 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)829 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)852 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 853 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 854 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 855 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 830 856 END DO 831 857 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
