Changeset 12489 for NEMO/trunk/src/OCE/ZDF
- Timestamp:
- 2020-02-28T16:55:11+01:00 (4 years ago)
- Location:
- NEMO/trunk/src/OCE/ZDF
- Files:
-
- 7 edited
-
zdfdrg.F90 (modified) (1 diff)
-
zdfgls.F90 (modified) (4 diffs)
-
zdfiwm.F90 (modified) (10 diffs)
-
zdfmxl.F90 (modified) (1 diff)
-
zdfosm.F90 (modified) (12 diffs)
-
zdfric.F90 (modified) (1 diff)
-
zdftke.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdfdrg.F90
r12377 r12489 165 165 !!--------------------------------------------------------------------- 166 166 ! 167 !!gm bug : time step is only r dt (not 2 rdt if euler start !)168 zm1_2dt = - 1._wp / ( 2._wp * r dt )167 !!gm bug : time step is only rn_Dt (not 2 rn_Dt if euler start !) 168 zm1_2dt = - 1._wp / ( 2._wp * rn_Dt ) 169 169 170 170 IF( l_trddyn ) THEN ! trends: store the input trends -
NEMO/trunk/src/OCE/ZDF/zdfgls.F90
r12377 r12489 170 170 ! 171 171 ! surface friction 172 ustar2_surf(ji,jj) = r1_r au0 * taum(ji,jj) * tmask(ji,jj,1)172 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) 173 173 ! 174 174 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... … … 267 267 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 268 268 ! ! diagonal 269 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + r dt * zdiss * wmask(ji,jj,jk)269 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 270 270 ! ! right hand side in en 271 en(ji,jj,jk) = en(ji,jj,jk) + r dt * zesh2 * wmask(ji,jj,jk)271 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 272 272 END_3D 273 273 ! … … 477 477 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 478 478 ! ! diagonal 479 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + r dt * zdiss * wmask(ji,jj,jk)479 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 480 480 ! ! right hand side in psi 481 psi(ji,jj,jk) = psi(ji,jj,jk) + r dt * zesh2 * wmask(ji,jj,jk)481 psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 482 482 END_3D 483 483 ! … … 1007 1007 rc04 = rc03 * rc0 1008 1008 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1009 rsbc_tke2 = r dt * rn_crban / rl_sf ! Neumann + Wave breaking1009 rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking 1010 1010 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1011 1011 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1012 1012 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1013 1013 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1014 rsbc_psi1 = -0.5_wp * r dt * rc0**(rpp-2._wp*rmm) / rsc_psi1015 rsbc_psi2 = -0.5_wp * r dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking1016 ! 1017 rfact_tke = -0.5_wp / rsc_tke * r dt ! Cst used for the Diffusion term of tke1018 rfact_psi = -0.5_wp / rsc_psi * r dt ! Cst used for the Diffusion term of tke1014 rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi 1015 rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1016 ! 1017 rfact_tke = -0.5_wp / rsc_tke * rn_Dt ! Cst used for the Diffusion term of tke 1018 rfact_psi = -0.5_wp / rsc_psi * rn_Dt ! Cst used for the Diffusion term of tke 1019 1019 ! 1020 1020 ! !* Wall proximity function -
NEMO/trunk/src/OCE/ZDF/zdfiwm.F90
r12377 r12489 87 87 !! This is divided into three components: 88 88 !! 1. Bottom-intensified low-mode dissipation at critical slopes 89 !! zemx_iwm(z) = ( ecri_iwm / r au0 ) * EXP( -(H-z)/hcri_iwm )89 !! zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm ) 90 90 !! / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 91 91 !! where hcri_iwm is the characteristic length scale of the bottom 92 92 !! intensification, ecri_iwm a map of available power, and H the ocean depth. 93 93 !! 2. Pycnocline-intensified low-mode dissipation 94 !! zemx_iwm(z) = ( epyc_iwm / r au0 ) * ( sqrt(rn2(z))^nn_zpyc )94 !! zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 95 95 !! / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 96 96 !! where epyc_iwm is a map of available power, and nn_zpyc … … 98 98 !! energy dissipation. 99 99 !! 3. WKB-height dependent high mode dissipation 100 !! zemx_iwm(z) = ( ebot_iwm / r au0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)100 !! zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 101 101 !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 102 102 !! where hbot_iwm is the characteristic length scale of the WKB bottom … … 151 151 DO_2D_11_11 152 152 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 153 zfact(ji,jj) = r au0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) )153 zfact(ji,jj) = rho0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) 154 154 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 155 155 END_2D … … 181 181 ! 182 182 DO_2D_11_11 183 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( r au0 * zfact(ji,jj) )183 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 184 184 END_2D 185 185 ! … … 196 196 ! 197 197 DO_2D_11_11 198 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( r au0 * zfact(ji,jj) )198 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 199 199 END_2D 200 200 ! … … 243 243 ! 244 244 DO_2D_11_11 245 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( r au0 * zfact(ji,jj) )245 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 246 246 END_2D 247 247 ! … … 255 255 ! Calculate molecular kinematic viscosity 256 256 znu_t(:,:,:) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm) & 257 & + 0.02305_wp * ts(:,:,:,jp_sal,Kmm) ) * tmask(:,:,:) * r1_r au0257 & + 0.02305_wp * ts(:,:,:,jp_sal,Kmm) ) * tmask(:,:,:) * r1_rho0 258 258 DO jk = 2, jpkm1 259 259 znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) … … 293 293 END_3D 294 294 CALL mpp_sum( 'zdfiwm', zztmp ) 295 zztmp = r au0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing295 zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 296 296 ! 297 297 IF(lwp) THEN … … 337 337 !* output useful diagnostics: Kz*N^2 , 338 338 !!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 339 ! vertical integral of r au0 * Kz * N^2 , energy density (zemx_iwm)339 ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 340 340 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 341 341 ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) … … 345 345 z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 346 346 END DO 347 z2d(:,:) = r au0 * z2d(:,:)347 z2d(:,:) = rho0 * z2d(:,:) 348 348 CALL iom_put( "bflx_iwm", z3d ) 349 349 CALL iom_put( "pcmap_iwm", z2d ) -
NEMO/trunk/src/OCE/ZDF/zdfmxl.F90
r12377 r12489 97 97 nmln(:,:) = nlb10 ! Initialization to the number of w ocean point 98 98 hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 99 zN2_c = grav * rho_c * r1_r au0 ! convert density criteria into N^2 criteria99 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 100 100 DO_3D_11_11( nlb10, jpkm1 ) 101 101 ikt = mbkt(ji,jj) -
NEMO/trunk/src/OCE/ZDF/zdfosm.F90
r12377 r12489 300 300 DO_2D_00_00 301 301 ! Surface downward irradiance (so always +ve) 302 zrad0(ji,jj) = qsr(ji,jj) * r1_r au0_rcp302 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 303 303 ! Downwards irradiance at base of boundary layer 304 304 zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) … … 312 312 zbeta = rab_n(ji,jj,1,jp_sal) 313 313 ! Upwards surface Temperature flux for non-local term 314 zwth0(ji,jj) = - qns(ji,jj) * r1_r au0_rcp * tmask(ji,jj,1)314 zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 315 315 ! Upwards surface salinity flux for non-local term 316 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_r au0 * tmask(ji,jj,1)316 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 317 317 ! Non radiative upwards surface buoyancy flux 318 318 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) … … 324 324 zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj) 325 325 ! Surface upward velocity fluxes 326 zuw0(ji,jj) = -utau(ji,jj) * r1_r au0 * tmask(ji,jj,1)327 zvw0(ji,jj) = -vtau(ji,jj) * r1_r au0 * tmask(ji,jj,1)326 zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 327 zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 328 328 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 329 329 zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) … … 441 441 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 442 442 443 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_ rdt / hbl(ji,jj) ) &443 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 444 444 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 445 445 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. … … 447 447 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 448 448 449 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_ rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / &449 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 450 450 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 451 451 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) … … 458 458 IF ( zzdhdt < 0._wp ) THEN 459 459 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 460 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_ rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)460 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 461 461 ELSE 462 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_ rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) &462 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 463 463 & + MAX( zdb_bl(ji,jj), 0.0 ) 464 464 ENDIF … … 472 472 ibld(:,:) = 3 473 473 474 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_ rdt ! certainly need wb here, so subtract it474 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 475 475 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 476 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_ rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom476 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 477 477 478 478 DO_3D_00_00( 4, jpkm1 ) … … 496 496 IF ( lconv(ji,jj) ) THEN 497 497 !unstable 498 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_ rdt / hbl(ji,jj) ) &498 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 499 499 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 500 500 … … 503 503 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 504 504 505 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_ rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) )505 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 506 506 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 507 507 … … 1250 1250 IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind ) ! x surface Stokes drift 1251 1251 IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind ) ! y surface Stokes drift 1252 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1252 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1253 1253 ! Stokes drift read in from sbcwave (=2). 1254 1254 CASE(2) 1255 1255 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1256 1256 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1257 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*r au0*tmask(:,:,1)*zustar**2* &1257 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 1258 1258 & SQRT(ut0sd**2 + vt0sd**2 ) ) 1259 1259 END SELECT … … 1271 1271 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1272 1272 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1273 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*r au0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine1274 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*r au0*tmask(:,:,1)*zustar**2*zustke )1273 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1274 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 1275 1275 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1276 1276 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine … … 1500 1500 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 1501 1501 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 1502 zN2_c = grav * rho_c * r1_r au0 ! convert density criteria into N^2 criteria1502 zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 1503 1503 ! 1504 1504 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 -
NEMO/trunk/src/OCE/ZDF/zdfric.F90
r12377 r12489 174 174 ! 175 175 DO_2D_00_00 176 zustar = SQRT( taum(ji,jj) * r1_r au0 )176 zustar = SQRT( taum(ji,jj) * r1_rho0 ) 177 177 zhek = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 178 178 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zhek , rn_mldmax ) ) ! set allowed range -
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r12377 r12489 206 206 !!-------------------------------------------------------------------- 207 207 ! 208 zbbrau = rn_ebb / r au0 ! Local constant initialisation209 zfact1 = -.5_wp * r dt210 zfact2 = 1.5_wp * r dt * rn_ediss208 zbbrau = rn_ebb / rho0 ! Local constant initialisation 209 zfact1 = -.5_wp * rn_Dt 210 zfact2 = 1.5_wp * rn_Dt * rn_ediss 211 211 zfact3 = 0.5_wp * rn_ediss 212 212 ! … … 228 228 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 229 229 ! 230 ! en(bot) = (ebb0/r au0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)230 ! en(bot) = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 231 231 ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 232 232 ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 … … 237 237 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 238 238 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 239 ! ! where 0.001875 = (rn_ebb0/r au0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)239 ! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 240 240 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 241 241 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) … … 246 246 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 247 247 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 248 ! ! where 0.001875 = (rn_ebb0/r au0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)248 ! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 249 249 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 250 250 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) … … 288 288 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i 289 289 ! ! TKE Langmuir circulation source term 290 en(ji,jj,jk) = en(ji,jj,jk) + r dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)290 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 291 291 ENDIF 292 292 ENDIF … … 325 325 ! 326 326 ! ! right hand side in en 327 en(ji,jj,jk) = en(ji,jj,jk) + r dt * ( p_sh2(ji,jj,jk)& ! shear327 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * ( p_sh2(ji,jj,jk) & ! shear 328 328 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 329 329 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation … … 439 439 zmxld(:,:,:) = rmxl_min 440 440 ! 441 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(r au0*g)442 zraug = vkarmn * 2.e5_wp / ( r au0 * grav )441 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 442 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 443 443 DO_2D_00_00 444 444 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
Note: See TracChangeset
for help on using the changeset viewer.
