Changeset 14554
- Timestamp:
- 2021-02-26T20:30:58+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
r14551 r14554 285 285 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 286 286 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 287 REAL(wp), DIMENSION(jpi,jpj) :: zddhdt ! correction to dhdt due to internal structure.288 287 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients 289 288 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients … … 299 298 REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient. 300 299 ! Flux-gradient relationship variables 301 REAL(wp), DIMENSION(jpi, jpj) :: zshear , zri_i ! Shear production and interfacial richardon number.300 REAL(wp), DIMENSION(jpi, jpj) :: zshear ! Shear production 302 301 303 302 REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep … … 370 369 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 371 370 372 zddhdt(:,:) = 0._wp373 371 ! hbl = MAX(hbl,epsln) 374 372 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 560 558 ! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1 561 559 CALL zdf_osm_vertical_average( Kbb, Kmm, & 562 & i bld, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, ibld-imld+1, &560 & imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, ibld-imld+1, & 563 561 & zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 564 562 ! Velocity components in frame aligned with surface stress. … … 566 564 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 567 565 ! Determine the state of the OSBL, stable/unstable, shear/no shear 568 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear , zri_i)566 CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 569 567 570 568 IF ( ln_osm_mle ) THEN … … 586 584 CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 587 585 CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 588 CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )586 CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 589 587 ELSE ! ln_osm_mle 590 588 ! FK not selected, Boundary Layer only. … … 622 620 & ibld-imld+1, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 623 621 ! Rate of change of hbl 624 CALL zdf_osm_calculate_dhdt( zdhdt , zddhdt)622 CALL zdf_osm_calculate_dhdt( zdhdt ) 625 623 DO_2D( 0, 0, 0, 0 ) 626 624 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 … … 685 683 ! Calculate non-gradient components of the flux-gradient relationships 686 684 ! -------------------------------------------------------------------- 687 CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, ibld_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear,&685 CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear, & 688 686 & zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla, & 689 687 & zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, & … … 928 926 zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj) 929 927 930 IF ( lshear(ji,jj) . and. j_ddh(ji,jj) == 1) THEN928 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 931 929 zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj) 932 930 ENDIF … … 938 936 zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj) 939 937 zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac 940 IF ( lshear(ji,jj) . and. j_ddh(ji,jj) == 1) THEN938 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN 941 939 zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj) 942 940 ENDIF … … 944 942 zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) ) 945 943 zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac 946 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) )944 zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) ) 947 945 948 946 zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third … … 1022 1020 END SUBROUTINE zdf_osm_diffusivity_viscosity 1023 1021 1024 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear , zri_i)1022 SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear ) 1025 1023 1026 1024 !!--------------------------------------------------------------------- … … 1039 1037 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer. 1040 1038 REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline 1041 REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number1042 1039 1043 1040 ! Local Variables … … 1046 1043 1047 1044 REAL(wp), DIMENSION(jpi,jpj) :: zekman 1048 REAL(wp) :: zri_p, zri_b ! Richardson numbers1045 REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers 1049 1046 REAL(wp) :: zshear_u, zshear_v, zwb_shr 1050 1047 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1051 1048 1052 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1 1053 REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59 1054 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04 1049 REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8 1050 REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03 1055 1051 REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15 1056 1052 REAL, PARAMETER :: rn_ri_p_thresh = 27.0 1053 REAL, PARAMETER :: zri_c = 0.25 1054 REAL, PARAMETER :: zek = 4.0 1057 1055 REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress. 1058 1056 … … 1067 1065 END_2D 1068 1066 1069 zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1070 1071 WHERE ( lconv ) 1072 zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 ) 1073 END WHERE 1067 zekman(:,:) = EXP( -1.0_wp * zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) ) 1074 1068 1075 1069 zshear(:,:) = 0._wp … … 1079 1073 IF ( lconv(ji,jj) ) THEN 1080 1074 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1081 zri_p = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 1082 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1083 1084 zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 ) 1085 1086 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 1075 zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 & 1076 & / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp ) 1077 1078 IF ( ff_t(ji,jj) >= 0.0_wp ) THEN 1079 ! Northern hemisphere 1080 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( -1.0_wp * zdv_ml(ji,jj), 1e-5_wp)**2 ) 1081 ELSE 1082 ! Southern hemisphere 1083 zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( zdv_ml(ji,jj), 1e-5_wp)**2 ) 1084 END IF 1085 zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) ) 1086 ! Stability dependence 1087 zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) ) 1087 1088 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1088 1089 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when ! 1089 1090 ! full code available ! 1090 1091 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1091 IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN 1092 IF ( zshear(ji,jj) > 1e-10 ) THEN 1093 IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN 1092 1094 ! Growing shear layer 1093 j_ddh(ji,jj) = 01094 lshear(ji,jj) = .TRUE.1095 ELSE1096 j_ddh(ji,jj) = 11097 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN1095 j_ddh(ji,jj) = 0 1096 lshear(ji,jj) = .TRUE. 1097 ELSE 1098 j_ddh(ji,jj) = 1 1099 ! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN 1098 1100 ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer. 1099 lshear(ji,jj) = .TRUE. 1100 ELSE 1101 lshear(ji,jj) = .TRUE. 1102 ! ELSE 1103 END IF 1104 ELSE 1105 j_ddh(ji,jj) = 2 1106 lshear(ji,jj) = .FALSE. 1107 END IF 1101 1108 ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline. 1102 zshear(ji,jj) = 0.5 * zshear(ji,jj) 1103 lshear(ji,jj) = .FALSE. 1104 ENDIF 1105 ENDIF 1109 ! zshear(ji,jj) = 0.5 * zshear(ji,jj) 1110 ! lshear(ji,jj) = .FALSE. 1111 ! ENDIF 1106 1112 ELSE ! zdb_bl test, note zshear set to zero 1107 1113 j_ddh(ji,jj) = 2 … … 1128 1134 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1129 1135 END IF 1130 zwb_ent(ji,jj) = - 2.0 * 0.2* zrf_conv * zwbav(ji,jj) &1131 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &1132 & + zr_stokes * ( 0.15 * EXP( -1.5* zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &1133 & - zrf_langmuir * 0.03* zwstrl(ji,jj)**3 ) / zhml(ji,jj)1136 zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * zwbav(ji,jj) & 1137 & - zalpha_s * zrf_shear * zustar(ji,jj)**3 / zhml(ji,jj) & 1138 & + zr_stokes * ( zalpha_s * EXP( -1.5_wp * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1139 & - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1134 1140 ! 1135 1141 ENDIF … … 1142 1148 IF ( lconv(ji,jj) ) THEN 1143 1149 ! Unstable OSBL 1144 zwb_shr = - za_wb_s* zshear(ji,jj)1150 zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * zshear(ji,jj) 1145 1151 IF ( j_ddh(ji,jj) == 0 ) THEN 1146 1152 1147 ! Developing shear layer, additional shear production possible. 1148 1149 zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) 1150 zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) ) 1151 zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1152 1153 zwb_shr = -za_wb_s * zshear(ji,jj) 1153 ! ! Developing shear layer, additional shear production possible. 1154 1155 ! zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp ) 1156 ! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 ) 1157 ! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u ) 1158 1159 ! zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 ) 1160 ! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u ) 1154 1161 1155 1162 ENDIF 1156 1163 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1157 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)1164 ! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1158 1165 ELSE ! IF ( lconv ) THEN - ENDIF 1159 1166 ! Stable OSBL - shear production not coded for first attempt. 1160 1167 ENDIF ! lconv 1161 E LSE! lshear1162 1168 END IF ! lshear 1169 IF ( lconv(ji,jj) ) THEN 1163 1170 ! Unstable OSBL 1164 zwb_shr = -za_wb_s * zshear(ji,jj) 1165 zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr 1166 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1167 ENDIF ! lconv 1168 ENDIF ! lshear 1171 zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj) 1172 END IF ! lconv 1169 1173 END_2D 1170 1174 IF( ln_timing ) CALL timing_stop('zdf_osm_os') … … 1196 1200 ztemp = zdu(ji,jj) 1197 1201 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1198 zdv(ji,jj) = zdv(ji,jj) * z sin_w(ji,jj) - ztemp * zsin_w(ji,jj)1202 zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1199 1203 END_2D 1200 1204 IF( ln_timing ) CALL timing_stop('zdf_osm_vr') … … 1221 1225 REAL(wp), DIMENSION(jpi,jpj) :: znd_param 1222 1226 REAL(wp) :: zbuoy, ztmp, zpe_mle_layer 1223 REAL(wp) :: zpe_mle_ref, z wb_ent, zdbdz_mle_int1227 REAL(wp) :: zpe_mle_ref, zdbdz_mle_int 1224 1228 1225 1229 IF( ln_timing ) CALL timing_start('zdf_osm_osf') … … 1258 1262 DO_2D( 0, 0, 0, 0 ) 1259 1263 IF ( lconv(ji,jj) ) THEN 1260 zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) & 1261 & - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) & 1262 & + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 & 1263 & - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1264 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN 1264 IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN 1265 1265 IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN 1266 1266 ! MLE layer growing … … 1335 1335 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1336 1336 zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & 1337 & / e3t(ji,jj, ibld(ji,jj),Kmm)1337 & / e3t(ji,jj,jkb,Kmm) 1338 1338 zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & 1339 & / e3t(ji,jj, ibld(ji,jj),Kmm)1339 & / e3t(ji,jj,jkb,Kmm) 1340 1340 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1341 1341 ELSE … … 1377 1377 zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj) 1378 1378 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1379 DO jk = 2, ibld(ji,jj) +ibld_ext1379 DO jk = 2, ibld(ji,jj) 1380 1380 znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp 1381 1381 IF ( znd <= zzeta_m ) THEN … … 1426 1426 END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 1427 1427 1428 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt , zddhdt)1428 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) 1429 1429 !!--------------------------------------------------------------------- 1430 1430 !! *** ROUTINE zdf_osm_calculate_dhdt *** … … 1436 1436 !!---------------------------------------------------------------------- 1437 1437 1438 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt , zddhdt! Rate of change of hbl1438 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl 1439 1439 1440 1440 INTEGER :: jj, ji 1441 1441 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 1442 REAL(wp) :: zvel_max !, zwb_min1443 REAL(wp) :: zzeta_m = 0.31444 REAL(wp) :: zgamma_c = 2.01445 REAL(wp) :: zdhoh = 0.11446 REAL(wp) :: alpha_bc = 0.51447 REAL , PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth1442 REAL(wp) :: zvel_max, zddhdt 1443 REAL(wp), PARAMETER :: zzeta_m = 0.3_wp 1444 REAL(wp), PARAMETER :: zgamma_c = 2.0_wp 1445 REAL(wp), PARAMETER :: zdhoh = 0.1_wp 1446 REAL(wp), PARAMETER :: zalpha_b = 0.3_wp 1447 REAL(wp), PARAMETER :: a_ddh = 2.5_wp, a_ddh_2 = 3.5 ! also in pycnocline_depth 1448 1448 1449 1449 IF( ln_timing ) CALL timing_start('zdf_osm_cd') … … 1466 1466 ! OSBL is deepening, entrainment > restratification 1467 1467 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 1468 ! *** Used for shear Needs to be changed to work stabily 1469 ! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml 1470 ! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd ) 1471 ! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) ) 1472 ! za_1 = 1.0 / zgamma_b**2 - 0.017 1473 ! za_2 = 1.0 / zgamma_b**3 - 0.0025 1474 ! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl ) 1475 zpsi = 0._wp 1476 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1477 zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj) 1468 zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1469 zpsi = -1.0_wp * zalpha_pyc(ji,jj) * ( zwb0(ji,jj) - ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj) 1470 zpsi = zpsi - 4.0_wp * ( 1.0_wp + zdh(ji,jj) / zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) 1471 zpsi = zalpha_b * MAX( zpsi, 0.0_wp ) 1472 zdhdt(ji,jj) = -1.0_wp * ( zwb_ent(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15_wp ) ) & 1473 & + zpsi / MAX( zdb_bl(ji,jj), 1e-15 ) 1478 1474 IF ( j_ddh(ji,jj) == 1 ) THEN 1479 1475 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN … … 1483 1479 ENDIF 1484 1480 ! Relaxation to dh_ref = zari * hbl 1485 zddhdt (ji,jj)= -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)1486 1487 ELSE ! j_ddh == 01481 zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1482 1483 ELSE IF ( j_ddh(ji,jj) == 0 ) THEN 1488 1484 ! Growing shear layer 1489 zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1490 ENDIF ! j_ddh 1491 zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj) 1485 zddhdt = -1.0_wp * a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1486 zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1e-8_wp ) ) * zddhdt 1487 ELSE 1488 zddhdt = 0.0_wp 1489 ENDIF ! j_ddh 1490 zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0_wp * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) & 1491 & * zddhdt / zdb_bl(ji,jj) 1492 1492 ELSE ! zdb_bl >0 1493 1493 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) … … 1516 1516 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1517 1517 ELSE 1518 zpert = MAX( 2.0 *zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )1518 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1519 1519 ENDIF 1520 1520 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) … … 1563 1563 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj) 1564 1564 ELSE 1565 zpert = MAX( 2.0 *zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )1565 zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1566 1566 ENDIF 1567 1567 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) … … 1695 1695 INTEGER :: jj, ji 1696 1696 INTEGER :: inhml 1697 REAL(wp) :: zari, ztau, zdh_ref 1698 REAL, PARAMETER :: a_ddh _2 = 3.5 ! also in pycnocline_depth1697 REAL(wp) :: zari, ztau, zdh_ref, zddhdt 1698 REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth 1699 1699 1700 1700 IF( ln_timing ) CALL timing_start('zdf_osm_pt') … … 1703 1703 IF ( lshear(ji,jj) ) THEN 1704 1704 IF ( lconv(ji,jj) ) THEN 1705 IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN 1705 1706 IF ( j_ddh(ji,jj) == 0 ) THEN 1706 1707 ! ddhdt for pycnocline determined in osm_calculate_dhdt 1707 dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt 1708 ! zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj) 1709 ! zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt 1710 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer 1711 ! dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) ) 1712 ztau = MAX( zdb_bl(ji,jj) * ( 0.625_wp * hbl(ji,jj) ) / ( a_ddh * MAX( -1.0_wp * zwb_ent(ji,jj), 1e-12_wp ) ), 2.0_wp * rn_Dt ) 1713 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.625_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 1714 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8_wp * zhbl(ji,jj) 1708 1715 ELSE 1709 1716 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt … … 1717 1724 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj) 1718 1725 ENDIF 1719 1726 ELSE 1727 ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5_wp * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt ) 1728 dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) + 0.2_wp * zhbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) ) 1729 IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj) 1730 END IF 1720 1731 ELSE ! lconv 1721 1732 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL … … 1910 1921 1911 1922 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1912 SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )1923 SUBROUTINE zdf_osm_mle_parameters( pmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle ) 1913 1924 !!---------------------------------------------------------------------- 1914 1925 !! *** ROUTINE zdf_osm_mle_parameters *** … … 1921 1932 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1922 1933 1934 REAL(wp), DIMENSION(jpi,jpj) :: pmld ! == estimated FK BLD used for MLE horiz gradients == ! 1923 1935 INTEGER, DIMENSION(jpi,jpj) :: mld_prof 1924 1936 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle … … 1961 1973 hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj)) 1962 1974 IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) ) 1975 ! For now try just set hmle to zmld 1976 hmle(ji,jj) = pmld(ji,jj) 1963 1977 END_2D 1964 1978 … … 2081 2095 END SUBROUTINE zdf_osm_vertical_average 2082 2096 2083 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, kbld_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear,&2097 SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear, & 2084 2098 & pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla, & 2085 2099 & pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml, & … … 2097 2111 INTEGER, DIMENSION(:,:), INTENT(in ) :: kmld ! ML base layer 2098 2112 INTEGER, DIMENSION(:,:), INTENT(in ) :: kp_ext ! Offset for external level 2099 INTEGER, INTENT(in ) :: kbld_ext ! Offset for external level2100 2113 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldconv ! BL stability flags 2101 2114 LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldpyc ! Pycnocline flags … … 2200 2213 IF ( jk <= kbld(ji,jj) ) THEN 2201 2214 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 2202 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) * &2215 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) * & 2203 2216 & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj) 2204 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) * &2217 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) * & 2205 2218 & ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj) 2206 2219 END IF … … 2251 2264 ! ---------------------------------------------------------------------- 2252 2265 WHERE ( ldconv(A2D(0)) ) 2253 zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) / &2266 zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / & 2254 2267 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2255 zsc_ws_1(:,:) = pwbav(A2D(0)) * pws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) / &2268 zsc_ws_1(:,:) = pwbav(A2D(0)) * pws0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) * phml(A2D(0)) / & 2256 2269 & ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln ) 2257 2270 ELSEWHERE … … 2263 2276 IF ( jk <= kmld(ji,jj) ) THEN 2264 2277 zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj) 2265 ! Calculate turbulent lengthscale2266 zl_c = 0.9_wp * ( 1.0_wp - EXP( - 7.0_wp * ( zznd_ml -zznd_ml**3 / 3.0_wp ) ) ) * &2267 & ( 1.0_wp - EXP( -15.0_wp * ( 1. 1_wp - zznd_ml ) ) )2268 zl_l = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml -zznd_ml**3 / 3.0_wp ) ) ) * &2269 & ( 1.0_wp - EXP( - 5.0_wp * ( 1.0_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) )2278 ! Calculate turbulent time scale 2279 zl_c = 0.9_wp * ( 1.0_wp - EXP( -5.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) * & 2280 & ( 1.0_wp - EXP( -15.0_wp * ( 1.2_wp - zznd_ml ) ) ) 2281 zl_l = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) * & 2282 & ( 1.0_wp - EXP( -8.0_wp * ( 1.15_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) ) 2270 2283 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp ) 2271 2284 ! Non-gradient buoyancy terms 2272 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0. 5_wp * zsc_wth_1(ji,jj) * zl_eps * phml(ji,jj)/ ( 0.15_wp + zznd_ml )2273 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0. 5_wp * zsc_ws_1(ji,jj) * zl_eps * phml(ji,jj)/ ( 0.15_wp + zznd_ml )2285 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_wth_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) 2286 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.4_wp * zsc_ws_1(ji,jj) * zl_eps / ( 0.15_wp + zznd_ml ) 2274 2287 END IF 2275 2288 ELSE ! Stable conditions … … 2288 2301 zws_ent(ji,jj) = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird * & 2289 2302 & ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj) 2290 zbuoy_pyc_sc = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 2291 zdelta_pyc = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird / & 2292 & SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 2293 zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) * & 2294 & zdelta_pyc**2 / pdh(ji,jj) 2295 zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) * & 2296 & zdelta_pyc**2 / pdh(ji,jj) 2297 zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 2303 IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 2304 zbuoy_pyc_sc = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj) 2305 zdelta_pyc = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird / & 2306 & SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) ) 2307 zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) * & 2308 & zdelta_pyc**2 / pdh(ji,jj) 2309 zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) * & 2310 & zdelta_pyc**2 / pdh(ji,jj) 2311 zzeta_pyc(ji,jj) = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) ) 2312 END IF 2298 2313 END IF 2299 2314 END_2D … … 2307 2322 & 0.045_wp * ( ( zws_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * & 2308 2323 & MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp ) 2309 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & 2310 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2311 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2312 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp * zws_pyc_sc_1(ji,jj) * & 2313 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2314 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2324 IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN 2325 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp * zwt_pyc_sc_1(ji,jj) * & 2326 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2327 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2328 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp * zws_pyc_sc_1(ji,jj) * & 2329 & EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) * & 2330 & pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird 2331 END IF 2315 2332 END IF ! End of pycnocline 2316 2333 END_3D … … 2348 2365 ! 2349 2366 DO_2D( 0, 0, 0, 0 ) 2350 IF ( ld pyc(ji,jj) ) THEN2367 IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN 2351 2368 IF ( k_ddh(ji,jj) == 0 ) THEN 2352 2369 ! Place holding code. Parametrization needs checking for these conditions. 2353 zomega 2370 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 2354 2371 zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 2355 2372 zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) 2356 2373 ELSE 2357 zomega 2374 zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird 2358 2375 zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj) 2359 2376 zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj) … … 2367 2384 END_2D 2368 2385 DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld ) ! Need ztau_sc_u to be available. Change to array. 2369 IF ( ld pyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN2386 IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN 2370 2387 zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 2371 2388 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) * & … … 2375 2392 & ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) * & 2376 2393 & ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk) 2377 END IF ! ld pyc2394 END IF ! ldconv .AND. ldpyc 2378 2395 END_3D 2379 2396 ! … … 2608 2625 ! 2609 2626 DO_2D( 0, 0, 0, 0 ) 2610 ghamt(ji,jj,kbld(ji,jj) +kbld_ext) = 0.0_wp2611 ghams(ji,jj,kbld(ji,jj) +kbld_ext) = 0.0_wp2612 ghamu(ji,jj,kbld(ji,jj) +kbld_ext) = 0.0_wp2613 ghamv(ji,jj,kbld(ji,jj) +kbld_ext) = 0.0_wp2627 ghamt(ji,jj,kbld(ji,jj)) = 0.0_wp 2628 ghams(ji,jj,kbld(ji,jj)) = 0.0_wp 2629 ghamu(ji,jj,kbld(ji,jj)) = 0.0_wp 2630 ghamv(ji,jj,kbld(ji,jj)) = 0.0_wp 2614 2631 END_2D 2615 2632 !
Note: See TracChangeset
for help on using the changeset viewer.
