- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdftke.F90
r13466 r13469 231 231 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 232 232 ! 233 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 234 DO ji = fs_2, fs_jpim1 ! vector opt. 233 DO_2D_00_00 235 234 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 236 235 !! one way around would be to increase zbbirau 237 236 !! en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 238 237 !! & fr_i(ji,jj) * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 239 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 240 END DO 241 END DO 238 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 239 END_2D 242 240 ! 243 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 251 249 IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE 252 250 ! 253 DO jj = 2, jpjm1 ! bottom friction 254 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 256 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 257 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 258 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 259 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 260 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 261 END DO 262 END DO 251 DO_2D_00_00 252 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 253 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 254 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 255 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Nnn)+uu(ji-1,jj,mbkt(ji,jj),Nnn) ) )**2 & 256 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Nnn)+vv(ji,jj-1,mbkt(ji,jj),Nnn) ) )**2 ) 257 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 258 END_2D 263 259 IF( ln_isfcav ) THEN ! top friction 264 DO jj = 2, jpjm1 265 DO ji = fs_2, fs_jpim1 ! vector opt. 266 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 267 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 268 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 269 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 270 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 271 en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) & 272 & + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 273 END DO 274 END DO 260 DO_2D_00_00 261 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 262 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 263 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 264 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Nnn)+uu(ji-1,jj,mikt(ji,jj),Nnn) ) )**2 & 265 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Nnn)+vv(ji,jj-1,mikt(ji,jj),Nnn) ) )**2 ) 266 en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) & 267 & + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 268 END_2D 275 269 ENDIF 276 270 ! … … 289 283 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 290 284 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 291 DO jk = jpkm1, 2, -1 292 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 293 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 294 zus = zcof * taum(ji,jj) 295 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 296 END DO 297 END DO 298 END DO 285 DO_3D_11_11( jpkm1, 2, -1 ) 286 zus = zcof * taum(ji,jj) 287 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 288 END_3D 299 289 ! ! finite LC depth 300 DO jj = 1, jpj 301 DO ji = 1, jpi 302 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 303 END DO 304 END DO 290 DO_2D_11_11 291 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 292 END_2D 305 293 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 306 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 ! vector opt. 308 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 309 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 310 END DO 311 END DO 312 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 313 DO jj = 2, jpjm1 314 DO ji = fs_2, fs_jpim1 ! vector opt. 315 IF ( zus3(ji,jj) /= 0._wp ) THEN 316 ! vertical velocity due to LC 317 IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 318 ! ! vertical velocity due to LC 319 zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 320 ! ! TKE Langmuir circulation source term 321 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 322 ENDIF 323 ENDIF 324 END DO 325 END DO 326 END DO 294 DO_2D_00_00 295 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 296 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 297 END_2D 298 DO_3D_00_00( 2, jpkm1 ) 299 IF ( zus3(ji,jj) /= 0._wp ) THEN 300 ! vertical velocity due to LC 301 IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 302 ! ! vertical velocity due to LC 303 zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 304 ! ! TKE Langmuir circulation source term 305 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 306 ENDIF 307 ENDIF 308 END_3D 327 309 ! 328 310 ENDIF … … 336 318 ! 337 319 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 338 DO jk = 2, jpkm1 339 DO jj = 2, jpjm1 340 DO ji = 2, jpim1 341 ! ! local Richardson number 342 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 343 ! ! inverse of Prandtl number 344 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 345 END DO 346 END DO 347 END DO 320 DO_3D_00_00( 2, jpkm1 ) 321 ! ! local Richardson number 322 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 323 ! ! inverse of Prandtl number 324 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 325 END_3D 348 326 ENDIF 349 327 ! 350 DO jk = 2, jpkm1 !* Matrix and right hand side in en 351 DO jj = 2, jpjm1 352 DO ji = fs_2, fs_jpim1 ! vector opt. 353 zcof = zfact1 * tmask(ji,jj,jk) 354 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 355 ! ! eddy coefficient (ensure numerical stability) 356 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 357 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 358 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 359 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 360 ! 361 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 362 zd_lw(ji,jj,jk) = zzd_lw 363 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 364 ! 365 ! ! right hand side in en 366 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 367 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 368 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 369 & ) * wmask(ji,jj,jk) 370 END DO 371 END DO 372 END DO 328 DO_3D_00_00( 2, jpkm1 ) 329 zcof = zfact1 * tmask(ji,jj,jk) 330 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 331 ! ! eddy coefficient (ensure numerical stability) 332 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 333 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 334 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 335 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 336 ! 337 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 338 zd_lw(ji,jj,jk) = zzd_lw 339 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 340 ! 341 ! ! right hand side in en 342 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 343 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 344 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 345 & ) * wmask(ji,jj,jk) 346 END_3D 373 347 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 374 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 375 DO jj = 2, jpjm1 376 DO ji = fs_2, fs_jpim1 ! vector opt. 377 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 378 END DO 379 END DO 380 END DO 381 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 384 END DO 385 END DO 386 DO jk = 3, jpkm1 387 DO jj = 2, jpjm1 388 DO ji = fs_2, fs_jpim1 ! vector opt. 389 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) 390 END DO 391 END DO 392 END DO 393 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 394 DO ji = fs_2, fs_jpim1 ! vector opt. 395 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 396 END DO 397 END DO 398 DO jk = jpk-2, 2, -1 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 401 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 402 END DO 403 END DO 404 END DO 405 DO jk = 2, jpkm1 ! set the minimum value of tke 406 DO jj = 2, jpjm1 407 DO ji = fs_2, fs_jpim1 ! vector opt. 408 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 409 END DO 410 END DO 411 END DO 348 DO_3D_00_00( 3, jpkm1 ) 349 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 350 END_3D 351 DO_2D_00_00 352 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 353 END_2D 354 DO_3D_00_00( 3, jpkm1 ) 355 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) 356 END_3D 357 DO_2D_00_00 358 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 359 END_2D 360 DO_3D_00_00( jpk-2, 2, -1 ) 361 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 362 END_3D 363 DO_3D_00_00( 2, jpkm1 ) 364 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 365 END_3D 412 366 ! 413 367 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 419 373 420 374 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 421 DO jk = 2, jpkm1 ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF 422 DO jj = 2, jpjm1 423 DO ji = fs_2, fs_jpim1 ! vector opt. 424 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 425 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 426 END DO 427 END DO 428 END DO 375 DO_3D_00_00( 2, jpkm1 ) 376 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 377 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 378 END_3D 429 379 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 430 DO jj = 2, jpjm1 431 DO ji = fs_2, fs_jpim1 ! vector opt. 432 jk = nmln(ji,jj) 433 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 434 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 435 END DO 436 END DO 380 DO_2D_00_00 381 jk = nmln(ji,jj) 382 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 383 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 384 END_2D 437 385 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 438 DO jk = 2, jpkm1 439 DO jj = 2, jpjm1 440 DO ji = fs_2, fs_jpim1 ! vector opt. 441 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 442 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 443 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 444 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 445 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 446 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 447 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 448 END DO 449 END DO 450 END DO 386 DO_3D_00_00( 2, jpkm1 ) 387 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 388 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 389 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 390 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 391 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 392 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 393 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 394 END_3D 451 395 ENDIF 452 396 ! … … 515 459 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 516 460 #if ! defined key_si3 && ! defined key_cice 517 DO jj = 2, jpjm1 ! No sea-ice 518 DO ji = fs_2, fs_jpim1 519 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 520 END DO 521 END DO 461 DO_2D_00_00 462 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 463 END_2D 522 464 #else 523 465 … … 525 467 ! 526 468 CASE( 0 ) ! No scaling under sea-ice 527 DO jj = 2, jpjm1 528 DO ji = fs_2, fs_jpim1 529 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 530 END DO 531 END DO 469 DO_2D_00_00 470 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 471 END_2D 532 472 ! 533 473 CASE( 1 ) ! scaling with constant sea-ice thickness 534 DO jj = 2, jpjm1 535 DO ji = fs_2, fs_jpim1 536 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 537 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 538 END DO 539 END DO 474 DO_2D_00_00 475 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 476 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 477 END_2D 540 478 ! 541 479 CASE( 2 ) ! scaling with mean sea-ice thickness 542 DO jj = 2, jpjm1 543 DO ji = fs_2, fs_jpim1 480 DO_2D_00_00 544 481 #if defined key_si3 545 546 482 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 483 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 547 484 #elif defined key_cice 548 549 550 485 zmaxice = MAXVAL( h_i(ji,jj,:) ) 486 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 487 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 551 488 #endif 552 END DO 553 END DO 489 END_2D 554 490 ! 555 491 CASE( 3 ) ! scaling with max sea-ice thickness 556 DO jj = 2, jpjm1 557 DO ji = fs_2, fs_jpim1 558 zmaxice = MAXVAL( h_i(ji,jj,:) ) 559 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 560 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 561 END DO 562 END DO 492 DO_2D_00_00 493 zmaxice = MAXVAL( h_i(ji,jj,:) ) 494 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 495 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 496 END_2D 563 497 ! 564 498 END SELECT 565 499 #endif 566 500 ! 567 DO jj = 2, jpjm1 568 DO ji = fs_2, fs_jpim1 569 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 570 END DO 571 END DO 501 DO_2D_00_00 502 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 503 END_2D 572 504 ! 573 505 ELSE … … 575 507 ENDIF 576 508 ! 577 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 578 DO jj = 2, jpjm1 579 DO ji = fs_2, fs_jpim1 ! vector opt. 580 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 581 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 582 END DO 583 END DO 584 END DO 509 DO_3D_00_00( 2, jpkm1 ) 510 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 511 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 512 END_3D 585 513 ! 586 514 ! !* Physical limits for the mixing length … … 594 522 ! where wmask = 0 set zmxlm == p_e3w 595 523 CASE ( 0 ) ! bounded by the distance to surface and bottom 596 DO jk = 2, jpkm1 597 DO jj = 2, jpjm1 598 DO ji = fs_2, fs_jpim1 ! vector opt. 599 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 600 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 601 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 602 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 603 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 604 END DO 605 END DO 606 END DO 524 DO_3D_00_00( 2, jpkm1 ) 525 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 526 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 527 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 528 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 529 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 530 END_3D 607 531 ! 608 532 CASE ( 1 ) ! bounded by the vertical scale factor 609 DO jk = 2, jpkm1 610 DO jj = 2, jpjm1 611 DO ji = fs_2, fs_jpim1 ! vector opt. 612 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 613 zmxlm(ji,jj,jk) = zemxl 614 zmxld(ji,jj,jk) = zemxl 615 END DO 616 END DO 617 END DO 533 DO_3D_00_00( 2, jpkm1 ) 534 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 535 zmxlm(ji,jj,jk) = zemxl 536 zmxld(ji,jj,jk) = zemxl 537 END_3D 618 538 ! 619 539 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 620 DO jk = 2, jpkm1 ! from the surface to the bottom : 621 DO jj = 2, jpjm1 622 DO ji = fs_2, fs_jpim1 ! vector opt. 623 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 624 END DO 625 END DO 626 END DO 627 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 628 DO jj = 2, jpjm1 629 DO ji = fs_2, fs_jpim1 ! vector opt. 630 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 631 zmxlm(ji,jj,jk) = zemxl 632 zmxld(ji,jj,jk) = zemxl 633 END DO 634 END DO 635 END DO 540 DO_3D_00_00( 2, jpkm1 ) 541 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 542 END_3D 543 DO_3D_00_00( jpkm1, 2, -1 ) 544 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 545 zmxlm(ji,jj,jk) = zemxl 546 zmxld(ji,jj,jk) = zemxl 547 END_3D 636 548 ! 637 549 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 638 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 639 DO jj = 2, jpjm1 640 DO ji = fs_2, fs_jpim1 ! vector opt. 641 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 642 END DO 643 END DO 644 END DO 645 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 646 DO jj = 2, jpjm1 647 DO ji = fs_2, fs_jpim1 ! vector opt. 648 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 649 END DO 650 END DO 651 END DO 652 DO jk = 2, jpkm1 653 DO jj = 2, jpjm1 654 DO ji = fs_2, fs_jpim1 ! vector opt. 655 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 656 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 657 zmxlm(ji,jj,jk) = zemlm 658 zmxld(ji,jj,jk) = zemlp 659 END DO 660 END DO 661 END DO 550 DO_3D_00_00( 2, jpkm1 ) 551 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 552 END_3D 553 DO_3D_00_00( jpkm1, 2, -1 ) 554 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 555 END_3D 556 DO_3D_00_00( 2, jpkm1 ) 557 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 558 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 559 zmxlm(ji,jj,jk) = zemlm 560 zmxld(ji,jj,jk) = zemlp 561 END_3D 662 562 ! 663 563 END SELECT … … 666 566 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 667 567 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 668 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 669 DO jj = 2, jpjm1 670 DO ji = fs_2, fs_jpim1 ! vector opt. 671 zsqen = SQRT( en(ji,jj,jk) ) 672 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 673 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 674 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 675 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 676 END DO 677 END DO 678 END DO 568 DO_3D_00_00( 1, jpkm1 ) 569 zsqen = SQRT( en(ji,jj,jk) ) 570 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 571 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 572 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 573 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 574 END_3D 679 575 ! 680 576 ! 681 577 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 682 DO jk = 2, jpkm1 683 DO jj = 2, jpjm1 684 DO ji = fs_2, fs_jpim1 ! vector opt. 685 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 686 END DO 687 END DO 688 END DO 578 DO_3D_00_00( 2, jpkm1 ) 579 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 580 END_3D 689 581 ENDIF 690 582 !
Note: See TracChangeset
for help on using the changeset viewer.
