Changeset 13402
- Timestamp:
- 2020-08-14T13:47:09+02:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715
- Files:
-
- 30 added
- 20 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/C1D_PAPA/cpp_C1D_PAPA.fcm
r9799 r13402 1 bld::tool::fppkeys key_c1d key_ mpp_mpi key_iomput1 bld::tool::fppkeys key_c1d key_iomput -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/GYRE_PISCES/EXPREF/namelist_cfg
r11536 r13402 142 142 rn_Ld = 100.e+3 ! lateral diffusive length [m] 143 143 / 144 !----------------------------------------------------------------------- 145 &namtra_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 146 !----------------------------------------------------------------------- 147 / 144 148 !!====================================================================== 145 149 !! *** Dynamics namelists *** !! … … 212 216 !----------------------------------------------------------------------- 213 217 nn_etau = 0 ! penetration of tke below the mixed layer (ML) due to internal & intertial waves 218 / 219 !----------------------------------------------------------------------- 220 &namzdf_osm ! OSM vertical diffusion (ln_zdfosm =T) 221 !----------------------------------------------------------------------- 222 / 223 !----------------------------------------------------------------------- 224 &namosm_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 225 !----------------------------------------------------------------------- 214 226 / 215 227 !!====================================================================== -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r11536 r13402 135 135 / 136 136 !----------------------------------------------------------------------- 137 &namtra_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 138 !----------------------------------------------------------------------- 139 / 140 !----------------------------------------------------------------------- 137 141 &namsbc_ssr ! surface boundary condition : sea surface restoring (ln_ssr =T) 138 142 !----------------------------------------------------------------------- … … 377 381 / 378 382 !----------------------------------------------------------------------- 383 &namzdf_osm ! OSM vertical diffusion (ln_zdfosm =T) 384 !----------------------------------------------------------------------- 385 / 386 !----------------------------------------------------------------------- 387 &namosm_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 388 !----------------------------------------------------------------------- 389 / 390 !----------------------------------------------------------------------- 379 391 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 380 392 !----------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/domain_def_nemo.xml
r11536 r13402 10 10 </domain> 11 11 12 <!-- Eq section --> 12 <domain id="1point" domain_ref="grid_T" > 13 <zoom_domain ibegin="1" jbegin="1" ni="1" nj="1"/> 14 </domain> 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 --> 13 28 <domain id="EqT" domain_ref="grid_T" > <zoom_domain id="EqT"/> </domain> 14 29 <!-- TAO : see example above --> -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/field_def_nemo-oce.xml
r11536 r13402 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 … … 245 300 <field id="emp_oce" long_name="Evap minus Precip over ocean" standard_name="evap_minus_precip_over_sea_water" unit="kg/m2/s" /> 246 301 <field id="emp_ice" long_name="Evap minus Precip over ice" standard_name="evap_minus_precip_over_sea_ice" unit="kg/m2/s" /> 247 <field id="saltflx" long_name=" Downward salt flux"unit="1e-3/m2/s" />248 <field id="fmmflx" long_name="Water flux due to freezing/melting" 302 <field id="saltflx" long_name="Total downward salinity flux" standard_name="total_downward_salinity_flux" unit="1e-3/m2/s" /> 303 <field id="fmmflx" long_name="Water flux due to freezing/melting" standard_name="water_flux_due_to_freezing/melting" unit="kg/m2/s" /> 249 304 <field id="snowpre" long_name="Snow precipitation" standard_name="snowfall_flux" unit="kg/m2/s" /> 250 305 <field id="runoffs" long_name="River Runoffs" standard_name="water_flux_into_sea_water_from_rivers" unit="kg/m2/s" /> … … 676 731 <field id="strd_zdfp" long_name="salinity -trend: pure vert. diffusion" unit="1e-3/s" /> 677 732 678 <!-- --> 733 <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 734 <field id="ttrd_osm" long_name="temperature-trend: OSM-OSBL non-local forcing" unit="degC/s" /> 735 <field id="strd_osm" long_name="salinity -trend: OSM-OSBL non-local forcing" unit="1e-3/s" /> 736 737 738 <!-- --> 679 739 <field id="ttrd_dmp" long_name="temperature-trend: interior restoring" unit="degC/s" /> 680 740 <field id="strd_dmp" long_name="salinity -trend: interior restoring" unit="1e-3/s" /> … … 712 772 <field id="strd_zdfp_e3t" unit="1e-3/s * m" > strd_zdfp * e3t </field> 713 773 774 <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 775 <field id="ttrd_osm_e3t" long_name="temperature-trend: OSM-OSBL non-local forcing" unit="degC/s * m" > ttrd_osm * e3t </field> 776 <field id="strd_osm_e3t" long_name="salinity -trend: OSM-OSBL non-local forcing" unit="1e-3/s * m" > strd_osm * e3t </field> 777 714 778 <!-- --> 715 779 <field id="ttrd_dmp_e3t" unit="degC/s * m" > ttrd_dmp * e3t </field> … … 720 784 <field id="strd_npc_e3t" unit="1e-3/s * m" > strd_npc * e3t </field> 721 785 <field id="ttrd_qns_e3t" unit="degC/s * m" > ttrd_qns * e3t_surf </field> 722 < field id="strd_cdt_e3t" unit="degC/s * m" > strd_cdt * e3t_surf </field>786 <!-- <field id="strd_cdt_e3t" unit="degC/s * m" > strd_cdt * e3t_surf </field> --> 723 787 <field id="ttrd_qsr_e3t" unit="degC/s * m" > ttrd_qsr * e3t </field> 724 788 <field id="ttrd_bbc_e3t" unit="degC/s * m" > ttrd_bbc * e3t </field> 725 789 726 790 <!-- OMIP layer-integrated trends --> 791 <field id="sfd" long_name="Total downward salt flux" unit="kg/(m^2 s)" > saltflux * 1026.0 * 0.001 </field> 792 <field id="wfo" long_name="Total downward FW mass flux" unit="kg/(m^2 s)" > -empmr </field> 793 727 794 <field id="ttrd_totad_li" long_name="layer integrated heat-trend: total advection" unit="W/m^2" > ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field> 728 795 <field id="strd_totad_li" long_name="layer integrated salt-trend: total advection" unit="kg/(m^2 s)" > strd_totad_e3t * 1026.0 * 0.001 </field> 796 <field id="ttrd_osm_li" long_name="layer integrated heat-trend: non-local OSM" unit="W/m^2" > ttrd_osm_e3t * 1026.0 * 3991.86795711963 </field> 797 <field id="strd_osm_li" long_name="layer integrated salt-trend: non-local OSM" unit="kg/(m^2 s)" > strd_osm_e3t * 1026.0 * 0.001 </field> 729 798 <field id="ttrd_evd_li" long_name="layer integrated heat-trend: EVD convection" unit="W/m^2" > ttrd_evd_e3t * 1026.0 * 3991.86795711963 </field> 730 799 <field id="strd_evd_li" long_name="layer integrated salt-trend: EVD convection" unit="kg/(m^2 s)" > strd_evd_e3t * 1026.0 * 0.001 </field> 800 <field id="ttrd_npc_li" long_name="layer integrated heat-trend: NPC convection" unit="W/m^2" > ttrd_npc_e3t * 1026.0 * 3991.86795711963 </field> 801 <field id="strd_npc_li" long_name="layer integrated salt-trend: NPC convection" unit="kg/(m^2 s)" > strd_npc_e3t * 1026.0 * 0.001 </field> 731 802 <field id="ttrd_iso_li" long_name="layer integrated heat-trend: isopycnal diffusion" unit="W/m^2" > ttrd_iso_e3t * 1026.0 * 3991.86795711963 </field> 732 803 <field id="strd_iso_li" long_name="layer integrated salt-trend: isopycnal diffusion" unit="kg/(m^2 s)" > strd_iso_e3t * 1026.0 * 0.001 </field> … … 735 806 <field id="ttrd_qns_li" long_name="layer integrated heat-trend: non-solar flux + runoff" unit="W/m^2" grid_ref="grid_T_2D"> ttrd_qns_e3t * 1026.0 * 3991.86795711963 </field> 736 807 <field id="ttrd_qsr_li" long_name="layer integrated heat-trend: solar flux" unit="W/m^2" grid_ref="grid_T_3D"> ttrd_qsr_e3t * 1026.0 * 3991.86795711963 </field> 808 <field id="ttrd_bbc_li" long_name="layer integrated heat-trend: geothermal heating " unit="W/m^2" > ttrd_bbc_e3t * 1026.0 * 3991.86795711963 </field> 737 809 <field id="ttrd_bbl_li" long_name="layer integrated heat-trend: bottom boundary layer " unit="W/m^2" > ttrd_bbl_e3t * 1026.0 * 3991.86795711963 </field> 738 810 <field id="strd_bbl_li" long_name="layer integrated salt-trend: bottom boundary layer " unit="kg/(m^2 s)" > strd_bbl_e3t * 1026.0 * 0.001 </field> -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/grid_def_nemo.xml
r11536 r13402 58 58 </grid> 59 59 60 </grid_definition> 60 <!-- --> 61 <grid id="grid_T_SFC"> 62 <domain domain_ref="grid_T" /> 63 <scalar> 64 <extract_axis position="0" /> 65 </scalar> 66 </grid> 67 </grid_definition> 61 68 -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/SHARED/namelist_ref
r11715 r13402 257 257 ln_COARE_3p5 = .false. ! "COARE 3.5" algorithm (Edson et al. 2013) 258 258 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 259 260 ln_humi_dpt = .false. ! Supply dewpoint tempearture instead of specific humidity (is true for ERA5) 261 259 262 ! 260 263 rn_zqt = 10. ! Air temperature & humidity reference height (m) … … 1086 1089 &namzdf_osm ! OSM vertical diffusion (ln_zdfosm =T) 1087 1090 !----------------------------------------------------------------------- 1088 ln_use_osm_la = .false. ! Use namelistrn_osm_la1091 ln_use_osm_la = .false. ! Use rn_osm_la 1089 1092 rn_osm_la = 0.3 ! Turbulent Langmuir number 1090 1093 rn_osm_dstokes = 5. ! Depth scale of Stokes drift (m) … … 1101 1104 ! ! = 1: Pierson Moskowitz wave spectrum 1102 1105 ! ! = 0: Constant La# = 0.3 1106 ln_osm_mle = .false. ! Use integrated FK-OSM model 1107 / 1108 !----------------------------------------------------------------------- 1109 &namosm_mle ! mixed layer eddy parametrisation (Fox-Kemper) (default: OFF) 1110 !----------------------------------------------------------------------- 1111 rn_osm_mle_ce = 0.06 ! magnitude of the MLE (typical value: 0.06 to 0.08) 1112 nn_osm_mle = 0 ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 1113 rn_osm_mle_lf = 5.e+3 ! typical scale of mixed layer front (meters) (case rn_osm_mle=0) 1114 rn_osm_mle_time = 172800. ! time scale for mixing momentum across the mixed layer (seconds) (case rn_osm_mle=0) 1115 rn_osm_mle_lat = 20. ! reference latitude (degrees) of MLE coef. (case rn_mle=1) 1116 rn_osm_mle_rho_c = 0.01 ! delta rho criterion used to calculate MLD for FK 1117 rn_osm_mle_thresh = 0.0005 ! delta b criterion used for FK MLD criterion 1118 rn_osm_mle_tau = 172800. ! time scale for FK-OSM (seconds) (case rn_osm_mle=0) 1103 1119 / 1104 1120 !----------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/cfgs/ref_cfgs.txt
r11715 r13402 9 9 ORCA2_ICE_PISCES OCE TOP ICE NST 10 10 SPITZ12 OCE ICE 11 eORCA1 OCE ICE -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/C1D/step_c1d.F90
r11715 r13402 58 58 59 59 indic = 0 ! reset to no error condition 60 IF( kstp == nit000 ) CALL iom_init( "nemo") ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 60 IF( kstp == nit000 ) THEN 61 CALL iom_init( "nemo") ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 62 CALL dia_hth_init ! extra ML depth diagnostics, thermocline depths 63 END IF 61 64 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) 62 65 CALL iom_setkt( kstp - nit000 + 1, "nemo" ) ! say to iom that we are at time step kstp … … 86 89 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 87 90 CALL dia_wri( kstp ) ! ocean model: outputs 88 IF( l k_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20°C)91 IF( ll_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20°C) 89 92 90 93 -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/DIA/diahth.F90
r11715 r13402 5 5 !!====================================================================== 6 6 !! History : OPA ! 1994-09 (J.-P. Boulanger) Original code 7 !! ! 1996-11 (E. Guilyardi) OPA8 7 !! ! 1996-11 (E. Guilyardi) OPA8 8 8 !! ! 1997-08 (G. Madec) optimization 9 !! ! 1999-07 (E. Guilyardi) hd28 + heat content 9 !! ! 1999-07 (E. Guilyardi) hd28 + heat content 10 10 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 11 11 !! 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag 12 !!---------------------------------------------------------------------- 13 #if defined key_diahth 14 !!---------------------------------------------------------------------- 15 !! 'key_diahth' : thermocline depth diag. 12 !! 4.0.1! 2020-01 (G. Nurser) remove need for key lk_diahth 16 13 !!---------------------------------------------------------------------- 17 14 !! dia_hth : Compute varius diagnostics associated with the mixed layer … … 24 21 USE lib_mpp ! MPP library 25 22 USE iom ! I/O library 26 USE timing ! preformance summary27 23 28 24 IMPLICIT NONE … … 30 26 31 27 PUBLIC dia_hth ! routine called by step.F90 32 PUBLIC dia_hth_alloc ! routine called by nemogcm.F90 33 34 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 35 36 ! note: following variables should move to local variables once iom_put is always used 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd20 !: depth of 20 C isotherm [m] 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd28 !: depth of 28 C isotherm [m] 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W] 28 PUBLIC dia_hth_init ! routine called by nemogcm.F90 29 30 LOGICAL, PUBLIC :: ll_diahth !: Compute further diagnostics of ML and thermocline depth 41 31 42 32 !!---------------------------------------------------------------------- … … 47 37 CONTAINS 48 38 49 FUNCTION dia_hth_alloc() 50 !!--------------------------------------------------------------------- 51 INTEGER :: dia_hth_alloc 52 !!--------------------------------------------------------------------- 39 40 SUBROUTINE dia_hth( kt ) 41 !!--------------------------------------------------------------------- 42 !! *** ROUTINE dia_hth *** 43 !! 44 !! ** Purpose : Computes 45 !! the mixing layer depth (turbocline): avt = 5.e-4 46 !! the depth of strongest vertical temperature gradient 47 !! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01) 48 !! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 49 !! the top of the thermochine: tn = tn(10m) - ztem2 50 !! the pycnocline depth with density criteria equivalent to a temperature variation 51 !! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 52 !! the barrier layer thickness 53 !! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 54 !! the depth of the 20 degree isotherm (linear interpolation) 55 !! the depth of the 28 degree isotherm (linear interpolation) 56 !! the heat content of first 300 m 57 !! 58 !! ** Method : 59 !!------------------------------------------------------------------- 60 INTEGER, INTENT( in ) :: kt ! ocean time-step index 61 !! 62 INTEGER :: ji, jj, jk ! dummy loop arguments 63 INTEGER :: iid, ilevel ! temporary integers 64 INTEGER, DIMENSION(jpi,jpj) :: ik20, ik28 ! levels 65 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 66 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 67 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 68 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 69 REAL(wp) :: zthick_0, zcoef ! temporary scalars 70 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 71 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 72 REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 73 REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 74 REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3 75 REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 76 REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion 77 REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion 78 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 79 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 80 REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz 81 REAL(wp), DIMENSION(jpi,jpj) :: zthick ! vertical integration thickness 82 REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2 83 ! note: following variables should move to local variables once iom_put is always used 84 REAL(wp), DIMENSION(jpi,jpj) :: zhth !: depth of the max vertical temperature gradient [m] 85 REAL(wp), DIMENSION(jpi,jpj) :: zhd20 !: depth of 20 C isotherm [m] 86 REAL(wp), DIMENSION(jpi,jpj) :: zhd28 !: depth of 28 C isotherm [m] 87 REAL(wp), DIMENSION(jpi,jpj) :: zhtc3 !: heat content of first 300 m [W] 88 89 IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN 90 ! ------------------------------------------------------------- ! 91 ! thermocline depth: strongest vertical gradient of temperature ! 92 ! turbocline depth (mixing layer depth): avt = zavt5 ! 93 ! MLD: rho = rho(1) + zrho3 ! 94 ! MLD: rho = rho(1) + zrho1 ! 95 ! ------------------------------------------------------------- ! 96 zmaxdzT(:,:) = 0._wp 97 IF( nla10 > 1 ) THEN 98 DO jj = 1, jpj 99 DO ji = 1, jpi 100 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 101 zrho0_3(ji,jj) = zztmp 102 zrho0_1(ji,jj) = zztmp 103 zhth(ji,jj) = zztmp 104 END DO 105 END DO 106 ELSE IF (iom_use("mlddzt")) THEN 107 DO jj = 1, jpj 108 DO ji = 1, jpi 109 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 110 zhth(ji,jj) = zztmp 111 END DO 112 END DO 113 ELSE 114 zhth(:,:) = 0._wp 115 116 ENDIF 117 118 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 119 DO jj = 1, jpj 120 DO ji = 1, jpi 121 ! 122 zzdep = gdepw_n(ji,jj,jk) 123 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 124 zzdep = zzdep * tmask(ji,jj,1) 125 126 IF( zztmp > zmaxdzT(ji,jj) ) THEN 127 zmaxdzT(ji,jj) = zztmp ; zhth (ji,jj) = zzdep ! max and depth of dT/dz 128 ENDIF 129 130 IF( nla10 > 1 ) THEN 131 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 132 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 133 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 134 ENDIF 135 136 END DO 137 END DO 138 END DO 139 140 IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) ) ! depth of the thermocline 141 IF( nla10 > 1 ) THEN 142 IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) ) ! MLD delta rho(surf) = 0.03 143 IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) ) ! MLD delta rho(surf) = 0.01 144 ENDIF 145 ENDIF 146 147 IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR. & 148 & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 152 zabs2 (ji,jj) = zztmp 153 ztm2 (ji,jj) = zztmp 154 zrho10_3(ji,jj) = zztmp 155 zpycn (ji,jj) = zztmp 156 END DO 157 END DO 158 ztinv (:,:) = 0._wp 159 zdepinv(:,:) = 0._wp 160 161 IF (iom_use("pycndep")) THEN 162 ! Preliminary computation 163 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 164 DO jj = 1, jpj 165 DO ji = 1, jpi 166 IF( tmask(ji,jj,nla10) == 1. ) THEN 167 zu = 1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80 * tsn(ji,jj,nla10,jp_sal) & 168 & - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) & 169 & - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 170 zv = 5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00 * tsn(ji,jj,nla10,jp_sal) & 171 & - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 172 zut = 11.25 - 0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01 * tsn(ji,jj,nla10,jp_sal) 173 zvt = 38.00 - 0.750 * tsn(ji,jj,nla10,jp_tem) 174 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 175 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 176 ELSE 177 zdelr(ji,jj) = 0._wp 178 ENDIF 179 END DO 180 END DO 181 ELSE 182 zdelr(:,:) = 0._wp 183 ENDIF 184 185 ! ------------------------------------------------------------- ! 186 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 187 ! Top of thermocline: tn = tn(10m) - ztem2 ! 188 ! MLD: rho = rho10m + zrho3 ! 189 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 190 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 191 ! depth of temperature inversion ! 192 ! ------------------------------------------------------------- ! 193 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 ! 197 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 198 ! 199 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) 200 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 201 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 202 zztmp = -zztmp ! delta T(10m) 203 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 204 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 205 ENDIF 206 207 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 208 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 209 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 210 ! 211 END DO 212 END DO 213 END DO 214 215 IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1) ) ! MLD abs(delta t) - 0.2 216 IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1) ) ! T(10) - 0.2 217 IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1) ) ! MLD delta rho(10m) = 0.03 218 IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1) ) ! MLD delta rho equi. delta T(10m) = 0.2 219 IF (iom_use("tinv")) CALL iom_put( "tinv" , ztinv*tmask(:,:,1) ) ! max. temp. inv. (t10 ref) 220 IF (iom_use("depti")) CALL iom_put( "depti" , zdepinv*tmask(:,:,1) ) ! depth of max. temp. inv. (t10 ref) 221 ENDIF 222 223 IF(iom_use("20d") .OR. iom_use("28d")) THEN 224 ! ----------------------------------- ! 225 ! search deepest level above 20C/28C ! 226 ! ----------------------------------- ! 227 ik20(:,:) = 1 228 ik28(:,:) = 1 229 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 zztmp = tsn(ji,jj,jk,jp_tem) 233 IF( zztmp >= 20. ) ik20(ji,jj) = jk 234 IF( zztmp >= 28. ) ik28(ji,jj) = jk 235 END DO 236 END DO 237 END DO 238 239 ! --------------------------- ! 240 ! Depth of 20C/28C isotherm ! 241 ! --------------------------- ! 242 DO jj = 1, jpj 243 DO ji = 1, jpi 244 ! 245 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 246 ! 247 iid = ik20(ji,jj) 248 IF( iid /= 1 ) THEN 249 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 250 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 251 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 252 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 253 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 254 ELSE 255 zhd20(ji,jj) = 0._wp 256 ENDIF 257 ! 258 iid = ik28(ji,jj) 259 IF( iid /= 1 ) THEN 260 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 261 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 262 & * ( 28.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 263 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 264 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth 265 ELSE 266 zhd28(ji,jj) = 0._wp 267 ENDIF 268 269 END DO 270 END DO 271 CALL iom_put( "20d", zhd20 ) ! depth of the 20 isotherm 272 CALL iom_put( "28d", zhd28 ) ! depth of the 28 isotherm 273 ENDIF 274 275 ! ----------------------------- ! 276 ! Heat content of first 300 m ! 277 ! ----------------------------- ! 278 IF (iom_use("hc300")) THEN 279 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 280 ilevel = 0 281 zthick_0 = 0._wp 282 DO jk = 1, jpkm1 283 zthick_0 = zthick_0 + e3t_1d(jk) 284 IF( zthick_0 < 300. ) ilevel = jk 285 END DO 286 ! surface boundary condition 287 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 288 ELSE ; zthick(:,:) = 0._wp ; zhtc3(:,:) = 0._wp 289 ENDIF 290 ! integration down to ilevel 291 DO jk = 1, ilevel 292 zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 293 zhtc3 (:,:) = zhtc3 (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 294 END DO 295 ! deepest layer 296 zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m 297 DO jj = 1, jpj 298 DO ji = 1, jpi 299 zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) & 300 & * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 301 END DO 302 END DO 303 ! from temperature to heat contain 304 zcoef = rau0 * rcp 305 zhtc3(:,:) = zcoef * zhtc3(:,:) 306 CALL iom_put( "hc300", zhtc3*tmask(:,:,1) ) ! first 300m heat content 307 ENDIF 308 ! 309 END SUBROUTINE dia_hth 310 311 312 SUBROUTINE dia_hth_init 313 !!--------------------------------------------------------------------------- 314 !! *** ROUTINE dia_hth_init *** 315 !! 316 !! ** Purpose: Initialization for ML and thermocline depths 317 !! 318 !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth 319 !!--------------------------------------------------------------------------- 320 ! 321 IF(lwp) THEN 322 WRITE(numout,*) 323 WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 324 WRITE(numout,*) '~~~~~~~~~~~~ ' 325 ENDIF 326 ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR. & 327 & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR. & 328 & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. & 329 & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300") 330 IF(lwp) THEN 331 WRITE(numout,*) ' output upper ocean diagnostics (T) or not (F) ll_diahth = ', ll_diahth 332 ENDIF 53 333 ! 54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 55 ! 56 CALL mpp_sum ( 'diahth', dia_hth_alloc ) 57 IF(dia_hth_alloc /= 0) CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' ) 58 ! 59 END FUNCTION dia_hth_alloc 60 61 62 SUBROUTINE dia_hth( kt ) 63 !!--------------------------------------------------------------------- 64 !! *** ROUTINE dia_hth *** 65 !! 66 !! ** Purpose : Computes 67 !! the mixing layer depth (turbocline): avt = 5.e-4 68 !! the depth of strongest vertical temperature gradient 69 !! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01) 70 !! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 71 !! the top of the thermochine: tn = tn(10m) - ztem2 72 !! the pycnocline depth with density criteria equivalent to a temperature variation 73 !! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 74 !! the barrier layer thickness 75 !! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 76 !! the depth of the 20 degree isotherm (linear interpolation) 77 !! the depth of the 28 degree isotherm (linear interpolation) 78 !! the heat content of first 300 m 79 !! 80 !! ** Method : 81 !!------------------------------------------------------------------- 82 INTEGER, INTENT( in ) :: kt ! ocean time-step index 83 !! 84 INTEGER :: ji, jj, jk ! dummy loop arguments 85 INTEGER :: iid, ilevel ! temporary integers 86 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ik20, ik28 ! levels 87 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 88 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 89 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 90 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 91 REAL(wp) :: zthick_0, zcoef ! temporary scalars 92 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 93 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho10_3 ! MLD: rho = rho10m + zrho3 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztinv ! max of temperature inversion 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdepinv ! depth of temperature inversion 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zmaxdzT ! max of dT/dz 103 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zthick ! vertical integration thickness 104 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdelr ! delta rho equivalent to deltaT = 0.2 105 !!---------------------------------------------------------------------- 106 IF( ln_timing ) CALL timing_start('dia_hth') 107 108 IF( kt == nit000 ) THEN 109 ! ! allocate dia_hth array 110 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 111 112 IF(.NOT. ALLOCATED(ik20) ) THEN 113 ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 114 & zabs2(jpi,jpj), & 115 & ztm2(jpi,jpj), & 116 & zrho10_3(jpi,jpj),& 117 & zpycn(jpi,jpj), & 118 & ztinv(jpi,jpj), & 119 & zdepinv(jpi,jpj), & 120 & zrho0_3(jpi,jpj), & 121 & zrho0_1(jpi,jpj), & 122 & zmaxdzT(jpi,jpj), & 123 & zthick(jpi,jpj), & 124 & zdelr(jpi,jpj), STAT=ji) 125 CALL mpp_sum('diahth', ji) 126 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 127 END IF 128 129 IF(lwp) WRITE(numout,*) 130 IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 131 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 132 IF(lwp) WRITE(numout,*) 133 ENDIF 134 135 ! initialization 136 ztinv (:,:) = 0._wp 137 zdepinv(:,:) = 0._wp 138 zmaxdzT(:,:) = 0._wp 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 142 hth (ji,jj) = zztmp 143 zabs2 (ji,jj) = zztmp 144 ztm2 (ji,jj) = zztmp 145 zrho10_3(ji,jj) = zztmp 146 zpycn (ji,jj) = zztmp 147 END DO 148 END DO 149 IF( nla10 > 1 ) THEN 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 153 zrho0_3(ji,jj) = zztmp 154 zrho0_1(ji,jj) = zztmp 155 END DO 156 END DO 157 ENDIF 158 159 ! Preliminary computation 160 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF( tmask(ji,jj,nla10) == 1. ) THEN 164 zu = 1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80 * tsn(ji,jj,nla10,jp_sal) & 165 & - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) & 166 & - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 167 zv = 5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00 * tsn(ji,jj,nla10,jp_sal) & 168 & - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 169 zut = 11.25 - 0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01 * tsn(ji,jj,nla10,jp_sal) 170 zvt = 38.00 - 0.750 * tsn(ji,jj,nla10,jp_tem) 171 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 172 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 173 ELSE 174 zdelr(ji,jj) = 0._wp 175 ENDIF 176 END DO 177 END DO 178 179 ! ------------------------------------------------------------- ! 180 ! thermocline depth: strongest vertical gradient of temperature ! 181 ! turbocline depth (mixing layer depth): avt = zavt5 ! 182 ! MLD: rho = rho(1) + zrho3 ! 183 ! MLD: rho = rho(1) + zrho1 ! 184 ! ------------------------------------------------------------- ! 185 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 186 DO jj = 1, jpj 187 DO ji = 1, jpi 188 ! 189 zzdep = gdepw_n(ji,jj,jk) 190 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 191 zzdep = zzdep * tmask(ji,jj,1) 192 193 IF( zztmp > zmaxdzT(ji,jj) ) THEN 194 zmaxdzT(ji,jj) = zztmp ; hth (ji,jj) = zzdep ! max and depth of dT/dz 195 ENDIF 196 197 IF( nla10 > 1 ) THEN 198 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 199 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 200 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 201 ENDIF 202 203 END DO 204 END DO 205 END DO 206 207 CALL iom_put( "mlddzt", hth ) ! depth of the thermocline 208 IF( nla10 > 1 ) THEN 209 CALL iom_put( "mldr0_3", zrho0_3 ) ! MLD delta rho(surf) = 0.03 210 CALL iom_put( "mldr0_1", zrho0_1 ) ! MLD delta rho(surf) = 0.01 211 ENDIF 212 213 ! ------------------------------------------------------------- ! 214 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 215 ! Top of thermocline: tn = tn(10m) - ztem2 ! 216 ! MLD: rho = rho10m + zrho3 ! 217 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 218 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 219 ! depth of temperature inversion ! 220 ! ------------------------------------------------------------- ! 221 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 222 DO jj = 1, jpj 223 DO ji = 1, jpi 224 ! 225 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 226 ! 227 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) 228 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 229 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 230 zztmp = -zztmp ! delta T(10m) 231 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 232 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 233 ENDIF 234 235 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 236 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 237 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 238 ! 239 END DO 240 END DO 241 END DO 242 243 CALL iom_put( "mld_dt02", zabs2 ) ! MLD abs(delta t) - 0.2 244 CALL iom_put( "topthdep", ztm2 ) ! T(10) - 0.2 245 CALL iom_put( "mldr10_3", zrho10_3 ) ! MLD delta rho(10m) = 0.03 246 CALL iom_put( "pycndep" , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2 247 CALL iom_put( "tinv" , ztinv ) ! max. temp. inv. (t10 ref) 248 CALL iom_put( "depti" , zdepinv ) ! depth of max. temp. inv. (t10 ref) 249 250 251 ! ----------------------------------- ! 252 ! search deepest level above 20C/28C ! 253 ! ----------------------------------- ! 254 ik20(:,:) = 1 255 ik28(:,:) = 1 256 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 257 DO jj = 1, jpj 258 DO ji = 1, jpi 259 zztmp = tsn(ji,jj,jk,jp_tem) 260 IF( zztmp >= 20. ) ik20(ji,jj) = jk 261 IF( zztmp >= 28. ) ik28(ji,jj) = jk 262 END DO 263 END DO 264 END DO 265 266 ! --------------------------- ! 267 ! Depth of 20C/28C isotherm ! 268 ! --------------------------- ! 269 DO jj = 1, jpj 270 DO ji = 1, jpi 271 ! 272 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 273 ! 274 iid = ik20(ji,jj) 275 IF( iid /= 1 ) THEN 276 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 277 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 278 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 279 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 280 hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 281 ELSE 282 hd20(ji,jj) = 0._wp 283 ENDIF 284 ! 285 iid = ik28(ji,jj) 286 IF( iid /= 1 ) THEN 287 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 288 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 289 & * ( 28.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 290 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 291 hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth 292 ELSE 293 hd28(ji,jj) = 0._wp 294 ENDIF 295 296 END DO 297 END DO 298 CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm 299 CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm 300 301 ! ----------------------------- ! 302 ! Heat content of first 300 m ! 303 ! ----------------------------- ! 304 305 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 306 ilevel = 0 307 zthick_0 = 0._wp 308 DO jk = 1, jpkm1 309 zthick_0 = zthick_0 + e3t_1d(jk) 310 IF( zthick_0 < 300. ) ilevel = jk 311 END DO 312 ! surface boundary condition 313 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 314 ELSE ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp 315 ENDIF 316 ! integration down to ilevel 317 DO jk = 1, ilevel 318 zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 319 htc3 (:,:) = htc3 (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 320 END DO 321 ! deepest layer 322 zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m 323 DO jj = 1, jpj 324 DO ji = 1, jpi 325 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) & 326 & * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 327 END DO 328 END DO 329 ! from temperature to heat contain 330 zcoef = rau0 * rcp 331 htc3(:,:) = zcoef * htc3(:,:) 332 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 333 ! 334 IF( ln_timing ) CALL timing_stop('dia_hth') 335 ! 336 END SUBROUTINE dia_hth 337 338 #else 339 !!---------------------------------------------------------------------- 340 !! Default option : Empty module 341 !!---------------------------------------------------------------------- 342 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .FALSE. !: thermocline-20d depths flag 343 CONTAINS 344 SUBROUTINE dia_hth( kt ) ! Empty routine 345 IMPLICIT NONE 346 INTEGER, INTENT( in ) :: kt 347 WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 348 END SUBROUTINE dia_hth 349 #endif 350 351 !!====================================================================== 334 END SUBROUTINE dia_hth_init 352 335 END MODULE diahth -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/DIA/diawri.F90
r11715 r13402 43 43 USE zdfdrg ! ocean vertical physics: top/bottom friction 44 44 USE zdfmxl ! mixed layer 45 USE zdfosm ! mixed layer 45 46 ! 46 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 926 927 CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity 927 928 ENDIF 929 930 IF( ln_zdfosm ) THEN 931 CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth 932 CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth 933 CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion 934 CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity 935 CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing 936 CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing 937 CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask ) ! non-local u forcing 938 CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask ) ! non-local v forcing 939 IF( ln_osm_mle ) THEN 940 CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth 941 END IF 942 ENDIF 928 943 929 944 #if defined key_si3 -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/IOM/iom.F90
r11715 r13402 856 856 ELSE ; llstop = .TRUE. 857 857 ENDIF 858 WRITE(numout,*)'iom_varid 1: entering routine for variable '//TRIM(cdvar) 859 WRITE(numout,*)'iom_varid 1: kiomid = ',kiomid 858 860 ! 859 861 IF( kiomid > 0 ) THEN … … 869 871 ll_fnd = ( TRIM(cdvar) == TRIM(iom_file(kiomid)%cn_var(iiv)) ) 870 872 END DO 873 874 WRITE(numout,*)'iom_varid 2: ll_fnd, iiv, nvars, llstop = ',ll_fnd, ' ', iiv, ' ', iom_file(kiomid)%nvars, ' ', llstop 875 871 876 ! 872 877 IF( .NOT.ll_fnd ) THEN … … 878 883 & 'increase the parameter jpmax_vars') 879 884 ENDIF 885 WRITE(numout,*)'iom_varid: cdvar = '//TRIM(cdvar) 886 IF( PRESENT(ldstop) ) WRITE(numout,*)'iom_varid: ldstop = ',ldstop 887 WRITE(numout,*)'iom_varid: llstop = ',llstop 888 WRITE(numout,*)'iom_varid: iom_varid = ',iom_varid 880 889 IF( llstop .AND. iom_varid == -1 ) CALL ctl_stop( TRIM(clinfo)//' not found' ) 881 890 ELSE -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/SBC/sbcblk.F90
r11715 r13402 124 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 125 125 126 LOGICAL :: ln_humi_dpt = .FALSE. ! calculate specific hunidity from dewpoint 127 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: qair ! specific humidity of air at input height 128 126 129 INTEGER :: nblk ! choice of the bulk algorithm 127 130 ! ! associated indices: … … 145 148 !!------------------------------------------------------------------- 146 149 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 147 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )150 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), qair(jpi,jpj), STAT=sbc_blk_alloc ) 148 151 ! 149 152 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) … … 171 174 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 172 175 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 173 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, 176 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, ln_humi_dpt,& ! bulk algorithm 174 177 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 178 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 … … 323 326 ! 324 327 ! ! compute the surface ocean fluxes using bulk formulea 328 ! ..... if dewpoint supplied instead of specific humidaity, calculate specific humidity 329 IF(ln_humi_dpt) THEN 330 qair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 331 ELSE 332 qair(:,:) = sf(jp_humi)%fnow(:,:,1) 333 END IF 334 325 335 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 326 336 … … 332 342 ENDIF 333 343 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 334 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)344 qatm_ice(:,:) = qair(:,:) 335 345 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 336 346 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 434 444 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 435 445 !! (since reanalysis products provide T at z, not theta !) 436 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt446 ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), qair(:,:) ) * rn_zqt 437 447 438 448 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 439 449 ! 440 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2450 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! NCAR-COREv2 441 451 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 442 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0452 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! COARE v3.0 443 453 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 444 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5454 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! COARE v3.5 445 455 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF456 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm, & ! ECMWF 447 457 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 458 CASE DEFAULT … … 454 464 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 455 465 ELSE ! At zt: 456 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )466 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 457 467 END IF 458 468 … … 495 505 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 496 506 !! q_air and t_air are given at 10m (wind reference height) 497 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed498 zqsb (:,:) = cp_air( sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed507 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - qair(:,:)) ) ! Evaporation, using bulk wind speed 508 zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 499 509 ELSE 500 510 !! q_air and t_air are not given at 10m (wind reference height) 501 511 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 502 512 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 503 zqsb (:,:) = cp_air( sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed513 zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 504 514 ENDIF 505 515 … … 742 752 ! local scalars ( place there for vector optimisation purposes) 743 753 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 744 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))754 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1)) 745 755 746 756 !!gm brutal.... … … 806 816 zcoef_dqla = -Ls * 11637800. * (-5897.8) 807 817 ! 808 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )818 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 809 819 ! 810 820 zztmp = 1. / ( 1. - albo ) … … 837 847 ! Latent Heat 838 848 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 839 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )849 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - qair(ji,jj) ) ) 840 850 ! Latent heat sensitivity for ice (Dqla/Dt) 841 851 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/TRA/tramle.F90
r11715 r13402 21 21 USE lbclnk ! lateral boundary condition / mpp link 22 22 23 ! where OSMOSIS_OBL is used with integrated FK 24 USE zdf_oce, ONLY : ln_zdfosm 25 USE zdfosm, ONLY : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 26 23 27 IMPLICIT NONE 24 28 PRIVATE … … 56 60 CONTAINS 57 61 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) & 62 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 63 !!---------------------------------------------------------------------- 64 !! *** ROUTINE tra_mle_trp *** 65 !! 66 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 67 !! 68 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 69 !! transport are computed as follows : 70 !! zu_mle = dk[ zpsi_uw ] 71 !! zv_mle = dk[ zpsi_vw ] 72 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 73 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 74 !! and added to the input velocity : 75 !! p.n = p.n + z._mle 76 !! 77 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 78 !! CAUTION, the transport is not updated at the last line/raw 79 !! this may be a problem for some advection schemes 80 !! 81 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 82 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 83 !!---------------------------------------------------------------------- 84 INTEGER , INTENT(in ) :: kt ! ocean time-step index 85 INTEGER , INTENT(in ) :: kit000 ! first time step index 86 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 88 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 89 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 90 ! 91 INTEGER :: ji, jj, jk ! dummy loop indices 92 INTEGER :: ii, ij, ik, ikmax ! local integers 93 REAL(wp) :: zcuw, zmuw, zc ! local scalar 94 REAL(wp) :: zcvw, zmvw ! - - 95 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 96 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 98 !!---------------------------------------------------------------------- 99 ! 100 ! 101 IF(ln_osm_mle.and.ln_zdfosm) THEN 102 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 103 ! 104 ! 105 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 106 CASE ( 0 ) != min of the 2 neighbour MLDs 107 DO jj = 1, jpjm1 108 DO ji = 1, fs_jpim1 ! vector opt. 109 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 110 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 111 END DO 112 END DO 113 CASE ( 1 ) != average of the 2 neighbour MLDs 114 DO jj = 1, jpjm1 115 DO ji = 1, fs_jpim1 ! vector opt. 116 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 117 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 118 END DO 119 END DO 120 CASE ( 2 ) != max of the 2 neighbour MLDs 121 DO jj = 1, jpjm1 122 DO ji = 1, fs_jpim1 ! vector opt. 123 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 124 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 125 END DO 126 END DO 127 END SELECT 128 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 129 DO jj = 1, jpjm1 130 DO ji = 1, fs_jpim1 ! vector opt. 131 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 132 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 133 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 134 ! 135 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 136 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 137 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 138 END DO 139 END DO 140 ! 141 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 142 DO jj = 1, jpjm1 143 DO ji = 1, fs_jpim1 ! vector opt. 144 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 145 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 146 ! 147 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 148 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 149 END DO 150 END DO 151 ENDIF 152 153 ELSE !do not use osn_mle 154 ! !== MLD used for MLE ==! 155 ! ! compute from the 10m density to deal with the diurnal cycle 156 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 157 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 158 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 159 DO jj = 1, jpj 160 DO ji = 1, jpi ! index of the w-level at the ML based 161 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 162 END DO 163 END DO 164 END DO 165 ENDIF 166 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 167 168 ! 169 ! 170 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 171 zbm (:,:) = 0._wp 172 zn2 (:,:) = 0._wp 173 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 174 DO jj = 1, jpj 175 DO ji = 1, jpi 176 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 177 zmld(ji,jj) = zmld(ji,jj) + zc 178 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 179 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 180 END DO 181 END DO 182 END DO 183 184 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 185 CASE ( 0 ) != min of the 2 neighbour MLDs 186 DO jj = 1, jpjm1 187 DO ji = 1, fs_jpim1 ! vector opt. 188 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 189 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 190 END DO 191 END DO 192 CASE ( 1 ) != average of the 2 neighbour MLDs 193 DO jj = 1, jpjm1 194 DO ji = 1, fs_jpim1 ! vector opt. 195 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 196 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 197 END DO 198 END DO 199 CASE ( 2 ) != max of the 2 neighbour MLDs 200 DO jj = 1, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 203 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 204 END DO 205 END DO 206 END SELECT 207 ! ! convert density into buoyancy 208 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 209 ! 210 ! 211 ! !== Magnitude of the MLE stream function ==! 212 ! 213 ! di[bm] Ds 214 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 215 ! e1u Lf fu and the e2u for the "transport" 216 ! (not *e3u as divided by e3u at the end) 217 ! 218 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 219 DO jj = 1, jpjm1 220 DO ji = 1, fs_jpim1 ! vector opt. 221 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 222 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 223 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 224 ! 225 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 226 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 227 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 228 END DO 229 END DO 230 ! 231 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 232 DO jj = 1, jpjm1 233 DO ji = 1, fs_jpim1 ! vector opt. 234 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 235 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 236 ! 237 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 238 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 239 END DO 240 END DO 241 ENDIF 242 ! 243 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 244 DO jj = 1, jpjm1 245 DO ji = 1, fs_jpim1 ! vector opt. 246 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 247 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 248 END DO 249 END DO 250 ENDIF 251 ! 252 ENDIF ! end of lm_osm_mle loop 253 ! !== structure function value at uw- and vw-points ==! 254 DO jj = 1, jpjm1 255 DO ji = 1, fs_jpim1 ! vector opt. 256 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 257 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 258 END DO 259 END DO 260 ! 261 zpsi_uw(:,:,:) = 0._wp 262 zpsi_vw(:,:,:) = 0._wp 263 ! 264 DO jk = 2, ikmax ! start from 2 : surface value = 0 265 DO jj = 1, jpjm1 266 DO ji = 1, fs_jpim1 ! vector opt. 267 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 268 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 269 zcuw = zcuw * zcuw 270 zcvw = zcvw * zcvw 271 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 272 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 273 ! 274 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 275 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 276 END DO 277 END DO 278 END DO 279 ! 280 ! !== transport increased by the MLE induced transport ==! 281 DO jk = 1, ikmax 282 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 283 DO ji = 1, fs_jpim1 ! vector opt. 284 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 285 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 286 END DO 287 END DO 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 231 291 & + 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 292 END DO 293 END DO 294 END DO 295 296 IF( cdtype == 'TRA') THEN !== outputs ==! 297 ! 298 IF (ln_osm_mle.and.ln_zdfosm) THEN 299 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 300 ELSE 301 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 302 END IF 303 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 304 ! 305 ! divide by cross distance to give streamfunction with dimensions m^2/s 306 DO jk = 1, ikmax+1 307 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 308 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 309 END DO 310 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 311 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 312 ENDIF 313 ! 314 END SUBROUTINE tra_mle_trp 251 315 252 316 -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/TRD/trd_oce.F90
r11715 r13402 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.1_FKOSM_m11715/src/OCE/TRD/trdtra.F90
r11715 r13402 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.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90
r11715 r13402 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 … … 105 141 !!---------------------------------------------------------------------- 106 142 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 107 !! $Id $143 !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $ 108 144 !! Software governed by the CeCILL license (see ./LICENSE) 109 145 !!---------------------------------------------------------------------- … … 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 2 levels deep. 409 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 410 ibld(:,:) = 3 411 DO jk = 4, jpkm1 412 DO jj = 2, jpjm1 413 DO ji = 2, jpim1 484 ! BL must be always 4 levels deep. 485 ! For calculation of lateral buoyancy gradients for FK in 486 ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 487 ! previously exist for hbl also. 488 hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 489 ibld(:,:) = 4 490 DO jk = 5, jpkm1 491 DO jj = 1, jpj 492 DO ji = 1, jpi 414 493 IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 415 494 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) … … 419 498 END DO 420 499 421 DO jj = 2, jpjm1 ! Vertical slab500 DO jj = 2, jpjm1 422 501 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 502 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 503 imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 )) 504 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 505 END DO 483 506 END DO 484 507 ! Averages over well-mixed and boundary layer 508 ! Alan: do we need zb_nl?, zb_ml? 509 CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl) 510 CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) 511 ! External gradient 512 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 513 514 515 IF ( ln_osm_mle ) THEN 516 CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 517 CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 518 ENDIF 519 520 ! Rate of change of hbl 521 CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 485 522 ! Calculate averages over depth of boundary layer 486 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 523 DO jj = 2, jpjm1 524 DO ji = 2, jpim1 525 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 526 ! adjustment to represent limiting by ocean bottom 527 zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:)) 528 END DO 529 END DO 530 531 imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index 532 ibld(:,:) = 4 492 533 493 534 DO jk = 4, jpkm1 … … 495 536 DO ji = 2, jpim1 496 537 IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 497 ibld(ji,jj) = MIN(mbkt(ji,jj), jk)538 ibld(ji,jj) = jk 498 539 ENDIF 499 540 END DO … … 504 545 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 505 546 ! 506 DO jj = 2, jpjm1507 DO ji = 2, jpim1508 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN547 CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 548 ! Alan: do we need zb_ml? 549 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 550 ! 510 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.511 551 ! 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 552 CALL zdf_osm_pycnocline_thickness( dh, zdh ) 567 553 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 ! 554 ! 555 ! Average over the depth of the mixed layer in the convective boundary layer 556 ! Alan: do we need zb_ml? 557 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 558 ! rotate mean currents and changes onto wind align co-ordinates 725 559 ! 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 560 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml ) 561 CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl ) 745 562 zuw_bse = 0._wp 746 563 zvw_bse = 0._wp 564 zwth_ent = 0._wp 565 zws_ent = 0._wp 747 566 DO jj = 2, jpjm1 748 567 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) 568 IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN 569 IF ( lconv(ji,jj) ) THEN 570 zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + & 571 & 1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ & 572 & ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln) 573 zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ & 574 & 2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj)) 575 IF ( zdb_ml(ji,jj) > 0._wp ) THEN 576 zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 577 zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 578 ENDIF 579 ELSE 580 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 581 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 754 582 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 583 ENDIF 759 584 END DO … … 764 589 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 765 590 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 591 CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext ) 592 CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc ) 593 CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc ) 847 594 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 848 595 ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 849 596 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 850 597 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 598 DO jj = 2, jpjm1 863 599 DO ji = 2, jpim1 … … 868 604 zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 869 605 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)606 zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln ) 871 607 ELSE 872 608 zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 873 609 zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 874 END IF875 END DO876 END DO610 END IF 611 END DO 612 END DO 877 613 ! 878 614 DO jj = 2, jpjm1 … … 889 625 ! pycnocline - if present linear profile 890 626 IF ( zdh(ji,jj) > 0._wp ) THEN 627 zgamma_b = 6.0 891 628 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 629 zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 893 630 ! 894 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 +zznd_pyc )631 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 895 632 ! 896 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 +zznd_pyc )633 zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc ) 897 634 END DO 635 IF ( ibld_ext == 0 ) THEN 636 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 637 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 638 ELSE 639 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 640 zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) ) 641 ENDIF 898 642 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)) 643 ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out 644 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6) 645 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5) 901 646 ! could be taken out, take account of entrainment represents as a diffusivity 902 647 ! should remove w from here, represents entrainment … … 908 653 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 909 654 END DO 655 656 IF ( ibld_ext == 0 ) THEN 657 zdiffut(ji,jj,ibld(ji,jj)) = 0._wp 658 zviscos(ji,jj,ibld(ji,jj)) = 0._wp 659 ELSE 660 zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 661 zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj)) 662 ENDIF 910 663 ENDIF ! end if ( lconv ) 911 !664 ! 912 665 END DO ! end of ji loop 913 666 END DO ! end of jj loop … … 952 705 END DO ! end of jj loop 953 706 954 955 707 ! 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 708 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 )709 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 ) 710 zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 ) 711 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 712 ELSEWHERE 961 713 zsc_uw_1 = zustar**2 962 zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)714 zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln) 963 715 ENDWHERE 964 716 IF(ln_dia_osm) THEN 717 IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu ) 718 IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv ) 719 END IF 965 720 DO jj = 2, jpjm1 966 721 DO ji = 2, jpim1 … … 1007 762 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & 1008 763 & * ( 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)764 zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 1010 765 ! non-gradient buoyancy terms 1011 766 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 775 END DO ! ji loop 1021 776 END DO ! jj loop 1022 1023 777 1024 778 WHERE ( lconv ) … … 1051 805 END DO ! jj loop 1052 806 807 IF(ln_dia_osm) THEN 808 IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu ) 809 IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 ) 810 END IF 1053 811 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 1054 812 … … 1089 847 END DO ! jj loop 1090 848 1091 1092 849 WHERE ( lconv ) 1093 850 zsc_uw_1 = zustar**2 … … 1134 891 END DO 1135 892 END DO 893 894 IF(ln_dia_osm) THEN 895 IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu ) 896 IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv ) 897 IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 ) 898 IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 ) 899 IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 ) 900 IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 ) 901 END IF 1136 902 ! 1137 903 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. … … 1165 931 END DO 1166 932 933 IF(ln_dia_osm) THEN 934 IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu ) 935 IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv ) 936 END IF 1167 937 ! pynocline contributions 1168 938 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small … … 1170 940 DO jj = 2, jpjm1 1171 941 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 942 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 943 DO jk= 2, ibld(ji,jj) 944 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 945 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 946 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 947 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 948 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) 949 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 950 END DO 951 END IF 952 END DO 1181 953 END DO 1182 954 … … 1185 957 DO jj=2, jpjm1 1186 958 DO ji = 2, jpim1 1187 IF ( lconv(ji,jj)) THEN959 IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN 1188 960 DO jk = 1, imld(ji,jj) - 1 1189 961 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) * znd962 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 963 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 1192 964 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 1193 965 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd … … 1195 967 DO jk = imld(ji,jj), ibld(ji,jj) 1196 968 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 )969 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 970 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 1199 971 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 1200 972 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 1201 973 END DO 1202 974 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 975 976 ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 977 ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 978 ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 979 ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp 1207 980 END DO ! ji loop 1208 981 END DO ! jj loop 1209 982 1210 983 IF(ln_dia_osm) THEN 984 IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu ) 985 IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv ) 986 IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse ) 987 IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse ) 988 IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc ) 989 IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc ) 990 IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos ) 991 END IF 1211 992 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1212 993 ! Need to put in code for contributions that are applied explicitly to … … 1232 1013 IF(ln_dia_osm) THEN 1233 1014 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1015 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1016 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1234 1017 END IF 1235 1018 … … 1286 1069 END IF ! ln_convmix = .true. 1287 1070 1071 1072 1073 IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing 1074 DO jj = 2 , jpjm1 1075 DO ji = 2, jpim1 1076 IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL. 1077 ! Calculate MLE flux profiles 1078 ! DO jk = 1, mld_prof(ji,jj) 1079 ! znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1080 ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + & 1081 ! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1082 ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + & 1083 ! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 ) 1084 ! END DO 1085 ! Calculate MLE flux contribution from surface fluxes 1086 DO jk = 1, ibld(ji,jj) 1087 znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln) 1088 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd ) 1089 ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd ) 1090 END DO 1091 DO jk = 1, mld_prof(ji,jj) 1092 znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln) 1093 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd ) 1094 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd ) 1095 END DO 1096 ! Viscosity for MLEs 1097 DO jk = ibld(ji,jj), mld_prof(ji,jj) 1098 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) ) 1099 END DO 1100 ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj) 1101 jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj)) 1102 jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj)) 1103 DO jk = mld_prof(ji,jj), jl 1104 zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * & 1105 & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / & 1106 & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln)) 1107 END DO 1108 ENDIF 1109 END DO 1110 END DO 1111 ENDIF 1112 1113 IF(ln_dia_osm) THEN 1114 IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc ) 1115 IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc ) 1116 IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc ) 1117 END IF 1118 1119 1288 1120 ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids 1289 CALL lbc_lnk( 'zdfosm',zviscos(:,:,:), 'W', 1. )1121 !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. ) 1290 1122 1291 1123 ! GN 25/8: need to change tmask --> wmask … … 1316 1148 END DO 1317 1149 END DO 1150 ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged) 1151 CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. ) 1318 1152 ! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged) 1319 1153 ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged) 1320 1154 CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., & 1321 & ghamu, 'U', 1. , ghamv, 'V',1. )1322 1323 1155 & ghamu, 'U', -1. , ghamv, 'V', -1. ) 1156 1157 IF(ln_dia_osm) THEN 1324 1158 SELECT CASE (nn_osm_wave) 1325 1159 ! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1). … … 1330 1164 ! Stokes drift read in from sbcwave (=2). 1331 1165 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 1166 IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift 1167 IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift 1168 IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period 1169 IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height 1170 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 1171 IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) ) ! significant wave height from NP spectrum 1172 IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10 1334 1173 IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 1335 1174 & SQRT(ut0sd**2 + vt0sd**2 ) ) … … 1342 1181 IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0> 1343 1182 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 1183 IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k 1184 IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base 1185 IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base 1186 IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base 1187 IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base 1188 IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base 1189 IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth 1190 IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth 1345 1191 IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth 1346 1192 IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points … … 1348 1194 IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale 1349 1195 IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale 1196 IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale 1197 IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir # 1350 1198 IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 1351 1199 IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 1352 1200 IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine 1353 1201 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 1202 IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine 1203 IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine 1355 1204 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 1205 IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux 1206 IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux 1207 IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux 1208 IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux 1209 IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML 1210 1211 IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth 1212 IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth 1213 IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux 1214 IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML 1215 IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k 1216 IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt 1217 IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt 1218 IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt 1219 IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt 1220 IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt 1221 IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt 1222 IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1223 IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt 1224 1359 1225 END IF 1360 ! Lateral boundary conditions on p_avt (sign unchanged) 1361 CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. ) 1226 1227 CONTAINS 1228 1229 1230 ! Alan: do we need zb? 1231 SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv ) 1232 !!--------------------------------------------------------------------- 1233 !! *** ROUTINE zdf_vertical_average *** 1234 !! 1235 !! ** Purpose : Determines vertical averages from surface to jnlev. 1236 !! 1237 !! ** Method : Averages are calculated from the surface to jnlev. 1238 !! The external level used to calculate differences is ibld+ibld_ext 1239 !! 1240 !!---------------------------------------------------------------------- 1241 1242 INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over. 1243 1244 ! Alan: do we need zb? 1245 REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity 1246 REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components 1247 REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL 1248 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components. 1249 1250 INTEGER :: jk, ji, jj 1251 REAL(wp) :: zthick, zthermal, zbeta 1252 1253 1254 zt = 0._wp 1255 zs = 0._wp 1256 zu = 0._wp 1257 zv = 0._wp 1258 DO jj = 2, jpjm1 ! Vertical slab 1259 DO ji = 2, jpim1 1260 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1261 zbeta = rab_n(ji,jj,1,jp_sal) 1262 ! average over depth of boundary layer 1263 zthick = epsln 1264 DO jk = 2, jnlev_av(ji,jj) 1265 zthick = zthick + e3t_n(ji,jj,jk) 1266 zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 1267 zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 1268 zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) & 1269 & * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) & 1270 & / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) ) 1271 zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) & 1272 & * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) & 1273 & / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) ) 1274 END DO 1275 zt(ji,jj) = zt(ji,jj) / zthick 1276 zs(ji,jj) = zs(ji,jj) / zthick 1277 zu(ji,jj) = zu(ji,jj) / zthick 1278 zv(ji,jj) = zv(ji,jj) / zthick 1279 ! Alan: do we need zb? 1280 zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem) 1281 zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal) 1282 zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) & 1283 & / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) 1284 zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) & 1285 & / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) 1286 zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj) 1287 END DO 1288 END DO 1289 END SUBROUTINE zdf_osm_vertical_average 1290 1291 SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv ) 1292 !!--------------------------------------------------------------------- 1293 !! *** ROUTINE zdf_velocity_rotation *** 1294 !! 1295 !! ** Purpose : Rotates frame of reference of averaged velocity components. 1296 !! 1297 !! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w 1298 !! 1299 !!---------------------------------------------------------------------- 1300 1301 REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle 1302 REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current 1303 REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline 1304 1305 INTEGER :: ji, jj 1306 REAL(wp) :: ztemp 1307 1308 DO jj = 2, jpjm1 1309 DO ji = 2, jpim1 1310 ztemp = zu(ji,jj) 1311 zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj) 1312 zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj) 1313 ztemp = zdu(ji,jj) 1314 zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj) 1315 zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj) 1316 END DO 1317 END DO 1318 END SUBROUTINE zdf_osm_velocity_rotation 1319 1320 SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz ) 1321 !!--------------------------------------------------------------------- 1322 !! *** ROUTINE zdf_osm_external_gradients *** 1323 !! 1324 !! ** Purpose : Calculates the gradients below the OSBL 1325 !! 1326 !! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient. 1327 !! 1328 !!---------------------------------------------------------------------- 1329 1330 REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy. 1331 1332 INTEGER :: jj, ji, jkb, jkb1 1333 REAL(wp) :: zthermal, zbeta 1334 1335 1336 DO jj = 2, jpjm1 1337 DO ji = 2, jpim1 1338 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1339 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 1340 zbeta = rab_n(ji,jj,1,jp_sal) 1341 jkb = ibld(ji,jj) + ibld_ext 1342 jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 1343 zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) & 1344 & / e3t_n(ji,jj,ibld(ji,jj)) 1345 zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) & 1346 & / e3t_n(ji,jj,ibld(ji,jj)) 1347 zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 1348 END IF 1349 END DO 1350 END DO 1351 END SUBROUTINE zdf_osm_external_gradients 1352 1353 SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz ) 1354 1355 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline 1356 1357 INTEGER :: jk, jj, ji 1358 REAL(wp) :: ztgrad, zsgrad, zbgrad 1359 REAL(wp) :: zgamma_b_nd, zgamma_c, znd 1360 REAL(wp) :: zzeta_s=0.3, ztmp 1361 1362 DO jj = 2, jpjm1 1363 DO ji = 2, jpim1 1364 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1365 IF ( lconv(ji,jj) ) THEN ! convective conditions 1366 IF ( zdbdz_ext(ji,jj) > 0._wp .AND. & 1367 & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh & 1368 & .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0 1369 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1370 ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj) 1371 zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj) 1372 zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj) 1373 zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln) 1374 zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /& 1375 & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check 1376 DO jk = 2, ibld(ji,jj)+ibld_ext 1377 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp 1378 IF ( znd <= zzeta_s ) THEN 1379 zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * & 1380 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1381 zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * & 1382 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1383 zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * & 1384 & EXP( -6.0 * ( znd -zzeta_s )**2 ) 1385 ELSE 1386 zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1387 zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1388 zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 ) 1389 ENDIF 1390 END DO 1391 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing 1392 ELSE 1393 ! stable conditions 1394 ! if pycnocline profile only defined when depth steady of increasing. 1395 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady. 1396 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1397 IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline 1398 ztmp = 1._wp/MAX(zhbl(ji,jj), epsln) 1399 ztgrad = zdt_bl(ji,jj) * ztmp 1400 zsgrad = zds_bl(ji,jj) * ztmp 1401 zbgrad = zdb_bl(ji,jj) * ztmp 1402 DO jk = 2, ibld(ji,jj) 1403 znd = gdepw_n(ji,jj,jk) * ztmp 1404 zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1405 zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1406 zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 1407 END DO 1408 ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 1409 ztmp = 1._wp/MAX(zdh(ji,jj), epsln) 1410 ztgrad = zdt_bl(ji,jj) * ztmp 1411 zsgrad = zds_bl(ji,jj) * ztmp 1412 zbgrad = zdb_bl(ji,jj) * ztmp 1413 DO jk = 2, ibld(ji,jj) 1414 znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp 1415 zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1416 zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1417 zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 1418 END DO 1419 ENDIF ! IF (zhol >=0.5) 1420 ENDIF ! IF (zdb_bl> 0.) 1421 ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero 1422 ENDIF ! IF (lconv) 1423 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1424 END DO 1425 END DO 1426 1427 END SUBROUTINE zdf_osm_pycnocline_scalar_profiles 1428 1429 SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz ) 1430 !!--------------------------------------------------------------------- 1431 !! *** ROUTINE zdf_osm_pycnocline_shear_profiles *** 1432 !! 1433 !! ** Purpose : Calculates velocity shear in the pycnocline 1434 !! 1435 !! ** Method : 1436 !! 1437 !!---------------------------------------------------------------------- 1438 1439 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz 1440 1441 INTEGER :: jk, jj, ji 1442 REAL(wp) :: zugrad, zvgrad, znd 1443 REAL(wp) :: zzeta_v = 0.45 1362 1444 ! 1363 END SUBROUTINE zdf_osm 1445 DO jj = 2, jpjm1 1446 DO ji = 2, jpim1 1447 ! 1448 IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN 1449 IF ( lconv (ji,jj) ) THEN 1450 ! Unstable conditions 1451 zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / & 1452 & ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * & 1453 & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 )) 1454 !Alan is this right? 1455 zvgrad = ( 0.7 * zdv_ml(ji,jj) + & 1456 & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / & 1457 & ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird + epsln ) & 1458 & )/ (zdh(ji,jj) + epsln ) 1459 DO jk = 2, ibld(ji,jj) - 1 + ibld_ext 1460 znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v 1461 IF ( znd <= 0.0 ) THEN 1462 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd ) 1463 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd ) 1464 ELSE 1465 zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd ) 1466 zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd ) 1467 ENDIF 1468 END DO 1469 ELSE 1470 ! stable conditions 1471 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 1472 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 1473 DO jk = 2, ibld(ji,jj) 1474 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 1475 IF ( znd < 1.0 ) THEN 1476 zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 1477 ELSE 1478 zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 1479 ENDIF 1480 zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 1481 END DO 1482 ENDIF 1483 ! 1484 END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) 1485 END DO 1486 END DO 1487 END SUBROUTINE zdf_osm_pycnocline_shear_profiles 1488 1489 SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 ) 1490 !!--------------------------------------------------------------------- 1491 !! *** ROUTINE zdf_osm_calculate_dhdt *** 1492 !! 1493 !! ** Purpose : Calculates the rate at which hbl changes. 1494 !! 1495 !! ** Method : 1496 !! 1497 !!---------------------------------------------------------------------- 1498 1499 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! Rate of change of hbl 1500 1501 INTEGER :: jj, ji 1502 REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert 1503 REAL(wp) :: zvel_max, zwb_min 1504 REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 1505 REAL(wp) :: zzeta_m = 0.3 1506 REAL(wp) :: zgamma_c = 2.0 1507 REAL(wp) :: zdhoh = 0.1 1508 REAL(wp) :: alpha_bc = 0.5 1509 1510 DO jj = 2, jpjm1 1511 DO ji = 2, jpim1 1512 IF ( lconv(ji,jj) ) THEN ! Convective 1513 ! Alan is this right? Yes, it's a bit different from the previous relationship 1514 ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) & 1515 ! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1516 zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln 1517 zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 ) 1518 zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 1519 zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 1520 zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 1521 & * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 1522 1523 zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 1524 & - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) & 1525 & + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 & 1526 & - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj) 1527 ! 1528 zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj) 1529 1530 IF ( ln_osm_mle ) THEN 1531 ! 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 + & 1532 ! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL 1533 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1534 zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * & 1535 (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 1536 ELSE 1537 zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 1538 ENDIF 1539 zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1540 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 1541 ! OSBL is deepening, entrainment > restratification 1542 IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1543 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1544 ELSE 1545 zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15) 1546 ENDIF 1547 ELSE 1548 ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 1549 zdhdt(ji,jj) = - zvel_mle(ji,jj) 1550 1551 1552 ENDIF 1553 1554 ELSE 1555 ! Fox-Kemper not used. 1556 1557 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) / & 1558 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1559 zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 1560 ! added ajgn 23 July as temporay fix 1561 1562 ENDIF 1563 1564 zdhdt_2(ji,jj) = 0._wp 1565 1566 ! commented out ajgn 23 July as temporay fix 1567 ! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN 1568 ! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary. 1569 ! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1570 ! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj) 1571 ! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) 1572 ! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj) 1573 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min ) 1574 ! ! Alan no idea what this should be? 1575 ! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) & 1576 ! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) & 1577 ! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj) 1578 ! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) ) 1579 ! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN 1580 ! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj) 1581 ! ENDIF 1582 ELSE ! Stable 1583 zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj) 1584 zdhdt_2(ji,jj) = 0._wp 1585 IF ( zdhdt(ji,jj) < 0._wp ) THEN 1586 ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 1587 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj) 1588 ELSE 1589 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) ) 1590 ENDIF 1591 zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert 1592 ENDIF 1593 END DO 1594 END DO 1595 END SUBROUTINE zdf_osm_calculate_dhdt 1596 1597 SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 ) 1598 !!--------------------------------------------------------------------- 1599 !! *** ROUTINE zdf_osm_timestep_hbl *** 1600 !! 1601 !! ** Purpose : Increments hbl. 1602 !! 1603 !! ** Method : If thechange in hbl exceeds one model level the change is 1604 !! is calculated by moving down the grid, changing the buoyancy 1605 !! jump. This is to ensure that the change in hbl does not 1606 !! overshoot a stable layer. 1607 !! 1608 !!---------------------------------------------------------------------- 1609 1610 1611 REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl. 1612 1613 INTEGER :: jk, jj, ji, jm 1614 REAL(wp) :: zhbl_s, zvel_max, zdb 1615 REAL(wp) :: zthermal, zbeta 1616 1617 DO jj = 2, jpjm1 1618 DO ji = 2, jpim1 1619 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 1620 ! 1621 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 1622 ! 1623 zhbl_s = hbl(ji,jj) 1624 jm = imld(ji,jj) 1625 zthermal = rab_n(ji,jj,1,jp_tem) 1626 zbeta = rab_n(ji,jj,1,jp_sal) 1627 1628 1629 IF ( lconv(ji,jj) ) THEN 1630 !unstable 1631 1632 IF( ln_osm_mle ) THEN 1633 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 1634 ELSE 1635 1636 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) / & 1637 & ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1638 1639 ENDIF 1640 1641 DO jk = imld(ji,jj), ibld(ji,jj) 1642 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 1643 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), & 1644 & 0.0 ) + zvel_max 1645 1646 1647 IF ( ln_osm_mle ) THEN 1648 zhbl_s = zhbl_s + MIN( & 1649 & 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) ), & 1650 & e3w_n(ji,jj,jm) ) 1651 ELSE 1652 zhbl_s = zhbl_s + MIN( & 1653 & rn_rdt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 1654 & e3w_n(ji,jj,jm) ) 1655 ENDIF 1656 1657 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1658 1659 IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 1660 END DO 1661 hbl(ji,jj) = zhbl_s 1662 ibld(ji,jj) = jm 1663 ELSE 1664 ! stable 1665 DO jk = imld(ji,jj), ibld(ji,jj) 1666 zdb = MAX( & 1667 & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )& 1668 & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),& 1669 & 0.0 ) + & 1670 & 2.0 * zvstr(ji,jj)**2 / zhbl_s 1671 1672 ! Alan is thuis right? I have simply changed hbli to hbl 1673 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) ) 1674 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) ) ) * & 1675 & zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) ) 1676 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj) 1677 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) ) 1678 1679 zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol) 1680 IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 1681 END DO 1682 ENDIF ! IF ( lconv ) 1683 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) ) 1684 ibld(ji,jj) = MAX(jm, 4 ) 1685 ELSE 1686 ! change zero or one model level. 1687 hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) ) 1688 ENDIF 1689 zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 1690 END DO 1691 END DO 1692 1693 END SUBROUTINE zdf_osm_timestep_hbl 1694 1695 SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 1696 !!--------------------------------------------------------------------- 1697 !! *** ROUTINE zdf_osm_pycnocline_thickness *** 1698 !! 1699 !! ** Purpose : Calculates thickness of the pycnocline 1700 !! 1701 !! ** Method : The thickness is calculated from a prognostic equation 1702 !! that relaxes the pycnocine thickness to a diagnostic 1703 !! value. The time change is calculated assuming the 1704 !! thickness relaxes exponentially. This is done to deal 1705 !! with large timesteps. 1706 !! 1707 !!---------------------------------------------------------------------- 1708 1709 REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness. 1710 ! 1711 INTEGER :: jj, ji 1712 INTEGER :: inhml 1713 REAL(wp) :: zari, ztau, zddhdt 1714 1715 1716 DO jj = 2, jpjm1 1717 DO ji = 2, jpim1 1718 IF ( lconv(ji,jj) ) THEN 1719 1720 IF( ln_osm_mle ) THEN 1721 IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN 1722 ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F 1723 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1724 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1725 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1726 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1727 ELSE ! unstable 1728 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1729 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1730 ENDIF 1731 ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1732 zddhdt = zari * hbl(ji,jj) 1733 ELSE 1734 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1735 zddhdt = 0.2 * hbl(ji,jj) 1736 ENDIF 1737 ELSE 1738 ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1739 zddhdt = 0.2 * hbl(ji,jj) 1740 ENDIF 1741 ELSE ! ln_osm_mle 1742 IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN 1743 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability 1744 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1745 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1746 ELSE ! unstable 1747 zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / & 1748 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 ) 1749 ENDIF 1750 ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird 1751 zddhdt = zari * hbl(ji,jj) 1752 ELSE 1753 ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 1754 zddhdt = 0.2 * hbl(ji,jj) 1755 ENDIF 1756 1757 END IF ! ln_osm_mle 1758 1759 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1760 ! Alan: this hml is never defined or used 1761 ELSE ! IF (lconv) 1762 ztau = hbl(ji,jj) / zvstr(ji,jj) 1763 IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here 1764 ! boundary layer deepening 1765 IF ( zdb_bl(ji,jj) > 0._wp ) THEN 1766 ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 1767 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 1768 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 1769 zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj) 1770 ELSE 1771 zddhdt = 0.2 * hbl(ji,jj) 1772 ENDIF 1773 ELSE ! IF(dhdt < 0) 1774 zddhdt = 0.2 * hbl(ji,jj) 1775 ENDIF ! IF (dhdt >= 0) 1776 dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) ) 1777 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 1778 ! Alan: this hml is never defined or used -- do we need it? 1779 ENDIF ! IF (lconv) 1780 1781 hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) 1782 inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 1783 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3) 1784 zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 1785 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 1786 END DO 1787 END DO 1788 1789 END SUBROUTINE zdf_osm_pycnocline_thickness 1790 1791 1792 SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle ) 1793 !!---------------------------------------------------------------------- 1794 !! *** ROUTINE zdf_osm_horizontal_gradients *** 1795 !! 1796 !! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization. 1797 !! 1798 !! ** Method : 1799 !! 1800 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1801 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1802 1803 1804 REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points 1805 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == ! 1806 REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy 1807 1808 INTEGER :: ji, jj, jk ! dummy loop indices 1809 INTEGER :: ii, ij, ik, ikmax ! local integers 1810 REAL(wp) :: zc 1811 REAL(wp) :: zN2_c ! local buoyancy difference from 10m value 1812 REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH 1813 REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv 1814 REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv 1815 !!---------------------------------------------------------------------- 1816 ! 1817 ! !== MLD used for MLE ==! 1818 1819 mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point 1820 zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 1821 zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria 1822 DO jk = nlb10, jpkm1 1823 DO jj = 1, jpj ! Mixed layer level: w-level 1824 DO ji = 1, jpi 1825 ikt = mbkt(ji,jj) 1826 zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 1827 IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1828 END DO 1829 END DO 1830 END DO 1831 DO jj = 1, jpj 1832 DO ji = 1, jpi 1833 mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj)) 1834 zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1835 END DO 1836 END DO 1837 ! ensure mld_prof .ge. ibld 1838 ! 1839 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 1840 ! 1841 ztm(:,:) = 0._wp 1842 zsm(:,:) = 0._wp 1843 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 1844 DO jj = 1, jpj 1845 DO ji = 1, jpi 1846 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 1847 ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) 1848 zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal) 1849 END DO 1850 END DO 1851 END DO 1852 ! average temperature and salinity. 1853 ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1854 zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 1855 ! calculate horizontal gradients at u & v points 1856 1857 DO jj = 2, jpjm1 1858 DO ji = 1, jpim1 1859 zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1860 zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj) 1861 zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj)) 1862 ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) ) 1863 ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) ) 1864 END DO 1865 END DO 1866 1867 DO jj = 1, jpjm1 1868 DO ji = 2, jpim1 1869 zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1870 zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 1871 zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj)) 1872 ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) ) 1873 ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) ) 1874 END DO 1875 END DO 1876 1877 CALL eos_rab(ztsm_midu, zmld_midu, zabu) 1878 CALL eos_rab(ztsm_midv, zmld_midv, zabv) 1879 1880 DO jj = 2, jpjm1 1881 DO ji = 1, jpim1 1882 dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal)) 1883 END DO 1884 END DO 1885 DO jj = 1, jpjm1 1886 DO ji = 2, jpim1 1887 dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal)) 1888 END DO 1889 END DO 1890 1891 END SUBROUTINE zdf_osm_zmld_horizontal_gradients 1892 SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle ) 1893 !!---------------------------------------------------------------------- 1894 !! *** ROUTINE zdf_osm_mle_parameters *** 1895 !! 1896 !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity. 1897 !! 1898 !! ** Method : 1899 !! 1900 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 1901 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 1902 1903 REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle 1904 INTEGER :: ji, jj, jk ! dummy loop indices 1905 INTEGER :: ii, ij, ik ! local integers 1906 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 1907 REAL(wp) :: zdb_mle, ztmp, zdbds_mle 1908 1909 mld_prof(:,:) = 4 1910 DO jk = 5, jpkm1 1911 DO jj = 2, jpjm1 1912 DO ji = 2, jpim1 1913 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1914 END DO 1915 END DO 1916 END DO 1917 ! DO jj = 2, jpjm1 1918 ! DO ji = 1, jpim1 1919 ! zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj)) 1920 ! END DO 1921 ! END DO 1922 ! Timestep mixed layer eddy depth. 1923 DO jj = 2, jpjm1 1924 DO ji = 2, jpim1 1925 zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG 1926 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 1927 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. 1928 ELSE 1929 ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10 1930 ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN 1931 ! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1932 ! ELSE 1933 ! 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 1934 ! ENDIF 1935 IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 1936 hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau 1937 ELSE 1938 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 1939 ENDIF 1940 ENDIF 1941 hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj)) 1942 hmle(ji,jj) = MIN(hmle(ji,jj), 1.2*zmld(ji,jj)) 1943 END DO 1944 END DO 1945 1946 mld_prof = 4 1947 DO jk = 5, jpkm1 1948 DO jj = 2, jpjm1 1949 DO ji = 2, jpim1 1950 IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 1951 END DO 1952 END DO 1953 END DO 1954 DO jj = 2, jpjm1 1955 DO ji = 2, jpim1 1956 zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj)) 1957 END DO 1958 END DO 1959 ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE. 1960 1961 DO jj = 2, jpjm1 1962 DO ji = 2, jpim1 1963 IF ( lconv(ji,jj) ) THEN 1964 ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 1965 ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) & 1966 ! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) ) 1967 ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) & 1968 ! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) ) 1969 zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) & 1970 & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) ) 1971 zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle 1972 ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt. 1973 zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 1974 zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf 1975 ENDIF 1976 END DO 1977 END DO 1978 END SUBROUTINE zdf_osm_mle_parameters 1979 1980 END SUBROUTINE zdf_osm 1364 1981 1365 1982 … … 1378 1995 INTEGER :: ios ! local integer 1379 1996 INTEGER :: ji, jj, jk ! dummy loop indices 1997 REAL z1_t2 1380 1998 !! 1381 1999 NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 1382 2000 & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 1383 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv 2001 & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 2002 & ,ln_osm_mle 2003 ! Namelist for Fox-Kemper parametrization. 2004 NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 2005 & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau 2006 1384 2007 !!---------------------------------------------------------------------- 1385 2008 ! … … 1397 2020 WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation' 1398 2021 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 2022 WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters' 2023 WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la 2024 WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle 1401 2025 WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la 1402 2026 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_dstokes2027 WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes 1404 2028 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 1405 2029 WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave … … 1423 2047 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1424 2048 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl 2049 2050 IF( ln_osm_mle ) THEN 2051 ! Initialise Fox-Kemper parametrization 2052 REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme 2053 READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903) 2054 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist') 2055 2056 REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme 2057 READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 ) 2058 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist') 2059 IF(lwm) WRITE ( numond, namosm_mle ) 2060 2061 IF(lwp) THEN ! Namelist print 2062 WRITE(numout,*) 2063 WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)' 2064 WRITE(numout,*) '~~~~~~~~~~~~~' 2065 WRITE(numout,*) ' Namelist namosm_mle : ' 2066 WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle 2067 WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce 2068 WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm' 2069 WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's' 2070 WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg' 2071 WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c 2072 WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s' 2073 WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's' 2074 ENDIF ! 2075 ENDIF 2076 ! 2077 IF(lwp) THEN 2078 WRITE(numout,*) 2079 IF( ln_osm_mle ) THEN 2080 WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation' 2081 IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation' 2082 IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation' 2083 ELSE 2084 WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation' 2085 ENDIF 2086 ENDIF 2087 ! 2088 IF( ln_osm_mle ) THEN ! MLE initialisation 2089 ! 2090 rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria 2091 IF(lwp) WRITE(numout,*) 2092 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' 2093 IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3' 2094 ! 2095 IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation ! 2096 ! 2097 ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation 2098 rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) ) 2099 ! 2100 ENDIF 2101 ! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 2102 z1_t2 = 2.e-5 2103 do jj=1,jpj 2104 do ji = 1,jpi 2105 r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2) 2106 end do 2107 end do 2108 ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj ) 2109 ! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 2110 ! 2111 ENDIF 2112 2113 call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle 2114 1426 2115 1427 2116 IF( ln_zdfddm) THEN … … 1512 2201 CALL iom_set_rstw_var_active('wn') 1513 2202 CALL iom_set_rstw_var_active('hbl') 1514 CALL iom_set_rstw_var_active('hbli') 2203 CALL iom_set_rstw_var_active('dh') 2204 IF( ln_osm_mle ) THEN 2205 CALL iom_set_rstw_var_active('hmle') 2206 END IF 1515 2207 ENDIF 1516 2208 END SUBROUTINE zdf_osm_init … … 1530 2222 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1531 2223 1532 INTEGER :: id1, id2 ! iom enquiry index2224 INTEGER :: id1, id2, id3 ! iom enquiry index 1533 2225 INTEGER :: ji, jj, jk ! dummy loop indices 1534 2226 INTEGER :: iiki, ikt ! local integer … … 1536 2228 REAL(wp) :: zN2_c ! local scalar 1537 2229 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)2230 INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top) 1539 2231 !!---------------------------------------------------------------------- 1540 2232 ! … … 1551 2243 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1552 2244 END IF 2245 1553 2246 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) 1554 id2 = iom_varid( numror, ' hbli' , ldstop = .FALSE. )2247 id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. ) 1555 2248 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1556 2249 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 & hbli read from restart file' 2250 CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios ) 2251 WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file' 2252 IF( ln_osm_mle ) THEN 2253 id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. ) 2254 IF( id3 > 0) THEN 2255 CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios ) 2256 WRITE(numout,*) ' ===>>>> : hmle read from restart file' 2257 ELSE 2258 WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl' 2259 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2260 END IF 2261 END IF 1559 2262 RETURN 1560 ELSE ! 'hbl' & ' hbli' not in restart file, recalculate2263 ELSE ! 'hbl' & 'dh' not in restart file, recalculate 1561 2264 WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification' 1562 2265 END IF … … 1568 2271 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 1569 2272 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1570 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn , ldxios = lwxios ) 1571 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1572 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) 2273 CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn, ldxios = lwxios ) 2274 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios ) 2275 CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios ) 2276 IF( ln_osm_mle ) THEN 2277 CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios ) 2278 END IF 1573 2279 RETURN 1574 2280 END IF … … 1578 2284 !!----------------------------------------------------------------------------- 1579 2285 IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification' 1580 ALLOCATE( imld_rst(jpi,jpj) )1581 2286 ! w-level of the mixing and mixed layers 1582 2287 CALL eos_rab( tsn, rab_n ) … … 1599 2304 DO jj = 1, jpj 1600 2305 DO ji = 1, jpi 1601 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 2306 iiki = MAX(4,imld_rst(ji,jj)) 2307 hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth 2308 dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth 1603 2309 END DO 1604 2310 END DO 1605 hbl = MAX(hbl,epsln) 1606 hbli(:,:) = hbl(:,:) 1607 DEALLOCATE( imld_rst ) 2311 1608 2312 WRITE(numout,*) ' ===>>>> : hbl computed from stratification' 2313 2314 IF( ln_osm_mle ) THEN 2315 hmle(:,:) = hbl(:,:) ! Initialise MLE depth. 2316 WRITE(numout,*) ' ===>>>> : hmle set = to hbl' 2317 END IF 2318 2319 wn(:,:,:) = 0._wp 2320 WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially' 1609 2321 END SUBROUTINE osm_rst 1610 2322 … … 1634 2346 ENDIF 1635 2347 1636 ! add non-local temperature and salinity flux1637 2348 DO jk = 1, jpkm1 1638 2349 DO jj = 2, jpjm1 … … 1648 2359 END DO 1649 2360 1650 1651 ! save the non-local tracer flux trends for diagnostic 2361 ! save the non-local tracer flux trends for diagnostics 1652 2362 IF( l_trdtra ) THEN 1653 2363 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 1654 2364 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 )2365 2366 CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt ) 2367 CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds ) 1658 2368 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1659 2369 ENDIF … … 1723 2433 1724 2434 !!====================================================================== 2435 1725 2436 END MODULE zdfosm -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfphy.F90
r11715 r13402 56 56 !!---------------------------------------------------------------------- 57 57 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 58 !! $Id $58 !! $Id: zdfphy.F90 12178 2019-12-11 11:02:38Z agn $ 59 59 !! Software governed by the CeCILL license (see ./LICENSE) 60 60 !!---------------------------------------------------------------------- … … 172 172 IF( ln_zdfosm .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' ) 173 173 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 174 IF( lk_top .AND. ln_zdfosm ) CALL ctl_ stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )174 IF( lk_top .AND. ln_zdfosm ) CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' ) 175 175 IF(lwp) THEN 176 176 WRITE(numout,*) -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/nemogcm.F90
r11715 r13402 394 394 ! !-----------------------------------------! 395 395 CALL mpp_init 396 WRITE(numout,*)'NEMOGCM 1: nstop = ',nstop 396 397 397 398 ! Now we know the dimensions of the grid and numout has been set: we can allocate arrays … … 405 406 ! 406 407 ! ! General initialization 408 WRITE(numout,*)'NEMOGCM 2: nstop = ',nstop 407 409 IF( ln_timing ) CALL timing_init ! timing 408 410 IF( ln_timing ) CALL timing_start( 'nemo_init') … … 411 413 CALL eos_init ! Equation of state 412 414 IF( lk_c1d ) CALL c1d_init ! 1D column configuration 415 WRITE(numout,*)'NEMOGCM 3: nstop = ',nstop 413 416 CALL wad_init ! Wetting and drying options 414 417 CALL dom_init("OPA") ! Domain 415 418 IF( ln_crs ) CALL crs_init ! coarsened grid: domain initialization 416 419 IF( ln_ctl ) CALL prt_ctl_init ! Print control 420 WRITE(numout,*)'NEMOGCM 4: nstop = ',nstop 417 421 418 422 CALL diurnal_sst_bulk_init ! diurnal sst … … 433 437 434 438 CALL istate_init ! ocean initial state (Dynamics and tracers) 439 WRITE(numout,*)'NEMOGCM 5: nstop = ',nstop 435 440 436 441 ! ! external forcing … … 438 443 CALL sbc_init ! surface boundary conditions (including sea-ice) 439 444 CALL bdy_init ! Open boundaries initialisation 445 WRITE(numout,*)'NEMOGCM 6: nstop = ',nstop 440 446 441 447 ! ! Ocean physics 442 448 CALL zdf_phy_init ! Vertical physics 449 WRITE(numout,*)'NEMOGCM 7: nstop = ',nstop 443 450 444 451 ! ! Lateral physics … … 446 453 CALL ldf_eiv_init ! eddy induced velocity param. 447 454 CALL ldf_dyn_init ! Lateral ocean momentum physics 455 WRITE(numout,*)'NEMOGCM 8: nstop = ',nstop 448 456 449 457 ! ! Active tracers … … 451 459 CALL tra_bbc_init ! bottom heat flux 452 460 CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme 461 WRITE(numout,*)'NEMOGCM 9: nstop = ',nstop 453 462 CALL tra_dmp_init ! internal tracer damping 454 463 CALL tra_adv_init ! horizontal & vertical advection 455 464 CALL tra_ldf_init ! lateral mixing 465 WRITE(numout,*)'NEMOGCM 10: nstop = ',nstop 456 466 457 467 ! ! Dynamics … … 459 469 CALL dyn_adv_init ! advection (vector or flux form) 460 470 CALL dyn_vor_init ! vorticity term including Coriolis 471 WRITE(numout,*)'NEMOGCM 11: nstop = ',nstop 461 472 CALL dyn_ldf_init ! lateral mixing 462 473 CALL dyn_hpg_init ! horizontal gradient of Hydrostatic pressure 463 474 CALL dyn_spg_init ! surface pressure gradient 475 WRITE(numout,*)'NEMOGCM 12: nstop = ',nstop 464 476 465 477 #if defined key_top … … 467 479 CALL trc_init 468 480 #endif 481 WRITE(numout,*)'NEMOGCM 13: nstop = ',nstop 469 482 IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing 470 483 471 484 ! ! Icebergs 472 485 CALL icb_init( rdt, nit000) ! initialise icebergs instance 486 WRITE(numout,*)'NEMOGCM 14: nstop = ',nstop 473 487 474 488 ! ! Misc. options 475 489 CALL sto_par_init ! Stochastic parametrization 476 490 IF( ln_sto_eos ) CALL sto_pts_init ! RRandom T/S fluctuations 491 WRITE(numout,*)'NEMOGCM 15: nstop = ',nstop 477 492 478 493 ! ! Diagnostics … … 480 495 IF( ln_diacfl ) CALL dia_cfl_init ! Initialise CFL diagnostics 481 496 CALL dia_ptr_init ! Poleward TRansports initialization 497 WRITE(numout,*)'NEMOGCM 16: nstop = ',nstop 482 498 CALL dia_dct_init ! Sections tranports 483 499 CALL dia_hsb_init ! heat content, salt content and volume budgets 484 500 CALL trd_init ! Mixed-layer/Vorticity/Integral constraints trends 485 501 CALL dia_obs_init ! Initialize observational data 502 WRITE(numout,*)'NEMOGCM 17: nstop = ',nstop 486 503 CALL dia_tmb_init ! TMB outputs 487 504 CALL dia_25h_init ! 25h mean outputs 488 505 CALL dia_harm_init ! tidal harmonics outputs 489 506 IF( ln_diaobs ) CALL dia_obs( nit000-1 ) ! Observation operator for restart 507 WRITE(numout,*)'NEMOGCM 18: nstop = ',nstop 490 508 491 509 ! ! Assimilation increments -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/step.F90
r11715 r13402 101 101 IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 102 102 CALL iom_init( cxios_context ) ! for model grid (including passible AGRIF zoom) 103 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid 103 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid 104 CALL dia_hth_init ! extra ML depth diagnostics, thermocline depths 104 105 ENDIF 105 106 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) … … 205 206 IF( ln_floats ) CALL flo_stp ( kstp ) ! drifting Floats 206 207 IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics 207 IF( l k_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth)208 IF( ll_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 208 209 IF( ln_diadct ) CALL dia_dct ( kstp ) ! Transports 209 210 CALL dia_ar5 ( kstp ) ! ar5 diag
Note: See TracChangeset
for help on using the changeset viewer.