- Timestamp:
- 2021-03-25T12:51:33+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90
r14574 r14636 400 400 zz0 = rn_abs ! surface equi-partition in 2-bands 401 401 zz1 = 1. - rn_abs 402 DO_2D( 0, 0, 0, 0 ) 402 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 403 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 403 404 ! Surface downward irradiance (so always +ve) 404 405 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp … … 410 411 END_2D 411 412 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 412 DO_2D( 0, 0, 0, 0 ) 413 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 414 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 413 415 zthermal = rab_n(ji,jj,1,jp_tem) 414 416 zbeta = rab_n(ji,jj,1,jp_sal) … … 437 439 ! Assume constant La#=0.3 438 440 CASE(0) 439 DO_2D( 0, 0, 0, 0 ) 441 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 442 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 440 443 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 441 444 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 … … 446 449 ! Assume Pierson-Moskovitz wind-wave spectrum 447 450 CASE(1) 448 DO_2D( 0, 0, 0, 0 ) 451 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 452 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 449 453 ! Use wind speed wndm included in sbc_oce module 450 454 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) … … 455 459 zfac = 2.0_wp * rpi / 16.0_wp 456 460 457 DO_2D( 0, 0, 0, 0 ) 461 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 462 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 458 463 IF (hsw(ji,jj) > 1.e-4) THEN 459 464 ! Use wave fields … … 472 477 IF (ln_zdfosm_ice_shelter) THEN 473 478 ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 474 DO_2D( 0, 0, 0, 0 ) 479 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 480 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 475 481 zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 476 482 dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) … … 494 500 z_two_thirds = 2.0_wp / 3.0_wp 495 501 496 DO_2D( 0, 0, 0, 0 ) 502 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 503 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 497 504 zthickness = rn_osm_hblfrac*hbl(ji,jj) 498 505 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) … … 509 516 zsqrtpi = SQRT(rpi) 510 517 511 DO_2D( 0, 0, 0, 0 ) 518 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 519 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 512 520 zthickness = rn_osm_hblfrac*hbl(ji,jj) 513 521 z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) … … 530 538 ! Langmuir velocity scale (zwstrl), La # (zla) 531 539 ! mixed scale (zvstr), convective velocity scale (zwstrc) 532 DO_2D( 0, 0, 0, 0 ) 540 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 541 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 533 542 ! Langmuir velocity scale (zwstrl), at T-point 534 543 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird … … 563 572 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 564 573 ibld(:,:) = 4 565 DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 574 ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 575 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 ) 566 576 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 567 577 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) … … 570 580 ! ########################################################################## 571 581 572 DO_2D( 0, 0, 0, 0 ) 582 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 583 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 573 584 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 574 585 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 )) … … 590 601 ! Fox-Kemper Scheme 591 602 mld_prof = 4 592 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 603 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 604 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 593 605 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 594 606 END_3D … … 596 608 CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle) 597 609 598 DO_2D( 0, 0, 0, 0 ) 610 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 611 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 599 612 zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 600 613 END_2D … … 611 624 lflux(:,:) = .FALSE. 612 625 lmle(:,:) = .FALSE. 613 DO_2D( 0, 0, 0, 0 ) 626 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 627 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 614 628 IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE. 615 629 END_2D … … 617 631 618 632 ! Test if pycnocline well resolved 619 DO_2D( 0, 0, 0, 0 ) 633 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 634 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 620 635 IF (lconv(ji,jj) ) THEN 621 636 ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm) … … 638 653 ! Rate of change of hbl 639 654 CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt ) 640 DO_2D( 0, 0, 0, 0 ) 655 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 656 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 641 657 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need ww here, so subtract it 642 658 ! adjustment to represent limiting by ocean bottom … … 650 666 ibld(:,:) = 4 651 667 652 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 668 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 669 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 ) 653 670 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 654 671 ibld(ji,jj) = jk … … 669 686 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 670 687 671 DO_2D( 0, 0, 0, 0 ) 688 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 689 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 672 690 IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE. 673 691 END_2D … … 709 727 710 728 711 DO_2D( 0, 0, 0, 0 ) 729 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 730 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 712 731 IF ( lconv(ji,jj) ) THEN 713 732 DO jk = 2, imld(ji,jj) … … 744 763 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 745 764 END IF 746 DO_2D( 0, 0, 0, 0 ) 765 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 766 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 747 767 IF ( lconv(ji,jj) ) THEN 748 768 DO jk = 2, imld(ji,jj) … … 776 796 ENDWHERE 777 797 778 DO_2D( 0, 0, 0, 0 ) 798 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 799 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 779 800 IF (lconv(ji,jj) ) THEN 780 801 DO jk = 2, imld(ji,jj) … … 838 859 ENDWHERE 839 860 840 DO_2D( 0, 0, 0, 0 ) 861 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 862 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 841 863 IF ( lconv(ji,jj) ) THEN 842 864 DO jk = 2 , imld(ji,jj) … … 856 878 END_2D 857 879 858 DO_2D( 0, 0, 0, 0 ) 880 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 881 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 859 882 IF ( lpyc(ji,jj) ) THEN 860 883 IF ( j_ddh(ji,jj) == 0 ) THEN … … 891 914 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 892 915 893 DO_2D( 1, 0, 1, 0 ) 916 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 917 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 894 918 895 919 IF ( lconv(ji,jj) ) THEN … … 907 931 END_2D 908 932 909 DO_2D( 0, 0, 0, 0 ) 933 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 934 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 910 935 IF ( lconv(ji,jj) ) THEN 911 936 DO jk = 2, imld(ji,jj) … … 954 979 ENDWHERE 955 980 956 DO_2D( 0, 0, 0, 0 ) 981 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 982 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 957 983 IF ( lconv(ji,jj) ) THEN 958 984 DO jk = 2, imld(ji,jj) … … 1001 1027 ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer. 1002 1028 1003 DO_2D( 0, 0, 0, 0 ) 1029 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1030 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1004 1031 IF ( .not. lconv(ji,jj) ) THEN 1005 1032 DO jk = 2, ibld(ji,jj) … … 1017 1044 1018 1045 ! pynocline contributions 1019 DO_2D( 0, 0, 0, 0 ) 1046 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1047 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1020 1048 IF ( .not. lconv(ji,jj) ) THEN 1021 1049 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 1035 1063 END IF 1036 1064 1037 DO_2D( 0, 0, 0, 0 ) 1065 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1066 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1038 1067 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1039 1068 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp … … 1060 1089 ! rotate non-gradient velocity terms back to model reference frame 1061 1090 1062 DO_2D( 0, 0, 0, 0 ) 1091 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1092 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1063 1093 DO jk = 2, ibld(ji,jj) 1064 1094 ztemp = ghamu(ji,jj,jk) … … 1076 1106 ! KPP-style Ri# mixing 1077 1107 IF( ln_kpprimix) THEN 1078 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) !* Shear production at uw- and vw-points (energy conserving form) 1108 1109 ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 2, jpkm1 ) !* Shear production at uw- and vw-points (energy conserving form) 1110 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 ) 1079 1111 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 1080 1112 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & … … 1085 1117 END_3D 1086 1118 ! 1087 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1119 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1120 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 1088 1121 ! ! shear prod. at w-point weightened by mask 1089 1122 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & … … 1096 1129 END_3D 1097 1130 1098 DO_2D( 0, 0, 0, 0 ) 1131 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1132 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1099 1133 DO jk = ibld(ji,jj) + 1, jpkm1 1100 1134 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri … … 1107 1141 ! KPP-style set diffusivity large if unstable below BL 1108 1142 IF( ln_convmix) THEN 1109 DO_2D( 0, 0, 0, 0 ) 1143 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1144 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1110 1145 DO jk = ibld(ji,jj) + 1, jpkm1 1111 1146 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv … … 1117 1152 1118 1153 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1119 DO_2D( 0, 0, 0, 0 ) 1154 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1155 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1120 1156 IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer 1121 1157 ! Calculate MLE flux contribution from surface fluxes … … 1158 1194 ! GN 25/8: need to change tmask --> wmask 1159 1195 1160 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1196 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1197 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 1161 1198 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1162 1199 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 1163 1200 END_3D 1164 1201 ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids 1165 CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, &1166 & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp )1202 IF (nn_hls.eq.1) CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & 1203 & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) 1167 1204 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1168 1205 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & … … 1176 1213 END_3D 1177 1214 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1178 CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1215 ! [comm_cleanup] ! no need lbc_lnk for output 1216 ! CALL lbc_lnk( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1179 1217 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1180 1218 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign changed) 1181 CALL lbc_lnk( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & 1182 & ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) 1219 ! [comm_cleanup] ! no need lbc_lnk for output 1220 ! CALL lbc_lnk( 'zdfosm', ghamt, 'W', 1.0_wp , ghams, 'W', 1.0_wp, & 1221 ! & ghamu, 'U', -1.0_wp , ghamv, 'V', -1.0_wp ) 1183 1222 1184 1223 IF(ln_dia_osm) THEN … … 1280 1319 REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15 1281 1320 1282 DO_2D( 0, 0, 0, 0 ) 1321 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1322 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1283 1323 IF ( lconv(ji,jj) ) THEN 1284 1324 … … 1323 1363 END_2D 1324 1364 ! 1325 DO_2D( 0, 0, 0, 0 ) 1365 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1366 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1326 1367 IF ( lconv(ji,jj) ) THEN 1327 1368 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity … … 1422 1463 1423 1464 ! Determins stability and set flag lconv 1424 DO_2D( 0, 0, 0, 0 ) 1465 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1466 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1425 1467 IF ( zhol(ji,jj) < 0._wp ) THEN 1426 1468 lconv(ji,jj) = .TRUE. … … 1439 1481 j_ddh(:,:) = 1 1440 1482 1441 DO_2D( 0, 0, 0, 0 ) 1483 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1484 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1442 1485 IF ( lconv(ji,jj) ) THEN 1443 1486 IF ( zdb_bl(ji,jj) > 0._wp ) THEN … … 1476 1519 ! Calculate entrainment buoyancy flux due to surface fluxes. 1477 1520 1478 DO_2D( 0, 0, 0, 0 ) 1521 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1522 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1479 1523 IF ( lconv(ji,jj) ) THEN 1480 1524 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln … … 1501 1545 zwb_min(:,:) = 0._wp 1502 1546 1503 DO_2D( 0, 0, 0, 0 ) 1547 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1548 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1504 1549 IF ( lshear(ji,jj) ) THEN 1505 1550 IF ( lconv(ji,jj) ) THEN … … 1562 1607 zu = 0._wp 1563 1608 zv = 0._wp 1564 DO_2D( 0, 0, 0, 0 ) 1609 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1610 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1565 1611 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1566 1612 zbeta = rab_n(ji,jj,1,jp_sal) … … 1619 1665 REAL(wp) :: ztemp 1620 1666 1621 DO_2D( 0, 0, 0, 0 ) 1667 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1668 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1622 1669 ztemp = zu(ji,jj) 1623 1670 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) … … 1653 1700 znd_param(:,:) = 0._wp 1654 1701 1655 DO_2D( 0, 0, 0, 0 ) 1702 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1703 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1656 1704 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1657 1705 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj) 1658 1706 END_2D 1659 DO_2D( 0, 0, 0, 0 ) 1707 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1708 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1660 1709 ! 1661 1710 IF ( lconv(ji,jj) ) THEN … … 1681 1730 1682 1731 ! Diagnosis 1683 DO_2D( 0, 0, 0, 0 ) 1732 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1733 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1684 1734 IF ( lconv(ji,jj) ) THEN 1685 1735 zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & … … 1751 1801 1752 1802 1753 DO_2D( 0, 0, 0, 0 ) 1803 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1804 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1754 1805 IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 1755 1806 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? … … 1781 1832 REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15 1782 1833 1783 DO_2D( 0, 0, 0, 0 ) 1834 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1835 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1784 1836 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 1785 1837 IF ( lconv(ji,jj) ) THEN ! convective conditions … … 1866 1918 REAL(wp) :: zzeta_v = 0.45 1867 1919 ! 1868 DO_2D( 0, 0, 0, 0 ) 1920 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1921 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1869 1922 ! 1870 1923 IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN … … 1929 1982 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 1930 1983 1931 DO_2D( 0, 0, 0, 0 ) 1984 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 1985 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1932 1986 1933 1987 IF ( lshear(ji,jj) ) THEN … … 2072 2126 REAL(wp) :: zthermal, zbeta 2073 2127 2074 DO_2D( 0, 0, 0, 0 ) 2128 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2129 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2075 2130 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 2076 2131 ! … … 2176 2231 REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth 2177 2232 2178 DO_2D( 0, 0, 0, 0 ) 2233 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2234 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2179 2235 2180 2236 IF ( lshear(ji,jj) ) THEN … … 2322 2378 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 2323 2379 zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria 2324 DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 2380 ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 2381 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) 2325 2382 ikt = mbkt(ji,jj) 2326 2383 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 2327 2384 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 2328 2385 END_3D 2329 DO_2D( 1, 1, 1, 1 ) 2386 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 2387 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2330 2388 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 2331 2389 zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) … … 2337 2395 ztm(:,:) = 0._wp 2338 2396 zsm(:,:) = 0._wp 2339 DO_3D( 1, 1, 1, 1, 1, ikmax ) 2397 ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, ikmax ) 2398 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax ) 2340 2399 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 2341 2400 ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) … … 2347 2406 ! calculate horizontal gradients at u & v points 2348 2407 2349 DO_2D( 1, 0, 0, 0 ) 2408 ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 ) 2409 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 2350 2410 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 2351 2411 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) … … 2355 2415 END_2D 2356 2416 2357 DO_2D( 0, 0, 1, 0 ) 2417 ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 ) 2418 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 2358 2419 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 2359 2420 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) … … 2366 2427 CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm) 2367 2428 2368 DO_2D( 1, 0, 0, 0 ) 2429 ! [comm_cleanup] ! DO_2D( 1, 0, 0, 0 ) 2430 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 2369 2431 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 2370 2432 END_2D 2371 DO_2D( 0, 0, 1, 0 ) 2433 ! [comm_cleanup] ! DO_2D( 0, 0, 1, 0 ) 2434 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 2372 2435 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 2373 2436 END_2D 2374 2437 2375 DO_2D( 0, 0, 0, 0 ) 2438 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2439 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2376 2440 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 2377 2441 zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & … … 2400 2464 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 2401 2465 2402 DO_2D( 0, 0, 0, 0 ) 2466 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2467 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2403 2468 IF ( lconv(ji,jj) ) THEN 2404 2469 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf … … 2409 2474 END_2D 2410 2475 ! Timestep mixed layer eddy depth. 2411 DO_2D( 0, 0, 0, 0 ) 2476 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2477 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2412 2478 IF ( lmle(ji,jj) ) THEN ! MLE layer growing. 2413 2479 ! Buoyancy gradient at base of MLE layer. … … 2433 2499 2434 2500 mld_prof = 4 2435 DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 2501 2502 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 2503 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 2436 2504 IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 2437 2505 END_3D 2438 DO_2D( 0, 0, 0, 0 ) 2506 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 2507 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 2439 2508 zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm) 2440 2509 END_2D … … 2585 2654 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2586 2655 z1_t2 = 2.e-5 2587 DO_2D( 1, 1, 1, 1 ) 2656 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 2657 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2588 2658 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2589 2659 END_2D … … 2630 2700 etmean(:,:,:) = 0.e0 2631 2701 2632 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2702 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2703 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 2633 2704 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 2634 2705 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & … … 2644 2715 etmean(:,:,:) = 0.e0 2645 2716 2646 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2717 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2718 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 2647 2719 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 2648 2720 & / MAX( 1., 2.* tmask(ji,jj,jk) & … … 2759 2831 ! 2760 2832 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 2761 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 2833 ! [comm_cleanup] ! DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 2834 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 2762 2835 ikt = mbkt(ji,jj) 2763 2836 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) … … 2765 2838 END_3D 2766 2839 ! 2767 DO_2D( 1, 1, 1, 1 ) 2840 ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 2841 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 2768 2842 iiki = MAX(4,imld_rst(ji,jj)) 2769 2843 hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth … … 2812 2886 ENDIF 2813 2887 2814 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2888 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 2889 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 2815 2890 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 2816 2891 & - ( ghamt(ji,jj,jk ) & … … 2879 2954 !code saving tracer trends removed, replace with trdmxl_oce 2880 2955 2881 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 2956 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 2957 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! add non-local u and v fluxes 2882 2958 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 2883 2959 & - ( ghamu(ji,jj,jk ) &
Note: See TracChangeset
for help on using the changeset viewer.
