Changeset 14900
- Timestamp:
- 2021-05-25T21:26:19+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF
- Files:
-
- 2 edited
-
zdfosm.F90 (modified) (89 diffs)
-
zdfphy.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14889 r14900 61 61 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 62 62 !! dyn_osm : compute and add to u & v trensd the non-local flux 63 !! zdf_osm_iomput : iom_put wrapper that accepts arrays without halo 64 !! zdf_osm_iomput_2d : iom_put wrapper for 2D fields 65 !! zdf_osm_iomput_3d : iom_put wrapper for 3D fields 63 66 !!---------------------------------------------------------------------- 64 67 USE oce ! Ocean dynamics and active tracers … … 114 117 MODULE PROCEDURE zdf_osm_velocity_rotation_2d 115 118 MODULE PROCEDURE zdf_osm_velocity_rotation_3d 119 END INTERFACE 120 ! 121 INTERFACE zdf_osm_iomput 122 !!--------------------------------------------------------------------- 123 !! *** INTERFACE zdf_osm_iomput *** 124 !!--------------------------------------------------------------------- 125 MODULE PROCEDURE zdf_osm_iomput_2d 126 MODULE PROCEDURE zdf_osm_iomput_3d 116 127 END INTERFACE 117 128 ! … … 370 381 ! 371 382 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zhmle ! MLE depth - grid 372 REAL(wp), DIMENSION( jpi,jpj) :: zmld! ML depth on grid383 REAL(wp), DIMENSION(A2D(nn_hls)) :: zmld ! ML depth on grid 373 384 ! 374 385 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdh ! Pycnocline depth - grid 375 386 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdhdt ! BL depth tendency 376 387 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ! External temperature/salinity and buoyancy gradients 377 REAL(wp), DIMENSION( jpi,jpj):: zdtdx, zdtdy, zdsdx, zdsdy ! Horizontal gradients for Fox-Kemper parametrization388 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdtdx, zdtdy, zdsdx, zdsdy ! Horizontal gradients for Fox-Kemper parametrization 378 389 ! 379 390 REAL(wp), DIMENSION(A2D((nn_hls-1))) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient … … 384 395 ! 385 396 ! For calculating Ri#-dependent mixing 386 REAL(wp), DIMENSION( jpi,jpj) :: z2du ! u-shear^2387 REAL(wp), DIMENSION( jpi,jpj) :: z2dv ! v-shear^2388 REAL(wp) :: zrimix ! Spatial form of ri#-induced diffusion397 REAL(wp), DIMENSION(A2D(nn_hls)) :: z2du ! u-shear^2 398 REAL(wp), DIMENSION(A2D(nn_hls)) :: z2dv ! v-shear^2 399 REAL(wp) :: zrimix ! Spatial form of ri#-induced diffusion 389 400 ! 390 401 ! Temporary variables … … 401 412 REAL(wp), PARAMETER :: pp_large = -1e10_wp 402 413 ! 403 nbld(:,:) = 0 404 nmld(:,:) = 0 405 sustke(:,:) = pp_large 406 l_pyc(:,:) = .FALSE. 407 l_flux(:,:) = .FALSE. 408 l_mle(:,:) = .FALSE. 414 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 415 nmld(ji,jj) = 0 416 sustke(ji,jj) = pp_large 417 l_pyc(ji,jj) = .FALSE. 418 l_flux(ji,jj) = .FALSE. 419 l_mle(ji,jj) = .FALSE. 420 END_2D 409 421 ! Mixed layer 410 422 ! No initialization of zhbl or zhml (or zdh?) … … 413 425 IF ( ln_osm_mle ) THEN ! Only initialise arrays if needed 414 426 zdtdx(:,:) = pp_large ; zdtdy(:,:) = pp_large ; zdsdx(:,:) = pp_large 415 zdsdy(:,:) = pp_large ; dbdx_mle(:,:) = pp_large ; dbdy_mle(:,:) = pp_large427 zdsdy(:,:) = pp_large 416 428 zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large 417 429 zhmle(:,:) = pp_large ; zmld(:,:) = pp_large 430 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 431 dbdx_mle(ji,jj) = pp_large 432 dbdy_mle(ji,jj) = pp_large 433 END_2D 418 434 ENDIF 419 435 zhbl_t(:,:) = pp_large … … 422 438 zviscos(:,:,:) = 0.0_wp 423 439 ! 424 ghamt(:,:,:) = pp_large ; ghams(:,:,:) = pp_large 425 ghamt(A2D((nn_hls-1)),:) = 0.0_wp ; ghams(A2D((nn_hls-1)),:) = 0.0_wp 426 ghamu(:,:,:) = pp_large ; ghamv(:,:,:) = pp_large 427 ghamu(A2D((nn_hls-1)),:) = 0.0_wp ; ghamv(A2D((nn_hls-1)),:) = 0.0_wp 440 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 441 ghamt(ji,jj,jk) = pp_large 442 ghams(ji,jj,jk) = pp_large 443 ghamu(ji,jj,jk) = pp_large 444 ghamv(ji,jj,jk) = pp_large 445 END_3D 446 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 447 ghamt(ji,jj,jk) = 0.0_wp 448 ghams(ji,jj,jk) = 0.0_wp 449 ghamu(ji,jj,jk) = 0.0_wp 450 ghamv(ji,jj,jk) = 0.0_wp 451 END_3D 428 452 ! 429 453 zdiff_mle(:,:) = 0.0_wp … … 431 455 ! Ensure only positive hbl values are accessed when using extended halo 432 456 ! (nn_hls==2) 433 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )457 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 434 458 hbl(ji,jj) = MAX( hbl(ji,jj), epsln ) 435 459 END_2D … … 442 466 zz0 = rn_abs ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands 443 467 zz1 = 1.0_wp - rn_abs 444 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )468 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 445 469 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp ! Surface downward irradiance (so always +ve) 446 470 zradh(ji,jj) = zrad0(ji,jj) * & ! Downwards irradiance at base of boundary layer … … 453 477 & zradav ) ! over depth of OSBL 454 478 END_2D 455 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )479 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 456 480 sws0(ji,jj) = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + & ! Upwards surface salinity flux 457 481 & sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) ! for non-local term … … 464 488 & grav * zbeta * swsav(ji,jj) ! OBSBL 465 489 END_2D 466 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )490 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 467 491 suw0(ji,jj) = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) ! Surface upward velocity fluxes 468 492 zvw0 = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) … … 476 500 ! Assume constant La#=0.3 477 501 CASE(0) 478 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )502 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 479 503 zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 480 504 zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 … … 485 509 ! Assume Pierson-Moskovitz wind-wave spectrum 486 510 CASE(1) 487 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )511 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 488 512 ! Use wind speed wndm included in sbc_oce module 489 513 sustke(ji,jj) = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp ) … … 494 518 zfac = 2.0_wp * rpi / 16.0_wp 495 519 ! 496 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )520 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 497 521 IF ( hsw(ji,jj) > 1e-4_wp ) THEN 498 522 ! Use wave fields … … 511 535 IF (ln_zdfosm_ice_shelter) THEN 512 536 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 513 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )537 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 514 538 sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) 515 539 dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) … … 524 548 ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested. 525 549 IF(nn_osm_wave > 0) THEN 526 sustke(A2D((nn_hls-1))) = rn_zdfosm_adjust_sd * sustke(A2D((nn_hls-1))) 550 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 551 sustke(ji,jj) = rn_zdfosm_adjust_sd * sustke(ji,jj) 552 END_2D 527 553 END IF 528 554 CASE(1) … … 531 557 zsqrtpi = SQRT(rpi) 532 558 z_two_thirds = 2.0_wp / 3.0_wp 533 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )559 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 534 560 zthickness = rn_osm_hblfrac*hbl(ji,jj) 535 561 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) … … 545 571 ! Assumes approximate depth profile of SD from Breivik (2016) 546 572 zsqrtpi = SQRT(rpi) 547 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )573 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 548 574 zthickness = rn_osm_hblfrac*hbl(ji,jj) 549 575 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) … … 567 593 ! Langmuir velocity scale (swstrl), La # (sla) 568 594 ! Mixed scale (svstr), convective velocity scale (swstrc) 569 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )595 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 570 596 ! Langmuir velocity scale (swstrl), at T-point 571 597 swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird … … 594 620 ! BL must be always 4 levels deep. 595 621 ! For calculation of lateral buoyancy gradients for FK in 596 ! zdf_osm_zmld_horizontal_gradients need halo values for nbld, so must 597 ! previously exist for hbl also. 622 ! zdf_osm_zmld_horizontal_gradients need halo values for nbld 598 623 ! 599 624 ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 600 625 ! ########################################################################## 601 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 602 nbld(:,:) = 4 603 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 ) 604 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 626 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 627 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,4,Kmm) ) 628 END_2D 629 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 630 nbld(ji,jj) = 4 631 END_2D 632 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 ) 633 IF ( MAX( hbl(ji,jj), gdepw(ji,jj,4,Kmm) ) >= gdepw(ji,jj,jk,Kmm) ) THEN 605 634 nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 606 635 ENDIF … … 608 637 ! ########################################################################## 609 638 ! 610 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )639 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 611 640 zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) 612 641 nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) ) … … 617 646 ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere 618 647 jk_ext(:,:) = 1 ! ag 19/03 619 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 620 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 621 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 622 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03 623 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 624 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & 625 & av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 648 CALL zdf_osm_vertical_average( Kbb, Kmm, & 649 & nbld(A2D((nn_hls-1))), av_t_bl(A2D((nn_hls-1))), & 650 & av_s_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), & 651 & av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), & 652 & jk_ext, av_dt_bl(A2D((nn_hls-1))), & 653 & av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), & 654 & av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) ) 655 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1 ! ag 19/03 656 CALL zdf_osm_vertical_average( Kbb, Kmm, & 657 & nmld - 1, av_t_ml(A2D((nn_hls-1))), & 658 & av_s_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), & 659 & av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), & 660 & jk_ext, av_dt_ml(A2D((nn_hls-1))), & 661 & av_ds_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), & 662 & av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) ) 626 663 ! Velocity components in frame aligned with surface stress 627 CALL zdf_osm_velocity_rotation( av_u_ml , av_v_ml)628 CALL zdf_osm_velocity_rotation( av_du_ml , av_dv_ml)629 CALL zdf_osm_velocity_rotation( av_u_bl , av_v_bl)630 CALL zdf_osm_velocity_rotation( av_du_bl , av_dv_bl)664 CALL zdf_osm_velocity_rotation( av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))) ) 665 CALL zdf_osm_velocity_rotation( av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) ) 666 CALL zdf_osm_velocity_rotation( av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))) ) 667 CALL zdf_osm_velocity_rotation( av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) ) 631 668 ! 632 669 ! Determine the state of the OSBL, stable/unstable, shear/no shear … … 636 673 IF ( ln_osm_mle ) THEN 637 674 ! Fox-Kemper Scheme 638 mld_prof = 4 639 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 675 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 676 mld_prof(ji,jj) = 4 677 END_2D 678 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 640 679 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 641 680 END_3D 642 CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D((nn_hls-1))), av_t_mle, av_s_mle, & 643 & av_b_mle, av_u_mle, av_v_mle ) 681 CALL zdf_osm_vertical_average( Kbb, Kmm, & 682 & mld_prof(A2D((nn_hls-1))), av_t_mle(A2D((nn_hls-1))), & 683 & av_s_mle(A2D((nn_hls-1))), av_b_mle(A2D((nn_hls-1))), & 684 & av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))) ) 644 685 ! 645 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )686 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 646 687 zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 647 688 END_2D … … 649 690 ! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients 650 691 CALL zdf_osm_zmld_horizontal_gradients( Kmm, zmld, zdtdx, zdtdy, zdsdx, & 651 & zdsdy, dbdx_mle, dbdy_mle,zdbds_mle )692 & zdsdy, zdbds_mle ) 652 693 ! Calculate max vertical FK flux zwb_fk & set logical descriptors 653 694 CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zhbl, zhmle, zwb_ent, & 654 695 & zdbds_mle ) 655 696 ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle 656 CALL zdf_osm_mle_parameters( Kmm, zmld , zhmle, zvel_mle, zdiff_mle, &697 CALL zdf_osm_mle_parameters( Kmm, zmld(A2D((nn_hls-1))), zhmle, zvel_mle, zdiff_mle, & 657 698 & zdbds_mle, zhbl, zwb0tot ) 658 699 ELSE ! ln_osm_mle 659 700 ! FK not selected, Boundary Layer only. 660 l_pyc(:,:) = .TRUE.661 l_flux(:,:) = .FALSE.662 l_mle(:,:)= .FALSE.663 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )701 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 702 l_pyc(ji,jj) = .TRUE. 703 l_flux(ji,jj) = .FALSE. 704 l_mle(ji,jj) = .FALSE. 664 705 IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 665 706 END_2D … … 688 729 ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here.. 689 730 jk_ext(:,:) = 1 ! ag 19/03 690 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 691 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 692 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 693 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1 ! ag 19/03 694 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 695 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & 696 & av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) ! ag 19/03 731 CALL zdf_osm_vertical_average( Kbb, Kmm, & 732 & nbld(A2D((nn_hls-1))), av_t_bl(A2D((nn_hls-1))), & 733 & av_s_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), & 734 & av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), & 735 & jk_ext, av_dt_bl(A2D((nn_hls-1))), & 736 & av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), & 737 & av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) ) 738 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1 ! ag 19/03 739 CALL zdf_osm_vertical_average( Kbb, Kmm, & 740 & nmld(A2D((nn_hls-1))) - 1, av_t_ml(A2D((nn_hls-1))), & 741 & av_s_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), & 742 & av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), & 743 & jk_ext, av_dt_ml(A2D((nn_hls-1))), & 744 & av_ds_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), & 745 & av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) ) ! ag 19/03 697 746 ! 698 747 ! Rate of change of hbl … … 700 749 & zdbdz_bl_ext, zwb_fk_b, zwb_fk, zvel_mle ) 701 750 ! Test if surface boundary layer coupled to bottom 702 l_coup(:,:) = .FALSE. ! ag 19/03703 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )751 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 752 l_coup(ji,jj) = .FALSE. ! ag 19/03 704 753 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 705 754 ! Adjustment to represent limiting by ocean bottom … … 713 762 END_2D 714 763 ! 715 nmld(:,:) = nbld(A2D((nn_hls-1))) ! use nmld to hold previous blayer index 716 nbld(:,:) = 4 717 ! 718 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 ) 764 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 765 nmld(ji,jj) = nbld(ji,jj) ! use nmld to hold previous blayer index 766 nbld(ji,jj) = 4 767 END_2D 768 ! 769 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 ) 719 770 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 720 771 nbld(ji,jj) = jk … … 731 782 ! Recalculate BL averages and differences using new BL depth 732 783 jk_ext(:,:) = 1 ! ag 19/03 733 CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl, & 734 & av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl, & 735 & av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 784 CALL zdf_osm_vertical_average( Kbb, Kmm, & 785 & nbld(A2D((nn_hls-1))), av_t_bl(A2D((nn_hls-1))), & 786 & av_s_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), & 787 & av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), & 788 & jk_ext, av_dt_bl(A2D((nn_hls-1))), & 789 & av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), & 790 & av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) ) 736 791 ! 737 792 CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zhbl, & … … 739 794 ! 740 795 ! Reset l_pyc before calculating terms in the flux-gradient relationship 741 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )796 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 742 797 IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR. & 743 798 & nbld(ji,jj) - nmld(ji,jj) == 1 .OR. zdhdt(ji,jj) < 0.0_wp ) THEN ! ag 19/03 … … 753 808 END_2D 754 809 ! 755 dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D((nn_hls-1)))/ 3.0_wp ) ! Limit delta for shallow boundary layers for calculating 756 ! ! flux-gradient terms 810 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Limit delta for shallow boundary layers for calculating 811 dstokes(ji,jj) = MIN ( dstokes(ji,jj), hbl(ji,jj) / 3.0_wp ) ! flux-gradient terms 812 END_2D 813 ! 757 814 ! 758 815 ! Average over the depth of the mixed layer in the convective boundary layer … … 760 817 ! Recalculate ML averages and differences using new ML depth 761 818 jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1 ! ag 19/03 762 CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml, & 763 & av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml, & 764 & av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 819 CALL zdf_osm_vertical_average( Kbb, Kmm, & 820 & nmld(A2D((nn_hls-1))) - 1, av_t_ml(A2D((nn_hls-1))), & 821 & av_s_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), & 822 & av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), & 823 & jk_ext, av_dt_ml(A2D((nn_hls-1))), & 824 & av_ds_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), & 825 & av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) ) 765 826 ! 766 827 CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 767 828 ! Rotate mean currents and changes onto wind aligned co-ordinates 768 CALL zdf_osm_velocity_rotation( av_u_ml , av_v_ml)769 CALL zdf_osm_velocity_rotation( av_du_ml , av_dv_ml)770 CALL zdf_osm_velocity_rotation( av_u_bl , av_v_bl)771 CALL zdf_osm_velocity_rotation( av_du_bl , av_dv_bl)829 CALL zdf_osm_velocity_rotation( av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))) ) 830 CALL zdf_osm_velocity_rotation( av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) ) 831 CALL zdf_osm_velocity_rotation( av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))) ) 832 CALL zdf_osm_velocity_rotation( av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) ) 772 833 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 773 834 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship … … 791 852 ! 792 853 ! Rotate non-gradient velocity terms back to model reference frame 793 CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:), .FALSE., 2, nbld(A2D((nn_hls-1))) ) 854 CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:), & 855 & .FALSE., 2, & 856 & nbld(A2D((nn_hls-1))) ) 794 857 ! 795 858 ! KPP-style Ri# mixing 796 859 IF ( ln_kpprimix ) THEN 797 860 jkflt = jpk 798 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )861 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 799 862 IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj) 800 863 END_2D … … 807 870 & wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 808 871 END_2D 809 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )872 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 810 873 IF ( jk > nbld(ji,jj) ) THEN 811 874 ! Shear prod. at w-point weightened by mask … … 826 889 ! KPP-style set diffusivity large if unstable below BL 827 890 IF ( ln_convmix) THEN 828 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )891 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 829 892 DO jk = nbld(ji,jj) + 1, jpkm1 830 893 IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) … … 834 897 ! 835 898 IF ( ln_osm_mle ) THEN ! Set up diffusivity and non-gradient mixing 836 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )899 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 837 900 IF ( l_flux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 838 901 ! Calculate MLE flux contribution from surface fluxes … … 867 930 ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 868 931 ! GN 25/8: need to change tmask --> wmask 869 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )932 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 870 933 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 871 934 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 872 935 END_3D 873 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and874 ! v grids875 IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp, &876 & ghamv, 'W', 1.0_wp )877 DO_3D( 0, 0, 0, 0, 2, jpkm1 )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) ) * &879 & umask(ji,jj,jk)880 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * &881 & vmask(ji,jj,jk)882 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)883 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)884 END_3D885 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)886 CALL lbc_lnk( 'zdfosm', hbl, 'T', 1.0_wp, &887 & hmle, 'T', 1.0_wp )888 936 ! 889 937 IF ( ln_dia_osm ) THEN … … 891 939 ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1) 892 940 CASE(0:1) 893 IF ( iom_use("us_x") ) THEN ! x surface Stokes drift 894 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * scos_wind 895 CALL iom_put( "us_x", osmdia2d ) 896 END IF 897 IF ( iom_use("us_y") ) THEN ! y surface Stokes drift 898 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * ssin_wind 899 CALL iom_put( "us_y", osmdia2d ) 900 END IF 901 IF ( iom_use("wind_wave_abs_power") ) THEN 902 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 903 CALL iom_put( "wind_wave_abs_power", osmdia2d ) 904 END IF 941 CALL zdf_osm_iomput( "us_x", tmask(A2D(0),1) * sustke(A2D(0)) * scos_wind(A2D(0)) ) ! x surface Stokes drift 942 CALL zdf_osm_iomput( "us_y", tmask(A2D(0),1) * sustke(A2D(0)) * scos_wind(A2D(0)) ) ! y surface Stokes drift 943 CALL zdf_osm_iomput( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar(A2D(0))**2 * sustke(A2D(0)) ) 905 944 ! Stokes drift read in from sbcwave (=2). 906 945 CASE(2:3) 907 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 908 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 909 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! Wave mean period 910 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! Significant wave height 911 IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp / & ! Wave mean period from NP 912 & ( 0.877_wp * grav ) ) * & ! spectrum 913 & wndm * tmask(:,:,1) ) 914 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 * & ! Significant wave 915 & tmask(:,:,1) ) ! height from NP spectrum 916 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 917 IF ( iom_use("wind_wave_abs_power") ) THEN 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 ) 919 CALL iom_put( "wind_wave_abs_power", osmdia2d ) 920 END IF 946 CALL zdf_osm_iomput( "us_x", ut0sd(A2D(0)) * umask(A2D(0),1) ) ! x surface Stokes drift 947 CALL zdf_osm_iomput( "us_y", vt0sd(A2D(0)) * vmask(A2D(0),1) ) ! y surface Stokes drift 948 CALL zdf_osm_iomput( "wmp", wmp(A2D(0)) * tmask(A2D(0),1) ) ! Wave mean period 949 CALL zdf_osm_iomput( "hsw", hsw(A2D(0)) * tmask(A2D(0),1) ) ! Significant wave height 950 CALL zdf_osm_iomput( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp / ( 0.877_wp * grav ) ) * & ! Wave mean period from NP 951 & wndm(A2D(0)) * tmask(A2D(0),1) ) ! spectrum 952 CALL zdf_osm_iomput( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 * tmask(A2D(0),1) ) ! Significant wave height from 953 ! ! NP spectrum 954 CALL zdf_osm_iomput( "wndm", wndm(A2D(0)) * tmask(A2D(0),1) ) ! U_10 955 CALL zdf_osm_iomput( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar(A2D(0))**2 * & 956 & SQRT( ut0sd(A2D(0))**2 + vt0sd(A2D(0))**2 ) ) 921 957 END SELECT 922 IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL> 923 IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL> 924 IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL> 925 IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL> 926 IF ( iom_use("zwth0") ) THEN ! <Tw_0> 927 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0; CALL iom_put( "zwth0", osmdia2d ) 928 END IF 929 IF ( iom_use("zws0") ) THEN ! <Sw_0> 930 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sws0; CALL iom_put( "zws0", osmdia2d ) 931 END IF 932 IF ( iom_use("zwb0") ) THEN ! <Sw_0> 933 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swb0; CALL iom_put( "zwb0", osmdia2d ) 934 END IF 935 IF ( iom_use("zwbav") ) THEN ! Upward BL-avged turb buoyancy flux 936 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0; CALL iom_put( "zwbav", osmdia2d ) 937 END IF 938 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! Boundary-layer depth 939 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld ) ! Boundary-layer max k 940 IF ( iom_use("zdt_bl") ) THEN ! dt at ml base 941 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_bl; CALL iom_put( "zdt_bl", osmdia2d ) 942 END IF 943 IF ( iom_use("zds_bl") ) THEN ! ds at ml base 944 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_bl; CALL iom_put( "zds_bl", osmdia2d ) 945 END IF 946 IF ( iom_use("zdb_bl") ) THEN ! db at ml base 947 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_bl; CALL iom_put( "zdb_bl", osmdia2d ) 948 END IF 949 IF ( iom_use("zdu_bl") ) THEN ! du at ml base 950 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_du_bl; CALL iom_put( "zdu_bl", osmdia2d ) 951 END IF 952 IF ( iom_use("zdv_bl") ) THEN ! dv at ml base 953 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dv_bl; CALL iom_put( "zdv_bl", osmdia2d ) 954 END IF 955 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 956 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 957 IF ( iom_use("zdt_ml") ) THEN ! dt at ml base 958 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_ml; CALL iom_put( "zdt_ml", osmdia2d ) 959 END IF 960 IF ( iom_use("zds_ml") ) THEN ! ds at ml base 961 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_ml; CALL iom_put( "zds_ml", osmdia2d ) 962 END IF 963 IF ( iom_use("zdb_ml") ) THEN ! db at ml base 964 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_ml; CALL iom_put( "zdb_ml", osmdia2d ) 965 END IF 966 IF ( iom_use("dstokes") ) THEN ! Stokes drift penetration depth 967 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * dstokes; CALL iom_put( "dstokes", osmdia2d ) 968 END IF 969 IF ( iom_use("zustke") ) THEN ! Stokes drift magnitude at T-points 970 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke; CALL iom_put( "zustke", osmdia2d ) 971 END IF 972 IF ( iom_use("zwstrc") ) THEN ! Convective velocity scale 973 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrc; CALL iom_put( "zwstrc", osmdia2d ) 974 END IF 975 IF ( iom_use("zwstrl") ) THEN ! Langmuir velocity scale 976 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrl; CALL iom_put( "zwstrl", osmdia2d ) 977 END IF 978 IF ( iom_use("zustar") ) THEN ! Friction velocity scale 979 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustar; CALL iom_put( "zustar", osmdia2d ) 980 END IF 981 IF ( iom_use("zvstr") ) THEN ! Mixed velocity scale 982 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * svstr; CALL iom_put( "zvstr", osmdia2d ) 983 END IF 984 IF ( iom_use("zla") ) THEN ! Langmuir # 985 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sla; CALL iom_put( "zla", osmdia2d ) 986 END IF 987 IF ( iom_use("wind_power") ) THEN ! BL depth internal to zdf_osm routine 988 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**3 989 CALL iom_put( "wind_power", osmdia2d ) 990 END IF 991 IF ( iom_use("wind_wave_power") ) THEN 992 osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 993 CALL iom_put( "wind_wave_power", osmdia2d ) 994 END IF 995 IF ( iom_use("zhbl") ) THEN ! BL depth internal to zdf_osm routine 996 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhbl; CALL iom_put( "zhbl", osmdia2d ) 997 END IF 998 IF ( iom_use("zhml") ) THEN ! ML depth internal to zdf_osm routine 999 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhml; CALL iom_put( "zhml", osmdia2d ) 1000 END IF 1001 IF ( iom_use("imld") ) THEN ! Index for ML depth internal to zdf_osm routine 1002 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * nmld; CALL iom_put( "imld", osmdia2d ) 1003 END IF 1004 IF ( iom_use("jp_ext") ) THEN ! =1 if pycnocline resolved internal to zdf_osm routine 1005 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * jk_ext; CALL iom_put( "jp_ext", osmdia2d ) 1006 END IF 1007 IF ( iom_use("j_ddh") ) THEN ! Index forpyc thicknessh internal to zdf_osm routine 1008 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * n_ddh; CALL iom_put( "j_ddh", osmdia2d ) 1009 END IF 1010 IF ( iom_use("zshear") ) THEN ! Shear production of TKE internal to zdf_osm routine 1011 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zshear; CALL iom_put( "zshear", osmdia2d ) 1012 END IF 1013 IF ( iom_use("zdh") ) THEN ! Pyc thicknessh internal to zdf_osm routine 1014 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdh; CALL iom_put( "zdh", osmdia2d ) 1015 END IF 1016 IF ( iom_use("zhol") ) THEN ! ML depth internal to zdf_osm routine 1017 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * shol; CALL iom_put( "zhol", osmdia2d ) 1018 END IF 1019 IF ( iom_use("zwb_ent") ) THEN ! Upward turb buoyancy entrainment flux 1020 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_ent; CALL iom_put( "zwb_ent", osmdia2d ) 1021 END IF 1022 IF ( iom_use("zt_ml") ) THEN ! Average T in ML 1023 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_t_ml; CALL iom_put( "zt_ml", osmdia2d ) 1024 END IF 1025 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1026 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1027 IF ( iom_use("zwb_fk") ) THEN ! FK b flux 1028 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk; CALL iom_put( "zwb_fk", osmdia2d ) 1029 END IF 1030 IF ( iom_use("zwb_fk_b") ) THEN ! FK b flux averaged over ML 1031 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk_b; CALL iom_put( "zwb_fk_b", osmdia2d ) 1032 END IF 1033 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof ) ! FK layer max k 1034 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1035 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1036 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1037 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1038 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1039 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1040 IF ( iom_use("zdiff_mle") ) THEN ! FK diff in MLE at t-pt 1041 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d ) 1042 END IF 1043 IF ( iom_use("zvel_mle") ) THEN ! FK diff in MLE at t-pt 1044 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zvel_mle", osmdia2d ) 1045 END IF 958 CALL zdf_osm_iomput( "zwth0", tmask(A2D(0),1) * swth0(A2D(0)) ) ! <Tw_0> 959 CALL zdf_osm_iomput( "zws0", tmask(A2D(0),1) * sws0(A2D(0)) ) ! <Sw_0> 960 CALL zdf_osm_iomput( "zwb0", tmask(A2D(0),1) * swb0(A2D(0)) ) ! <Sw_0> 961 CALL zdf_osm_iomput( "zwbav", tmask(A2D(0),1) * swth0(A2D(0)) ) ! Upward BL-avged turb buoyancy flux 962 CALL zdf_osm_iomput( "ibld", tmask(A2D(0),1) * nbld(A2D(0)) ) ! Boundary-layer max k 963 CALL zdf_osm_iomput( "zdt_bl", tmask(A2D(0),1) * av_dt_bl(A2D(0)) ) ! dt at ml base 964 CALL zdf_osm_iomput( "zds_bl", tmask(A2D(0),1) * av_ds_bl(A2D(0)) ) ! ds at ml base 965 CALL zdf_osm_iomput( "zdb_bl", tmask(A2D(0),1) * av_db_bl(A2D(0)) ) ! db at ml base 966 CALL zdf_osm_iomput( "zdu_bl", tmask(A2D(0),1) * av_du_bl(A2D(0)) ) ! du at ml base 967 CALL zdf_osm_iomput( "zdv_bl", tmask(A2D(0),1) * av_dv_bl(A2D(0)) ) ! dv at ml base 968 CALL zdf_osm_iomput( "dh", tmask(A2D(0),1) * dh(A2D(0)) ) ! Initial boundary-layer depth 969 CALL zdf_osm_iomput( "hml", tmask(A2D(0),1) * hml(A2D(0)) ) ! Initial boundary-layer depth 970 CALL zdf_osm_iomput( "zdt_ml", tmask(A2D(0),1) * av_dt_ml(A2D(0)) ) ! dt at ml base 971 CALL zdf_osm_iomput( "zds_ml", tmask(A2D(0),1) * av_ds_ml(A2D(0)) ) ! ds at ml base 972 CALL zdf_osm_iomput( "zdb_ml", tmask(A2D(0),1) * av_db_ml(A2D(0)) ) ! db at ml base 973 CALL zdf_osm_iomput( "dstokes", tmask(A2D(0),1) * dstokes(A2D(0)) ) ! Stokes drift penetration depth 974 CALL zdf_osm_iomput( "zustke", tmask(A2D(0),1) * sustke(A2D(0)) ) ! Stokes drift magnitude at T-points 975 CALL zdf_osm_iomput( "zwstrc", tmask(A2D(0),1) * swstrc(A2D(0)) ) ! Convective velocity scale 976 CALL zdf_osm_iomput( "zwstrl", tmask(A2D(0),1) * swstrl(A2D(0)) ) ! Langmuir velocity scale 977 CALL zdf_osm_iomput( "zustar", tmask(A2D(0),1) * sustar(A2D(0)) ) ! Friction velocity scale 978 CALL zdf_osm_iomput( "zvstr", tmask(A2D(0),1) * svstr(A2D(0)) ) ! Mixed velocity scale 979 CALL zdf_osm_iomput( "zla", tmask(A2D(0),1) * sla(A2D(0)) ) ! Langmuir # 980 CALL zdf_osm_iomput( "wind_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * & ! BL depth internal to zdf_osm routine 981 & sustar(A2D(0))**3 ) 982 CALL zdf_osm_iomput( "wind_wave_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * & 983 & sustar(A2D(0))**2 * sustke(A2D(0)) ) 984 CALL zdf_osm_iomput( "zhbl", tmask(A2D(0),1) * zhbl(A2D(0)) ) ! BL depth internal to zdf_osm routine 985 CALL zdf_osm_iomput( "zhml", tmask(A2D(0),1) * zhml(A2D(0)) ) ! ML depth internal to zdf_osm routine 986 CALL zdf_osm_iomput( "imld", tmask(A2D(0),1) * nmld(A2D(0)) ) ! Index for ML depth internal to zdf_osm routine 987 CALL zdf_osm_iomput( "jp_ext", tmask(A2D(0),1) * jk_ext(A2D(0)) ) ! =1 if pycnocline resolved internal to zdf_osm routine 988 CALL zdf_osm_iomput( "j_ddh", tmask(A2D(0),1) * n_ddh(A2D(0)) ) ! Index forpyc thicknessh internal to zdf_osm routine 989 CALL zdf_osm_iomput( "zshear", tmask(A2D(0),1) * zshear(A2D(0)) ) ! Shear production of TKE internal to zdf_osm routine 990 CALL zdf_osm_iomput( "zdh", tmask(A2D(0),1) * zdh(A2D(0)) ) ! Pyc thicknessh internal to zdf_osm routine 991 CALL zdf_osm_iomput( "zhol", tmask(A2D(0),1) * shol(A2D(0)) ) ! ML depth internal to zdf_osm routine 992 CALL zdf_osm_iomput( "zwb_ent", tmask(A2D(0),1) * zwb_ent(A2D(0)) ) ! Upward turb buoyancy entrainment flux 993 CALL zdf_osm_iomput( "zt_ml", tmask(A2D(0),1) * av_t_ml(A2D(0)) ) ! Average T in ML 994 CALL zdf_osm_iomput( "zmld", tmask(A2D(0),1) * zmld(A2D(0)) ) ! FK target layer depth 995 CALL zdf_osm_iomput( "zwb_fk", tmask(A2D(0),1) * zwb_fk(A2D(0)) ) ! FK b flux 996 CALL zdf_osm_iomput( "zwb_fk_b", tmask(A2D(0),1) * zwb_fk_b(A2D(0)) ) ! FK b flux averaged over ML 997 CALL zdf_osm_iomput( "mld_prof", tmask(A2D(0),1) * mld_prof(A2D(0)) ) ! FK layer max k 998 CALL zdf_osm_iomput( "zdtdx", umask(A2D(0),1) * zdtdx(A2D(0)) ) ! FK dtdx at u-pt 999 CALL zdf_osm_iomput( "zdtdy", vmask(A2D(0),1) * zdtdy(A2D(0)) ) ! FK dtdy at v-pt 1000 CALL zdf_osm_iomput( "zdsdx", umask(A2D(0),1) * zdsdx(A2D(0)) ) ! FK dtdx at u-pt 1001 CALL zdf_osm_iomput( "zdsdy", vmask(A2D(0),1) * zdsdy(A2D(0)) ) ! FK dsdy at v-pt 1002 CALL zdf_osm_iomput( "dbdx_mle", umask(A2D(0),1) * dbdx_mle(A2D(0)) ) ! FK dbdx at u-pt 1003 CALL zdf_osm_iomput( "dbdy_mle", vmask(A2D(0),1) * dbdy_mle(A2D(0)) ) ! FK dbdy at v-pt 1004 CALL zdf_osm_iomput( "zdiff_mle", tmask(A2D(0),1) * zdiff_mle(A2D(0)) ) ! FK diff in MLE at t-pt 1005 CALL zdf_osm_iomput( "zvel_mle", tmask(A2D(0),1) * zdiff_mle(A2D(0)) ) ! FK diff in MLE at t-pt 1006 END IF 1007 ! 1008 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and 1009 ! v grids 1010 IF ( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Finalise ghamu, ghamv, hbl, and hmle only after full domain has been processed 1011 IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp, & 1012 & ghamv, 'W', 1.0_wp ) 1013 DO jk = 2, jpkm1 1014 DO jj = Njs0, Nje0 1015 DO ji = Nis0, Nie0 1016 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) ) * & 1017 & umask(ji,jj,jk) 1018 ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * & 1019 & vmask(ji,jj,jk) 1020 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk) 1021 ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk) 1022 END DO 1023 END DO 1024 END DO 1025 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1026 CALL lbc_lnk( 'zdfosm', hbl, 'T', 1.0_wp, & 1027 & hmle, 'T', 1.0_wp ) 1028 ! 1029 CALL zdf_osm_iomput( "ghamt", tmask*ghamt ) ! <Tw_NL> 1030 CALL zdf_osm_iomput( "ghams", tmask*ghams ) ! <Sw_NL> 1031 CALL zdf_osm_iomput( "ghamu", umask*ghamu ) ! <uw_NL> 1032 CALL zdf_osm_iomput( "ghamv", vmask*ghamv ) ! <vw_NL> 1033 CALL zdf_osm_iomput( "hbl", tmask(:,:,1) * hbl ) ! Boundary-layer depth 1034 CALL zdf_osm_iomput( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1046 1035 END IF 1047 1036 ! … … 1079 1068 ! 1080 1069 ! Averages over depth of boundary layer 1081 pt(:,:) = 0.0_wp 1082 ps(:,:) = 0.0_wp 1083 pu(:,:) = 0.0_wp 1084 pv(:,:) = 0.0_wp 1070 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1071 pt(ji,jj) = 0.0_wp 1072 ps(ji,jj) = 0.0_wp 1073 pu(ji,jj) = 0.0_wp 1074 pv(ji,jj) = 0.0_wp 1075 END_2D 1085 1076 zthick(:,:) = epsln 1086 1077 jkflt = jpk 1087 1078 jkmax = 0 1088 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1079 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1089 1080 IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj) 1090 1081 IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 1091 1082 END_2D 1092 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt ) ! Upper, flat part of layer1083 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt ) ! Upper, flat part of layer 1093 1084 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 1094 1085 pt(ji,jj) = pt(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) … … 1101 1092 & MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1102 1093 END_3D 1103 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax ) ! Lower, non-flat part of layer1094 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax ) ! Lower, non-flat part of layer 1104 1095 IF ( knlev(ji,jj) >= jk ) THEN 1105 1096 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) … … 1114 1105 END IF 1115 1106 END_3D 1116 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1107 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1117 1108 pt(ji,jj) = pt(ji,jj) / zthick(ji,jj) 1118 1109 ps(ji,jj) = ps(ji,jj) / zthick(ji,jj) … … 1126 1117 ! Differences between vertical averages and values at an external layer 1127 1118 IF ( PRESENT( kp_ext ) ) THEN 1128 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1119 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1129 1120 ibld_ext = knlev(ji,jj) + kp_ext(ji,jj) 1130 1121 IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN ! ag 09/03 … … 1171 1162 zfwd = 1.0_wp 1172 1163 IF( PRESENT(fwd) .AND. ( .NOT. fwd ) ) zfwd = -1.0_wp 1173 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1164 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1174 1165 ztmp = pu(ji,jj) 1175 1166 pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj) … … 1208 1199 IF( PRESENT(knlev) ) THEN 1209 1200 jkmax = 0 1210 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1201 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1211 1202 IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 1212 1203 END_2D … … 1216 1207 llkbot = .TRUE. 1217 1208 END IF 1218 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax )1209 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax ) 1219 1210 IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN 1220 1211 ztmp = pu(ji,jj,jk) … … 1264 1255 ! 1265 1256 ! Initialise arrays 1266 l_conv(:,:) = .FALSE. 1267 l_shear(:,:) = .FALSE. 1268 n_ddh(:,:) = 1 1257 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1258 l_conv(ji,jj) = .FALSE. 1259 l_shear(ji,jj) = .FALSE. 1260 n_ddh(ji,jj) = 1 1261 END_2D 1269 1262 ! Initialise INTENT( out) arrays 1270 pwb_ent(:,:) = pp_large 1271 pwb_min(:,:) = pp_large 1263 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1264 pwb_ent(ji,jj) = pp_large 1265 pwb_min(ji,jj) = pp_large 1266 END_2D 1272 1267 ! 1273 1268 ! Determins stability and set flag l_conv 1274 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1269 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1275 1270 IF ( shol(ji,jj) < 0.0_wp ) THEN 1276 1271 l_conv(ji,jj) = .TRUE. … … 1280 1275 END_2D 1281 1276 ! 1282 pshear(:,:) = 0.0_wp 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 ) 1277 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1278 pshear(ji,jj) = 0.0_wp 1279 END_2D 1280 zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D((nn_hls-1))) ) * phbl(A2D((nn_hls-1))) / MAX( sustar(A2D((nn_hls-1))), 1.e-8 ) ) 1281 ! 1282 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1286 1283 IF ( l_conv(ji,jj) ) THEN 1287 1284 IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN … … 1335 1332 ! 1336 1333 ! Calculate entrainment buoyancy flux due to surface fluxes. 1337 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1334 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1338 1335 IF ( l_conv(ji,jj) ) THEN 1339 1336 zwcor = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln … … 1356 1353 END_2D 1357 1354 ! 1358 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1355 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1359 1356 IF ( l_shear(ji,jj) ) THEN 1360 1357 IF ( l_conv(ji,jj) ) THEN … … 1404 1401 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1405 1402 ! 1406 pdtdz(:,:) = pp_large 1407 pdsdz(:,:) = pp_large 1408 pdbdz(:,:) = pp_large 1409 ! 1410 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1403 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1404 pdtdz(ji,jj) = pp_large 1405 pdsdz(ji,jj) = pp_large 1406 pdbdz(ji,jj) = pp_large 1407 END_2D 1408 ! 1409 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1411 1410 IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1412 1411 zthermal = rab_n(ji,jj,1,jp_tem) ! Ideally use nbld not 1?? … … 1455 1454 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1456 1455 ! 1457 pdhdt(:,:) = pp_large 1458 pwb_fk_b(:,:) = pp_large 1459 ! 1460 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1456 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1457 pdhdt(ji,jj) = pp_large 1458 pwb_fk_b(ji,jj) = pp_large 1459 END_2D 1460 ! 1461 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1461 1462 ! 1462 1463 IF ( l_shear(ji,jj) ) THEN … … 1612 1613 REAL(wp) :: zthermal, zbeta 1613 1614 ! 1614 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1615 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1615 1616 IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN 1616 1617 ! … … 1718 1719 REAL, PARAMETER :: pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp ! Also in pycnocline_depth 1719 1720 ! 1720 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )1721 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1721 1722 ! 1722 1723 IF ( l_shear(ji,jj) ) THEN … … 1881 1882 REAL(wp), PARAMETER :: pp_large = -1e10_wp 1882 1883 ! 1883 pdbdz(:,:,:) = pp_large 1884 palpha(:,:) = pp_large 1885 ! 1886 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1884 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1885 pdbdz(ji,jj,:) = pp_large 1886 palpha(ji,jj) = pp_large 1887 END_2D 1888 ! 1889 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1887 1890 ! 1888 1891 IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 1952 1955 ! 1953 1956 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 1954 IF ( iom_use("zdbdz_pyc") ) THEN 1955 osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d ) 1956 END IF 1957 CALL zdf_osm_iomput( "zdbdz_pyc", wmask(A2D(0),:) * pdbdz(A2D(0),:) ) 1957 1958 END IF 1958 1959 ! … … 2004 2005 zb_coup(:,:) = 0.0_wp 2005 2006 ! 2006 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2007 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2007 2008 IF ( l_conv(ji,jj) ) THEN 2008 2009 ! … … 2083 2084 END_2D 2084 2085 ! 2085 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2086 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2086 2087 IF ( l_conv(ji,jj) ) THEN 2087 2088 DO jk = 2, nmld(ji,jj) ! Mixed layer diffusivity … … 2152 2153 ! 2153 2154 END_2D 2154 IF( iom_use("pb_coup") ) THEN 2155 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d ) ! BBL-coupling velocity scale 2156 END IF 2155 CALL zdf_osm_iomput( "pb_coup", tmask(A2D(0),1) * zb_coup(A2D(0)) ) ! BBL-coupling velocity scale 2157 2156 ! 2158 2157 END SUBROUTINE zdf_osm_diffusivity_viscosity … … 2223 2222 jkf_mld = jpk 2224 2223 jkm_mld = 0 2225 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2224 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2226 2225 IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj) 2227 2226 IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj) … … 2238 2237 zsc_ws_1(:,:) = 2.0_wp * swsav(A2D((nn_hls-1))) 2239 2238 ENDWHERE 2240 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )2239 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2241 2240 IF ( l_conv(ji,jj) ) THEN 2242 2241 IF ( jk <= nmld(ji,jj) ) THEN … … 2259 2258 ! 2260 2259 IF ( ln_dia_osm ) THEN 2261 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask * ghamu)2262 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask * ghamv)2260 CALL zdf_osm_iomput( "ghamu_00", wmask(A2D(0),:) * ghamu(A2D(0),:) ) 2261 CALL zdf_osm_iomput( "ghamv_00", wmask(A2D(0),:) * ghamv(A2D(0),:) ) 2263 2262 END IF 2264 2263 ! … … 2278 2277 & ( svstr(A2D((nn_hls-1)))**2 + epsln ) 2279 2278 ENDWHERE 2280 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )2279 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2281 2280 IF ( l_conv(ji,jj) ) THEN 2282 2281 IF ( jk <= nmld(ji,jj) ) THEN … … 2309 2308 zsc_ws_1(:,:) = 0.0_wp 2310 2309 ENDWHERE 2311 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )2310 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2312 2311 IF ( l_conv(ji,jj) ) THEN 2313 2312 IF ( jk <= nmld(ji,jj) ) THEN … … 2330 2329 END IF 2331 2330 END_3D 2332 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2331 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2333 2332 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 2334 2333 ztau_sc_u(ji,jj) = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird * & … … 2350 2349 END IF 2351 2350 END_2D 2352 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )2351 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2353 2352 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2354 2353 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) … … 2371 2370 ! 2372 2371 IF ( ln_dia_osm ) THEN 2373 IF ( iom_use("zwth_ent") ) THEN ! Upward turb. temperature entrainment flux 2374 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d ) 2375 END IF 2376 IF ( iom_use("zws_ent") ) THEN ! Upward turb. salinity entrainment flux 2377 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zws_ent(:,:); CALL iom_put( "zws_ent", osmdia2d ) 2378 END IF 2372 CALL iom_put( "zwth_ent", tmask(A2D(0),1) * zwth_ent(A2D(0)) ) ! Upward turb. temperature entrainment flux 2373 CALL iom_put( "zws_ent", tmask(A2D(0),1) * zws_ent(A2D(0)) ) ! Upward turb. salinity entrainment flux 2379 2374 END IF 2380 2375 ! … … 2388 2383 zsc_uw_1(:,:) = 0.0_wp 2389 2384 ENDWHERE 2390 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )2385 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2391 2386 IF ( l_conv(ji,jj) ) THEN 2392 2387 IF ( jk <= nmld(ji,jj) ) THEN … … 2405 2400 END_3D 2406 2401 ! 2407 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2402 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2408 2403 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 2409 2404 IF ( n_ddh(ji,jj) == 0 ) THEN … … 2424 2419 END IF 2425 2420 END_2D 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.2421 DO_3D_OVR( 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. 2427 2422 IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2428 2423 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) … … 2437 2432 ! 2438 2433 IF ( ln_dia_osm ) THEN 2439 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask * ghamu ) 2440 IF ( iom_use("zsc_uw_1_0") ) THEN 2441 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d ) 2442 END IF 2434 CALL iom_put( "ghamu_0", wmask(A2D(0),:) * ghamu(A2D(0),:) ) 2435 CALL iom_put( "zsc_uw_1_0", tmask(A2D(0),1) * zsc_uw_1(A2D(0)) ) 2443 2436 END IF 2444 2437 ! … … 2457 2450 zsc_ws_1(:,:) = sws0(A2D((nn_hls-1))) 2458 2451 END WHERE 2459 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) )2452 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) ) 2460 2453 IF ( l_conv(ji,jj) ) THEN 2461 2454 IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN … … 2498 2491 & zsc_vw_1(:,:) 2499 2492 ENDWHERE 2500 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )2493 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 2501 2494 IF ( l_conv(ji,jj) ) THEN 2502 2495 IF ( jk <= nmld(ji,jj) ) THEN … … 2530 2523 ! 2531 2524 IF ( ln_dia_osm ) THEN 2532 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask * ghamu ) 2533 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask * ghamv ) 2534 IF ( iom_use("zsc_uw_1_f") ) THEN 2535 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d ) 2536 END IF 2537 IF ( iom_use("zsc_vw_1_f") ) THEN 2538 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d ) 2539 END IF 2540 IF ( iom_use("zsc_uw_2_f") ) THEN 2541 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d ) 2542 END IF 2543 IF ( iom_use("zsc_vw_2_f") ) THEN 2544 osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d ) 2545 END IF 2525 CALL zdf_osm_iomput( "ghamu_f", wmask(A2D(0),:) * ghamu(A2D(0),:) ) 2526 CALL zdf_osm_iomput( "ghamv_f", wmask(A2D(0),:) * ghamv(A2D(0),:) ) 2527 CALL zdf_osm_iomput( "zsc_uw_1_f", tmask(A2D(0),1) * zsc_uw_1(A2D(0)) ) 2528 CALL zdf_osm_iomput( "zsc_vw_1_f", tmask(A2D(0),1) * zsc_vw_1(A2D(0)) ) 2529 CALL zdf_osm_iomput( "zsc_uw_2_f", tmask(A2D(0),1) * zsc_uw_2(A2D(0)) ) 2530 CALL zdf_osm_iomput( "zsc_vw_2_f", tmask(A2D(0),1) * zsc_vw_2(A2D(0)) ) 2546 2531 END IF 2547 2532 ! … … 2551 2536 ! Make surface forced velocity non-gradient terms go to zero at the base 2552 2537 ! of the boundary layer. 2553 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )2538 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2554 2539 IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 2555 2540 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj) ! ALMG to think about … … 2567 2552 ! 2568 2553 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Allocate arrays for output of pycnocline gradient/shear profiles 2569 ALLOCATE( z3ddz_pyc_1( jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat )2554 ALLOCATE( z3ddz_pyc_1(A2D(nn_hls),jpk), z3ddz_pyc_2(A2D(nn_hls),jpk), STAT=istat ) 2570 2555 IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' ) 2571 2556 z3ddz_pyc_1(:,:,:) = 0.0_wp 2572 2557 z3ddz_pyc_2(:,:,:) = 0.0_wp 2573 2558 END IF 2574 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )2559 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2575 2560 IF ( l_conv (ji,jj) ) THEN 2576 2561 ! Unstable conditions. Shouldn;t be needed with no pycnocline code. … … 2637 2622 END_3D 2638 2623 IF ( ln_dia_pyc_scl ) THEN ! Output of pycnocline gradient profiles 2639 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )2640 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )2624 CALL zdf_osm_iomput( "zdtdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_1(A2D(0),:) ) 2625 CALL zdf_osm_iomput( "zdsdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_2(A2D(0),:) ) 2641 2626 END IF 2642 DO_3D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )2627 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 2643 2628 IF ( .NOT. l_conv (ji,jj) ) THEN 2644 2629 IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 2664 2649 END_3D 2665 2650 IF ( ln_dia_pyc_shr ) THEN ! Output of pycnocline shear profiles 2666 IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )2667 IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )2651 CALL zdf_osm_iomput( "zdudz_pyc", wmask(A2D(0),:) * z3ddz_pyc_1(A2D(0),:) ) 2652 CALL zdf_osm_iomput( "zdvdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_2(A2D(0),:) ) 2668 2653 END IF 2669 2654 IF ( ln_dia_osm ) THEN 2670 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask * ghamu)2671 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask * ghamv)2655 CALL zdf_osm_iomput( "ghamu_b", wmask(A2D(0),:) * ghamu(A2D(0),:) ) 2656 CALL zdf_osm_iomput( "ghamv_b", wmask(A2D(0),:) * ghamv(A2D(0),:) ) 2672 2657 END IF 2673 2658 IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN ! Deallocate arrays used for output of pycnocline gradient/shear profiles … … 2675 2660 END IF 2676 2661 ! 2677 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2662 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2678 2663 ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp 2679 2664 ghams(ji,jj,nbld(ji,jj)) = 0.0_wp … … 2683 2668 ! 2684 2669 IF ( ln_dia_osm ) THEN 2685 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask * ghamu ) 2686 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask * ghamv ) 2687 IF ( iom_use("zviscos") ) THEN 2688 osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pviscos; CALL iom_put( "zviscos", osmdia3d ) 2689 END IF 2670 CALL zdf_osm_iomput( "ghamu_1", wmask(A2D(0),:) * ghamu(A2D(0),:) ) 2671 CALL zdf_osm_iomput( "ghamv_1", wmask(A2D(0),:) * ghamv(A2D(0),:) ) 2672 CALL zdf_osm_iomput( "zviscos", wmask(A2D(0),:) * pviscos(A2D(0),:) ) 2690 2673 END IF 2691 2674 ! … … 2693 2676 2694 2677 SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, pmld, pdtdx, pdtdy, pdsdx, & 2695 & pdsdy, pdbd x_mle, pdbdy_mle, pdbds_mle )2678 & pdsdy, pdbds_mle ) 2696 2679 !!---------------------------------------------------------------------- 2697 2680 !! *** ROUTINE zdf_osm_zmld_horizontal_gradients *** … … 2707 2690 !!---------------------------------------------------------------------- 2708 2691 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2709 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pmld ! == Estimated FK BLD used for MLE horizontal gradients == ! 2710 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdtdx ! Horizontal gradient for Fox-Kemper parametrization 2711 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdtdy ! Horizontal gradient for Fox-Kemper parametrization 2712 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdsdx ! Horizontal gradient for Fox-Kemper parametrization 2713 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdsdy ! Horizontal gradient for Fox-Kemper parametrization 2714 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdx_mle ! MLE horiz gradients at u points 2715 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pdbdy_mle ! MLE horiz gradients at v points 2716 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 2692 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT( out) :: pmld ! == Estimated FK BLD used for MLE horizontal gradients == ! 2693 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(inout) :: pdtdx ! Horizontal gradient for Fox-Kemper parametrization 2694 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(inout) :: pdtdy ! Horizontal gradient for Fox-Kemper parametrization 2695 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(inout) :: pdsdx ! Horizontal gradient for Fox-Kemper parametrization 2696 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(inout) :: pdsdy ! Horizontal gradient for Fox-Kemper parametrization 2697 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pdbds_mle ! Magnitude of horizontal buoyancy gradient 2717 2698 ! 2718 2699 ! Local variables 2719 2700 INTEGER :: ji, jj, jk ! Dummy loop indices 2701 INTEGER, DIMENSION(A2D(nn_hls)) :: jk_mld_prof ! Base level of MLE layer 2720 2702 INTEGER :: ikt, ikmax ! Local integers 2721 2703 REAL(wp) :: zc … … 2731 2713 ! 2732 2714 ! == MLD used for MLE ==! 2733 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 2734 pmld(:,:) = 0.0_wp ! Here hmlp used as a dummy variable, integrating vertically N^2 2715 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2716 jk_mld_prof(ji,jj) = nlb10 ! Initialization to the number of w ocean point 2717 pmld(ji,jj) = 0.0_wp ! Here hmlp used as a dummy variable, integrating vertically N^2 2718 END_2D 2735 2719 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! Convert density criteria into N^2 criteria 2736 2720 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) 2737 2721 ikt = mbkt(ji,jj) 2738 2722 pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm) 2739 IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level2723 IF( pmld(ji,jj) < zN2_c ) jk_mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2740 2724 END_3D 2741 2725 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2742 mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) ) ! Ensure mld_prof .ge. nbld 2743 pmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 2744 END_2D 2745 ! 2746 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! Max level of the computation 2726 jk_mld_prof(ji,jj) = MAX( jk_mld_prof(ji,jj), nbld(ji,jj) ) ! Ensure jk_mld_prof .ge. nbld 2727 pmld(ji,jj) = gdepw(ji,jj,jk_mld_prof(ji,jj),Kmm) 2728 END_2D 2729 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2730 mld_prof(ji,jj) = jk_mld_prof(ji,jj) 2731 END_2D 2732 ! 2733 ikmax = MIN( MAXVAL( jk_mld_prof(A2D(nn_hls)) ), jpkm1 ) ! Max level of the computation 2747 2734 ztm(:,:) = 0.0_wp 2748 2735 zsm(:,:) = 0.0_wp 2749 2736 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax ) 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-points2737 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, jk_mld_prof(ji,jj) - jk ), 1 ), KIND=wp ) ! zc being 0 outside the ML t-points 2751 2738 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) 2752 2739 zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm) 2753 2740 END_3D 2754 2741 ! Average temperature and salinity 2755 ztm(:,:) = ztm(:,:) / MAX( e3t( :,:,1,Kmm), pmld(:,:) )2756 zsm(:,:) = zsm(:,:) / MAX( e3t( :,:,1,Kmm), pmld(:,:) )2742 ztm(:,:) = ztm(:,:) / MAX( e3t(A2D(nn_hls),1,Kmm), pmld(A2D(nn_hls)) ) 2743 zsm(:,:) = zsm(:,:) / MAX( e3t(A2D(nn_hls),1,Kmm), pmld(A2D(nn_hls)) ) 2757 2744 ! Calculate horizontal gradients at u & v points 2758 2745 zmld_midu(:,:) = 0.0_wp … … 2776 2763 CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm ) 2777 2764 CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm ) 2778 DO_2D ( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )2779 pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) )2780 END_2D 2781 DO_2D ( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )2782 pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) )2765 DO_2D_OVR( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 2766 dbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) ) 2767 END_2D 2768 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 2769 dbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) ) 2783 2770 END_2D 2784 2771 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 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 ) + &2786 & pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) )2772 pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji, jj) * dbdx_mle(ji, jj) + dbdy_mle(ji,jj ) * dbdy_mle(ji,jj ) + & 2773 & dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 2787 2774 END_2D 2788 2775 ! … … 2826 2813 REAL(wp) :: zdbdz_mle_int 2827 2814 ! 2828 znd_param( A2D((nn_hls-1))) = 0.0_wp2829 ! 2830 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2815 znd_param(:,:) = 0.0_wp 2816 ! 2817 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2831 2818 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2832 2819 pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj) 2833 2820 END_2D 2834 2821 ! 2835 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2822 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2836 2823 ! 2837 2824 IF ( l_conv(ji,jj) ) THEN … … 2861 2848 ! 2862 2849 ! Diagnosis 2863 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2850 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2864 2851 ! 2865 2852 IF ( l_conv(ji,jj) ) THEN … … 2922 2909 !!---------------------------------------------------------------------- 2923 2910 INTEGER, INTENT(in ) :: Kmm ! Time-level index 2924 REAL(wp), DIMENSION( :,:),INTENT(in ) :: pmld ! == Estimated FK BLD used for MLE horiz gradients == !2911 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in ) :: pmld ! == Estimated FK BLD used for MLE horiz gradients == ! 2925 2912 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: phmle ! MLE depth 2926 2913 REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) :: pvel_mle ! Velocity scale for dhdt with stable ML and FK … … 2942 2929 ! 2943 2930 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE 2944 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2931 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2945 2932 IF ( l_conv(ji,jj) ) THEN 2946 2933 ztmp = r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf … … 2951 2938 END_2D 2952 2939 ! Timestep mixed layer eddy depth 2953 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2940 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2954 2941 IF ( l_mle(ji,jj) ) THEN ! MLE layer growing 2955 2942 ! Buoyancy gradient at base of MLE layer … … 2973 2960 END_2D 2974 2961 ! 2975 mld_prof(:,:) = 4 2976 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 2962 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 2977 2963 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk ) 2978 2964 END_3D 2979 DO_2D ( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )2965 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2980 2966 phmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 2981 2967 END_2D … … 3432 3418 END SUBROUTINE dyn_osm 3433 3419 3420 SUBROUTINE zdf_osm_iomput_2d( cdname, posmdia2d ) 3421 !!---------------------------------------------------------------------- 3422 !! *** ROUTINE zdf_osm_iomput_2d *** 3423 !! 3424 !! ** Purpose : Wrapper for subroutine iom_put that accepts 2D arrays 3425 !! with and without halo 3426 !! 3427 !!---------------------------------------------------------------------- 3428 CHARACTER(LEN=*), INTENT(in ) :: cdname 3429 REAL(wp), DIMENSION(:,:), INTENT(in ) :: posmdia2d 3430 ! 3431 IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN 3432 IF ( SIZE( posmdia2d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia2d, 2 ) == ntej-ntsj+1 ) THEN ! Halo absent 3433 osmdia2d(A2D(0)) = posmdia2d(:,:) 3434 CALL iom_put( cdname, osmdia2d(A2D(nn_hls)) ) 3435 ELSE ! Halo present 3436 CALL iom_put( cdname, osmdia2d ) 3437 END IF 3438 END IF 3439 ! 3440 END SUBROUTINE zdf_osm_iomput_2d 3441 3442 SUBROUTINE zdf_osm_iomput_3d( cdname, posmdia3d ) 3443 !!---------------------------------------------------------------------- 3444 !! *** ROUTINE zdf_osm_iomput_3d *** 3445 !! 3446 !! ** Purpose : Wrapper for subroutine iom_put that accepts 3D arrays 3447 !! with and without halo 3448 !! 3449 !!---------------------------------------------------------------------- 3450 CHARACTER(LEN=*), INTENT(in ) :: cdname 3451 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: posmdia3d 3452 ! 3453 IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN 3454 IF ( SIZE( posmdia3d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia3d, 2 ) == ntej-ntsj+1 ) THEN ! Halo absent 3455 osmdia3d(A2D(0),:) = posmdia3d(:,:,:) 3456 CALL iom_put( cdname, osmdia3d(A2D(nn_hls),:) ) 3457 ELSE ! Halo present 3458 CALL iom_put( cdname, osmdia3d ) 3459 END IF 3460 END IF 3461 ! 3462 END SUBROUTINE zdf_osm_iomput_3d 3463 3434 3464 !!====================================================================== 3435 3465 -
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfphy.F90
r14856 r14900 186 186 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 187 187 IF( lk_top .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 188 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled)189 IF( ln_tile .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis does not yet work with tiling' )190 188 IF( lk_top .AND. ln_zdfmfc ) CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' ) 191 189 IF(lwp) THEN … … 256 254 INTEGER :: ji, jj, jk ! dummy loop indice 257 255 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zsh2 ! shear production 258 ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled)259 LOGICAL :: lskip260 256 !! --------------------------------------------------------------------- 261 257 ! 262 258 IF( ln_timing ) CALL timing_start('zdf_phy') 263 264 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled)265 lskip = .FALSE.266 267 IF( ln_tile .AND. nzdf_phy == np_OSM ) THEN268 IF( ntile == 1 ) THEN269 CALL dom_tile_stop( ldhold=.TRUE. )270 ELSE271 lskip = .TRUE.272 ENDIF273 ENDIF274 259 ! 275 260 IF( l_zdfdrg ) THEN !== update top/bottom drag ==! (non-linear cases) … … 301 286 ! 302 287 CALL zdf_mxl( kt, Kmm ) !* mixed layer depth, and level 303 304 ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled) 305 IF( .NOT. lskip ) THEN 306 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 307 ! 308 ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2. 309 IF( l_zdfsh2 ) THEN !* shear production at w-points (energy conserving form) 310 IF( ln_tile ) THEN 311 IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:) ! Preserve "now" avm_k for calculation of zsh2 312 CALL zdf_sh2( Kbb, Kmm, avm_k_n, & ! <<== in 313 & zsh2 ) ! ==>> out : shear production 314 ELSE 315 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 316 & zsh2 ) ! ==>> out : shear production 317 ENDIF 318 ENDIF 319 ! 320 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 321 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 322 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 323 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 324 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 288 ! 289 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 290 ! 291 ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2. 292 IF( l_zdfsh2 ) THEN !* shear production at w-points (energy conserving form) 293 IF( ln_tile ) THEN 294 IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:) ! Preserve "now" avm_k for calculation of zsh2 295 CALL zdf_sh2( Kbb, Kmm, avm_k_n, & ! <<== in 296 & zsh2 ) ! ==>> out : shear production 297 ELSE 298 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 299 & zsh2 ) ! ==>> out : shear production 300 ENDIF 301 ENDIF 302 ! 303 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 304 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 305 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 306 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 307 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 325 308 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 326 309 ! ! avt_k and avm_k set one for all at initialisation phase 327 310 !!gm avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 328 311 !!gm avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 329 END SELECT 330 331 IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. ) 332 ENDIF 312 END SELECT 333 313 ! 334 314 ! !== ocean Kz ==! (avt, avs, avm)
Note: See TracChangeset
for help on using the changeset viewer.
