Changeset 5951 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2015-11-30T12:48:01+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r5950 r5951 28 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 29 29 !!---------------------------------------------------------------------- 30 #if defined key_zdftke || defined key_esopa30 #if defined key_zdftke 31 31 !!---------------------------------------------------------------------- 32 32 !! 'key_zdftke' TKE vertical physics … … 53 53 USE timing ! Timing 54 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 55 59 56 60 IMPLICIT NONE … … 85 89 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 86 90 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2]88 91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 89 92 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k ! not enhanced Kz 91 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 92 94 #if defined key_c1d 93 95 ! !!** 1D cfg only ** ('key_c1d') … … 100 102 # include "vectopt_loop_substitute.h90" 101 103 !!---------------------------------------------------------------------- 102 !! NEMO/OPA 4.0 , NEMO Consortium (2011)104 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 103 105 !! $Id$ 104 106 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 115 117 & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & 116 118 #endif 117 & en (jpi,jpj,jpk) , htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 118 & avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk), & 119 & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc ) 119 & htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 120 & apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 120 121 ! 121 122 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) … … 173 174 !!---------------------------------------------------------------------- 174 175 ! 176 #if defined key_agrif 177 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 178 IF( .NOT.Agrif_Root() ) CALL Agrif_Tke 179 #endif 180 ! 175 181 IF( kt /= nit000 ) THEN ! restore before value to compute tke 176 182 avt (:,:,:) = avt_k (:,:,:) … … 189 195 avmv_k(:,:,:) = avmv(:,:,:) 190 196 ! 191 END SUBROUTINE zdf_tke 197 #if defined key_agrif 198 ! Update child grid f => parent grid 199 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 200 #endif 201 ! 202 END SUBROUTINE zdf_tke 192 203 193 204 … … 221 232 REAL(wp) :: zzd_up, zzd_lw ! - - 222 233 !!bfr REAL(wp) :: zebot ! - - 223 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 224 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 225 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 234 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 235 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 236 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 237 REAL(wp) :: zri ! local Richardson number 226 238 !!-------------------------------------------------------------------- 227 239 ! 228 240 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 229 241 ! 230 CALL wrk_alloc( jpi,jpj, imlc ) ! integer231 CALL wrk_alloc( jpi,jpj, zhlc )232 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw)242 CALL wrk_alloc( jpi,jpj, imlc ) ! integer 243 CALL wrk_alloc( jpi,jpj, zhlc ) 244 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 233 245 ! 234 246 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 244 256 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 245 257 DO ji = fs_2, fs_jpim1 ! vector opt. 246 en(ji,jj,mikt(ji,jj)) =rn_emin * tmask(ji,jj,1)258 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 247 259 END DO 248 260 END DO … … 265 277 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 266 278 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 267 !CDIR NOVERRCHK268 279 !! DO jj = 2, jpjm1 269 !CDIR NOVERRCHK270 280 !! DO ji = fs_2, fs_jpim1 ! vector opt. 271 281 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & … … 306 316 END DO 307 317 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 308 !CDIR NOVERRCHK309 318 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 310 !CDIR NOVERRCHK 311 DO jj = 2, jpjm1 312 !CDIR NOVERRCHK 319 DO jj = 2, jpjm1 313 320 DO ji = fs_2, fs_jpim1 ! vector opt. 314 321 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 332 339 ! 333 340 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 334 DO jj = 1, jpj ! here avmu, avmv used as workspace 335 DO ji = 1, jpi 336 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 337 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 338 & / ( fse3uw_n(ji,jj,jk) & 339 & * fse3uw_b(ji,jj,jk) ) 340 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 341 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 342 & / ( fse3vw_n(ji,jj,jk) & 343 & * fse3vw_b(ji,jj,jk) ) 344 END DO 345 END DO 346 END DO 347 ! 341 DO jj = 1, jpjm1 342 DO ji = 1, fs_jpim1 ! vector opt. 343 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 344 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 345 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 346 & / ( fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) ) 347 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 348 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 349 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 350 & / ( fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) ) 351 END DO 352 END DO 353 END DO 354 ! 355 IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr 356 ! Note that zesh2 is also computed in the next loop. 357 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 358 DO jk = 2, jpkm1 359 DO jj = 2, jpjm1 360 DO ji = fs_2, fs_jpim1 ! vector opt. 361 ! ! shear prod. at w-point weightened by mask 362 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 363 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 364 ! ! local Richardson number 365 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 366 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 367 368 END DO 369 END DO 370 END DO 371 ! 372 ENDIF 373 ! 348 374 DO jk = 2, jpkm1 !* Matrix and right hand side in en 349 375 DO jj = 2, jpjm1 … … 354 380 zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal 355 381 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 356 !! shear prod. at w-point weightened by mask357 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &358 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )359 382 ! ! shear prod. at w-point weightened by mask 383 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 384 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 385 ! 360 386 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 361 387 zd_lw(ji,jj,jk) = zzd_lw … … 377 403 END DO 378 404 END DO 379 ! 380 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 381 DO jj = 2, jpjm1 405 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 382 406 DO ji = fs_2, fs_jpim1 ! vector opt. 383 407 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke … … 391 415 END DO 392 416 END DO 393 ! 394 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 395 DO jj = 2, jpjm1 417 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 396 418 DO ji = fs_2, fs_jpim1 ! vector opt. 397 419 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) … … 434 456 END DO 435 457 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 436 !CDIR NOVERRCHK437 458 DO jk = 2, jpkm1 438 !CDIR NOVERRCHK 439 DO jj = 2, jpjm1 440 !CDIR NOVERRCHK 459 DO jj = 2, jpjm1 441 460 DO ji = fs_2, fs_jpim1 ! vector opt. 442 461 ztx2 = utau(ji-1,jj ) + utau(ji,jj) … … 453 472 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 454 473 ! 455 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer456 CALL wrk_dealloc( jpi,jpj, zhlc )457 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw)474 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 475 CALL wrk_dealloc( jpi,jpj, zhlc ) 476 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 458 477 ! 459 478 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 499 518 INTEGER :: ji, jj, jk ! dummy loop indices 500 519 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 501 REAL(wp) :: zdku, z pdlr, zri, zsqen! - -520 REAL(wp) :: zdku, zri, zsqen ! - - 502 521 REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - - 503 522 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld … … 529 548 ENDIF 530 549 ! 531 !CDIR NOVERRCHK532 550 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 533 !CDIR NOVERRCHK 534 DO jj = 2, jpjm1 535 !CDIR NOVERRCHK 551 DO jj = 2, jpjm1 536 552 DO ji = fs_2, fs_jpim1 ! vector opt. 537 553 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 538 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )554 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 539 555 END DO 540 556 END DO … … 543 559 ! !* Physical limits for the mixing length 544 560 ! 545 zmxld(:,:, 1) = zmxlm(:,:,1) ! surface set to the minimum value561 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 546 562 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 547 563 ! 548 564 SELECT CASE ( nn_mxl ) 549 565 ! 566 !!gm Not sure of that coding for ISF.... 550 567 ! where wmask = 0 set zmxlm == fse3w 551 568 CASE ( 0 ) ! bounded by the distance to surface and bottom … … 606 623 END DO 607 624 END DO 608 !CDIR NOVERRCHK609 625 DO jk = 2, jpkm1 610 !CDIR NOVERRCHK 611 DO jj = 2, jpjm1 612 !CDIR NOVERRCHK 626 DO jj = 2, jpjm1 613 627 DO ji = fs_2, fs_jpim1 ! vector opt. 614 628 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) … … 630 644 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 631 645 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 632 !CDIR NOVERRCHK633 646 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 634 !CDIR NOVERRCHK 635 DO jj = 2, jpjm1 636 !CDIR NOVERRCHK 647 DO jj = 2, jpjm1 637 648 DO ji = fs_2, fs_jpim1 ! vector opt. 638 649 zsqen = SQRT( en(ji,jj,jk) ) … … 660 671 DO jj = 2, jpjm1 661 672 DO ji = fs_2, fs_jpim1 ! vector opt. 662 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 663 ! ! shear 664 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) ) & 665 & + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) ) 666 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) ) & 667 & + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) ) 668 ! ! local Richardson number 669 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear ) 670 zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) ) 671 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 672 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 673 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 673 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 674 674 # if defined key_c1d 675 e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 676 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 675 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 676 !!gm bug NO zri here.... 677 !!gm remove the specific diag for c1d ! 678 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 677 679 # endif 678 680 END DO
Note: See TracChangeset
for help on using the changeset viewer.