- Timestamp:
- 2021-04-30T13:27:48+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14760 r14775 39 39 !! 'ln_zdfosm' OSMOSIS scheme 40 40 !!---------------------------------------------------------------------- 41 !! zdf_osm : update momentum and tracer Kz from osm scheme 42 !! zdf_osm_vertical_average : compute vertical averages over boundary layers 43 !! zdf_osm_init : initialization, namelist read, and parameters control 44 !! osm_rst : read (or initialize) and write osmosis restart fields 45 !! tra_osm : compute and add to the T & S trend the non-local flux 46 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 47 !! dyn_osm : compute and add to u & v trensd the non-local flux 48 !! 49 !! Subroutines in revised code. 41 !! zdf_osm : update momentum and tracer Kz from osm scheme 42 !! zdf_osm_vertical_average : compute vertical averages over boundary layers 43 !! zdf_osm_velocity_rotation : rotate velocity components 44 !! zdf_osm_velocity_rotation_2d : rotation of 2d fields 45 !! zdf_osm_velocity_rotation_3d : rotation of 3d fields 46 !! zdf_osm_osbl_state : determine the state of the OSBL 47 !! zdf_osm_diffusivity_viscosity : compute eddy diffusivity and viscosity profiles 48 !! zdf_osm_fgr_terms : compute flux-gradient relationship terms 49 !! zdf_osm_init : initialization, namelist read, and parameters control 50 !! osm_rst : read (or initialize) and write osmosis restart fields 51 !! tra_osm : compute and add to the T & S trend the non-local flux 52 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 53 !! dyn_osm : compute and add to u & v trensd the non-local flux 50 54 !!---------------------------------------------------------------------- 51 55 USE oce ! ocean dynamics and active tracers … … 83 87 PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90 84 88 85 ! Scales86 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swth0 ! Surface heat flux (Kinematic)87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sws0 ! Surface freshwater flux88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swb0 ! Surface buoyancy flux89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: suw0 ! Surface u-momentum flux90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustar ! Friction velocity91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: scos_wind ! Cos angle of surface stress92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssin_wind ! Sin angle of surface stress93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swthav ! Heat flux - bl average94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swsav ! Freshwater flux - bl average95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swbav ! Buoyancy flux - bl average96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustke ! Surface Stokes drift97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes ! Penetration depth of the Stokes drift98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrl ! Langmuir velocity scale99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrc ! Convective velocity scale100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sla ! Trubulent Langmuir number101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: svstr ! Velocity scale that tends to sustar for large Langmuir number102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: shol ! Stability parameter for boundary layer103 104 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux 105 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux … … 111 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth 112 97 113 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts114 98 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 115 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML 116 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML 117 101 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer. 102 103 INTERFACE zdf_osm_velocity_rotation 104 !!--------------------------------------------------------------------- 105 !! *** INTERFACE zdf_velocity_rotation *** 106 !!--------------------------------------------------------------------- 107 MODULE PROCEDURE zdf_osm_velocity_rotation_2d 108 MODULE PROCEDURE zdf_osm_velocity_rotation_3d 109 END INTERFACE 110 111 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts 112 113 ! Scales 114 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swth0 ! Surface heat flux (Kinematic) 115 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sws0 ! Surface freshwater flux 116 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swb0 ! Surface buoyancy flux 117 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: suw0 ! Surface u-momentum flux 118 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustar ! Friction velocity 119 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: scos_wind ! Cos angle of surface stress 120 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssin_wind ! Sin angle of surface stress 121 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swthav ! Heat flux - bl average 122 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swsav ! Freshwater flux - bl average 123 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swbav ! Buoyancy flux - bl average 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sustke ! Surface Stokes drift 125 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes ! Penetration depth of the Stokes drift 126 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrl ! Langmuir velocity scale 127 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: swstrc ! Convective velocity scale 128 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sla ! Trubulent Langmuir number 129 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: svstr ! Velocity scale that tends to sustar for large Langmuir number 130 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: shol ! Stability parameter for boundary layer 118 131 119 132 ! !!** Namelist namzdf_osm ** … … 174 187 INTEGER :: idebug = 236 175 188 INTEGER :: jdebug = 228 176 177 INTERFACE zdf_osm_velocity_rotation178 !!---------------------------------------------------------------------179 !! *** INTERFACE zdf_velocity_rotation ***180 !!---------------------------------------------------------------------181 MODULE PROCEDURE zdf_osm_velocity_rotation_2d182 MODULE PROCEDURE zdf_osm_velocity_rotation_3d183 END INTERFACE184 189 185 190 !! * Substitutions … … 287 292 REAL(wp), DIMENSION(jpi,jpj) :: zwb_min 288 293 289 290 294 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 291 295 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux … … 366 370 swstrl(:,:) = zlarge ; svstr(:,:) = zlarge ; swstrc(:,:) = zlarge ; suw0(:,:) = zlarge 367 371 swth0(:,:) = zlarge ; sws0(:,:) = zlarge ; swb0(:,:) = zlarge 368 swthav(:,:) = zlarge ; swsav(:,:) = zlarge ; swbav(:,:) = zlarge ; zwb_ent(:,:) = zlarge372 swthav(:,:) = zlarge ; swsav(:,:) = zlarge ; swbav(:,:) = zlarge 369 373 sustke(:,:) = zlarge ; sla(:,:) = zlarge ; scos_wind(:,:) = zlarge ; ssin_wind(:,:) = zlarge 370 374 shol(:,:) = zlarge ; zalpha_pyc(:,:) = zlarge 371 l conv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE.375 lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE. 372 376 ! mixed layer 373 377 ! no initialization of zhbl or zhml (or zdh?) … … 705 709 706 710 ! Determine the state of the OSBL, stable/unstable, shear/no shear 707 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 711 CALL zdf_osm_osbl_state( Kmm, lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, & 712 & zhbl, zhml, zdh, zdb_bl, zdu_bl, zdv_bl, zdb_ml, zdu_ml, zdv_ml ) 708 713 709 714 #ifdef key_osm_debug … … 965 970 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 966 971 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 967 CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos ) 972 CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, ibld, imld, lconv, lshear, lpyc, lcoup, j_ddh, zdiffut, zviscos, & 973 & zhbl, zhml, zdh, zdhdt, zshear, zb_bl, zu_bl, zv_bl, zb_ml, zu_ml, zv_ml, zwb_ent, & 974 & zwb_min ) 968 975 #ifdef key_osm_debug 969 976 IF(narea==nn_narea_db) THEN … … 1212 1219 1213 1220 CONTAINS 1214 ! subroutine code changed, needs syntax checking.1215 SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )1216 1217 !!---------------------------------------------------------------------1218 !! *** ROUTINE zdf_osm_diffusivity_viscosity ***1219 !!1220 !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.1221 !!1222 !! ** Method :1223 !!1224 !! !!----------------------------------------------------------------------1225 REAL(wp), DIMENSION(:,:,:) :: zdiffut1226 REAL(wp), DIMENSION(:,:,:) :: zviscos1227 ! local1228 1229 ! Scales used to calculate eddy diffusivity and viscosity profiles1230 REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc1231 REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr1232 REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr1233 REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc1234 REAL(wp), DIMENSION(jpi,jpj) :: zb_coup, zc_coup_vis, zc_coup_dif1235 !1236 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b1237 REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! Coefficients in cubic polynomial specifying diffusivity in pycnocline1238 REAL(wp) :: zznd_ml, zznd_pyc1239 REAL(wp) :: zmsku, zmskv1240 1241 REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.3751242 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.1421243 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.151244 1245 IF( ln_timing ) CALL timing_start('zdf_osm_dv')1246 1247 zb_coup(:,:) = 0.0_wp1248 1249 DO_2D( 0, 0, 0, 0 )1250 IF ( lconv(ji,jj) ) THEN1251 1252 zvel_sc_pyc = ( 0.15 * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird1253 zvel_sc_ml = ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird1254 zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-shol(ji,jj) ) ) )**1.25 ) )**21255 1256 zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml1257 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)1258 #ifdef key_osm_debug1259 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1260 WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)1261 WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac1262 FLUSH(narea+100)1263 END IF1264 #endif1265 1266 IF ( lpyc(ji,jj) ) THEN1267 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)1268 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)1269 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac1270 #ifdef key_osm_debug1271 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1272 WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)1273 FLUSH(narea+100)1274 END IF1275 #endif1276 1277 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN1278 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj)1279 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj)1280 ENDIF1281 #ifdef key_osm_debug1282 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1283 WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)1284 FLUSH(narea+100)1285 END IF1286 #endif1287 1288 zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) )1289 zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) )1290 #ifdef key_osm_debug1291 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1292 WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)1293 FLUSH(narea+100)1294 END IF1295 #endif1296 zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac1297 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac1298 #ifdef key_osm_debug1299 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1300 WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)1301 FLUSH(narea+100)1302 END IF1303 #endif1304 1305 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) )1306 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) )1307 1308 zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third1309 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )1310 ELSE1311 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/031312 zdifpyc_s_sc(ji,jj) = 0.0_wp ! ag 19/031313 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) ! ag 19/031314 zvispyc_s_sc(ji,jj) = 0.0_wp ! ag 19/031315 IF(lcoup(ji,jj) ) THEN ! ag 19/031316 ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub|1317 ! already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F901318 ! Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)1319 ! wet-cell averaging ..1320 zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )1321 zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )1322 zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) * &1323 & SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 &1324 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) )1325 1326 zz_b = -gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/031327 zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / zhml(ji,jj) - zb_coup(ji,jj) ) / &1328 & ( zhml(ji,jj) + zz_b ) ! ag 19/031329 #ifdef key_osm_debug1330 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1331 WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj), &1332 & ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj)1333 WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv1334 FLUSH(narea+100)1335 END IF1336 #endif1337 !#ifdef key_osm_debug1338 ! WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj)1339 ! WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', zz_b, 'pb_coup', zb_coup(ji,jj), &1340 ! & ' pc_coup_vis', zc_coup_vis(ji,jj)1341 ! FLUSH(narea+400)1342 !#endif1343 zz_b = -zhml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/031344 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 ) / &1345 & zvisml_sc(ji,jj) ! ag 19/031346 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / &1347 & zdifml_sc(ji,jj) )**p2third1348 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 + &1349 & 1.5_wp * ( zdifml_sc(ji,jj) / zhml(ji,jj) ) * zbeta_d_sc(ji,jj) * &1350 & SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b ! ag 19/031351 #ifdef key_osm_debug1352 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1353 WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj)1354 FLUSH(narea+100)1355 END IF1356 #endif1357 !#ifdef key_osm_debug1358 ! WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', zc_coup_dif(ji,jj)1359 ! FLUSH(narea+400)1360 !#endif1361 ELSE ! ag 19/031362 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) / &1363 & ( zdifml_sc(ji,jj) + epsln ) )**p2third ! ag 19/031364 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / &1365 & ( zvisml_sc(ji,jj) + epsln ) ! ag 19/031366 ENDIF ! ag 19/031367 ENDIF ! ag 19/031368 #ifdef key_osm_debug1369 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1370 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)1371 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)1372 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)1373 FLUSH(narea+100)1374 END IF1375 #endif1376 ELSE1377 zdifml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)1378 zvisml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)1379 #ifdef key_osm_debug1380 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1381 WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)1382 FLUSH(narea+100)1383 END IF1384 #endif1385 END IF1386 END_2D1387 !1388 DO_2D( 0, 0, 0, 0 )1389 IF ( lconv(ji,jj) ) THEN1390 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity1391 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)1392 !1393 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.51394 !1395 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &1396 & * ( 1.0 - 0.5 * zznd_ml**2 )1397 END DO1398 1399 ! Coupling to bottom1400 1401 IF ( lcoup(ji,jj) ) THEN ! ag 19/031402 DO jk = mbkt(ji,jj), imld(ji,jj), -1 ! ag 19/031403 zz_b = - ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) ! ag 19/031404 zviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ! ag 19/031405 zdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2 ! ag 19/031406 END DO ! ag 19/031407 ENDIF ! ag 19/031408 ! pycnocline1409 IF ( lpyc(ji,jj) ) THEN1410 ! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed.1411 za_cubic = 0.51412 zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)1413 zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) &1414 & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8)1415 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )1416 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic1417 DO jk = imld(ji,jj) , ibld(ji,jj)1418 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)1419 !1420 zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )1421 1422 zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 )1423 END DO1424 ! viscosity profiles.1425 za_cubic = 0.51426 zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)1427 zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj) ) / MAX(zvispyc_n_sc(ji,jj), 1.e-8)1428 zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )1429 zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic1430 DO jk = imld(ji,jj) , ibld(ji,jj)1431 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)1432 zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )1433 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 )1434 END DO1435 ! IF ( zdhdt(ji,jj) > 0._wp ) THEN1436 ! zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )1437 ! zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )1438 ! ELSE1439 ! zdiffut(ji,jj,ibld(ji,jj)) = 0._wp1440 ! zviscos(ji,jj,ibld(ji,jj)) = 0._wp1441 ! ENDIF1442 ENDIF1443 ELSE1444 ! stable conditions1445 DO jk = 2, ibld(ji,jj)1446 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)1447 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.51448 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )1449 END DO1450 1451 IF ( zdhdt(ji,jj) > 0._wp ) THEN1452 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)1453 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)1454 ENDIF1455 ENDIF ! end if ( lconv )1456 !1457 END_2D1458 IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) ) ! BBL-coupling velocity scale1459 IF( ln_timing ) CALL timing_stop('zdf_osm_dv')1460 1461 END SUBROUTINE zdf_osm_diffusivity_viscosity1462 1463 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )1464 1465 !!---------------------------------------------------------------------1466 !! *** ROUTINE zdf_osm_osbl_state ***1467 !!1468 !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number1469 !!1470 !! ** Method :1471 !!1472 !! !!----------------------------------------------------------------------1473 1474 INTEGER, DIMENSION(jpi,jpj) :: j_ddh ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low.1475 1476 LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear1477 1478 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.1479 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline1480 1481 ! Local Variables1482 1483 INTEGER :: jj, ji1484 1485 REAL(wp), DIMENSION(jpi,jpj) :: zekman1486 REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers1487 REAL(wp) :: zshear_u, zshear_v, zwb_shr1488 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes1489 1490 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.81491 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.031492 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.151493 REAL, PARAMETER :: rn_ri_p_thresh = 27.01494 REAL, PARAMETER :: zri_c = 0.251495 REAL, PARAMETER :: zek = 4.01496 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress.1497 1498 IF( ln_timing ) CALL timing_start('zdf_osm_os')1499 ! Determins stability and set flag lconv1500 DO_2D( 0, 0, 0, 0 )1501 IF ( shol(ji,jj) < 0._wp ) THEN1502 lconv(ji,jj) = .TRUE.1503 ELSE1504 lconv(ji,jj) = .FALSE.1505 ENDIF1506 END_2D1507 1508 zekman(A2D(0)) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * zhbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) )1509 1510 zshear(:,:) = zlarge1511 zshear(A2D(0)) = 0._wp1512 #ifdef key_osm_debug1513 IF(narea==nn_narea_db) THEN1514 ji=iloc_db; jj=jloc_db1515 WRITE(narea+100,'(a,g11.3)') &1516 & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj)1517 FLUSH(narea+100)1518 END IF1519 #endif1520 j_ddh(A2D(0)) = 11521 1522 DO_2D( 0, 0, 0, 0 )1523 IF ( lconv(ji,jj) ) THEN1524 IF ( zdb_bl(ji,jj) > 0._wp ) THEN1525 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) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1.e-6 ) )**2 &1526 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp )1527 1528 IF ( ff_t(ji,jj) >= 0.0_wp ) THEN1529 ! Northern hemisphere1530 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( -1.0_wp * zdv_ml(ji,jj), 1e-5_wp)**2 )1531 ELSE1532 ! Southern hemisphere1533 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( zdv_ml(ji,jj), 1e-5_wp)**2 )1534 END IF1535 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( sustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )1536 #ifdef key_osm_debug1537 ! IF(narea==nn_narea_db)THEN1538 ! WRITE(narea+100,'(2(a,i10.4))')'ji',ji,'jj',jj1539 ! WRITE(narea+100,'(2(a,i10.4))')'iloc_db',iloc_db,'jloc_db',jloc_db1540 ! WRITE(narea+100,'(2(a,i10.4))')'iloc_db+',mi0(nn_idb),'jloc_db+',mj0(nn_jdb)1541 ! FLUSH(narea+100)1542 ! END IF1543 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1544 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',zshear(ji,jj)1545 WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj)1546 FLUSH(narea+100)1547 END IF1548 #endif1549 ! Stability dependence1550 zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) )1551 #ifdef key_osm_debug1552 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1553 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',zshear(ji,jj)1554 FLUSH(narea+100)1555 END IF1556 #endif1557 1558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1559 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when !1560 ! full code available !1561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1562 IF ( zshear(ji,jj) > 1e-10 ) THEN1563 IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN1564 ! Growing shear layer1565 j_ddh(ji,jj) = 01566 lshear(ji,jj) = .TRUE.1567 ELSE1568 j_ddh(ji,jj) = 11569 ! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN1570 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.1571 lshear(ji,jj) = .TRUE.1572 ! ELSE1573 END IF1574 ELSE1575 j_ddh(ji,jj) = 21576 lshear(ji,jj) = .FALSE.1577 END IF1578 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.1579 ! zshear(ji,jj) = 0.5 * zshear(ji,jj)1580 ! lshear(ji,jj) = .FALSE.1581 ! ENDIF1582 ELSE ! zdb_bl test, note zshear set to zero1583 j_ddh(ji,jj) = 21584 lshear(ji,jj) = .FALSE.1585 ENDIF1586 ENDIF1587 END_2D1588 1589 ! Calculate entrainment buoyancy flux due to surface fluxes.1590 1591 DO_2D( 0, 0, 0, 0 )1592 IF ( lconv(ji,jj) ) THEN1593 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln1594 zrf_conv = TANH( ( swstrc(ji,jj) / zwcor )**0.69 )1595 zrf_shear = TANH( ( sustar(ji,jj) / zwcor )**0.69 )1596 zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69 )1597 IF (nn_osm_SD_reduce > 0 ) THEN1598 ! Effective Stokes drift already reduced from surface value1599 zr_stokes = 1.0_wp1600 ELSE1601 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,1602 ! requires further reduction where BL is deep1603 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &1604 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )1605 END IF1606 zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) &1607 & - zalpha_s * zrf_shear * sustar(ji,jj)**3 / zhml(ji,jj) &1608 & + zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 &1609 & - zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / zhml(ji,jj)1610 !1611 #ifdef key_osm_debug1612 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1613 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',zwb_ent(ji,jj)1614 FLUSH(narea+100)1615 END IF1616 #endif1617 1618 ENDIF1619 END_2D1620 1621 zwb_min(:,:) = zlarge1622 1623 DO_2D( 0, 0, 0, 0 )1624 IF ( lshear(ji,jj) ) THEN1625 IF ( lconv(ji,jj) ) THEN1626 ! Unstable OSBL1627 zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * zshear(ji,jj)1628 #ifdef key_osm_debug1629 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1630 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr1631 FLUSH(narea+100)1632 END IF1633 #endif1634 IF ( j_ddh(ji,jj) == 0 ) THEN1635 1636 ! ! Developing shear layer, additional shear production possible.1637 1638 ! zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp )1639 ! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )1640 ! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )1641 1642 ! 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 )1643 ! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u )1644 #ifdef key_osm_debug1645 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1646 WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, &1647 & ' zshear=',zshear(ji,jj),' zshear_u=', zshear_u1648 FLUSH(narea+100)1649 END IF1650 #endif1651 1652 ENDIF1653 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr1654 ! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)1655 ELSE ! IF ( lconv ) THEN - ENDIF1656 ! Stable OSBL - shear production not coded for first attempt.1657 ENDIF ! lconv1658 END IF ! lshear1659 IF ( lconv(ji,jj) ) THEN1660 ! Unstable OSBL1661 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2.0_wp * swbav(ji,jj)1662 END IF ! lconv1663 #ifdef key_osm_debug1664 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1665 WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), &1666 & ' zwb_min=',zwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' swbav= ', swbav(ji,jj)1667 FLUSH(narea+100)1668 END IF1669 #endif1670 END_2D1671 IF( ln_timing ) CALL timing_stop('zdf_osm_os')1672 END SUBROUTINE zdf_osm_osbl_state1673 1674 1221 1675 1222 SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) … … 2752 2299 ! 2753 2300 END SUBROUTINE zdf_osm_velocity_rotation_3d 2301 2302 SUBROUTINE zdf_osm_osbl_state( Kmm, ldconv, ldshear, k_ddh, pwb_ent, pwb_min, pshear, & 2303 & phbl, phml, pdh, pdb_bl, pdu_bl, pdv_bl, pdb_ml, pdu_ml, pdv_ml ) 2304 !!--------------------------------------------------------------------- 2305 !! *** ROUTINE zdf_osm_osbl_state *** 2306 !! 2307 !! ** Purpose : Determines the state of the OSBL, stable/unstable, 2308 !! shear/ noshear. Also determines shear production, 2309 !! entrainment buoyancy flux and interfacial Richardson 2310 !! number 2311 !! 2312 !! ** Method : 2313 !! 2314 !!---------------------------------------------------------------------- 2315 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 2316 LOGICAL, DIMENSION(:,:), INTENT( out) :: ldconv ! BL stability flags 2317 LOGICAL, DIMENSION(:,:), INTENT( out) :: ldshear ! Shear layers 2318 INTEGER, DIMENSION(:,:), INTENT( out) :: k_ddh ! k_ddh=0: active shear layer 2319 ! ! k_ddh=1: shear layer not active 2320 ! ! k_ddh=2: shear production low 2321 REAL(wp), DIMENSION(:,:), INTENT( out) :: pwb_ent ! Buoyancy fluxes at base 2322 REAL(wp), DIMENSION(:,:), INTENT( out) :: pwb_min ! of well-mixed layer 2323 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear ! Production of TKE due to shear across the pycnocline 2324 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 2325 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phml ! ML depth 2326 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline depth 2327 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value 2328 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_bl ! Velocity diff. (u) between BL average and basal value 2329 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_bl ! Velocity diff. (v) between BL average and basal value 2330 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_ml ! Buoyancy diff. between mixed-layer average and basal value 2331 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdu_ml ! Velocity diff. (u) between mixed-layer average and basal value 2332 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_ml ! Velocity diff. (v) between mixed-layer average and basal value 2333 ! 2334 ! Local Variables 2335 INTEGER :: jj, ji ! Loop indices 2336 ! 2337 REAL(wp), DIMENSION(A2D(0)) :: zekman 2338 REAL(wp), DIMENSION(A2D(0)) :: zri_p, zri_b ! Richardson numbers 2339 REAL(wp) :: zshear_u, zshear_v, zwb_shr 2340 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 2341 ! 2342 REAL(wp), PARAMETER :: za_shr = 0.4_wp, zb_shr = 6.5_wp, za_wb_s = 0.8_wp 2343 REAL(wp), PARAMETER :: zalpha_c = 0.2_wp, zalpha_lc = 0.03_wp 2344 REAL(wp), PARAMETER :: zalpha_ls = 0.06_wp, zalpha_s = 0.15_wp 2345 REAL(wp), PARAMETER :: rn_ri_p_thresh = 27.0_wp 2346 REAL(wp), PARAMETER :: zri_c = 0.25_wp 2347 REAL(wp), PARAMETER :: zek = 4.0_wp 2348 REAL(wp), PARAMETER :: zrot = 0.0_wp ! Dummy rotation rate of surface stress 2349 REAL(wp), PARAMETER :: zlarge = -1e10_wp 2350 ! 2351 IF( ln_timing ) CALL timing_start('zdf_osm_os') 2352 ! 2353 ! Initialise INTENT( out) arrays 2354 ldconv(:,:) = .FALSE. 2355 ldshear(:,:) = .FALSE. 2356 k_ddh(:,:) = 1 2357 pwb_ent(:,:) = zlarge 2358 pwb_min(:,:) = zlarge 2359 pshear(:,:) = zlarge 2360 ! 2361 ! Determins stability and set flag ldconv 2362 DO_2D( 0, 0, 0, 0 ) 2363 IF ( shol(ji,jj) < 0._wp ) THEN 2364 ldconv(ji,jj) = .TRUE. 2365 ELSE 2366 ldconv(ji,jj) = .FALSE. 2367 ENDIF 2368 END_2D 2369 ! 2370 pshear(A2D(0)) = 0._wp 2371 zekman(:,:) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * phbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) ) 2372 ! 2373 #ifdef key_osm_debug 2374 IF(narea==nn_narea_db) THEN 2375 ji=iloc_db; jj=jloc_db 2376 WRITE(narea+100,'(a,g11.3)') & 2377 & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj) 2378 FLUSH(narea+100) 2379 END IF 2380 #endif 2381 ! 2382 DO_2D( 0, 0, 0, 0 ) 2383 IF ( ldconv(ji,jj) ) THEN 2384 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 2385 zri_p(ji,jj) = MAX ( SQRT( pdb_bl(ji,jj) * pdh(ji,jj) / MAX( pdu_bl(ji,jj)**2 + pdv_bl(ji,jj)**2, 1e-8_wp ) ) * & 2386 & ( phbl(ji,jj) / pdh(ji,jj) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 / & 2387 & MAX( zekman(ji,jj), 1.0e-6_wp ) , 5.0_wp ) 2388 IF ( ff_t(ji,jj) >= 0.0_wp ) THEN ! Northern hemisphere 2389 zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 + & 2390 & MAX( -1.0_wp * pdv_ml(ji,jj), 1e-5_wp)**2 ) 2391 ELSE ! Southern hemisphere 2392 zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 + & 2393 & MAX( pdv_ml(ji,jj), 1e-5_wp)**2 ) 2394 END IF 2395 pshear(ji,jj) = za_shr * zekman(ji,jj) * & 2396 & ( MAX( sustar(ji,jj)**2 * pdu_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) + & 2397 & zb_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * & 2398 & pdv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) ) 2399 #ifdef key_osm_debug 2400 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2401 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',pshear(ji,jj) 2402 WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj) 2403 FLUSH(narea+100) 2404 END IF 2405 #endif 2406 ! Stability dependence 2407 pshear(ji,jj) = pshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) ) 2408 #ifdef key_osm_debug 2409 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2410 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',pshear(ji,jj) 2411 FLUSH(narea+100) 2412 END IF 2413 #endif 2414 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2415 ! Test ensures k_ddh=0 is not selected. Change to zri_p<27 when ! 2416 ! full code available ! 2417 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2418 IF ( pshear(ji,jj) > 1e-10 ) THEN 2419 IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN 2420 ! Growing shear layer 2421 k_ddh(ji,jj) = 0 2422 ldshear(ji,jj) = .TRUE. 2423 ELSE 2424 k_ddh(ji,jj) = 1 2425 ! IF ( zri_b <= 1.5 .and. pshear(ji,jj) > 0._wp ) THEN 2426 ! Shear production large enough to determine layer charcteristics, but can't maintain a shear layer 2427 ldshear(ji,jj) = .TRUE. 2428 ! ELSE 2429 END IF 2430 ELSE 2431 k_ddh(ji,jj) = 2 2432 ldshear(ji,jj) = .FALSE. 2433 END IF 2434 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline 2435 ! pshear(ji,jj) = 0.5 * pshear(ji,jj) 2436 ! ldshear(ji,jj) = .FALSE. 2437 ! ENDIF 2438 ELSE ! pdb_bl test, note pshear set to zero 2439 k_ddh(ji,jj) = 2 2440 ldshear(ji,jj) = .FALSE. 2441 ENDIF 2442 ENDIF 2443 END_2D 2444 ! 2445 ! Calculate entrainment buoyancy flux due to surface fluxes. 2446 DO_2D( 0, 0, 0, 0 ) 2447 IF ( ldconv(ji,jj) ) THEN 2448 zwcor = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln 2449 zrf_conv = TANH( ( swstrc(ji,jj) / zwcor )**0.69_wp ) 2450 zrf_shear = TANH( ( sustar(ji,jj) / zwcor )**0.69_wp ) 2451 zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69_wp ) 2452 IF ( nn_osm_SD_reduce > 0 ) THEN 2453 ! Effective Stokes drift already reduced from surface value 2454 zr_stokes = 1.0_wp 2455 ELSE 2456 ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 2457 ! requires further reduction where BL is deep 2458 zr_stokes = 1.0 - EXP( -25.0_wp * dstokes(ji,jj) / hbl(ji,jj) * ( 1.0_wp + 4.0_wp * dstokes(ji,jj) / hbl(ji,jj) ) ) 2459 END IF 2460 pwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) - & 2461 & zalpha_s * zrf_shear * sustar(ji,jj)**3 / phml(ji,jj) + & 2462 & zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 - & 2463 & zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / phml(ji,jj) 2464 #ifdef key_osm_debug 2465 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2466 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',pwb_ent(ji,jj) 2467 FLUSH(narea+100) 2468 END IF 2469 #endif 2470 ENDIF 2471 END_2D 2472 ! 2473 DO_2D( 0, 0, 0, 0 ) 2474 IF ( ldshear(ji,jj) ) THEN 2475 IF ( ldconv(ji,jj) ) THEN 2476 ! Unstable OSBL 2477 zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * pshear(ji,jj) 2478 #ifdef key_osm_debug 2479 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2480 WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr 2481 FLUSH(narea+100) 2482 END IF 2483 #endif 2484 IF ( k_ddh(ji,jj) == 0 ) THEN 2485 ! Developing shear layer, additional shear production possible. 2486 2487 ! pshear_u = MAX( zustar(ji,jj)**2 * MAX( pdu_ml(ji,jj), 0._wp ) / phbl(ji,jj), 0._wp ) 2488 ! pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) 2489 ! pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u ) 2490 2491 ! zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) 2492 ! zwb_shr = MAX( zwb_shr, -0.25 * pshear_u ) 2493 #ifdef key_osm_debug 2494 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2495 WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state k_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, & 2496 & ' zshear=',pshear(ji,jj),' zshear_u=', pshear_u 2497 FLUSH(narea+100) 2498 END IF 2499 #endif 2500 ENDIF 2501 pwb_ent(ji,jj) = pwb_ent(ji,jj) + zwb_shr 2502 ! pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * zwb0(ji,jj) 2503 ELSE ! IF ( ldconv ) THEN - ENDIF 2504 ! Stable OSBL - shear production not coded for first attempt. 2505 ENDIF ! ldconv 2506 END IF ! ldshear 2507 IF ( ldconv(ji,jj) ) THEN 2508 ! Unstable OSBL 2509 pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * 2.0_wp * swbav(ji,jj) 2510 END IF ! ldconv 2511 #ifdef key_osm_debug 2512 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2513 WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',pwb_ent(ji,jj), & 2514 & ' zwb_min=',pwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' swbav= ', swbav(ji,jj) 2515 FLUSH(narea+100) 2516 END IF 2517 #endif 2518 END_2D 2519 ! 2520 IF( ln_timing ) CALL timing_stop('zdf_osm_os') 2521 ! 2522 END SUBROUTINE zdf_osm_osbl_state 2523 2524 SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, kbld, kmld, ldconv, ldshear, ldpyc, ldcoup, k_ddh, pdiffut, pviscos, & 2525 & phbl, phml, pdh, pdhdt, pshear, pb_bl, pu_bl, pv_bl, pb_ml, pu_ml, pv_ml, pwb_ent, & 2526 & pwb_min ) 2527 !!--------------------------------------------------------------------- 2528 !! *** ROUTINE zdf_osm_diffusivity_viscosity *** 2529 !! 2530 !! ** Purpose : Determines the eddy diffusivity and eddy viscosity 2531 !! profiles in the mixed layer and the pycnocline. 2532 !! 2533 !! ** Method : 2534 !! 2535 !!---------------------------------------------------------------------- 2536 INTEGER, INTENT(in ) :: Kbb, Kmm ! Ocean time-level indices 2537 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 2538 INTEGER, DIMENSION(:,:), INTENT(in ) :: kmld ! ML base layer 2539 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 2540 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldshear ! Shear layers 2541 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldpyc ! Pycnocline flags 2542 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldcoup ! Coupling to bottom 2543 INTEGER, DIMENSION(:,:), INTENT(in ) :: k_ddh ! Type of shear layer 2544 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdiffut ! t-diffusivity 2545 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pviscos ! Viscosity 2546 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 2547 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phml ! ML depth 2548 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline depth 2549 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdhdt ! BL depth tendency 2550 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pshear ! Shear production 2551 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the boundary layer 2552 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pu_bl ! Velocity average over the depth of the boundary layer 2553 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pv_bl ! Velocity average over the depth of the boundary layer 2554 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_ml ! Buoyancy average over the depth of the mixed layer 2555 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pu_ml ! Velocity average over the depth of the mixed layer 2556 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pv_ml ! Velocity average over the depth of the mixed layer 2557 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 2558 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_min 2559 ! 2560 ! Local variables 2561 INTEGER :: ji, jj, jk ! Loop indices 2562 ! 2563 ! Scales used to calculate eddy diffusivity and viscosity profiles 2564 REAL(wp), DIMENSION(A2D(0)) :: zdifml_sc, zvisml_sc 2565 REAL(wp), DIMENSION(A2D(0)) :: zdifpyc_n_sc, zdifpyc_s_sc 2566 REAL(wp), DIMENSION(A2D(0)) :: zvispyc_n_sc, zvispyc_s_sc 2567 REAL(wp), DIMENSION(A2D(0)) :: zbeta_d_sc, zbeta_v_sc 2568 REAL(wp), DIMENSION(A2D(0)) :: zb_coup, zc_coup_vis, zc_coup_dif 2569 ! 2570 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b 2571 REAL(wp) :: za_cubic, zb_d_cubic, zc_d_cubic, zd_d_cubic, & ! Coefficients in cubic polynomial specifying diffusivity 2572 & zb_v_cubic, zc_v_cubic, zd_v_cubic ! and viscosity in pycnocline 2573 REAL(wp) :: zznd_ml, zznd_pyc, ztmp 2574 REAL(wp) :: zmsku, zmskv 2575 ! 2576 REAL(wp), PARAMETER :: rn_dif_ml = 0.8_wp, rn_vis_ml = 0.375_wp 2577 REAL(wp), PARAMETER :: rn_dif_pyc = 0.15_wp, rn_vis_pyc = 0.142_wp 2578 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15_wp 2579 ! 2580 IF( ln_timing ) CALL timing_start('zdf_osm_dv') 2581 ! 2582 zb_coup(:,:) = 0.0_wp 2583 ! 2584 DO_2D( 0, 0, 0, 0 ) 2585 IF ( ldconv(ji,jj) ) THEN 2586 ! 2587 zvel_sc_pyc = ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25_wp * pshear(ji,jj) * phbl(ji,jj) )**pthird 2588 zvel_sc_ml = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird 2589 zstab_fac = ( phml(ji,jj) / zvel_sc_ml * & 2590 & ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10(-shol(ji,jj) ) ) )**1.25_wp ) )**2 2591 ! 2592 zdifml_sc(ji,jj) = rn_dif_ml * phml(ji,jj) * zvel_sc_ml 2593 zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj) 2594 #ifdef key_osm_debug 2595 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2596 WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, ldconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) 2597 WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac 2598 FLUSH(narea+100) 2599 END IF 2600 #endif 2601 ! 2602 IF ( ldpyc(ji,jj) ) THEN 2603 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj) 2604 zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 * & 2605 & ( 0.005 * ( pu_ml(ji,jj)-pu_bl(ji,jj) )**2 + 0.0075 * ( pv_ml(ji,jj)-pv_bl(ji,jj) )**2 ) / & 2606 & pdh(ji,jj) 2607 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 2608 #ifdef key_osm_debug 2609 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2610 WRITE(narea+100,'(2(a,g11.3))')' lpyc=ldconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 2611 FLUSH(narea+100) 2612 END IF 2613 #endif 2614 ! 2615 IF ( ldshear(ji,jj) .AND. k_ddh(ji,jj) /= 2 ) THEN 2616 ztmp = rn_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj) 2617 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + ztmp 2618 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + ztmp 2619 ENDIF 2620 #ifdef key_osm_debug 2621 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2622 WRITE(narea+100,'(2(a,g11.3))')' lpyc=ldconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 2623 FLUSH(narea+100) 2624 END IF 2625 #endif 2626 ! 2627 zdifpyc_s_sc(ji,jj) = pwb_ent(ji,jj) + 0.0025_wp * zvel_sc_pyc * ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) * & 2628 & ( pb_ml(ji,jj) - pb_bl(ji,jj) ) 2629 zvispyc_s_sc(ji,jj) = 0.09_wp * ( pwb_min(ji,jj) + 0.0025_wp * zvel_sc_pyc * & 2630 & ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) * & 2631 & ( pb_ml(ji,jj) - pb_bl(ji,jj) ) ) 2632 #ifdef key_osm_debug 2633 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2634 WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 2635 FLUSH(narea+100) 2636 END IF 2637 #endif 2638 zdifpyc_s_sc(ji,jj) = 0.09_wp * zdifpyc_s_sc(ji,jj) * zstab_fac 2639 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 2640 #ifdef key_osm_debug 2641 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2642 WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 2643 FLUSH(narea+100) 2644 END IF 2645 #endif 2646 ! 2647 zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5_wp * zdifpyc_n_sc(ji,jj) ) 2648 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) ) 2649 2650 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) / & 2651 & ( zdifml_sc(ji,jj) + epsln ) )**p2third 2652 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 2653 ELSE 2654 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/03 2655 zdifpyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 2656 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj) ! ag 19/03 2657 zvispyc_s_sc(ji,jj) = 0.0_wp ! ag 19/03 2658 IF(ldcoup(ji,jj) ) THEN ! ag 19/03 2659 ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub| 2660 ! already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90 2661 ! Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub) 2662 ! wet-cell averaging .. 2663 zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 2664 zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 2665 zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) * & 2666 & SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 2667 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) ) 2668 2669 zz_b = -1.0_wp * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 2670 zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / phml(ji,jj) - zb_coup(ji,jj) ) / & 2671 & ( phml(ji,jj) + zz_b ) ! ag 19/03 2672 #ifdef key_osm_debug 2673 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2674 WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj), & 2675 & ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj) 2676 WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv 2677 FLUSH(narea+100) 2678 END IF 2679 #endif 2680 zz_b = -1.0_wp * phml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! ag 19/03 2681 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 ) / & 2682 & zvisml_sc(ji,jj) ! ag 19/03 2683 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) / & 2684 & zdifml_sc(ji,jj) )**p2third 2685 zc_coup_dif(ji,jj) = 0.5_wp * ( -zdifml_sc(ji,jj) / phml(ji,jj) * ( 1.0_wp - zbeta_d_sc(ji,jj) )**1.5_wp + & 2686 & 1.5_wp * ( zdifml_sc(ji,jj) / phml(ji,jj) ) * zbeta_d_sc(ji,jj) * & 2687 & SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b ! ag 19/03 2688 #ifdef key_osm_debug 2689 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2690 WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj) 2691 FLUSH(narea+100) 2692 END IF 2693 #endif 2694 ELSE ! ag 19/03 2695 zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) / & 2696 & ( zdifml_sc(ji,jj) + epsln ) )**p2third ! ag 19/03 2697 zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / & 2698 & ( zvisml_sc(ji,jj) + epsln ) ! ag 19/03 2699 ENDIF ! ag 19/03 2700 ENDIF ! ag 19/03 2701 #ifdef key_osm_debug 2702 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2703 WRITE(narea+100,'(2(a,g11.3))')'ldconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj) 2704 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj) 2705 WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj) 2706 FLUSH(narea+100) 2707 END IF 2708 #endif 2709 ELSE 2710 zdifml_sc(ji,jj) = svstr(ji,jj) * phbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp) 2711 zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 2712 #ifdef key_osm_debug 2713 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2714 WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, ldconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj) 2715 FLUSH(narea+100) 2716 END IF 2717 #endif 2718 END IF 2719 END_2D 2720 ! 2721 DO_2D( 0, 0, 0, 0 ) 2722 IF ( ldconv(ji,jj) ) THEN 2723 DO jk = 2, kmld(ji,jj) ! Mixed layer diffusivity 2724 zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 2725 pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 2726 pviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_v_sc(ji,jj) * zznd_ml ) * & 2727 & ( 1.0_wp - 0.5_wp * zznd_ml**2 ) 2728 END DO 2729 ! 2730 ! Coupling to bottom 2731 ! 2732 IF ( ldcoup(ji,jj) ) THEN ! ag 19/03 2733 DO jk = mbkt(ji,jj), kmld(ji,jj), -1 ! ag 19/03 2734 zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) ! ag 19/03 2735 pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ! ag 19/03 2736 pdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2 ! ag 19/03 2737 END DO ! ag 19/03 2738 ENDIF ! ag 19/03 2739 ! Pycnocline 2740 IF ( ldpyc(ji,jj) ) THEN 2741 ! Diffusivity and viscosity profiles in the pycnocline given by 2742 ! cubic polynomial. Note, if ldpyc TRUE can't be coupled to seabed. 2743 za_cubic = 0.5_wp 2744 zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj) 2745 zd_d_cubic = ( pdh(ji,jj) * zdifml_sc(ji,jj) / phml(ji,jj) * SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) * & 2746 & ( 2.5_wp * zbeta_d_sc(ji,jj) - 1.0_wp ) - 0.85_wp * zdifpyc_s_sc(ji,jj) ) / & 2747 & MAX( zdifpyc_n_sc(ji,jj), 1.0e-8_wp ) 2748 zd_d_cubic = zd_d_cubic - zb_d_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_d_cubic ) 2749 zc_d_cubic = 1.0_wp - za_cubic - zb_d_cubic - zd_d_cubic 2750 zb_v_cubic = -1.75_wp * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 2751 zd_v_cubic = ( 0.5_wp * zvisml_sc(ji,jj) * pdh(ji,jj) / phml(ji,jj) - 0.85_wp * zvispyc_s_sc(ji,jj) ) / & 2752 & MAX( zvispyc_n_sc(ji,jj), 1.0e-8_wp ) 2753 zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic ) 2754 zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic 2755 DO jk = kmld(ji,jj) , kbld(ji,jj) 2756 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp ) 2757 ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ) 2758 ! 2759 pdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * & 2760 & ( za_cubic + zb_d_cubic * zznd_pyc + zc_d_cubic * zznd_pyc**2 + zd_d_cubic * zznd_pyc**3 ) 2761 ! 2762 pdiffut(ji,jj,jk) = pdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ztmp 2763 pviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * & 2764 & ( za_cubic + zb_v_cubic * zznd_pyc + zc_v_cubic * zznd_pyc**2 + zd_v_cubic * zznd_pyc**3 ) 2765 pviscos(ji,jj,jk) = pviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ztmp 2766 END DO 2767 ! IF ( pdhdt(ji,jj) > 0._wp ) THEN 2768 ! zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 2769 ! zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 ) 2770 ! ELSE 2771 ! zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 2772 ! zviscos(ji,jj,ibld(ji,jj)) = 0._wp 2773 ! ENDIF 2774 ENDIF 2775 ELSE 2776 ! Stable conditions 2777 DO jk = 2, kbld(ji,jj) 2778 zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj) 2779 pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp 2780 pviscos(ji,jj,jk) = 0.375_wp * zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml ) * ( 1.0_wp - zznd_ml**2 ) 2781 END DO 2782 ! 2783 IF ( pdhdt(ji,jj) > 0.0_wp ) THEN 2784 pdiffut(ji,jj,kbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, kbld(ji,jj), Kmm) 2785 pviscos(ji,jj,kbld(ji,jj)) = pdiffut(ji,jj,kbld(ji,jj)) 2786 ENDIF 2787 ENDIF ! End if ( ldconv ) 2788 ! 2789 END_2D 2790 IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) ) ! BBL-coupling velocity scale 2791 IF( ln_timing ) CALL timing_stop('zdf_osm_dv') 2792 ! 2793 END SUBROUTINE zdf_osm_diffusivity_viscosity 2754 2794 2755 2795 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, pshear, &
Note: See TracChangeset
for help on using the changeset viewer.
