Changeset 12218 for NEMO/branches/2019
- Timestamp:
- 2019-12-12T17:01:31+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90
r12216 r12218 76 76 PUBLIC tra_osm ! routine called by step.F90 77 77 PUBLIC trc_osm ! routine called by trcstp.F90 78 PUBLIC dyn_osm ! routine called by 'step.F90' 78 PUBLIC dyn_osm ! routine called by step.F90 79 80 PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90 79 81 80 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux … … 84 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt 85 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbli !: intial boundary layer depth for stable blayer 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth 87 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift. 91 92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML 95 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML 96 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer. 88 97 89 98 ! !!** Namelist namzdf_osm ** 90 99 LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la 100 101 LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation 102 91 103 REAL(wp) :: rn_osm_la ! Turbulent Langmuir number 92 104 REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift … … 102 114 LOGICAL :: ln_convmix = .true. ! Convective instability mixing 103 115 REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s) 116 117 ! OSMOSIS mixed layer eddy parametrization constants 118 INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt 119 REAL(wp) :: rn_osm_mle_ce ! MLE coefficient 120 ! ! parameters used in nn_osm_mle = 0 case 121 REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front 122 REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer 123 ! ! parameters used in nn_osm_mle = 1 case 124 REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front 125 REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK 126 REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation 127 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 128 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case 129 REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base. 130 REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE. 131 104 132 105 133 ! !!! ** General constants ** … … 122 150 !! *** FUNCTION zdf_osm_alloc *** 123 151 !!---------------------------------------------------------------------- 124 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 125 & hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , & 126 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 152 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & 153 & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 154 & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 155 156 ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), & 157 & mld_prof(jpi,jpj), STAT= zdf_osm_alloc ) 158 159 ! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ? 160 ! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), & 161 ! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc ) 162 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 163 ! IF ( ln_osm_mle ) THEN 164 ! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc ) 165 ! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays') 166 ! ENDIF 167 127 168 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 128 169 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) … … 169 210 !! 170 211 INTEGER :: ji, jj, jk ! dummy loop indices 212 213 INTEGER :: jl ! dummy loop indices 214 171 215 INTEGER :: ikbot, jkmax, jkm1, jkp2 ! 172 216 … … 200 244 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 201 245 REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux 246 247 248 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL 249 REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux 250 REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff 251 REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK 252 202 253 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift 203 254 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number … … 217 268 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 218 269 REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid 270 271 REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid 272 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid 273 219 274 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 220 275 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 276 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2 ! correction to dhdt due to internal structure. 277 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext ! external temperature/salinity and buoyancy gradients 278 279 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization. 280 221 281 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 222 282 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer … … 289 349 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 290 350 ! 351 zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp 352 353 IF ( ln_osm_mle ) THEN ! only initialise arrays if needed 354 zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp 355 zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp 356 zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp 357 zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp 358 ENDIF 359 zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 360 291 361 ! Flux-Gradient arrays. 292 362 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp … … 432 502 END DO 433 503 END DO 434 504 ! Averages over well-mixed and boundary layer 505 ! Alan: do we need zb_nl?, zb_ml? 506 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 507 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 508 ! External gradient 509 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 510 511 512 IF ( ln_osm_mle ) THEN 513 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 514 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 515 ENDIF 516 517 ! Rate of change of hbl 518 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 435 519 ! Calculate averages over depth of boundary layer 520 DO jj = 2, jpjm1 521 DO ji = 2, jpim1 522 zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it 523 ! adjustment to represent limiting by ocean bottom 524 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 525 END DO 526 END DO 527 436 528 imld = ibld ! use imld to hold previous blayer index 437 ibld(:,:) = 3 438 439 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 440 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 441 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 529 ibld(:,:) = 4 442 530 443 531 DO jk = 4, jpkm1 … … 550 638 ENDIF 551 639 ENDIF 552 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 553 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 640 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 641 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 642 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 554 643 ! could be taken out, take account of entrainment represents as a diffusivity 555 644 ! should remove w from here, represents entrainment … … 919 1008 END DO 920 1009 921 IF(ln_dia_osm) THEN922 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )923 END IF924 925 1010 ! KPP-style Ri# mixing 926 1011 IF( ln_kpprimix) THEN … … 974 1059 END DO 975 1060 END IF ! ln_convmix = .true. 1061 1062 1063 1064 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1065 DO jj = 2 , jpjm1 1066 DO ji = 2, jpim1 1067 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1068 ! Calculate MLE flux profiles 1069 ! DO jk = 1, mld_prof(ji,jj) 1070 ! znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1071 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1072 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1073 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1074 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1075 ! END DO 1076 ! Calculate MLE flux contribution from surface fluxes 1077 DO jk = 1, ibld(ji,jj) 1078 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 1079 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1080 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1081 END DO 1082 DO jk = 1, mld_prof(ji,jj) 1083 znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1084 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1085 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1086 END DO 1087 ! Viscosity for MLEs 1088 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1089 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1090 END DO 1091 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1092 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1093 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1094 DO jk = mld_prof(ji,jj), jl 1095 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1096 & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 1097 & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 1098 END DO 1099 ENDIF 1100 END DO 1101 END DO 1102 ENDIF 1103 1104 IF(ln_dia_osm) THEN 1105 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1106 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1107 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1108 END IF 1109 976 1110 977 1111 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids … … 1058 1192 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1059 1193 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1060 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine 1061 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine 1062 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1194 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1195 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1196 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1197 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux 1198 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1199 1200 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1201 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1202 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux 1203 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML 1204 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 1205 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1206 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1207 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1208 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1209 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1210 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1211 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1212 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1213 1063 1214 END IF 1064 ! Lateral boundary conditions on p_avt (sign unchanged) 1065 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 1215 1216 CONTAINS 1217 1218 1219 ! Alan: do we need zb? 1220 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 1221 !!--------------------------------------------------------------------- 1222 !! *** ROUTINE zdf_vertical_average *** 1223 !! 1224 !! ** Purpose : Determines vertical averages from surface to jnlev. 1225 !! 1226 !! ** Method : Averages are calculated from the surface to jnlev. 1227 !! The external level used to calculate differences is ibld+ibld_ext 1228 !! 1229 !!---------------------------------------------------------------------- 1230 1231 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1232 1233 ! Alan: do we need zb? 1234 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity 1235 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1236 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1237 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1238 1239 INTEGER :: jk, ji, jj 1240 REAL(wp) :: zthick, zthermal, zbeta 1241 1242 1243 zt = 0._wp 1244 zs = 0._wp 1245 zu = 0._wp 1246 zv = 0._wp 1247 DO jj = 2, jpjm1 ! Vertical slab 1248 DO ji = 2, jpim1 1249 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1250 zbeta = rab_n(ji,jj,1,jp_sal) 1251 ! average over depth of boundary layer 1252 zthick = epsln 1253 DO jk = 2, jnlev_av(ji,jj) 1254 zthick = zthick + e3t_n(ji,jj,jk) 1255 zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 1256 zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 1257 zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & 1258 & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 1259 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1260 zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & 1261 & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 1262 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1263 END DO 1264 zt(ji,jj) = zt(ji,jj) / zthick 1265 zs(ji,jj) = zs(ji,jj) / zthick 1266 zu(ji,jj) = zu(ji,jj) / zthick 1267 zv(ji,jj) = zv(ji,jj) / zthick 1268 ! Alan: do we need zb? 1269 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1270 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1271 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1272 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1273 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1274 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1275 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1276 END DO 1277 END DO 1278 END SUBROUTINE zdf_osm_vertical_average 1279 1280 SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 1281 !!--------------------------------------------------------------------- 1282 !! *** ROUTINE zdf_velocity_rotation *** 1283 !! 1284 !! ** Purpose : Rotates frame of reference of averaged velocity components. 1285 !! 1286 !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w 1287 !! 1288 !!---------------------------------------------------------------------- 1289 1290 REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle 1291 REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current 1292 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline 1293 1294 INTEGER :: ji, jj 1295 REAL(wp) :: ztemp 1296 1297 DO jj = 2, jpjm1 1298 DO ji = 2, jpim1 1299 ztemp = zu(ji,jj) 1300 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1301 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1302 ztemp = zdu(ji,jj) 1303 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1304 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1305 END DO 1306 END DO 1307 END SUBROUTINE zdf_osm_velocity_rotation 1308 1309 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1310 !!--------------------------------------------------------------------- 1311 !! *** ROUTINE zdf_osm_external_gradients *** 1312 !! 1313 !! ** Purpose : Calculates the gradients below the OSBL 1314 !! 1315 !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. 1316 !! 1317 !!---------------------------------------------------------------------- 1318 1319 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1320 1321 INTEGER :: jj, ji, jkb, jkb1 1322 REAL(wp) :: zthermal, zbeta 1323 1324 1325 DO jj = 2, jpjm1 1326 DO ji = 2, jpim1 1327 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1328 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1329 zbeta = rab_n(ji,jj,1,jp_sal) 1330 jkb = ibld(ji,jj) + ibld_ext 1331 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1332 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1333 & / e3t_n(ji,jj,ibld(ji,jj)) 1334 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1335 & / e3t_n(ji,jj,ibld(ji,jj)) 1336 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1337 END IF 1338 END DO 1339 END DO 1340 END SUBROUTINE zdf_osm_external_gradients 1341 1342 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 1343 1344 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1345 1346 INTEGER :: jk, jj, ji 1347 REAL(wp) :: ztgrad, zsgrad, zbgrad 1348 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1349 REAL(wp) :: zzeta_s=0.3 1350 1351 DO jj = 2, jpjm1 1352 DO ji = 2, jpim1 1353 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1354 IF ( lconv(ji,jj) ) THEN ! convective conditions 1355 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1356 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1357 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1358 ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj) 1359 zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj) 1360 zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj) 1361 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1362 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1363 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1364 DO jk = 2, ibld(ji,jj)+ibld_ext 1365 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 1366 IF ( znd <= zzeta_s ) THEN 1367 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * & 1368 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1369 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * & 1370 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1371 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * & 1372 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1373 ELSE 1374 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1375 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1376 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1377 ENDIF 1378 END DO 1379 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1380 ELSE 1381 ! stable conditions 1382 ! if pycnocline profile only defined when depth steady of increasing. 1383 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1384 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1385 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1386 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 1387 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 1388 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 1389 DO jk = 2, ibld(ji,jj) 1390 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1391 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1392 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1393 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1394 END DO 1395 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1396 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 1397 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 1398 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 1399 DO jk = 2, ibld(ji,jj) 1400 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 1401 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1402 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1403 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1404 END DO 1405 ENDIF ! IF (zhol >=0.5) 1406 ENDIF ! IF (zdb_bl> 0.) 1407 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1408 ENDIF ! IF (lconv) 1409 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1410 END DO 1411 END DO 1412 1413 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1414 1415 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1416 !!--------------------------------------------------------------------- 1417 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1418 !! 1419 !! ** Purpose : Calculates velocity shear in the pycnocline 1420 !! 1421 !! ** Method : 1422 !! 1423 !!---------------------------------------------------------------------- 1424 1425 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1426 1427 INTEGER :: jk, jj, ji 1428 REAL(wp) :: zugrad, zvgrad, znd 1429 REAL(wp) :: zzeta_v = 0.45 1066 1430 ! 1067 END SUBROUTINE zdf_osm 1431 DO jj = 2, jpjm1 1432 DO ji = 2, jpim1 1433 ! 1434 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1435 IF ( lconv (ji,jj) ) THEN 1436 ! Unstable conditions 1437 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1438 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1439 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1440 !Alan is this right? 1441 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1442 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1443 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1444 & )/ (zdh(ji,jj) + epsln ) 1445 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1446 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1447 IF ( znd <= 0.0 ) THEN 1448 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1449 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1450 ELSE 1451 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1452 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1453 ENDIF 1454 END DO 1455 ELSE 1456 ! stable conditions 1457 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1458 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1459 DO jk = 2, ibld(ji,jj) 1460 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1461 IF ( znd < 1.0 ) THEN 1462 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1463 ELSE 1464 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1465 ENDIF 1466 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1467 END DO 1468 ENDIF 1469 ! 1470 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1471 END DO 1472 END DO 1473 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1474 1475 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 1476 !!--------------------------------------------------------------------- 1477 !! *** ROUTINE zdf_osm_calculate_dhdt *** 1478 !! 1479 !! ** Purpose : Calculates the rate at which hbl changes. 1480 !! 1481 !! ** Method : 1482 !! 1483 !!---------------------------------------------------------------------- 1484 1485 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! Rate of change of hbl 1486 1487 INTEGER :: jj, ji 1488 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 1489 REAL(wp) :: zvel_max, zwb_min 1490 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1491 REAL(wp) :: zzeta_m = 0.3 1492 REAL(wp) :: zgamma_c = 2.0 1493 REAL(wp) :: zdhoh = 0.1 1494 REAL(wp) :: alpha_bc = 0.5 1495 1496 DO jj = 2, jpjm1 1497 DO ji = 2, jpim1 1498 IF ( lconv(ji,jj) ) THEN ! Convective 1499 ! Alan is this right? Yes, it's a bit different from the previous relationship 1500 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1501 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1502 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1503 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1504 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1505 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1506 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1507 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1508 1509 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1510 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1511 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1512 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1513 ! 1514 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1515 1516 IF ( ln_osm_mle ) THEN 1517 ! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + & 1518 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 1519 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1520 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1521 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1522 ELSE 1523 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1524 ENDIF 1525 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1526 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1527 ! OSBL is deepening, entrainment > restratification 1528 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1529 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1530 ELSE 1531 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 1532 ENDIF 1533 ELSE 1534 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 1535 zdhdt(ji,jj) = - zvel_mle(ji,jj) 1536 1537 1538 ENDIF 1539 1540 ELSE 1541 ! Fox-Kemper not used. 1542 1543 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1544 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1545 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1546 ! added ajgn 23 July as temporay fix 1547 1548 ENDIF 1549 1550 zdhdt_2(ji,jj) = 0._wp 1551 1552 ! commented out ajgn 23 July as temporay fix 1553 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1554 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 1555 ! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1556 ! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 1557 ! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1558 ! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 1559 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 1560 ! ! Alan no idea what this should be? 1561 ! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 1562 ! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 1563 ! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 1564 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 1565 ! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 1566 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1567 ! ENDIF 1568 ELSE ! Stable 1569 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1570 zdhdt_2(ji,jj) = 0._wp 1571 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1572 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1573 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1574 ELSE 1575 zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 1576 ENDIF 1577 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1578 ENDIF 1579 END DO 1580 END DO 1581 END SUBROUTINE zdf_osm_calculate_dhdt 1582 1583 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 1584 !!--------------------------------------------------------------------- 1585 !! *** ROUTINE zdf_osm_timestep_hbl *** 1586 !! 1587 !! ** Purpose : Increments hbl. 1588 !! 1589 !! ** Method : If thechange in hbl exceeds one model level the change is 1590 !! is calculated by moving down the grid, changing the buoyancy 1591 !! jump. This is to ensure that the change in hbl does not 1592 !! overshoot a stable layer. 1593 !! 1594 !!---------------------------------------------------------------------- 1595 1596 1597 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl. 1598 1599 INTEGER :: jk, jj, ji, jm 1600 REAL(wp) :: zhbl_s, zvel_max, zdb 1601 REAL(wp) :: zthermal, zbeta 1602 1603 DO jj = 2, jpjm1 1604 DO ji = 2, jpim1 1605 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1606 ! 1607 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 1608 ! 1609 zhbl_s = hbl(ji,jj) 1610 jm = imld(ji,jj) 1611 zthermal = rab_n(ji,jj,1,jp_tem) 1612 zbeta = rab_n(ji,jj,1,jp_sal) 1613 1614 1615 IF ( lconv(ji,jj) ) THEN 1616 !unstable 1617 1618 IF( ln_osm_mle ) THEN 1619 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1620 ELSE 1621 1622 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 1623 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1624 1625 ENDIF 1626 1627 DO jk = imld(ji,jj), ibld(ji,jj) 1628 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 1629 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 1630 & 0.0 ) + zvel_max 1631 1632 1633 IF ( ln_osm_mle ) THEN 1634 zhbl_s = zhbl_s + MIN( & 1635 & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1636 & e3w_n(ji,jj,jm) ) 1637 ELSE 1638 zhbl_s = zhbl_s + MIN( & 1639 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1640 & e3w_n(ji,jj,jm) ) 1641 ENDIF 1642 1643 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1644 1645 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1646 END DO 1647 hbl(ji,jj) = zhbl_s 1648 ibld(ji,jj) = jm 1649 ELSE 1650 ! stable 1651 DO jk = imld(ji,jj), ibld(ji,jj) 1652 zdb = MAX( & 1653 & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 1654 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 1655 & 0.0 ) + & 1656 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 1657 1658 ! Alan is thuis right? I have simply changed hbli to hbl 1659 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 1660 zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * & 1661 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 1662 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 1663 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1664 1665 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1666 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1667 END DO 1668 ENDIF ! IF ( lconv ) 1669 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 1670 ibld(ji,jj) = MAX(jm, 4 ) 1671 ELSE 1672 ! change zero or one model level. 1673 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 1674 ENDIF 1675 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 1676 END DO 1677 END DO 1678 1679 END SUBROUTINE zdf_osm_timestep_hbl 1680 1681 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 1682 !!--------------------------------------------------------------------- 1683 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 1684 !! 1685 !! ** Purpose : Calculates thickness of the pycnocline 1686 !! 1687 !! ** Method : The thickness is calculated from a prognostic equation 1688 !! that relaxes the pycnocine thickness to a diagnostic 1689 !! value. The time change is calculated assuming the 1690 !! thickness relaxes exponentially. This is done to deal 1691 !! with large timesteps. 1692 !! 1693 !!---------------------------------------------------------------------- 1694 1695 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1696 ! 1697 INTEGER :: jj, ji 1698 INTEGER :: inhml 1699 REAL(wp) :: zari, ztau, zddhdt 1700 1701 1702 DO jj = 2, jpjm1 1703 DO ji = 2, jpim1 1704 IF ( lconv(ji,jj) ) THEN 1705 1706 IF( ln_osm_mle ) THEN 1707 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1708 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1709 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1710 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1711 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1712 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1713 ELSE ! unstable 1714 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1715 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1716 ENDIF 1717 ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1718 zddhdt = zari * hbl(ji,jj) 1719 ELSE 1720 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1721 zddhdt = 0.2 * hbl(ji,jj) 1722 ENDIF 1723 ELSE 1724 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1725 zddhdt = 0.2 * hbl(ji,jj) 1726 ENDIF 1727 ELSE ! ln_osm_mle 1728 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1729 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1730 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1731 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1732 ELSE ! unstable 1733 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1734 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1735 ENDIF 1736 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1737 zddhdt = zari * hbl(ji,jj) 1738 ELSE 1739 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1740 zddhdt = 0.2 * hbl(ji,jj) 1741 ENDIF 1742 1743 END IF ! ln_osm_mle 1744 1745 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1746 ! Alan: this hml is never defined or used 1747 ELSE ! IF (lconv) 1748 ztau = hbl(ji,jj) / zvstr(ji,jj) 1749 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1750 ! boundary layer deepening 1751 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1752 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1753 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1754 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1755 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1756 ELSE 1757 zddhdt = 0.2 * hbl(ji,jj) 1758 ENDIF 1759 ELSE ! IF(dhdt < 0) 1760 zddhdt = 0.2 * hbl(ji,jj) 1761 ENDIF ! IF (dhdt >= 0) 1762 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1763 IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse 1764 ! Alan: this hml is never defined or used -- do we need it? 1765 ENDIF ! IF (lconv) 1766 1767 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1768 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1769 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1770 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1771 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1772 END DO 1773 END DO 1774 1775 END SUBROUTINE zdf_osm_pycnocline_thickness 1776 1777 1778 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 1779 !!---------------------------------------------------------------------- 1780 !! *** ROUTINE zdf_osm_horizontal_gradients *** 1781 !! 1782 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 1783 !! 1784 !! ** Method : 1785 !! 1786 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1787 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1788 1789 1790 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 1791 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 1792 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy 1793 1794 INTEGER :: ji, jj, jk ! dummy loop indices 1795 INTEGER :: ii, ij, ik, ikmax ! local integers 1796 REAL(wp) :: zc 1797 REAL(wp) :: zN2_c ! local buoyancy difference from 10m value 1798 REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH 1799 REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 1800 REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv 1801 !!---------------------------------------------------------------------- 1802 ! 1803 ! !== MLD used for MLE ==! 1804 1805 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 1806 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 1807 zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria 1808 DO jk = nlb10, jpkm1 1809 DO jj = 1, jpj ! Mixed layer level: w-level 1810 DO ji = 1, jpi 1811 ikt = mbkt(ji,jj) 1812 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 1813 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1814 END DO 1815 END DO 1816 END DO 1817 DO jj = 1, jpj 1818 DO ji = 1, jpi 1819 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 1820 zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1821 END DO 1822 END DO 1823 ! ensure mld_prof .ge. ibld 1824 ! 1825 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 1826 ! 1827 ztm(:,:) = 0._wp 1828 zsm(:,:) = 0._wp 1829 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 1830 DO jj = 1, jpj 1831 DO ji = 1, jpi 1832 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 1833 ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 1834 zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 1835 END DO 1836 END DO 1837 END DO 1838 ! average temperature and salinity. 1839 ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1840 zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1841 ! calculate horizontal gradients at u & v points 1842 1843 DO jj = 2, jpjm1 1844 DO ji = 1, jpim1 1845 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1846 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1847 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 1848 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 1849 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 1850 END DO 1851 END DO 1852 1853 DO jj = 1, jpjm1 1854 DO ji = 2, jpim1 1855 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1856 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1857 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 1858 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 1859 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 1860 END DO 1861 END DO 1862 1863 CALL eos_rab(ztsm_midu, zmld_midu, zabu) 1864 CALL eos_rab(ztsm_midv, zmld_midv, zabv) 1865 1866 DO jj = 2, jpjm1 1867 DO ji = 1, jpim1 1868 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 1869 END DO 1870 END DO 1871 DO jj = 1, jpjm1 1872 DO ji = 2, jpim1 1873 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 1874 END DO 1875 END DO 1876 1877 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1878 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 1879 !!---------------------------------------------------------------------- 1880 !! *** ROUTINE zdf_osm_mle_parameters *** 1881 !! 1882 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 1883 !! 1884 !! ** Method : 1885 !! 1886 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1887 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1888 1889 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 1890 INTEGER :: ji, jj, jk ! dummy loop indices 1891 INTEGER :: ii, ij, ik ! local integers 1892 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 1893 REAL(wp) :: zdb_mle, ztmp, zdbds_mle 1894 1895 mld_prof(:,:) = 4 1896 DO jk = 5, jpkm1 1897 DO jj = 2, jpjm1 1898 DO ji = 2, jpim1 1899 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1900 END DO 1901 END DO 1902 END DO 1903 ! DO jj = 2, jpjm1 1904 ! DO ji = 1, jpim1 1905 ! zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1906 ! END DO 1907 ! END DO 1908 ! Timestep mixed layer eddy depth. 1909 DO jj = 2, jpjm1 1910 DO ji = 2, jpim1 1911 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 1912 IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN 1913 hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy. 1914 ELSE 1915 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 1916 IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 1917 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1918 ELSE 1919 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1920 ENDIF 1921 ENDIF 1922 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 1923 END DO 1924 END DO 1925 1926 mld_prof = 4 1927 DO jk = 5, jpkm1 1928 DO jj = 2, jpjm1 1929 DO ji = 2, jpim1 1930 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1931 END DO 1932 END DO 1933 END DO 1934 DO jj = 2, jpjm1 1935 DO ji = 2, jpim1 1936 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 1937 END DO 1938 END DO 1939 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 1940 1941 DO jj = 2, jpjm1 1942 DO ji = 2, jpim1 1943 IF ( lconv(ji,jj) ) THEN 1944 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1945 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1946 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1947 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1948 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1949 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1950 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1951 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1952 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1953 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1954 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1955 ENDIF 1956 END DO 1957 END DO 1958 END SUBROUTINE zdf_osm_mle_parameters 1959 1960 END SUBROUTINE zdf_osm 1068 1961 1069 1962 … … 1086 1979 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1087 1980 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1088 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 1981 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 1982 & ,ln_osm_mle 1983 ! Namelist for Fox-Kemper parametrization. 1984 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 1985 & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 1986 1089 1987 !!---------------------------------------------------------------------- 1090 1988 ! … … 1102 2000 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1103 2001 WRITE(numout,*) '~~~~~~~~~~~~' 1104 WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters' 1105 WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2002 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 2003 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2004 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 1106 2005 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1107 2006 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 … … 1128 2027 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1129 2028 1130 call osm_rst( nit000, 'READ' ) !* read or initialize hbl 2029 2030 IF( ln_osm_mle ) THEN 2031 ! Initialise Fox-Kemper parametrization 2032 REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme 2033 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2034 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2035 2036 REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 2037 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2038 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 2039 IF(lwm) WRITE ( numond, namosm_mle ) 2040 2041 IF(lwp) THEN ! Namelist print 2042 WRITE(numout,*) 2043 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 2044 WRITE(numout,*) '~~~~~~~~~~~~~' 2045 WRITE(numout,*) ' Namelist namosm_mle : ' 2046 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 2047 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 2048 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 2049 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 2050 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2051 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2052 WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2053 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2054 ENDIF ! 2055 ENDIF 2056 ! 2057 IF(lwp) THEN 2058 WRITE(numout,*) 2059 IF( ln_osm_mle ) THEN 2060 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 2061 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' 2062 IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' 2063 ELSE 2064 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 2065 ENDIF 2066 ENDIF 2067 ! 2068 IF( ln_osm_mle ) THEN ! MLE initialisation 2069 ! 2070 rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria 2071 IF(lwp) WRITE(numout,*) 2072 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 2073 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 2074 ! 2075 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 2076 ! 2077 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 2078 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 2079 ! 2080 ENDIF 2081 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2082 z1_t2 = 2.e-5 2083 do jj=1,jpj 2084 do ji = 1,jpi 2085 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2086 end do 2087 end do 2088 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 2089 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 2090 ! 2091 ENDIF 2092 2093 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 2094 1131 2095 1132 2096 IF( ln_zdfddm) THEN … … 1309 2273 END DO 1310 2274 END DO 1311 hbl = MAX(hbl,epsln) 1312 hbli(:,:) = hbl(:,:)1313 DEALLOCATE( imld_rst ) 2275 2276 IF( ln_osm_mle ) hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2277 1314 2278 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 1315 2279 wn(:,:,:) = 0._wp
Note: See TracChangeset
for help on using the changeset viewer.
