Changeset 14408
- Timestamp:
- 2021-02-05T13:46:32+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90
r14407 r14408 308 308 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 309 309 ! Flux-gradient relationship variables 310 REAL(wp), DIMENSION(jpi, jpj) :: zshear , zri_i ! Shear production and interfacial richardon number.310 REAL(wp), DIMENSION(jpi, jpj) :: zshear ! Shear production. 311 311 312 312 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale. … … 611 611 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 612 612 ! Determine the state of the OSBL, stable/unstable, shear/no shear 613 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear , zri_i)613 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 614 614 615 615 IF ( ln_osm_mle ) THEN … … 766 766 DO jk = 2, ibld(ji,jj) 767 767 zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj) 768 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9* zznd_d ) &768 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & 769 769 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 770 770 ! 771 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9* zznd_d ) &771 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) & 772 772 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 773 773 END DO … … 817 817 818 818 WHERE ( lconv ) 819 zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )820 zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )819 zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 820 zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 821 821 ELSEWHERE 822 822 zsc_wth_1 = 0._wp … … 829 829 DO jk = 2, imld(ji,jj) 830 830 zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 831 ! calculate turbulent lengthscale832 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml -zznd_ml**3 / 3.0 ) ) ) &833 & * ( 1.0 - EXP ( -15.0 * ( 1. 1- zznd_ml ) ) )834 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml -zznd_ml**3 / 3.0 ) ) ) &835 & * ( 1.0 - EXP ( - 5.0 * ( 1.0- zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )831 ! calculate turbulent time scale 832 zl_c = 0.9 * ( 1.0 - EXP ( - 5.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & 833 & * ( 1.0 - EXP ( -15.0 * ( 1.2 - zznd_ml ) ) ) 834 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) & 835 & * ( 1.0 - EXP ( - 8.0 * ( 1.15 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 836 836 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 837 837 ! non-gradient buoyancy terms 838 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0. 5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) * zznd_ml/ ( 0.15 + zznd_ml )839 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0. 5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) * zznd_ml/ ( 0.15 + zznd_ml )838 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.4 * zsc_wth_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) 839 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 * zsc_ws_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml ) 840 840 END DO 841 841 … … 846 846 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) 847 847 ! Cubic profile used for buoyancy term 848 za_cubic = 0.755 * ztau_sc_u(ji,jj)849 zb_cubic = 0.25 * ztau_sc_u(ji,jj)850 848 DO jk = 2, ibld(ji,jj) 851 849 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) … … 1129 1127 DO jj=2, jpjm1 1130 1128 DO ji = 2, jpim1 1131 ghamt(ji,jj,ibld(ji,jj) +ibld_ext) = 0._wp1132 ghams(ji,jj,ibld(ji,jj) +ibld_ext) = 0._wp1133 ghamu(ji,jj,ibld(ji,jj) +ibld_ext) = 0._wp1134 ghamv(ji,jj,ibld(ji,jj) +ibld_ext) = 0._wp1129 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1130 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1131 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1132 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 1135 1133 END DO ! ji loop 1136 1134 END DO ! jj loop … … 1435 1433 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) ) ) 1436 1434 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 1437 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_ s_sc(ji,jj) )1435 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_n_sc(ji,jj) ) 1438 1436 1439 1437 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 … … 1515 1513 END SUBROUTINE zdf_osm_diffusivity_viscosity 1516 1514 1517 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear , zri_i)1515 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 1518 1516 1519 1517 !!--------------------------------------------------------------------- … … 1532 1530 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 1533 1531 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline 1534 REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number1535 1532 1536 1533 ! Local Variables … … 1564 1561 zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1565 1562 1566 WHERE ( lconv )1567 zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 )1568 END WHERE1569 1570 1563 zshear(:,:) = 0._wp 1571 1564 j_ddh(:,:) = 1 … … 1578 1571 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1579 1572 1580 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 ) 1581 1573 IF ( ff_t(ji,jj) >= 0._wp ) THEN 1574 ! Northern Hemisphere 1575 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 ) 1576 ELSE 1577 ! Southern 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 ENDIF 1582 1580 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 ) ) 1583 1581 ! Stability Dependence … … 1796 1794 REAL(wp), DIMENSION(jpi,jpj) :: znd_param 1797 1795 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer 1798 REAL(wp) :: zpe_mle_ref, z wb_ent, zdbdz_mle_int1796 REAL(wp) :: zpe_mle_ref, zdbdz_mle_int 1799 1797 1800 1798 znd_param(:,:) = 0._wp … … 1838 1836 DO ji = 2, jpim1 1839 1837 IF ( lconv(ji,jj) ) THEN 1840 zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 1841 & - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 1842 & + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 1843 & - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1844 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 1838 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN 1845 1839 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1846 1840 ! MLE layer growing … … 1915 1909 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1916 1910 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1917 & / e3t_n(ji,jj, ibld(ji,jj))1911 & / e3t_n(ji,jj,jkb) 1918 1912 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1919 & / e3t_n(ji,jj, ibld(ji,jj))1913 & / e3t_n(ji,jj,jkb) 1920 1914 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1921 1915 ELSE … … 1957 1951 zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1958 1952 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1959 DO jk = 2, ibld(ji,jj) +ibld_ext1953 DO jk = 2, ibld(ji,jj) 1960 1954 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1961 1955 IF ( znd <= zzeta_m ) THEN … … 2158 2152 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 2159 2153 ELSE 2160 zpert = MAX( 2.0 *zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )2154 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2161 2155 ENDIF 2162 2156 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) … … 2205 2199 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 2206 2200 ELSE 2207 zpert = MAX( 2.0 *zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )2201 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 2208 2202 ENDIF 2209 2203 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
Note: See TracChangeset
for help on using the changeset viewer.
