Changeset 14405
- Timestamp:
- 2021-02-05T13:46:18+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
r14404 r14405 288 288 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 289 289 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 290 REAL(wp), DIMENSION(jpi,jpj) :: zddhdt ! correction to dhdt due to internal structure.291 290 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients 292 291 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients … … 395 394 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 396 395 397 zddhdt(:,:) = 0._wp 398 ! hbl = MAX(hbl,epsln) 396 ! hbl = MAX(hbl,epsln) 399 397 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 400 398 ! Calculate boundary layer scales … … 675 673 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) 676 674 ! Rate of change of hbl 677 CALL zdf_osm_calculate_dhdt( zdhdt , zddhdt)675 CALL zdf_osm_calculate_dhdt( zdhdt ) 678 676 DO jj = 2, jpjm1 679 677 DO ji = 2, jpim1 … … 856 854 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 855 END DO 856 ! 857 IF ( dh(ji,jj) < 0.2*hbl(ji,jj) ) THEN 858 zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 859 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 ) ) 858 860 ! 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 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) 861 862 ! 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 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) 863 864 ! 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 zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 866 DO jk = 2, ibld(ji,jj) 867 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 868 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 865 869 ! 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 870 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 871 END DO 872 END IF 873 873 ENDIF ! End of pycnocline 874 874 ELSE ! lconv test - stable conditions … … 912 912 DO jj = 2, jpjm1 913 913 DO ji = 2, jpim1 914 IF ( lpyc(ji,jj) ) THEN 915 IF ( j_ddh(ji,jj) == 0 ) THEN 914 IF( lconv(ji,jj) ) THEN 915 IF ( lpyc(ji,jj) ) THEN 916 IF ( j_ddh(ji,jj) == 0 ) THEN 916 917 ! Place holding code. Parametrization needs checking for these conditions. 917 918 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird … … 940 941 END DO 941 942 ENDIF ! lpyc 942 END DO ! ji loop 943 ENDIF ! lconv 944 END DO ! ji loop 943 945 END DO ! jj loop 944 946 … … 1029 1031 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 1030 1032 END DO 1033 1031 1034 ELSE 1032 1035 DO jk = 2, ibld(ji,jj) … … 1052 1055 END DO 1053 1056 1054 IF(ln_dia_osm) THEN 1057 ! DO jj = 1, jpjm1 1058 ! DO ji = 1, jpim1 1059 ! IF ( lconv(ji,jj) ) THEN 1060 ! IF ( lpyc(ji,jj) ) THEN 1061 ! zd_cubic = ( 0.948 - 2.13 * zdh(ji,jj) / zhml(ji,jj) ) * zustar(ji,jj)**2 1062 ! zc_cubic = -0.474 * zustar(ji,jj)**2 - zd_cubic 1063 ! DO jk = imld(ji,jj), ibld(ji,jj) 1064 ! zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 1065 ! ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1066 ! END DO 1067 ! zc_cubic= 3.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj) 1068 ! zd_cubic = -2.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj) 1069 ! DO jk = imld(ji,jj), ibld(ji,jj) 1070 ! zznd_pyc = -( gdepw_n(ji,jj,jk)-zhbl(ji,jj) ) / zdh(ji,jj) 1071 ! ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) 1072 ! END DO 1073 ! ENDIF 1074 ! ENDIF 1075 ! END DO 1076 ! END DO 1077 1078 IF(ln_dSia_osm) THEN 1055 1079 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 1056 1080 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) … … 1388 1412 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 1389 1413 1390 IF ( lshear(ji,jj) . and. j_ddh(ji,jj) == 1) THEN1414 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 1391 1415 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 1392 1416 ENDIF … … 1398 1422 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) 1399 1423 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 1400 IF ( lshear(ji,jj) . and. j_ddh(ji,jj) == 1) THEN1424 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 1401 1425 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 1402 1426 ENDIF … … 1508 1532 1509 1533 REAL(wp), DIMENSION(jpi,jpj) :: zekman 1510 REAL(wp) :: zri_p, zri_b ! Richardson numbers1534 REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers 1511 1535 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1512 1536 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1513 1537 1514 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 1515 REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 1516 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 1538 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8 1539 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03 1517 1540 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 1518 1541 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 1542 REAL, PARAMETER :: zri_c = 0.25 1543 REAL, PARAMETER :: zek = 4.0 1519 1544 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. 1520 1545 … … 1530 1555 END DO 1531 1556 1532 zekman(:,:) = EXP( - 4.0* ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )1557 zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1533 1558 1534 1559 WHERE ( lconv ) … … 1543 1568 IF ( lconv(ji,jj) ) THEN 1544 1569 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1545 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 &1570 zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 1546 1571 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1547 1572 1548 zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8)1573 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 ) 1549 1574 1550 1575 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 ) ) 1576 ! Stability Dependence 1577 zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) ) 1551 1578 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1552 1579 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! 1553 1580 ! full code available ! 1554 1581 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1555 IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 1582 IF ( zshear(ji,jj) > 1.e-10 ) THEN 1583 IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN 1556 1584 ! Growing shear layer 1557 1585 j_ddh(ji,jj) = 0 … … 1559 1587 ELSE 1560 1588 j_ddh(ji,jj) = 1 1561 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN1589 ! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 1562 1590 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 1563 lshear(ji,jj) = .TRUE. 1564 ELSE 1591 lshear(ji,jj) = .TRUE. 1592 ! ELSE 1593 ENDIF 1594 ELSE 1595 j_ddh(ji,jj) = 2 1596 lshear(ji,jj) = .FALSE. 1597 ENDIF 1565 1598 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 1566 zshear(ji,jj) = 0.5 * zshear(ji,jj) 1567 lshear(ji,jj) = .FALSE. 1568 ENDIF 1569 ENDIF 1599 ! zshear(ji,jj) = 0.5 * zshear(ji,jj) 1600 ! lshear(ji,jj) = .FALSE. 1601 ! ENDIF 1570 1602 ELSE ! zdb_bl test, note zshear set to zero 1571 1603 j_ddh(ji,jj) = 2 … … 1594 1626 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1595 1627 END IF 1596 zwb_ent(ji,jj) = - 2.0 * 0.2* zrf_conv * zwbav(ji,jj) &1597 & - 0.15* zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &1598 & + zr_stokes * ( 0.15* EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &1599 & - zrf_langmuir * 0.03* zwstrl(ji,jj)**3 ) / zhml(ji,jj)1628 zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) & 1629 & - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1630 & + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1631 & - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1600 1632 ! 1601 1633 ENDIF … … 1610 1642 IF ( lconv(ji,jj) ) THEN 1611 1643 ! Unstable OSBL 1612 zwb_shr = -za_wb_s * z shear(ji,jj)1644 zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj) 1613 1645 IF ( j_ddh(ji,jj) == 0 ) THEN 1614 1646 1615 ! Developing shear layer, additional shear production possible.1616 1617 zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp )1618 zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ))1619 zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )1647 ! ! Developing shear layer, additional shear production possible. 1648 1649 ! zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp ) 1650 ! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) 1651 ! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1620 1652 1621 zwb_shr = -za_wb_s * zshear(ji,jj) 1653 ! zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) 1654 ! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u ) 1622 1655 1623 1656 ENDIF 1624 1657 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1625 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)1658 ! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1626 1659 ELSE ! IF ( lconv ) THEN - ENDIF 1627 1660 ! Stable OSBL - shear production not coded for first attempt. 1628 1661 ENDIF ! lconv 1629 E LSE! lshear1630 1662 ENDIF ! lshear 1663 IF ( lconv(ji,jj) ) THEN 1631 1664 ! Unstable OSBL 1632 zwb_shr = -za_wb_s * zshear(ji,jj) 1633 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1634 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1635 ENDIF ! lconv 1636 ENDIF ! lshear 1665 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1666 ENDIF ! lconv 1637 1667 END DO ! ji 1638 1668 END DO ! jj … … 2033 2063 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 2034 2064 2035 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt , zddhdt)2065 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) 2036 2066 !!--------------------------------------------------------------------- 2037 2067 !! *** ROUTINE zdf_osm_calculate_dhdt *** … … 2043 2073 !!---------------------------------------------------------------------- 2044 2074 2045 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt , zddhdt! Rate of change of hbl2075 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl 2046 2076 2047 2077 INTEGER :: jj, ji 2048 2078 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 2049 REAL(wp) :: zvel_max !, zwb_min2050 REAL(wp) :: zzeta_m = 0.32051 REAL(wp) :: zgamma_c = 2.02052 REAL(wp) :: zdhoh = 0.12053 REAL(wp) :: alpha_bc = 0.52079 REAL(wp) :: zvel_max,zddhdt 2080 REAL(wp), PARAMETER :: zzeta_m = 0.3 2081 REAL(wp), PARAMETER :: zgamma_c = 2.0 2082 REAL(wp), PARAMETER :: zdhoh = 0.1 2083 REAL(wp), PARAMETER :: zalpha_b = 0.3 2054 2084 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 2055 2085 … … 2073 2103 ! OSBL is deepening, entrainment > restratification 2074 2104 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 2075 ! *** Used for shear Needs to be changed to work stabily 2076 ! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 2077 ! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 2078 ! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 2079 ! za_1 = 1.0 / zgamma_b**2 - 0.017 2080 ! za_2 = 1.0 / zgamma_b**3 - 0.0025 2081 ! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 2082 zpsi = 0._wp 2083 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 2084 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) 2105 zgamma_b_nd = zdbdz_bl_ext(ji,jj) + zdh(ji,jj) / zdb_ml(ji,jj) 2106 zpsi = -zalpha_pyc(ji,jj) * ( zwb0(ji,jj) - ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj) 2107 zpsi = zpsi - 4.0 * ( 1.0 + zdh(ji,jj) /zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) 2108 zpsi = zalpha_b * MAX ( zpsi, 0._wp ) 2109 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / MAX( zdb_bl(ji,jj), 1.e-15 ) 2085 2110 IF ( j_ddh(ji,jj) == 1 ) THEN 2086 2111 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN … … 2090 2115 ENDIF 2091 2116 ! Relaxation to dh_ref = zari * hbl 2092 zddhdt (ji,jj)= -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)2117 zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2093 2118 2094 ELSE ! j_ddh == 02119 ELSE IF ( j_ddh(ji,jj) == 0 ) THEN 2095 2120 ! Growing shear layer 2096 zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2121 zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2122 zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 2123 ELSE 2124 zddhdt = 0._wp 2097 2125 ENDIF ! j_ddh 2098 zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj)2126 zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0 * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) * zddhdt / zdb_bl(ji,jj) 2099 2127 ELSE ! zdb_bl >0 2100 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)2128 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 2101 2129 ENDIF 2102 2130 ELSE ! zwb_min + 2*zwb_fk_b < 0 … … 2302 2330 INTEGER :: jj, ji 2303 2331 INTEGER :: inhml 2304 REAL(wp) :: zari, ztau, zdh_ref 2305 REAL, PARAMETER :: a_ddh _2 = 3.5 ! also in pycnocline_depth2332 REAL(wp) :: zari, ztau, zdh_ref, zddhdt 2333 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 2306 2334 2307 2335 DO jj = 2, jpjm1 … … 2310 2338 IF ( lshear(ji,jj) ) THEN 2311 2339 IF ( lconv(ji,jj) ) THEN 2340 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 2312 2341 IF ( j_ddh(ji,jj) == 0 ) THEN 2313 2342 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2314 dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt 2343 ! zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 2344 ! zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 2345 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer 2346 ! dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) ) 2347 ztau = MAX( zdb_bl(ji,jj) * ( 0.625 * hbl(ji,jj) ) / ( a_ddh * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt ) 2348 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.625 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2349 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8 * zhbl(ji,jj) 2315 2350 ELSE 2316 2351 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt … … 2324 2359 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 2325 2360 ENDIF 2326 2361 ELSE 2362 ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0 * rn_rdt ) 2363 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) ) 2364 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(ji,jj) 2365 ENDIF 2327 2366 ELSE ! lconv 2328 2367 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL
Note: See TracChangeset
for help on using the changeset viewer.
