- Timestamp:
- 2021-04-20T15:46:13+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining
- Files:
-
- 2 edited
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (1 diff)
-
src/OCE/ZDF/zdfosm.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/SHARED/field_def_nemo-oce.xml
r14571 r14734 314 314 <field id="zds_ml" long_name="salinity jump at base of ML" unit="10^-3" /> 315 315 <field id="zdb_ml" long_name="buoyancy jump at base of ML" unit="m/s^2" /> 316 <field id="pb_coup" long_name="bottom coupling velocity" unit="m/s" /> 316 317 <!-- extra OSMOSIS diagnostics for debugging --> 317 318 <field id="zsc_uw_1_0" long_name="zsc u-momentum flux on T after Stokes" unit="m^2/s^2" /> -
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14732 r14734 58 58 USE eosbn2 ! equation of state 59 59 USE traqsr ! details of solar radiation absorption 60 USE zdfdrg, ONLY : rCdU_bot ! bottom friction velocity 60 61 USE zdfddm ! double diffusion mixing (avs array) 61 62 USE iom ! I/O library … … 272 273 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 273 274 LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers 275 LOGICAL, DIMENSION(jpi,jpj) :: lcoup ! Coupling to bottom 274 276 LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present 275 277 LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer. … … 580 582 zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 581 583 ELSE 584 zwstrc(ji,jj) = 0.0_wp 582 585 zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln ) 583 586 ENDIF 584 587 #ifdef key_osm_debug 585 588 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 586 WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/, 2(a,g11.3),/)') &589 WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,3(a,g11.3),/)') & 587 590 & 'After reduction: zustke=', zustke(ji,jj), ' dstokes=', dstokes(ji,jj), & 588 591 & ' zustar =', zustar(ji,jj), ' zwstrl=', zwstrl(ji,jj), ' zwstrc=', zwstrc(ji,jj),& 589 & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj) 592 & ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj), ' zvstr=', zvstr(ji,jj) 590 593 FLUSH(narea+100) 591 594 END IF … … 607 610 DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 608 611 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 609 ibld(ji,jj) = MIN(mbkt(ji,jj) , jk)612 ibld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 610 613 ENDIF 611 614 END_3D … … 641 644 642 645 ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere 643 jp_ext(:,:) = 2646 jp_ext(:,:) = 1 ! ag 19/03 644 647 CALL zdf_osm_vertical_average( Kbb, Kmm, & 645 648 & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & 646 649 & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 647 ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1648 CALL zdf_osm_vertical_average( Kbb, Kmm, 649 & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, ibld-imld+1, &650 jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 651 CALL zdf_osm_vertical_average( Kbb, Kmm, & 652 & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, jp_ext, & 650 653 & zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 651 654 #ifdef key_osm_debug … … 747 750 748 751 !! External gradient below BL needed both with and w/o FK 749 CALL zdf_osm_external_gradients( ibld+ 2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )752 CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03 750 753 751 754 ! Test if pycnocline well resolved 752 DO_2D( 0, 0, 0, 0 ) 753 IF (lconv(ji,jj) ) THEN 754 ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm)755 IF ( ztmp > 6 ) THEN756 ! pycnocline well resolved757 jp_ext(ji,jj) = 1758 ELSE759 ! pycnocline poorly resolved760 jp_ext(ji,jj) = 0761 ENDIF762 ELSE763 ! Stable conditions764 jp_ext(ji,jj) = 0765 ENDIF766 END_2D755 ! DO_2D( 0, 0, 0, 0 ) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 756 ! IF (lconv(ji,jj) ) THEN should account for this. 757 ! ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) 758 ! IF ( ztmp > 6 ) THEN 759 ! ! pycnocline well resolved 760 ! jp_ext(ji,jj) = 1 761 ! ELSE 762 ! ! pycnocline poorly resolved 763 ! jp_ext(ji,jj) = 0 764 ! ENDIF 765 ! ELSE 766 ! ! Stable conditions 767 ! jp_ext(ji,jj) = 0 768 ! ENDIF 769 ! END_2D 767 770 #ifdef key_osm_debug 768 771 IF(narea==nn_narea_db) THEN … … 777 780 778 781 ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here.. 782 jp_ext(:,:) = 1 ! ag 19/03 779 783 CALL zdf_osm_vertical_average( Kbb, Kmm, & 780 784 & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & 781 785 & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 782 ! jp_ext = ibld-imld+1786 jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 783 787 CALL zdf_osm_vertical_average( Kbb, Kmm, & 784 788 & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, & 785 & ibld-imld+1, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )789 & jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) ! ag 19/03 786 790 #ifdef key_osm_debug 787 791 IF(narea==nn_narea_db) THEN … … 802 806 ! Rate of change of hbl 803 807 CALL zdf_osm_calculate_dhdt( zdhdt ) 808 ! Test if surface boundary layer coupled to bottom 809 lcoup(:,:) = .FALSE. ! ag 19/03 804 810 DO_2D( 0, 0, 0, 0 ) 805 811 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it 806 812 ! adjustment to represent limiting by ocean bottom 807 IF ( zhbl_t(ji,jj) >= gdepw(ji, jj, mbkt(ji,jj) + 1, Kmm ) ) THEN 808 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm) - depth_tol)! ht(:,:)) 809 lpyc(ji,jj) = .FALSE. 810 ENDIF 813 IF ( mbkt(ji,jj) > 2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access 814 IF ( zhbl_t(ji,jj) > gdepw(ji, jj,mbkt(ji,jj)-2,Kmm) ) THEN 815 zhbl_t(ji,jj) = MIN( zhbl_t(ji,jj), gdepw(ji,jj,mbkt(ji,jj)-2,Kmm) ) ! ht(:,:)) 816 lpyc(ji,jj) = .FALSE. 817 lcoup(ji,jj) = .TRUE. ! ag 19/03 818 END IF 819 END IF 811 820 #ifdef key_osm_debug 812 821 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 813 WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3))')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),& 814 & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_Dt,' delta hbl from w ', ww(ji,jj,ibld(ji,jj))*rn_Dt 822 WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3)),2(a,l7)')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),& 823 & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_Dt,' delta hbl from w ', ww(ji,jj,ibld(ji,jj))*rn_Dt, & 824 & ' lcoup= ', lcoup(ji,jj), ' lpyc= ', lpyc(ji,jj) 815 825 FLUSH(narea+100) 816 826 END IF … … 834 844 835 845 ! Recalculate BL averages and differences using new BL depth 846 jp_ext(:,:) = 1 ! ag 19/03 836 847 CALL zdf_osm_vertical_average( Kbb, Kmm, & 837 848 & ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, & 838 849 & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 839 !840 !841 ! Check to see if lpyc needs to be changed842 850 843 851 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 844 852 853 ! Reset l_pyc before calculating terms in the flux-gradient relationship 854 845 855 DO_2D( 0, 0, 0, 0 ) 846 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE. 856 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) >= mbkt(ji,jj) - 2 .or. & 857 & ibld(ji,jj)-imld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0.0_wp ) THEN ! ag 19/03 858 lpyc(ji,jj) = .FALSE. ! ag 19/03 859 IF ( ibld(ji,jj) >= mbkt(ji,jj) -2 ) THEN 860 imld(ji,jj) = ibld(ji,jj) - 1 ! ag 19/03 861 zdh(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) - gdepw(ji,jj,imld(ji,jj),Kmm) ! ag 19/03 862 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) ! ag 19/03 863 dh(ji,jj) = zdh(ji,jj) ! ag 19/03 864 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) ! ag 19/03 865 #ifdef key_osm_debug 866 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 867 WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= F5: ibld(ji,jj) >= mbkt(ji,jj) -2' 868 WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj) 869 WRITE(narea+100,'(2(a,g11.3))')'dh=',dh(ji,jj),' hml=',hml(ji,jj) 870 FLUSH(narea+100) 871 END IF 872 #endif 873 ENDIF 874 ENDIF ! ag 19/03 847 875 END_2D 848 876 … … 852 880 ! jp_ext = ibld - imld +1 853 881 ! Recalculate ML averages and differences using new ML depth 882 jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03 854 883 CALL zdf_osm_vertical_average( Kbb, Kmm, & 855 884 & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, & 856 & ibld-imld+1, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )857 ! rotate mean currents and changes onto wind align co-ordinates 858 !885 & jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 886 887 CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 859 888 #ifdef key_osm_debug 860 889 IF(narea==nn_narea_db) THEN … … 871 900 END IF 872 901 #endif 902 903 ! rotate mean currents and changes onto wind align co-ordinates 904 873 905 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 874 906 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) … … 887 919 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 888 920 889 CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )921 jp_ext(:,:) = 1 ! ag 19/03 890 922 CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_pyc ) 891 923 #ifdef key_osm_debug … … 978 1010 zfri = ( 1.0_wp - zfri * zfri ) 979 1011 zrimix = zfri * zfri * zfri * wmask(ji, jj, jk) 980 zdiffut(ji,jj,jk) = zrimix*rn_difri981 zviscos(ji,jj,jk) = zrimix*rn_difri1012 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zrimix*rn_difri ) 1013 zviscos(ji,jj,jk) = MAX( zviscos(ji,jj,jk), zrimix*rn_difri ) 982 1014 END IF 983 1015 END_2D … … 989 1021 DO_2D( 0, 0, 0, 0 ) 990 1022 DO jk = ibld(ji,jj) + 1, jpkm1 991 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv1023 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) 992 1024 END DO 993 1025 END_2D … … 1190 1222 REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr 1191 1223 REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc 1224 REAL(wp), DIMENSION(jpi,jpj) :: zb_coup, zc_coup_vis, zc_coup_dif 1192 1225 ! 1193 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac 1226 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b 1194 1227 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! Coefficients in cubic polynomial specifying diffusivity in pycnocline 1228 REAL(wp) :: zznd_ml, zznd_pyc 1229 REAL(wp) :: zmsku, zmskv 1195 1230 1196 1231 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375 … … 1199 1234 1200 1235 IF( ln_timing ) CALL timing_start('zdf_osm_dv') 1236 1237 zb_coup(:,:) = 0.0_wp 1238 1201 1239 DO_2D( 0, 0, 0, 0 ) 1202 1240 IF ( lconv(ji,jj) ) THEN … … 1257 1295 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) ) 1258 1296 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) ) 1259 #ifdef key_osm_debug1260 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1261 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)1262 FLUSH(narea+100)1263 END IF1264 #endif1265 1297 1266 1298 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 1267 1299 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 1268 1300 ELSE 1269 zbeta_d_sc(ji,jj) = 1.0 1270 zbeta_v_sc(ji,jj) = 1.0 1271 ENDIF 1301 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/03 1302 zdifpyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 1303 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/03 1304 zvispyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 1305 IF(lcoup(ji,jj) ) THEN ! ag 19/03 1306 ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub| 1307 ! already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90 1308 ! Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub) 1309 ! wet-cell averaging .. 1310 zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 1311 zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 1312 zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) * & 1313 & SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 1314 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) ) 1315 1316 zz_b = -gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 1317 zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / zhml(ji,jj) - zb_coup(ji,jj) ) / & 1318 & ( zhml(ji,jj) + zz_b ) ! ag 19/03 1319 #ifdef key_osm_debug 1320 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 1321 WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj), & 1322 & ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj) 1323 WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv 1324 FLUSH(narea+100) 1325 END IF 1326 #endif 1327 !#ifdef key_osm_debug 1328 ! WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj) 1329 ! WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', zz_b, 'pb_coup', zb_coup(ji,jj), & 1330 ! & ' pc_coup_vis', zc_coup_vis(ji,jj) 1331 ! FLUSH(narea+400) 1332 !#endif 1333 zz_b = -zhml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 1334 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / & 1335 & zvisml_sc(ji,jj) ! ag 19/03 1336 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / & 1337 & zdifml_sc(ji,jj) )**p2third 1338 zc_coup_dif(ji,jj) = 0.5_wp * ( -zdifml_sc(ji,jj) / zhml(ji,jj) * ( 1.0_wp - zbeta_d_sc(ji,jj) )**1.5_wp + & 1339 & 1.5_wp * ( zdifml_sc(ji,jj) / zhml(ji,jj) ) * zbeta_d_sc(ji,jj) * & 1340 & SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b ! ag 19/03 1341 #ifdef key_osm_debug 1342 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 1343 WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj) 1344 FLUSH(narea+100) 1345 END IF 1346 #endif 1347 !#ifdef key_osm_debug 1348 ! WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', zc_coup_dif(ji,jj) 1349 ! FLUSH(narea+400) 1350 !#endif 1351 ELSE ! ag 19/03 1352 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) / & 1353 & ( zdifml_sc(ji,jj) + epsln ) )**p2third ! ag 19/03 1354 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / & 1355 & ( zvisml_sc(ji,jj) + epsln ) ! ag 19/03 1356 ENDIF ! ag 19/03 1357 ENDIF ! ag 19/03 1272 1358 #ifdef key_osm_debug 1273 1359 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 1274 1360 WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj) 1361 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 1362 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 1275 1363 FLUSH(narea+100) 1276 1364 END IF … … 1298 1386 & * ( 1.0 - 0.5 * zznd_ml**2 ) 1299 1387 END DO 1388 1389 ! Coupling to bottom 1390 1391 IF ( lcoup(ji,jj) ) THEN ! ag 19/03 1392 DO jk = mbkt(ji,jj), imld(ji,jj), -1 ! ag 19/03 1393 zz_b = - ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) ! ag 19/03 1394 zviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ! ag 19/03 1395 zdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2 ! ag 19/03 1396 END DO ! ag 19/03 1397 ENDIF ! ag 19/03 1300 1398 ! pycnocline 1301 IF ( lpyc(ji,jj) ) THEN 1302 ! Diffusivity profile in the pycnocline given by cubic polynomial. 1399 IF ( lpyc(ji,jj) ) THEN 1400 ! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed. 1303 1401 za_cubic = 0.5 1304 1402 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) … … 1325 1423 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 ) 1326 1424 END DO 1327 IF ( zdhdt(ji,jj) > 0._wp ) THEN1328 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )1329 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )1330 ELSE1331 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp1332 zviscos(ji,jj,ibld(ji,jj)) = 0._wp1333 ENDIF1425 ! IF ( zdhdt(ji,jj) > 0._wp ) THEN 1426 ! zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 1427 ! zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 1428 ! ELSE 1429 ! zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 1430 ! zviscos(ji,jj,ibld(ji,jj)) = 0._wp 1431 ! ENDIF 1334 1432 ENDIF 1335 1433 ELSE … … 1348 1446 ! 1349 1447 END_2D 1448 IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) ) ! BBL-coupling velocity scale 1350 1449 IF( ln_timing ) CALL timing_stop('zdf_osm_dv') 1351 1450 … … 2564 2663 DO_2D( 0, 0, 0, 0 ) 2565 2664 ibld_ext = knlev(ji,jj) + kp_ext(ji,jj) 2566 IF ( ibld_ext < mbkt(ji,jj) ) THEN 2665 IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN ! ag 09/03 2666 ! Two external levels are available 2567 2667 pdt(ji,jj) = pt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm) 2568 2668 pds(ji,jj) = ps(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm) … … 2861 2961 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 2862 2962 #endif 2863 IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN2963 IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) .AND. kbld(ji,jj) - kmld(ji,jj) > 3 ) THEN 2864 2964 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & 2865 2965 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & … … 3574 3674 DO_2D( 1, 1, 1, 1 ) 3575 3675 iiki = MAX(4,imld_rst(ji,jj)) 3576 hbl 3676 hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth 3577 3677 dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth 3678 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 3578 3679 END_2D 3579 3680
Note: See TracChangeset
for help on using the changeset viewer.
