Changeset 13472 for NEMO/trunk/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-09-16T15:05:19+02:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r13461 r13472 52 52 #endif 53 53 ! 54 #if defined key_si3 55 USE ice, ONLY: hm_i, h_i 56 #endif 57 #if defined key_cice 58 USE sbc_ice, ONLY: h_i 59 #endif 54 60 USE in_out_manager ! I/O manager 55 61 USE iom ! I/O manager library … … 68 74 ! !!** Namelist namzdf_tke ** 69 75 LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not 76 INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 77 REAL(wp) :: rn_mxlice ! ice thickness value when scaling under sea-ice 70 78 INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) 71 79 REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] 72 INTEGER :: nn_mxlice ! type of scaling under sea-ice73 REAL(wp) :: rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1)74 80 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 75 81 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) … … 82 88 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 83 89 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 84 REAL(wp) :: rn_eice ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/485 90 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 86 91 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 92 INTEGER :: nn_eice ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3) 87 93 88 94 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 199 205 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 200 206 ! 201 INTEGER :: ji, jj, jk ! dummy loop arguments207 INTEGER :: ji, jj, jk ! dummy loop arguments 202 208 REAL(wp) :: zetop, zebot, zmsku, zmskv ! local scalars 203 209 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 204 210 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 205 REAL(wp) :: zbbrau, z ri! local scalars206 REAL(wp) :: zfact1, zfact2, zfact3 ! - 207 REAL(wp) :: ztx2 , zty2 , zcof ! - 208 REAL(wp) :: ztau , zdif ! - 209 REAL(wp) :: zus , zwlc , zind ! - 210 REAL(wp) :: zzd_up, zzd_lw ! - 211 REAL(wp) :: zbbrau, zbbirau, zri ! local scalars 212 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 213 REAL(wp) :: ztx2 , zty2 , zcof ! - - 214 REAL(wp) :: ztau , zdif ! - - 215 REAL(wp) :: zus , zwlc , zind ! - - 216 REAL(wp) :: zzd_up, zzd_lw ! - - 211 217 INTEGER , DIMENSION(jpi,jpj) :: imlc 212 REAL(wp), DIMENSION(jpi,jpj) :: z hlc, zfr_i218 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3 213 219 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 214 220 !!-------------------------------------------------------------------- 215 221 ! 216 zbbrau = rn_ebb / rho0 ! Local constant initialisation 217 zfact1 = -.5_wp * rn_Dt 218 zfact2 = 1.5_wp * rn_Dt * rn_ediss 219 zfact3 = 0.5_wp * rn_ediss 222 zbbrau = rn_ebb / rho0 ! Local constant initialisation 223 zbbirau = 3.75_wp / rho0 224 zfact1 = -.5_wp * rn_Dt 225 zfact2 = 1.5_wp * rn_Dt * rn_ediss 226 zfact3 = 0.5_wp * rn_ediss 227 ! 228 ! ice fraction considered for attenuation of langmuir & wave breaking 229 SELECT CASE ( nn_eice ) 230 CASE( 0 ) ; zice_fra(:,:) = 0._wp 231 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(:,:) * 10._wp ) 232 CASE( 2 ) ; zice_fra(:,:) = fr_i(:,:) 233 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 234 END SELECT 220 235 ! 221 236 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 222 237 ! ! Surface/top/bottom boundary condition on tke 223 238 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 224 ! 239 ! 225 240 DO_2D( 0, 0, 0, 0 ) 241 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 242 !! one way around would be to increase zbbirau 243 !! en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 244 !! & fr_i(ji,jj) * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 226 245 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 227 246 END_2D … … 273 292 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 274 293 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 275 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 276 zus = zcof * taum(ji,jj)294 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! Last w-level at which zpelc>=0.5*us*us 295 zus = zcof * taum(ji,jj) ! with us=0.016*wind(starting from jpk-1) 277 296 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 278 297 END_3D … … 284 303 DO_2D( 0, 0, 0, 0 ) 285 304 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 286 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 287 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 305 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 288 306 END_2D 289 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 290 IF ( zfr_i(ji,jj) /= 0. ) THEN 291 ! vertical velocity due to LC 307 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 308 IF ( zus3(ji,jj) /= 0._wp ) THEN 292 309 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 293 310 ! ! vertical velocity due to LC 294 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i311 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) 295 312 ! ! TKE Langmuir circulation source term 296 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * z fr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)313 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 297 314 ENDIF 298 315 ENDIF … … 326 343 ! ! eddy coefficient (ensure numerical stability) 327 344 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 328 & / ( e3t(ji,jj,jk ,Kmm) & 329 & * e3w(ji,jj,jk ,Kmm) ) 345 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) ) 330 346 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 331 & / ( e3t(ji,jj,jk-1,Kmm) & 332 & * e3w(ji,jj,jk ,Kmm) ) 347 & / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) ) 333 348 ! 334 349 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 370 385 371 386 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 372 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 387 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF 373 388 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 374 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) )* wmask(ji,jj,jk) * tmask(ji,jj,1)389 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 375 390 END_3D 376 391 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) … … 378 393 jk = nmln(ji,jj) 379 394 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 380 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) )* wmask(ji,jj,jk) * tmask(ji,jj,1)395 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 381 396 END_2D 382 397 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) … … 388 403 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 389 404 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 390 & * MAX( 0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)405 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 391 406 END_3D 392 407 ENDIF … … 450 465 zmxlm(:,:,:) = rmxl_min 451 466 zmxld(:,:,:) = rmxl_min 452 ! 467 ! 453 468 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 454 469 ! … … 468 483 CASE( 1 ) ! scaling with constant sea-ice thickness 469 484 DO_2D( 0, 0, 0, 0 ) 470 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 485 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 486 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 471 487 END_2D 472 488 ! … … 474 490 DO_2D( 0, 0, 0, 0 ) 475 491 #if defined key_si3 476 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 492 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 493 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 477 494 #elif defined key_cice 478 495 zmaxice = MAXVAL( h_i(ji,jj,:) ) 479 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 496 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 497 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 480 498 #endif 481 499 END_2D … … 484 502 DO_2D( 0, 0, 0, 0 ) 485 503 zmaxice = MAXVAL( h_i(ji,jj,:) ) 486 zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 504 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 505 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 487 506 END_2D 488 507 ! … … 610 629 & rn_mxl0 , nn_mxlice, rn_mxlice, & 611 630 & nn_pdl , ln_lc , rn_lc , & 612 & nn_etau , nn_htau , rn_efr , rn_eice631 & nn_etau , nn_htau , rn_efr , nn_eice 613 632 !!---------------------------------------------------------------------- 614 633 ! … … 642 661 ENDIF 643 662 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 663 IF( ln_mxl0 ) THEN 664 WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice 665 IF( nn_mxlice == 1 ) & 666 WRITE(numout,*) ' ice thickness when scaling under sea-ice rn_mxlice = ', rn_mxlice 667 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 668 CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice' 669 CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness' 670 CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness' 671 CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness' 672 CASE DEFAULT 673 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 674 END SELECT 675 ENDIF 644 676 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 645 677 WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc … … 647 679 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau 648 680 WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr 649 WRITE(numout,*) ' below sea-ice: =0 ON rn_eice = ', rn_eice 650 WRITE(numout,*) ' =4 OFF when ice fraction > 1/4 ' 681 WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice 682 SELECT CASE( nn_eice ) 683 CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on langmuir & surface wave breaking' 684 CASE( 1 ) ; WRITE(numout,*) ' ==>>> weigthed by 1-TANH( fr_i(:,:) * 10 )' 685 CASE( 2 ) ; WRITE(numout,*) ' ==>>> weighted by 1-fr_i(:,:)' 686 CASE( 3 ) ; WRITE(numout,*) ' ==>>> weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 687 CASE DEFAULT 688 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 689 END SELECT 651 690 IF( .NOT.ln_drg_OFF ) THEN 652 691 WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.
