Changeset 12245
- Timestamp:
- 2019-12-13T14:31:46+01:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0_FKOSM
- Files:
-
- 10 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/C1D_PAPA/cpp_C1D_PAPA.fcm
r9799 r12245 1 bld::tool::fppkeys key_c1d key_ mpp_mpi key_iomput1 bld::tool::fppkeys key_c1d key_iomput -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/domain_def_nemo.xml
r9930 r12245 13 13 <zoom_domain ibegin="1" jbegin="1" ni="1" nj="1"/> 14 14 </domain> 15 <!-- Eq section --> 15 <domain id="danger_T" domain_ref="grid_T" > 16 <zoom_domain id="danger_T" ibegin="45" jbegin="196" ni="5" nj="5"/> 17 </domain> 18 <domain id="danger_U" domain_ref="grid_U" > 19 <zoom_domain id="danger_U" ibegin="45" jbegin="196" ni="5" nj="5"/> 20 </domain> 21 <domain id="danger_V" domain_ref="grid_V" > 22 <zoom_domain id="danger_V" ibegin="45" jbegin="196" ni="5" nj="5"/> 23 </domain> 24 <domain id="danger_W" domain_ref="grid_W" > 25 <zoom_domain id="danger_W" ibegin="45" jbegin="196" ni="5" nj="5"/> 26 </domain> 27 <!-- Eq section --> 16 28 <domain id="EqT" domain_ref="grid_T" > <zoom_domain id="EqT"/> </domain> 17 29 <!-- TAO : see example above --> -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml
r10823 r12245 197 197 198 198 <field_group id="OSMOSIS_T" grid_ref="grid_T_2D"> 199 <field id="hml" long_name="mixed layr depth" unit="m" /> 200 <field id="hbl" long_name="boundary layer depth" unit="m" /> 201 <field id="dh" long_name="Pycnocline thickness" unit=" m" /> 202 <field id="ibld" long_name="index of boundary layer depth" unit="#" /> 203 <field id="imld" long_name="index of mixed layer depth" unit="#" /> 204 <field id="zhbl" long_name="boundary layer depth -grid" unit="m" /> 205 <field id="zhml" long_name="mixed layer depth - grid" unit="m" /> 206 <field id="zdh" long_name="Pycnocline depth - grid" unit=" m" /> 207 <field id="zustke" long_name="magnitude of stokes drift at T-points" unit="m/s" /> 208 <field id="us_x" long_name="i component of active Stokes drift" unit="m/s" /> 209 <field id="us_y" long_name="j component of active Stokes drift" unit="m/s" /> 210 <field id="dstokes" long_name="stokes drift depth scale" unit="m" /> 199 211 <field id="zwth0" long_name="surface non-local temperature flux" unit="deg m/s" /> 200 212 <field id="zws0" long_name="surface non-local salinity flux" unit="psu m/s" /> 201 <field id="hbl" long_name="boundary layer depth" unit="m" />202 <field id="hbli" long_name="initial boundary layer depth" unit="m" />203 <field id="dstokes" long_name="stokes drift depth scale" unit="m" />204 <field id="zustke" long_name="magnitude of stokes drift at T-points" unit="m/s" />205 213 <field id="zwstrc" long_name="convective velocity scale" unit="m/s" /> 214 <field id="zustar" long_name="friction velocity" unit="m/s" /> 206 215 <field id="zwstrl" long_name="langmuir velocity scale" unit="m/s" /> 207 <field id="zustar" long_name="friction velocity" unit="m/s" /> 208 <field id="zhbl" long_name="boundary layer depth" unit="m" /> 209 <field id="zhml" long_name="mixed layer depth" unit="m" /> 216 <field id="zvstr" long_name="mixed velocity scale" unit="m/s" /> 217 <field id="zla" long_name="langmuir number" unit="m/s" /> 210 218 <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2" unit="mW" /> 211 219 <field id="wind_wave_power" long_name="U_s \dot tau" unit="mW" /> 212 220 <field id="wind_power" long_name="\rho u*^3" unit="mW" /> 213 221 214 <!-- extraOSMOSIS diagnostics -->222 <!-- interior BL OSMOSIS diagnostics --> 215 223 <field id="zwthav" long_name="av turb flux of T in ml" unit="deg m/s" /> 216 224 <field id="zt_ml" long_name="av T in ml" unit="deg" /> 225 <field id="zhol" long_name="Hoenekker number" unit="#" /> 226 <field id="zws_ent" long_name="entrainment turb flux of S" unit="10^-3 m/s" /> 217 227 <field id="zwth_ent" long_name="entrainment turb flux of T" unit="deg m/s" /> 218 <field id="zhol" long_name="Hoenekker number" unit="#" /> 219 <field id="zdh" long_name="Pycnocline depth - grid" unit=" m" /> 220 </field_group> 221 222 <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" operation="instant" > 228 <field id="zwb_ent" long_name="entrainment turb flux of buoyancy" unit="m^2/s^-3" /> 229 230 <field id="zdt_bl" long_name="temperature jump at base of BL" unit="deg" /> 231 <field id="zds_bl" long_name="salinity jump at base of BL" unit="10^-3" /> 232 <field id="zdb_bl" long_name="buoyancy jump at base of BL" unit="m/s^2" /> 233 <field id="zdu_bl" long_name="u jump at base of BL" unit="m/s" /> 234 <field id="zdv_bl" long_name="v jump at base of BL" unit="m/s" /> 235 236 <!-- extra OSMOSIS diagnostics for debugging --> 237 <field id="zsc_uw_1_0" long_name="zsc u-momentum flux on T after Stokes" unit="m^2/s^2" /> 238 <field id="zsc_uw_1_f" long_name="zsc u-momentum flux on T after Coriolis" unit="m^2/s^2" /> 239 <field id="zsc_vw_1_f" long_name="zsc v-momentum flux on T after Coriolis" unit="m^2/s^2" /> 240 <field id="zsc_uw_2_f" long_name="2nd zsc u-momentum flux on T after Coriolis" unit="m^2/s^2" /> 241 <field id="zsc_vw_2_f" long_name="2nd zsc v-momentum flux on T after Coriolis" unit="m^2/s^2" /> 242 <field id="zuw_bse" long_name="base u-flux T-points" unit="m^2/s^2" /> 243 <field id="zvw_bse" long_name="base v-flux T-points" unit="m^2/s^2" /> 244 245 <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 246 <field id="hmle" long_name="OBL FK-layer thickness" unit="m" /> 247 <field id="mld_prof" long_name="FK-layer depth index" unit="#" /> 248 <field id="zmld" long_name="target FK-layer thickness" unit="m" /> 249 <field id="zwb_fk" long_name="FK b-flux" unit="m^2 s^-3" /> 250 <field id="zwb_fk_b" long_name="layer averaged FK b-flux" unit="m^2 s^-3" /> 251 <field id="zdiff_mle" long_name="max FK diffusivity in MLE" unit=" 10^-4 m^2 s^-1" /> 252 <field id="zvel_mle" long_name="FK velocity scale in MLE" unit=" m s^-1" /> 253 </field_group> 254 255 <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" > 256 <field id="zviscos" long_name="BL viscosity" unit="m^2/s" /> 223 257 <field id="ghamt" long_name="non-local temperature flux" unit="deg m/s" /> 224 258 <field id="ghams" long_name="non-local salinity flux" unit="psu m/s" /> 225 259 <field id="zdtdz_pyc" long_name="Pycnocline temperature gradient" unit=" deg/m" /> 226 </field_group> 260 <field id="zdsdz_pyc" long_name="Pycnocline salinity gradient" unit=" 10^-3/m" /> 261 <field id="zdbdz_pyc" long_name="Pycnocline buoyancy gradient" unit=" s^-2" /> 262 <field id="zdudz_pyc" long_name="Pycnocline u gradient" unit=" s^-2" /> 263 <field id="zdvdz_pyc" long_name="Pycnocline v gradient" unit=" s^-2" /> 264 265 <!-- extra OSMOSIS diagnostics for debugging --> 266 <field id="ghamu_00" long_name="initial non-local u-momentum flux" unit="m^2/s^2" /> 267 <field id="ghamv_00" long_name="initial non-local v-momentum flux" unit="m^2/s^2" /> 268 <field id="ghamu_0" long_name="after dstokes non-local u-momentum flux" unit="m^2/s^2" /> 269 <field id="ghamu_f" long_name="after Coriolis non-local u-momentum flux" unit="m^2/s^2" /> 270 <field id="ghamv_f" long_name="after Coriolis non-local v-momentum flux" unit="m^2/s^2" /> 271 <field id="ghamu_b" long_name="after buoyancy added non-local u-momentum flux" unit="m^2/s^2" /> 272 <field id="ghamv_b" long_name="after buoyancy added non-local v-momentum flux" unit="m^2/s^2" /> 273 <field id="ghamu_1" long_name="after entrainment non-local u-momentum flux" unit="m^2/s^2" /> 274 <field id="ghamv_1" long_name="after entrainment non-local v-momentum flux" unit="m^2/s^2" /> 275 </field_group> 227 276 228 277 <field_group id="OSMOSIS_U" grid_ref="grid_U_2D" > 229 278 <field id="ghamu" long_name="non-local u-momentum flux" grid_ref="grid_U_3D" unit="m^2/s^2" /> 230 <field id="us_x" long_name="i component of Stokes drift" unit="m/s" /> 231 </field_group> 279 <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 280 <field id="zdtdx" long_name="FK T x-gradient" unit=" deg C m^-1" /> 281 <field id="zdsdx" long_name="FK S x-gradient" unit=" 10^-3 m^-1" /> 282 <field id="dbdx_mle" long_name="FK B x-gradient" unit=" s^-2" /> 283 </field_group> 232 284 233 285 <field_group id="OSMOSIS_V" grid_ref="grid_V_2D" > 234 286 <field id="ghamv" long_name="non-local v-momentum flux" grid_ref="grid_V_3D" unit="m^2/s^2" /> 235 <field id="us_y" long_name="j component of Stokes drift" unit="m/s" /> 287 <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 288 <field id="zdtdy" long_name="FK T y-gradient" unit=" deg C m^-1" /> 289 <field id="zdsdy" long_name="FK S y-gradient" unit=" 10^-3 m^-1" /> 290 <field id="dbdy_mle" long_name="FK B y-gradient" unit=" s^-2" /> 236 291 </field_group> 237 292 … … 679 734 <field id="strd_zdfp" long_name="salinity -trend: pure vert. diffusion" unit="1e-3/s" /> 680 735 681 <!-- --> 736 <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 737 <field id="ttrd_osm" long_name="temperature-trend: OSM-OSBL non-local forcing" unit="degC/s" /> 738 <field id="strd_osm" long_name="salinity -trend: OSM-OSBL non-local forcing" unit="1e-3/s" /> 739 740 741 <!-- --> 682 742 <field id="ttrd_dmp" long_name="temperature-trend: interior restoring" unit="degC/s" /> 683 743 <field id="strd_dmp" long_name="salinity -trend: interior restoring" unit="1e-3/s" /> … … 715 775 <field id="strd_zdfp_e3t" unit="1e-3/s * m" > strd_zdfp * e3t </field> 716 776 777 <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 778 <field id="ttrd_osm" long_name="temperature-trend: OSM-OSBL non-local forcing" unit="degC/s * m" > ttrd_osm * e3t </field> 779 <field id="strd_osm" long_name="salinity -trend: OSM-OSBL non-local forcing" unit="1e-3/s * m" > strd_osm * e3t </field> 780 717 781 <!-- --> 718 782 <field id="ttrd_dmp_e3t" unit="degC/s * m" > ttrd_dmp * e3t </field> -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/grid_def_nemo.xml
r10226 r12245 48 48 <axis id="depthw" /> 49 49 </grid> 50 50 <!-- --> 51 51 <grid id="grid_1point" > 52 52 <domain domain_ref="1point"/> … … 57 57 <axis id="nfloat" /> 58 58 </grid> 59 <grid id="danger_T_2D" > 60 <domain id="danger_T" /> 61 </grid> 62 <!-- --> 63 <grid id="danger_T_3D" > 64 <domain id="danger_T" /> 65 <axis id="deptht" /> 66 </grid> 67 <!-- --> 68 <grid id="danger_U_2D" > 69 <domain id="danger_U" /> 70 </grid> 71 <!-- --> 72 <grid id="danger_U_3D" > 73 <domain id="danger_U" /> 74 <axis id="depthu" /> 75 </grid> 76 <!-- --> 77 <grid id="danger_V_2D" > 78 <domain id="danger_V" /> 79 </grid> 80 <!-- --> 81 <grid id="danger_V_3D" > 82 <domain id="danger_V" /> 83 <axis id="depthv" /> 84 </grid> 85 <!-- --> 86 <grid id="danger_W_2D" > 87 <domain id="danger_W" /> 88 </grid> 89 <!-- --> 90 <grid id="danger_W_3D" > 91 <domain id="danger_W" /> 92 <axis id="depthw" /> 93 </grid> 94 <!-- --> 95 <grid id="danger_T_2D" > 96 <domain id="danger_T" /> 97 </grid> 98 <!-- --> 99 <grid id="danger_T_3D" > 100 <domain id="danger_T" /> 101 <axis axis_ref="deptht" /> 102 </grid> 103 <!-- --> 104 <grid id="danger_U_2D" > 105 <domain id="danger_U" /> 106 </grid> 107 <!-- --> 108 <grid id="danger_U_3D" > 109 <domain id="danger_U" /> 110 <axis axis_ref="depthu" /> 111 </grid> 112 <!-- --> 113 <grid id="danger_V_2D" > 114 <domain id="danger_V" /> 115 </grid> 116 <!-- --> 117 <grid id="danger_V_3D" > 118 <domain id="danger_V" /> 119 <axis axis_ref="depthv" /> 120 </grid> 121 <!-- --> 122 <grid id="danger_W_2D" > 123 <domain id="danger_W" /> 124 </grid> 125 <!-- --> 126 <grid id="danger_W_3D" > 127 <domain id="danger_W" /> 128 <axis axis_ref="depthw" /> 129 </grid> 130 </grid_definition> 59 131 60 </grid_definition>61 -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/namelist_ref
r10888 r12245 1070 1070 &namzdf_osm ! OSM vertical diffusion (ln_zdfosm =T) 1071 1071 !----------------------------------------------------------------------- 1072 ln_use_osm_la = .false. ! Use namelistrn_osm_la1072 ln_use_osm_la = .false. ! Use rn_osm_la 1073 1073 rn_osm_la = 0.3 ! Turbulent Langmuir number 1074 1074 rn_osm_dstokes = 5. ! Depth scale of Stokes drift (m) … … 1085 1085 ! ! = 1: Pierson Moskowitz wave spectrum 1086 1086 ! ! = 0: Constant La# = 0.3 1087 ln_osm_mle = .false. ! Use integrated FK-OSM model 1088 / 1089 !----------------------------------------------------------------------- 1090 &namosm_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 1091 !----------------------------------------------------------------------- 1092 rn_osm_mle_ce = 0.06 ! magnitude of the MLE (typical value: 0.06 to 0.08) 1093 nn_osm_mle = 0 ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 1094 rn_osm_mle_lf = 5.e+3 ! typical scale of mixed layer front (meters) (case rn_osm_mle=0) 1095 rn_osm_mle_time = 172800. ! time scale for mixing momentum across the mixed layer (seconds) (case rn_osm_mle=0) 1096 rn_osm_mle_lat = 20. ! reference latitude (degrees) of MLE coef. (case rn_mle=1) 1097 rn_osm_mle_rho_c = 0.01 ! delta rho criterion used to calculate MLD for FK 1098 rn_osm_mle_thresh = 0.0005 ! delta b criterion used for FK MLD criterion 1099 rn_osm_mle_tau = 172800. ! time scale for FK-OSM (seconds) (case rn_osm_mle=0) 1087 1100 / 1088 1101 !----------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/ref_cfgs.txt
r10888 r12245 9 9 ORCA2_ICE_PISCES OCE TOP ICE NST 10 10 SPITZ12 OCE ICE 11 eORCA1 OCE ICE -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRA/tramle.F90
r10888 r12245 21 21 USE lbclnk ! lateral boundary condition / mpp link 22 22 23 USE zdfosm, ONLY : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 24 23 25 IMPLICIT NONE 24 26 PRIVATE … … 56 58 CONTAINS 57 59 58 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 59 !!---------------------------------------------------------------------- 60 !! *** ROUTINE tra_mle_trp *** 61 !! 62 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 63 !! 64 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 65 !! transport are computed as follows : 66 !! zu_mle = dk[ zpsi_uw ] 67 !! zv_mle = dk[ zpsi_vw ] 68 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 69 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 70 !! and added to the input velocity : 71 !! p.n = p.n + z._mle 72 !! 73 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 74 !! CAUTION, the transport is not updated at the last line/raw 75 !! this may be a problem for some advection schemes 76 !! 77 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 78 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 79 !!---------------------------------------------------------------------- 80 INTEGER , INTENT(in ) :: kt ! ocean time-step index 81 INTEGER , INTENT(in ) :: kit000 ! first time step index 82 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 83 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 84 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 85 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 86 ! 87 INTEGER :: ji, jj, jk ! dummy loop indices 88 INTEGER :: ii, ij, ik, ikmax ! local integers 89 REAL(wp) :: zcuw, zmuw, zc ! local scalar 90 REAL(wp) :: zcvw, zmvw ! - - 91 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 92 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 94 !!---------------------------------------------------------------------- 95 ! 96 ! !== MLD used for MLE ==! 97 ! ! compute from the 10m density to deal with the diurnal cycle 98 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 99 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 100 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 101 DO jj = 1, jpj 102 DO ji = 1, jpi ! index of the w-level at the ML based 103 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 104 END DO 105 END DO 106 END DO 107 ENDIF 108 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 109 ! 110 ! 111 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 112 zbm (:,:) = 0._wp 113 zn2 (:,:) = 0._wp 114 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 115 DO jj = 1, jpj 116 DO ji = 1, jpi 117 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 118 zmld(ji,jj) = zmld(ji,jj) + zc 119 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 120 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 121 END DO 122 END DO 123 END DO 124 125 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 126 CASE ( 0 ) != min of the 2 neighbour MLDs 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 ! vector opt. 129 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 130 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 131 END DO 132 END DO 133 CASE ( 1 ) != average of the 2 neighbour MLDs 134 DO jj = 1, jpjm1 135 DO ji = 1, fs_jpim1 ! vector opt. 136 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 137 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 138 END DO 139 END DO 140 CASE ( 2 ) != max of the 2 neighbour MLDs 141 DO jj = 1, jpjm1 142 DO ji = 1, fs_jpim1 ! vector opt. 143 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 144 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 145 END DO 146 END DO 147 END SELECT 148 ! ! convert density into buoyancy 149 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 150 ! 151 ! 152 ! !== Magnitude of the MLE stream function ==! 153 ! 154 ! di[bm] Ds 155 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 156 ! e1u Lf fu and the e2u for the "transport" 157 ! (not *e3u as divided by e3u at the end) 158 ! 159 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 160 DO jj = 1, jpjm1 161 DO ji = 1, fs_jpim1 ! vector opt. 162 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 163 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 164 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 165 ! 166 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 167 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 168 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 169 END DO 170 END DO 171 ! 172 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 173 DO jj = 1, jpjm1 174 DO ji = 1, fs_jpim1 ! vector opt. 175 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 176 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 177 ! 178 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 179 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 180 END DO 181 END DO 182 ENDIF 183 ! 184 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 185 DO jj = 1, jpjm1 186 DO ji = 1, fs_jpim1 ! vector opt. 187 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 188 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 189 END DO 190 END DO 191 ENDIF 192 ! 193 ! !== structure function value at uw- and vw-points ==! 194 DO jj = 1, jpjm1 195 DO ji = 1, fs_jpim1 ! vector opt. 196 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 197 zhv(ji,jj) = 1._wp / zhv(ji,jj) 198 END DO 199 END DO 200 ! 201 zpsi_uw(:,:,:) = 0._wp 202 zpsi_vw(:,:,:) = 0._wp 203 ! 204 DO jk = 2, ikmax ! start from 2 : surface value = 0 205 DO jj = 1, jpjm1 206 DO ji = 1, fs_jpim1 ! vector opt. 207 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 208 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 209 zcuw = zcuw * zcuw 210 zcvw = zcvw * zcvw 211 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 212 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 213 ! 214 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 215 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ! 220 ! !== transport increased by the MLE induced transport ==! 221 DO jk = 1, ikmax 222 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 223 DO ji = 1, fs_jpim1 ! vector opt. 224 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 225 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 226 END DO 227 END DO 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 60 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 61 !!---------------------------------------------------------------------- 62 !! *** ROUTINE tra_mle_trp *** 63 !! 64 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 65 !! 66 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 67 !! transport are computed as follows : 68 !! zu_mle = dk[ zpsi_uw ] 69 !! zv_mle = dk[ zpsi_vw ] 70 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 71 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 72 !! and added to the input velocity : 73 !! p.n = p.n + z._mle 74 !! 75 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 76 !! CAUTION, the transport is not updated at the last line/raw 77 !! this may be a problem for some advection schemes 78 !! 79 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 80 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 81 !!---------------------------------------------------------------------- 82 INTEGER , INTENT(in ) :: kt ! ocean time-step index 83 INTEGER , INTENT(in ) :: kit000 ! first time step index 84 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 85 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 86 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 88 ! 89 INTEGER :: ji, jj, jk ! dummy loop indices 90 INTEGER :: ii, ij, ik, ikmax ! local integers 91 REAL(wp) :: zcuw, zmuw, zc ! local scalar 92 REAL(wp) :: zcvw, zmvw ! - - 93 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 94 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 95 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 96 !!---------------------------------------------------------------------- 97 ! 98 ! 99 IF(ln_osm_mle) THEN 100 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 101 ! 102 ! 103 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 104 CASE ( 0 ) != min of the 2 neighbour MLDs 105 DO jj = 1, jpjm1 106 DO ji = 1, fs_jpim1 ! vector opt. 107 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 108 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 109 END DO 110 END DO 111 CASE ( 1 ) != average of the 2 neighbour MLDs 112 DO jj = 1, jpjm1 113 DO ji = 1, fs_jpim1 ! vector opt. 114 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 115 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 116 END DO 117 END DO 118 CASE ( 2 ) != max of the 2 neighbour MLDs 119 DO jj = 1, jpjm1 120 DO ji = 1, fs_jpim1 ! vector opt. 121 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 122 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 123 END DO 124 END DO 125 END SELECT 126 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 ! vector opt. 129 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 130 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 131 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 132 ! 133 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 134 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 135 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 136 END DO 137 END DO 138 ! 139 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 140 DO jj = 1, jpjm1 141 DO ji = 1, fs_jpim1 ! vector opt. 142 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 143 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 144 ! 145 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 146 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 147 END DO 148 END DO 149 ENDIF 150 151 ELSE !do not use osn_mle 152 ! !== MLD used for MLE ==! 153 ! ! convert density into buoyancy 154 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 155 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 156 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 157 DO jj = 1, jpj 158 DO ji = 1, jpi ! index of the w-level at the ML based 159 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 160 END DO 161 END DO 162 END DO 163 ENDIF 164 165 ! 166 ! 167 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 168 zbm (:,:) = 0._wp 169 zn2 (:,:) = 0._wp 170 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 171 DO jj = 1, jpj 172 DO ji = 1, jpi 173 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 174 zmld(ji,jj) = zmld(ji,jj) + zc 175 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 176 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 177 END DO 178 END DO 179 END DO 180 181 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 182 CASE ( 0 ) != min of the 2 neighbour MLDs 183 DO jj = 1, jpjm1 184 DO ji = 1, fs_jpim1 ! vector opt. 185 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 186 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 187 END DO 188 END DO 189 CASE ( 1 ) != average of the 2 neighbour MLDs 190 DO jj = 1, jpjm1 191 DO ji = 1, fs_jpim1 ! vector opt. 192 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 193 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 194 END DO 195 END DO 196 CASE ( 2 ) != max of the 2 neighbour MLDs 197 DO jj = 1, jpjm1 198 DO ji = 1, fs_jpim1 ! vector opt. 199 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 200 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 201 END DO 202 END DO 203 END SELECT 204 ! ! convert density into buoyancy 205 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 206 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 207 ! 208 ! 209 ! !== Magnitude of the MLE stream function ==! 210 ! 211 ! di[bm] Ds 212 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 213 ! e1u Lf fu and the e2u for the "transport" 214 ! (not *e3u as divided by e3u at the end) 215 ! 216 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 217 DO jj = 1, jpjm1 218 DO ji = 1, fs_jpim1 ! vector opt. 219 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 220 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 221 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 222 ! 223 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 224 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 225 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 226 END DO 227 END DO 228 ! 229 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 230 DO jj = 1, jpjm1 231 DO ji = 1, fs_jpim1 ! vector opt. 232 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 233 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 234 ! 235 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 236 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 237 END DO 238 END DO 239 ENDIF 240 ! 241 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 242 DO jj = 1, jpjm1 243 DO ji = 1, fs_jpim1 ! vector opt. 244 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 245 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 246 END DO 247 END DO 248 ENDIF 249 ! 250 ENDIF ! end of lm_osm_mle loop 251 ! !== structure function value at uw- and vw-points ==! 252 DO jj = 1, jpjm1 253 DO ji = 1, fs_jpim1 ! vector opt. 254 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 255 zhv(ji,jj) = 1._wp / zhv(ji,jj) 256 END DO 257 END DO 258 ! 259 zpsi_uw(:,:,:) = 0._wp 260 zpsi_vw(:,:,:) = 0._wp 261 ! 262 DO jk = 2, ikmax ! start from 2 : surface value = 0 263 DO jj = 1, jpjm1 264 DO ji = 1, fs_jpim1 ! vector opt. 265 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 266 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 267 zcuw = zcuw * zcuw 268 zcvw = zcvw * zcvw 269 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 270 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 271 ! 272 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 273 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 274 END DO 275 END DO 276 END DO 277 ! 278 ! !== transport increased by the MLE induced transport ==! 279 DO jk = 1, ikmax 280 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 281 DO ji = 1, fs_jpim1 ! vector opt. 282 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 283 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 284 END DO 285 END DO 286 DO jj = 2, jpjm1 287 DO ji = fs_2, fs_jpim1 ! vector opt. 288 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 231 289 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 232 END DO 233 END DO 234 END DO 235 236 IF( cdtype == 'TRA') THEN !== outputs ==! 237 ! 238 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 239 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 240 ! 241 ! divide by cross distance to give streamfunction with dimensions m^2/s 242 DO jk = 1, ikmax+1 243 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 244 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 245 END DO 246 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 247 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 248 ENDIF 249 ! 250 END SUBROUTINE tra_mle_trp 290 END DO 291 END DO 292 END DO 293 294 IF( cdtype == 'TRA') THEN !== outputs ==! 295 ! 296 IF (ln_osm_mle) THEN 297 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 298 ELSE 299 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 300 END IF 301 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 302 ! 303 ! divide by cross distance to give streamfunction with dimensions m^2/s 304 DO jk = 1, ikmax+1 305 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 306 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 307 END DO 308 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 309 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 310 ENDIF 311 ! 312 END SUBROUTINE tra_mle_trp 251 313 252 314 … … 290 352 WRITE(numout,*) ' Density difference used to define ML for FK rn_rho_c_mle = ', rn_rho_c_mle 291 353 ENDIF 292 ! 293 IF(lwp) THEN 354 355 356 IF( ln_osm_mle .AND. ln_mle ) THEN 357 WRITE(numout,*) 'WARNING: You are running with both OSM-FK and default FK' 358 ! CALL ctl_stop('STOP in zdf_osm_init: Cannot run with both OSM-FK and default FK') 359 END IF 360 361 IF(lwp) THEN 294 362 WRITE(numout,*) 295 363 IF( ln_mle ) THEN -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRD/trd_oce.F90
r10888 r12245 33 33 # endif 34 34 ! !!!* Active tracers trends indexes 35 INTEGER, PUBLIC, PARAMETER :: jptot_tra = 2 0!: Total trend nb: change it when adding/removing one indice below35 INTEGER, PUBLIC, PARAMETER :: jptot_tra = 21 !: Total trend nb: change it when adding/removing one indice below 36 36 ! =============== ! 37 37 INTEGER, PUBLIC, PARAMETER :: jptra_xad = 1 !: x- horizontal advection … … 46 46 INTEGER, PUBLIC, PARAMETER :: jptra_bbc = 10 !: Bottom Boundary Condition (geoth. heating) 47 47 INTEGER, PUBLIC, PARAMETER :: jptra_bbl = 11 !: Bottom Boundary Layer (diffusive and/or advective) 48 INTEGER, PUBLIC, PARAMETER :: jptra_osm = 21 !: Non-local terms from OSMOSIS OBL model 48 49 INTEGER, PUBLIC, PARAMETER :: jptra_npc = 12 !: non-penetrative convection treatment 49 50 INTEGER, PUBLIC, PARAMETER :: jptra_dmp = 13 !: internal restoring (damping) -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRD/trdtra.F90
r10888 r12245 346 346 CASE( jptra_bbl ) ; CALL iom_put( "ttrd_bbl" , ptrdx ) ! bottom boundary layer 347 347 CALL iom_put( "strd_bbl" , ptrdy ) 348 CASE( jptra_osm ) ; CALL iom_put( "ttrd_osm" , ptrdx ) ! OSMOSIS non-local forcing 349 CALL iom_put( "strd_osm" , ptrdy ) 348 350 CASE( jptra_npc ) ; CALL iom_put( "ttrd_npc" , ptrdx ) ! static instability mixing 349 351 CALL iom_put( "strd_npc" , ptrdy ) -
NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90
r10888 r12245 25 25 !! (12) Replace zwstrl with zvstr in calculation of eddy viscosity. 26 26 !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information 27 !! (14) B ouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).27 !! (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence. 28 28 !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added. 29 29 !! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out) 30 30 !! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out) 31 !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made, 32 !! (a) Pycnocline temperature and salinity profies changed for unstable layers 33 !! (b) The stable OSBL depth parametrization changed. 34 !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code. 35 !! 23/05/19 (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1 31 36 !!---------------------------------------------------------------------- 32 37 … … 40 45 !! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD) 41 46 !! dyn_osm : compute and add to u & v trensd the non-local flux 47 !! 48 !! Subroutines in revised code. 42 49 !!---------------------------------------------------------------------- 43 50 USE oce ! ocean dynamics and active tracers … … 69 76 PUBLIC tra_osm ! routine called by step.F90 70 77 PUBLIC trc_osm ! routine called by trcstp.F90 71 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 72 81 73 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux … … 77 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt 78 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth 79 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 80 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. 81 97 82 98 ! !!** Namelist namzdf_osm ** 83 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 84 103 REAL(wp) :: rn_osm_la ! Turbulent Langmuir number 85 104 REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift … … 96 115 REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s) 97 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 132 98 133 ! !!! ** General constants ** 99 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number 134 REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero 135 REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw 100 136 REAL(wp) :: pthird = 1._wp/3._wp ! 1/3 101 137 REAL(wp) :: p2third = 2._wp/3._wp ! 2/3 … … 114 150 !! *** FUNCTION zdf_osm_alloc *** 115 151 !!---------------------------------------------------------------------- 116 ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), & 117 & hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , & 118 & 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 119 168 IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays') 120 169 CALL mpp_sum ( 'zdfosm', zdf_osm_alloc ) … … 161 210 !! 162 211 INTEGER :: ji, jj, jk ! dummy loop indices 212 213 INTEGER :: jl ! dummy loop indices 214 163 215 INTEGER :: ikbot, jkmax, jkm1, jkp2 ! 164 216 … … 166 218 REAL(wp) :: zbeta, zthermal ! 167 219 REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales 168 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 220 REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm ! 221 169 222 REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density 170 223 INTEGER :: jm ! dummy loop indices … … 191 244 REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average 192 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 193 253 REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift 194 254 REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number … … 196 256 REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress 197 257 REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer 198 LOGICAL, DIMENSION( :,:), ALLOCATABLE :: lconv! unstable/stable bl258 LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl 199 259 200 260 ! mixed-layer variables … … 208 268 REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid 209 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 210 274 REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid 211 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 212 281 REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer 213 282 REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer … … 238 307 ! Temporary variables 239 308 INTEGER :: inhml 240 INTEGER :: i_lconv_alloc241 309 REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines 242 310 REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables … … 248 316 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity 249 317 318 INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme 319 REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt 320 REAL(wp) :: zzeta_s = 0._wp 321 REAL(wp) :: zzeta_v = 0.46 322 REAL(wp) :: zabsstke 323 250 324 ! For debugging 251 325 INTEGER :: ikt 252 326 !!-------------------------------------------------------------------- 253 327 ! 254 ALLOCATE( lconv(jpi,jpj), STAT= i_lconv_alloc )255 IF( i_lconv_alloc /= 0 ) CALL ctl_warn('zdf_osm: failed to allocate lconv')256 257 328 ibld(:,:) = 0 ; imld(:,:) = 0 258 329 zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp … … 268 339 zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp 269 340 zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp 341 270 342 zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp 271 343 zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp … … 277 349 zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp 278 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 279 361 ! Flux-Gradient arrays. 280 362 zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp … … 287 369 ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp 288 370 371 zdhdt_2(:,:) = 0._wp 289 372 ! hbl = MAX(hbl,epsln) 290 373 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 350 433 ! Use wind speed wndm included in sbc_oce module 351 434 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 352 dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav435 dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 353 436 END DO 354 437 END DO … … 362 445 ! It could represent the effects of the spread of wave directions 363 446 ! around the mean wind. The effect of this adjustment needs to be tested. 364 z ustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), &365 & zustar(ji,jj) / ( 0.45 * 0.45 ))366 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(z ustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !447 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 448 zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8) 449 dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 367 450 END DO 368 451 END DO … … 375 458 ! Langmuir velocity scale (zwstrl), at T-point 376 459 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 377 ! Modify zwstrl to allow for small and large values of dstokes/hbl. 378 ! Intended as a possible test. Doesn't affect LES results for entrainment, 379 ! but hasn't been shown to be correct as dstokes/h becomes large or small. 380 zwstrl(ji,jj) = zwstrl(ji,jj) * & 381 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 382 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 383 ! define La this way so effects of Stokes penetration depth on velocity scale are included 384 zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 460 zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2) 461 IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj)) 385 462 ! Velocity scale that tends to zustar for large Langmuir numbers 386 463 zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + & … … 389 466 ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 390 467 ! Note zustke and zwstrl are not amended. 391 IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45392 468 ! 393 469 ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv … … 406 482 ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth 407 483 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 408 ! BL must be always 2levels deep.409 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:, 3) )410 ibld(:,:) = 3411 DO jk = 4, jpkm1484 ! BL must be always 4 levels deep. 485 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 486 ibld(:,:) = 4 487 DO jk = 5, jpkm1 412 488 DO jj = 2, jpjm1 413 489 DO ji = 2, jpim1 … … 419 495 END DO 420 496 421 DO jj = 2, jpjm1 ! Vertical slab497 DO jj = 2, jpjm1 422 498 DO ji = 2, jpim1 423 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 424 zbeta = rab_n(ji,jj,1,jp_sal) 425 zt = 0._wp 426 zs = 0._wp 427 zu = 0._wp 428 zv = 0._wp 429 ! average over depth of boundary layer 430 zthick=0._wp 431 DO jm = 2, ibld(ji,jj) 432 zthick=zthick+e3t_n(ji,jj,jm) 433 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 434 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 435 zu = zu + e3t_n(ji,jj,jm) & 436 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 437 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 438 zv = zv + e3t_n(ji,jj,jm) & 439 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 440 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 441 END DO 442 zt_bl(ji,jj) = zt / zthick 443 zs_bl(ji,jj) = zs / zthick 444 zu_bl(ji,jj) = zu / zthick 445 zv_bl(ji,jj) = zv / zthick 446 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 447 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 448 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 449 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 450 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 451 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 452 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 453 IF ( lconv(ji,jj) ) THEN ! Convective 454 zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 455 & + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 456 457 zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 458 & * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 459 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 460 ! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 461 ! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 462 463 ! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 464 ! & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 465 zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 466 ELSE ! Stable 467 zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 468 & + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 469 & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 470 & * zwstrl(ji,jj)**3 / hbli(ji,jj) 471 zzdhdt = zzdhdt + zwbav(ji,jj) 472 IF ( zzdhdt < 0._wp ) THEN 473 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 474 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 475 ELSE 476 zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 477 & + MAX( zdb_bl(ji,jj), 0.0 ) 478 ENDIF 479 zzdhdt = 2.0 * zzdhdt / zpert 480 ENDIF 481 zdhdt(ji,jj) = zzdhdt 482 END DO 499 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 500 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 501 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 502 END DO 483 503 END DO 484 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 ) 485 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 486 528 imld = ibld ! use imld to hold previous blayer index 487 ibld(:,:) = 3 488 489 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 490 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 491 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 529 ibld(:,:) = 4 492 530 493 531 DO jk = 4, jpkm1 … … 495 533 DO ji = 2, jpim1 496 534 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk)535 ibld(ji,jj) = jk 498 536 ENDIF 499 537 END DO … … 504 542 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 505 543 ! 506 DO jj = 2, jpjm1507 DO ji = 2, jpim1508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN544 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 545 ! Alan: do we need zb_ml? 546 CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl ) 509 547 ! 510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.511 548 ! 512 zhbl_s = hbl(ji,jj) 513 jm = imld(ji,jj) 514 zthermal = rab_n(ji,jj,1,jp_tem) 515 zbeta = rab_n(ji,jj,1,jp_sal) 516 IF ( lconv(ji,jj) ) THEN 517 !unstable 518 zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 519 & * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 520 521 DO jk = imld(ji,jj), ibld(ji,jj) 522 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 523 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 524 525 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 526 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 527 528 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 529 END DO 530 hbl(ji,jj) = zhbl_s 531 ibld(ji,jj) = jm 532 hbli(ji,jj) = hbl(ji,jj) 533 ELSE 534 ! stable 535 DO jk = imld(ji,jj), ibld(ji,jj) 536 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 537 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 538 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 539 540 zhbl_s = zhbl_s + ( & 541 & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) & 542 & * zwstrl(ji,jj)**3 / hbli(ji,jj) & 543 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 544 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 545 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wn here 546 547 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 548 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 549 END DO 550 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 551 ibld(ji,jj) = MAX(jm, 3 ) 552 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 553 ENDIF ! IF ( lconv ) 554 ELSE 555 ! change zero or one model level. 556 hbl(ji,jj) = zhbl_t(ji,jj) 557 IF ( lconv(ji,jj) ) THEN 558 hbli(ji,jj) = hbl(ji,jj) 559 ELSE 560 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 561 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 562 ENDIF 563 ENDIF 564 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 565 END DO 566 END DO 549 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 567 550 dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms. 568 569 ! Recalculate averages over boundary layer after depth updated 570 ! Consider later combining this into the loop above and looking for columns 571 ! where the index for base of the boundary layer have changed 572 DO jj = 2, jpjm1 ! Vertical slab 573 DO ji = 2, jpim1 574 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 575 zbeta = rab_n(ji,jj,1,jp_sal) 576 zt = 0._wp 577 zs = 0._wp 578 zu = 0._wp 579 zv = 0._wp 580 ! average over depth of boundary layer 581 zthick=0._wp 582 DO jm = 2, ibld(ji,jj) 583 zthick=zthick+e3t_n(ji,jj,jm) 584 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 585 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 586 zu = zu + e3t_n(ji,jj,jm) & 587 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 588 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 589 zv = zv + e3t_n(ji,jj,jm) & 590 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 591 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 592 END DO 593 zt_bl(ji,jj) = zt / zthick 594 zs_bl(ji,jj) = zs / zthick 595 zu_bl(ji,jj) = zu / zthick 596 zv_bl(ji,jj) = zv / zthick 597 zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 598 zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 599 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 600 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 601 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 602 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 603 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 604 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 605 IF ( lconv(ji,jj) ) THEN 606 IF ( zdb_bl(ji,jj) > 0._wp )THEN 607 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 608 zari = 4.5 * ( zvstr(ji,jj)**2 ) & 609 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 610 ELSE ! unstable 611 zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 612 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 613 ENDIF 614 IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline 615 zari = 0.2 616 zwb_ent(ji,jj) = 0._wp 617 ENDIF 618 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 619 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 620 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 621 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 622 ELSE ! IF (zdb_bl) 623 imld(ji,jj) = ibld(ji,jj) - 1 624 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 625 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 626 ENDIF 627 ELSE ! IF (lconv) 628 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 629 ! boundary layer deepening 630 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 631 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 632 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 633 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 634 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 635 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 636 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 637 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 638 ELSE 639 imld(ji,jj) = ibld(ji,jj) - 1 640 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 641 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 642 ENDIF ! IF (zdb_bl > 0.0) 643 ELSE ! IF(dhdt >= 0) 644 ! boundary layer collapsing. 645 imld(ji,jj) = ibld(ji,jj) 646 zhml(ji,jj) = zhbl(ji,jj) 647 zdh(ji,jj) = 0._wp 648 ENDIF ! IF (dhdt >= 0) 649 ENDIF ! IF (lconv) 650 END DO 651 END DO 652 653 ! Average over the depth of the mixed layer in the convective boundary layer 654 ! Also calculate entrainment fluxes for temperature and salinity 655 DO jj = 2, jpjm1 ! Vertical slab 656 DO ji = 2, jpim1 657 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 658 zbeta = rab_n(ji,jj,1,jp_sal) 659 IF ( lconv(ji,jj) ) THEN 660 zt = 0._wp 661 zs = 0._wp 662 zu = 0._wp 663 zv = 0._wp 664 ! average over depth of boundary layer 665 zthick=0._wp 666 DO jm = 2, imld(ji,jj) 667 zthick=zthick+e3t_n(ji,jj,jm) 668 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 669 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 670 zu = zu + e3t_n(ji,jj,jm) & 671 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 672 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 673 zv = zv + e3t_n(ji,jj,jm) & 674 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 675 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 676 END DO 677 zt_ml(ji,jj) = zt / zthick 678 zs_ml(ji,jj) = zs / zthick 679 zu_ml(ji,jj) = zu / zthick 680 zv_ml(ji,jj) = zv / zthick 681 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 682 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 683 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 684 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 685 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 686 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 687 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 688 ELSE 689 ! stable, if entraining calulate average below interface layer. 690 IF ( zdhdt(ji,jj) >= 0._wp ) THEN 691 zt = 0._wp 692 zs = 0._wp 693 zu = 0._wp 694 zv = 0._wp 695 ! average over depth of boundary layer 696 zthick=0._wp 697 DO jm = 2, imld(ji,jj) 698 zthick=zthick+e3t_n(ji,jj,jm) 699 zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 700 zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 701 zu = zu + e3t_n(ji,jj,jm) & 702 & * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 703 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 704 zv = zv + e3t_n(ji,jj,jm) & 705 & * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 706 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 707 END DO 708 zt_ml(ji,jj) = zt / zthick 709 zs_ml(ji,jj) = zs / zthick 710 zu_ml(ji,jj) = zu / zthick 711 zv_ml(ji,jj) = zv / zthick 712 zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 713 zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 714 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 715 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 716 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 717 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 718 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 719 ENDIF 720 ENDIF 721 END DO 722 END DO 723 ! 551 ! 552 ! Average over the depth of the mixed layer in the convective boundary layer 553 ! Alan: do we need zb_ml? 554 CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 724 555 ! rotate mean currents and changes onto wind align co-ordinates 725 556 ! 726 727 DO jj = 2, jpjm1 728 DO ji = 2, jpim1 729 ztemp = zu_ml(ji,jj) 730 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 731 zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 732 ztemp = zdu_ml(ji,jj) 733 zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 734 zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 735 ! 736 ztemp = zu_bl(ji,jj) 737 zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 738 zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 739 ztemp = zdu_bl(ji,jj) 740 zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 741 zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 742 END DO 743 END DO 744 557 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 558 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 745 559 zuw_bse = 0._wp 746 560 zvw_bse = 0._wp 561 zwth_ent = 0._wp 562 zws_ent = 0._wp 747 563 DO jj = 2, jpjm1 748 564 DO ji = 2, jpim1 749 750 IF ( lconv(ji,jj) ) THEN 751 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 752 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 753 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 565 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 566 IF ( lconv(ji,jj) ) THEN 567 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 568 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 569 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 570 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 571 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 572 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 573 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 574 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 575 ENDIF 576 ELSE 577 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 578 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 754 579 ENDIF 755 ELSE756 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )757 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )758 580 ENDIF 759 581 END DO … … 764 586 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 765 587 766 DO jj = 2, jpjm1 767 DO ji = 2, jpim1 768 ! 769 IF ( lconv (ji,jj) ) THEN 770 ! Unstable conditions 771 IF( zdb_bl(ji,jj) > 0._wp ) THEN 772 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 773 ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 774 zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 775 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 776 DO jk = 2 , ibld(ji,jj) 777 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 778 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 779 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 780 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 781 END DO 782 ENDIF 783 ELSE 784 ! stable conditions 785 ! if pycnocline profile only defined when depth steady of increasing. 786 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 787 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 788 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 789 ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 790 zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 791 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 792 DO jk = 2, ibld(ji,jj) 793 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 794 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 795 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 796 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 797 END DO 798 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 799 ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 800 zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 801 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 802 DO jk = 2, ibld(ji,jj) 803 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 804 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 805 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 806 zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 807 END DO 808 ENDIF ! IF (zhol >=0.5) 809 ENDIF ! IF (zdb_bl> 0.) 810 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 811 ENDIF ! IF (lconv) 812 ! 813 END DO 814 END DO 815 ! 816 DO jj = 2, jpjm1 817 DO ji = 2, jpim1 818 ! 819 IF ( lconv (ji,jj) ) THEN 820 ! Unstable conditions 821 zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 822 & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 823 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 824 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 825 DO jk = 2 , ibld(ji,jj)-1 826 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 827 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 828 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 829 END DO 830 ELSE 831 ! stable conditions 832 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 833 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 834 DO jk = 2, ibld(ji,jj) 835 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 836 IF ( znd < 1.0 ) THEN 837 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 838 ELSE 839 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 840 ENDIF 841 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 842 END DO 843 ENDIF 844 ! 845 END DO 846 END DO 588 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 589 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 590 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 847 591 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 848 592 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 849 593 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 850 594 851 ! WHERE ( lconv )852 ! zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird853 ! zvisml_sc = zdifml_sc854 ! zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml )855 ! zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml )856 ! zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third857 ! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml858 ! ELSEWHERE859 ! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )860 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )861 ! ENDWHERE862 595 DO jj = 2, jpjm1 863 596 DO ji = 2, jpim1 … … 868 601 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 869 602 zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 870 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj)603 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 871 604 ELSE 872 605 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 873 606 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 874 END IF875 END DO876 END DO607 END IF 608 END DO 609 END DO 877 610 ! 878 611 DO jj = 2, jpjm1 … … 889 622 ! pycnocline - if present linear profile 890 623 IF ( zdh(ji,jj) > 0._wp ) THEN 624 zgamma_b = 6.0 891 625 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 626 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 893 627 ! 894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 +zznd_pyc )628 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 895 629 ! 896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 +zznd_pyc )630 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 897 631 END DO 632 IF ( ibld_ext == 0 ) THEN 633 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 634 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 635 ELSE 636 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 637 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 638 ENDIF 898 639 ENDIF 899 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 900 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) 901 643 ! could be taken out, take account of entrainment represents as a diffusivity 902 644 ! should remove w from here, represents entrainment … … 908 650 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 909 651 END DO 652 653 IF ( ibld_ext == 0 ) THEN 654 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 655 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 656 ELSE 657 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 658 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 659 ENDIF 910 660 ENDIF ! end if ( lconv ) 911 !661 ! 912 662 END DO ! end of ji loop 913 663 END DO ! end of jj loop … … 952 702 END DO ! end of jj loop 953 703 954 955 704 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero) 956 705 WHERE ( lconv ) 957 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0))958 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln)959 zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )706 zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 ) 707 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 708 zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln ) 960 709 ELSEWHERE 961 710 zsc_uw_1 = zustar**2 962 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)711 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 963 712 ENDWHERE 964 713 IF(ln_dia_osm) THEN 714 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 715 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 716 END IF 965 717 DO jj = 2, jpjm1 966 718 DO ji = 2, jpim1 … … 1007 759 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1008 760 & * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 1009 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)761 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 1010 762 ! non-gradient buoyancy terms 1011 763 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) … … 1020 772 END DO ! ji loop 1021 773 END DO ! jj loop 1022 1023 774 1024 775 WHERE ( lconv ) … … 1051 802 END DO ! jj loop 1052 803 804 IF(ln_dia_osm) THEN 805 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 806 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 807 END IF 1053 808 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 1054 809 … … 1089 844 END DO ! jj loop 1090 845 1091 1092 846 WHERE ( lconv ) 1093 847 zsc_uw_1 = zustar**2 … … 1134 888 END DO 1135 889 END DO 890 891 IF(ln_dia_osm) THEN 892 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 893 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 894 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 895 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 896 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 897 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 898 END IF 1136 899 ! 1137 900 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. … … 1165 928 END DO 1166 929 930 IF(ln_dia_osm) THEN 931 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 932 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 933 END IF 1167 934 ! pynocline contributions 1168 935 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small … … 1170 937 DO jj = 2, jpjm1 1171 938 DO ji = 2, jpim1 1172 DO jk= 2, ibld(ji,jj) 1173 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1174 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1175 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 1176 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 1177 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 1178 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 1179 END DO 1180 END DO 939 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 940 DO jk= 2, ibld(ji,jj) 941 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 942 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 943 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 944 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 945 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 946 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 947 END DO 948 END IF 949 END DO 1181 950 END DO 1182 951 … … 1185 954 DO jj=2, jpjm1 1186 955 DO ji = 2, jpim1 1187 IF ( lconv(ji,jj)) THEN956 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 1188 957 DO jk = 1, imld(ji,jj) - 1 1189 958 znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 1190 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd1191 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd959 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 960 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1192 961 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1193 962 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd … … 1195 964 DO jk = imld(ji,jj), ibld(ji,jj) 1196 965 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 1197 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )1198 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )966 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 967 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1199 968 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1200 969 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1201 970 END DO 1202 971 ENDIF 1203 ghamt(ji,jj,ibld(ji,jj)) = 0._wp 1204 ghams(ji,jj,ibld(ji,jj)) = 0._wp 1205 ghamu(ji,jj,ibld(ji,jj)) = 0._wp 1206 ghamv(ji,jj,ibld(ji,jj)) = 0._wp 972 973 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 974 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 975 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 976 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1207 977 END DO ! ji loop 1208 978 END DO ! jj loop 1209 979 1210 980 IF(ln_dia_osm) THEN 981 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 982 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 983 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 984 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 985 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 986 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 987 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 988 END IF 1211 989 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1212 990 ! Need to put in code for contributions that are applied explicitly to … … 1229 1007 END DO 1230 1008 END DO 1231 1232 IF(ln_dia_osm) THEN1233 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )1234 END IF1235 1009 1236 1010 ! KPP-style Ri# mixing … … 1286 1060 END IF ! ln_convmix = .true. 1287 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 1110 1288 1111 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1289 CALL lbc_lnk( 'zdfosm',zviscos(:,:,:), 'W', 1. )1112 !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 1290 1113 1291 1114 ! GN 25/8: need to change tmask --> wmask … … 1319 1142 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1320 1143 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & 1321 & ghamu, 'U', 1. , ghamv, 'V',1. )1322 1323 1144 & ghamu, 'U', -1. , ghamv, 'V', -1. ) 1145 1146 IF(ln_dia_osm) THEN 1324 1147 SELECT CASE (nn_osm_wave) 1325 1148 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1330 1153 ! Stokes drift read in from sbcwave (=2). 1331 1154 CASE(2) 1332 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift 1333 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift 1155 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1156 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1157 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1158 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1159 IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum 1160 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1161 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1334 1162 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 1335 1163 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1342 1170 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1343 1171 IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth 1344 IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth 1172 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1173 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1174 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1175 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1176 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1177 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1178 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1179 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1345 1180 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1346 1181 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1348 1183 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1349 1184 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1185 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1186 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1350 1187 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1351 1188 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 1352 1189 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1353 1190 IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine 1354 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine 1191 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1192 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1355 1193 IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine 1356 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine 1357 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine 1358 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 1359 1214 END IF 1360 ! Lateral boundary conditions on p_avt (sign unchanged) 1361 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 1362 1430 ! 1363 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 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1922 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1923 ELSE 1924 hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD 1925 ENDIF 1926 ENDIF 1927 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 1928 END DO 1929 END DO 1930 1931 mld_prof = 4 1932 DO jk = 5, jpkm1 1933 DO jj = 2, jpjm1 1934 DO ji = 2, jpim1 1935 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1936 END DO 1937 END DO 1938 END DO 1939 DO jj = 2, jpjm1 1940 DO ji = 2, jpim1 1941 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 1942 END DO 1943 END DO 1944 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 1945 1946 DO jj = 2, jpjm1 1947 DO ji = 2, jpim1 1948 IF ( lconv(ji,jj) ) THEN 1949 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1950 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1951 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1952 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1953 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1954 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1955 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1956 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1957 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1958 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1959 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1960 ENDIF 1961 END DO 1962 END DO 1963 END SUBROUTINE zdf_osm_mle_parameters 1964 1965 END SUBROUTINE zdf_osm 1364 1966 1365 1967 … … 1378 1980 INTEGER :: ios ! local integer 1379 1981 INTEGER :: ji, jj, jk ! dummy loop indices 1982 REAL z1_t2 1380 1983 !! 1381 1984 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1382 1985 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1383 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 1986 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 1987 & ,ln_osm_mle 1988 ! Namelist for Fox-Kemper parametrization. 1989 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 1990 & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 1991 1384 1992 !!---------------------------------------------------------------------- 1385 1993 ! … … 1397 2005 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1398 2006 WRITE(numout,*) '~~~~~~~~~~~~' 1399 WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters' 1400 WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2007 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 2008 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2009 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 1401 2010 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1402 2011 WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0 1403 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes2012 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1404 2013 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1405 2014 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1423 2032 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1424 2033 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl 2034 2035 IF( ln_osm_mle ) THEN 2036 ! Initialise Fox-Kemper parametrization 2037 REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme 2038 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2039 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2040 2041 REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 2042 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2043 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 2044 IF(lwm) WRITE ( numond, namosm_mle ) 2045 2046 IF(lwp) THEN ! Namelist print 2047 WRITE(numout,*) 2048 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 2049 WRITE(numout,*) '~~~~~~~~~~~~~' 2050 WRITE(numout,*) ' Namelist namosm_mle : ' 2051 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 2052 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 2053 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 2054 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 2055 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2056 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2057 WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2058 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2059 ENDIF ! 2060 ENDIF 2061 ! 2062 IF(lwp) THEN 2063 WRITE(numout,*) 2064 IF( ln_osm_mle ) THEN 2065 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 2066 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' 2067 IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' 2068 ELSE 2069 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 2070 ENDIF 2071 ENDIF 2072 ! 2073 IF( ln_osm_mle ) THEN ! MLE initialisation 2074 ! 2075 rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria 2076 IF(lwp) WRITE(numout,*) 2077 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 2078 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 2079 ! 2080 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 2081 ! 2082 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 2083 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 2084 ! 2085 ENDIF 2086 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2087 z1_t2 = 2.e-5 2088 do jj=1,jpj 2089 do ji = 1,jpi 2090 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2091 end do 2092 end do 2093 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 2094 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 2095 ! 2096 ENDIF 2097 2098 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 2099 1426 2100 1427 2101 IF( ln_zdfddm) THEN … … 1536 2210 REAL(wp) :: zN2_c ! local scalar 1537 2211 REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth 1538 INTEGER, DIMENSION( :,:), ALLOCATABLE:: imld_rst ! level of mixed-layer depth (pycnocline top)2212 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 1539 2213 !!---------------------------------------------------------------------- 1540 2214 ! … … 1551 2225 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1552 2226 END IF 2227 1553 2228 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1554 id2 = iom_varid( numror, ' hbli' , ldstop = .FALSE. )2229 id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) 1555 2230 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1556 2231 CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios ) 1557 CALL iom_get( numror, jpdom_autoglo, ' hbli', hbli, ldxios = lrxios )1558 WRITE(numout,*) ' ===>>>> : hbl & hbliread from restart file'2232 CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios ) 2233 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 1559 2234 RETURN 1560 ELSE ! 'hbl' & ' hbli' not in restart file, recalculate2235 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 1561 2236 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1562 2237 END IF … … 1570 2245 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn , ldxios = lwxios ) 1571 2246 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1572 CALL iom_rstput( kt, nitrst, numrow, ' hbli' , hbli, ldxios = lwxios )2247 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) 1573 2248 RETURN 1574 2249 END IF … … 1578 2253 !!----------------------------------------------------------------------------- 1579 2254 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1580 ALLOCATE( imld_rst(jpi,jpj) )1581 2255 ! w-level of the mixing and mixed layers 1582 2256 CALL eos_rab( tsn, rab_n ) … … 1599 2273 DO jj = 1, jpj 1600 2274 DO ji = 1, jpi 1601 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 2275 iiki = MAX(4,imld_rst(ji,jj)) 2276 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 2277 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 1603 2278 END DO 1604 2279 END DO 1605 hbl = MAX(hbl,epsln) 1606 hbli(:,:) = hbl(:,:)1607 DEALLOCATE( imld_rst ) 2280 2281 IF( ln_osm_mle ) hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2282 1608 2283 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 2284 wn(:,:,:) = 0._wp 2285 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1609 2286 END SUBROUTINE osm_rst 1610 2287 … … 1634 2311 ENDIF 1635 2312 1636 ! add non-local temperature and salinity flux1637 2313 DO jk = 1, jpkm1 1638 2314 DO jj = 2, jpjm1 … … 1648 2324 END DO 1649 2325 1650 1651 ! save the non-local tracer flux trends for diagnostic 2326 ! save the non-local tracer flux trends for diagnostics 1652 2327 IF( l_trdtra ) THEN 1653 2328 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 1654 2329 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 1655 !!bug gm jpttdzdf ==> jpttosm 1656 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ zdf, ztrdt )1657 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ zdf, ztrds )2330 2331 CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 2332 CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 1658 2333 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1659 2334 ENDIF … … 1723 2398 1724 2399 !!====================================================================== 2400 1725 2401 END MODULE zdfosm
Note: See TracChangeset
for help on using the changeset viewer.