- Timestamp:
- 2017-05-20T13:49:38+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7991 r8055 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 TKE46 45 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 46 #if defined key_agrif 56 47 USE agrif_opa_interp 57 48 USE agrif_opa_update 58 49 #endif 50 ! 51 USE in_out_manager ! I/O manager 52 USE iom ! I/O manager library 53 USE lib_mpp ! MPP library 54 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 55 USE prtctl ! Print control 56 USE wrk_nemo ! work arrays 57 USE timing ! Timing 58 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 59 59 60 60 IMPLICIT NONE … … 87 87 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 88 88 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 dissipation89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 92 92 93 93 !! * Substitutions 94 # include " vectopt_loop_substitute.h90"94 # include "domain_substitute.h90" 95 95 !!---------------------------------------------------------------------- 96 96 !! NEMO/OPA 3.7 , NEMO Consortium (2015) … … 112 112 113 113 114 SUBROUTINE zdf_tke( kt )114 SUBROUTINE zdf_tke( ARG_2D, kt, p_sh2, p_avm, p_avt ) 115 115 !!---------------------------------------------------------------------- 116 116 !! *** ROUTINE zdf_tke *** … … 157 157 !! Bruchard OM 2002 158 158 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: kt ! ocean time step 159 INTEGER , INTENT(in ) :: ARG_2D ! inner domain start-end i-indices 160 INTEGER , INTENT(in ) :: kt ! ocean time step 161 REAL(wp), DIMENSION(WRK_3D), INTENT(in ) :: p_sh2 ! shear production term 162 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 160 163 !!---------------------------------------------------------------------- 161 164 ! … … 165 168 #endif 166 169 ! 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(:,:,:) 170 CALL tke_tke( ARG_2D, gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en) 171 ! 172 173 !!gm not sure we need this lbc .... ==>>> certainly NOT !!! 174 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 175 ! 176 CALL tke_avn( ARG_2D, gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl 176 177 ! 177 178 #if defined key_agrif … … 179 180 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 180 181 #endif 181 !182 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' )183 182 ! 184 183 END SUBROUTINE zdf_tke 185 184 186 185 187 SUBROUTINE tke_tke 186 SUBROUTINE tke_tke( ARG_2D, pdepw, p_e3t, p_e3w, p_sh2 & 187 & , p_avm, p_avt ) 188 188 !!---------------------------------------------------------------------- 189 189 !! *** ROUTINE tke_tke *** … … 200 200 !! ** Action : - en : now turbulent kinetic energy) 201 201 !! --------------------------------------------------------------------- 202 INTEGER :: ji, jj, jk ! dummy loop arguments 202 INTEGER , INTENT(in ) :: ARG_2D ! inner domain start-end i-indices 203 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! vertical eddy viscosity (w-points) 204 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! vertical eddy viscosity (w-points) 205 REAL(wp), DIMENSION(WRK_3D), INTENT(in ) :: p_sh2 ! shear production term 206 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity (w-points) 207 ! 208 INTEGER :: ji, jj, jk ! dummy loop arguments 203 209 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! local integers 204 210 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! - - 205 211 !!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 ! - - 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 212 REAL(wp):: zrhoa = 1.22 ! Air density kg/m3 213 REAL(wp):: zcdrag = 1.5e-3 ! drag coefficient 214 REAL(wp):: zbbrau, zri ! local scalars 215 REAL(wp):: zfact1, zfact2, zfact3 ! - - 216 REAL(wp):: ztx2 , zty2 , zcof ! - - 217 REAL(wp):: ztau , zdif ! - - 218 REAL(wp):: zus , zwlc , zind ! - - 219 REAL(wp):: zzd_up, zzd_lw ! - - 220 INTEGER , DIMENSION(WRK_2D) :: imlc 221 REAL(wp), DIMENSION(WRK_2D) :: zhlc 222 REAL(wp), DIMENSION(WRK_3D) :: zpelc, zdiag, zd_up, zd_lw 218 223 !!-------------------------------------------------------------------- 219 224 ! 220 225 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 226 ! 225 227 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 233 235 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 234 236 IF ( ln_isfcav ) THEN 235 DO jj = 2, jpjm1! en(mikt(ji,jj)) = rn_emin236 DO ji = fs_2, fs_jpim1 ! vector opt.237 DO jj = k_Jstr, k_Jend ! en(mikt(ji,jj)) = rn_emin 238 DO ji = k_Istr, k_Iend 237 239 en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 238 240 END DO 239 241 END DO 240 242 END IF 241 DO jj = 2, jpjm1! en(1) = rn_ebb taum / rau0 (min value rn_emin0)242 DO ji = fs_2, fs_jpim1 ! vector opt.243 DO jj = k_Jstr, k_Jend ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 244 DO ji = k_Istr, k_Iend 243 245 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 244 246 END DO … … 256 258 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 257 259 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 258 !! DO jj = 2, jpjm1259 !! DO ji = fs_2, fs_jpim1 ! vector opt.260 !! DO jj = k_Jstr, k_Jend 261 !! DO ji = k_Jstr, k_Iend 260 262 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 261 263 !! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) ) … … 269 271 ! 270 272 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 271 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke 273 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke ! (Axell JGR 2002) 272 274 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 273 275 ! 274 276 ! !* total energy produce by LC : cumulative sum over jk 275 zpelc( :,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)277 zpelc(WRK_2D,1) = MAX( rn2b(WRK_2D,1), 0._wp ) * pdepw(WRK_2D,1) * p_e3w(WRK_2D,1) 276 278 DO jk = 2, jpk 277 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk) 279 zpelc(WRK_2D,jk) = zpelc(WRK_2D,jk-1) & 280 & + MAX( rn2b(WRK_2D,jk), 0._wp ) & 281 & * pdepw(WRK_2D,jk) * p_e3w(WRK_2D,jk) 278 282 END DO 279 283 ! !* finite Langmuir Circulation depth 280 284 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 281 imlc( :,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land)285 imlc(WRK_2D) = mbkt(WRK_2D) + 1 ! Initialization to the number of w ocean point (=2 over land) 282 286 DO jk = jpkm1, 2, -1 283 DO jj = 1, jpj! Last w-level at which zpelc>=0.5*us*us284 DO ji = 1, jpi! with us=0.016*wind(starting from jpk-1)287 DO jj = k_Jstr, k_Jend ! Last w-level at which zpelc>=0.5*us*us 288 DO ji = k_Istr, k_Iend ! with us=0.016*wind(starting from jpk-1) 285 289 zus = zcof * taum(ji,jj) 286 290 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk … … 289 293 END DO 290 294 ! ! finite LC depth 291 DO jj = 1, jpj292 DO ji = 1, jpi293 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))295 DO jj = k_Jstr, k_Jend 296 DO ji = k_Istr, k_Iend 297 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 294 298 END DO 295 299 END DO 296 300 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 297 301 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 298 DO jj = 2, jpjm1299 DO ji = fs_2, fs_jpim1 ! vector opt.302 DO jj = k_Jstr, k_Jend 303 DO ji = k_Istr, k_Iend 300 304 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 301 305 ! ! 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) )306 zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 307 zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 304 308 ! ! TKE Langmuir circulation source term 305 309 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) & … … 318 322 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 319 323 ! 320 ! !* Shear production at uw- and vw-points (energy conserving form)321 CALL zdf_sh2( zsh2 )322 !323 324 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 324 325 DO jk = 2, jpkm1 325 DO jj = 2, jpjm1326 DO ji = 2, jpim1326 DO jj = k_Jstr, k_Jend 327 DO ji = k_Istr, k_Iend 327 328 ! ! local Richardson number 328 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear )329 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 329 330 ! ! inverse of Prandtl number 330 331 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) … … 335 336 ! 336 337 DO jk = 2, jpkm1 !* Matrix and right hand side in en 337 DO jj = 2, jpjm1338 DO ji = fs_2, fs_jpim1 ! vector opt.338 DO jj = k_Jstr, k_Jend 339 DO ji = k_Istr, k_Iend 339 340 zcof = zfact1 * tmask(ji,jj,jk) 340 341 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 341 342 ! ! 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 ) )343 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 344 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 345 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 346 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 346 347 ! 347 348 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 350 351 ! 351 352 ! ! 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)& ! stratification353 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 354 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 354 355 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 355 356 & ) * wmask(ji,jj,jk) … … 359 360 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 360 361 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 361 DO jj = 2, jpjm1362 DO ji = fs_2, fs_jpim1 ! vector opt.362 DO jj = k_Jstr, k_Jend 363 DO ji = k_Istr, k_Iend 363 364 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 364 365 END DO 365 366 END DO 366 367 END DO 367 DO jj = 2, jpjm1! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1368 DO ji = fs_2, fs_jpim1 ! vector opt.368 DO jj = k_Jstr, k_Jend ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 369 DO ji = k_Istr, k_Iend 369 370 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 370 371 END DO 371 372 END DO 372 373 DO jk = 3, jpkm1 373 DO jj = 2, jpjm1374 DO ji = fs_2, fs_jpim1 ! vector opt.374 DO jj = k_Jstr, k_Jend 375 DO ji = k_Istr, k_Iend 375 376 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) 376 377 END DO 377 378 END DO 378 379 END DO 379 DO jj = 2, jpjm1! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk380 DO ji = fs_2, fs_jpim1 ! vector opt.380 DO jj = k_Jstr, k_Jend ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 381 DO ji = k_Istr, k_Iend 381 382 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 382 383 END DO 383 384 END DO 384 385 DO jk = jpk-2, 2, -1 385 DO jj = 2, jpjm1386 DO ji = fs_2, fs_jpim1 ! vector opt.386 DO jj = k_Jstr, k_Jend 387 DO ji = k_Istr, k_Iend 387 388 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 388 389 END DO … … 390 391 END DO 391 392 DO jk = 2, jpkm1 ! set the minimum value of tke 392 DO jj = 2, jpjm1393 DO ji = fs_2, fs_jpim1 ! vector opt.393 DO jj = k_Jstr, k_Jend 394 DO ji = k_Istr, k_Iend 394 395 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 395 396 END DO … … 405 406 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 406 407 DO jk = 2, jpkm1 407 DO jj = 2, jpjm1408 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) ) &408 DO jj = k_Jstr, k_Jend 409 DO ji = k_Istr, k_Iend 410 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 410 411 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 411 412 END DO … … 413 414 END DO 414 415 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 415 DO jj = 2, jpjm1416 DO ji = fs_2, fs_jpim1 ! vector opt.416 DO jj = k_Jstr, k_Jend 417 DO ji = k_Istr, k_Iend 417 418 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) ) &419 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 419 420 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 420 421 END DO … … 422 423 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 423 424 DO jk = 2, jpkm1 424 DO jj = 2, jpjm1425 DO ji = fs_2, fs_jpim1 ! vector opt.425 DO jj = k_Jstr, k_Jend 426 DO ji = k_Istr, k_Iend 426 427 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 427 428 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) … … 429 430 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 430 431 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) ) &432 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 432 433 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 433 434 END DO … … 435 436 END DO 436 437 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 438 ! 443 439 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 446 442 447 443 448 SUBROUTINE tke_avn 444 SUBROUTINE tke_avn( ARG_2D, pdepw, p_e3t, p_e3w, p_avm, p_avt ) 449 445 !!---------------------------------------------------------------------- 450 446 !! *** ROUTINE tke_avn *** … … 480 476 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 481 477 !!---------------------------------------------------------------------- 478 INTEGER , INTENT(in ) :: ARG_2D ! inner domain start-end i-indices 479 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! vertical eddy viscosity (w-points) 480 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! vertical eddy viscosity (w-points) 481 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity (w-points) 482 ! 482 483 INTEGER :: ji, jj, jk ! dummy loop indices 483 484 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 484 485 REAL(wp) :: zdku, zdkv, zsqen ! - - 485 486 REAL(wp) :: zemxl, zemlm, zemlp ! - - 486 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl,zmxlm, zmxld487 REAL(wp), DIMENSION(WRK_3D) :: zmxlm, zmxld 487 488 !!-------------------------------------------------------------------- 488 489 ! 489 490 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 490 491 491 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )492 493 492 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 494 493 ! ! Mixing length … … 498 497 ! 499 498 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 500 zmxlm( :,:,:) = rmxl_min501 zmxld( :,:,:) = rmxl_min499 zmxlm(WRK_3D) = rmxl_min 500 zmxld(WRK_3D) = rmxl_min 502 501 ! 503 502 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 504 DO jj = 2, jpjm1505 DO ji = fs_2, fs_jpim1506 zraug = vkarmn * 2.e5_wp / ( rau0 * grav )503 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 504 DO jj = k_Jstr, k_Jend 505 DO ji = k_Istr, k_Iend 507 506 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 508 507 END DO 509 508 END DO 510 509 ELSE 511 zmxlm( :,:,1) = rn_mxl0510 zmxlm(WRK_2D,1) = rn_mxl0 512 511 ENDIF 513 512 ! 514 513 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 515 DO jj = 2, jpjm1516 DO ji = fs_2, fs_jpim1 ! vector opt.514 DO jj = k_Jstr, k_Jend 515 DO ji = k_Istr, k_Iend 517 516 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 518 517 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 523 522 ! !* Physical limits for the mixing length 524 523 ! 525 zmxld( :,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value526 zmxld( :,:,jpk) = rmxl_min! last level set to the minimum value524 zmxld(WRK_2D, 1 ) = zmxlm(WRK_2D,1) ! surface set to the minimum value 525 zmxld(WRK_2D,jpk) = rmxl_min ! last level set to the minimum value 527 526 ! 528 527 SELECT CASE ( nn_mxl ) 529 528 ! 530 529 !!gm Not sure of that coding for ISF.... 531 ! where wmask = 0 set zmxlm == e3w_n530 ! where wmask = 0 set zmxlm == p_e3w 532 531 CASE ( 0 ) ! bounded by the distance to surface and bottom 533 532 DO jk = 2, jpkm1 534 DO jj = 2, jpjm1535 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) )533 DO jj = k_Jstr, k_Jend 534 DO ji = k_Istr, k_Iend 535 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 536 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 538 537 ! 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))538 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 539 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 540 END DO 542 541 END DO … … 545 544 CASE ( 1 ) ! bounded by the vertical scale factor 546 545 DO jk = 2, jpkm1 547 DO jj = 2, jpjm1548 DO ji = fs_2, fs_jpim1 ! vector opt.549 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )546 DO jj = k_Jstr, k_Jend 547 DO ji = k_Istr, k_Iend 548 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 550 549 zmxlm(ji,jj,jk) = zemxl 551 550 zmxld(ji,jj,jk) = zemxl … … 556 555 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 557 556 DO jk = 2, jpkm1 ! from the surface to the bottom : 558 DO jj = 2, jpjm1559 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) )557 DO jj = k_Jstr, k_Jend 558 DO ji = k_Istr, k_Iend 559 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 561 560 END DO 562 561 END DO 563 562 END DO 564 563 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 565 DO jj = 2, jpjm1566 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) )564 DO jj = k_Jstr, k_Jend 565 DO ji = k_Istr, k_Iend 566 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 568 567 zmxlm(ji,jj,jk) = zemxl 569 568 zmxld(ji,jj,jk) = zemxl … … 574 573 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 575 574 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 576 DO jj = 2, jpjm1577 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) )575 DO jj = k_Jstr, k_Jend 576 DO ji = k_Istr, k_Iend 577 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 579 578 END DO 580 579 END DO 581 580 END DO 582 581 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 583 DO jj = 2, jpjm1584 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) )582 DO jj = k_Jstr, k_Jend 583 DO ji = k_Istr, k_Iend 584 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 586 585 END DO 587 586 END DO 588 587 END DO 589 588 DO jk = 2, jpkm1 590 DO jj = 2, jpjm1591 DO ji = fs_2, fs_jpim1 ! vector opt.589 DO jj = k_Jstr, k_Jend 590 DO ji = k_Istr, k_Iend 592 591 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 593 592 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) … … 605 604 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 606 605 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 607 DO jj = 2, jpjm1608 DO ji = fs_2, fs_jpim1 ! vector opt.606 DO jj = k_Jstr, k_Jend 607 DO ji = k_Istr, k_Iend 609 608 zsqen = SQRT( en(ji,jj,jk) ) 610 609 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)610 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 611 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 613 612 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 614 613 END DO … … 619 618 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 620 619 DO jk = 2, jpkm1 621 DO jj = 2, jpjm1622 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)620 DO jj = k_Jstr, k_Jend 621 DO ji = k_Istr, k_Iend 622 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 623 END DO 625 624 END DO 626 625 END DO 627 626 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 627 632 628 IF(ln_ctl) THEN … … 634 630 CALL prt_ctl( tab3d_1=avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 635 631 ENDIF 636 !637 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )638 632 ! 639 633 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 737 731 ! !* set vertical eddy coef. to the background value 738 732 DO jk = 1, jpk 739 avt (:,:,jk) = avtb(jk) * wmask(:,:,jk)740 avm (:,:,jk) = avmb(jk) * wmask(:,:,jk)733 avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 734 avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 741 735 END DO 742 736 dissl(:,:,:) = 1.e-12_wp … … 748 742 749 743 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 ! 744 !!--------------------------------------------------------------------- 745 !! *** ROUTINE tke_rst *** 746 !! 747 !! ** Purpose : Read or write TKE file (en) in restart file 748 !! 749 !! ** Method : use of IOM library 750 !! if the restart does not contain TKE, en is either 751 !! set to rn_emin or recomputed 752 !!---------------------------------------------------------------------- 753 INTEGER , INTENT(in) :: kt ! ocean time-step 754 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 755 ! 756 INTEGER :: jit, jk ! dummy loop indices 757 INTEGER :: id1, id2, id3, id4 ! local integers 758 !!---------------------------------------------------------------------- 759 ! 760 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 761 ! ! --------------- 762 IF( ln_rstart ) THEN !* Read the restart file 763 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 764 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 765 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 766 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 767 ! 768 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 769 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 770 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 771 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 772 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 773 ELSE ! start TKE from rest 774 IF(lwp) WRITE(numout,*) ' ==>> previous run without TKE scheme, set en to background values' 775 en(:,:,:) = rn_emin * wmask(:,:,:) 776 ! avt_k, avm_k already set to the background value in zdf_phy_init 777 ENDIF 778 ELSE !* Start from rest 779 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set en to the background value' 780 en(:,:,:) = rn_emin * wmask(:,:,:) 781 ! avt_k, avm_k already set to the background value in zdf_phy_init 782 ENDIF 783 ! 784 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 785 ! ! ------------------- 786 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 787 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 788 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 789 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 790 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 791 ! 792 ENDIF 793 ! 812 794 END SUBROUTINE tke_rst 813 795
Note: See TracChangeset
for help on using the changeset viewer.
