Changeset 13858
- Timestamp:
- 2020-11-24T10:45:09+01:00 (3 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/cfgs/SHARED/namelist_ref
r13400 r13858 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 … … 1122 1123 rn_osm_mle_lat = 20. ! reference latitude (degrees) of MLE coef. (case rn_mle=1) 1123 1124 rn_osm_mle_rho_c = 0.01 ! delta rho criterion used to calculate MLD for FK 1124 rn_osm_mle_thresh = 0.0005 ! delta b criterion used for FK MLDcriterion1125 rn_osm_mle_thresh = 0.0005 ! delta b criterion used for FK MLE criterion 1125 1126 rn_osm_mle_tau = 172800. ! time scale for FK-OSM (seconds) (case rn_osm_mle=0) 1126 1127 ln_osm_hmle_limit = .false. ! limit hmle to rn_osm_hmle_limit*hbl -
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90
r13684 r13858 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 … … 249 251 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 250 252 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 253 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 251 254 252 255 … … 262 265 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 263 266 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 267 LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers 268 LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present 269 LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer. 270 LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness. 264 271 265 272 ! mixed-layer variables … … 267 274 INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base 268 275 INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top) 276 INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level 277 INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 269 278 270 279 REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients … … 279 288 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 280 289 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 281 REAL(wp), DIMENSION(jpi,jpj) :: zd hdt_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 gradients283 290 REAL(wp), DIMENSION(jpi,jpj) :: zddhdt ! correction to dhdt due to internal structure. 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, zri_i ! Shear production and interfacial richardon number. 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 … … 340 357 zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp 341 358 zhol(:,:) = 0._wp 342 lconv(:,:) = .FALSE. 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 zd hdt_2(:,:) = 0._wp397 zddhdt(:,:) = 0._wp 379 398 ! hbl = MAX(hbl,epsln) 380 399 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 548 567 zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 549 568 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 569 ELSE 552 570 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 553 lconv(ji,jj) = .FALSE.554 571 ENDIF 555 572 END DO … … 584 601 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 585 602 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 603 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 586 604 END DO 587 605 END DO 588 606 ! 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 607 jp_ext(:,:) = 2 608 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) 609 ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 610 CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 611 ! Velocity components in frame aligned with surface stress. 612 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 613 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 614 ! Determine the state of the OSBL, stable/unstable, shear/no shear 615 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 595 616 596 617 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 618 ! Fox-Kemper Scheme 619 mld_prof = 4 620 DO jk = 5, jpkm1 621 DO jj = 2, jpjm1 622 DO ji = 2, jpim1 623 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 624 END DO 625 END DO 626 END DO 627 jp_ext_mle(:,:) = 2 628 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) 629 630 DO jj = 2, jpjm1 631 DO ji = 2, jpim1 632 zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 633 END DO 634 END DO 635 636 !! External gradient 637 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 638 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 639 CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 640 CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 641 CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 642 ELSE ! ln_osm_mle 643 ! FK not selected, Boundary Layer only. 644 lpyc(:,:) = .TRUE. 645 lflux(:,:) = .FALSE. 646 lmle(:,:) = .FALSE. 647 DO jj = 2, jpjm1 648 DO ji = 2, jpim1 649 IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 650 END DO 651 END DO 652 ENDIF ! ln_osm_mle 653 654 ! Test if pycnocline well resolved 655 DO jj = 2, jpjm1 656 DO ji = 2,jpim1 657 IF (lconv(ji,jj) ) THEN 658 ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj)) 659 IF ( ztmp > 6 ) THEN 660 ! pycnocline well resolved 661 jp_ext(ji,jj) = 1 662 ELSE 663 ! pycnocline poorly resolved 664 jp_ext(ji,jj) = 0 665 ENDIF 666 ELSE 667 ! Stable conditions 668 jp_ext(ji,jj) = 0 669 ENDIF 670 END DO 671 END DO 672 673 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 ) 674 ! jp_ext = ibld-imld+1 675 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 676 ! 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 677 CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 678 DO jj = 2, jpjm1 679 DO ji = 2, jpim1 680 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 681 ! 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 682 IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN 683 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 684 lpyc(ji,jj) = .FALSE. 685 ENDIF 610 686 END DO 687 END DO 611 688 612 689 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index … … 626 703 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 627 704 ! 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 ) 705 CALL zdf_osm_timestep_hbl( zdhdt ) 706 ! is external level in bounds? 707 708 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 709 ! 632 710 ! 711 ! Check to see if lpyc needs to be changed 712 633 713 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 714 715 DO jj = 2, jpjm1 716 DO ji = 2, jpim1 717 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(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE. 718 END DO 719 END DO 720 634 721 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 635 722 ! 636 723 ! 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 )724 ! jp_ext = ibld - imld +1 725 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 726 ! rotate mean currents and changes onto wind align co-ordinates 640 727 ! 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 728 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 729 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 668 730 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 669 731 ! Pycnocline gradients for scalars and velocity 670 732 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 671 733 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 )734 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 735 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc ) 674 736 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 675 737 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 676 738 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 677 739 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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 740 CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 748 741 749 742 ! … … 843 836 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 844 837 & * ( 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)838 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 846 839 ! non-gradient buoyancy terms 847 840 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 841 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 ) 849 842 END DO 850 ELSE 843 844 IF ( lpyc(ji,jj) ) THEN 845 ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 846 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 ) 847 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) 848 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) 849 ! Cubic profile used for buoyancy term 850 za_cubic = 0.755 * ztau_sc_u(ji,jj) 851 zb_cubic = 0.25 * ztau_sc_u(ji,jj) 852 DO jk = 2, ibld(ji,jj) 853 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 854 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 ) 855 856 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 ) 857 END DO 858 ! 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 ENDIF ! End of pycnocline 874 ELSE ! lconv test - stable conditions 851 875 DO jk = 2, ibld(ji,jj) 852 876 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) … … 886 910 END DO ! jj loop 887 911 912 DO jj = 2, jpjm1 913 DO ji = 2, jpim1 914 IF ( lpyc(ji,jj) ) THEN 915 IF ( j_ddh(ji,jj) == 0 ) THEN 916 ! Place holding code. Parametrization needs checking for these conditions. 917 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 918 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 919 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 920 ELSE 921 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 922 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 923 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 924 ENDIF 925 zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 926 zc_cubic = zuw_bse - zd_cubic 927 ! need ztau_sc_u to be available. Change to array. 928 DO jk = imld(ji,jj), ibld(ji,jj) 929 zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 930 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 931 END DO 932 zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 933 zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 934 zc_cubic = zvw_bse - zd_cubic 935 DO jk = imld(ji,jj), ibld(ji,jj) 936 zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj) 937 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 938 END DO 939 ENDIF ! lpyc 940 END DO ! ji loop 941 END DO ! jj loop 942 888 943 IF(ln_dia_osm) THEN 889 944 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) … … 892 947 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 893 948 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 949 DO jj = 1, jpjm1 950 DO ji = 1, jpim1 951 952 IF ( lconv(ji,jj) ) THEN 953 zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 954 zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 955 IF ( lpyc(ji,jj) ) THEN 956 ! Pycnocline scales 957 zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 958 zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 959 ENDIF 960 ELSE 961 zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 962 zsc_ws_1(ji,jj) = zws0(ji,jj) 963 ENDIF 964 END DO 965 END DO 901 966 902 967 DO jj = 2, jpjm1 … … 915 980 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 916 981 END DO 982 ! 983 IF ( lpyc(ji,jj) ) THEN 984 ! pycnocline 985 DO jk = imld(ji,jj), ibld(ji,jj) 986 zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 987 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 ) ) 988 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 ) ) 989 END DO 990 ENDIF 917 991 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 992 IF( zdhdt(ji,jj) > 0. ) THEN 993 DO jk = 2, ibld(ji,jj) 994 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 995 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 996 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 997 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 998 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 999 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 1000 END DO 1001 ENDIF 926 1002 ENDIF 927 1003 ENDDO ! ji loop … … 984 1060 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 985 1061 1062 1063 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 1064 986 1065 DO jj = 2, jpjm1 987 1066 DO ji = 2, jpim1 988 IF ( lconv(ji,jj) ) THEN1067 IF ( .not. lconv(ji,jj) ) THEN 989 1068 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 1069 znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 1002 1070 IF ( znd >= 0.0 ) THEN 1003 1071 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) … … 1012 1080 END DO 1013 1081 1014 IF(ln_dia_osm) THEN1015 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )1016 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )1017 END IF1018 1082 ! pynocline contributions 1019 ! 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 1083 DO jj = 2, jpjm1 1022 1084 DO ji = 2, jpim1 1023 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1085 IF ( .not. lconv(ji,jj) ) THEN 1086 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1024 1087 DO jk= 2, ibld(ji,jj) 1025 1088 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) … … 1027 1090 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1028 1091 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 1092 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1031 1093 END DO 1032 1094 END IF 1095 END IF 1033 1096 END DO 1034 1097 END DO 1035 1036 ! Entrainment contribution. 1098 IF(ln_dia_osm) THEN 1099 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1100 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1101 END IF 1037 1102 1038 1103 DO jj=2, jpjm1 1039 1104 DO ji = 2, jpim1 1040 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN1041 DO jk = 1, imld(ji,jj) - 11042 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)1043 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd1044 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd1045 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd1046 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd1047 END DO1048 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 DO1055 ENDIF1056 1057 1105 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1058 1106 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp … … 1065 1113 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 1066 1114 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 1115 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 1070 1116 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) … … 1155 1201 DO jj = 2 , jpjm1 1156 1202 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 1203 IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 1204 ! Calculate MLE flux contribution from surface fluxes 1167 1205 DO jk = 1, ibld(ji,jj) 1168 1206 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) … … 1176 1214 END DO 1177 1215 ! 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) ) 1216 DO jk = 1, mld_prof(ji,jj) 1217 znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1218 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 1219 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)) 1220 ELSE 1221 ! Surface transports limited to OSBL. 1222 ! Viscosity for MLEs 1223 DO jk = 1, mld_prof(ji,jj) 1224 znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1225 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 1226 END DO 1189 ENDIF1227 ENDIF 1190 1228 END DO 1191 1229 END DO … … 1307 1345 1308 1346 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 ) 1347 ! subroutine code changed, needs syntax checking. 1348 SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 1349 1350 !!--------------------------------------------------------------------- 1351 !! *** ROUTINE zdf_osm_diffusivity_viscosity *** 1352 !! 1353 !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline. 1354 !! 1355 !! ** Method : 1356 !! 1357 !! !!---------------------------------------------------------------------- 1358 REAL(wp), DIMENSION(:,:,:) :: zdiffut 1359 REAL(wp), DIMENSION(:,:,:) :: zviscos 1360 ! local 1361 1362 ! Scales used to calculate eddy diffusivity and viscosity profiles 1363 REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc 1364 REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr 1365 REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 1366 REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 1367 ! 1368 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 1369 1370 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 1371 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142 1372 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 1373 1374 DO jj = 2, jpjm1 1375 DO ji = 2, jpim1 1376 IF ( lconv(ji,jj) ) THEN 1377 1378 zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird 1379 zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1380 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2 1381 1382 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml 1383 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 1384 1385 IF ( lpyc(ji,jj) ) THEN 1386 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 1387 1388 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 1389 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 1390 ENDIF 1391 1392 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) ) 1393 zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac 1394 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 1395 1396 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) 1397 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 1398 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN 1399 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 1400 ENDIF 1401 1402 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) ) ) 1403 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 1404 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) ) 1405 1406 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 1407 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 1408 ELSE 1409 zbeta_d_sc = 1.0 1410 zbeta_v_sc = 1.0 1411 ENDIF 1412 ELSE 1413 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1414 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 1415 END IF 1416 END DO 1417 END DO 1418 ! 1419 DO jj = 2, jpjm1 1420 DO ji = 2, jpim1 1421 IF ( lconv(ji,jj) ) THEN 1422 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 1423 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 1424 ! 1425 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 1426 ! 1427 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) & 1428 & * ( 1.0 - 0.5 * zznd_ml**2 ) 1429 END DO 1430 ! pycnocline 1431 IF ( lpyc(ji,jj) ) THEN 1432 ! Diffusivity profile in the pycnocline given by cubic polynomial. 1433 za_cubic = 0.5 1434 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 1435 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 ) & 1436 & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8) 1437 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 1438 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1439 DO jk = imld(ji,jj) , ibld(ji,jj) 1440 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1441 ! 1442 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 ) 1443 1444 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 ) 1445 END DO 1446 ! viscosity profiles. 1447 za_cubic = 0.5 1448 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 1449 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) 1450 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) 1451 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 1452 DO jk = imld(ji,jj) , ibld(ji,jj) 1453 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6) 1454 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 ) 1455 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 ) 1456 END DO 1457 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1458 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 1459 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 ) 1460 ELSE 1461 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 1462 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 1463 ENDIF 1464 ENDIF 1465 ELSE 1466 ! stable conditions 1467 DO jk = 2, ibld(ji,jj) 1468 zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1469 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 1470 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 1471 END DO 1472 1473 IF ( zdhdt(ji,jj) > 0._wp ) THEN 1474 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 1475 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj)) 1476 ENDIF 1477 ENDIF ! end if ( lconv ) 1478 ! 1479 END DO ! end of ji loop 1480 END DO ! end of jj loop 1481 1482 END SUBROUTINE zdf_osm_diffusivity_viscosity 1483 1484 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i ) 1485 1486 !!--------------------------------------------------------------------- 1487 !! *** ROUTINE zdf_osm_osbl_state *** 1488 !! 1489 !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number 1490 !! 1491 !! ** Method : 1492 !! 1493 !! !!---------------------------------------------------------------------- 1494 1495 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. 1496 1497 LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear 1498 1499 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 1500 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline 1501 REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number 1502 1503 ! Local Variables 1504 1505 INTEGER :: jj, ji 1506 1507 REAL(wp), DIMENSION(jpi,jpj) :: zekman 1508 REAL(wp) :: zri_p, zri_b ! Richardson numbers 1509 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1510 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1511 1512 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 1513 REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 1514 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 1515 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 1516 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 1517 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. 1518 1519 ! Determins stability and set flag lconv 1520 DO jj = 2, jpjm1 1521 DO ji = 2, jpim1 1522 IF ( zhol(ji,jj) < 0._wp ) THEN 1523 lconv(ji,jj) = .TRUE. 1524 ELSE 1525 lconv(ji,jj) = .FALSE. 1526 ENDIF 1527 END DO 1528 END DO 1529 1530 zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1531 1532 WHERE ( lconv ) 1533 zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 1534 END WHERE 1535 1536 zshear(:,:) = 0._wp 1537 j_ddh(:,:) = 1 1538 1539 DO jj = 2, jpjm1 1540 DO ji = 2, jpim1 1541 IF ( lconv(ji,jj) ) THEN 1542 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1543 zri_p = 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 & 1544 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1545 1546 zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 1547 1548 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 ) ) 1549 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1550 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! 1551 ! full code available ! 1552 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1553 IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 1554 ! Growing shear layer 1555 j_ddh(ji,jj) = 0 1556 lshear(ji,jj) = .TRUE. 1557 ELSE 1558 j_ddh(ji,jj) = 1 1559 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 1560 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 1561 lshear(ji,jj) = .TRUE. 1562 ELSE 1563 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 1564 zshear(ji,jj) = 0.5 * zshear(ji,jj) 1565 lshear(ji,jj) = .FALSE. 1566 ENDIF 1567 ENDIF 1568 ELSE ! zdb_bl test, note zshear set to zero 1569 j_ddh(ji,jj) = 2 1570 lshear(ji,jj) = .FALSE. 1571 ENDIF 1572 ENDIF 1573 END DO 1574 END DO 1575 1576 ! Calculate entrainment buoyancy flux due to surface fluxes. 1577 1578 DO jj = 2, jpjm1 1579 DO ji = 2, jpim1 1580 IF ( lconv(ji,jj) ) THEN 1581 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1582 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1583 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1584 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1585 IF (nn_osm_SD_reduce > 0 ) THEN 1586 ! Effective Stokes drift already reduced from surface value 1587 zr_stokes = 1.0_wp 1588 ELSE 1589 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 1590 ! requires further reduction where BL is deep 1591 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1592 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1593 END IF 1594 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1595 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1596 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1597 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1598 ! 1599 ENDIF 1600 END DO ! ji loop 1601 END DO ! jj loop 1602 1603 zwb_min(:,:) = 0._wp 1604 1605 DO jj = 2, jpjm1 1606 DO ji = 2, jpjm1 1607 IF ( lshear(ji,jj) ) THEN 1608 IF ( lconv(ji,jj) ) THEN 1609 ! Unstable OSBL 1610 zwb_shr = -za_wb_s * zshear(ji,jj) 1611 IF ( j_ddh(ji,jj) == 0 ) THEN 1612 1613 ! Developing shear layer, additional shear production possible. 1614 1615 zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) 1616 zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 1617 zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1618 1619 zwb_shr = -za_wb_s * zshear(ji,jj) 1620 1621 ENDIF 1622 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1623 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1624 ELSE ! IF ( lconv ) THEN - ENDIF 1625 ! Stable OSBL - shear production not coded for first attempt. 1626 ENDIF ! lconv 1627 ELSE ! lshear 1628 IF ( lconv(ji,jj) ) THEN 1629 ! Unstable OSBL 1630 zwb_shr = -za_wb_s * zshear(ji,jj) 1631 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1632 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1633 ENDIF ! lconv 1634 ENDIF ! lshear 1635 END DO ! ji 1636 END DO ! jj 1637 END SUBROUTINE zdf_osm_osbl_state 1638 1639 1640 SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv ) 1313 1641 !!--------------------------------------------------------------------- 1314 1642 !! *** ROUTINE zdf_vertical_average *** … … 1322 1650 1323 1651 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1652 INTEGER, DIMENSION(jpi,jpj) :: jp_ext 1324 1653 1325 1654 ! Alan: do we need zb? 1326 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity1655 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity 1327 1656 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1328 1657 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1329 1658 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1330 1659 1331 INTEGER :: jk, ji, jj 1660 INTEGER :: jk, ji, jj, ibld_ext 1332 1661 REAL(wp) :: zthick, zthermal, zbeta 1333 1662 … … 1358 1687 zu(ji,jj) = zu(ji,jj) / zthick 1359 1688 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) 1689 zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj) 1690 ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj) 1691 IF ( ibld_ext < mbkt(ji,jj) ) THEN 1692 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem) 1693 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal) 1694 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) & 1695 & / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) ) 1696 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) & 1697 & / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) ) 1698 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1699 ELSE 1700 zdt(ji,jj) = 0._wp 1701 zds(ji,jj) = 0._wp 1702 zdu(ji,jj) = 0._wp 1703 zdv(ji,jj) = 0._wp 1704 zdb(ji,jj) = 0._wp 1705 ENDIF 1368 1706 END DO 1369 1707 END DO … … 1399 1737 END SUBROUTINE zdf_osm_velocity_rotation 1400 1738 1401 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1739 SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 1740 !!--------------------------------------------------------------------- 1741 !! *** ROUTINE zdf_osm_osbl_state_fk *** 1742 !! 1743 !! ** 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. 1744 !! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined 1745 !! lflux :: determines whether effects of surface flux extend below the base of the OSBL 1746 !! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl. 1747 !! 1748 !! ** Method : 1749 !! 1750 !! 1751 !!---------------------------------------------------------------------- 1752 1753 ! Outputs 1754 LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle 1755 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk 1756 ! 1757 REAL(wp), DIMENSION(jpi,jpj) :: znd_param 1758 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer 1759 REAL(wp) :: zpe_mle_ref, zwb_ent, zdbdz_mle_int 1760 1761 znd_param(:,:) = 0._wp 1762 1763 DO jj = 2, jpjm1 1764 DO ji = 2, jpim1 1765 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1766 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 1767 END DO 1768 END DO 1769 DO jj = 2, jpjm1 1770 DO ji = 2, jpim1 1771 ! 1772 IF ( lconv(ji,jj) ) THEN 1773 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1774 zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1775 zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1776 zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1777 zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) ) 1778 ! Calculate potential energies of actual profile and reference profile. 1779 zpe_mle_layer = 0._wp 1780 zpe_mle_ref = 0._wp 1781 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1782 zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) ) 1783 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 1784 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) 1785 END DO 1786 ! Non-dimensional parameter to diagnose the presence of thermocline 1787 1788 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) ) 1789 ENDIF 1790 ENDIF 1791 END DO 1792 END DO 1793 1794 ! Diagnosis 1795 DO jj = 2, jpjm1 1796 DO ji = 2, jpim1 1797 IF ( lconv(ji,jj) ) THEN 1798 zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 1799 & - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 1800 & + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 1801 & - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1802 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 1803 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1804 ! MLE layer growing 1805 IF ( znd_param (ji,jj) > 100. ) THEN 1806 ! Thermocline present 1807 lflux(ji,jj) = .FALSE. 1808 lmle(ji,jj) =.FALSE. 1809 ELSE 1810 ! Thermocline not present 1811 lflux(ji,jj) = .TRUE. 1812 lmle(ji,jj) = .TRUE. 1813 ENDIF ! znd_param > 100 1814 ! 1815 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1816 lpyc(ji,jj) = .FALSE. 1817 ELSE 1818 lpyc = .TRUE. 1819 ENDIF 1820 ELSE 1821 ! MLE layer restricted to OSBL or just below. 1822 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 1823 ! Weak stratification MLE layer can grow. 1824 lpyc(ji,jj) = .FALSE. 1825 lflux(ji,jj) = .TRUE. 1826 lmle(ji,jj) = .TRUE. 1827 ELSE 1828 ! Strong stratification 1829 lpyc(ji,jj) = .TRUE. 1830 lflux(ji,jj) = .FALSE. 1831 lmle(ji,jj) = .FALSE. 1832 ENDIF ! zdb_bl < rn_mle_thresh_bl and 1833 ENDIF ! zhmle > 1.2 zhbl 1834 ELSE 1835 lpyc(ji,jj) = .TRUE. 1836 lflux(ji,jj) = .FALSE. 1837 lmle(ji,jj) = .FALSE. 1838 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 1839 ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 1840 ELSE 1841 ! Stable Boundary Layer 1842 lpyc(ji,jj) = .FALSE. 1843 lflux(ji,jj) = .FALSE. 1844 lmle(ji,jj) = .FALSE. 1845 ENDIF ! lconv 1846 END DO 1847 END DO 1848 END SUBROUTINE zdf_osm_osbl_state_fk 1849 1850 SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 1402 1851 !!--------------------------------------------------------------------- 1403 1852 !! *** ROUTINE zdf_osm_external_gradients *** … … 1409 1858 !!---------------------------------------------------------------------- 1410 1859 1860 INTEGER, DIMENSION(jpi,jpj) :: jbase 1411 1861 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1412 1862 … … 1417 1867 DO jj = 2, jpjm1 1418 1868 DO ji = 2, jpim1 1419 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN1869 IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1420 1870 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1421 1871 zbeta = rab_n(ji,jj,1,jp_sal) 1422 jkb = ibld(ji,jj) + ibld_ext1872 jkb = jbase(ji,jj) 1423 1873 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1424 1874 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & … … 1427 1877 & / e3t_n(ji,jj,ibld(ji,jj)) 1428 1878 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1879 ELSE 1880 zdtdz(ji,jj) = 0._wp 1881 zdsdz(ji,jj) = 0._wp 1882 zdbdz(ji,jj) = 0._wp 1429 1883 END IF 1430 1884 END DO … … 1432 1886 END SUBROUTINE zdf_osm_external_gradients 1433 1887 1434 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )1888 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha ) 1435 1889 1436 1890 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1891 REAL(wp), DIMENSION(jpi,jpj) :: zalpha 1437 1892 1438 1893 INTEGER :: jk, jj, ji 1439 1894 REAL(wp) :: ztgrad, zsgrad, zbgrad 1440 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1441 REAL(wp) :: zzeta_s=0.3, ztmp 1895 REAL(wp) :: zgamma_b_nd, znd 1896 REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc 1897 REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 1442 1898 1443 1899 DO jj = 2, jpjm1 1444 1900 DO ji = 2, jpim1 1445 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN1901 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1446 1902 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 1903 IF ( lpyc(ji,jj) ) THEN 1904 zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 1905 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 ) ) 1906 zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp ) 1907 1908 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1909 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1910 ! Commented lines in this section are not needed in new code, once tested ! 1911 ! can be removed ! 1912 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1913 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 1914 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 1915 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1916 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1917 DO jk = 2, ibld(ji,jj)+ibld_ext 1918 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1919 IF ( znd <= zzeta_m ) THEN 1920 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 1921 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1922 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 1923 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1924 zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * & 1925 & EXP( -6.0 * ( znd -zzeta_m )**2 ) 1926 ELSE 1927 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1928 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1929 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 1930 ENDIF 1931 END DO 1932 ENDIF ! if no pycnocline pycnocline gradients set to zero 1473 1933 ELSE 1474 1934 ! stable conditions 1475 1935 ! if pycnocline profile only defined when depth steady of increasing. 1476 IF ( zdhdt(ji,jj) > =0.0 ) THEN ! Depth increasing, or steady.1936 IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady. 1477 1937 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1478 1938 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline … … 1502 1962 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1503 1963 ENDIF ! IF (lconv) 1504 END IF ! IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) )1964 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) ) 1505 1965 END DO 1506 1966 END DO … … 1527 1987 DO ji = 2, jpim1 1528 1988 ! 1529 IF ( ibld(ji,jj) + ibld_ext< mbkt(ji,jj) ) THEN1989 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1530 1990 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 1991 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 1992 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1993 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1994 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1535 1995 !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 DO1996 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1997 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1998 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1999 ! & )/ (zdh(ji,jj) + epsln ) 2000 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 2001 ! znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 2002 ! IF ( znd <= 0.0 ) THEN 2003 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 2004 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 2005 ! ELSE 2006 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 2007 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 2008 ! ENDIF 2009 ! END DO 1550 2010 ELSE 1551 2011 ! stable conditions … … 1568 2028 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1569 2029 1570 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2)2030 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 1571 2031 !!--------------------------------------------------------------------- 1572 2032 !! *** ROUTINE zdf_osm_calculate_dhdt *** … … 1578 2038 !!---------------------------------------------------------------------- 1579 2039 1580 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zd hdt_2! Rate of change of hbl2040 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt ! Rate of change of hbl 1581 2041 1582 2042 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 2043 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 2044 REAL(wp) :: zvel_max!, zwb_min 1586 2045 REAL(wp) :: zzeta_m = 0.3 1587 2046 REAL(wp) :: zgamma_c = 2.0 1588 2047 REAL(wp) :: zdhoh = 0.1 1589 2048 REAL(wp) :: alpha_bc = 0.5 1590 1591 DO jj = 2, jpjm1 1592 DO ji = 2, jpim1 2049 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 2050 2051 DO jj = 2, jpjm1 2052 DO ji = 2, jpim1 2053 2054 IF ( lshear(ji,jj) ) THEN 1593 2055 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 2056 1618 2057 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 2058 1621 2059 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2060 ! Fox-Kemper buoyancy flux average over OSBL 1622 2061 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1623 2062 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) … … 1628 2067 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1629 2068 ! OSBL is deepening, entrainment > restratification 1630 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 2069 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 2070 ! *** Used for shear Needs to be changed to work stabily 2071 ! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 2072 ! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 2073 ! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 2074 ! za_1 = 1.0 / zgamma_b**2 - 0.017 2075 ! za_2 = 1.0 / zgamma_b**3 - 0.0025 2076 ! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 2077 zpsi = 0._wp 2078 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2079 zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) 2080 IF ( j_ddh(ji,jj) == 1 ) THEN 2081 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 2082 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 ) 2083 ELSE 2084 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 ) 2085 ENDIF 2086 ! Relaxation to dh_ref = zari * hbl 2087 zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2088 2089 ELSE ! j_ddh == 0 2090 ! Growing shear layer 2091 zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2092 ENDIF ! j_ddh 2093 zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 2094 ELSE ! zdb_bl >0 2095 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 2096 ENDIF 2097 ELSE ! zwb_min + 2*zwb_fk_b < 0 2098 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 2099 zdhdt(ji,jj) = - zvel_mle(ji,jj) 2100 2101 2102 ENDIF 2103 2104 ELSE 2105 ! Fox-Kemper not used. 2106 2107 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) / & 2108 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 2109 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2110 ! added ajgn 23 July as temporay fix 2111 2112 ENDIF ! ln_osm_mle 2113 2114 ELSE ! lconv - Stable 2115 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 2116 IF ( zdhdt(ji,jj) < 0._wp ) THEN 2117 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 2118 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 2119 ELSE 2120 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2121 ENDIF 2122 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 2123 ENDIF ! lconv 2124 ELSE ! lshear 2125 IF ( lconv(ji,jj) ) THEN ! Convective 2126 2127 IF ( ln_osm_mle ) THEN 2128 2129 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 2130 ! Fox-Kemper buoyancy flux average over OSBL 2131 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 2132 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 2133 ELSE 2134 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 2135 ENDIF 2136 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2137 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 2138 ! OSBL is deepening, entrainment > restratification 2139 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 1631 2140 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 2141 ELSE … … 1643 2152 ! Fox-Kemper not used. 1644 2153 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) / &2154 zvel_max = -zwb_ent(ji,jj) / & 1646 2155 & MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln) 1647 2156 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1648 2157 ! added ajgn 23 July as temporay fix 1649 2158 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 2159 ENDIF ! ln_osm_mle 2160 1670 2161 ELSE ! Stable 1671 2162 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 2163 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1674 2164 ! 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)2165 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 1676 2166 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) )2167 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1678 2168 ENDIF 1679 2169 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 1680 ENDIF 1681 END DO 1682 END DO 2170 ENDIF ! lconv 2171 ENDIF ! lshear 2172 END DO 2173 END DO 1683 2174 END SUBROUTINE zdf_osm_calculate_dhdt 1684 2175 1685 SUBROUTINE zdf_osm_timestep_hbl( zdhdt , zdhdt_2)2176 SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 1686 2177 !!--------------------------------------------------------------------- 1687 2178 !! *** ROUTINE zdf_osm_timestep_hbl *** … … 1697 2188 1698 2189 1699 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt , zdhdt_2! rates of change of hbl.2190 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl. 1700 2191 1701 2192 INTEGER :: jk, jj, ji, jm … … 1735 2226 IF ( ln_osm_mle ) THEN 1736 2227 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) ), &2228 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1738 2229 & e3w_n(ji,jj,jm) ) 1739 2230 ELSE 1740 2231 zhbl_s = zhbl_s + MIN( & 1741 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj)) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &2232 & rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1742 2233 & e3w_n(ji,jj,jm) ) 1743 2234 ENDIF 1744 2235 1745 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1746 2236 ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2237 IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN 2238 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2239 lpyc(ji,jj) = .FALSE. 2240 ENDIF 1747 2241 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1748 2242 END DO … … 1765 2259 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1766 2260 1767 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2261 ! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2262 IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 2263 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 2264 lpyc(ji,jj) = .FALSE. 2265 ENDIF 1768 2266 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1769 2267 END DO … … 1796 2294 1797 2295 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1798 !2296 ! 1799 2297 INTEGER :: jj, ji 1800 2298 INTEGER :: inhml 1801 REAL(wp) :: zari, ztau, zddhdt 1802 1803 1804 DO jj = 2, jpjm1 1805 DO ji = 2, jpim1 2299 REAL(wp) :: zari, ztau, zdh_ref 2300 REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 2301 2302 DO jj = 2, jpjm1 2303 DO ji = 2, jpim1 2304 2305 IF ( lshear(ji,jj) ) THEN 1806 2306 IF ( lconv(ji,jj) ) THEN 2307 IF ( j_ddh(ji,jj) == 0 ) THEN 2308 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2309 dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt 2310 ELSE 2311 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 2312 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN 2313 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 ) 2314 ELSE 2315 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 ) 2316 ENDIF 2317 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 ) 2318 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2319 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 2320 ENDIF 2321 2322 ELSE ! lconv 2323 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 2324 2325 ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln) 2326 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 2327 ! boundary layer deepening 2328 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2329 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2330 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 2331 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 2332 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 2333 ELSE 2334 zdh_ref = 0.2 * hbl(ji,jj) 2335 ENDIF 2336 ELSE ! IF(dhdt < 0) 2337 zdh_ref = 0.2 * hbl(ji,jj) 2338 ENDIF ! IF (dhdt >= 0) 2339 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2340 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 2341 ! Alan: this hml is never defined or used -- do we need it? 2342 ENDIF 2343 2344 ELSE ! lshear 2345 ! for lshear = .FALSE. calculate ddhdt here 2346 2347 IF ( lconv(ji,jj) ) THEN 1807 2348 1808 2349 IF( ln_osm_mle ) THEN 1809 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN2350 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1810 2351 ! 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)THEN2352 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 1812 2353 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 ) 2354 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 2355 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 ) 2356 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 2357 ENDIF 1819 2358 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)2359 zdh_ref = zari * hbl(ji,jj) 1821 2360 ELSE 1822 2361 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)2362 zdh_ref = 0.2 * hbl(ji,jj) 1824 2363 ENDIF 1825 2364 ELSE 1826 2365 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)2366 zdh_ref = 0.2 * hbl(ji,jj) 1828 2367 ENDIF 1829 2368 ELSE ! ln_osm_mle 1830 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ ext(ji,jj) > 0._wp)THEN2369 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN 1831 2370 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 ) 2371 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 2372 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 ) 2373 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 2374 ENDIF 1838 2375 ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird) 1839 zd dhdt= zari * hbl(ji,jj)2376 zdh_ref = zari * hbl(ji,jj) 1840 2377 ELSE 1841 2378 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)2379 zdh_ref = 0.2 * hbl(ji,jj) 1843 2380 ENDIF 1844 2381 1845 2382 END IF ! ln_osm_mle 1846 2383 1847 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2384 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2385 ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2386 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 1848 2387 ! Alan: this hml is never defined or used 1849 2388 ELSE ! IF (lconv) … … 1855 2394 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1856 2395 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 ) 1857 zd dhdt= MIN( zari, 0.2 ) * hbl(ji,jj)2396 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj) 1858 2397 ELSE 1859 zd dhdt= 0.2 * hbl(ji,jj)2398 zdh_ref = 0.2 * hbl(ji,jj) 1860 2399 ENDIF 1861 2400 ELSE ! IF(dhdt < 0) 1862 zd dhdt= 0.2 * hbl(ji,jj)2401 zdh_ref = 0.2 * hbl(ji,jj) 1863 2402 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? 2403 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2404 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 2405 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 2406 ENDIF ! lshear 2407 2408 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 2409 inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)), 1.e-3) ) , 1 ) 2410 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 2411 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 2412 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 2413 END DO 2414 END DO 1876 2415 1877 2416 END SUBROUTINE zdf_osm_pycnocline_thickness 1878 2417 1879 2418 1880 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )2419 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 1881 2420 !!---------------------------------------------------------------------- 1882 2421 !! *** ROUTINE zdf_osm_horizontal_gradients *** … … 1891 2430 1892 2431 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 2432 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 1893 2433 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 1894 2434 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy … … 1977 2517 END DO 1978 2518 2519 DO jj = 2, jpjm1 2520 DO ji = 2, jpim1 2521 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2522 zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 2523 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 2524 END DO 2525 END DO 2526 1979 2527 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1980 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )2528 SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 1981 2529 !!---------------------------------------------------------------------- 1982 2530 !! *** ROUTINE zdf_osm_mle_parameters *** … … 1989 2537 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1990 2538 1991 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 2539 INTEGER, DIMENSION(jpi,jpj) :: mld_prof 2540 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle 1992 2541 INTEGER :: ji, jj, jk ! dummy loop indices 1993 INTEGER :: ii, ij, ik ! local integers2542 INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers 1994 2543 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 2544 REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle 2545 2546 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 2547 2548 DO jj = 2, jpjm1 2549 DO ji = 2, jpim1 2550 IF ( lconv(ji,jj) ) THEN 2551 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2552 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 2553 zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 2554 zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2 2555 ENDIF 2003 2556 END DO 2004 2557 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 2558 ! Timestep mixed layer eddy depth. 2011 2559 DO jj = 2, jpjm1 2012 2560 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)) 2561 IF ( lmle(ji,jj) ) THEN ! MLE layer growing. 2562 ! Buoyancy gradient at base of MLE layer. 2563 zthermal = rab_n(ji,jj,1,jp_tem) 2564 zbeta = rab_n(ji,jj,1,jp_sal) 2565 jkb = mld_prof(ji,jj) 2566 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 2567 ! 2568 zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) ) 2569 zdb_mle = zb_bl(ji,jj) - zbuoy 2570 ! Timestep hmle. 2571 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / zdb_mle 2572 ELSE 2573 IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN 2574 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 2575 ELSE 2576 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau 2577 ENDIF 2578 ENDIF 2579 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 2580 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 2031 2581 END DO 2032 2582 END DO … … 2045 2595 END DO 2046 2596 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 2597 END SUBROUTINE zdf_osm_mle_parameters 2067 2598 … … 2088 2619 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 2089 2620 & ,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_shelter2621 & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter 2091 2622 ! Namelist for Fox-Kemper parametrization. 2092 2623 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& … … 2138 2669 WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 2139 2670 WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm 2671 WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s' 2140 2672 WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix 2141 2673 WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty … … 2144 2676 WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv 2145 2677 ENDIF 2678 2679 2680 ! ! Check wave coupling settings ! 2681 ! ! Further work needed - see ticket #2447 ! 2682 IF( nn_osm_wave == 2 ) THEN 2683 IF (.NOT. ( ln_wave .AND. ln_sdw )) & 2684 & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 2685 END IF 2146 2686 2147 2687 ! ! allocate zdfosm arrays … … 2171 2711 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2172 2712 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'2713 WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2174 2714 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2175 2715 WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
Note: See TracChangeset
for help on using the changeset viewer.