Changeset 14889
- Timestamp:
- 2021-05-19T17:20:25+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE
- Files:
-
- 2 edited
-
TRA/tramle.F90 (modified) (1 diff)
-
ZDF/zdfosm.F90 (modified) (108 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/TRA/tramle.F90
r14856 r14889 366 366 r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 367 367 ! 368 ! Specifically, dbdx_mle, dbdy_mle and mld_prof need to be defined for nn_hls = 2369 IF( nn_hls == 2 .AND. ln_osm_mle .AND. ln_zdfosm ) THEN370 CALL ctl_stop('nn_hls = 2 cannot be used with ln_mle = ln_osm_mle = ln_zdfosm = T (zdfosm not updated for nn_hls = 2)')371 ENDIF372 368 ENDIF 373 369 ! -
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14868 r14889 262 262 zdf_osm_alloc = zdf_osm_alloc + ierr 263 263 ! 264 ALLOCATE( etmean(A2D( 0),jpk), dh(jpi,jpj), r1_ft(A2D(0)), STAT=ierr )264 ALLOCATE( etmean(A2D((nn_hls-1)),jpk), dh(jpi,jpj), r1_ft(A2D((nn_hls-1))), STAT=ierr ) 265 265 zdf_osm_alloc = zdf_osm_alloc + ierr 266 266 ! 267 ALLOCATE( nbld(jpi,jpj), nmld(A2D( 0)), STAT=ierr )267 ALLOCATE( nbld(jpi,jpj), nmld(A2D((nn_hls-1))), STAT=ierr ) 268 268 zdf_osm_alloc = zdf_osm_alloc + ierr 269 269 ! 270 ALLOCATE( n_ddh(A2D( 0)), STAT=ierr )270 ALLOCATE( n_ddh(A2D((nn_hls-1))), STAT=ierr ) 271 271 zdf_osm_alloc = zdf_osm_alloc + ierr 272 272 ! 273 ALLOCATE( l_conv(A2D( 0)), l_shear(A2D(0)), l_coup(A2D(0)), l_pyc(A2D(0)), l_flux(A2D(0)), l_mle(A2D(0)), STAT=ierr )273 ALLOCATE( l_conv(A2D((nn_hls-1))), l_shear(A2D((nn_hls-1))), l_coup(A2D((nn_hls-1))), l_pyc(A2D((nn_hls-1))), l_flux(A2D((nn_hls-1))), l_mle(A2D((nn_hls-1))), STAT=ierr ) 274 274 zdf_osm_alloc = zdf_osm_alloc + ierr 275 275 ! 276 ALLOCATE( swth0(A2D( 0)), sws0(A2D(0)), swb0(A2D(0)), suw0(A2D(0)), sustar(A2D(0)), scos_wind(A2D(0)), &277 & ssin_wind(A2D( 0)), swthav(A2D(0)), swsav(A2D(0)), swbav(A2D(0)), sustke(A2D(0)), dstokes(A2D(0)), &278 & swstrl(A2D( 0)), swstrc(A2D(0)), sla(A2D(0)), svstr(A2D(0)), shol(A2D(0)), STAT=ierr )276 ALLOCATE( swth0(A2D((nn_hls-1))), sws0(A2D((nn_hls-1))), swb0(A2D((nn_hls-1))), suw0(A2D((nn_hls-1))), sustar(A2D((nn_hls-1))), scos_wind(A2D((nn_hls-1))), & 277 & ssin_wind(A2D((nn_hls-1))), swthav(A2D((nn_hls-1))), swsav(A2D((nn_hls-1))), swbav(A2D((nn_hls-1))), sustke(A2D((nn_hls-1))), dstokes(A2D((nn_hls-1))), & 278 & swstrl(A2D((nn_hls-1))), swstrc(A2D((nn_hls-1))), sla(A2D((nn_hls-1))), svstr(A2D((nn_hls-1))), shol(A2D((nn_hls-1))), STAT=ierr ) 279 279 zdf_osm_alloc = zdf_osm_alloc + ierr 280 280 ! 281 ALLOCATE( av_t_bl(A2D( 0)), av_s_bl(A2D(0)), av_u_bl(A2D(0)), av_v_bl(A2D(0)), av_b_bl(A2D(0)), STAT=ierr)281 ALLOCATE( av_t_bl(A2D((nn_hls-1))), av_s_bl(A2D((nn_hls-1))), av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), STAT=ierr) 282 282 zdf_osm_alloc = zdf_osm_alloc + ierr 283 283 ! 284 ALLOCATE( av_dt_bl(A2D( 0)), av_ds_bl(A2D(0)), av_du_bl(A2D(0)), av_dv_bl(A2D(0)), av_db_bl(A2D(0)), STAT=ierr)284 ALLOCATE( av_dt_bl(A2D((nn_hls-1))), av_ds_bl(A2D((nn_hls-1))), av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), STAT=ierr) 285 285 zdf_osm_alloc = zdf_osm_alloc + ierr 286 286 ! 287 ALLOCATE( av_t_ml(A2D( 0)), av_s_ml(A2D(0)), av_u_ml(A2D(0)), av_v_ml(A2D(0)), av_b_ml(A2D(0)), STAT=ierr)287 ALLOCATE( av_t_ml(A2D((nn_hls-1))), av_s_ml(A2D((nn_hls-1))), av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), STAT=ierr) 288 288 zdf_osm_alloc = zdf_osm_alloc + ierr 289 289 ! 290 ALLOCATE( av_dt_ml(A2D( 0)), av_ds_ml(A2D(0)), av_du_ml(A2D(0)), av_dv_ml(A2D(0)), av_db_ml(A2D(0)), STAT=ierr)290 ALLOCATE( av_dt_ml(A2D((nn_hls-1))), av_ds_ml(A2D((nn_hls-1))), av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), STAT=ierr) 291 291 zdf_osm_alloc = zdf_osm_alloc + ierr 292 292 ! 293 ALLOCATE( av_t_mle(A2D( 0)), av_s_mle(A2D(0)), av_u_mle(A2D(0)), av_v_mle(A2D(0)), av_b_mle(A2D(0)), STAT=ierr)293 ALLOCATE( av_t_mle(A2D((nn_hls-1))), av_s_mle(A2D((nn_hls-1))), av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))), av_b_mle(A2D((nn_hls-1))), STAT=ierr) 294 294 zdf_osm_alloc = zdf_osm_alloc + ierr 295 295 ! … … 351 351 ! 352 352 ! Scales 353 REAL(wp), DIMENSION(A2D( 0)) :: zrad0 ! Surface solar temperature flux (deg m/s)354 REAL(wp), DIMENSION(A2D( 0)) :: zradh ! Radiative flux at bl base (Buoyancy units)353 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zrad0 ! Surface solar temperature flux (deg m/s) 354 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zradh ! Radiative flux at bl base (Buoyancy units) 355 355 REAL(wp) :: zradav ! Radiative flux, bl average (Buoyancy Units) 356 356 REAL(wp) :: zvw0 ! Surface v-momentum flux 357 REAL(wp), DIMENSION(A2D( 0)) :: zwb0tot ! Total surface buoyancy flux including insolation358 REAL(wp), DIMENSION(A2D( 0)) :: zwb_ent ! Buoyancy entrainment flux359 REAL(wp), DIMENSION(A2D( 0)) :: zwb_min360 REAL(wp), DIMENSION(A2D( 0)) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL361 REAL(wp), DIMENSION(A2D( 0)) :: zwb_fk ! Max MLE buoyancy flux362 REAL(wp), DIMENSION(A2D( 0)) :: zdiff_mle ! Extra MLE vertical diff363 REAL(wp), DIMENSION(A2D( 0)) :: zvel_mle ! Velocity scale for dhdt with stable ML and FK357 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwb0tot ! Total surface buoyancy flux including insolation 358 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwb_ent ! Buoyancy entrainment flux 359 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwb_min 360 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 361 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwb_fk ! Max MLE buoyancy flux 362 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdiff_mle ! Extra MLE vertical diff 363 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zvel_mle ! Velocity scale for dhdt with stable ML and FK 364 364 ! 365 365 ! mixed-layer variables 366 INTEGER, DIMENSION(A2D( 0)) :: jk_ext ! Offset for external level367 ! 368 REAL(wp), DIMENSION(A2D( 0)) :: zhbl ! BL depth - grid369 REAL(wp), DIMENSION(A2D( 0)) :: zhml ! ML depth - grid370 ! 371 REAL(wp), DIMENSION(A2D( 0)) :: zhmle ! MLE depth - grid366 INTEGER, DIMENSION(A2D((nn_hls-1))) :: jk_ext ! Offset for external level 367 ! 368 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zhbl ! BL depth - grid 369 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zhml ! ML depth - grid 370 ! 371 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zhmle ! MLE depth - grid 372 372 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 373 373 ! 374 REAL(wp), DIMENSION(A2D( 0)) :: zdh ! Pycnocline depth - grid375 REAL(wp), DIMENSION(A2D( 0)) :: zdhdt ! BL depth tendency376 REAL(wp), DIMENSION(A2D( 0)) :: zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ! External temperature/salinity and buoyancy gradients374 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdh ! Pycnocline depth - grid 375 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdhdt ! BL depth tendency 376 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ! External temperature/salinity and buoyancy gradients 377 377 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! Horizontal gradients for Fox-Kemper parametrization 378 378 ! 379 REAL(wp), DIMENSION(A2D( 0)) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient379 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient 380 380 ! Flux-gradient relationship variables 381 REAL(wp), DIMENSION(A2D( 0)) :: zshear ! Shear production382 ! 383 REAL(wp), DIMENSION(A2D( 0)) :: zhbl_t ! Holds boundary layer depth updated by full timestep381 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zshear ! Shear production 382 ! 383 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zhbl_t ! Holds boundary layer depth updated by full timestep 384 384 ! 385 385 ! For calculating Ri#-dependent mixing … … 392 392 REAL(wp) :: zz0, zz1, zfac 393 393 REAL(wp) :: zus_x, zus_y ! Temporary Stokes drift 394 REAL(wp), DIMENSION(A2D( 0),jpk) :: zviscos ! Viscosity395 REAL(wp), DIMENSION(A2D( 0),jpk) :: zdiffut ! t-diffusivity394 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) :: zviscos ! Viscosity 395 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) :: zdiffut ! t-diffusivity 396 396 REAL(wp) :: zabsstke 397 397 REAL(wp) :: zsqrtpi, z_two_thirds, zthickness … … 423 423 ! 424 424 ghamt(:,:,:) = pp_large ; ghams(:,:,:) = pp_large 425 ghamt(A2D( 0),:) = 0.0_wp ; ghams(A2D(0),:) = 0.0_wp425 ghamt(A2D((nn_hls-1)),:) = 0.0_wp ; ghams(A2D((nn_hls-1)),:) = 0.0_wp 426 426 ghamu(:,:,:) = pp_large ; ghamv(:,:,:) = pp_large 427 ghamu(A2D( 0),:) = 0.0_wp ; ghamv(A2D(0),:) = 0.0_wp427 ghamu(A2D((nn_hls-1)),:) = 0.0_wp ; ghamv(A2D((nn_hls-1)),:) = 0.0_wp 428 428 ! 429 429 zdiff_mle(:,:) = 0.0_wp 430 430 ! 431 ! hbl = MAX(hbl,epsln) 431 ! Ensure only positive hbl values are accessed when using extended halo 432 ! (nn_hls==2) 433 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 434 hbl(ji,jj) = MAX( hbl(ji,jj), epsln ) 435 END_2D 436 ! 432 437 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 433 438 ! Calculate boundary layer scales … … 437 442 zz0 = rn_abs ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands 438 443 zz1 = 1.0_wp - rn_abs 439 DO_2D( 0, 0, 0, 0)444 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 440 445 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp ! Surface downward irradiance (so always +ve) 441 446 zradh(ji,jj) = zrad0(ji,jj) * & ! Downwards irradiance at base of boundary layer … … 448 453 & zradav ) ! over depth of OSBL 449 454 END_2D 450 DO_2D( 0, 0, 0, 0)455 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 451 456 sws0(ji,jj) = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + & ! Upwards surface salinity flux 452 457 & sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) ! for non-local term … … 459 464 & grav * zbeta * swsav(ji,jj) ! OBSBL 460 465 END_2D 461 DO_2D( 0, 0, 0, 0)466 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 462 467 suw0(ji,jj) = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) ! Surface upward velocity fluxes 463 468 zvw0 = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) … … 471 476 ! Assume constant La#=0.3 472 477 CASE(0) 473 DO_2D( 0, 0, 0, 0)478 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 474 479 zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 475 480 zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 … … 480 485 ! Assume Pierson-Moskovitz wind-wave spectrum 481 486 CASE(1) 482 DO_2D( 0, 0, 0, 0)487 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 483 488 ! Use wind speed wndm included in sbc_oce module 484 489 sustke(ji,jj) = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp ) … … 489 494 zfac = 2.0_wp * rpi / 16.0_wp 490 495 ! 491 DO_2D( 0, 0, 0, 0)496 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 492 497 IF ( hsw(ji,jj) > 1e-4_wp ) THEN 493 498 ! Use wave fields … … 506 511 IF (ln_zdfosm_ice_shelter) THEN 507 512 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 508 DO_2D( 0, 0, 0, 0)513 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 509 514 sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) 510 515 dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) … … 519 524 ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested. 520 525 IF(nn_osm_wave > 0) THEN 521 sustke(A2D( 0)) = rn_zdfosm_adjust_sd * sustke(A2D(0))526 sustke(A2D((nn_hls-1))) = rn_zdfosm_adjust_sd * sustke(A2D((nn_hls-1))) 522 527 END IF 523 528 CASE(1) … … 526 531 zsqrtpi = SQRT(rpi) 527 532 z_two_thirds = 2.0_wp / 3.0_wp 528 DO_2D( 0, 0, 0, 0)533 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 529 534 zthickness = rn_osm_hblfrac*hbl(ji,jj) 530 535 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) … … 540 545 ! Assumes approximate depth profile of SD from Breivik (2016) 541 546 zsqrtpi = SQRT(rpi) 542 DO_2D( 0, 0, 0, 0)547 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 543 548 zthickness = rn_osm_hblfrac*hbl(ji,jj) 544 549 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) … … 562 567 ! Langmuir velocity scale (swstrl), La # (sla) 563 568 ! Mixed scale (svstr), convective velocity scale (swstrc) 564 DO_2D( 0, 0, 0, 0)569 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 565 570 ! Langmuir velocity scale (swstrl), at T-point 566 571 swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird … … 596 601 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 597 602 nbld(:,:) = 4 598 DO_3D( 1, 1, 1, 1, 5, jpkm1 )603 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 ) 599 604 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 600 605 nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) … … 603 608 ! ########################################################################## 604 609 ! 605 DO_2D( 0, 0, 0, 0)610 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 606 611 zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) 607 612 nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) ) … … 612 617 ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere 613 618 jk_ext(:,:) = 1 ! ag 19/03 614 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D( 0)), av_t_bl, av_s_bl, &619 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 615 620 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 616 621 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 617 jk_ext(:,:) = nbld(A2D( 0)) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03622 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03 618 623 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 619 624 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & … … 632 637 ! Fox-Kemper Scheme 633 638 mld_prof = 4 634 DO_3D( 0, 0, 0, 0, 5, jpkm1 )639 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 635 640 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 636 641 END_3D 637 CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D( 0)), av_t_mle, av_s_mle, &642 CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D((nn_hls-1))), av_t_mle, av_s_mle, & 638 643 & av_b_mle, av_u_mle, av_v_mle ) 639 644 ! 640 DO_2D( 0, 0, 0, 0)645 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 641 646 zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 642 647 END_2D … … 656 661 l_flux(:,:) = .FALSE. 657 662 l_mle(:,:) = .FALSE. 658 DO_2D( 0, 0, 0, 0)663 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 659 664 IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 660 665 END_2D … … 662 667 ! 663 668 !! External gradient below BL needed both with and w/o FK 664 CALL zdf_osm_external_gradients( Kmm, nbld(A2D( 0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03669 CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03 665 670 ! 666 671 ! Test if pycnocline well resolved 667 ! DO_2D( 0, 0, 0, 0) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity672 ! DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 668 673 ! IF (l_conv(ji,jj) ) THEN should account for this. 669 674 ! ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm) … … 683 688 ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here.. 684 689 jk_ext(:,:) = 1 ! ag 19/03 685 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D( 0)), av_t_bl, av_s_bl, &690 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 686 691 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 687 692 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 688 jk_ext(:,:) = nbld(A2D( 0)) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03693 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03 689 694 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 690 695 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & … … 696 701 ! Test if surface boundary layer coupled to bottom 697 702 l_coup(:,:) = .FALSE. ! ag 19/03 698 DO_2D( 0, 0, 0, 0)703 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 699 704 zhbl_t(ji,jj) = hbl(ji,jj) + ( zdhdt(ji,jj) - ww(ji,jj,nbld(ji,jj)) ) * rn_Dt ! Certainly need ww here, so subtract it 700 705 ! Adjustment to represent limiting by ocean bottom … … 708 713 END_2D 709 714 ! 710 nmld(:,:) = nbld(A2D( 0)) ! use nmld to hold previous blayer index715 nmld(:,:) = nbld(A2D((nn_hls-1))) ! use nmld to hold previous blayer index 711 716 nbld(:,:) = 4 712 717 ! 713 DO_3D( 0, 0, 0, 0, 4, jpkm1 )718 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 ) 714 719 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 715 720 nbld(ji,jj) = jk … … 726 731 ! Recalculate BL averages and differences using new BL depth 727 732 jk_ext(:,:) = 1 ! ag 19/03 728 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D( 0)), av_t_bl, av_s_bl, &733 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 729 734 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 730 735 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) … … 734 739 ! 735 740 ! Reset l_pyc before calculating terms in the flux-gradient relationship 736 DO_2D( 0, 0, 0, 0)741 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 737 742 IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR. & 738 743 & nbld(ji,jj) - nmld(ji,jj) == 1 .OR. zdhdt(ji,jj) < 0.0_wp ) THEN ! ag 19/03 … … 748 753 END_2D 749 754 ! 750 dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D( 0))/ 3.0_wp ) ! Limit delta for shallow boundary layers for calculating755 dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D((nn_hls-1)))/ 3.0_wp ) ! Limit delta for shallow boundary layers for calculating 751 756 ! ! flux-gradient terms 752 757 ! … … 754 759 ! jk_ext = nbld - nmld + 1 755 760 ! Recalculate ML averages and differences using new ML depth 756 jk_ext(:,:) = nbld(A2D( 0)) - nmld(A2D(0)) + jk_ext(:,:) + 1 ! ag 19/03761 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1 ! ag 19/03 757 762 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 758 763 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & 759 764 & av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 760 765 ! 761 CALL zdf_osm_external_gradients( Kmm, nbld(A2D( 0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )766 CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 762 767 ! Rotate mean currents and changes onto wind aligned co-ordinates 763 768 CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml ) … … 786 791 ! 787 792 ! Rotate non-gradient velocity terms back to model reference frame 788 CALL zdf_osm_velocity_rotation( ghamu(A2D( 0),:), ghamv(A2D(0),:), .FALSE., 2, nbld(A2D(0)) )793 CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:), .FALSE., 2, nbld(A2D((nn_hls-1))) ) 789 794 ! 790 795 ! KPP-style Ri# mixing 791 796 IF ( ln_kpprimix ) THEN 792 797 jkflt = jpk 793 DO_2D( 0, 0, 0, 0)798 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 794 799 IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj) 795 800 END_2D 796 801 DO jk = jkflt+1, jpkm1 797 802 ! Shear production at uw- and vw-points (energy conserving form) 798 DO_2D( 1, 0, 1, 0)803 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 799 804 z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * & 800 805 & wumask(ji,jj,jk) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) … … 802 807 & wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 803 808 END_2D 804 DO_2D( 0, 0, 0, 0)809 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 805 810 IF ( jk > nbld(ji,jj) ) THEN 806 811 ! Shear prod. at w-point weightened by mask … … 821 826 ! KPP-style set diffusivity large if unstable below BL 822 827 IF ( ln_convmix) THEN 823 DO_2D( 0, 0, 0, 0)828 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 824 829 DO jk = nbld(ji,jj) + 1, jpkm1 825 830 IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) … … 829 834 ! 830 835 IF ( ln_osm_mle ) THEN ! Set up diffusivity and non-gradient mixing 831 DO_2D( 0, 0, 0, 0)836 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 832 837 IF ( l_flux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 833 838 ! Calculate MLE flux contribution from surface fluxes … … 862 867 ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 863 868 ! GN 25/8: need to change tmask --> wmask 864 DO_3D( 0, 0, 0, 0, 2, jpkm1 )869 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 865 870 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 866 871 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) … … 868 873 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and 869 874 ! v grids 870 CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp, & 871 & p_avm, 'W', 1.0_wp, & 872 & ghamu, 'W', 1.0_wp, & 873 & ghamv, 'W', 1.0_wp ) 875 IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp, & 876 & ghamv, 'W', 1.0_wp ) 874 877 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 875 878 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) * & … … 889 892 CASE(0:1) 890 893 IF ( iom_use("us_x") ) THEN ! x surface Stokes drift 891 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sustke * scos_wind894 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * scos_wind 892 895 CALL iom_put( "us_x", osmdia2d ) 893 896 END IF 894 897 IF ( iom_use("us_y") ) THEN ! y surface Stokes drift 895 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sustke * ssin_wind898 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * ssin_wind 896 899 CALL iom_put( "us_y", osmdia2d ) 897 900 END IF 898 901 IF ( iom_use("wind_wave_abs_power") ) THEN 899 osmdia2d(A2D( 0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke902 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 900 903 CALL iom_put( "wind_wave_abs_power", osmdia2d ) 901 904 END IF … … 913 916 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 914 917 IF ( iom_use("wind_wave_abs_power") ) THEN 915 osmdia2d(A2D( 0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * SQRT( ut0sd(A2D(0))**2 + vt0sd(A2D(0))**2 )918 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * SQRT( ut0sd(A2D((nn_hls-1)))**2 + vt0sd(A2D((nn_hls-1)))**2 ) 916 919 CALL iom_put( "wind_wave_abs_power", osmdia2d ) 917 920 END IF … … 922 925 IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL> 923 926 IF ( iom_use("zwth0") ) THEN ! <Tw_0> 924 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * swth0; CALL iom_put( "zwth0", osmdia2d )927 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0; CALL iom_put( "zwth0", osmdia2d ) 925 928 END IF 926 929 IF ( iom_use("zws0") ) THEN ! <Sw_0> 927 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sws0; CALL iom_put( "zws0", osmdia2d )930 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sws0; CALL iom_put( "zws0", osmdia2d ) 928 931 END IF 929 932 IF ( iom_use("zwb0") ) THEN ! <Sw_0> 930 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * swb0; CALL iom_put( "zwb0", osmdia2d )933 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swb0; CALL iom_put( "zwb0", osmdia2d ) 931 934 END IF 932 935 IF ( iom_use("zwbav") ) THEN ! Upward BL-avged turb buoyancy flux 933 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * swth0; CALL iom_put( "zwbav", osmdia2d )936 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0; CALL iom_put( "zwbav", osmdia2d ) 934 937 END IF 935 938 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! Boundary-layer depth 936 939 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld ) ! Boundary-layer max k 937 940 IF ( iom_use("zdt_bl") ) THEN ! dt at ml base 938 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_dt_bl; CALL iom_put( "zdt_bl", osmdia2d )941 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_bl; CALL iom_put( "zdt_bl", osmdia2d ) 939 942 END IF 940 943 IF ( iom_use("zds_bl") ) THEN ! ds at ml base 941 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_ds_bl; CALL iom_put( "zds_bl", osmdia2d )944 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_bl; CALL iom_put( "zds_bl", osmdia2d ) 942 945 END IF 943 946 IF ( iom_use("zdb_bl") ) THEN ! db at ml base 944 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_db_bl; CALL iom_put( "zdb_bl", osmdia2d )947 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_bl; CALL iom_put( "zdb_bl", osmdia2d ) 945 948 END IF 946 949 IF ( iom_use("zdu_bl") ) THEN ! du at ml base 947 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_du_bl; CALL iom_put( "zdu_bl", osmdia2d )950 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_du_bl; CALL iom_put( "zdu_bl", osmdia2d ) 948 951 END IF 949 952 IF ( iom_use("zdv_bl") ) THEN ! dv at ml base 950 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_dv_bl; CALL iom_put( "zdv_bl", osmdia2d )953 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dv_bl; CALL iom_put( "zdv_bl", osmdia2d ) 951 954 END IF 952 955 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 953 956 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 954 957 IF ( iom_use("zdt_ml") ) THEN ! dt at ml base 955 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_dt_ml; CALL iom_put( "zdt_ml", osmdia2d )958 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_ml; CALL iom_put( "zdt_ml", osmdia2d ) 956 959 END IF 957 960 IF ( iom_use("zds_ml") ) THEN ! ds at ml base 958 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_ds_ml; CALL iom_put( "zds_ml", osmdia2d )961 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_ml; CALL iom_put( "zds_ml", osmdia2d ) 959 962 END IF 960 963 IF ( iom_use("zdb_ml") ) THEN ! db at ml base 961 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_db_ml; CALL iom_put( "zdb_ml", osmdia2d )964 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_ml; CALL iom_put( "zdb_ml", osmdia2d ) 962 965 END IF 963 966 IF ( iom_use("dstokes") ) THEN ! Stokes drift penetration depth 964 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * dstokes; CALL iom_put( "dstokes", osmdia2d )967 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * dstokes; CALL iom_put( "dstokes", osmdia2d ) 965 968 END IF 966 969 IF ( iom_use("zustke") ) THEN ! Stokes drift magnitude at T-points 967 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sustke; CALL iom_put( "zustke", osmdia2d )970 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke; CALL iom_put( "zustke", osmdia2d ) 968 971 END IF 969 972 IF ( iom_use("zwstrc") ) THEN ! Convective velocity scale 970 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * swstrc; CALL iom_put( "zwstrc", osmdia2d )973 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrc; CALL iom_put( "zwstrc", osmdia2d ) 971 974 END IF 972 975 IF ( iom_use("zwstrl") ) THEN ! Langmuir velocity scale 973 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * swstrl; CALL iom_put( "zwstrl", osmdia2d )976 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrl; CALL iom_put( "zwstrl", osmdia2d ) 974 977 END IF 975 978 IF ( iom_use("zustar") ) THEN ! Friction velocity scale 976 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sustar; CALL iom_put( "zustar", osmdia2d )979 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustar; CALL iom_put( "zustar", osmdia2d ) 977 980 END IF 978 981 IF ( iom_use("zvstr") ) THEN ! Mixed velocity scale 979 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * svstr; CALL iom_put( "zvstr", osmdia2d )982 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * svstr; CALL iom_put( "zvstr", osmdia2d ) 980 983 END IF 981 984 IF ( iom_use("zla") ) THEN ! Langmuir # 982 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * sla; CALL iom_put( "zla", osmdia2d )985 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sla; CALL iom_put( "zla", osmdia2d ) 983 986 END IF 984 987 IF ( iom_use("wind_power") ) THEN ! BL depth internal to zdf_osm routine 985 osmdia2d(A2D( 0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**3988 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**3 986 989 CALL iom_put( "wind_power", osmdia2d ) 987 990 END IF 988 991 IF ( iom_use("wind_wave_power") ) THEN 989 osmdia2d(A2D( 0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke992 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 990 993 CALL iom_put( "wind_wave_power", osmdia2d ) 991 994 END IF 992 995 IF ( iom_use("zhbl") ) THEN ! BL depth internal to zdf_osm routine 993 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zhbl; CALL iom_put( "zhbl", osmdia2d )996 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhbl; CALL iom_put( "zhbl", osmdia2d ) 994 997 END IF 995 998 IF ( iom_use("zhml") ) THEN ! ML depth internal to zdf_osm routine 996 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zhml; CALL iom_put( "zhml", osmdia2d )999 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhml; CALL iom_put( "zhml", osmdia2d ) 997 1000 END IF 998 1001 IF ( iom_use("imld") ) THEN ! Index for ML depth internal to zdf_osm routine 999 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * nmld; CALL iom_put( "imld", osmdia2d )1002 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * nmld; CALL iom_put( "imld", osmdia2d ) 1000 1003 END IF 1001 1004 IF ( iom_use("jp_ext") ) THEN ! =1 if pycnocline resolved internal to zdf_osm routine 1002 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * jk_ext; CALL iom_put( "jp_ext", osmdia2d )1005 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * jk_ext; CALL iom_put( "jp_ext", osmdia2d ) 1003 1006 END IF 1004 1007 IF ( iom_use("j_ddh") ) THEN ! Index forpyc thicknessh internal to zdf_osm routine 1005 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * n_ddh; CALL iom_put( "j_ddh", osmdia2d )1008 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * n_ddh; CALL iom_put( "j_ddh", osmdia2d ) 1006 1009 END IF 1007 1010 IF ( iom_use("zshear") ) THEN ! Shear production of TKE internal to zdf_osm routine 1008 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zshear; CALL iom_put( "zshear", osmdia2d )1011 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zshear; CALL iom_put( "zshear", osmdia2d ) 1009 1012 END IF 1010 1013 IF ( iom_use("zdh") ) THEN ! Pyc thicknessh internal to zdf_osm routine 1011 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zdh; CALL iom_put( "zdh", osmdia2d )1014 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdh; CALL iom_put( "zdh", osmdia2d ) 1012 1015 END IF 1013 1016 IF ( iom_use("zhol") ) THEN ! ML depth internal to zdf_osm routine 1014 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * shol; CALL iom_put( "zhol", osmdia2d )1017 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * shol; CALL iom_put( "zhol", osmdia2d ) 1015 1018 END IF 1016 1019 IF ( iom_use("zwb_ent") ) THEN ! Upward turb buoyancy entrainment flux 1017 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zwb_ent; CALL iom_put( "zwb_ent", osmdia2d )1020 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_ent; CALL iom_put( "zwb_ent", osmdia2d ) 1018 1021 END IF 1019 1022 IF ( iom_use("zt_ml") ) THEN ! Average T in ML 1020 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * av_t_ml; CALL iom_put( "zt_ml", osmdia2d )1023 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_t_ml; CALL iom_put( "zt_ml", osmdia2d ) 1021 1024 END IF 1022 1025 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1023 1026 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1024 1027 IF ( iom_use("zwb_fk") ) THEN ! FK b flux 1025 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zwb_fk; CALL iom_put( "zwb_fk", osmdia2d )1028 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk; CALL iom_put( "zwb_fk", osmdia2d ) 1026 1029 END IF 1027 1030 IF ( iom_use("zwb_fk_b") ) THEN ! FK b flux averaged over ML 1028 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zwb_fk_b; CALL iom_put( "zwb_fk_b", osmdia2d )1031 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk_b; CALL iom_put( "zwb_fk_b", osmdia2d ) 1029 1032 END IF 1030 1033 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof ) ! FK layer max k … … 1036 1039 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1037 1040 IF ( iom_use("zdiff_mle") ) THEN ! FK diff in MLE at t-pt 1038 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d )1041 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d ) 1039 1042 END IF 1040 1043 IF ( iom_use("zvel_mle") ) THEN ! FK diff in MLE at t-pt 1041 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zvel_mle", osmdia2d )1044 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zvel_mle", osmdia2d ) 1042 1045 END IF 1043 1046 END IF … … 1060 1063 !!---------------------------------------------------------------------- 1061 1064 INTEGER, INTENT(in ) :: Kbb, Kmm ! Ocean time-level indices 1062 INTEGER, DIMENSION(A2D( 0)), INTENT(in ) :: knlev ! Number of levels to average over.1063 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pt, ps ! Average temperature and salinity1064 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pb ! Average buoyancy1065 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pu, pv ! Average current components1066 INTEGER, DIMENSION(A2D( 0)), INTENT(in ), OPTIONAL :: kp_ext ! External-level offsets1067 REAL(wp), DIMENSION(A2D( 0)), INTENT( out), OPTIONAL :: pdt, pds ! Difference between average temperature, salinity,1068 REAL(wp), DIMENSION(A2D( 0)), INTENT( out), OPTIONAL :: pdb ! buoyancy,1069 REAL(wp), DIMENSION(A2D( 0)), INTENT( out), OPTIONAL :: pdu, pdv ! velocity components and the OSBL1065 INTEGER, DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: knlev ! Number of levels to average over. 1066 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pt, ps ! Average temperature and salinity 1067 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pb ! Average buoyancy 1068 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pu, pv ! Average current components 1069 INTEGER, DIMENSION(A2D((nn_hls-1))), INTENT(in ), OPTIONAL :: kp_ext ! External-level offsets 1070 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out), OPTIONAL :: pdt, pds ! Difference between average temperature, salinity, 1071 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out), OPTIONAL :: pdb ! buoyancy, 1072 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out), OPTIONAL :: pdu, pdv ! velocity components and the OSBL 1070 1073 ! 1071 1074 INTEGER :: jk, jkflt, jkmax, ji, jj ! Loop indices 1072 1075 INTEGER :: ibld_ext ! External-layer index 1073 REAL(wp), DIMENSION(A2D( 0)) :: zthick ! Layer thickness1076 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zthick ! Layer thickness 1074 1077 REAL(wp) :: zthermal, zbeta ! Thermal/haline expansion/contraction coefficients 1075 1078 !!---------------------------------------------------------------------- … … 1083 1086 jkflt = jpk 1084 1087 jkmax = 0 1085 DO_2D( 0, 0, 0, 0)1088 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1086 1089 IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj) 1087 1090 IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 1088 1091 END_2D 1089 DO_3D( 0, 0, 0, 0, 2, jkflt ) ! Upper, flat part of layer1092 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt ) ! Upper, flat part of layer 1090 1093 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 1091 1094 pt(ji,jj) = pt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) … … 1098 1101 & MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1099 1102 END_3D 1100 DO_3D( 0, 0, 0, 0, jkflt+1, jkmax ) ! Lower, non-flat part of layer1103 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax ) ! Lower, non-flat part of layer 1101 1104 IF ( knlev(ji,jj) >= jk ) THEN 1102 1105 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) … … 1111 1114 END IF 1112 1115 END_3D 1113 DO_2D( 0, 0, 0, 0)1116 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1114 1117 pt(ji,jj) = pt(ji,jj) / zthick(ji,jj) 1115 1118 ps(ji,jj) = ps(ji,jj) / zthick(ji,jj) … … 1123 1126 ! Differences between vertical averages and values at an external layer 1124 1127 IF ( PRESENT( kp_ext ) ) THEN 1125 DO_2D( 0, 0, 0, 0)1128 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1126 1129 ibld_ext = knlev(ji,jj) + kp_ext(ji,jj) 1127 1130 IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN ! ag 09/03 … … 1160 1163 !! 1161 1164 !!---------------------------------------------------------------------- 1162 REAL(wp), INTENT(inout), DIMENSION(A2D( 0)) :: pu, pv ! Components of current1165 REAL(wp), INTENT(inout), DIMENSION(A2D((nn_hls-1))) :: pu, pv ! Components of current 1163 1166 LOGICAL, OPTIONAL, INTENT(in ) :: fwd ! Forward (default) or reverse rotation 1164 1167 ! … … 1168 1171 zfwd = 1.0_wp 1169 1172 IF( PRESENT(fwd) .AND. ( .NOT. fwd ) ) zfwd = -1.0_wp 1170 DO_2D( 0, 0, 0, 0)1173 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1171 1174 ztmp = pu(ji,jj) 1172 1175 pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj) … … 1190 1193 !! 1191 1194 !!---------------------------------------------------------------------- 1192 REAL(wp), INTENT(inout), DIMENSION(A2D( 0),jpk) :: pu, pv ! Components of current1195 REAL(wp), INTENT(inout), DIMENSION(A2D((nn_hls-1)),jpk) :: pu, pv ! Components of current 1193 1196 LOGICAL, OPTIONAL, INTENT(in ) :: fwd ! Forward (default) or reverse rotation 1194 1197 INTEGER, OPTIONAL, INTENT(in ) :: ktop ! Minimum depth index 1195 INTEGER, OPTIONAL, INTENT(in ), DIMENSION(A2D( 0)) :: knlev ! Array of maximum depth indices1198 INTEGER, OPTIONAL, INTENT(in ), DIMENSION(A2D((nn_hls-1))) :: knlev ! Array of maximum depth indices 1196 1199 ! 1197 1200 INTEGER :: ji, jj, jk, jktop, jkmax ! Loop indices … … 1205 1208 IF( PRESENT(knlev) ) THEN 1206 1209 jkmax = 0 1207 DO_2D( 0, 0, 0, 0)1210 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1208 1211 IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 1209 1212 END_2D … … 1213 1216 llkbot = .TRUE. 1214 1217 END IF 1215 DO_3D( 0, 0, 0, 0, jktop, jkmax )1218 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax ) 1216 1219 IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN 1217 1220 ztmp = pu(ji,jj,jk) … … 1237 1240 !!---------------------------------------------------------------------- 1238 1241 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 1239 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pwb_ent ! Buoyancy fluxes at base1240 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pwb_min ! of well-mixed layer1241 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pshear ! Production of TKE due to shear across the pycnocline1242 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth1243 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phml ! ML depth1244 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdh ! Pycnocline depth1242 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pwb_ent ! Buoyancy fluxes at base 1243 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pwb_min ! of well-mixed layer 1244 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pshear ! Production of TKE due to shear across the pycnocline 1245 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 1246 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phml ! ML depth 1247 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdh ! Pycnocline depth 1245 1248 ! 1246 1249 ! Local variables 1247 1250 INTEGER :: jj, ji ! Loop indices 1248 1251 ! 1249 REAL(wp), DIMENSION(A2D( 0)) :: zekman1250 REAL(wp), DIMENSION(A2D( 0)) :: zri_p, zri_b ! Richardson numbers1252 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zekman 1253 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zri_p, zri_b ! Richardson numbers 1251 1254 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1252 1255 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes … … 1269 1272 ! 1270 1273 ! Determins stability and set flag l_conv 1271 DO_2D( 0, 0, 0, 0)1274 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1272 1275 IF ( shol(ji,jj) < 0.0_wp ) THEN 1273 1276 l_conv(ji,jj) = .TRUE. … … 1278 1281 ! 1279 1282 pshear(:,:) = 0.0_wp 1280 zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D( 0)) ) * phbl(:,:) / MAX( sustar(A2D(0)), 1.e-8 ) )1281 ! 1282 DO_2D( 0, 0, 0, 0)1283 zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D((nn_hls-1))) ) * phbl(:,:) / MAX( sustar(A2D((nn_hls-1))), 1.e-8 ) ) 1284 ! 1285 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1283 1286 IF ( l_conv(ji,jj) ) THEN 1284 1287 IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN … … 1332 1335 ! 1333 1336 ! Calculate entrainment buoyancy flux due to surface fluxes. 1334 DO_2D( 0, 0, 0, 0)1337 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1335 1338 IF ( l_conv(ji,jj) ) THEN 1336 1339 zwcor = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln … … 1353 1356 END_2D 1354 1357 ! 1355 DO_2D( 0, 0, 0, 0)1358 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1356 1359 IF ( l_shear(ji,jj) ) THEN 1357 1360 IF ( l_conv(ji,jj) ) THEN … … 1392 1395 !!---------------------------------------------------------------------- 1393 1396 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 1394 INTEGER, DIMENSION(A2D( 0)), INTENT(in ) :: kbase ! OSBL base layer index1395 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pdtdz, pdsdz, pdbdz ! External gradients of temperature, salinity and buoyancy1397 INTEGER, DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: kbase ! OSBL base layer index 1398 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pdtdz, pdsdz, pdbdz ! External gradients of temperature, salinity and buoyancy 1396 1399 ! 1397 1400 ! Local variables … … 1405 1408 pdbdz(:,:) = pp_large 1406 1409 ! 1407 DO_2D( 0, 0, 0, 0)1410 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1408 1411 IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1409 1412 zthermal = rab_n(ji,jj,1,jp_tem) ! Ideally use nbld not 1?? … … 1433 1436 !! 1434 1437 !!---------------------------------------------------------------------- 1435 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pdhdt ! Rate of change of hbl1436 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth1437 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdh ! Pycnocline depth1438 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux1439 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_min1440 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients1441 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL1442 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_fk ! Max MLE buoyancy flux1443 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pvel_mle ! Vvelocity scale for dhdt with stable ML and FK1438 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pdhdt ! Rate of change of hbl 1439 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 1440 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdh ! Pycnocline depth 1441 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 1442 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_min 1443 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 1444 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL 1445 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_fk ! Max MLE buoyancy flux 1446 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pvel_mle ! Vvelocity scale for dhdt with stable ML and FK 1444 1447 ! 1445 1448 ! Local variables … … 1455 1458 pwb_fk_b(:,:) = pp_large 1456 1459 ! 1457 DO_2D( 0, 0, 0, 0)1460 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1458 1461 ! 1459 1462 IF ( l_shear(ji,jj) ) THEN … … 1598 1601 !!---------------------------------------------------------------------- 1599 1602 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 1600 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pdhdt ! Rates of change of hbl1601 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: phbl ! BL depth1602 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl_t ! BL depth1603 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux1604 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL1603 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdhdt ! Rates of change of hbl 1604 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: phbl ! BL depth 1605 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl_t ! BL depth 1606 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 1607 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL 1605 1608 ! 1606 1609 ! Local variables … … 1609 1612 REAL(wp) :: zthermal, zbeta 1610 1613 ! 1611 DO_2D( 0, 0, 0, 0)1614 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1612 1615 IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN 1613 1616 ! … … 1698 1701 !!---------------------------------------------------------------------- 1699 1702 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 1700 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pdh ! Pycnocline thickness1701 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: phml ! ML depth1702 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdhdt ! BL depth tendency1703 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth1704 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux1705 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients1706 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL1703 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdh ! Pycnocline thickness 1704 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: phml ! ML depth 1705 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdhdt ! BL depth tendency 1706 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 1707 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 1708 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 1709 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_fk_b ! MLE buoyancy flux averaged over OSBL 1707 1710 1708 1711 ! … … 1715 1718 REAL, PARAMETER :: pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp ! Also in pycnocline_depth 1716 1719 ! 1717 DO_2D( 0, 0, 0, 0)1720 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1718 1721 ! 1719 1722 IF ( l_shear(ji,jj) ) THEN … … 1859 1862 !!---------------------------------------------------------------------- 1860 1863 INTEGER, INTENT(in ) :: Kmm ! Ocean time-level index 1861 INTEGER, DIMENSION(A2D( 0)), INTENT(in ) :: kp_ext ! External-level offsets1862 REAL(wp), DIMENSION(A2D( 0),jpk), INTENT( out) :: pdbdz ! Gradients in the pycnocline1863 REAL(wp), DIMENSION(A2D( 0)), INTENT( out) :: palpha1864 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdh ! Pycnocline thickness1865 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth1866 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients1867 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phml ! ML depth1868 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdhdt ! Rates of change of hbl1864 INTEGER, DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: kp_ext ! External-level offsets 1865 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk), INTENT( out) :: pdbdz ! Gradients in the pycnocline 1866 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT( out) :: palpha 1867 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdh ! Pycnocline thickness 1868 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 1869 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 1870 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phml ! ML depth 1871 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdhdt ! Rates of change of hbl 1869 1872 ! 1870 1873 ! Local variables … … 1881 1884 palpha(:,:) = pp_large 1882 1885 ! 1883 DO_2D( 0, 0, 0, 0)1886 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1884 1887 ! 1885 1888 IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 1950 1953 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 1951 1954 IF ( iom_use("zdbdz_pyc") ) THEN 1952 osmdia3d(A2D( 0),:) = wmask(A2D(0),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d )1955 osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d ) 1953 1956 END IF 1954 1957 END IF … … 1969 1972 !!---------------------------------------------------------------------- 1970 1973 INTEGER, INTENT(in ) :: Kbb, Kmm ! Ocean time-level indices 1971 REAL(wp), DIMENSION(A2D( 0),jpk), INTENT(inout) :: pdiffut ! t-diffusivity1972 REAL(wp), DIMENSION(A2D( 0),jpk), INTENT(inout) :: pviscos ! Viscosity1973 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth1974 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phml ! ML depth1975 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdh ! Pycnocline depth1976 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdhdt ! BL depth tendency1977 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pshear ! Shear production1978 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux1979 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_min1974 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk), INTENT(inout) :: pdiffut ! t-diffusivity 1975 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk), INTENT(inout) :: pviscos ! Viscosity 1976 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 1977 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phml ! ML depth 1978 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdh ! Pycnocline depth 1979 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdhdt ! BL depth tendency 1980 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pshear ! Shear production 1981 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 1982 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_min 1980 1983 ! 1981 1984 ! Local variables … … 1983 1986 ! 1984 1987 ! Scales used to calculate eddy diffusivity and viscosity profiles 1985 REAL(wp), DIMENSION(A2D( 0)) :: zdifml_sc, zvisml_sc1986 REAL(wp), DIMENSION(A2D( 0)) :: zdifpyc_n_sc, zdifpyc_s_sc1987 REAL(wp), DIMENSION(A2D( 0)) :: zvispyc_n_sc, zvispyc_s_sc1988 REAL(wp), DIMENSION(A2D( 0)) :: zbeta_d_sc, zbeta_v_sc1989 REAL(wp), DIMENSION(A2D( 0)) :: zb_coup, zc_coup_vis, zc_coup_dif1988 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdifml_sc, zvisml_sc 1989 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdifpyc_n_sc, zdifpyc_s_sc 1990 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zvispyc_n_sc, zvispyc_s_sc 1991 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zbeta_d_sc, zbeta_v_sc 1992 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zb_coup, zc_coup_vis, zc_coup_dif 1990 1993 ! 1991 1994 REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b … … 2001 2004 zb_coup(:,:) = 0.0_wp 2002 2005 ! 2003 DO_2D( 0, 0, 0, 0)2006 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2004 2007 IF ( l_conv(ji,jj) ) THEN 2005 2008 ! … … 2080 2083 END_2D 2081 2084 ! 2082 DO_2D( 0, 0, 0, 0)2085 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2083 2086 IF ( l_conv(ji,jj) ) THEN 2084 2087 DO jk = 2, nmld(ji,jj) ! Mixed layer diffusivity … … 2150 2153 END_2D 2151 2154 IF( iom_use("pb_coup") ) THEN 2152 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d ) ! BBL-coupling velocity scale2155 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d ) ! BBL-coupling velocity scale 2153 2156 END IF 2154 2157 ! … … 2167 2170 !!---------------------------------------------------------------------- 2168 2171 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2169 INTEGER, DIMENSION(A2D( 0)), INTENT(in ) :: kp_ext ! Offset for external level2170 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth2171 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phml ! ML depth2172 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdh ! Pycnocline depth2173 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdhdt ! BL depth tendency2174 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pshear ! Shear production2175 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdtdz_bl_ext ! External temperature gradients2176 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdsdz_bl_ext ! External salinity gradients2177 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients2178 REAL(wp), DIMENSION(A2D( 0),jpk), INTENT(in ) :: pdiffut ! t-diffusivity2179 REAL(wp), DIMENSION(A2D( 0),jpk), INTENT(in ) :: pviscos ! Viscosity2180 ! 2181 REAL(wp), DIMENSION(A2D( 0)) :: zalpha_pyc !2182 REAL(wp), DIMENSION(A2D( 0),jpk) :: zdbdz_pyc ! Parametrised gradient of buoyancy in the pycnocline2172 INTEGER, DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: kp_ext ! Offset for external level 2173 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 2174 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phml ! ML depth 2175 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdh ! Pycnocline depth 2176 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdhdt ! BL depth tendency 2177 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pshear ! Shear production 2178 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdtdz_bl_ext ! External temperature gradients 2179 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdsdz_bl_ext ! External salinity gradients 2180 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbdz_bl_ext ! External buoyancy gradients 2181 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk), INTENT(in ) :: pdiffut ! t-diffusivity 2182 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk), INTENT(in ) :: pviscos ! Viscosity 2183 ! 2184 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zalpha_pyc ! 2185 REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) :: zdbdz_pyc ! Parametrised gradient of buoyancy in the pycnocline 2183 2186 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z3ddz_pyc_1, z3ddz_pyc_2 ! Pycnocline gradient/shear profiles 2184 2187 ! … … 2186 2189 INTEGER :: istat ! Memory allocation status 2187 2190 REAL(wp) :: zznd_d, zznd_ml, zznd_pyc, znd ! Temporary non-dimensional depths 2188 REAL(wp), DIMENSION(A2D( 0)) :: zsc_wth_1,zsc_ws_1 ! Temporary scales2189 REAL(wp), DIMENSION(A2D( 0)) :: zsc_uw_1, zsc_uw_2 ! Temporary scales2190 REAL(wp), DIMENSION(A2D( 0)) :: zsc_vw_1, zsc_vw_2 ! Temporary scales2191 REAL(wp), DIMENSION(A2D( 0)) :: ztau_sc_u ! Dissipation timescale at base of WML2191 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zsc_wth_1,zsc_ws_1 ! Temporary scales 2192 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zsc_uw_1, zsc_uw_2 ! Temporary scales 2193 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zsc_vw_1, zsc_vw_2 ! Temporary scales 2194 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: ztau_sc_u ! Dissipation timescale at base of WML 2192 2195 REAL(wp) :: zbuoy_pyc_sc, zdelta_pyc ! 2193 2196 REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale 2194 REAL(wp), DIMENSION(A2D( 0)) :: za_cubic, zb_cubic ! Coefficients in cubic polynomial specifying2195 REAL(wp), DIMENSION(A2D( 0)) :: zc_cubic, zd_cubic ! diffusivity in pycnocline2196 REAL(wp), DIMENSION(A2D( 0)) :: zwt_pyc_sc_1, zws_pyc_sc_1 !2197 REAL(wp), DIMENSION(A2D( 0)) :: zzeta_pyc !2197 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: za_cubic, zb_cubic ! Coefficients in cubic polynomial specifying 2198 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zc_cubic, zd_cubic ! diffusivity in pycnocline 2199 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwt_pyc_sc_1, zws_pyc_sc_1 ! 2200 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zzeta_pyc ! 2198 2201 REAL(wp) :: zomega, zvw_max ! 2199 REAL(wp), DIMENSION(A2D( 0)) :: zuw_bse,zvw_bse ! Momentum, heat, and salinity fluxes2200 REAL(wp), DIMENSION(A2D( 0)) :: zwth_ent,zws_ent ! at the top of the pycnocline2201 REAL(wp), DIMENSION(A2D( 0)) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term2202 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zuw_bse,zvw_bse ! Momentum, heat, and salinity fluxes 2203 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zwth_ent,zws_ent ! at the top of the pycnocline 2204 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term 2202 2205 REAL(wp) :: ztmp ! 2203 2206 REAL(wp) :: ztgrad, zsgrad, zbgrad ! Variables used to calculate pycnocline gradients … … 2220 2223 jkf_mld = jpk 2221 2224 jkm_mld = 0 2222 DO_2D( 0, 0, 0, 0)2225 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2223 2226 IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj) 2224 2227 IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj) … … 2228 2231 ! Stokes term in scalar flux, flux-gradient relationship 2229 2232 ! ------------------------------------------------------ 2230 WHERE ( l_conv(A2D( 0)) )2231 zsc_wth_1(:,:) = swstrl(A2D( 0))**3 * swth0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )2232 zsc_ws_1(:,:) = swstrl(A2D( 0))**3 * sws0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )2233 WHERE ( l_conv(A2D((nn_hls-1))) ) 2234 zsc_wth_1(:,:) = swstrl(A2D((nn_hls-1)))**3 * swth0(A2D((nn_hls-1))) / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 2235 zsc_ws_1(:,:) = swstrl(A2D((nn_hls-1)))**3 * sws0(A2D((nn_hls-1))) / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 2233 2236 ELSEWHERE 2234 zsc_wth_1(:,:) = 2.0_wp * swthav(A2D( 0))2235 zsc_ws_1(:,:) = 2.0_wp * swsav(A2D( 0))2237 zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1))) 2238 zsc_ws_1(:,:) = 2.0_wp * swsav(A2D((nn_hls-1))) 2236 2239 ENDWHERE 2237 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )2240 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2238 2241 IF ( l_conv(ji,jj) ) THEN 2239 2242 IF ( jk <= nmld(ji,jj) ) THEN … … 2263 2266 ! svstr since term needs to go to zero as swstrl goes to zero) 2264 2267 ! --------------------------------------------------------------------- 2265 WHERE ( l_conv(A2D( 0)) )2266 zsc_uw_1(:,:) = ( swstrl(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) / &2267 & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D( 0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )2268 zsc_uw_2(:,:) = ( swstrl(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) / &2269 & MIN( sla(A2D( 0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )2270 zsc_vw_1(:,:) = ff_t(A2D( 0)) * phml(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / &2271 & ( ( svstr(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln )2268 WHERE ( l_conv(A2D((nn_hls-1))) ) 2269 zsc_uw_1(:,:) = ( swstrl(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) / & 2270 & MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 2271 zsc_uw_2(:,:) = ( swstrl(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) / & 2272 & MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 2273 zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 * MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2274 & ( ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 2272 2275 ELSEWHERE 2273 zsc_uw_1(:,:) = sustar(A2D( 0))**22274 zsc_vw_1(:,:) = ff_t(A2D( 0)) * phbl(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / &2275 & ( svstr(A2D( 0))**2 + epsln )2276 zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 2277 zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 * MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / & 2278 & ( svstr(A2D((nn_hls-1)))**2 + epsln ) 2276 2279 ENDWHERE 2277 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )2280 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2278 2281 IF ( l_conv(ji,jj) ) THEN 2279 2282 IF ( jk <= nmld(ji,jj) ) THEN … … 2297 2300 ! (X0.3) and pressure (X0.5)] 2298 2301 ! ---------------------------------------------------------------------- 2299 WHERE ( l_conv(A2D( 0)) )2300 zsc_wth_1(:,:) = swbav(A2D( 0)) * swth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) / &2301 & ( svstr(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )2302 zsc_ws_1(:,:) = swbav(A2D( 0)) * sws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) / &2303 & ( svstr(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )2302 WHERE ( l_conv(A2D((nn_hls-1))) ) 2303 zsc_wth_1(:,:) = swbav(A2D((nn_hls-1))) * swth0(A2D((nn_hls-1))) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) * phml(A2D((nn_hls-1))) / & 2304 & ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 2305 zsc_ws_1(:,:) = swbav(A2D((nn_hls-1))) * sws0(A2D((nn_hls-1))) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) * phml(A2D((nn_hls-1))) / & 2306 & ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 2304 2307 ELSEWHERE 2305 2308 zsc_wth_1(:,:) = 0.0_wp 2306 2309 zsc_ws_1(:,:) = 0.0_wp 2307 2310 ENDWHERE 2308 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )2311 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2309 2312 IF ( l_conv(ji,jj) ) THEN 2310 2313 IF ( jk <= nmld(ji,jj) ) THEN … … 2327 2330 END IF 2328 2331 END_3D 2329 DO_2D( 0, 0, 0, 0)2332 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2330 2333 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 2331 2334 ztau_sc_u(ji,jj) = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird * & … … 2347 2350 END IF 2348 2351 END_2D 2349 DO_3D( 0, 0, 0, 0, 2, jkm_bld )2352 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2350 2353 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2351 2354 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) … … 2369 2372 IF ( ln_dia_osm ) THEN 2370 2373 IF ( iom_use("zwth_ent") ) THEN ! Upward turb. temperature entrainment flux 2371 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d )2374 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d ) 2372 2375 END IF 2373 2376 IF ( iom_use("zws_ent") ) THEN ! Upward turb. salinity entrainment flux 2374 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zws_ent(:,:); CALL iom_put( "zws_ent", osmdia2d )2377 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zws_ent(:,:); CALL iom_put( "zws_ent", osmdia2d ) 2375 2378 END IF 2376 2379 END IF 2377 2380 ! 2378 2381 zsc_vw_1(:,:) = 0.0_wp 2379 WHERE ( l_conv(A2D( 0)) )2380 zsc_uw_1(:,:) = -1.0_wp * swb0(A2D( 0)) * sustar(A2D(0))**2 * phml(A2D(0)) / &2381 & ( svstr(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )2382 zsc_uw_2(:,:) = swb0(A2D( 0)) * sustke(A2D(0)) * phml(A2D(0)) / &2383 & ( svstr(A2D( 0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp )2382 WHERE ( l_conv(A2D((nn_hls-1))) ) 2383 zsc_uw_1(:,:) = -1.0_wp * swb0(A2D((nn_hls-1))) * sustar(A2D((nn_hls-1)))**2 * phml(A2D((nn_hls-1))) / & 2384 & ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 2385 zsc_uw_2(:,:) = swb0(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) / & 2386 & ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )**( 2.0_wp / 3.0_wp ) 2384 2387 ELSEWHERE 2385 2388 zsc_uw_1(:,:) = 0.0_wp 2386 2389 ENDWHERE 2387 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )2390 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2388 2391 IF ( l_conv(ji,jj) ) THEN 2389 2392 IF ( jk <= nmld(ji,jj) ) THEN … … 2402 2405 END_3D 2403 2406 ! 2404 DO_2D( 0, 0, 0, 0)2407 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2405 2408 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 2406 2409 IF ( n_ddh(ji,jj) == 0 ) THEN … … 2421 2424 END IF 2422 2425 END_2D 2423 DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld ) ! Need ztau_sc_u to be available. Change to array.2426 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkf_mld, jkm_bld ) ! Need ztau_sc_u to be available. Change to array. 2424 2427 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2425 2428 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) … … 2436 2439 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask * ghamu ) 2437 2440 IF ( iom_use("zsc_uw_1_0") ) THEN 2438 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d )2441 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d ) 2439 2442 END IF 2440 2443 END IF … … 2443 2446 ! (X0.3) ] 2444 2447 ! ----------------------------------------------------------------------- 2445 WHERE ( l_conv(A2D( 0)) )2446 zsc_wth_1(:,:) = swth0(A2D( 0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) )2447 zsc_ws_1(:,:) = sws0(A2D( 0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) )2448 WHERE ( l_pyc(A2D( 0)) ) ! Pycnocline scales2449 zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D( 0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_dt_ml(A2D(0))2450 zsc_ws_pyc(:,:) = -0.003_wp * swstrc(A2D( 0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_ds_ml(A2D(0))2448 WHERE ( l_conv(A2D((nn_hls-1))) ) 2449 zsc_wth_1(:,:) = swth0(A2D((nn_hls-1))) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) ) 2450 zsc_ws_1(:,:) = sws0(A2D((nn_hls-1))) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) ) 2451 WHERE ( l_pyc(A2D((nn_hls-1))) ) ! Pycnocline scales 2452 zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) * av_dt_ml(A2D((nn_hls-1))) 2453 zsc_ws_pyc(:,:) = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) * av_ds_ml(A2D((nn_hls-1))) 2451 2454 END WHERE 2452 2455 ELSEWHERE 2453 zsc_wth_1(:,:) = 2.0_wp * swthav(A2D( 0))2454 zsc_ws_1(:,:) = sws0(A2D( 0))2456 zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1))) 2457 zsc_ws_1(:,:) = sws0(A2D((nn_hls-1))) 2455 2458 END WHERE 2456 DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) )2459 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) ) 2457 2460 IF ( l_conv(ji,jj) ) THEN 2458 2461 IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN … … 2484 2487 END_3D 2485 2488 ! 2486 WHERE ( l_conv(A2D( 0)) )2487 zsc_uw_1(:,:) = sustar(A2D( 0))**22488 zsc_vw_1(:,:) = ff_t(A2D( 0)) * sustke(A2D(0)) * phml(A2D(0))2489 WHERE ( l_conv(A2D((nn_hls-1))) ) 2490 zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 2491 zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) 2489 2492 ELSEWHERE 2490 zsc_uw_1(:,:) = sustar(A2D( 0))**22493 zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 2491 2494 zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) * & 2492 2495 & zsc_uw_1(:,:) 2493 zsc_vw_1(:,:) = ff_t(A2D( 0)) * sustke(A2D(0)) * phbl(A2D(0))2496 zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1))) 2494 2497 zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) * & 2495 2498 & zsc_vw_1(:,:) 2496 2499 ENDWHERE 2497 DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )2500 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2498 2501 IF ( l_conv(ji,jj) ) THEN 2499 2502 IF ( jk <= nmld(ji,jj) ) THEN … … 2530 2533 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask * ghamv ) 2531 2534 IF ( iom_use("zsc_uw_1_f") ) THEN 2532 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d )2535 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d ) 2533 2536 END IF 2534 2537 IF ( iom_use("zsc_vw_1_f") ) THEN 2535 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d )2538 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d ) 2536 2539 END IF 2537 2540 IF ( iom_use("zsc_uw_2_f") ) THEN 2538 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d )2541 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d ) 2539 2542 END IF 2540 2543 IF ( iom_use("zsc_vw_2_f") ) THEN 2541 osmdia2d(A2D( 0)) = tmask(A2D(0),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d )2544 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d ) 2542 2545 END IF 2543 2546 END IF … … 2548 2551 ! Make surface forced velocity non-gradient terms go to zero at the base 2549 2552 ! of the boundary layer. 2550 DO_3D( 0, 0, 0, 0, 2, jkm_bld )2553 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2551 2554 IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2552 2555 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj) ! ALMG to think about … … 2569 2572 z3ddz_pyc_2(:,:,:) = 0.0_wp 2570 2573 END IF 2571 DO_3D( 0, 0, 0, 0, 2, jkm_bld )2574 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2572 2575 IF ( l_conv (ji,jj) ) THEN 2573 2576 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. … … 2637 2640 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 2638 2641 END IF 2639 DO_3D( 0, 0, 0, 0, 2, jkm_bld )2642 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2640 2643 IF ( .NOT. l_conv (ji,jj) ) THEN 2641 2644 IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 2672 2675 END IF 2673 2676 ! 2674 DO_2D( 0, 0, 0, 0)2677 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2675 2678 ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp 2676 2679 ghams(ji,jj,nbld(ji,jj)) = 0.0_wp … … 2683 2686 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask * ghamv ) 2684 2687 IF ( iom_use("zviscos") ) THEN 2685 osmdia3d(A2D( 0),:) = wmask(A2D(0),:) * pviscos; CALL iom_put( "zviscos", osmdia3d )2688 osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pviscos; CALL iom_put( "zviscos", osmdia3d ) 2686 2689 END IF 2687 2690 END IF … … 2711 2714 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdx_mle ! MLE horiz gradients at u points 2712 2715 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdy_mle ! MLE horiz gradients at v points 2713 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient2716 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 2714 2717 ! 2715 2718 ! Local variables … … 2718 2721 REAL(wp) :: zc 2719 2722 REAL(wp) :: zN2_c ! Local buoyancy difference from 10m value 2720 REAL(wp), DIMENSION(A2D( 1)) :: ztm2721 REAL(wp), DIMENSION(A2D( 1)) :: zsm2722 REAL(wp), DIMENSION(A2D( 1),jpts) :: ztsm_midu2723 REAL(wp), DIMENSION(A2D( 1),jpts) :: ztsm_midv2724 REAL(wp), DIMENSION(A2D( 1),jpts) :: zabu2725 REAL(wp), DIMENSION(A2D( 1),jpts) :: zabv2726 REAL(wp), DIMENSION(A2D( 1)) :: zmld_midu2727 REAL(wp), DIMENSION(A2D( 1)) :: zmld_midv2723 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztm 2724 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsm 2725 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: ztsm_midu 2726 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: ztsm_midv 2727 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zabu 2728 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zabv 2729 REAL(wp), DIMENSION(A2D(nn_hls)) :: zmld_midu 2730 REAL(wp), DIMENSION(A2D(nn_hls)) :: zmld_midv 2728 2731 ! 2729 2732 ! == MLD used for MLE ==! … … 2731 2734 pmld(:,:) = 0.0_wp ! Here hmlp used as a dummy variable, integrating vertically N^2 2732 2735 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! Convert density criteria into N^2 criteria 2733 DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )2736 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) 2734 2737 ikt = mbkt(ji,jj) 2735 2738 pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm) 2736 2739 IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2737 2740 END_3D 2738 DO_2D( 1, 1, 1, 1)2741 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2739 2742 mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) ) ! Ensure mld_prof .ge. nbld 2740 2743 pmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) … … 2744 2747 ztm(:,:) = 0.0_wp 2745 2748 zsm(:,:) = 0.0_wp 2746 DO_3D( 1, 1, 1, 1, 1, ikmax )2749 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax ) 2747 2750 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 2748 2751 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) … … 2755 2758 zmld_midu(:,:) = 0.0_wp 2756 2759 ztsm_midu(:,:,:) = 10.0_wp 2757 DO_2D( 1, 0, 0, 0)2760 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 2758 2761 pdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm(ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 2759 2762 pdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm(ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) … … 2764 2767 zmld_midv(:,:) = 0.0_wp 2765 2768 ztsm_midv(:,:,:) = 10.0_wp 2766 DO_2D( 0, 0, 1, 0)2769 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 2767 2770 pdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 2768 2771 pdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) … … 2773 2776 CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm ) 2774 2777 CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm ) 2775 DO_2D( 1, 0, 0, 0)2778 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 2776 2779 pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) ) 2777 2780 END_2D 2778 DO_2D( 0, 0, 1, 0)2781 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 2779 2782 pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) ) 2780 2783 END_2D 2781 DO_2D( 0, 0, 0, 0)2784 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2782 2785 pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( pdbdx_mle(ji, jj) * pdbdx_mle(ji, jj) + pdbdy_mle(ji,jj ) * pdbdy_mle(ji,jj ) + & 2783 2786 & pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) ) … … 2807 2810 ! Outputs 2808 2811 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2809 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pwb_fk2810 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth2811 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phmle ! MLE depth2812 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux2813 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient2812 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pwb_fk 2813 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 2814 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phmle ! MLE depth 2815 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb_ent ! Buoyancy entrainment flux 2816 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 2814 2817 ! 2815 2818 ! Local variables 2816 2819 INTEGER :: ji, jj, jk ! Dummy loop indices 2817 REAL(wp), DIMENSION(A2D( 0)) :: znd_param2820 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: znd_param 2818 2821 REAL(wp) :: zthermal, zbeta 2819 2822 REAL(wp) :: zbuoy … … 2823 2826 REAL(wp) :: zdbdz_mle_int 2824 2827 ! 2825 znd_param(A2D( 0)) = 0.0_wp2826 ! 2827 DO_2D( 0, 0, 0, 0)2828 znd_param(A2D((nn_hls-1))) = 0.0_wp 2829 ! 2830 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2828 2831 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2829 2832 pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj) 2830 2833 END_2D 2831 2834 ! 2832 DO_2D( 0, 0, 0, 0)2835 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2833 2836 ! 2834 2837 IF ( l_conv(ji,jj) ) THEN … … 2858 2861 ! 2859 2862 ! Diagnosis 2860 DO_2D( 0, 0, 0, 0)2863 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2861 2864 ! 2862 2865 IF ( l_conv(ji,jj) ) THEN … … 2920 2923 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2921 2924 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pmld ! == Estimated FK BLD used for MLE horiz gradients == ! 2922 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: phmle ! MLE depth2923 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pvel_mle ! Velocity scale for dhdt with stable ML and FK2924 REAL(wp), DIMENSION(A2D( 0)), INTENT(inout) :: pdiff_mle ! Extra MLE vertical diff2925 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient2926 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: phbl ! BL depth2927 REAL(wp), DIMENSION(A2D( 0)), INTENT(in ) :: pwb0tot ! Total surface buoyancy flux including insolation2925 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: phmle ! MLE depth 2926 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pvel_mle ! Velocity scale for dhdt with stable ML and FK 2927 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdiff_mle ! Extra MLE vertical diff 2928 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 2929 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: phbl ! BL depth 2930 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pwb0tot ! Total surface buoyancy flux including insolation 2928 2931 ! 2929 2932 ! Local variables … … 2939 2942 ! 2940 2943 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE 2941 DO_2D( 0, 0, 0, 0)2944 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2942 2945 IF ( l_conv(ji,jj) ) THEN 2943 2946 ztmp = r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf … … 2948 2951 END_2D 2949 2952 ! Timestep mixed layer eddy depth 2950 DO_2D( 0, 0, 0, 0)2953 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2951 2954 IF ( l_mle(ji,jj) ) THEN ! MLE layer growing 2952 2955 ! Buoyancy gradient at base of MLE layer … … 2971 2974 ! 2972 2975 mld_prof(:,:) = 4 2973 DO_3D( 0, 0, 0, 0, 5, jpkm1 )2976 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 2974 2977 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk ) 2975 2978 END_3D 2976 DO_2D( 0, 0, 0, 0)2979 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2977 2980 phmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 2978 2981 END_2D … … 3126 3129 ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 3127 3130 z1_t2 = 2e-5_wp 3128 DO_2D( 0, 0, 0, 0)3131 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 3129 3132 r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 ) 3130 3133 END_2D … … 3167 3170 etmean(:,:,:) = 0.0_wp 3168 3171 ! 3169 DO_3D( 0, 0, 0, 0, 1, jpkm1 )3172 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 3170 3173 etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.0_wp, umask(ji-1,jj, jk) + umask(ji,jj,jk) + & 3171 3174 & vmask(ji, jj-1,jk) + vmask(ji,jj,jk) ) … … 3179 3182 etmean(:,:,:) = 0.0_wp 3180 3183 ! 3181 DO_3D( 0, 0, 0, 0, 1, jpkm1 )3184 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 3182 3185 etmean(ji,jj,jk) = tmask(ji, jj,jk) / MAX( 1.0_wp, 2.0_wp * tmask(ji,jj,jk) + & 3183 3186 & 0.5_wp * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) + & … … 3296 3299 ! 3297 3300 hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2 3298 DO_3D( 1, 1, 1, 1, 1, jpkm1 )3301 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 3299 3302 ikt = mbkt(ji,jj) 3300 3303 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0.0_wp ) * e3w(ji,jj,jk,Kmm) … … 3302 3305 END_3D 3303 3306 ! 3304 DO_2D( 1, 1, 1, 1)3307 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 3305 3308 iiki = MAX( 4, imld_rst(ji,jj) ) 3306 3309 hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth
Note: See TracChangeset
for help on using the changeset viewer.
