Changeset 14785
- Timestamp:
- 2021-05-04T21:46:50+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
r14783 r14785 52 52 !! zdf_osm_diffusivity_viscosity : compute eddy diffusivity and viscosity profiles 53 53 !! zdf_osm_fgr_terms : compute flux-gradient relationship terms 54 !! zdf_osm_zmld_horizontal_gradients : calculate horizontal buoyancy gradients for use with Fox-Kemper parametrization 55 !! zdf_osm_osbl_state_fk : determine state of OSBL and MLE layers 56 !! zdf_osm_mle_parameters : timestep MLE depth and calculate MLE fluxes 54 57 !! zdf_osm_init : initialization, namelist read, and parameters control 55 58 !! osm_rst : read (or initialize) and write osmosis restart fields … … 746 749 747 750 !! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients 748 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )751 CALL zdf_osm_zmld_horizontal_gradients( Kmm, ibld, zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 749 752 !! Calculate max vertical FK flux zwb_fk & set logical descriptors 750 CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 753 CALL zdf_osm_osbl_state_fk( Kmm, ibld, lconv, lpyc, lflux, lmle, zwb_fk, zt_mle, zs_mle, zb_mle, zhbl, & 754 & zhmle, zwb_ent, zt_bl, zs_bl, zb_bl, zdb_bl, zdbds_mle ) 751 755 !! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 752 CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 756 CALL zdf_osm_mle_parameters( Kmm, lconv, lmle, mld_prof, zmld, zhmle, & 757 & zvel_mle, zdiff_mle, zdbds_mle, zhbl, zb_bl, zwb0tot ) 753 758 #ifdef key_osm_debug 754 759 IF(narea==nn_narea_db) THEN … … 1217 1222 END IF 1218 1223 IF( ln_timing ) CALL timing_stop('zdf_osm') 1219 1220 CONTAINS1221 1222 SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )1223 !!---------------------------------------------------------------------1224 !! *** ROUTINE zdf_osm_osbl_state_fk ***1225 !!1226 !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme.1227 !! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined1228 !! lflux :: determines whether effects of surface flux extend below the base of the OSBL1229 !! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.1230 !!1231 !! ** Method :1232 !!1233 !!1234 !!----------------------------------------------------------------------1235 1236 ! Outputs1237 LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle1238 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk1239 !1240 REAL(wp), DIMENSION(jpi,jpj) :: znd_param1241 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer1242 REAL(wp) :: zpe_mle_ref, zdbdz_mle_int1243 1244 IF( ln_timing ) CALL timing_start('zdf_osm_osf')1245 znd_param(A2D(0)) = 0.0_wp1246 1247 DO_2D( 0, 0, 0, 0 )1248 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf1249 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)1250 END_2D1251 DO_2D( 0, 0, 0, 0 )1252 !1253 IF ( lconv(ji,jj) ) THEN1254 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN1255 zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )1256 zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )1257 zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )1258 zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )1259 ! Calculate potential energies of actual profile and reference profile.1260 zpe_mle_layer = 0._wp1261 zpe_mle_ref = 0._wp1262 zthermal = rab_n(ji,jj,1,jp_tem)1263 zbeta = rab_n(ji,jj,1,jp_sal)1264 DO jk = ibld(ji,jj), mld_prof(ji,jj)1265 zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) )1266 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)1267 zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)1268 END DO1269 ! Non-dimensional parameter to diagnose the presence of thermocline1270 1271 znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) )1272 ENDIF1273 ENDIF1274 #ifdef key_osm_debug1275 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1276 WRITE(narea+100,'(4(a,g11.3))')'start of zdf_osm_osbl_state_fk: zwb_fk=',zwb_fk(ji,jj), &1277 & ' znd_param=',znd_param(ji,jj), ' zpe_mle_ref=', zpe_mle_ref, ' zpe_mle_layer=', zpe_mle_layer1278 FLUSH(narea+100)1279 END IF1280 #endif1281 END_2D1282 1283 ! Diagnosis1284 DO_2D( 0, 0, 0, 0 )1285 IF ( lconv(ji,jj) ) THEN1286 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN1287 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN1288 ! MLE layer growing1289 IF ( znd_param (ji,jj) > 100. ) THEN1290 ! Thermocline present1291 lflux(ji,jj) = .FALSE.1292 lmle(ji,jj) =.FALSE.1293 ELSE1294 ! Thermocline not present1295 lflux(ji,jj) = .TRUE.1296 lmle(ji,jj) = .TRUE.1297 ENDIF ! znd_param > 1001298 !1299 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN1300 lpyc(ji,jj) = .FALSE.1301 ELSE1302 lpyc(ji,jj) = .TRUE.1303 ENDIF1304 ELSE1305 ! MLE layer restricted to OSBL or just below.1306 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN1307 ! Weak stratification MLE layer can grow.1308 lpyc(ji,jj) = .FALSE.1309 lflux(ji,jj) = .TRUE.1310 lmle(ji,jj) = .TRUE.1311 ELSE1312 ! Strong stratification1313 lpyc(ji,jj) = .TRUE.1314 lflux(ji,jj) = .FALSE.1315 lmle(ji,jj) = .FALSE.1316 ENDIF ! zdb_bl < rn_mle_thresh_bl and1317 ENDIF ! zhmle > 1.2 zhbl1318 ELSE1319 lpyc(ji,jj) = .TRUE.1320 lflux(ji,jj) = .FALSE.1321 lmle(ji,jj) = .FALSE.1322 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.1323 ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.51324 ELSE1325 ! Stable Boundary Layer1326 lpyc(ji,jj) = .FALSE.1327 lflux(ji,jj) = .FALSE.1328 lmle(ji,jj) = .FALSE.1329 ENDIF ! lconv1330 #ifdef key_osm_debug1331 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1332 WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',zwb_ent(ji,jj), &1333 & ' zhmle=',zhmle(ji,jj), ' zhbl=', zhbl(ji,jj), &1334 & ' lpyc= ', lpyc(ji,jj), ' lflux= ', lflux(ji,jj), ' lmle= ', lmle(ji,jj), ' lconv= ', lconv(ji,jj)1335 FLUSH(narea+100)1336 END IF1337 #endif1338 END_2D1339 IF( ln_timing ) CALL timing_stop('zdf_osm_osf')1340 END SUBROUTINE zdf_osm_osbl_state_fk1341 1342 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )1343 !!----------------------------------------------------------------------1344 !! *** ROUTINE zdf_osm_horizontal_gradients ***1345 !!1346 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.1347 !!1348 !! ** Method :1349 !!1350 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 20081351 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 20081352 1353 1354 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points1355 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.1356 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !1357 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy1358 1359 INTEGER :: ji, jj, jk ! dummy loop indices1360 INTEGER :: ii, ij, ik, ikmax ! local integers1361 REAL(wp) :: zc1362 REAL(wp) :: zN2_c ! local buoyancy difference from 10m value1363 REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH1364 REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv1365 REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv1366 !!----------------------------------------------------------------------1367 !1368 IF( ln_timing ) CALL timing_start('zdf_osm_zhg')1369 ! !== MLD used for MLE ==!1370 1371 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point1372 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^21373 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria1374 DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )1375 ikt = mbkt(ji,jj)1376 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)1377 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level1378 END_3D1379 DO_2D( 1, 1, 1, 1 )1380 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))1381 zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)1382 END_2D1383 ! ensure mld_prof .ge. ibld1384 !1385 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation1386 !1387 ztm(:,:) = 0._wp1388 zsm(:,:) = 0._wp1389 DO_3D( 1, 1, 1, 1, 1, ikmax )1390 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points1391 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm)1392 zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm)1393 END_3D1394 ! average temperature and salinity.1395 ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )1396 zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )1397 ! calculate horizontal gradients at u & v points1398 1399 zmld_midu(:,:) = 0.0_wp1400 ztsm_midu(:,:,:) = 10.0_wp1401 DO_2D( 0, 0, 1, 0 )1402 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)1403 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)1404 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))1405 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )1406 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )1407 END_2D1408 1409 zmld_midv(:,:) = 0.0_wp1410 ztsm_midv(:,:,:) = 10.0_wp1411 DO_2D( 1, 0, 0, 0 )1412 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)1413 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)1414 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))1415 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )1416 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )1417 END_2D1418 1419 CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm)1420 CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm)1421 1422 DO_2D( 0, 0, 1, 0 )1423 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))1424 END_2D1425 DO_2D( 1, 0, 0, 0 )1426 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))1427 END_2D1428 1429 DO_2D( 0, 0, 0, 0 )1430 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf1431 zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &1432 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )1433 END_2D1434 IF( ln_timing ) CALL timing_stop('zdf_osm_zhg')1435 1436 END SUBROUTINE zdf_osm_zmld_horizontal_gradients1437 SUBROUTINE zdf_osm_mle_parameters( pmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )1438 !!----------------------------------------------------------------------1439 !! *** ROUTINE zdf_osm_mle_parameters ***1440 !!1441 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.1442 !!1443 !! ** Method :1444 !!1445 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 20081446 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 20081447 1448 REAL(wp), DIMENSION(jpi,jpj) :: pmld ! == estimated FK BLD used for MLE horiz gradients == !1449 INTEGER, DIMENSION(jpi,jpj) :: mld_prof1450 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle1451 INTEGER :: ji, jj, jk ! dummy loop indices1452 INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers1453 INTEGER , DIMENSION(jpi,jpj) :: inml_mle1454 REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle1455 1456 IF( ln_timing ) CALL timing_start('zdf_osm_mp')1457 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.1458 1459 DO_2D( 0, 0, 0, 0 )1460 IF ( lconv(ji,jj) ) THEN1461 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf1462 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.1463 zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)1464 zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**21465 ENDIF1466 END_2D1467 ! Timestep mixed layer eddy depth.1468 DO_2D( 0, 0, 0, 0 )1469 IF ( lmle(ji,jj) ) THEN ! MLE layer growing.1470 ! Buoyancy gradient at base of MLE layer.1471 zthermal = rab_n(ji,jj,1,jp_tem)1472 zbeta = rab_n(ji,jj,1,jp_sal)1473 jkb = mld_prof(ji,jj)1474 jkb1 = MIN(jkb + 1, mbkt(ji,jj))1475 !1476 zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) )1477 zdb_mle = zb_bl(ji,jj) - zbuoy1478 ! Timestep hmle.1479 hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_Dt / zdb_mle1480 ELSE1481 IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN1482 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau1483 ELSE1484 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau1485 ENDIF1486 ENDIF1487 hmle(ji,jj) = MAX( MIN( hmle(ji,jj), ht(ji,jj) ), gdepw(ji,jj,4,Kmm) )1488 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN( hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) )1489 ! For now try just set hmle to zmld1490 hmle(ji,jj) = pmld(ji,jj)1491 END_2D1492 1493 mld_prof = 41494 DO_3D( 0, 0, 0, 0, 5, jpkm1 )1495 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)1496 END_3D1497 DO_2D( 0, 0, 0, 0 )1498 zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm)1499 END_2D1500 IF( ln_timing ) CALL timing_stop('zdf_osm_mp')1501 END SUBROUTINE zdf_osm_mle_parameters1502 1224 1503 1225 END SUBROUTINE zdf_osm … … 3514 3236 END SUBROUTINE zdf_osm_fgr_terms 3515 3237 3238 SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, kbld, pmld, pdtdx, pdtdy, pdsdx, pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle ) 3239 !!---------------------------------------------------------------------- 3240 !! *** ROUTINE zdf_osm_zmld_horizontal_gradients *** 3241 !! 3242 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with 3243 !! Fox-Kemper parametrization 3244 !! 3245 !! ** Method : 3246 !! 3247 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 3248 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 3249 !! 3250 !!---------------------------------------------------------------------- 3251 INTEGER, INTENT(in ) :: Kmm ! Time-level index 3252 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 3253 REAL(wp), DIMENSION(:,:), INTENT( out) :: pmld ! == Estimated FK BLD used for MLE horizontal gradients == ! 3254 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdtdx ! Horizontal gradient for Fox-Kemper parametrization 3255 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdtdy ! Horizontal gradient for Fox-Kemper parametrization 3256 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdsdx ! Horizontal gradient for Fox-Kemper parametrization 3257 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdsdy ! Horizontal gradient for Fox-Kemper parametrization 3258 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdx_mle ! MLE horiz gradients at u points 3259 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdy_mle ! MLE horiz gradients at v points 3260 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 3261 ! 3262 ! Local variables 3263 INTEGER :: ji, jj, jk ! Dummy loop indices 3264 INTEGER :: ikt, ikmax ! Local integers 3265 REAL(wp) :: zc 3266 REAL(wp) :: zN2_c ! Local buoyancy difference from 10m value 3267 REAL(wp), DIMENSION(A2D(1)) :: ztm 3268 REAL(wp), DIMENSION(A2D(1)) :: zsm 3269 REAL(wp), DIMENSION(A2D(1),jpts) :: ztsm_midu 3270 REAL(wp), DIMENSION(A2D(1),jpts) :: ztsm_midv 3271 REAL(wp), DIMENSION(A2D(1),jpts) :: zabu 3272 REAL(wp), DIMENSION(A2D(1),jpts) :: zabv 3273 REAL(wp), DIMENSION(A2D(1)) :: zmld_midu 3274 REAL(wp), DIMENSION(A2D(1)) :: zmld_midv 3275 ! 3276 IF( ln_timing ) CALL timing_start('zdf_osm_zhg') 3277 ! 3278 ! == MLD used for MLE ==! 3279 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 3280 pmld(:,:) = 0.0_wp ! Here hmlp used as a dummy variable, integrating vertically N^2 3281 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! Convert density criteria into N^2 criteria 3282 DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 3283 ikt = mbkt(ji,jj) 3284 pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm) 3285 IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 3286 END_3D 3287 DO_2D( 1, 1, 1, 1 ) 3288 mld_prof(ji,jj) = MAX( mld_prof(ji,jj), kbld(ji,jj) ) ! Ensure mld_prof .ge. kbld 3289 pmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 3290 END_2D 3291 ! 3292 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! Max level of the computation 3293 ztm(:,:) = 0.0_wp 3294 zsm(:,:) = 0.0_wp 3295 DO_3D( 1, 1, 1, 1, 1, ikmax ) 3296 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj) - jk ), 1 ), KIND=wp ) ! zc being 0 outside the ML t-points 3297 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) 3298 zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) 3299 END_3D 3300 ! Average temperature and salinity 3301 ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) ) 3302 zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) ) 3303 ! Calculate horizontal gradients at u & v points 3304 zmld_midu(:,:) = 0.0_wp 3305 ztsm_midu(:,:,:) = 10.0_wp 3306 DO_2D( 0, 0, 1, 0 ) 3307 pdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm(ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 3308 pdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm(ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 3309 zmld_midu(ji,jj) = 0.25_wp * ( pmld(ji+1,jj) + pmld(ji,jj)) 3310 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm( ji+1,jj) + ztm( ji,jj) ) 3311 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm( ji+1,jj) + zsm( ji,jj) ) 3312 END_2D 3313 zmld_midv(:,:) = 0.0_wp 3314 ztsm_midv(:,:,:) = 10.0_wp 3315 DO_2D( 1, 0, 0, 0 ) 3316 pdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 3317 pdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 3318 zmld_midv(ji,jj) = 0.25_wp * ( pmld(ji,jj+1) + pmld( ji,jj) ) 3319 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm( ji,jj+1) + ztm( ji,jj) ) 3320 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm( ji,jj+1) + zsm( ji,jj) ) 3321 END_2D 3322 CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm ) 3323 CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm ) 3324 DO_2D( 0, 0, 1, 0 ) 3325 pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) ) 3326 END_2D 3327 DO_2D( 1, 0, 0, 0 ) 3328 pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) ) 3329 END_2D 3330 DO_2D( 0, 0, 0, 0 ) 3331 pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( pdbdx_mle(ji, jj) * pdbdx_mle(ji, jj) + pdbdy_mle(ji,jj ) * pdbdy_mle(ji,jj ) + & 3332 & pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) ) 3333 END_2D 3334 ! 3335 IF( ln_timing ) CALL timing_stop('zdf_osm_zhg') 3336 ! 3337 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 3338 3339 SUBROUTINE zdf_osm_osbl_state_fk( Kmm, kbld, ldconv, ldpyc, ldflux, ldmle, pwb_fk, pt_mle, ps_mle, pb_mle, phbl, & 3340 & phmle, pwb_ent, pt_bl, ps_bl, pb_bl, pdb_bl, pdbds_mle ) 3341 !!--------------------------------------------------------------------- 3342 !! *** ROUTINE zdf_osm_osbl_state_fk *** 3343 !! 3344 !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is 3345 !! returned in the logicals ldpyc, ldflux and ldmle. Used 3346 !! with Fox-Kemper scheme. 3347 !! ldpyc :: determines whether pycnocline flux-grad 3348 !! relationship needs to be determined 3349 !! ldflux :: determines whether effects of surface flux 3350 !! extend below the base of the OSBL 3351 !! ldmle :: determines whether the layer with MLE is 3352 !! increasing with time or if base is relaxing 3353 !! towards hbl 3354 !! 3355 !! ** Method : 3356 !! 3357 !!---------------------------------------------------------------------- 3358 ! Outputs 3359 INTEGER, INTENT(in ) :: Kmm ! Time-level index 3360 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 3361 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 3362 LOGICAL, DIMENSION(:,:), INTENT(inout) :: ldpyc 3363 LOGICAL, DIMENSION(:,:), INTENT(inout) :: ldflux ! Surface flux extends below OSBL into MLE layer 3364 LOGICAL, DIMENSION(:,:), INTENT(inout) :: ldmle ! MLE layer increases in thickness 3365 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pwb_fk 3366 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_mle ! Temperature average over the depth of the MLE layer 3367 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ps_mle ! Salinity average over the depth of the MLE layer 3368 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pb_mle ! Buoyancy average over the depth of the MLE layer 3369 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 3370 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phmle ! MLE depth 3371 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 3372 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pt_bl ! Temperature average over the depth of the blayer 3373 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ps_bl ! Salinity average over the depth of the blayer 3374 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the blayer 3375 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Difference between blayer average and parameter at base of blayer 3376 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 3377 ! 3378 ! Local variables 3379 INTEGER :: ji, jj, jk ! Dummy loop indices 3380 REAL(wp), DIMENSION(A2D(0)) :: znd_param 3381 REAL(wp) :: zthermal, zbeta 3382 REAL(wp) :: zbuoy 3383 REAL(wp) :: ztmp 3384 REAL(wp) :: zpe_mle_layer 3385 REAL(wp) :: zpe_mle_ref 3386 REAL(wp) :: zdbdz_mle_int 3387 ! 3388 IF( ln_timing ) CALL timing_start('zdf_osm_osf') 3389 ! 3390 znd_param(A2D(0)) = 0.0_wp 3391 ! 3392 DO_2D( 0, 0, 0, 0 ) 3393 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 3394 pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj) 3395 END_2D 3396 ! 3397 DO_2D( 0, 0, 0, 0 ) 3398 ! 3399 IF ( ldconv(ji,jj) ) THEN 3400 IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN 3401 pt_mle(ji,jj) = ( pt_mle(ji,jj) * phmle(ji,jj) - pt_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3402 ps_mle(ji,jj) = ( ps_mle(ji,jj) * phmle(ji,jj) - ps_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3403 pb_mle(ji,jj) = ( pb_mle(ji,jj) * phmle(ji,jj) - pb_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3404 zdbdz_mle_int = ( pb_bl(ji,jj) - ( 2.0_wp * pb_mle(ji,jj) - pb_bl(ji,jj) ) ) / ( phmle(ji,jj) - phbl(ji,jj) ) 3405 ! Calculate potential energies of actual profile and reference profile 3406 zpe_mle_layer = 0.0_wp 3407 zpe_mle_ref = 0.0_wp 3408 zthermal = rab_n(ji,jj,1,jp_tem) 3409 zbeta = rab_n(ji,jj,1,jp_sal) 3410 DO jk = kbld(ji,jj), mld_prof(ji,jj) 3411 zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) ) 3412 zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 3413 zpe_mle_ref = zpe_mle_ref + ( pb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) ) * & 3414 & gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 3415 END DO 3416 ! Non-dimensional parameter to diagnose the presence of thermocline 3417 znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / & 3418 & ( MAX( pwb_fk(ji,jj), 1e-10 ) * phmle(ji,jj) ) 3419 END IF 3420 END IF 3421 #ifdef key_osm_debug 3422 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 3423 WRITE(narea+100,'(4(a,g11.3))')'start of zdf_osm_osbl_state_fk: zwb_fk=',pwb_fk(ji,jj), & 3424 & ' znd_param=',znd_param(ji,jj), ' zpe_mle_ref=', zpe_mle_ref, ' zpe_mle_layer=', zpe_mle_layer 3425 FLUSH(narea+100) 3426 END IF 3427 #endif 3428 ! 3429 END_2D 3430 ! 3431 ! Diagnosis 3432 DO_2D( 0, 0, 0, 0 ) 3433 ! 3434 IF ( ldconv(ji,jj) ) THEN 3435 IF ( -2.0_wp * pwb_fk(ji,jj) / pwb_ent(ji,jj) > 0.5_wp ) THEN 3436 IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN ! MLE layer growing 3437 IF ( znd_param (ji,jj) > 100.0_wp ) THEN ! Thermocline present 3438 ldflux(ji,jj) = .FALSE. 3439 ldmle(ji,jj) = .FALSE. 3440 ELSE ! Thermocline not present 3441 ldflux(ji,jj) = .TRUE. 3442 ldmle(ji,jj) = .TRUE. 3443 ENDIF ! znd_param > 100 3444 ! 3445 IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN 3446 ldpyc(ji,jj) = .FALSE. 3447 ELSE 3448 ldpyc(ji,jj) = .TRUE. 3449 ENDIF 3450 ELSE ! MLE layer restricted to OSBL or just below 3451 IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN ! Weak stratification MLE layer can grow 3452 ldpyc(ji,jj) = .FALSE. 3453 ldflux(ji,jj) = .TRUE. 3454 ldmle(ji,jj) = .TRUE. 3455 ELSE ! Strong stratification 3456 ldpyc(ji,jj) = .TRUE. 3457 ldflux(ji,jj) = .FALSE. 3458 ldmle(ji,jj) = .FALSE. 3459 END IF ! pdb_bl < rn_mle_thresh_bl and 3460 END IF ! phmle > 1.2 phbl 3461 ELSE 3462 ldpyc(ji,jj) = .TRUE. 3463 ldflux(ji,jj) = .FALSE. 3464 ldmle(ji,jj) = .FALSE. 3465 IF ( pdb_bl(ji,jj) < rn_osm_bl_thresh ) ldpyc(ji,jj) = .FALSE. 3466 END IF ! -2.0 * pwb_fk(ji,jj) / pwb_ent > 0.5 3467 ELSE ! Stable Boundary Layer 3468 ldpyc(ji,jj) = .FALSE. 3469 ldflux(ji,jj) = .FALSE. 3470 ldmle(ji,jj) = .FALSE. 3471 END IF ! ldconv 3472 #ifdef key_osm_debug 3473 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 3474 WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',pwb_ent(ji,jj), & 3475 & ' zhmle=',phmle(ji,jj), ' zhbl=', phbl(ji,jj), & 3476 & ' lpyc= ', ldpyc(ji,jj), ' lflux= ', ldflux(ji,jj), ' lmle= ', ldmle(ji,jj), ' lconv= ', ldconv(ji,jj) 3477 FLUSH(narea+100) 3478 END IF 3479 #endif 3480 ! 3481 END_2D 3482 ! 3483 IF( ln_timing ) CALL timing_stop('zdf_osm_osf') 3484 ! 3485 END SUBROUTINE zdf_osm_osbl_state_fk 3486 3487 SUBROUTINE zdf_osm_mle_parameters( Kmm, ldconv, ldmle, kmld_prof, pmld, phmle, pvel_mle, pdiff_mle, pdbds_mle, phbl, pb_bl, pwb0tot ) 3488 !!---------------------------------------------------------------------- 3489 !! *** ROUTINE zdf_osm_mle_parameters *** 3490 !! 3491 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates 3492 !! the mixed layer eddy fluxes for buoyancy, heat and 3493 !! salinity. 3494 !! 3495 !! ** Method : 3496 !! 3497 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 3498 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 3499 !! 3500 !!---------------------------------------------------------------------- 3501 INTEGER, INTENT(in ) :: Kmm ! Time-level index 3502 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 3503 LOGICAL, DIMENSION(:,:), INTENT(inout) :: ldmle ! MLE layer increases in thickness 3504 INTEGER, DIMENSION(:,:), INTENT(inout) :: kmld_prof 3505 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pmld ! == Estimated FK BLD used for MLE horiz gradients == ! 3506 REAL(wp), DIMENSION(:,:), INTENT(inout) :: phmle ! MLE depth 3507 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pvel_mle ! Velocity scale for dhdt with stable ML and FK 3508 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdiff_mle ! Extra MLE vertical diff 3509 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 3510 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 3511 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pb_bl ! Buoyancy average over the depth of the blayer 3512 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pwb0tot ! Total surface buoyancy flux including insolation 3513 ! 3514 ! Local variables 3515 INTEGER :: ji, jj, jk ! Dummy loop indices 3516 REAL(wp) :: ztmp 3517 REAL(wp) :: zdbdz 3518 REAL(wp) :: zdtdz 3519 REAL(wp) :: zdsdz 3520 REAL(wp) :: zthermal 3521 REAL(wp) :: zbeta 3522 REAL(wp) :: zbuoy 3523 REAL(wp) :: zdb_mle 3524 ! 3525 IF( ln_timing ) CALL timing_start('zdf_osm_mp') 3526 ! 3527 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE 3528 DO_2D( 0, 0, 0, 0 ) 3529 IF ( ldconv(ji,jj) ) THEN 3530 ztmp = r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf 3531 ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt 3532 pvel_mle(ji,jj) = pdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 3533 pdiff_mle(ji,jj) = 5e-4_wp * rn_osm_mle_ce * ztmp * pdbds_mle(ji,jj) * phmle(ji,jj)**2 3534 END IF 3535 END_2D 3536 ! Timestep mixed layer eddy depth 3537 DO_2D( 0, 0, 0, 0 ) 3538 IF ( ldmle(ji,jj) ) THEN ! MLE layer growing 3539 ! Buoyancy gradient at base of MLE layer 3540 zthermal = rab_n(ji,jj,1,jp_tem) 3541 zbeta = rab_n(ji,jj,1,jp_sal) 3542 zbuoy = grav * ( zthermal * ts(ji,jj,kmld_prof(ji,jj)+2,jp_tem,Kmm) - & 3543 & zbeta * ts(ji,jj,kmld_prof(ji,jj)+2,jp_sal,Kmm) ) 3544 zdb_mle = pb_bl(ji,jj) - zbuoy 3545 ! Timestep hmle 3546 hmle(ji,jj) = hmle(ji,jj) + pwb0tot(ji,jj) * rn_Dt / zdb_mle 3547 ELSE 3548 IF ( phmle(ji,jj) > phbl(ji,jj) ) THEN 3549 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau 3550 ELSE 3551 hmle(ji,jj) = hmle(ji,jj) - 10.0_wp * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau 3552 END IF 3553 END IF 3554 hmle(ji,jj) = MAX( MIN( hmle(ji,jj), ht(ji,jj) ), gdepw(ji,jj,4,Kmm) ) 3555 IF ( ln_osm_hmle_limit ) hmle(ji,jj) = MIN( hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) ) 3556 hmle(ji,jj) = pmld(ji,jj) ! For now try just set hmle to pmld 3557 END_2D 3558 ! 3559 kmld_prof(:,:) = 4 3560 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 3561 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) kmld_prof(ji,jj) = MIN( mbkt(ji,jj), jk ) 3562 END_3D 3563 DO_2D( 0, 0, 0, 0 ) 3564 phmle(ji,jj) = gdepw(ji,jj,kmld_prof(ji,jj),Kmm) 3565 END_2D 3566 ! 3567 IF( ln_timing ) CALL timing_stop('zdf_osm_mp') 3568 ! 3569 END SUBROUTINE zdf_osm_mle_parameters 3570 3516 3571 SUBROUTINE zdf_osm_init( Kmm ) 3517 3572 !!---------------------------------------------------------------------- … … 3764 3819 END SUBROUTINE zdf_osm_init 3765 3820 3766 3767 3821 SUBROUTINE osm_rst( kt, Kmm, cdrw ) 3768 3822 !!---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
