- Timestamp:
- 2017-05-30T10:13:14+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7991 r8093 43 43 USE sbc_oce ! surface boundary condition: ocean 44 44 USE zdf_oce ! vertical physics: ocean variables 45 USE zdfsh2 ! vertical physics: shear production term of TKE 45 !!gm new 46 USE zdfdrg ! vertical physics: top/bottom drag coef. 47 !!gm 46 48 USE zdfmxl ! vertical physics: mixed layer 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link)48 USE prtctl ! Print control49 USE in_out_manager ! I/O manager50 USE iom ! I/O manager library51 USE lib_mpp ! MPP library52 USE wrk_nemo ! work arrays53 USE timing ! Timing54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)55 49 #if defined key_agrif 56 50 USE agrif_opa_interp 57 51 USE agrif_opa_update 58 52 #endif 53 ! 54 USE in_out_manager ! I/O manager 55 USE iom ! I/O manager library 56 USE lib_mpp ! MPP library 57 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 58 USE prtctl ! Print control 59 USE wrk_nemo ! work arrays 60 USE timing ! Timing 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 59 62 60 63 IMPLICIT NONE … … 87 90 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 88 91 89 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau! depth of tke penetration (nn_htau)90 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl! now mixing lenght of dissipation91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr! now mixing lenght of dissipation92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 92 95 93 96 !! * Substitutions … … 112 115 113 116 114 SUBROUTINE zdf_tke( kt )117 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 115 118 !!---------------------------------------------------------------------- 116 119 !! *** ROUTINE zdf_tke *** … … 157 160 !! Bruchard OM 2002 158 161 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: kt ! ocean time step 162 INTEGER , INTENT(in ) :: kt ! ocean time step 163 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 164 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 160 165 !!---------------------------------------------------------------------- 161 166 ! … … 165 170 #endif 166 171 ! 167 avt(:,:,:) = avt_k(:,:,:) ! restore before Kz computed with TKE alone 168 avm(:,:,:) = avm_k(:,:,:) 169 ! 170 CALL tke_tke ! now tke (en) 171 ! 172 CALL tke_avn ! now avt, avm, dissl 173 ! 174 avt_k(:,:,:) = avt(:,:,:) ! store Kz computed with TKE alone 175 avm_k(:,:,:) = avm(:,:,:) 172 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en) 173 ! 174 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl 176 175 ! 177 176 #if defined key_agrif … … 179 178 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 180 179 #endif 181 !182 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' )183 180 ! 184 181 END SUBROUTINE zdf_tke 185 182 186 183 187 SUBROUTINE tke_tke 184 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2 & 185 & , p_avm, p_avt ) 188 186 !!---------------------------------------------------------------------- 189 187 !! *** ROUTINE tke_tke *** … … 200 198 !! ** Action : - en : now turbulent kinetic energy) 201 199 !! --------------------------------------------------------------------- 202 INTEGER :: ji, jj, jk ! dummy loop arguments 203 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! local integers 204 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! - - 205 !!bfr REAL(wp) :: zebot ! local scalars 200 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 201 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 202 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 203 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 204 ! 205 INTEGER :: ji, jj, jk ! dummy loop arguments 206 !!bfr REAL(wp) :: zebot, zmshu, zmskv ! local scalars 206 207 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 207 208 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient … … 212 213 REAL(wp) :: zus , zwlc , zind ! - - 213 214 REAL(wp) :: zzd_up, zzd_lw ! - - 214 INTEGER , DIMENSION(jpi,jpj) :: imlc 215 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 216 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 215 INTEGER , DIMENSION(jpi,jpj) :: imlc 216 REAL(wp), DIMENSION(jpi,jpj) :: zhlc 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 218 218 !!-------------------------------------------------------------------- 219 219 ! 220 220 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 221 !222 CALL wrk_alloc( jpi,jpj, zhlc )223 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )224 221 ! 225 222 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 256 253 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 257 254 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 255 !!gm old 258 256 !! DO jj = 2, jpjm1 259 257 !! DO ji = fs_2, fs_jpim1 ! vector opt. … … 263 261 !! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 264 262 !! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 265 !! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)263 !! en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 266 264 !! END DO 267 265 !! END DO 266 !!gm new 267 !! 268 !! ====>>>> add below an wet-only calculation of u and v at t-point like in zdfsh2: 269 !! zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 270 !! zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 271 !! 272 !! 273 !! DO jj = 2, jpjm1 274 !! DO ji = fs_2, fs_jpim1 ! vector opt. 275 !! zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 276 !! zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 277 !! ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 278 !! zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 279 !! & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 280 !! en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 281 !! END DO 282 !! END DO 283 !! IF( ln_isfcav ) THEN !top friction 284 !! DO jj = 2, jpjm1 285 !! DO ji = fs_2, fs_jpim1 ! vector opt. 286 !! zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 287 !! zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 288 !! ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 289 !! zebot = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 290 !! & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 291 !! en(ji,jj,mikt(ji,jj)+1) = MAX( zebot, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 292 !! END DO 293 !! END DO 294 !! ENDIF 295 !! 268 296 !!bfr - end commented area 269 297 ! 270 298 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 271 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke 299 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke ! (Axell JGR 2002) 272 300 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 273 301 ! 274 302 ! !* total energy produce by LC : cumulative sum over jk 275 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)303 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 276 304 DO jk = 2, jpk 277 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)305 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 278 306 END DO 279 307 ! !* finite Langmuir Circulation depth … … 291 319 DO jj = 1, jpj 292 320 DO ji = 1, jpi 293 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))321 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 294 322 END DO 295 323 END DO … … 300 328 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 301 329 ! ! vertical velocity due to LC 302 zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )303 zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )330 zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 331 zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 304 332 ! ! TKE Langmuir circulation source term 305 333 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) & … … 318 346 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 319 347 ! 320 ! !* Shear production at uw- and vw-points (energy conserving form)321 CALL zdf_sh2( zsh2 )322 !323 348 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 324 349 DO jk = 2, jpkm1 … … 326 351 DO ji = 2, jpim1 327 352 ! ! local Richardson number 328 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear )353 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 329 354 ! ! inverse of Prandtl number 330 355 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) … … 340 365 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 341 366 ! ! eddy coefficient (ensure numerical stability) 342 zzd_up = zcof * MAX( avm(ji,jj,jk+1) +avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal343 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) )344 zzd_lw = zcof * MAX( avm(ji,jj,jk ) +avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal345 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) )367 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 368 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 369 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 370 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 346 371 ! 347 372 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 350 375 ! 351 376 ! ! right hand side in en 352 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zsh2(ji,jj,jk)& ! shear353 & - avt(ji,jj,jk) * rn2(ji,jj,jk)& ! stratification377 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 378 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 354 379 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 355 380 & ) * wmask(ji,jj,jk) … … 407 432 DO jj = 2, jpjm1 408 433 DO ji = fs_2, fs_jpim1 ! vector opt. 409 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &434 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 410 435 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 411 436 END DO … … 416 441 DO ji = fs_2, fs_jpim1 ! vector opt. 417 442 jk = nmln(ji,jj) 418 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &443 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 419 444 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 420 445 END DO … … 429 454 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 430 455 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 431 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &456 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 432 457 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 433 458 END DO … … 435 460 END DO 436 461 ENDIF 437 !!gm not sure we need this lbc ....438 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)439 !440 CALL wrk_dealloc( jpi,jpj, zhlc )441 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )442 462 ! 443 463 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 446 466 447 467 448 SUBROUTINE tke_avn 468 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 449 469 !!---------------------------------------------------------------------- 450 470 !! *** ROUTINE tke_avn *** … … 480 500 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 481 501 !!---------------------------------------------------------------------- 502 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 503 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 504 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 505 ! 482 506 INTEGER :: ji, jj, jk ! dummy loop indices 483 507 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 484 508 REAL(wp) :: zdku, zdkv, zsqen ! - - 485 509 REAL(wp) :: zemxl, zemlm, zemlp ! - - 486 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl,zmxlm, zmxld510 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxlm, zmxld 487 511 !!-------------------------------------------------------------------- 488 512 ! 489 513 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 490 491 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )492 514 493 515 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 502 524 ! 503 525 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 526 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 504 527 DO jj = 2, jpjm1 505 528 DO ji = fs_2, fs_jpim1 506 zraug = vkarmn * 2.e5_wp / ( rau0 * grav )507 529 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 508 530 END DO … … 529 551 ! 530 552 !!gm Not sure of that coding for ISF.... 531 ! where wmask = 0 set zmxlm == e3w_n553 ! where wmask = 0 set zmxlm == p_e3w 532 554 CASE ( 0 ) ! bounded by the distance to surface and bottom 533 555 DO jk = 2, jpkm1 534 556 DO jj = 2, jpjm1 535 557 DO ji = fs_2, fs_jpim1 ! vector opt. 536 zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &537 & gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )558 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 559 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 538 560 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 539 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))540 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))561 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 562 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 541 563 END DO 542 564 END DO … … 547 569 DO jj = 2, jpjm1 548 570 DO ji = fs_2, fs_jpim1 ! vector opt. 549 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )571 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 550 572 zmxlm(ji,jj,jk) = zemxl 551 573 zmxld(ji,jj,jk) = zemxl … … 558 580 DO jj = 2, jpjm1 559 581 DO ji = fs_2, fs_jpim1 ! vector opt. 560 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )582 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 561 583 END DO 562 584 END DO … … 565 587 DO jj = 2, jpjm1 566 588 DO ji = fs_2, fs_jpim1 ! vector opt. 567 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )589 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 568 590 zmxlm(ji,jj,jk) = zemxl 569 591 zmxld(ji,jj,jk) = zemxl … … 576 598 DO jj = 2, jpjm1 577 599 DO ji = fs_2, fs_jpim1 ! vector opt. 578 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )600 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 579 601 END DO 580 602 END DO … … 583 605 DO jj = 2, jpjm1 584 606 DO ji = fs_2, fs_jpim1 ! vector opt. 585 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )607 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 586 608 END DO 587 609 END DO … … 609 631 zsqen = SQRT( en(ji,jj,jk) ) 610 632 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 611 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)612 avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)633 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 634 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 613 635 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 614 636 END DO … … 621 643 DO jj = 2, jpjm1 622 644 DO ji = fs_2, fs_jpim1 ! vector opt. 623 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) *avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)645 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 624 646 END DO 625 647 END DO 626 648 END DO 627 649 ENDIF 628 !!gm I'm sure this is needed to compute the shear term at next time-step629 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)630 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)631 650 632 651 IF(ln_ctl) THEN … … 634 653 CALL prt_ctl( tab3d_1=avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 635 654 ENDIF 636 !637 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )638 655 ! 639 656 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 737 754 ! !* set vertical eddy coef. to the background value 738 755 DO jk = 1, jpk 739 avt (:,:,jk) = avtb(jk) * wmask(:,:,jk)740 avm (:,:,jk) = avmb(jk) * wmask(:,:,jk)756 avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 757 avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 741 758 END DO 742 759 dissl(:,:,:) = 1.e-12_wp … … 748 765 749 766 SUBROUTINE tke_rst( kt, cdrw ) 750 !!--------------------------------------------------------------------- 751 !! *** ROUTINE tke_rst *** 752 !! 753 !! ** Purpose : Read or write TKE file (en) in restart file 754 !! 755 !! ** Method : use of IOM library 756 !! if the restart does not contain TKE, en is either 757 !! set to rn_emin or recomputed 758 !!---------------------------------------------------------------------- 759 INTEGER , INTENT(in) :: kt ! ocean time-step 760 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 761 ! 762 INTEGER :: jit, jk ! dummy loop indices 763 INTEGER :: id1, id2, id3, id4 ! local integers 764 !!---------------------------------------------------------------------- 765 ! 766 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 767 ! ! --------------- 768 IF( ln_rstart ) THEN !* Read the restart file 769 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 770 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 771 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 772 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 773 ! 774 IF( id1 > 0 ) THEN ! 'en' exists 775 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 776 IF( MIN( id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 777 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 778 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 779 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 780 ELSE ! one at least array is missing 781 CALL tke_avn ! compute avt, avm, and dissl (approximation) 782 ENDIF 783 ELSE ! No TKE array found: initialisation 784 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 785 en (:,:,:) = rn_emin * wmask(:,:,:) 786 CALL tke_avn ! recompute avt, avm, and dissl (approximation) 787 avt_k(:,:,:) = avt(:,:,:) 788 avm_k(:,:,:) = avm(:,:,:) 789 ! 790 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 791 avt_k(:,:,:) = avt(:,:,:) 792 avm_k(:,:,:) = avm(:,:,:) 793 ENDIF 794 ELSE !* Start from rest 795 en(:,:,:) = rn_emin * tmask(:,:,:) 796 DO jk = 1, jpk ! set the Kz to the background value 797 avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk) 798 avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk) 799 END DO 800 ENDIF 801 ! 802 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 803 ! ! ------------------- 804 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 805 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 806 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 807 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 808 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 809 ! 810 ENDIF 811 ! 767 !!--------------------------------------------------------------------- 768 !! *** ROUTINE tke_rst *** 769 !! 770 !! ** Purpose : Read or write TKE file (en) in restart file 771 !! 772 !! ** Method : use of IOM library 773 !! if the restart does not contain TKE, en is either 774 !! set to rn_emin or recomputed 775 !!---------------------------------------------------------------------- 776 INTEGER , INTENT(in) :: kt ! ocean time-step 777 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 778 ! 779 INTEGER :: jit, jk ! dummy loop indices 780 INTEGER :: id1, id2, id3, id4 ! local integers 781 !!---------------------------------------------------------------------- 782 ! 783 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 784 ! ! --------------- 785 IF( ln_rstart ) THEN !* Read the restart file 786 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 787 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 788 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 789 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 790 ! 791 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 792 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 793 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 794 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 795 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 796 ELSE ! start TKE from rest 797 IF(lwp) WRITE(numout,*) ' ==>> previous run without TKE scheme, set en to background values' 798 en(:,:,:) = rn_emin * wmask(:,:,:) 799 ! avt_k, avm_k already set to the background value in zdf_phy_init 800 ENDIF 801 ELSE !* Start from rest 802 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set en to the background value' 803 en(:,:,:) = rn_emin * wmask(:,:,:) 804 ! avt_k, avm_k already set to the background value in zdf_phy_init 805 ENDIF 806 ! 807 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 808 ! ! ------------------- 809 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 810 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 811 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 812 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 813 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 814 ! 815 ENDIF 816 ! 812 817 END SUBROUTINE tke_rst 813 818
Note: See TracChangeset
for help on using the changeset viewer.
