- Timestamp:
- 2021-05-04T18:05:47+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
r14779 r14783 40 40 !!---------------------------------------------------------------------- 41 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_external_gradients : calculate gradients below the OSBL 48 !! zdf_osm_calculate_dhdt : calculate rate of change of hbl 49 !! zdf_osm_timestep_hbl : hbl timestep 50 !! zdf_osm_diffusivity_viscosity : compute eddy diffusivity and viscosity profiles 51 !! zdf_osm_fgr_terms : compute flux-gradient relationship terms 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_external_gradients : calculate gradients below the OSBL 48 !! zdf_osm_calculate_dhdt : calculate rate of change of hbl 49 !! zdf_osm_timestep_hbl : hbl timestep 50 !! zdf_osm_pycnocline_thickness : calculate thickness of pycnocline 51 !! zdf_osm_pycnocline_buoyancy_profiles : calculate pycnocline buoyancy profiles 52 !! zdf_osm_diffusivity_viscosity : compute eddy diffusivity and viscosity profiles 53 !! zdf_osm_fgr_terms : compute flux-gradient relationship terms 52 54 !! zdf_osm_init : initialization, namelist read, and parameters control 53 55 !! osm_rst : read (or initialize) and write osmosis restart fields … … 874 876 & jp_ext, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 875 877 876 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 878 CALL zdf_osm_pycnocline_thickness( Kmm, ibld, imld, lshear, lconv, j_ddh, zdh, zhml, zdhdt, zdb_bl, & 879 & zhbl, zwb_ent, zdbdz_bl_ext, zwb_fk_b ) 877 880 878 881 ! Reset l_pyc before calculating terms in the flux-gradient relationship … … 945 948 946 949 jp_ext(:,:) = 1 ! ag 19/03 947 CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_pyc ) 950 CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, ibld, jp_ext, lconv, lpyc, zdbdz_pyc, zalpha_pyc, zdh, & 951 & zdb_bl, zhbl, zdbdz_bl_ext, zhml, zdb_ml, zdhdt ) 948 952 #ifdef key_osm_debug 949 953 IF(narea==nn_narea_db) THEN … … 1336 1340 END SUBROUTINE zdf_osm_osbl_state_fk 1337 1341 1338 SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha )1339 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: pdbdz ! Gradients in the pycnocline1340 REAL(wp), DIMENSION(:,:), INTENT( inout ) :: palpha1341 INTEGER :: jk, jj, ji1342 REAL(wp) :: zbgrad1343 REAL(wp) :: zgamma_b_nd, znd1344 REAL(wp) :: zzeta_m1345 REAL(wp), PARAMETER :: ppgamma_b = 2.25_wp1346 !1347 IF( ln_timing ) CALL timing_start('zdf_osm_pscp')1348 !1349 DO_2D( 0, 0, 0, 0 )1350 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN1351 IF ( lconv(ji,jj) ) THEN ! convective conditions1352 IF ( lpyc(ji,jj) ) THEN1353 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )1354 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) * &1355 & zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / &1356 & ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) )1357 palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp )1358 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln )1359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1360 ! Commented lines in this section are not needed in new code, once tested !1361 ! can be removed !1362 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1363 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)1364 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)1365 zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)1366 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)1367 DO jk = 2, ibld(ji,jj)1368 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp1369 IF ( znd <= zzeta_m ) THEN1370 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &1371 ! & EXP( -6.0 * ( znd -zzeta_m )**2 )1372 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &1373 ! & EXP( -6.0 * ( znd -zzeta_m )**2 )1374 pdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + palpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &1375 & EXP( -6.0_wp * ( znd -zzeta_m )**2 )1376 ELSE1377 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )1378 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )1379 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * ppgamma_b * ( znd - zzeta_m )**2 )1380 ENDIF1381 END DO1382 #ifdef key_osm_debug1383 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1384 WRITE(narea+100,'(a,/,3(a,g11.3),/,2(a,g11.3),/)')'end of zdf_osm_pycnocline_buoyancy_profiles:lconv=lpyc=T',&1385 & 'zzeta_m=', zzeta_m, ' zalpha=', palpha(ji,jj), ' ztmp=', ztmp,&1386 & ' zbgrad=', zbgrad, ' zgamma_b_nd=', zgamma_b_nd1387 FLUSH(narea+100)1388 END IF1389 #endif1390 ENDIF ! If no pycnocline pycnocline gradients set to zero1391 ELSE ! Stable conditions1392 ! If pycnocline profile only defined when depth steady of increasing.1393 IF ( zdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady.1394 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN1395 IF ( shol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline1396 ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln )1397 zbgrad = zdb_bl(ji,jj) * ztmp1398 DO jk = 2, ibld(ji,jj)1399 znd = gdepw(ji,jj,jk,Kmm) * ztmp1400 pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )1401 END DO1402 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.1403 ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln )1404 zbgrad = zdb_bl(ji,jj) * ztmp1405 DO jk = 2, ibld(ji,jj)1406 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp1407 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )1408 END DO1409 ENDIF ! IF (shol >=0.5)1410 #ifdef key_osm_debug1411 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1412 WRITE(narea+100,'(1(a,g11.3))')'end of zdf_osm_pycnocline_buoyancy_profiles:lconv=F zbgrad=', zbgrad1413 ! WRITE(narea+100,'(1(a,g11.3))')'end of zdf_osm_pycnocline_scalar_profiles:lconv=F ztgrad=',&1414 ! & ztgrad, ' zsgrad=', zsgrad, ' zbgrad=', zbgrad1415 FLUSH(narea+100)1416 END IF1417 #endif1418 ENDIF ! IF (zdb_bl> 0.)1419 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero1420 ENDIF ! IF (lconv)1421 ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) )1422 END_2D1423 !1424 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles1425 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) )1426 END IF1427 !1428 IF( ln_timing ) CALL timing_stop('zdf_osm_pscp')1429 !1430 END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles1431 1432 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )1433 !!---------------------------------------------------------------------1434 !! *** ROUTINE zdf_osm_pycnocline_thickness ***1435 !!1436 !! ** Purpose : Calculates thickness of the pycnocline1437 !!1438 !! ** Method : The thickness is calculated from a prognostic equation1439 !! that relaxes the pycnocine thickness to a diagnostic1440 !! value. The time change is calculated assuming the1441 !! thickness relaxes exponentially. This is done to deal1442 !! with large timesteps.1443 !!1444 !!----------------------------------------------------------------------1445 1446 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness.1447 !1448 INTEGER :: jj, ji1449 INTEGER :: inhml1450 REAL(wp) :: zari, ztau, zdh_ref, zddhdt, zvel_max1451 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth1452 1453 IF( ln_timing ) CALL timing_start('zdf_osm_pt')1454 DO_2D( 0, 0, 0, 0 )1455 1456 IF ( lshear(ji,jj) ) THEN1457 IF ( lconv(ji,jj) ) THEN1458 IF ( zdb_bl(ji,jj) > 1e-15_wp ) THEN1459 IF ( j_ddh(ji,jj) == 0 ) THEN1460 zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj)1461 ! ddhdt for pycnocline determined in osm_calculate_dhdt1462 zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) )1463 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt1464 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer1465 dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_Dt, 0.625_wp * hbl(ji,jj) )1466 ELSE1467 ! Need to recalculate because hbl has been updated.1468 IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5 ) THEN1469 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1470 ELSE1471 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1472 ENDIF1473 ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt )1474 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) )1475 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)1476 ENDIF1477 ELSE1478 ztau = MAX( MAX( hbl(ji,jj) / ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt )1479 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.2_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) )1480 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj)1481 END IF1482 ELSE ! lconv1483 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL1484 1485 ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln)1486 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here1487 ! boundary layer deepening1488 IF ( zdb_bl(ji,jj) > 0._wp ) THEN1489 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.1490 zari = MIN( 4.5 * ( svstr(ji,jj)**2 ) &1491 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )1492 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)1493 ELSE1494 zdh_ref = 0.2 * hbl(ji,jj)1495 ENDIF1496 ELSE ! IF(dhdt < 0)1497 zdh_ref = 0.2 * hbl(ji,jj)1498 ENDIF ! IF (dhdt >= 0)1499 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )1500 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse1501 ENDIF1502 1503 ELSE ! lshear1504 ! for lshear = .FALSE. calculate ddhdt here1505 1506 IF ( lconv(ji,jj) ) THEN1507 1508 IF( ln_osm_mle ) THEN1509 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN1510 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F1511 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN1512 IF ( ( swstrc(ji,jj) / MAX(svstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability1513 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1514 ELSE ! unstable1515 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1516 ENDIF1517 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird)1518 zdh_ref = zari * hbl(ji,jj)1519 ELSE1520 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird)1521 zdh_ref = 0.2 * hbl(ji,jj)1522 ENDIF1523 ELSE1524 ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird)1525 zdh_ref = 0.2 * hbl(ji,jj)1526 ENDIF1527 ELSE ! ln_osm_mle1528 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN1529 IF ( ( swstrc(ji,jj) / MAX(svstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability1530 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1531 ELSE ! unstable1532 zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )1533 ENDIF1534 ztau = hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird)1535 zdh_ref = zari * hbl(ji,jj)1536 ELSE1537 ztau = hbl(ji,jj) / MAX(epsln, (svstr(ji,jj)**3 + 0.5 *swstrc(ji,jj)**3)**pthird)1538 zdh_ref = 0.2 * hbl(ji,jj)1539 ENDIF1540 1541 END IF ! ln_osm_mle1542 1543 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )1544 ! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref1545 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref1546 ! Alan: this hml is never defined or used1547 ELSE ! IF (lconv)1548 ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln)1549 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here1550 ! boundary layer deepening1551 IF ( zdb_bl(ji,jj) > 0._wp ) THEN1552 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.1553 zari = MIN( 4.5 * ( svstr(ji,jj)**2 ) &1554 & / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )1555 zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)1556 ELSE1557 zdh_ref = 0.2 * hbl(ji,jj)1558 ENDIF1559 ELSE ! IF(dhdt < 0)1560 zdh_ref = 0.2 * hbl(ji,jj)1561 ENDIF ! IF (dhdt >= 0)1562 dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )1563 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse1564 ENDIF ! IF (lconv)1565 ENDIF ! lshear1566 1567 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)1568 inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj) - 1,Kmm), 1.e-3) ) , 1 )1569 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)1570 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)1571 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)1572 #ifdef key_osm_debug1573 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN1574 WRITE(narea+100,'(4(a,g11.3),2(a,i7),/,5(a,g11.3),/)') 'end of zdf_osm_pycnocline_thickness:hml=',hml(ji,jj), &1575 & ' zhml=',zhml(ji,jj),' zdh=', zdh(ji,jj), ' dh=', dh(ji,jj), ' imld=', imld(ji,jj), ' inhml=', inhml, &1576 & 'zvel_max=', zvel_max, ' ztau=', ztau,' zdh_ref=', zdh_ref, ' zar=', zari, ' zddhdt=', zddhdt1577 FLUSH(narea+100)1578 END IF1579 #endif1580 END_2D1581 IF( ln_timing ) CALL timing_stop('zdf_osm_pt')1582 1583 END SUBROUTINE zdf_osm_pycnocline_thickness1584 1585 1586 1342 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 1587 1343 !!---------------------------------------------------------------------- … … 1969 1725 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdv_ml ! Velocity diff. (v) between mixed-layer average and basal value 1970 1726 ! 1971 ! Local Variables1727 ! Local variables 1972 1728 INTEGER :: jj, ji ! Loop indices 1973 1729 ! … … 2225 1981 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_min 2226 1982 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value 2227 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdbdz_bl_ext ! external buoyancy gradients1983 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 2228 1984 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_ml ! Buoyancy diff. between mixed-layer average and basal value 2229 1985 REAL(wp), DIMENSION(:,:), INTENT( out) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL … … 2240 1996 REAL(wp), PARAMETER :: zdhoh = 0.1_wp 2241 1997 REAL(wp), PARAMETER :: zalpha_b = 0.3_wp 2242 REAL(wp), PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5 ! Also in pycnocline_depth1998 REAL(wp), PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5_wp ! Also in pycnocline_depth 2243 1999 REAL(wp), PARAMETER :: zlarge = -1e10_wp 2244 2000 ! … … 2429 2185 ! 2430 2186 ! Local variables 2431 INTEGER :: jk, jj, ji, jm2432 REAL(wp) :: zhbl_s, zvel_max, zdb2433 REAL(wp) :: zthermal, zbeta2187 INTEGER :: jk, jj, ji, jm 2188 REAL(wp) :: zhbl_s, zvel_max, zdb 2189 REAL(wp) :: zthermal, zbeta 2434 2190 ! 2435 2191 IF( ln_timing ) CALL timing_start('zdf_osm_th') … … 2548 2304 ! 2549 2305 END SUBROUTINE zdf_osm_timestep_hbl 2306 2307 SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, kbld, kmld, ldshear, ldconv, k_ddh, pdh, phml, pdhdt, pdb_bl, phbl, pwb_ent, pdbdz_bl_ext, pwb_fk_b ) 2308 !!--------------------------------------------------------------------- 2309 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 2310 !! 2311 !! ** Purpose : Calculates thickness of the pycnocline 2312 !! 2313 !! ** Method : The thickness is calculated from a prognostic equation 2314 !! that relaxes the pycnocine thickness to a diagnostic 2315 !! value. The time change is calculated assuming the 2316 !! thickness relaxes exponentially. This is done to deal 2317 !! with large timesteps. 2318 !! 2319 !!---------------------------------------------------------------------- 2320 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 2321 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 2322 INTEGER, DIMENSION(:,:), INTENT(inout) :: kmld ! ML base layer 2323 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldshear ! Shear layers 2324 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 2325 INTEGER, DIMENSION(:,:), INTENT(in ) :: k_ddh ! k_ddh=0: active shear layer 2326 ! ! k_ddh=1: shear layer not active 2327 ! ! k_ddh=2: shear production low 2328 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdh ! Pycnocline thickness 2329 REAL(wp), DIMENSION(:,:), INTENT(inout) :: phml ! ML depth 2330 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdhdt ! BL depth tendency 2331 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value 2332 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 2333 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 2334 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 2335 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL 2336 2337 ! 2338 ! Local variables 2339 INTEGER :: jj, ji 2340 INTEGER :: inhml 2341 REAL(wp) :: zari, ztau, zdh_ref, zddhdt, zvel_max 2342 REAL(wp) :: ztmp ! Auxiliary variable 2343 ! 2344 REAL, PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5_wp ! Also in pycnocline_depth 2345 ! 2346 IF( ln_timing ) CALL timing_start('zdf_osm_pt') 2347 ! 2348 DO_2D( 0, 0, 0, 0 ) 2349 ! 2350 IF ( ldshear(ji,jj) ) THEN 2351 ! 2352 IF ( ldconv(ji,jj) ) THEN 2353 ! 2354 IF ( pdb_bl(ji,jj) > 1e-15_wp ) THEN 2355 IF ( k_ddh(ji,jj) == 0 ) THEN 2356 zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 2357 ! ddhdt for pycnocline determined in osm_calculate_dhdt 2358 zddhdt = -1.0_wp * a_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) / & 2359 & ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15 ) ) 2360 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt 2361 ! Maximum limit for how thick the shear layer can grow relative to the thickness of the boundary layer 2362 dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_Dt, 0.625_wp * hbl(ji,jj) ) 2363 ELSE ! Need to recalculate because hbl has been updated 2364 IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN 2365 ztmp = svstr(ji,jj) 2366 ELSE 2367 ztmp = swstrc(ji,jj) 2368 END IF 2369 zari = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) + & 2370 & pdb_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2, & 2371 & 1e-12_wp ) ) ), 0.2_wp ) 2372 ztau = MAX( pdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ), & 2373 & 2.0_wp * rn_Dt ) 2374 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + & 2375 & zari * phbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2376 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * phbl(ji,jj) 2377 END IF 2378 ELSE 2379 ztau = MAX( MAX( hbl(ji,jj) / ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt ) 2380 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + & 2381 & 0.2_wp * phbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2382 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj) 2383 END IF 2384 ! 2385 ELSE ! ldconv 2386 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 2387 ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln) 2388 IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN ! Probably shouldn't include wm here 2389 ! Boundary layer deepening 2390 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 2391 ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions 2392 zari = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( pdb_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp, 0.2_wp ) 2393 zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj) 2394 ELSE 2395 zdh_ref = 0.2_wp * hbl(ji,jj) 2396 ENDIF 2397 ELSE ! IF(dhdt < 0) 2398 zdh_ref = 0.2_wp * hbl(ji,jj) 2399 ENDIF ! IF (dhdt >= 0) 2400 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + zdh_ref * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2401 IF ( pdhdt(ji,jj) < 0.0_wp .AND. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! Can be a problem with dh>hbl for 2402 ! ! rapid collapse 2403 ENDIF 2404 ! 2405 ELSE ! ldshear = .FALSE., calculate ddhdt here 2406 ! 2407 IF ( ldconv(ji,jj) ) THEN 2408 ! 2409 IF( ln_osm_mle ) THEN 2410 IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN ! OSBL is deepening. Note wb_fk_b is zero if 2411 ! ! ln_osm_mle=F 2412 IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 2413 IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln) )**3 <= 0.5_wp ) THEN ! Near neutral stability 2414 ztmp = svstr(ji,jj) 2415 ELSE ! Unstable 2416 ztmp = swstrc(ji,jj) 2417 END IF 2418 zari = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) * & 2419 & ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp) + & 2420 & pdb_bl(ji,jj)**2 / MAX(4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 2421 ELSE 2422 zari = 0.2_wp 2423 END IF 2424 ELSE 2425 zari = 0.2_wp 2426 END IF 2427 ztau = 0.2_wp * hbl(ji,jj) / MAX( epsln, ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird ) 2428 zdh_ref = zari * hbl(ji,jj) 2429 ELSE ! ln_osm_mle 2430 IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 2431 IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln ) )**3 <= 0.5_wp ) THEN ! Near neutral stability 2432 ztmp = svstr(ji,jj) 2433 ELSE ! Unstable 2434 ztmp = swstrc(ji,jj) 2435 END IF 2436 zari = MIN( 1.5_wp * pdb_bl(ji,jj) / ( phbl(ji,jj) * & 2437 & ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) + & 2438 & pdb_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp ) 2439 ELSE 2440 zari = 0.2_wp 2441 END IF 2442 ztau = hbl(ji,jj) / MAX( epsln, ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird ) 2443 zdh_ref = zari * hbl(ji,jj) 2444 END IF ! ln_osm_mle 2445 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + zdh_ref * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2446 ! IF ( pdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2447 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref 2448 ! Alan: this hml is never defined or used 2449 ELSE ! IF (ldconv) 2450 ! 2451 ztau = hbl(ji,jj) / MAX( svstr(ji,jj), epsln ) 2452 IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN ! Probably shouldn't include wm here 2453 ! Boundary layer deepening 2454 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 2455 ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions. 2456 zari = MIN( 4.5_wp * ( svstr(ji,jj)**2 ) / MAX( pdb_bl(ji,jj) * phbl(ji,jj), epsln ) + 0.01_wp , 0.2_wp ) 2457 zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj) 2458 ELSE 2459 zdh_ref = 0.2_wp * hbl(ji,jj) 2460 END IF 2461 ELSE ! IF(dhdt < 0) 2462 zdh_ref = 0.2_wp * hbl(ji,jj) 2463 END IF ! IF (dhdt >= 0) 2464 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + zdh_ref * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 2465 IF ( pdhdt(ji,jj) < 0.0_wp .AND. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! Can be a problem with dh>hbl for 2466 ! ! rapid collapse 2467 END IF ! IF (ldconv) 2468 ! 2469 END IF ! ldshear 2470 ! 2471 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 2472 inhml = MAX( INT( dh(ji,jj) / MAX( e3t(ji,jj,kbld(ji,jj)-1,Kmm), 1e-3_wp ) ), 1 ) 2473 kmld(ji,jj) = MAX( kbld(ji,jj) - inhml, 3 ) 2474 phml(ji,jj) = gdepw(ji,jj,kmld(ji,jj),Kmm) 2475 pdh(ji,jj) = phbl(ji,jj) - phml(ji,jj) 2476 #ifdef key_osm_debug 2477 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2478 WRITE(narea+100,'(4(a,g11.3),2(a,i7),/,5(a,g11.3),/)') 'end of zdf_osm_pycnocline_thickness:hml=',hml(ji,jj), & 2479 & ' zhml=',phml(ji,jj),' zdh=', pdh(ji,jj), ' dh=', dh(ji,jj), ' imld=', kmld(ji,jj), ' inhml=', inhml, & 2480 & 'zvel_max=', zvel_max, ' ztau=', ztau,' zdh_ref=', zdh_ref, ' zar=', zari, ' zddhdt=', zddhdt 2481 FLUSH(narea+100) 2482 END IF 2483 #endif 2484 ! 2485 END_2D 2486 ! 2487 IF( ln_timing ) CALL timing_stop('zdf_osm_pt') 2488 ! 2489 END SUBROUTINE zdf_osm_pycnocline_thickness 2490 2491 SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kbld, kp_ext, ldconv, ldpyc, pdbdz, palpha, pdh, pdb_bl, phbl, pdbdz_bl_ext, phml, pdb_ml, pdhdt ) 2492 !!--------------------------------------------------------------------- 2493 !! *** ROUTINE zdf_osm_pycnocline_buoyancy_profiles *** 2494 !! 2495 !! ** Purpose : calculate pycnocline buoyancy profiles 2496 !! 2497 !! ** Method : 2498 !! 2499 !!---------------------------------------------------------------------- 2500 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 2501 INTEGER, DIMENSION(:,:), INTENT(in ) :: kbld ! BL base layer 2502 INTEGER, DIMENSION(:,:), INTENT(in ) :: kp_ext ! External-level offsets 2503 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 2504 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldpyc ! Pycnocline flags 2505 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdbdz ! Gradients in the pycnocline 2506 REAL(wp), DIMENSION(:,:), INTENT(inout) :: palpha 2507 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdh ! Pycnocline thickness 2508 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_bl ! Buoyancy diff. between BL average and basal value 2509 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phbl ! BL depth 2510 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 2511 REAL(wp), DIMENSION(:,:), INTENT(in ) :: phml ! ML depth 2512 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pdb_ml ! Buoyancy diff. between mixed-layer average and basal value 2513 REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pdhdt ! Rates of change of hbl 2514 ! 2515 ! Local variables 2516 INTEGER :: jk, jj, ji 2517 REAL(wp) :: zbgrad 2518 REAL(wp) :: zgamma_b_nd, znd 2519 REAL(wp) :: zzeta_m 2520 REAL(wp) :: ztmp ! Auxiliary variable 2521 ! 2522 REAL(wp), PARAMETER :: ppgamma_b = 2.25_wp 2523 ! 2524 IF( ln_timing ) CALL timing_start('zdf_osm_pscp') 2525 ! 2526 DO_2D( 0, 0, 0, 0 ) 2527 ! 2528 IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 2529 ! 2530 IF ( ldconv(ji,jj) ) THEN ! Convective conditions 2531 ! 2532 IF ( ldpyc(ji,jj) ) THEN 2533 ! 2534 zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) ) 2535 palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) * & 2536 & pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / pdb_ml(ji,jj) ) / & 2537 & ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) ) 2538 palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp ) 2539 ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 2540 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2541 ! Commented lines in this section are not needed in new code, once tested ! 2542 ! can be removed ! 2543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2544 ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj) 2545 ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj) 2546 zbgrad = palpha(ji,jj) * pdb_ml(ji,jj) * ztmp + pdbdz_bl_ext(ji,jj) 2547 zgamma_b_nd = pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / MAX( pdb_ml(ji,jj), epsln ) 2548 DO jk = 2, kbld(ji,jj) 2549 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) * ztmp 2550 IF ( znd <= zzeta_m ) THEN 2551 ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * & 2552 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 2553 ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * & 2554 ! & EXP( -6.0 * ( znd -zzeta_m )**2 ) 2555 pdbdz(ji,jj,jk) = pdbdz_bl_ext(ji,jj) + palpha(ji,jj) * pdb_ml(ji,jj) * ztmp * & 2556 & EXP( -6.0_wp * ( znd -zzeta_m )**2 ) 2557 ELSE 2558 ! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 2559 ! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 ) 2560 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.0_wp * ppgamma_b * ( znd - zzeta_m )**2 ) 2561 END IF 2562 END DO 2563 #ifdef key_osm_debug 2564 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2565 WRITE(narea+100,'(a,/,3(a,g11.3),/,2(a,g11.3),/)')'end of zdf_osm_pycnocline_buoyancy_profiles:lconv=lpyc=T',& 2566 & 'zzeta_m=', zzeta_m, ' zalpha=', palpha(ji,jj), ' ztmp=', ztmp,& 2567 & ' zbgrad=', zbgrad, ' zgamma_b_nd=', zgamma_b_nd 2568 FLUSH(narea+100) 2569 END IF 2570 #endif 2571 END IF ! If no pycnocline pycnocline gradients set to zero 2572 ! 2573 ELSE ! Stable conditions 2574 ! If pycnocline profile only defined when depth steady of increasing. 2575 IF ( pdhdt(ji,jj) > 0.0_wp ) THEN ! Depth increasing, or steady. 2576 IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN 2577 IF ( shol(ji,jj) >= 0.5_wp ) THEN ! Very stable - 'thick' pycnocline 2578 ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln ) 2579 zbgrad = pdb_bl(ji,jj) * ztmp 2580 DO jk = 2, kbld(ji,jj) 2581 znd = gdepw(ji,jj,jk,Kmm) * ztmp 2582 pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 ) 2583 END DO 2584 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 2585 ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln ) 2586 zbgrad = pdb_bl(ji,jj) * ztmp 2587 DO jk = 2, kbld(ji,jj) 2588 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp 2589 pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 ) 2590 END DO 2591 END IF ! IF (shol >=0.5) 2592 #ifdef key_osm_debug 2593 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 2594 WRITE(narea+100,'(1(a,g11.3))')'end of zdf_osm_pycnocline_buoyancy_profiles:lconv=F zbgrad=', zbgrad 2595 ! WRITE(narea+100,'(1(a,g11.3))')'end of zdf_osm_pycnocline_scalar_profiles:lconv=F ztgrad=',& 2596 ! & ztgrad, ' zsgrad=', zsgrad, ' zbgrad=', zbgrad 2597 FLUSH(narea+100) 2598 END IF 2599 #endif 2600 END IF ! IF (pdb_bl> 0.) 2601 END IF ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 2602 ! 2603 END IF ! IF (ldconv) 2604 ! 2605 END IF ! IF ( kbld(ji,jj) < mbkt(ji,jj) ) 2606 ! 2607 END_2D 2608 ! 2609 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 2610 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) ) 2611 END IF 2612 ! 2613 IF( ln_timing ) CALL timing_stop('zdf_osm_pscp') 2614 ! 2615 END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 2550 2616 2551 2617 SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, kbld, kmld, ldconv, ldshear, ldpyc, ldcoup, k_ddh, pdiffut, pviscos, &
Note: See TracChangeset
for help on using the changeset viewer.
