Changeset 14413
- Timestamp:
- 2021-02-05T17:18:45+01:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715
- Files:
-
- 5 edited
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (5 diffs)
-
cfgs/SHARED/namelist_ref (modified) (2 diffs)
-
src/OCE/DIA/diawri.F90 (modified) (1 diff)
-
src/OCE/TRA/tranxt.F90 (modified) (1 diff)
-
src/OCE/ZDF/zdfosm.F90 (modified) (57 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/field_def_nemo-oce.xml
r13403 r14413 202 202 <field id="ibld" long_name="index of boundary layer depth" unit="#" /> 203 203 <field id="imld" long_name="index of mixed layer depth" unit="#" /> 204 <field id="jp_ext" long_name="flag =1 if pycnocline well resolved" unit="#" /> 205 <field id="j_ddh" long_name="index of mixed layer depth" unit="#" /> 206 <field id="zshear" long_name="shear production of TKE " unit="m^3/s^3" /> 204 207 <field id="zhbl" long_name="boundary layer depth -grid" unit="m" /> 205 208 <field id="zhml" long_name="mixed layer depth - grid" unit="m" /> … … 211 214 <field id="zwth0" long_name="surface non-local temperature flux" unit="deg m/s" /> 212 215 <field id="zws0" long_name="surface non-local salinity flux" unit="psu m/s" /> 216 <field id="zwb0" long_name="surface non-local buoyancy flux" unit="m^2/s^3" /> 213 217 <field id="zwstrc" long_name="convective velocity scale" unit="m/s" /> 214 218 <field id="zustar" long_name="friction velocity" unit="m/s" /> … … 221 225 222 226 <!-- interior BL OSMOSIS diagnostics --> 223 <field id="zw thav" long_name="av turb flux of T in ml" unit="deg m/s" />227 <field id="zwbav" long_name="av turb flux of buoyancy in ml" unit="m^2/s^3" /> 224 228 <field id="zt_ml" long_name="av T in ml" unit="deg" /> 225 229 <field id="zhol" long_name="Hoenekker number" unit="#" /> … … 233 237 <field id="zdu_bl" long_name="u jump at base of BL" unit="m/s" /> 234 238 <field id="zdv_bl" long_name="v jump at base of BL" unit="m/s" /> 235 239 <field id="zdt_ml" long_name="temperature jump at base of ML" unit="deg" /> 240 <field id="zds_ml" long_name="salinity jump at base of ML" unit="10^-3" /> 241 <field id="zdb_ml" long_name="buoyancy jump at base of ML" unit="m/s^2" /> 236 242 <!-- extra OSMOSIS diagnostics for debugging --> 237 243 <field id="zsc_uw_1_0" long_name="zsc u-momentum flux on T after Stokes" unit="m^2/s^2" /> … … 240 246 <field id="zsc_uw_2_f" long_name="2nd zsc u-momentum flux on T after Coriolis" unit="m^2/s^2" /> 241 247 <field id="zsc_vw_2_f" long_name="2nd zsc v-momentum flux on T after Coriolis" unit="m^2/s^2" /> 242 <field id="zuw_bse" long_name="base u-flux T-points" unit="m^2/s^2" />243 <field id="zvw_bse" long_name="base v-flux T-points" unit="m^2/s^2" />244 248 245 249 <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/namelist_ref
r13478 r14413 1093 1093 rn_zdfosm_adjust_sd = 1.0 ! Stokes drift reduction factor 1094 1094 rn_osm_hblfrac = 0.1 ! specify top part of hbl for nn_osm_wave = 3 or 4 1095 rn_osm_bl_thresh = 5.e-5 !Threshold buoyancy for deepening of OSBL base 1095 1096 nn_ave = 0 ! choice of horizontal averaging on avt, avmu, avmv 1096 1097 ln_dia_osm = .true. ! output OSMOSIS-OBL variables … … 1119 1120 nn_osm_mle = 0 ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 1120 1121 rn_osm_mle_lf = 5.e+3 ! typical scale of mixed layer front (meters) (case rn_osm_mle=0) 1121 rn_osm_mle_time = 172800. ! time scale for mixing momentum across the mixed layer (seconds) (case rn_osm_mle=0)1122 rn_osm_mle_time = 43200. ! time scale for mixing momentum across the mixed layer (seconds) (case rn_osm_mle=0) 1122 1123 rn_osm_mle_lat = 20. ! reference latitude (degrees) of MLE coef. (case rn_mle=1) 1123 rn_osm_mle_rho_c = 0.01! delta rho criterion used to calculate MLD for FK1124 rn_osm_mle_thresh = 0.0005! delta b criterion used for FK MLD criterion1125 rn_osm_mle_tau = 172800.! time scale for FK-OSM (seconds) (case rn_osm_mle=0)1126 ln_osm_hmle_limit = . false. !limit hmle to rn_osm_hmle_limit*hbl1127 rn_osm_hmle_limit = 1. 21128 1124 rn_osm_mle_rho_c = 0.03 ! delta rho criterion used to calculate MLD for FK 1125 rn_osm_mle_thresh = 0.0001 ! delta b criterion used for FK MLD criterion 1126 rn_osm_mle_tau = 172800. ! time scale for FK-OSM (seconds) (case rn_osm_mle=0) 1127 ln_osm_hmle_limit = .true. ! If true, limit hmle to rn_osm_hmle_limit*hbl 1128 rn_osm_hmle_limit = 1.5 1129 / 1129 1130 !----------------------------------------------------------------------- 1130 1131 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/DIA/diawri.F90
r13402 r14413 935 935 CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing 936 936 CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing 937 CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu* wmask ) ! non-local u forcing938 CALL iom_rstput( 0, 0, inum, 'ghamv', gham u*wmask ) ! non-local v forcing937 CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask ) ! non-local u forcing 938 CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask ) ! non-local v forcing 939 939 IF( ln_osm_mle ) THEN 940 940 CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/TRA/tranxt.F90
r11715 r14413 119 119 IF( l_trdtra ) THEN 120 120 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 121 ztrdt(:,:, jpk) = 0._wp122 ztrds(:,:, jpk) = 0._wp121 ztrdt(:,:,:) = 0._wp 122 ztrds(:,:,:) = 0._wp 123 123 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 124 124 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt ) -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90
r13685 r14413 58 58 USE traqsr ! details of solar radiation absorption 59 59 USE zdfddm ! double diffusion mixing (avs array) 60 ! USE zdfmxl ! mixed layer depth 60 61 USE iom ! I/O library 61 62 USE lib_mpp ! MPP library … … 133 134 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 134 135 REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. 136 REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base. 135 137 REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. 136 138 … … 245 247 REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux 246 248 REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux 249 REAL(wp), DIMENSION(jpi,jpj) :: zwb0tot ! Total surface buoyancy flux including insolation 247 250 REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average 248 251 REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average 249 252 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 250 253 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 254 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 251 255 252 256 … … 262 266 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 263 267 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 268 LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers 269 LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present 270 LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer. 271 LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness. 264 272 265 273 ! mixed-layer variables … … 267 275 INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 268 276 INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 277 INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level 278 INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 269 279 270 280 REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients … … 279 289 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 280 290 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 281 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2 ! correction to dhdt due to internal structure. 282 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext ! external temperature/salinity and buoyancy gradients 283 291 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients 292 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients 284 293 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. 285 294 286 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 287 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer 288 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 289 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 290 REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 291 REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline 295 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer 296 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer 297 REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer 298 REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer 299 REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer 300 REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer 301 ! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 302 REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 303 REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline 292 304 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline 293 305 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline … … 295 307 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline 296 308 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline 297 309 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 298 310 ! Flux-gradient relationship variables 311 REAL(wp), DIMENSION(jpi, jpj) :: zshear ! Shear production. 299 312 300 313 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale. 301 314 302 REAL(wp) , DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity315 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline. 303 316 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms. 317 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/ 304 318 REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms. 305 319 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep … … 316 330 REAL(wp) :: zthick, zz0, zz1 ! temporary variables 317 331 REAL(wp) :: zvel_max, zhbl_s ! temporary variables 318 REAL(wp) :: zfac 332 REAL(wp) :: zfac, ztmp ! temporary variable 319 333 REAL(wp) :: zus_x, zus_y ! temporary Stokes drift 320 334 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity 321 335 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 322 336 REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 337 REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML. 338 REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc 339 REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max 323 340 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 324 REAL(wp) :: z wb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt341 REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau 325 342 REAL(wp) :: zzeta_s = 0._wp 326 343 REAL(wp) :: zzeta_v = 0.46 … … 339 356 zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp 340 357 zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 341 zhol(:,:) = 0._wp 342 lconv(:,:) = .FALSE. 358 zhol(:,:) = 0._wp ; zwb0tot(:,:) = 0._wp 359 lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. 343 360 ! mixed layer 344 361 ! no initialization of zhbl or zhml (or zdh?) 345 362 zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp 346 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 347 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 348 349 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 350 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp 363 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp 364 zv_bl(:,:) = 0._wp ; zb_bl(:,:) = 0._wp 365 zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 366 zt_mle(:,:) = 0._wp ; zs_mle(:,:) = 0._wp ; zu_mle(:,:) = 0._wp 367 zb_mle(:,:) = 0._wp 368 zv_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 369 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp 351 370 zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp 352 zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:) = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp 353 zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp 371 zdb_ml(:,:) = 0._wp 372 zdt_mle(:,:) = 0._wp ; zds_mle(:,:) = 0._wp ; zdu_mle(:,:) = 0._wp 373 zdv_mle(:,:) = 0._wp ; zdb_mle(:,:) = 0._wp 374 zwth_ent = 0._wp ; zws_ent = 0._wp 354 375 ! 355 376 zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp 356 377 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 357 378 ! 358 zdtdz_ ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp379 zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp 359 380 360 381 IF ( ln_osm_mle ) THEN ! only initialise arrays if needed … … 367 388 368 389 ! Flux-Gradient arrays. 369 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp370 zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp371 390 zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp 372 391 zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp … … 376 395 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 377 396 378 zdhdt_2(:,:) = 0._wp 379 ! hbl = MAX(hbl,epsln) 397 ! hbl = MAX(hbl,epsln) 380 398 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 381 399 ! Calculate boundary layer scales … … 407 425 ! Non radiative upwards surface buoyancy flux 408 426 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) 427 ! Total upwards surface buoyancy flux 428 zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * zrad0(ji,jj) 409 429 ! turbulent heat flux averaged over depth of OSBL 410 430 zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) … … 548 568 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 549 569 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 550 lconv(ji,jj) = .TRUE. 551 ELSE 570 ELSE 552 571 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 553 lconv(ji,jj) = .FALSE.554 572 ENDIF 555 573 END DO … … 584 602 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 585 603 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 604 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 586 605 END DO 587 606 END DO 588 607 ! Averages over well-mixed and boundary layer 589 ! Alan: do we need zb_nl?, zb_ml? 590 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 591 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 592 ! External gradient 593 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 594 608 jp_ext(:,:) = 2 609 CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 610 ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 611 CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 612 ! Velocity components in frame aligned with surface stress. 613 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 614 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 615 ! Determine the state of the OSBL, stable/unstable, shear/no shear 616 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 595 617 596 618 IF ( ln_osm_mle ) THEN 597 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 598 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 599 ENDIF 600 619 ! Fox-Kemper Scheme 620 mld_prof = 4 621 DO jk = 5, jpkm1 622 DO jj = 2, jpjm1 623 DO ji = 2, jpim1 624 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 625 END DO 626 END DO 627 END DO 628 jp_ext_mle(:,:) = 2 629 CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) 630 631 DO jj = 2, jpjm1 632 DO ji = 2, jpim1 633 zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 634 END DO 635 END DO 636 637 !! External gradient 638 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 639 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 640 CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 641 CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 642 CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 643 ELSE ! ln_osm_mle 644 ! FK not selected, Boundary Layer only. 645 lpyc(:,:) = .TRUE. 646 lflux(:,:) = .FALSE. 647 lmle(:,:) = .FALSE. 648 DO jj = 2, jpjm1 649 DO ji = 2, jpim1 650 IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 651 END DO 652 END DO 653 ENDIF ! ln_osm_mle 654 655 ! Test if pycnocline well resolved 656 DO jj = 2, jpjm1 657 DO ji = 2,jpim1 658 IF (lconv(ji,jj) ) THEN 659 ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 660 IF ( ztmp > 6 ) THEN 661 ! pycnocline well resolved 662 jp_ext(ji,jj) = 1 663 ELSE 664 ! pycnocline poorly resolved 665 jp_ext(ji,jj) = 0 666 ENDIF 667 ELSE 668 ! Stable conditions 669 jp_ext(ji,jj) = 0 670 ENDIF 671 END DO 672 END DO 673 674 CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 675 ! jp_ext = ibld-imld+1 676 CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 601 677 ! Rate of change of hbl 602 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 603 ! Calculate averages over depth of boundary layer 604 DO jj = 2, jpjm1 605 DO ji = 2, jpim1 606 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 678 CALL zdf_osm_calculate_dhdt( zdhdt ) 679 DO jj = 2, jpjm1 680 DO ji = 2, jpim1 681 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 607 682 ! adjustment to represent limiting by ocean bottom 608 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 609 END DO 683 IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN 684 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 685 lpyc(ji,jj) = .FALSE. 686 ENDIF 610 687 END DO 688 END DO 611 689 612 690 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index … … 626 704 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 627 705 ! 628 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 629 ! Alan: do we need zb_ml? 630 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 706 CALL zdf_osm_timestep_hbl( zdhdt ) 707 ! is external level in bounds? 708 709 CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 631 710 ! 632 711 ! 712 ! Check to see if lpyc needs to be changed 713 633 714 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 715 716 DO jj = 2, jpjm1 717 DO ji = 2, jpim1 718 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE. 719 END DO 720 END DO 721 634 722 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 635 723 ! 636 724 ! Average over the depth of the mixed layer in the convective boundary layer 637 ! Alan: do we need zb_ml? 638 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )725 ! jp_ext = ibld - imld +1 726 CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 639 727 ! rotate mean currents and changes onto wind align co-ordinates 640 728 ! 641 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 642 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 643 zuw_bse = 0._wp 644 zvw_bse = 0._wp 645 zwth_ent = 0._wp 646 zws_ent = 0._wp 647 DO jj = 2, jpjm1 648 DO ji = 2, jpim1 649 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 650 IF ( lconv(ji,jj) ) THEN 651 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 652 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 653 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 654 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 655 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 656 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 657 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 658 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 659 ENDIF 660 ELSE 661 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 662 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 663 ENDIF 664 ENDIF 665 END DO 666 END DO 667 729 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 730 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 668 731 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 669 732 ! Pycnocline gradients for scalars and velocity 670 733 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 671 734 672 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )673 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc )735 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 736 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 674 737 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 675 738 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 676 739 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 677 740 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 678 679 DO jj = 2, jpjm1 680 DO ji = 2, jpim1 681 IF ( lconv(ji,jj) ) THEN 682 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 683 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 684 zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 685 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 686 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 687 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 688 ELSE 689 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 690 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 691 END IF 692 END DO 693 END DO 694 ! 695 DO jj = 2, jpjm1 696 DO ji = 2, jpim1 697 IF ( lconv(ji,jj) ) THEN 698 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 699 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 700 ! 701 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 702 ! 703 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 704 & * ( 1.0 - 0.5 * zznd_ml**2 ) 705 END DO 706 ! pycnocline - if present linear profile 707 IF ( zdh(ji,jj) > 0._wp ) THEN 708 zgamma_b = 6.0 709 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 710 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 711 ! 712 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 713 ! 714 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 715 END DO 716 IF ( ibld_ext == 0 ) THEN 717 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 718 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 719 ELSE 720 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 721 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 722 ENDIF 723 ENDIF 724 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 725 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 726 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 727 ! could be taken out, take account of entrainment represents as a diffusivity 728 ! should remove w from here, represents entrainment 729 ELSE 730 ! stable conditions 731 DO jk = 2, ibld(ji,jj) 732 zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 733 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 734 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 735 END DO 736 737 IF ( ibld_ext == 0 ) THEN 738 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 739 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 740 ELSE 741 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 742 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 743 ENDIF 744 ENDIF ! end if ( lconv ) 745 ! 746 END DO ! end of ji loop 747 END DO ! end of jj loop 741 CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 748 742 749 743 ! … … 775 769 DO jk = 2, ibld(ji,jj) 776 770 zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj) 777 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9* zznd_d ) &771 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & 778 772 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 779 773 ! 780 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9* zznd_d ) &774 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & 781 775 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 782 776 END DO … … 826 820 827 821 WHERE ( lconv ) 828 zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )829 zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )822 zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 823 zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 830 824 ELSEWHERE 831 825 zsc_wth_1 = 0._wp … … 838 832 DO jk = 2, imld(ji,jj) 839 833 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 840 ! calculate turbulent lengthscale841 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml -zznd_ml**3 / 3.0 ) ) ) &842 & * ( 1.0 - EXP ( -15.0 * ( 1. 1- zznd_ml ) ) )843 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml -zznd_ml**3 / 3.0 ) ) ) &844 & * ( 1.0 - EXP ( - 5.0 * ( 1.0- zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )845 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 /2.0)834 ! calculate turbulent time scale 835 zl_c = 0.9 * ( 1.0 - EXP ( - 5.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & 836 & * ( 1.0 - EXP ( -15.0 * ( 1.2 - zznd_ml ) ) ) 837 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & 838 & * ( 1.0 - EXP ( - 8.0 * ( 1.15 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 839 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 846 840 ! non-gradient buoyancy terms 847 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0. 5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj)/ ( 0.15 + zznd_ml )848 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0. 5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj)/ ( 0.15 + zznd_ml )841 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.4 * zsc_wth_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) 842 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 * zsc_ws_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) 849 843 END DO 850 ELSE 844 845 IF ( lpyc(ji,jj) ) THEN 846 ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 847 ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 ) 848 zwth_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) 849 zws_ent = -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 850 ! Cubic profile used for buoyancy term 851 DO jk = 2, ibld(ji,jj) 852 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 853 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 854 855 ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 ) 856 END DO 857 ! 858 IF ( dh(ji,jj) < 0.2*hbl(ji,jj) ) THEN 859 zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 860 zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) ) 861 ! 862 zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 863 ! 864 zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj) 865 ! 866 zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 867 DO jk = 2, ibld(ji,jj) 868 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 869 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 870 ! 871 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 872 END DO 873 END IF 874 ENDIF ! End of pycnocline 875 ELSE ! lconv test - stable conditions 851 876 DO jk = 2, ibld(ji,jj) 852 877 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) … … 886 911 END DO ! jj loop 887 912 913 DO jj = 2, jpjm1 914 DO ji = 2, jpim1 915 IF( lconv(ji,jj) ) THEN 916 IF ( lpyc(ji,jj) ) THEN 917 IF ( j_ddh(ji,jj) == 0 ) THEN 918 ! Place holding code. Parametrization needs checking for these conditions. 919 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird 920 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 921 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 922 ELSE 923 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird 924 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 925 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 926 ENDIF 927 zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 928 zc_cubic = zuw_bse - zd_cubic 929 ! need ztau_sc_u to be available. Change to array. 930 DO jk = imld(ji,jj), ibld(ji,jj) 931 zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 932 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse * & 933 & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 934 END DO 935 zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 936 zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 937 zc_cubic = zvw_bse - zd_cubic 938 DO jk = imld(ji,jj), ibld(ji,jj) 939 zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj) 940 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse * & 941 & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 942 END DO 943 ENDIF ! lpyc 944 ENDIF ! lconv 945 END DO ! ji loop 946 END DO ! jj loop 947 888 948 IF(ln_dia_osm) THEN 889 949 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) … … 892 952 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 893 953 894 WHERE ( lconv ) 895 zsc_wth_1 = zwth0 896 zsc_ws_1 = zws0 897 ELSEWHERE 898 zsc_wth_1 = 2.0 * zwthav 899 zsc_ws_1 = zws0 900 ENDWHERE 954 DO jj = 1, jpjm1 955 DO ji = 1, jpim1 956 957 IF ( lconv(ji,jj) ) THEN 958 zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 959 zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 960 IF ( lpyc(ji,jj) ) THEN 961 ! Pycnocline scales 962 zsc_wth_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj) 963 zsc_ws_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj) 964 ENDIF 965 ELSE 966 zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 967 zsc_ws_1(ji,jj) = zws0(ji,jj) 968 ENDIF 969 END DO 970 END DO 901 971 902 972 DO jj = 2, jpjm1 … … 915 985 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 916 986 END DO 987 ! 988 ! may need to comment out lpyc block 989 IF ( lpyc(ji,jj) ) THEN 990 ! pycnocline 991 DO jk = imld(ji,jj), ibld(ji,jj) 992 zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 993 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 994 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 995 END DO 996 ENDIF 917 997 ELSE 918 DO jk = 2, ibld(ji,jj) 919 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 920 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 921 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 922 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 923 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 924 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 925 END DO 998 IF( zdhdt(ji,jj) > 0. ) THEN 999 DO jk = 2, ibld(ji,jj) 1000 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 1001 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1002 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1003 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 1004 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1005 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 1006 END DO 1007 ENDIF 926 1008 ENDIF 927 1009 ENDDO ! ji loop … … 950 1032 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 951 1033 END DO 1034 952 1035 ELSE 953 1036 DO jk = 2, ibld(ji,jj) … … 973 1056 END DO 974 1057 975 IF(ln_dia_osm) THEN 1058 ! DO jj = 1, jpjm1 1059 ! DO ji = 1, jpim1 1060 ! IF ( lconv(ji,jj) ) THEN 1061 ! IF ( lpyc(ji,jj) ) THEN 1062 ! zd_cubic = ( 0.948 - 2.13 * zdh(ji,jj) / zhml(ji,jj) ) * zustar(ji,jj)**2 1063 ! zc_cubic = -0.474 * zustar(ji,jj)**2 - zd_cubic 1064 ! DO jk = imld(ji,jj), ibld(ji,jj) 1065 ! zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 1066 ! ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1067 ! END DO 1068 ! zc_cubic= 3.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj) 1069 ! zd_cubic = -2.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj) 1070 ! DO jk = imld(ji,jj), ibld(ji,jj) 1071 ! zznd_pyc = -( gdepw_n(ji,jj,jk)-zhbl(ji,jj) ) / zdh(ji,jj) 1072 ! ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1073 ! END DO 1074 ! ENDIF 1075 ! ENDIF 1076 ! END DO 1077 ! END DO 1078 1079 IF(ln_dSia_osm) THEN 976 1080 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 977 1081 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) … … 984 1088 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 985 1089 1090 1091 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 1092 986 1093 DO jj = 2, jpjm1 987 1094 DO ji = 2, jpim1 988 IF ( lconv(ji,jj) ) THEN1095 IF ( .not. lconv(ji,jj) ) THEN 989 1096 DO jk = 2, ibld(ji,jj) 990 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 991 IF ( znd >= 0.0 ) THEN 992 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 993 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 994 ELSE 995 ghamu(ji,jj,jk) = 0._wp 996 ghamv(ji,jj,jk) = 0._wp 997 ENDIF 998 END DO 999 ELSE 1000 DO jk = 2, ibld(ji,jj) 1001 znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1097 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 1002 1098 IF ( znd >= 0.0 ) THEN 1003 1099 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) … … 1012 1108 END DO 1013 1109 1014 IF(ln_dia_osm) THEN 1110 ! pynocline contributions 1111 DO jj = 2, jpjm1 1112 DO ji = 2, jpim1 1113 IF ( .not. lconv(ji,jj) ) THEN 1114 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1115 DO jk= 2, ibld(ji,jj) 1116 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1117 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1118 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1119 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1120 END DO 1121 END IF 1122 END IF 1123 END DO 1124 END DO 1125 IF(ln_dia_osm) THEN 1015 1126 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1016 1127 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1017 1128 END IF 1018 ! pynocline contributions1019 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small1020 zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )1021 DO jj = 2, jpjm11022 DO ji = 2, jpim11023 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN1024 DO jk= 2, ibld(ji,jj)1025 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)1026 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)1027 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)1028 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)1029 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk)1030 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)1031 END DO1032 END IF1033 END DO1034 END DO1035 1036 ! Entrainment contribution.1037 1129 1038 1130 DO jj=2, jpjm1 1039 1131 DO ji = 2, jpim1 1040 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 1041 DO jk = 1, imld(ji,jj) - 1 1042 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1043 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 1044 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1045 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1046 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 1047 END DO 1048 DO jk = imld(ji,jj), ibld(ji,jj) 1049 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 1050 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1051 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1052 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1053 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1054 END DO 1055 ENDIF 1056 1057 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1058 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1059 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1060 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1132 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1133 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1134 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1135 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 1061 1136 END DO ! ji loop 1062 1137 END DO ! jj loop … … 1065 1140 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 1066 1141 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 1067 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse )1068 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse )1069 1142 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 1070 1143 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) … … 1155 1228 DO jj = 2 , jpjm1 1156 1229 DO ji = 2, jpim1 1157 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1158 ! Calculate MLE flux profiles 1159 ! DO jk = 1, mld_prof(ji,jj) 1160 ! znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1161 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1162 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1163 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1164 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1165 ! END DO 1166 ! Calculate MLE flux contribution from surface fluxes 1230 IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 1231 ! Calculate MLE flux contribution from surface fluxes 1167 1232 DO jk = 1, ibld(ji,jj) 1168 1233 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 1169 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd )1234 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - zrad0(ji,jj) ) * ( 1.0 - znd ) 1170 1235 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1171 1236 END DO 1172 1237 DO jk = 1, mld_prof(ji,jj) 1173 1238 znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1174 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd )1239 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( zwth0(ji,jj) - zrad0(ji,jj) ) * ( 1.0 - znd ) 1175 1240 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1176 1241 END DO 1177 1242 ! Viscosity for MLEs 1178 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1179 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1243 DO jk = 1, mld_prof(ji,jj) 1244 znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1245 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1180 1246 END DO 1181 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1182 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1183 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1184 DO jk = mld_prof(ji,jj), jl 1185 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1186 & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 1187 & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 1247 ELSE 1248 ! Surface transports limited to OSBL. 1249 ! Viscosity for MLEs 1250 DO jk = 1, mld_prof(ji,jj) 1251 znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1252 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 ) 1188 1253 END DO 1189 ENDIF1254 ENDIF 1190 1255 END DO 1191 1256 END DO … … 1261 1326 IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0> 1262 1327 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1328 IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*zwb0 ) ! <Sw_0> 1329 IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb buoyancy flux 1263 1330 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1264 1331 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k … … 1270 1337 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1271 1338 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1339 IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base 1340 IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base 1341 IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base 1272 1342 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1273 1343 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1282 1352 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1283 1353 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1354 IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolved internal to zdf_osm routine 1355 IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*j_ddh ) ! index forpyc thicknessh internal to zdf_osm routine 1356 IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKE internal to zdf_osm routine 1284 1357 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1285 1358 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1286 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux1287 1359 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1288 1360 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux … … 1307 1379 1308 1380 CONTAINS 1309 1310 1311 ! Alan: do we need zb? 1312 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 1381 ! subroutine code changed, needs syntax checking. 1382 SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 1383 1384 !!--------------------------------------------------------------------- 1385 !! *** ROUTINE zdf_osm_diffusivity_viscosity *** 1386 !! 1387 !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. 1388 !! 1389 !! ** Method : 1390 !! 1391 !! !!---------------------------------------------------------------------- 1392 REAL(wp), DIMENSION(:,:,:) :: zdiffut 1393 REAL(wp), DIMENSION(:,:,:) :: zviscos 1394 ! local 1395 1396 ! Scales used to calculate eddy diffusivity and viscosity profiles 1397 REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 1398 REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 1399 REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 1400 REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 1401 ! 1402 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 1403 1404 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 1405 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 1406 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 1407 1408 DO jj = 2, jpjm1 1409 DO ji = 2, jpim1 1410 IF ( lconv(ji,jj) ) THEN 1411 1412 zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 1413 zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1414 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 1415 1416 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 1417 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 1418 1419 IF ( lpyc(ji,jj) ) THEN 1420 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 1421 1422 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 1423 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 1424 ENDIF 1425 1426 zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) 1427 zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 1428 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 1429 1430 zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 1431 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 1432 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 1433 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 1434 ENDIF 1435 1436 zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 1437 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 1438 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_n_sc(ji,jj) ) 1439 1440 zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third 1441 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 1442 ELSE 1443 zbeta_d_sc(ji,jj) = 1.0 1444 zbeta_v_sc(ji,jj) = 1.0 1445 ENDIF 1446 ELSE 1447 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1448 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1449 END IF 1450 END DO 1451 END DO 1452 ! 1453 DO jj = 2, jpjm1 1454 DO ji = 2, jpim1 1455 IF ( lconv(ji,jj) ) THEN 1456 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 1457 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 1458 ! 1459 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 1460 ! 1461 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 1462 & * ( 1.0 - 0.5 * zznd_ml**2 ) 1463 END DO 1464 ! pycnocline 1465 IF ( lpyc(ji,jj) ) THEN 1466 ! Diffusivity profile in the pycnocline given by cubic polynomial. 1467 za_cubic = 0.5 1468 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 1469 zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) & 1470 & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 1471 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 1472 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1473 DO jk = imld(ji,jj) , ibld(ji,jj) 1474 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1475 ! 1476 zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1477 1478 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ) 1479 END DO 1480 ! viscosity profiles. 1481 za_cubic = 0.5 1482 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 1483 zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj) ) / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 1484 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 1485 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1486 DO jk = imld(ji,jj) , ibld(ji,jj) 1487 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1488 zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1489 zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 ) 1490 END DO 1491 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1492 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 1493 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 1494 ELSE 1495 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 1496 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 1497 ENDIF 1498 ENDIF 1499 ELSE 1500 ! stable conditions 1501 DO jk = 2, ibld(ji,jj) 1502 zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1503 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 1504 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 1505 END DO 1506 1507 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1508 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 1509 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 1510 ENDIF 1511 ENDIF ! end if ( lconv ) 1512 ! 1513 END DO ! end of ji loop 1514 END DO ! end of jj loop 1515 1516 END SUBROUTINE zdf_osm_diffusivity_viscosity 1517 1518 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 1519 1520 !!--------------------------------------------------------------------- 1521 !! *** ROUTINE zdf_osm_osbl_state *** 1522 !! 1523 !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number 1524 !! 1525 !! ** Method : 1526 !! 1527 !! !!---------------------------------------------------------------------- 1528 1529 INTEGER, DIMENSION(jpi,jpj) :: j_ddh ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low. 1530 1531 LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear 1532 1533 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 1534 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline 1535 1536 ! Local Variables 1537 1538 INTEGER :: jj, ji 1539 1540 REAL(wp), DIMENSION(jpi,jpj) :: zekman 1541 REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers 1542 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1543 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1544 1545 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8 1546 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03 1547 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 1548 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 1549 REAL, PARAMETER :: zri_c = 0.25 1550 REAL, PARAMETER :: zek = 4.0 1551 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. 1552 1553 ! Determins stability and set flag lconv 1554 DO jj = 2, jpjm1 1555 DO ji = 2, jpim1 1556 IF ( zhol(ji,jj) < 0._wp ) THEN 1557 lconv(ji,jj) = .TRUE. 1558 ELSE 1559 lconv(ji,jj) = .FALSE. 1560 ENDIF 1561 END DO 1562 END DO 1563 1564 zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1565 1566 zshear(:,:) = 0._wp 1567 j_ddh(:,:) = 1 1568 1569 DO jj = 2, jpjm1 1570 DO ji = 2, jpim1 1571 IF ( lconv(ji,jj) ) THEN 1572 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1573 zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 1574 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1575 1576 IF ( ff_t(ji,jj) >= 0._wp ) THEN 1577 ! Northern Hemisphere 1578 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )**2 + MAX( -zdv_ml(ji,jj), 1.e-5)**2 ) 1579 ELSE 1580 ! Southern Hemisphere 1581 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )**2 + MAX( zdv_ml(ji,jj), 1.e-5)**2 ) 1582 ENDIF 1583 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 1584 ! Stability Dependence 1585 zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) ) 1586 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1587 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! 1588 ! full code available ! 1589 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1590 IF ( zshear(ji,jj) > 1.e-10 ) THEN 1591 IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN 1592 ! Growing shear layer 1593 j_ddh(ji,jj) = 0 1594 lshear(ji,jj) = .TRUE. 1595 ELSE 1596 j_ddh(ji,jj) = 1 1597 ! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 1598 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 1599 lshear(ji,jj) = .TRUE. 1600 ! ELSE 1601 ENDIF 1602 ELSE 1603 j_ddh(ji,jj) = 2 1604 lshear(ji,jj) = .FALSE. 1605 ENDIF 1606 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 1607 ! zshear(ji,jj) = 0.5 * zshear(ji,jj) 1608 ! lshear(ji,jj) = .FALSE. 1609 ! ENDIF 1610 ELSE ! zdb_bl test, note zshear set to zero 1611 j_ddh(ji,jj) = 2 1612 lshear(ji,jj) = .FALSE. 1613 ENDIF 1614 ENDIF 1615 END DO 1616 END DO 1617 1618 ! Calculate entrainment buoyancy flux due to surface fluxes. 1619 1620 DO jj = 2, jpjm1 1621 DO ji = 2, jpim1 1622 IF ( lconv(ji,jj) ) THEN 1623 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1624 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1625 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1626 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1627 IF (nn_osm_SD_reduce > 0 ) THEN 1628 ! Effective Stokes drift already reduced from surface value 1629 zr_stokes = 1.0_wp 1630 ELSE 1631 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 1632 ! requires further reduction where BL is deep 1633 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1634 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1635 END IF 1636 zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) & 1637 & - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1638 & + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1639 & - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1640 ! 1641 ENDIF 1642 END DO ! ji loop 1643 END DO ! jj loop 1644 1645 zwb_min(:,:) = 0._wp 1646 1647 DO jj = 2, jpjm1 1648 DO ji = 2, jpim1 1649 IF ( lshear(ji,jj) ) THEN 1650 IF ( lconv(ji,jj) ) THEN 1651 ! Unstable OSBL 1652 zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj) 1653 IF ( j_ddh(ji,jj) == 0 ) THEN 1654 1655 ! ! Developing shear layer, additional shear production possible. 1656 1657 ! zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp ) 1658 ! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) 1659 ! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1660 1661 ! zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) 1662 ! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u ) 1663 1664 ENDIF 1665 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1666 ! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1667 ELSE ! IF ( lconv ) THEN - ENDIF 1668 ! Stable OSBL - shear production not coded for first attempt. 1669 ENDIF ! lconv 1670 ENDIF ! lshear 1671 IF ( lconv(ji,jj) ) THEN 1672 ! Unstable OSBL 1673 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0tot(ji,jj) 1674 ENDIF ! lconv 1675 END DO ! ji 1676 END DO ! jj 1677 END SUBROUTINE zdf_osm_osbl_state 1678 1679 1680 SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) 1313 1681 !!--------------------------------------------------------------------- 1314 1682 !! *** ROUTINE zdf_vertical_average *** … … 1322 1690 1323 1691 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1692 INTEGER, DIMENSION(jpi,jpj) :: jp_ext 1324 1693 1325 1694 ! Alan: do we need zb? 1326 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity1695 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity 1327 1696 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1328 1697 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1329 1698 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1330 1699 1331 INTEGER :: jk, ji, jj 1700 INTEGER :: jk, ji, jj, ibld_ext 1332 1701 REAL(wp) :: zthick, zthermal, zbeta 1333 1702 … … 1358 1727 zu(ji,jj) = zu(ji,jj) / zthick 1359 1728 zv(ji,jj) = zv(ji,jj) / zthick 1360 ! Alan: do we need zb? 1361 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1362 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1363 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1364 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1365 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1366 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1367 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1729 zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 1730 ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 1731 IF ( ibld_ext < mbkt(ji,jj) ) THEN 1732 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem) 1733 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal) 1734 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) & 1735 & / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) 1736 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) & 1737 & / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 1738 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1739 ELSE 1740 zdt(ji,jj) = 0._wp 1741 zds(ji,jj) = 0._wp 1742 zdu(ji,jj) = 0._wp 1743 zdv(ji,jj) = 0._wp 1744 zdb(ji,jj) = 0._wp 1745 ENDIF 1368 1746 END DO 1369 1747 END DO … … 1394 1772 ztemp = zdu(ji,jj) 1395 1773 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1396 zdv(ji,jj) = zdv(ji,jj) * z sin_w(ji,jj) - ztemp * zsin_w(ji,jj)1774 zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1397 1775 END DO 1398 1776 END DO 1399 1777 END SUBROUTINE zdf_osm_velocity_rotation 1400 1778 1401 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1779 SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 1780 !!--------------------------------------------------------------------- 1781 !! *** ROUTINE zdf_osm_osbl_state_fk *** 1782 !! 1783 !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme. 1784 !! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined 1785 !! lflux :: determines whether effects of surface flux extend below the base of the OSBL 1786 !! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl. 1787 !! 1788 !! ** Method : 1789 !! 1790 !! 1791 !!---------------------------------------------------------------------- 1792 1793 ! Outputs 1794 LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle 1795 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk 1796 ! 1797 REAL(wp), DIMENSION(jpi,jpj) :: znd_param 1798 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer 1799 REAL(wp) :: zpe_mle_ref, zdbdz_mle_int 1800 1801 znd_param(:,:) = 0._wp 1802 1803 DO jj = 2, jpjm1 1804 DO ji = 2, jpim1 1805 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1806 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 1807 END DO 1808 END DO 1809 DO jj = 2, jpjm1 1810 DO ji = 2, jpim1 1811 ! 1812 IF ( lconv(ji,jj) ) THEN 1813 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1814 zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1815 zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1816 zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1817 zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1818 ! Calculate potential energies of actual profile and reference profile. 1819 zpe_mle_layer = 0._wp 1820 zpe_mle_ref = 0._wp 1821 zthermal = rab_n(ji,jj,1,jp_tem) 1822 zbeta = rab_n(ji,jj,1,jp_sal) 1823 1824 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1825 zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) ) 1826 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 1827 zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 1828 END DO 1829 ! Non-dimensional parameter to diagnose the presence of thermocline 1830 1831 znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) ) 1832 ENDIF 1833 ENDIF 1834 END DO 1835 END DO 1836 1837 ! Diagnosis 1838 DO jj = 2, jpjm1 1839 DO ji = 2, jpim1 1840 IF ( lconv(ji,jj) ) THEN 1841 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN 1842 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1843 ! MLE layer growing 1844 IF ( znd_param (ji,jj) > 100. ) THEN 1845 ! Thermocline present 1846 lflux(ji,jj) = .FALSE. 1847 lmle(ji,jj) =.FALSE. 1848 ELSE 1849 ! Thermocline not present 1850 lflux(ji,jj) = .TRUE. 1851 lmle(ji,jj) = .TRUE. 1852 ENDIF ! znd_param > 100 1853 ! 1854 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1855 lpyc(ji,jj) = .FALSE. 1856 ELSE 1857 lpyc(ji,jj) = .TRUE. 1858 ENDIF 1859 ELSE 1860 ! MLE layer restricted to OSBL or just below. 1861 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1862 ! Weak stratification MLE layer can grow. 1863 lpyc(ji,jj) = .FALSE. 1864 lflux(ji,jj) = .TRUE. 1865 lmle(ji,jj) = .TRUE. 1866 ELSE 1867 ! Strong stratification 1868 lpyc(ji,jj) = .TRUE. 1869 lflux(ji,jj) = .FALSE. 1870 lmle(ji,jj) = .FALSE. 1871 ENDIF ! zdb_bl < rn_mle_thresh_bl and 1872 ENDIF ! zhmle > 1.2 zhbl 1873 ELSE 1874 lpyc(ji,jj) = .TRUE. 1875 lflux(ji,jj) = .FALSE. 1876 lmle(ji,jj) = .FALSE. 1877 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 1878 ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 1879 ELSE 1880 ! Stable Boundary Layer 1881 lpyc(ji,jj) = .FALSE. 1882 lflux(ji,jj) = .FALSE. 1883 lmle(ji,jj) = .FALSE. 1884 ENDIF ! lconv 1885 END DO 1886 END DO 1887 END SUBROUTINE zdf_osm_osbl_state_fk 1888 1889 SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 1402 1890 !!--------------------------------------------------------------------- 1403 1891 !! *** ROUTINE zdf_osm_external_gradients *** … … 1409 1897 !!---------------------------------------------------------------------- 1410 1898 1899 INTEGER, DIMENSION(jpi,jpj) :: jbase 1411 1900 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1412 1901 … … 1417 1906 DO jj = 2, jpjm1 1418 1907 DO ji = 2, jpim1 1419 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN1908 IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1420 1909 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1421 1910 zbeta = rab_n(ji,jj,1,jp_sal) 1422 jkb = ibld(ji,jj) + ibld_ext1911 jkb = jbase(ji,jj) 1423 1912 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1424 1913 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1425 & / e3t_n(ji,jj, ibld(ji,jj))1914 & / e3t_n(ji,jj,jkb) 1426 1915 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1427 & / e3t_n(ji,jj, ibld(ji,jj))1916 & / e3t_n(ji,jj,jkb) 1428 1917 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1918 ELSE 1919 zdtdz(ji,jj) = 0._wp 1920 zdsdz(ji,jj) = 0._wp 1921 zdbdz(ji,jj) = 0._wp 1429 1922 END IF 1430 1923 END DO … … 1432 1925 END SUBROUTINE zdf_osm_external_gradients 1433 1926 1434 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )1927 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 1435 1928 1436 1929 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1930 REAL(wp), DIMENSION(jpi,jpj) :: zalpha 1437 1931 1438 1932 INTEGER :: jk, jj, ji 1439 1933 REAL(wp) :: ztgrad, zsgrad, zbgrad 1440 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1441 REAL(wp) :: zzeta_s=0.3, ztmp 1934 REAL(wp) :: zgamma_b_nd, znd 1935 REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 1936 REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 1442 1937 1443 1938 DO jj = 2, jpjm1 1444 1939 DO ji = 2, jpim1 1445 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN1940 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1446 1941 IF ( lconv(ji,jj) ) THEN ! convective conditions 1447 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1448 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1449 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1450 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1451 ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 1452 zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 1453 zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 1454 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1455 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1456 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1457 DO jk = 2, ibld(ji,jj)+ibld_ext 1458 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1459 IF ( znd <= zzeta_s ) THEN 1460 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 1461 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1462 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 1463 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1464 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 1465 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1466 ELSE 1467 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1468 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1469 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1470 ENDIF 1471 END DO 1472 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1942 IF ( lpyc(ji,jj) ) THEN 1943 zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 1944 zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) ) 1945 zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 1946 1947 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1948 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1949 ! Commented lines in this section are not needed in new code, once tested ! 1950 ! can be removed ! 1951 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1952 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 1953 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 1954 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1955 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1956 DO jk = 2, ibld(ji,jj) 1957 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1958 IF ( znd <= zzeta_m ) THEN 1959 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 1960 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1961 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 1962 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1963 zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 1964 & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1965 ELSE 1966 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1967 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1968 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1969 ENDIF 1970 END DO 1971 ENDIF ! if no pycnocline pycnocline gradients set to zero 1473 1972 ELSE 1474 1973 ! stable conditions 1475 1974 ! if pycnocline profile only defined when depth steady of increasing. 1476 IF ( zdhdt(ji,jj) > =0.0 ) THEN ! Depth increasing, or steady.1975 IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady. 1477 1976 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1478 1977 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline … … 1502 2001 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1503 2002 ENDIF ! IF (lconv) 1504 END IF ! IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) )2003 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 1505 2004 END DO 1506 2005 END DO … … 1527 2026 DO ji = 2, jpim1 1528 2027 ! 1529 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN2028 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1530 2029 IF ( lconv (ji,jj) ) THEN 1531 ! Unstable conditions 1532 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &1533 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &1534 2030 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 2031 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 2032 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 2033 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1535 2034 !Alan is this right? 1536 zvgrad = ( 0.7 * zdv_ml(ji,jj) + &1537 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &1538 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) &1539 & )/ (zdh(ji,jj) + epsln )1540 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext1541 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v1542 IF ( znd <= 0.0 ) THEN1543 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )1544 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )1545 ELSE1546 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )1547 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )1548 ENDIF1549 END DO2035 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 2036 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 2037 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 2038 ! & )/ (zdh(ji,jj) + epsln ) 2039 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 2040 ! znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 2041 ! IF ( znd <= 0.0 ) THEN 2042 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 2043 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 2044 ! ELSE 2045 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 2046 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 2047 ! ENDIF 2048 ! END DO 1550 2049 ELSE 1551 2050 ! stable conditions … … 1568 2067 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1569 2068 1570 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2)2069 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) 1571 2070 !!--------------------------------------------------------------------- 1572 2071 !! *** ROUTINE zdf_osm_calculate_dhdt *** … … 1578 2077 !!---------------------------------------------------------------------- 1579 2078 1580 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt , zdhdt_2! Rate of change of hbl2079 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl 1581 2080 1582 2081 INTEGER :: jj, ji 1583 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 1584 REAL(wp) :: zvel_max, zwb_min 1585 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1586 REAL(wp) :: zzeta_m = 0.3 1587 REAL(wp) :: zgamma_c = 2.0 1588 REAL(wp) :: zdhoh = 0.1 1589 REAL(wp) :: alpha_bc = 0.5 1590 1591 DO jj = 2, jpjm1 1592 DO ji = 2, jpim1 2082 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 2083 REAL(wp) :: zvel_max,zddhdt 2084 REAL(wp), PARAMETER :: zzeta_m = 0.3 2085 REAL(wp), PARAMETER :: zgamma_c = 2.0 2086 REAL(wp), PARAMETER :: zdhoh = 0.1 2087 REAL(wp), PARAMETER :: zalpha_b = 0.3 2088 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 2089 2090 DO jj = 2, jpjm1 2091 DO ji = 2, jpim1 2092 2093 IF ( lshear(ji,jj) ) THEN 1593 2094 IF ( lconv(ji,jj) ) THEN ! Convective 1594 ! Alan is this right? Yes, it's a bit different from the previous relationship1595 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &1596 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)1597 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln1598 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )1599 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )1600 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )1601 IF (nn_osm_SD_reduce > 0 ) THEN1602 ! Effective Stokes drift already reduced from surface value1603 zr_stokes = 1.0_wp1604 ELSE1605 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,1606 ! requires further reduction where BL is deep1607 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &1608 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )1609 END IF1610 1611 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &1612 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &1613 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &1614 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)1615 !1616 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)1617 2095 1618 2096 IF ( ln_osm_mle ) THEN 1619 ! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 1620 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 2097 1621 2098 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2099 ! Fox-Kemper buoyancy flux average over OSBL 1622 2100 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1623 2101 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) … … 1628 2106 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1629 2107 ! OSBL is deepening, entrainment > restratification 1630 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 2108 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 2109 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 2110 zpsi = -zalpha_pyc(ji,jj) * ( zwb0tot(ji,jj) - ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj) 2111 zpsi = zpsi - 4.0 * ( 1.0 + zdh(ji,jj) /zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) 2112 zpsi = zalpha_b * MAX ( zpsi, 0._wp ) 2113 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / MAX( zdb_bl(ji,jj), 1.e-15 ) 2114 IF ( j_ddh(ji,jj) == 1 ) THEN 2115 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 2116 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2117 ELSE 2118 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2119 ENDIF 2120 ! Relaxation to dh_ref = zari * hbl 2121 zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2122 2123 ELSE IF ( j_ddh(ji,jj) == 0 ) THEN 2124 ! Growing shear layer 2125 zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2126 zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 2127 ELSE 2128 zddhdt = 0._wp 2129 ENDIF ! j_ddh 2130 zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0 * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) * zddhdt / zdb_bl(ji,jj) 2131 ELSE ! zdb_bl >0 2132 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 2133 ENDIF 2134 ELSE ! zwb_min + 2*zwb_fk_b < 0 2135 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 2136 zdhdt(ji,jj) = - zvel_mle(ji,jj) 2137 2138 2139 ENDIF 2140 2141 ELSE 2142 ! Fox-Kemper not used. 2143 2144 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 2145 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 2146 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2147 ! added ajgn 23 July as temporay fix 2148 2149 ENDIF ! ln_osm_mle 2150 2151 ELSE ! lconv - Stable 2152 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 2153 IF ( zdhdt(ji,jj) < 0._wp ) THEN 2154 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 2155 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 2156 ELSE 2157 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2158 ENDIF 2159 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 2160 ENDIF ! lconv 2161 ELSE ! lshear 2162 IF ( lconv(ji,jj) ) THEN ! Convective 2163 2164 IF ( ln_osm_mle ) THEN 2165 2166 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2167 ! Fox-Kemper buoyancy flux average over OSBL 2168 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 2169 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 2170 ELSE 2171 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 2172 ENDIF 2173 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2174 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 2175 ! OSBL is deepening, entrainment > restratification 2176 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 1631 2177 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1632 2178 ELSE … … 1643 2189 ! Fox-Kemper not used. 1644 2190 1645 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) *zwb_ent(ji,jj) / &2191 zvel_max = -zwb_ent(ji,jj) / & 1646 2192 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 1647 2193 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1648 2194 ! added ajgn 23 July as temporay fix 1649 2195 1650 ENDIF 1651 1652 zdhdt_2(ji,jj) = 0._wp 1653 1654 ! commented out ajgn 23 July as temporay fix 1655 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1656 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 1657 ! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1658 ! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 1659 ! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1660 ! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 1661 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 1662 ! ! Alan no idea what this should be? 1663 ! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 1664 ! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 1665 ! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 1666 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 1667 ! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 1668 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1669 ! ENDIF 2196 ENDIF ! ln_osm_mle 2197 1670 2198 ELSE ! Stable 1671 2199 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1672 zdhdt_2(ji,jj) = 0._wp1673 2200 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1674 2201 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1675 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) *zvstr(ji,jj)**2 / hbl(ji,jj)2202 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 1676 2203 ELSE 1677 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) *zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )2204 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1678 2205 ENDIF 1679 2206 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 1680 ENDIF 1681 END DO 1682 END DO 2207 ENDIF ! lconv 2208 ENDIF ! lshear 2209 END DO 2210 END DO 1683 2211 END SUBROUTINE zdf_osm_calculate_dhdt 1684 2212 1685 SUBROUTINE zdf_osm_timestep_hbl( zdhdt , zdhdt_2)2213 SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 1686 2214 !!--------------------------------------------------------------------- 1687 2215 !! *** ROUTINE zdf_osm_timestep_hbl *** … … 1697 2225 1698 2226 1699 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt , zdhdt_2! rates of change of hbl.2227 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl. 1700 2228 1701 2229 INTEGER :: jk, jj, ji, jm … … 1735 2263 IF ( ln_osm_mle ) THEN 1736 2264 zhbl_s = zhbl_s + MIN( & 1737 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj)) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &2265 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1738 2266 & e3w_n(ji,jj,jm) ) 1739 2267 ELSE 1740 2268 zhbl_s = zhbl_s + MIN( & 1741 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj)) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &2269 & rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1742 2270 & e3w_n(ji,jj,jm) ) 1743 2271 ENDIF 1744 2272 1745 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1746 2273 ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2274 IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN 2275 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2276 lpyc(ji,jj) = .FALSE. 2277 ENDIF 1747 2278 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1748 2279 END DO … … 1765 2296 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1766 2297 1767 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2298 ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2299 IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 2300 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2301 lpyc(ji,jj) = .FALSE. 2302 ENDIF 1768 2303 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1769 2304 END DO … … 1796 2331 1797 2332 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1798 !2333 ! 1799 2334 INTEGER :: jj, ji 1800 2335 INTEGER :: inhml 1801 REAL(wp) :: zari, ztau, zddhdt 1802 1803 1804 DO jj = 2, jpjm1 1805 DO ji = 2, jpim1 2336 REAL(wp) :: zari, ztau, zdh_ref, zddhdt 2337 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 2338 2339 DO jj = 2, jpjm1 2340 DO ji = 2, jpim1 2341 2342 IF ( lshear(ji,jj) ) THEN 1806 2343 IF ( lconv(ji,jj) ) THEN 2344 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2345 IF ( j_ddh(ji,jj) == 0 ) THEN 2346 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2347 ! zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2348 ! zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 2349 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer 2350 ! dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) ) 2351 ztau = MAX( zdb_bl(ji,jj) * ( 0.625 * hbl(ji,jj) ) / ( a_ddh * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt ) 2352 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.625 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2353 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8 * zhbl(ji,jj) 2354 ELSE 2355 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 2356 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 2357 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2358 ELSE 2359 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 2360 ENDIF 2361 ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt ) 2362 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2363 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 2364 ENDIF 2365 ELSE 2366 ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0 * rn_rdt ) 2367 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2368 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(ji,jj) 2369 ENDIF 2370 ELSE ! lconv 2371 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 2372 2373 ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 2374 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2375 ! boundary layer deepening 2376 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2377 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2378 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 2379 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2380 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 2381 ELSE 2382 zdh_ref = 0.2 * hbl(ji,jj) 2383 ENDIF 2384 ELSE ! IF(dhdt < 0) 2385 zdh_ref = 0.2 * hbl(ji,jj) 2386 ENDIF ! IF (dhdt >= 0) 2387 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2388 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse 2389 ! Alan: this hml is never defined or used -- do we need it? 2390 ENDIF 2391 2392 ELSE ! lshear 2393 ! for lshear = .FALSE. calculate ddhdt here 2394 2395 IF ( lconv(ji,jj) ) THEN 1807 2396 1808 2397 IF( ln_osm_mle ) THEN 1809 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN2398 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1810 2399 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1811 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ ext(ji,jj) > 0._wp)THEN2400 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 1812 2401 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 1813 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1814 (1.0 + zdb_bl(ji,jj)**2 / MAX( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 2402 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1815 2403 ELSE ! unstable 1816 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1817 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 2404 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1818 2405 ENDIF 1819 2406 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1820 zd dhdt= zari * hbl(ji,jj)2407 zdh_ref = zari * hbl(ji,jj) 1821 2408 ELSE 1822 2409 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1823 zd dhdt= 0.2 * hbl(ji,jj)2410 zdh_ref = 0.2 * hbl(ji,jj) 1824 2411 ENDIF 1825 2412 ELSE 1826 2413 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1827 zd dhdt= 0.2 * hbl(ji,jj)2414 zdh_ref = 0.2 * hbl(ji,jj) 1828 2415 ENDIF 1829 2416 ELSE ! ln_osm_mle 1830 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ ext(ji,jj) > 0._wp)THEN2417 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 1831 2418 IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability 1832 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 1833 (1.0 + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 2419 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1834 2420 ELSE ! unstable 1835 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / & 1836 (1.0 + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 ) 2421 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 1837 2422 ENDIF 1838 2423 ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1839 zd dhdt= zari * hbl(ji,jj)2424 zdh_ref = zari * hbl(ji,jj) 1840 2425 ELSE 1841 2426 ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1842 zd dhdt= 0.2 * hbl(ji,jj)2427 zdh_ref = 0.2 * hbl(ji,jj) 1843 2428 ENDIF 1844 2429 1845 2430 END IF ! ln_osm_mle 1846 2431 1847 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2432 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2433 ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2434 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 1848 2435 ! Alan: this hml is never defined or used 1849 2436 ELSE ! IF (lconv) … … 1855 2442 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1856 2443 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 1857 zd dhdt= MIN( zari, 0.2 ) * hbl(ji,jj)2444 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 1858 2445 ELSE 1859 zd dhdt= 0.2 * hbl(ji,jj)2446 zdh_ref = 0.2 * hbl(ji,jj) 1860 2447 ENDIF 1861 2448 ELSE ! IF(dhdt < 0) 1862 zd dhdt= 0.2 * hbl(ji,jj)2449 zdh_ref = 0.2 * hbl(ji,jj) 1863 2450 ENDIF ! IF (dhdt >= 0) 1864 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1865 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1866 ! Alan: this hml is never defined or used -- do we need it? 2451 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2452 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse 1867 2453 ENDIF ! IF (lconv) 1868 1869 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1870 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1871 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1872 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1873 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1874 END DO 1875 END DO 2454 ENDIF ! lshear 2455 2456 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 2457 inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)), 1.e-3) ) , 1 ) 2458 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 2459 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 2460 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 2461 END DO 2462 END DO 1876 2463 1877 2464 END SUBROUTINE zdf_osm_pycnocline_thickness 1878 2465 1879 2466 1880 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )2467 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 1881 2468 !!---------------------------------------------------------------------- 1882 2469 !! *** ROUTINE zdf_osm_horizontal_gradients *** … … 1891 2478 1892 2479 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 2480 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 1893 2481 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 1894 2482 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy … … 1963 2551 END DO 1964 2552 2553 ! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on 2554 ztsm_midu(:,jpj,jp_sal) = 10. 2555 ztsm_midv(jpi,:,jp_sal) = 10. 2556 1965 2557 CALL eos_rab(ztsm_midu, zmld_midu, zabu) 1966 2558 CALL eos_rab(ztsm_midv, zmld_midv, zabv) … … 1977 2569 END DO 1978 2570 2571 DO jj = 2, jpjm1 2572 DO ji = 2, jpim1 2573 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2574 zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 2575 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 2576 END DO 2577 END DO 2578 1979 2579 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1980 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )2580 SUBROUTINE zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 1981 2581 !!---------------------------------------------------------------------- 1982 2582 !! *** ROUTINE zdf_osm_mle_parameters *** … … 1989 2589 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1990 2590 1991 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 2591 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 2592 INTEGER, DIMENSION(jpi,jpj) :: mld_prof 2593 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 1992 2594 INTEGER :: ji, jj, jk ! dummy loop indices 1993 INTEGER :: ii, ij, ik ! local integers2595 INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers 1994 2596 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 1995 REAL(wp) :: zdb_mle, ztmp, zdbds_mle 1996 1997 mld_prof(:,:) = 4 1998 DO jk = 5, jpkm1 1999 DO jj = 2, jpjm1 2000 DO ji = 2, jpim1 2001 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 2002 END DO 2597 REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle 2598 2599 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 2600 2601 DO jj = 2, jpjm1 2602 DO ji = 2, jpim1 2603 IF ( lconv(ji,jj) ) THEN 2604 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2605 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 2606 zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 2607 zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 2608 ENDIF 2003 2609 END DO 2004 2610 END DO 2005 ! DO jj = 2, jpjm12006 ! DO ji = 1, jpim12007 ! zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))2008 ! END DO2009 ! END DO2010 2611 ! Timestep mixed layer eddy depth. 2011 2612 DO jj = 2, jpjm1 2012 2613 DO ji = 2, jpim1 2013 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 2014 IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 2015 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 2016 ELSE 2017 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 2018 ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 2019 ! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 2020 ! ELSE 2021 ! hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 2022 ! ENDIF 2023 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2024 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 2025 ELSE 2026 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 2027 ENDIF 2028 ENDIF 2029 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 2030 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*zmld(ji,jj)) 2614 IF ( lmle(ji,jj) ) THEN ! MLE layer growing. 2615 ! Buoyancy gradient at base of MLE layer. 2616 zthermal = rab_n(ji,jj,1,jp_tem) 2617 zbeta = rab_n(ji,jj,1,jp_sal) 2618 jkb = mld_prof(ji,jj) 2619 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 2620 ! 2621 zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) ) 2622 zdb_mle = zb_bl(ji,jj) - zbuoy 2623 ! Timestep hmle. 2624 hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_rdt / zdb_mle 2625 ELSE 2626 IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN 2627 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 2628 ELSE 2629 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau 2630 ENDIF 2631 ENDIF 2632 hmle(ji,jj) = MAX(MIN(hmle(ji,jj), ht_n(ji,jj)), gdepw_n(ji,jj,4)) 2633 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) ) 2634 hmle(ji,jj) = zmld(ji,jj) 2031 2635 END DO 2032 2636 END DO … … 2045 2649 END DO 2046 2650 END DO 2047 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.2048 2049 DO jj = 2, jpjm12050 DO ji = 2, jpim12051 IF ( lconv(ji,jj) ) THEN2052 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf2053 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) &2054 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) )2055 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) &2056 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) )2057 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &2058 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )2059 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle2060 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.2061 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1)2062 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf2063 ENDIF2064 END DO2065 END DO2066 2651 END SUBROUTINE zdf_osm_mle_parameters 2067 2652 … … 2088 2673 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 2089 2674 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 2090 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, ln_zdfosm_ice_shelter2675 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 2091 2676 ! Namelist for Fox-Kemper parametrization. 2092 2677 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& … … 2138 2723 WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 2139 2724 WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm 2725 WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' 2140 2726 WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix 2141 2727 WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty … … 2144 2730 WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv 2145 2731 ENDIF 2732 2733 2734 ! ! Check wave coupling settings ! 2735 ! ! Further work needed - see ticket #2447 ! 2736 IF( nn_osm_wave == 2 ) THEN 2737 IF (.NOT. ( ln_wave .AND. ln_sdw )) & 2738 & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 2739 END IF 2146 2740 2147 2741 ! ! allocate zdfosm arrays … … 2171 2765 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2172 2766 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2173 WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s'2767 WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2174 2768 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2175 2769 WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
Note: See TracChangeset
for help on using the changeset viewer.
