- Timestamp:
- 2017-05-01T16:21:42+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
r7990 r7991 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 46 USE zdfmxl ! vertical physics: mixed layer 46 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 89 90 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 90 91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 91 #if defined key_c1d92 ! !!** 1D cfg only ** ('key_c1d')93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers95 #endif96 92 97 93 !! * Substitutions … … 108 104 !! *** FUNCTION zdf_tke_alloc *** 109 105 !!---------------------------------------------------------------------- 110 ALLOCATE( & 111 #if defined key_c1d 112 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & 113 & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & 114 #endif 115 & htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 116 & apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 106 ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 117 107 ! 118 108 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) … … 211 201 !! --------------------------------------------------------------------- 212 202 INTEGER :: ji, jj, jk ! dummy loop arguments 213 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar214 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar215 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 216 REAL(wp) :: z cdrag = 1.5e-3 ! drag coefficient217 REAL(wp) :: z bbrau, zesh2 ! temporary scalars218 REAL(wp) :: z fact1, zfact2, zfact3 ! - -219 REAL(wp) :: z tx2 , zty2 , zcof !- -220 REAL(wp) :: zt au , zdif !- -221 REAL(wp) :: z us , zwlc , zind !- -222 REAL(wp) :: z zd_up, zzd_lw !- -223 !!bfr REAL(wp) :: zebot !- -203 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! local integers 204 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! - - 205 !!bfr REAL(wp) :: zebot ! local scalars 206 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 207 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 208 REAL(wp) :: zbbrau, zri ! local scalars 209 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 210 REAL(wp) :: ztx2 , zty2 , zcof ! - - 211 REAL(wp) :: ztau , zdif ! - - 212 REAL(wp) :: zus , zwlc , zind ! - - 213 REAL(wp) :: zzd_up, zzd_lw ! - - 224 214 INTEGER , DIMENSION(jpi,jpj) :: imlc 225 215 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 226 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw , z3du, z3dv227 REAL(wp) :: zri ! local Richardson number216 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 228 218 !!-------------------------------------------------------------------- 229 219 ! … … 231 221 ! 232 222 CALL wrk_alloc( jpi,jpj, zhlc ) 233 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw , z3du, z3dv)223 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 234 224 ! 235 225 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 328 318 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 329 319 ! 330 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 331 DO jj = 1, jpjm1 332 DO ji = 1, fs_jpim1 ! vector opt. 333 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 334 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 335 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 336 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 337 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 338 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 339 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 340 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 341 END DO 342 END DO 343 END DO 344 ! 345 IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr 346 ! Note that zesh2 is also computed in the next loop. 347 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 320 ! !* Shear production at uw- and vw-points (energy conserving form) 321 CALL zdf_sh2( zsh2 ) 322 ! 323 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 348 324 DO jk = 2, jpkm1 349 325 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 351 ! ! shear prod. at w-point weightened by mask 352 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 353 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 354 ! ! local Richardson number 355 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 326 DO ji = 2, jpim1 327 ! ! local Richardson number 328 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear ) 329 ! ! inverse of Prandtl number 356 330 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 357 358 END DO 359 END DO 360 END DO 361 ! 331 END DO 332 END DO 333 END DO 362 334 ENDIF 363 335 ! … … 373 345 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 374 346 ! 375 ! ! shear prod. at w-point weightened by mask376 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &377 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )378 !379 347 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 380 348 zd_lw(ji,jj,jk) = zzd_lw 381 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)349 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 382 350 ! 383 351 ! ! right hand side in en 384 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 385 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 386 & * wmask(ji,jj,jk) 352 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zsh2(ji,jj,jk) & ! shear 353 & - avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 354 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 355 & ) * wmask(ji,jj,jk) 387 356 END DO 388 357 END DO … … 470 439 ! 471 440 CALL wrk_dealloc( jpi,jpj, zhlc ) 472 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw , z3du, z3dv)441 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 473 442 ! 474 443 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 631 600 END SELECT 632 601 ! 633 # if defined key_c1d634 e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales635 e_mix(:,:,:) = zmxlm(:,:,:)636 # endif637 602 638 603 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 657 622 DO ji = fs_2, fs_jpim1 ! vector opt. 658 623 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 659 # if defined key_c1d660 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number661 # endif662 624 END DO 663 625 END DO
Note: See TracChangeset
for help on using the changeset viewer.
