Changeset 14316
- Timestamp:
- 2021-01-19T19:57:12+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14305 r14316 277 277 INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer 278 278 279 REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients280 REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear281 282 279 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 283 280 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid … … 299 296 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 297 ! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline 301 REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline302 REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline303 REAL(wp) :: zdtdz_pyc ! Parametrized gradient of temperature in pycnocline304 REAL(wp) :: zdsdz_pyc ! Parametrised gradient of salinity in pycnocline305 298 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline 306 REAL(wp) :: zdudz_pyc ! u-shear across the pycnocline307 REAL(wp) :: zdvdz_pyc ! v-shear across the pycnocline308 299 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 309 300 ! Flux-gradient relationship variables 310 301 REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number. 311 302 312 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale.313 314 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.315 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.316 REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/317 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.318 303 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep 319 304 … … 334 319 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 335 320 REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc 336 REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML.337 REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc338 REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max339 321 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 340 322 REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau … … 344 326 REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 345 327 REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 346 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z3ddz_pyc_1, z3ddz_pyc_2 ! Pycnocline gradient/shear profiles347 INTEGER :: istat ! Memory allocation status348 328 349 329 ! For debugging … … 354 334 ibld(:,:) = 0 ; imld(:,:) = 0 355 335 zustar(:,:) = 0.0_wp 356 zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp336 zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp 357 337 zwth0(:,:) = 0.0_wp ; zws0(:,:) = 0.0_wp ; zwb0(:,:) = 0.0_wp 358 338 zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp … … 372 352 zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp 373 353 zdb_ml(:,:) = 0._wp 374 zwth_ent = 0._wp ; zws_ent = 0._wp375 354 ! 376 355 zdbdz_pyc(:,:,:) = 0.0_wp … … 386 365 zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 387 366 388 ! Flux-Gradient arrays. 389 zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp 390 zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp 391 zhbl_t(:,:) = 0._wp ; zdhdt(:,:) = 0._wp 367 zhbl_t(:,:) = 0._wp 392 368 393 369 zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp … … 706 682 CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 707 683 708 ! 709 ! calculate non-gradient components of the flux-gradient relationships 710 ! 711 ! Stokes term in scalar flux, flux-gradient relationship 712 WHERE ( lconv ) 713 zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln) 714 ! 715 zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 716 ELSEWHERE 717 zsc_wth_1 = 2.0 * zwthav 718 ! 719 zsc_ws_1 = 2.0 * zwsav 720 ENDWHERE 721 722 723 DO_2D( 0, 0, 0, 0 ) 724 IF ( lconv(ji,jj) ) THEN 725 DO jk = 2, imld(ji,jj) 726 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 727 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 728 ! 729 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 730 END DO ! end jk loop 731 ELSE ! else for if (lconv) 732 ! Stable conditions 733 DO jk = 2, ibld(ji,jj) 734 zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 735 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 736 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 737 ! 738 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 739 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj) 740 END DO 741 ENDIF ! endif for check on lconv 742 743 END_2D 744 745 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 746 WHERE ( lconv ) 747 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 ) 748 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 749 zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 750 ELSEWHERE 751 zsc_uw_1 = zustar**2 752 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 753 ENDWHERE 754 IF(ln_dia_osm) THEN 755 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 756 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 757 END IF 758 DO_2D( 0, 0, 0, 0 ) 759 IF ( lconv(ji,jj) ) THEN 760 DO jk = 2, imld(ji,jj) 761 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 762 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 763 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & 764 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) 765 ! 766 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) & 767 & * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 768 END DO ! end jk loop 769 ELSE 770 ! Stable conditions 771 DO jk = 2, ibld(ji,jj) ! corrected to ibld 772 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 773 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 774 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 775 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 776 END DO ! end jk loop 777 ENDIF 778 END_2D 779 780 ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] 781 782 WHERE ( lconv ) 783 zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 784 zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 785 ELSEWHERE 786 zsc_wth_1 = 0._wp 787 zsc_ws_1 = 0._wp 788 ENDWHERE 789 790 DO_2D( 0, 0, 0, 0 ) 791 IF (lconv(ji,jj) ) THEN 792 DO jk = 2, imld(ji,jj) 793 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 794 ! calculate turbulent length scale 795 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 796 & * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) ) 797 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 798 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 799 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0) 800 ! non-gradient buoyancy terms 801 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 ) 802 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 ) 803 END DO 804 805 IF ( lpyc(ji,jj) ) THEN 806 ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird 807 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 ) 808 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) 809 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) 810 ! Cubic profile used for buoyancy term 811 za_cubic = 0.755 * ztau_sc_u(ji,jj) 812 zb_cubic = 0.25 * ztau_sc_u(ji,jj) 813 DO jk = 2, ibld(ji,jj) 814 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 815 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 ) 816 817 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 ) 818 END DO 819 ! 820 zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj) 821 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 ) ) 822 ! 823 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) 824 ! 825 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) 826 ! 827 zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 828 DO jk = 2, ibld(ji,jj) 829 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 830 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 831 ! 832 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 833 END DO 834 ENDIF ! End of pycnocline 835 ELSE ! lconv test - stable conditions 836 DO jk = 2, ibld(ji,jj) 837 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 838 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 839 END DO 840 ENDIF 841 END_2D 842 843 WHERE ( lconv ) 844 zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln ) 845 zsc_uw_2 = zwb0 * zustke * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0) 846 zsc_vw_1 = 0._wp 847 ELSEWHERE 848 zsc_uw_1 = 0._wp 849 zsc_vw_1 = 0._wp 850 ENDWHERE 851 852 DO_2D( 0, 0, 0, 0 ) 853 IF ( lconv(ji,jj) ) THEN 854 DO jk = 2 , imld(ji,jj) 855 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 856 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 857 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & 858 & * zsc_uw_2(ji,jj) ) 859 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 860 END DO ! jk loop 861 ELSE 862 ! stable conditions 863 DO jk = 2, ibld(ji,jj) 864 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 865 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 866 END DO 867 ENDIF 868 END_2D 869 870 DO_2D( 0, 0, 0, 0 ) 871 IF ( lpyc(ji,jj) ) THEN 872 IF ( j_ddh(ji,jj) == 0 ) THEN 873 ! Place holding code. Parametrization needs checking for these conditions. 874 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 875 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 876 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 877 ELSE 878 zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 879 zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 880 zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 881 ENDIF 882 zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse 883 zc_cubic = zuw_bse - zd_cubic 884 ! need ztau_sc_u to be available. Change to array. 885 DO jk = imld(ji,jj), ibld(ji,jj) 886 zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 887 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) 888 END DO 889 zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 890 zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse 891 zc_cubic = zvw_bse - zd_cubic 892 DO jk = imld(ji,jj), ibld(ji,jj) 893 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj) 894 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) 895 END DO 896 ENDIF ! lpyc 897 END_2D 898 899 IF(ln_dia_osm) THEN 900 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 901 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 902 END IF 903 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 904 905 DO_2D( 1, 0, 1, 0 ) 906 907 IF ( lconv(ji,jj) ) THEN 908 zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) ) 909 zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) ) 910 IF ( lpyc(ji,jj) ) THEN 911 ! Pycnocline scales 912 zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj) 913 zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj) 914 ENDIF 915 ELSE 916 zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj) 917 zsc_ws_1(ji,jj) = zws0(ji,jj) 918 ENDIF 919 END_2D 920 921 DO_2D( 0, 0, 0, 0 ) 922 IF ( lconv(ji,jj) ) THEN 923 DO jk = 2, imld(ji,jj) 924 zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 925 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) & 926 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 927 & - EXP( - 6.0 * zznd_ml ) ) ) & 928 & * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) ) 929 ! 930 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) & 931 & * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) & 932 & - EXP( - 6.0 * zznd_ml ) ) ) & 933 & * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) ) 934 END DO 935 ! 936 IF ( lpyc(ji,jj) ) THEN 937 ! pycnocline 938 DO jk = imld(ji,jj), ibld(ji,jj) 939 zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj) 940 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 ) ) 941 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 ) ) 942 END DO 943 ENDIF 944 ELSE 945 IF( zdhdt(ji,jj) > 0. ) THEN 946 DO jk = 2, ibld(ji,jj) 947 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 948 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 949 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 950 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 951 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 952 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 953 END DO 954 ENDIF 955 ENDIF 956 END_2D 957 958 WHERE ( lconv ) 959 zsc_uw_1 = zustar**2 960 zsc_vw_1 = ff_t * zustke * zhml 961 ELSEWHERE 962 zsc_uw_1 = zustar**2 963 zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1 964 zsc_vw_1 = ff_t * zustke * zhbl 965 zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1 966 ENDWHERE 967 968 DO_2D( 0, 0, 0, 0 ) 969 IF ( lconv(ji,jj) ) THEN 970 DO jk = 2, imld(ji,jj) 971 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 972 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 973 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 974 & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 975 ! 976 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 977 & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 978 END DO 979 ELSE 980 DO jk = 2, ibld(ji,jj) 981 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 982 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 983 IF ( zznd_d <= 2.0 ) THEN 984 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 985 &* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 986 ! 987 ELSE 988 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 989 & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 990 ! 991 ENDIF 992 993 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 994 & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 995 ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 996 & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 997 END DO 998 ENDIF 999 END_2D 1000 1001 IF(ln_dia_osm) THEN 1002 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 1003 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 1004 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 1005 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 1006 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 1007 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 1008 END IF 1009 ! 1010 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 1011 1012 1013 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 1014 1015 DO_2D( 0, 0, 0, 0 ) 1016 IF ( .not. lconv(ji,jj) ) THEN 1017 DO jk = 2, ibld(ji,jj) 1018 znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 1019 IF ( znd >= 0.0 ) THEN 1020 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1021 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 1022 ELSE 1023 ghamu(ji,jj,jk) = 0._wp 1024 ghamv(ji,jj,jk) = 0._wp 1025 ENDIF 1026 END DO 1027 ENDIF 1028 END_2D 1029 1030 ! Pynocline contributions 1031 ! 1032 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Allocate arrays used for output of pycnocline gradient/shear profiles 1033 ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) 1034 CALL mpp_sum( 'zdfosm', istat ) 1035 IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 1036 z3ddz_pyc_1(:,:,:) = 0.0_wp 1037 z3ddz_pyc_2(:,:,:) = 0.0_wp 1038 END IF 1039 DO_2D( 0, 0, 0, 0 ) 1040 IF ( lconv (ji,jj) ) THEN 1041 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 1042 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1043 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1044 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1045 !Alan is this right? 1046 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1047 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1048 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1049 ! & )/ (zdh(ji,jj) + epsln ) 1050 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1051 ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1052 ! IF ( znd <= 0.0 ) THEN 1053 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1054 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1055 ! ELSE 1056 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1057 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1058 ! ENDIF 1059 ! END DO 1060 ELSE ! Stable conditions 1061 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1062 ! If pycnocline profile only defined when depth steady of increasing. 1063 IF ( zdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 1064 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 1065 IF ( zhol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 1066 ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln ) 1067 ztgrad = zdt_bl(ji,jj) * ztmp 1068 zsgrad = zds_bl(ji,jj) * ztmp 1069 zbgrad = zdb_bl(ji,jj) * ztmp 1070 DO jk = 2, ibld(ji,jj) 1071 znd = gdepw(ji,jj,jk,Kmm) * ztmp 1072 zdtdz_pyc = ztgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 1073 zdsdz_pyc = zsgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 ) 1074 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 1075 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 1076 IF ( ln_dia_pyc_scl ) THEN 1077 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 1078 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 1079 END IF 1080 END DO 1081 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1082 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln ) 1083 ztgrad = zdt_bl(ji,jj) * ztmp 1084 zsgrad = zds_bl(ji,jj) * ztmp 1085 zbgrad = zdb_bl(ji,jj) * ztmp 1086 DO jk = 2, ibld(ji,jj) 1087 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp 1088 zdtdz_pyc = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1089 zdsdz_pyc = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1090 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc 1091 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc 1092 IF ( ln_dia_pyc_scl ) THEN 1093 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 1094 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 1095 END IF 1096 END DO 1097 ENDIF ! IF (zhol >=0.5) 1098 ENDIF ! IF (zdb_bl> 0.) 1099 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1100 END IF 1101 END IF 1102 END_2D 1103 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 1104 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 1105 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 1106 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * zdbdz_pyc(:,:,:) ) 1107 END IF 1108 DO_2D( 0, 0, 0, 0 ) 1109 IF ( .NOT. lconv (ji,jj) ) THEN 1110 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1111 zugrad = 3.25_wp * zdu_bl(ji,jj) / zhbl(ji,jj) 1112 zvgrad = 2.75_wp * zdv_bl(ji,jj) / zhbl(ji,jj) 1113 DO jk = 2, ibld(ji,jj) 1114 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1115 IF ( znd < 1.0 ) THEN 1116 zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) 1117 ELSE 1118 zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) 1119 ENDIF 1120 zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) 1121 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc 1122 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc 1123 IF ( ln_dia_pyc_shr ) THEN 1124 z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc 1125 z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc 1126 END IF 1127 END DO 1128 END IF 1129 END IF 1130 END_2D 1131 IF ( ln_dia_pyc_shr ) THEN ! Output of pycnocline shear profiles 1132 IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 1133 IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 1134 END IF 1135 IF(ln_dia_osm) THEN 1136 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 1137 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 1138 END IF 1139 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Deallocate arrays used for output of pycnocline gradient/shear profiles 1140 DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) 1141 END IF 1142 1143 DO_2D( 0, 0, 0, 0 ) 1144 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1145 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1146 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1147 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1148 END_2D 1149 1150 IF(ln_dia_osm) THEN 1151 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 1152 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 1153 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 1154 END IF 684 ! 685 ! Calculate non-gradient components of the flux-gradient relationships 686 ! -------------------------------------------------------------------- 687 CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, ibld_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear, & 688 & zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla, & 689 & zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, & 690 & zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 1155 691 1156 692 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 1329 865 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1330 866 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1331 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux1332 867 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1333 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux1334 868 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1335 869 … … 1374 908 ! 1375 909 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 910 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! Coefficients in cubic polynomial specifying diffusivity in pycnocline 1376 911 1377 912 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 … … 1883 1418 END_2D 1884 1419 ! 1420 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 1421 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) ) 1422 END IF 1423 ! 1885 1424 IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 1886 1425 ! … … 2542 2081 END SUBROUTINE zdf_osm_vertical_average 2543 2082 2083 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, kbld_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear, & 2084 & pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla, & 2085 & pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml, & 2086 & pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) 2087 !!--------------------------------------------------------------------- 2088 !! *** ROUTINE zdf_osm_fgr_terms *** 2089 !! 2090 !! ** Purpose : Compute non-gradient terms in flux-gradient relationship 2091 !! 2092 !! ** Method : 2093 !! 2094 !!---------------------------------------------------------------------- 2095 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2096 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 2097 INTEGER, DIMENSION(:,:), INTENT(in ) :: kmld ! ML base layer 2098 INTEGER, DIMENSION(:,:), INTENT(in ) :: kp_ext ! Offset for external level 2099 INTEGER, INTENT(in ) :: kbld_ext ! Offset for external level 2100 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 2101 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldpyc ! Pycnocline flags 2102 INTEGER, DIMENSION(:,:), INTENT(in ) :: k_ddh ! Type of shear layer 2103 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 2104 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phml ! ML depth 2105 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline depth 2106 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdhdt ! BL depth tendency 2107 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phol ! Stability parameter for boundary layer 2108 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pshear ! Shear production 2109 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustar ! Friction velocity 2110 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrl ! Langmuir velocity scale 2111 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pvstr ! Velocity scale (approaches zustar for large Langmuir number) 2112 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwstrc ! Convective velocity scale 2113 REAL(wp), DIMENSION(:,:), INTENT(in ) :: puw0 ! Surface u-momentum flux 2114 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwth0 ! Surface heat flux 2115 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pws0 ! Surface freshwater flux 2116 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb0 ! Surface buoyancy flux 2117 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwthav ! BL average heat flux 2118 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwsav ! BL average freshwater flux 2119 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwbav ! BL average buoyancy flux 2120 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pustke ! Surface Stokes drift 2121 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pla ! Langmuir number 2122 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdt_bl ! Temperature diff. between BL average and basal value 2123 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pds_bl ! Salinity diff. between BL average and basal value 2124 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value 2125 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_bl ! Velocity diff. (u) between BL average and basal value 2126 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_bl ! Velocity diff. (u) between BL average and basal value 2127 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdt_ml ! Temperature diff. between mixed-layer average and basal value 2128 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pds_ml ! Salinity diff. between mixed-layer average and basal value 2129 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_ml ! Buoyancy diff. between mixed-layer average and basal value 2130 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_ml ! Velocity diff. (u) between mixed-layer average and basal value 2131 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_ml ! Velocity diff. (v) between mixed-layer average and basal value 2132 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdtdz_bl_ext ! External temperature gradients 2133 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdsdz_bl_ext ! External salinity gradients 2134 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 2135 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdbdz_pyc ! Pycnocline buoyancy gradients 2136 REAL(wp), DIMENSION(:,:), INTENT(in ) :: palpha_pyc ! 2137 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdiffut ! t-diffusivity 2138 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pviscos ! Viscosity 2139 ! 2140 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z3ddz_pyc_1, z3ddz_pyc_2 ! Pycnocline gradient/shear profiles 2141 ! 2142 INTEGER :: ji, jj, jk, jkm_bld, jkf_mld, jkm_mld ! Loop indices 2143 INTEGER :: istat ! Memory allocation status 2144 REAL(wp) :: zznd_d, zznd_ml, zznd_pyc, znd ! Temporary non-dimensional depths 2145 REAL(wp), DIMENSION(A2D(0)) :: zsc_wth_1,zsc_ws_1 ! Temporary scales 2146 REAL(wp), DIMENSION(A2D(0)) :: zsc_uw_1, zsc_uw_2 ! Temporary scales 2147 REAL(wp), DIMENSION(A2D(0)) :: zsc_vw_1, zsc_vw_2 ! Temporary scales 2148 REAL(wp), DIMENSION(A2D(0)) :: ztau_sc_u ! Dissipation timescale at base of WML 2149 REAL(wp) :: zbuoy_pyc_sc, zdelta_pyc ! 2150 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale 2151 REAL(wp), DIMENSION(A2D(0)) :: za_cubic, zb_cubic ! Coefficients in cubic polynomial specifying 2152 REAL(wp), DIMENSION(A2D(0)) :: zc_cubic, zd_cubic ! diffusivity in pycnocline 2153 REAL(wp), DIMENSION(A2D(0)) :: zwt_pyc_sc_1, zws_pyc_sc_1 ! 2154 REAL(wp), DIMENSION(A2D(0)) :: zzeta_pyc ! 2155 REAL(wp) :: zomega, zvw_max ! 2156 REAL(wp) :: zuw_bse,zvw_bse ! Momentum, heat, and salinity fluxes 2157 REAL(wp), DIMENSION(A2D(0)) :: zwth_ent,zws_ent ! at the top of the pycnocline 2158 REAL(wp), DIMENSION(A2D(0)) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term 2159 REAL(wp) :: ztmp ! 2160 REAL(wp) :: ztgrad, zsgrad, zbgrad ! Variables used to calculate pycnocline gradients 2161 REAL(wp) :: zugrad, zvgrad ! Variables for calculating pycnocline shear 2162 REAL(wp) :: zdtdz_pyc ! Parametrized gradient of temperature in pycnocline 2163 REAL(wp) :: zdsdz_pyc ! Parametrised gradient of salinity in pycnocline 2164 REAL(wp) :: zdudz_pyc ! u-shear across the pycnocline 2165 REAL(wp) :: zdvdz_pyc ! v-shear across the pycnocline 2166 !!---------------------------------------------------------------------- 2167 ! 2168 IF( ln_timing ) CALL timing_start('zdf_osm_ft') 2169 ! 2170 ! Auxiliary indices 2171 ! ----------------- 2172 jkm_bld = 0 2173 jkf_mld = jpk 2174 jkm_mld = 0 2175 DO_2D( 0, 0, 0, 0 ) 2176 IF ( kbld(ji,jj) > jkm_bld ) jkm_bld = kbld(ji,jj) 2177 IF ( kbld(ji,jj) < jkf_mld ) jkf_mld = kbld(ji,jj) 2178 IF ( kmld(ji,jj) > jkm_mld ) jkm_mld = kmld(ji,jj) 2179 END_2D 2180 ! 2181 ! Stokes term in scalar flux, flux-gradient relationship 2182 ! ------------------------------------------------------ 2183 WHERE ( ldconv(A2D(0)) ) 2184 zsc_wth_1(:,:) = pwstrl(A2D(0))**3 * pwth0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2185 zsc_ws_1(:,:) = pwstrl(A2D(0))**3 * pws0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2186 ELSEWHERE 2187 zsc_wth_1(:,:) = 2.0_wp * pwthav(A2D(0)) 2188 zsc_ws_1(:,:) = 2.0_wp * pwsav(A2D(0)) 2189 ENDWHERE 2190 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 2191 IF ( ldconv(ji,jj) ) THEN 2192 IF ( jk <= kmld(ji,jj) ) THEN 2193 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2194 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) * & 2195 & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 2196 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) * & 2197 & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 2198 END IF 2199 ELSE ! Stable conditions 2200 IF ( jk <= kbld(ji,jj) ) THEN 2201 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2202 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) * & 2203 & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 2204 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) * & 2205 & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 2206 END IF 2207 END IF ! Check on ldconv 2208 END_3D 2209 ! 2210 IF ( ln_dia_osm ) THEN 2211 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 2212 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 2213 END IF 2214 ! 2215 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use 2216 ! zvstr since term needs to go to zero as zwstrl goes to zero) 2217 ! --------------------------------------------------------------------- 2218 WHERE ( ldconv(A2D(0)) ) 2219 zsc_uw_1(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / & 2220 & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * pla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 2221 zsc_uw_2(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) / & 2222 & MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 2223 zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2224 & ( ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 2225 ELSEWHERE 2226 zsc_uw_1(:,:) = pustar(A2D(0))**2 2227 zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2228 & ( pvstr(A2D(0))**2 + epsln ) 2229 ENDWHERE 2230 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 2231 IF ( ldconv(ji,jj) ) THEN 2232 IF ( jk <= kmld(ji,jj) ) THEN 2233 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2234 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) + & 2235 & 0.00125_wp * EXP( -1.0_wp * zznd_d ) * zsc_uw_2(ji,jj) ) * & 2236 & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) 2237 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65_wp * 0.15_wp * EXP( -1.0_wp * zznd_d ) * & 2238 & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_vw_1(ji,jj) 2239 END IF 2240 ELSE ! Stable conditions 2241 IF ( jk <= kbld(ji,jj) ) THEN ! Corrected to ibld 2242 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2243 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) * & 2244 & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_uw_1(ji,jj) 2245 END IF 2246 END IF 2247 END_3D 2248 ! 2249 ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio 2250 ! (X0.3) and pressure (X0.5)] 2251 ! ---------------------------------------------------------------------- 2252 WHERE ( ldconv(A2D(0)) ) 2253 zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) / & 2254 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2255 zsc_ws_1(:,:) = pwbav(A2D(0)) * pws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) / & 2256 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2257 ELSEWHERE 2258 zsc_wth_1(:,:) = 0.0_wp 2259 zsc_ws_1(:,:) = 0.0_wp 2260 ENDWHERE 2261 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 2262 IF ( ldconv(ji,jj) ) THEN 2263 IF ( jk <= kmld(ji,jj) ) THEN 2264 zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 2265 ! Calculate turbulent length scale 2266 zl_c = 0.9_wp * ( 1.0_wp - EXP( -7.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) * & 2267 & ( 1.0_wp - EXP( -15.0_wp * ( 1.1_wp - zznd_ml ) ) ) 2268 zl_l = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) * & 2269 & ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 2270 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) 2271 ! Non-gradient buoyancy terms 2272 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.5_wp * zsc_wth_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 2273 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.5_wp * zsc_ws_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml ) 2274 END IF 2275 ELSE ! Stable conditions 2276 IF ( jk <= kbld(ji,jj) ) THEN 2277 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 2278 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj) 2279 END IF 2280 END IF 2281 END_3D 2282 DO_2D( 0, 0, 0, 0 ) 2283 IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 2284 ztau_sc_u(ji,jj) = phml(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & 2285 & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**1.5_wp ) 2286 zwth_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & 2287 & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj) 2288 zws_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & 2289 & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 2290 zbuoy_pyc_sc = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 2291 zdelta_pyc = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird / & 2292 & SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 2293 zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) * & 2294 & zdelta_pyc**2 / pdh(ji,jj) 2295 zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) * & 2296 & zdelta_pyc**2 / pdh(ji,jj) 2297 zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 2298 END IF 2299 END_2D 2300 DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 2301 IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN 2302 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 2303 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - & 2304 & 0.045_wp * ( ( zwth_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * & 2305 & MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 2306 ghams(ji,jj,jk) = ghams(ji,jj,jk) - & 2307 & 0.045_wp * ( ( zws_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * & 2308 & MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 2309 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & 2310 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2311 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2312 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp * zws_pyc_sc_1(ji,jj) * & 2313 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2314 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2315 END IF ! End of pycnocline 2316 END_3D 2317 ! 2318 IF(ln_dia_osm) THEN 2319 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! Upward turb. temperature entrainment flux 2320 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! Upward turb. salinity entrainment flux 2321 END IF 2322 ! 2323 zsc_vw_1(:,:) = 0.0_wp 2324 WHERE ( ldconv(A2D(0)) ) 2325 zsc_uw_1(:,:) = -1.0_wp * pwb0(A2D(0)) * pustar(A2D(0))**2 * phml(A2D(0)) / & 2326 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2327 zsc_uw_2(:,:) = pwb0(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) / & 2328 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp ) 2329 ELSEWHERE 2330 zsc_uw_1(:,:) = 0.0_wp 2331 ENDWHERE 2332 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 2333 IF ( ldconv(ji,jj) ) THEN 2334 IF ( jk <= kmld(ji,jj) ) THEN 2335 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2336 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp * & 2337 & ( zsc_uw_1(ji,jj) + 0.125_wp * EXP( -0.5_wp * zznd_d ) * & 2338 & ( 1.0_wp - EXP( -0.5_wp * zznd_d ) ) * zsc_uw_2(ji,jj) ) 2339 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 2340 END IF 2341 ELSE ! Stable conditions 2342 IF ( jk <= kbld(ji,jj) ) THEN 2343 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 2344 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 2345 END IF 2346 ENDIF 2347 END_3D 2348 ! 2349 DO_2D( 0, 0, 0, 0 ) 2350 IF ( ldpyc(ji,jj) ) THEN 2351 IF ( k_ddh(ji,jj) == 0 ) THEN 2352 ! Place holding code. Parametrization needs checking for these conditions. 2353 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird 2354 zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 2355 zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 2356 ELSE 2357 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird 2358 zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 2359 zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 2360 ENDIF 2361 zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse 2362 za_cubic(ji,jj) = zuw_bse - zb_cubic(ji,jj) 2363 zvw_max = 0.7_wp * ff_t(ji,jj) * ( pustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * pustar(ji,jj) * phml(ji,jj) ) 2364 zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse 2365 zc_cubic(ji,jj) = zvw_bse - zd_cubic(ji,jj) 2366 END IF 2367 END_2D 2368 DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld ) ! Need ztau_sc_u to be available. Change to array. 2369 IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 2370 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 2371 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 * & 2372 & ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) * & 2373 & ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 2374 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 * & 2375 & ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) * & 2376 & ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 2377 END IF ! ldpyc 2378 END_3D 2379 ! 2380 IF(ln_dia_osm) THEN 2381 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 2382 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 2383 END IF 2384 ! 2385 ! Transport term in flux-gradient relationship [note : includes ROI ratio 2386 ! (X0.3) ] 2387 ! ----------------------------------------------------------------------- 2388 WHERE ( ldconv(A2D(0)) ) 2389 zsc_wth_1(:,:) = pwth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) 2390 zsc_ws_1(:,:) = pws0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) ) 2391 WHERE ( ldpyc(A2D(0)) ) ! Pycnocline scales 2392 zsc_wth_pyc(:,:) = -0.2_wp * pwb0(A2D(0)) * pdt_bl(A2D(0)) / pdb_bl(A2D(0)) 2393 zsc_ws_pyc(:,:) = -0.2_wp * pwb0(A2D(0)) * pds_bl(A2D(0)) / pdb_bl(A2D(0)) 2394 END WHERE 2395 ELSEWHERE 2396 zsc_wth_1(:,:) = 2.0 * pwthav(A2D(0)) 2397 zsc_ws_1(:,:) = pws0(A2D(0)) 2398 END WHERE 2399 DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) 2400 IF ( ldconv(ji,jj) ) THEN 2401 IF ( ( jk > 1 ) .AND. ( jk <= kmld(ji,jj) ) ) THEN 2402 zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 2403 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) * & 2404 & ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) 2405 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) * & 2406 & ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) ) 2407 END IF 2408 IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN ! Pycnocline 2409 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 2410 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) 2411 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) ) 2412 END IF 2413 ELSE 2414 IF( pdhdt(ji,jj) > 0. ) THEN 2415 IF ( ( jk > 1 ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 2416 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2417 znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 2418 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) + & 2419 7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_wth_1(ji,jj) 2420 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) + & 2421 7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_ws_1(ji,jj) 2422 END IF 2423 ENDIF 2424 ENDIF 2425 END_3D 2426 ! 2427 WHERE ( ldconv(A2D(0)) ) 2428 zsc_uw_1(:,:) = pustar(A2D(0))**2 2429 zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phml(A2D(0)) 2430 ELSEWHERE 2431 zsc_uw_1(:,:) = pustar(A2D(0))**2 2432 zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) * & 2433 & zsc_uw_1(:,:) 2434 zsc_vw_1(:,:) = ff_t * pustke(A2D(0)) * phbl(A2D(0)) 2435 zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) * & 2436 & zsc_vw_1(:,:) 2437 ENDWHERE 2438 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 2439 IF ( ldconv(ji,jj) ) THEN 2440 IF ( jk <= kmld(ji,jj) ) THEN 2441 zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 2442 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2443 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + & 2444 & 0.3_wp * ( -2.0_wp + 2.5_wp * ( 1.0_wp + 0.1_wp * zznd_ml**4 ) - EXP( -8.0_wp * zznd_ml ) ) * & 2445 & zsc_uw_1(ji,jj) 2446 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + & 2447 & 0.3_wp * 0.1_wp * ( EXP( -1.0_wp * zznd_d ) + EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) * & 2448 & zsc_vw_1(ji,jj) 2449 END IF 2450 ELSE 2451 IF ( jk <= kbld(ji,jj) ) THEN 2452 znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 2453 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2454 IF ( zznd_d <= 2.0 ) THEN 2455 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp * & 2456 & ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) * & 2457 & ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ) * zsc_uw_1(ji,jj) 2458 ELSE 2459 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp * & 2460 & ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - znd ) ) ) * zsc_uw_2(ji,jj) 2461 ENDIF 2462 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) * & 2463 & EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj) 2464 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 2465 END IF 2466 END IF 2467 END_3D 2468 ! 2469 IF(ln_dia_osm) THEN 2470 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask *ghamu ) 2471 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask *ghamv ) 2472 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 2473 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 2474 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 2475 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 2476 END IF 2477 ! 2478 ! Make surface forced velocity non-gradient terms go to zero at the base 2479 ! of the mixed layer. 2480 ! 2481 ! Make surface forced velocity non-gradient terms go to zero at the base 2482 ! of the boundary layer. 2483 DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 2484 IF ( ( .NOT. ldconv(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 2485 znd = ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj) ! ALMG to think about 2486 IF ( znd >= 0.0_wp ) THEN 2487 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) 2488 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) ) 2489 ELSE 2490 ghamu(ji,jj,jk) = 0.0_wp 2491 ghamv(ji,jj,jk) = 0.0_wp 2492 ENDIF 2493 END IF 2494 END_3D 2495 ! 2496 ! Pynocline contributions 2497 ! 2498 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Allocate arrays for output of pycnocline gradient/shear profiles 2499 ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat ) 2500 IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 2501 z3ddz_pyc_1(:,:,:) = 0.0_wp 2502 z3ddz_pyc_2(:,:,:) = 0.0_wp 2503 END IF 2504 DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 2505 IF ( ldconv (ji,jj) ) THEN 2506 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 2507 ! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 2508 ! & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 2509 ! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 2510 !Alan is this right? 2511 ! zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 2512 ! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 2513 ! & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 2514 ! & )/ (zdh(ji,jj) + epsln ) 2515 ! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 2516 ! znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 2517 ! IF ( znd <= 0.0 ) THEN 2518 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 2519 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 2520 ! ELSE 2521 ! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 2522 ! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 2523 ! ENDIF 2524 ! END DO 2525 ELSE ! Stable conditions 2526 IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 2527 ! Pycnocline profile only defined when depth steady of increasing. 2528 IF ( pdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 2529 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 2530 IF ( phol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 2531 ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 2532 ztgrad = pdt_bl(ji,jj) * ztmp 2533 zsgrad = pds_bl(ji,jj) * ztmp 2534 zbgrad = pdb_bl(ji,jj) * ztmp 2535 IF ( jk <= kbld(ji,jj) ) THEN 2536 znd = gdepw(ji,jj,jk,Kmm) * ztmp 2537 zdtdz_pyc = ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 2538 zdsdz_pyc = zsgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 2539 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc 2540 ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc 2541 IF ( ln_dia_pyc_scl ) THEN 2542 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 2543 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 2544 END IF 2545 END IF 2546 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 2547 ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 2548 ztgrad = pdt_bl(ji,jj) * ztmp 2549 zsgrad = pds_bl(ji,jj) * ztmp 2550 zbgrad = pdb_bl(ji,jj) * ztmp 2551 IF ( jk <= kbld(ji,jj) ) THEN 2552 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 2553 zdtdz_pyc = ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 2554 zdsdz_pyc = zsgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 2555 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc 2556 ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc 2557 IF ( ln_dia_pyc_scl ) THEN 2558 z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc 2559 z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc 2560 END IF 2561 END IF 2562 ENDIF ! IF (zhol >=0.5) 2563 ENDIF ! IF (zdb_bl> 0.) 2564 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 2565 END IF 2566 END IF 2567 END_3D 2568 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 2569 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 2570 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 2571 END IF 2572 DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 2573 IF ( .NOT. ldconv (ji,jj) ) THEN 2574 IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 2575 zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj) 2576 zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj) 2577 IF ( jk <= kbld(ji,jj) ) THEN 2578 znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 2579 IF ( znd < 1.0 ) THEN 2580 zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 ) 2581 ELSE 2582 zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 ) 2583 ENDIF 2584 zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 ) 2585 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + pviscos(ji,jj,jk) * zdudz_pyc 2586 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + pviscos(ji,jj,jk) * zdvdz_pyc 2587 IF ( ln_dia_pyc_shr ) THEN 2588 z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc 2589 z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc 2590 END IF 2591 END IF 2592 END IF 2593 END IF 2594 END_3D 2595 IF ( ln_dia_pyc_shr ) THEN ! Output of pycnocline shear profiles 2596 IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) ) 2597 IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 2598 END IF 2599 IF(ln_dia_osm) THEN 2600 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 2601 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 2602 END IF 2603 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Deallocate arrays used for output of pycnocline gradient/shear profiles 2604 DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 ) 2605 END IF 2606 ! 2607 DO_2D( 0, 0, 0, 0 ) 2608 ghamt(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 2609 ghams(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 2610 ghamu(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 2611 ghamv(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp 2612 END_2D 2613 ! 2614 IF(ln_dia_osm) THEN 2615 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 2616 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 2617 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*pviscos ) 2618 END IF 2619 ! 2620 IF( ln_timing ) CALL timing_stop('zdf_osm_ft') 2621 ! 2622 END SUBROUTINE zdf_osm_fgr_terms 2623 2544 2624 SUBROUTINE zdf_osm_init( Kmm ) 2545 2625 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.