Changeset 7990
- Timestamp:
- 2017-04-29T17:24:54+02:00 (7 years ago)
- Location:
- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM
- Files:
-
- 27 edited
- 2 moved
-
CONFIG/AMM12/EXP00/namelist_cfg (modified) (2 diffs)
-
CONFIG/C1D_PAPA/EXP00/namelist_cfg (modified) (2 diffs)
-
CONFIG/GYRE_BFM/EXP00/namelist_cfg (modified) (2 diffs)
-
CONFIG/GYRE_PISCES/EXP00/namelist_cfg (modified) (2 diffs)
-
CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg (modified) (2 diffs)
-
CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg (modified) (3 diffs)
-
CONFIG/SHARED/field_def_nemo-opa.xml (modified) (1 diff)
-
CONFIG/SHARED/namelist_ref (modified) (5 diffs)
-
CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg (modified) (1 diff)
-
CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg (modified) (1 diff)
-
CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg (modified) (2 diffs)
-
CONFIG/cfg.txt (modified) (1 diff)
-
NEMO/OPA_SRC/DIA/dia25h.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/DYN/dynzdf.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/dynzdf_imp.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/SBC/sbcwave.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/ZDF/zdf_oce.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/ZDF/zdfddm.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/ZDF/zdfevd.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/ZDF/zdfgls.F90 (modified) (24 diffs)
-
NEMO/OPA_SRC/ZDF/zdfiwm.F90 (moved) (moved from branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90) (21 diffs)
-
NEMO/OPA_SRC/ZDF/zdfphy.F90 (modified) (8 diffs)
-
NEMO/OPA_SRC/ZDF/zdfric.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/ZDF/zdfswm.F90 (moved) (moved from branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90) (1 diff)
-
NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (19 diffs)
-
NEMO/OPA_SRC/module_example (modified) (5 diffs)
-
NEMO/OPA_SRC/step_oce.F90 (modified) (1 diff)
-
NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
r7954 r7990 324 324 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 325 325 ! 326 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)327 !328 ln_zdf qiao = .false. ! surface wave-induced mixing(T => ln_wave=ln_sdw=T )326 ! ! gravity wave-driven vertical mixing 327 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 328 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 329 329 ! 330 330 ! ! time-stepping … … 353 353 / 354 354 !----------------------------------------------------------------------- 355 &namzdf_ tmx ! internal wave-driven mixing parameterization (ln_zdftmx=T)355 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 356 356 !----------------------------------------------------------------------- 357 357 / -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
r7954 r7990 275 275 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 276 276 ! 277 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)278 !279 ln_zdf qiao = .false. ! surface wave-induced mixing(T => ln_wave=ln_sdw=T )277 ! ! gravity wave-driven vertical mixing 278 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 279 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 280 280 ! 281 281 ! ! time-stepping 282 282 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 283 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T 283 284 ! 284 285 ! ! coefficients … … 301 302 / 302 303 !----------------------------------------------------------------------- 303 &namzdf_ tmx ! internal wave-driven mixing parameterization (ln_zdftmx=T)304 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 304 305 !----------------------------------------------------------------------- 305 306 / -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
r7954 r7990 264 264 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 265 265 ! 266 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)267 !268 ln_zdf qiao = .false. ! surface wave-induced mixing(T => ln_wave=ln_sdw=T )266 ! ! gravity wave-driven vertical mixing 267 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 268 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 269 269 ! 270 270 ! ! time-stepping … … 292 292 / 293 293 !----------------------------------------------------------------------- 294 &namzdf_ tmx ! internal wave-driven mixing parameterization (ln_zdftmx=T)294 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 295 295 !----------------------------------------------------------------------- 296 296 / -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
r7954 r7990 218 218 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 219 219 ! 220 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)221 !222 ln_zdf qiao = .false. ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)220 ! ! gravity wave-driven vertical mixing 221 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 222 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 223 223 ! 224 224 ! ! time-stepping 225 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) time steppingscheme226 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T225 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 226 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T 227 227 ! 228 228 ! ! Coefficients … … 251 251 / 252 252 !----------------------------------------------------------------------- 253 &namzdf_tmx ! tidal mixing parameterization (ln_zdftmx=T)254 !-----------------------------------------------------------------------255 ln_tmx_itf = .false. ! ITF specific parameterisation256 /257 !-----------------------------------------------------------------------258 253 &nammpp ! Massively Parallel Processing ("key_mpp_mpi) 259 254 !----------------------------------------------------------------------- -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg
r7954 r7990 226 226 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 227 227 ! 228 ln_zdftmx = .true. ! tidal mixing parameterization (T => fill namzdf_tmx)229 !230 ln_zdf qiao = .false. ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)228 ! ! gravity wave-driven vertical mixing 229 ln_zdfiwm = .true. ! internal wave-induced mixing (T => fill namzdf_iwm) 230 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 231 231 ! 232 232 ! ! time-stepping 233 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) time steppingscheme234 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T233 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 234 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T 235 235 ! 236 236 ! ! Coefficients … … 245 245 / 246 246 !----------------------------------------------------------------------- 247 &namzdf_ tmx ! internal wave-driven mixing parameterization (ln_zdftmx=T)247 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 248 248 !----------------------------------------------------------------------- 249 249 nn_zpyc = 1 ! pycnocline-intensified dissipation scales as N (=1) or N^2 (=2) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg
r7954 r7990 27 27 / 28 28 !----------------------------------------------------------------------- 29 &namcrs ! Grid coarsening for dynamics output and/or 30 ! passive tracer coarsened online simulations 29 &namcrs ! coarsened grid (for outputs and/or TOP) (ln_crs =T) 31 30 !----------------------------------------------------------------------- 32 31 / … … 255 254 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 256 255 ! 257 ln_zdftmx = .true. ! tidal mixing parameterization (T => fill namzdf_tmx)258 !259 ln_zdf qiao = .false. ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)256 ! ! gravity wave-driven vertical mixing 257 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 258 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 260 259 ! 261 260 ! ! time-stepping 262 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) time steppingscheme263 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T261 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 262 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T 264 263 ! 265 264 ! ! Coefficients … … 274 273 / 275 274 !----------------------------------------------------------------------- 276 &namzdf_ tmx ! tidal mixing parameterization (ln_zdftmx=T)275 &namzdf_iwm ! tidal mixing parameterization (ln_zdfiwm =T) 277 276 !----------------------------------------------------------------------- 278 277 nn_zpyc = 2 ! pycnocline-intensified dissipation scales as N (=1) or N^2 (=2) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/field_def_nemo-opa.xml
r7954 r7990 398 398 <field id="avm_evd" long_name="convective enhancement of vertical viscosity" standard_name="ocean_vertical_momentum_diffusivity_due_to_convection" unit="m2/s" /> 399 399 400 <!-- variables available with ln_zdf tmx=T -->400 <!-- variables available with ln_zdfiwm =T --> 401 401 <field id="av_ratio" long_name="S over T diffusivity ratio" standard_name="salinity_over_temperature_diffusivity_ratio" unit="1" /> 402 402 <field id="av_wave" long_name="wave-induced vertical diffusivity" standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves" unit="m2/s" /> 403 <field id="bflx_ tmx" long_name="wave-induced buoyancy flux" standard_name="buoyancy_flux_due_to_internal_waves" unit="W/kg" />404 <field id="pcmap_ tmx" long_name="power consumed by wave-driven mixing" standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing" unit="W/m2" grid_ref="grid_W_2D" />405 <field id="emix_ tmx" long_name="power density available for mixing" standard_name="power_available_for_mixing_from_breaking_internal_waves" unit="W/kg" />403 <field id="bflx_iwm" long_name="wave-induced buoyancy flux" standard_name="buoyancy_flux_due_to_internal_waves" unit="W/kg" /> 404 <field id="pcmap_iwm" long_name="power consumed by wave-driven mixing" standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing" unit="W/m2" grid_ref="grid_W_2D" /> 405 <field id="emix_iwm" long_name="power density available for mixing" standard_name="power_available_for_mixing_from_breaking_internal_waves" unit="W/kg" /> 406 406 407 407 <!-- variables available with diaar5 --> -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
r7954 r7990 11 11 !! 6 - Tracer (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp) 12 12 !! 7 - dynamics (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 13 !! 8 - Verical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_ tmx, namzdf_tmx_new)13 !! 8 - Verical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm) 14 14 !! 9 - miscellaneous (nammpp, namctl) 15 15 !! 10 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb, namsto) … … 875 875 !! namzdf_tke TKE vertical mixing (ln_zdftke=T) 876 876 !! namzdf_gls GLS vertical mixing (ln_zdfgls=T) 877 !! namzdf_ tmx tidal mixing parameterization (ln_zdftmx=T)877 !! namzdf_iwm tidal mixing parameterization (ln_zdfiwm=T) 878 878 !!====================================================================== 879 879 ! … … 899 899 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 900 900 ! 901 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)902 !903 ln_zdf qiao = .false. ! surface wave-induced mixing(T => ln_wave=ln_sdw=T )901 ! ! gravity wave-driven vertical mixing 902 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 903 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 904 904 ! 905 905 ! ! time-stepping … … 919 919 rn_alp = 5. ! coefficient of the parameterization 920 920 nn_ric = 2 ! coefficient of the parameterization 921 rn_ekmfc = 0.7 ! Factor in the Ekman depth Equation922 rn_mldmin = 1.0 ! minimum allowable mixed-layer depth estimate (m)923 rn_mldmax = 1000.0 ! maximum allowable mixed-layer depth estimate (m)924 rn_wtmix = 10.0 ! vertical eddy viscosity coeff [m2/s] in the mixed-layer925 rn_wvmix = 10.0 ! vertical eddy diffusioncoeff [m2/s] in the mixed-layer926 ln_mldw = .true. ! Flag to use or not the mixed layer depth param.921 ln_mldw = .false. ! enhanced mixing in the Ekman layer 922 rn_ekmfc = 0.7 ! Factor in the Ekman depth Equation 923 rn_mldmin = 1.0 ! minimum allowable mixed-layer depth estimate (m) 924 rn_mldmax = 1000.0 ! maximum allowable mixed-layer depth estimate (m) 925 rn_wtmix = 10.0 ! vertical eddy viscosity coeff [m2/s] in the mixed-layer 926 rn_wvmix = 10.0 ! vertical eddy diffusion coeff [m2/s] in the mixed-layer 927 927 / 928 928 !----------------------------------------------------------------------- … … 974 974 / 975 975 !----------------------------------------------------------------------- 976 &namzdf_ tmx ! internal wave-driven mixing parameterization (ln_zdftmx=T)976 &namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T) 977 977 !----------------------------------------------------------------------- 978 978 nn_zpyc = 1 ! pycnocline-intensified dissipation scales as N (=1) or N^2 (=2) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg
r7954 r7990 344 344 ln_zdfddm = .false. ! double diffusive mixing 345 345 ! 346 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)347 !348 ln_zdf qiao = .false. ! surface wave-induced mixing(T => ln_wave=ln_sdw=T )346 ! ! gravity wave-driven vertical mixing 347 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 348 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 349 349 ! 350 350 ! ! time-stepping -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg
r7954 r7990 223 223 / 224 224 !----------------------------------------------------------------------- 225 &namzdf_ tmx ! tidal mixing parameterization (ln_zdftmx=T)225 &namzdf_iwm ! tidal mixing parameterization (ln_zdfiwm =T) 226 226 !----------------------------------------------------------------------- 227 227 / -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg
r7954 r7990 328 328 !! namzdf_tke TKE vertical mixing (ln_zdftke=T) 329 329 !! namzdf_gls GLS vertical mixing (ln_zdfgls=T) 330 !! namzdf_ tmx tidal mixing parameterization (ln_zdftmx=T)330 !! namzdf_iwm tidal mixing parameterization (ln_zdfiwm=T) 331 331 !!====================================================================== 332 332 !----------------------------------------------------------------------- … … 351 351 rn_hsbfr = 1.6 ! heat/salt buoyancy flux ratio 352 352 ! 353 ln_zdftmx = .false. ! tidal mixing parameterization (T => fill namzdf_tmx)354 !355 ln_zdf qiao = .false. ! surface wave-induced mixing (Qiao et al. 2010) (T => ln_wave=ln_sdw=T. & fill namsbc_wave)353 ! ! gravity wave-driven vertical mixing 354 ln_zdfiwm = .false. ! internal wave-induced mixing (T => fill namzdf_iwm) 355 ln_zdfswm = .false. ! surface wave-induced mixing (T => ln_wave=ln_sdw=T ) 356 356 ! 357 357 ! ! time-stepping 358 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 358 ln_zdfexp = .false. ! split-explicit (T) or implicit (F) scheme 359 nn_zdfexp= 3 ! number of sub-timestep for ln_zdfexp=T 359 360 ! 360 361 ! ! coefficients -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt
r7954 r7990 6 6 ORCA2_OFF_TRC OPA_SRC OFF_SRC TOP_SRC 7 7 ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC 8 GYRE_PISCES OPA_SRC TOP_SRC 8 9 GYRE_PISCES_XIOS OPA_SRC TOP_SRC 9 GYRE_PISCES OPA_SRC TOP_SRC -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
r7953 r7990 9 9 USE dom_oce ! ocean space and time domain 10 10 USE zdf_oce ! ocean vertical physics 11 USE zdfgls , ONLY: mxln11 USE zdfgls , ONLY : hmxn 12 12 USE in_out_manager ! I/O units 13 13 USE iom ! I/0 library … … 107 107 IF( ln_zdfgls ) THEN 108 108 en_25h(:,:,:) = en(:,:,:) 109 rmxln_25h(:,:,:) = mxln(:,:,:)109 rmxln_25h(:,:,:) = hmxn(:,:,:) 110 110 ENDIF 111 111 #if defined key_lim3 || defined key_lim2 … … 169 169 IF( ln_zdfgls ) THEN 170 170 en_25h(:,:,:) = en_25h(:,:,:) + en(:,:,:) 171 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + mxln(:,:,:)171 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxn(:,:,:) 172 172 ENDIF 173 173 cnt_25h = cnt_25h + 1 … … 217 217 zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 218 218 CALL iom_put("vomecrty25h", zw3d ) ! j-current 219 zw3d(:,:,:) = wn_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))219 zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 220 220 CALL iom_put("vomecrtz25h", zw3d ) ! k-current 221 221 ! Write vertical physics 222 zw3d(:,:,:) = avt_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))222 zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 223 223 CALL iom_put("avt25h", zw3d ) ! diffusivity 224 zw3d(:,:,:) = avm_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))224 zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 225 225 CALL iom_put("avm25h", zw3d) ! viscosity 226 226 IF( ln_zdftke ) THEN 227 zw3d(:,:,:) = en_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))227 zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 228 228 CALL iom_put("tke25h", zw3d) ! tke 229 229 ENDIF 230 230 IF( ln_zdfgls ) THEN 231 zw3d(:,:,:) = en_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))231 zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 232 232 CALL iom_put("tke25h", zw3d) ! tke 233 zw3d(:,:,:) = rmxln_25h(:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))233 zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 234 234 CALL iom_put( "mxln25h",zw3d) 235 235 ENDIF … … 249 249 IF( ln_zdfgls ) THEN 250 250 en_25h(:,:,:) = en(:,:,:) 251 rmxln_25h(:,:,:) = mxln(:,:,:)251 rmxln_25h(:,:,:) = hmxn(:,:,:) 252 252 ENDIF 253 253 cnt_25h = 1 -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r7953 r7990 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- 44 45 44 CONTAINS 46 45 … … 51 50 !! ** Purpose : compute the vertical ocean dynamics physics. 52 51 !!--------------------------------------------------------------------- 53 !!54 52 INTEGER, INTENT( in ) :: kt ! ocean time-step index 55 53 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r7753 r7990 12 12 13 13 !!---------------------------------------------------------------------- 14 !! dyn_zdf_imp : compute the vertical diffusion using a implicit scheme14 !! dyn_zdf_imp : compute the vertical diffusion using a implicit time-stepping 15 15 !! together with the Leap-Frog time integration. 16 16 !!---------------------------------------------------------------------- … … 39 39 # include "vectopt_loop_substitute.h90" 40 40 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 3.3 , NEMO Consortium (2010)41 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 42 42 !! $Id$ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r7815 r7990 9 9 !! 3.4 ! 2011_11 (C. Harris) more flexibility + multi-category fields 10 10 !!---------------------------------------------------------------------- 11 11 12 !!---------------------------------------------------------------------- 12 13 !! namsbc_cpl : coupled formulation namlist … … 974 975 !! emp upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 975 976 !!---------------------------------------------------------------------- 976 USE zdf_oce, ONLY : ln_zdf qiao977 USE zdf_oce, ONLY : ln_zdfswm 977 978 978 979 IMPLICIT NONE … … 1159 1160 ! ! Wave mean period ! 1160 1161 ! ! ========================= ! 1161 IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)1162 IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1) 1162 1163 ! 1163 1164 ! ! ========================= ! 1164 1165 ! ! Significant wave height ! 1165 1166 ! ! ========================= ! 1166 IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)1167 IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 1167 1168 ! 1168 1169 ! ! ========================= ! 1169 ! ! Vertical mixing Qiao!1170 ! ! surface wave mixing ! 1170 1171 ! ! ========================= ! 1171 IF( srcv(jpr_wnum)%laction .AND. ln_zdf qiao )wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)1172 IF( srcv(jpr_wnum)%laction .AND. ln_zdfswm ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1) 1172 1173 1173 1174 ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7864 r7990 19 19 USE oce ! ocean variables 20 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE zdf_oce, ONLY : ln_zdf qiao21 USE zdf_oce, ONLY : ln_zdfswm 22 22 USE bdy_oce ! open boundary condition variables 23 23 USE domvvl ! domain: variable volume layers … … 227 227 ! 228 228 ! Read also wave number if needed, so that it is available in coupling routines 229 IF( ln_zdf qiao.AND. .NOT.cpl_wnum ) THEN229 IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN 230 230 CALL fld_read( kt, nn_fsbc, sf_wn ) ! read wave parameters from external forcing 231 231 wnum(:,:) = sf_wn(1)%fnow(:,:,1) … … 345 345 vsd(:,:,:) = 0._wp 346 346 wsd(:,:,:) = 0._wp 347 ! Wave number needed only if ln_zdf qiao=T347 ! Wave number needed only if ln_zdfswm=T 348 348 IF( .NOT. cpl_wnum ) THEN 349 349 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
r7953 r7990 16 16 PUBLIC zdf_oce_alloc ! Called in nemogcm.F90 17 17 18 ! !!* namelist namzdf: vertical diffusion * 19 !!gm 18 ! !!* namelist namzdf: vertical physics * 20 19 ! ! vertical closure scheme flags 21 20 LOGICAL , PUBLIC :: ln_zdfcst !: constant coefficients … … 23 22 LOGICAL , PUBLIC :: ln_zdftke !: Turbulent Kinetic Energy closure 24 23 LOGICAL , PUBLIC :: ln_zdfgls !: Generic Length Sclare closure 25 ! ! tidal induced mixing 26 LOGICAL , PUBLIC :: ln_zdftmx !: tidal mixing parameterization flag 27 ! ! double diffusion 28 LOGICAL , PUBLIC :: ln_zdfddm !: double diffusive mixing flag 29 REAL(wp), PUBLIC :: rn_avts !: maximum value of avs for salt fingering 30 REAL(wp), PUBLIC :: rn_hsbfr !: heat/salt buoyancy flux ratio 31 !!gm 32 LOGICAL , PUBLIC :: ln_zdfexp !: explicit vertical diffusion scheme flag 33 INTEGER , PUBLIC :: nn_zdfexp !: number of sub-time step (explicit time stepping) 24 ! ! convection 34 25 LOGICAL , PUBLIC :: ln_zdfevd !: convection: enhanced vertical diffusion flag 35 26 INTEGER , PUBLIC :: nn_evdm !: =0/1 flag to apply enhanced avm or not … … 38 29 INTEGER , PUBLIC :: nn_npc !: non penetrative convective scheme call frequency 39 30 INTEGER , PUBLIC :: nn_npcp !: non penetrative convective scheme print frequency 40 ! ! Surface wave-induced mixing 41 LOGICAL , PUBLIC :: ln_zdfqiao !: Enhanced wave vertical mixing Qiao(2010) formulation flag 31 ! ! double diffusion 32 LOGICAL , PUBLIC :: ln_zdfddm !: double diffusive mixing flag 33 REAL(wp), PUBLIC :: rn_avts !: maximum value of avs for salt fingering 34 REAL(wp), PUBLIC :: rn_hsbfr !: heat/salt buoyancy flux ratio 35 ! ! gravity wave-induced vertical mixing 36 LOGICAL , PUBLIC :: ln_zdfswm !: surface wave-induced mixing flag 37 LOGICAL , PUBLIC :: ln_zdfiwm !: internal wave-induced mixing flag 38 ! ! time-stepping 39 LOGICAL , PUBLIC :: ln_zdfexp !: explicit vertical diffusion scheme flag 40 INTEGER , PUBLIC :: nn_zdfexp !: number of sub-time step (explicit time stepping) 42 41 ! ! coefficients 43 42 REAL(wp), PUBLIC :: rn_avm0 !: vertical eddy viscosity (m2/s) … … 47 46 48 47 49 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avt , avs !: vertical eddy diffusivity coef at w-pt [m2/s] 50 48 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avt , avs !: vertical T & S diffusivity coef at w-pt [m2/s] 49 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm !: vertical momentum viscosity coef at w-pt [m2/s] 50 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k !: avt, avs computed by turbulent closure alone 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 51 53 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: avmb , avtb !: background profile of avm and avt 52 54 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 53 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: bfrua, bfrva !: Bottom friction coefficients set in zdfbfr 54 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: tfrua, tfrva !: top friction coefficients set in zdfbfr 55 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 56 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm !: vertical viscosity & diffusivity coef at w-pt [m2/s] 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k ! not enhanced Kz 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmu_k, avmv_k ! not enhanced Kz 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 55 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: bfrua, bfrva !: bottom friction coefficients 56 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: tfrua, tfrva !: top friction coefficients 60 57 61 58 !!---------------------------------------------------------------------- … … 71 68 !!---------------------------------------------------------------------- 72 69 ! 73 ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) , & 74 & avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) , & 75 & tfrua(jpi, jpj), tfrva(jpi, jpj) , & 76 & avmu (jpi,jpj,jpk), avm (jpi,jpj,jpk) , & 77 & avmv (jpi,jpj,jpk), avt (jpi,jpj,jpk) , avs (jpi,jpj,jpk), & 78 & avt_k (jpi,jpj,jpk), avm_k (jpi,jpj,jpk) , & 79 & avmu_k(jpi,jpj,jpk), avmv_k(jpi,jpj,jpk) , & 80 & en (jpi,jpj,jpk), STAT = zdf_oce_alloc ) 70 ALLOCATE( avm (jpi,jpj,jpk) , avt (jpi,jpj,jpk) , avs(jpi,jpj,jpk) , & 71 & avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) , & 72 & avmb(jpk) , bfrua(jpi,jpj) , tfrua(jpi, jpj) , & 73 & avtb(jpk) , bfrva(jpi,jpj) , tfrva(jpi, jpj) , avtb_2d(jpi,jpj) , & 74 & avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk) , STAT = zdf_oce_alloc ) 81 75 ! 82 76 IF( zdf_oce_alloc /= 0 ) CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays') -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r7931 r7990 8 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 9 9 !! 3.6 ! 2013-04 (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 10 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 10 11 !!---------------------------------------------------------------------- 11 12 … … 22 23 USE prtctl ! Print control 23 24 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays25 25 USE timing ! Timing 26 26 … … 63 63 !! avt = avt + zavft + zavdt 64 64 !! avs = avs + zavfs + zavds 65 !! avm u, avmv arerequired to remain at least above avt and avs.65 !! avm is required to remain at least above avt and avs. 66 66 !! 67 67 !! ** Action : avt, avs : updated vertical eddy diffusivity coef. for T & S … … 77 77 REAL(wp) :: zavft, zavfs ! - - 78 78 REAL(wp) :: zavdt, zavds ! - - 79 REAL(wp), POINTER, DIMENSION(:,:) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd379 REAL(wp), DIMENSION(jpi,jpj) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 80 80 !!---------------------------------------------------------------------- 81 81 ! 82 IF( nn_timing == 1 ) CALL timing_start('zdf_ddm') 83 ! 84 CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 82 IF( nn_timing == 1 ) CALL timing_start('zdf_ddm') 85 83 ! 86 84 ! ! =============== … … 89 87 ! Define the mask 90 88 ! --------------- 89 !!gm WORK to be done: change the code from vector optimisation to scalar one. 90 !!gm ==>>> test in the loop instead of use of mask arrays 91 !!gm and many acces in memory 92 91 93 DO jj = 1, jpj ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 92 94 DO ji = 1, jpi 93 95 zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) & 96 !!gm please, use e3w_n below 94 97 & / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 95 98 ! … … 156 159 END DO 157 160 END DO 158 159 160 ! Increase avmu, avmv if necessary161 ! --------------------------------162 !!gm to be changed following the definition of avm.163 DO jj = 1, jpjm1164 DO ji = 1, fs_jpim1 ! vector opt.165 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), &166 & avt(ji,jj,jk), avt(ji+1,jj,jk), &167 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * wumask(ji,jj,jk)168 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), &169 & avt(ji,jj,jk), avt(ji,jj+1,jk), &170 & avs(ji,jj,jk), avs(ji,jj+1,jk) ) * wvmask(ji,jj,jk)171 END DO172 END DO173 161 ! ! =============== 174 162 END DO ! End of slab 175 163 ! ! =============== 176 164 ! 177 CALL lbc_lnk( avt , 'W', 1._wp ) ! Lateral boundary conditions (unchanged sign)178 CALL lbc_lnk( avs , 'W', 1._wp )179 CALL lbc_lnk( avm , 'W', 1._wp )180 CALL lbc_lnk( avmu, 'U', 1._wp )181 CALL lbc_lnk( avmv, 'V', 1._wp )182 183 165 IF(ln_ctl) THEN 184 166 CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk) 185 CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm - u: ', mask1=umask, &186 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk)187 167 ENDIF 188 !189 CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )190 168 ! 191 169 IF( nn_timing == 1 ) CALL timing_stop('zdf_ddm') -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r7931 r7990 8 8 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 9 9 !! 3.2 ! 2009-03 (M. Leclair, G. Madec, R. Benshila) test on both before & after 10 !! 4.0 ! 2017-04 (G. Madec) evd applied on avm (at t-point) 10 11 !!---------------------------------------------------------------------- 11 12 … … 23 24 USE iom ! for iom_put 24 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 USE wrk_nemo ! work arrays26 26 USE timing ! Timing 27 27 … … 45 45 !! sivity coefficients when a static instability is encountered. 46 46 !! 47 !! ** Method : avt, avm, and the 4 neighbouring avmu, avmv coefficients 48 !! are set to avevd (namelist parameter) if the water column is 49 !! statically unstable (i.e. if rn2 < -1.e-12 ) 47 !! ** Method : tracer (and momentum if nn_evdm=1) vertical mixing 48 !! coefficients are set to rn_evd (namelist parameter) 49 !! if the water column is statically unstable. 50 !! The test of static instability is performed using 51 !! Brunt-Vaisala frequency (rn2 < -1.e-12) of to successive 52 !! time-step (Leap-Frog environnement): before and 53 !! now time-step. 50 54 !! 51 !! ** Action : avt, avm, avmu, avmv updted in static instability cases 52 !! 53 !! References : Lazar, A., these de l'universite Paris VI, France, 1997 55 !! ** Action : avt, avm enhanced where static instability occurs 54 56 !!---------------------------------------------------------------------- 55 INTEGER, INTENT( in) :: kt ! ocean time-step indexocean time step57 INTEGER, INTENT(in) :: kt ! ocean time-step indexocean time step 56 58 ! 57 59 INTEGER :: ji, jj, jk ! dummy loop indices 58 REAL(wp), POINTER, DIMENSION(:,:,:) :: zavt_evd, zavm_evd60 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zavt_evd, zavm_evd 59 61 !!---------------------------------------------------------------------- 60 62 ! … … 68 70 ENDIF 69 71 ! 70 CALL wrk_alloc( jpi,jpj,jpk, zavt_evd, zavm_evd )71 72 ! 72 73 zavt_evd(:,:,:) = avt(:,:,:) ! set avt prior to evd application … … 74 75 SELECT CASE ( nn_evdm ) 75 76 ! 76 CASE ( 1 ) ! enhance vertical eddy viscosity and diffusivity(if rn2<-1.e-12)77 CASE ( 1 ) !== enhance tracer & momentum Kz ==! (if rn2<-1.e-12) 77 78 ! 78 79 zavm_evd(:,:,:) = avm(:,:,:) ! set avm prior to evd application 79 80 ! 81 !! change last digits results 82 ! WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) ) <= -1.e-12 ) THEN 83 ! avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 84 ! avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 85 ! END WHERE 86 ! 80 87 DO jk = 1, jpkm1 81 DO jj = 2, jpj ! no vector opt.82 DO ji = 2, jpi 88 DO jj = 2, jpjm1 89 DO ji = 2, jpim1 83 90 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 84 avt (ji ,jj ,jk) = rn_evd * tmask(ji ,jj ,jk) 85 avm (ji ,jj ,jk) = rn_evd * tmask(ji ,jj ,jk) 86 avmu(ji ,jj ,jk) = rn_evd * umask(ji ,jj ,jk) 87 avmu(ji-1,jj ,jk) = rn_evd * umask(ji-1,jj ,jk) 88 avmv(ji ,jj ,jk) = rn_evd * vmask(ji ,jj ,jk) 89 avmv(ji ,jj-1,jk) = rn_evd * vmask(ji ,jj-1,jk) 91 avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 92 avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 90 93 ENDIF 91 94 END DO 92 95 END DO 93 96 END DO 94 CALL lbc_lnk( avt , 'W', 1. ) ; CALL lbc_lnk( avm , 'W', 1. ) ! Lateral boundary conditions95 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )96 97 ! 97 98 zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:) ! change in avm due to evd 98 99 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 99 100 ! 100 CASE DEFAULT ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 101 CASE DEFAULT !== enhance tracer Kz ==! (if rn2<-1.e-12) 102 !! change last digits results 103 ! WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) ) <= -1.e-12 ) 104 ! avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 105 ! END WHERE 106 101 107 DO jk = 1, jpkm1 102 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 103 DO jj = 1, jpj ! loop over the whole domain (no lbc_lnk call) 104 DO ji = 1, jpi 108 DO jj = 2, jpjm1 109 DO ji = 2, jpim1 105 110 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 106 avt(ji,jj,jk) = rn_evd * tmask(ji,jj,jk)111 avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 107 112 END DO 108 113 END DO … … 110 115 ! 111 116 END SELECT 112 117 ! 113 118 zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd 114 119 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 115 120 IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 116 !117 CALL wrk_dealloc( jpi,jpj,jpk, zavt_evd, zavm_evd )118 121 ! 119 122 IF( nn_timing == 1 ) CALL timing_stop('zdf_evd') -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7953 r7990 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 9 10 !!---------------------------------------------------------------------- 10 11 … … 36 37 PRIVATE 37 38 38 PUBLIC zdf_gls ! routine called in step module39 PUBLIC zdf_gls_init ! routine called in zdfphy module40 PUBLIC gls_rst ! routine called in zdfphy module39 PUBLIC zdf_gls ! called in zdfphy 40 PUBLIC zdf_gls_init ! called in zdfphy 41 PUBLIC gls_rst ! called in zdfphy 41 42 42 43 ! 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxn !: now mixing length 44 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 45 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points … … 113 114 !! *** FUNCTION zdf_gls_alloc *** 114 115 !!---------------------------------------------------------------------- 115 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , &116 ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 116 117 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 117 118 ! … … 148 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 149 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3du, z3dv ! u- and v-shears 150 152 !!-------------------------------------------------------------------- 151 153 ! 152 IF( nn_timing == 1 ) CALL timing_start('zdf_gls')154 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 153 155 ! 154 156 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 155 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 156 157 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 158 CALL wrk_alloc( jpi,jpj,jpk, z3du, z3dv ) 159 157 160 ! Preliminary computing 158 161 159 162 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 160 163 161 IF( kt /= nit000 ) THEN ! restore before value to compute tke 162 avt (:,:,:) = avt_k (:,:,:) 163 avm (:,:,:) = avm_k (:,:,:) 164 avmu(:,:,:) = avmu_k(:,:,:) 165 avmv(:,:,:) = avmv_k(:,:,:) 166 ENDIF 164 ! restore before values computed GLS alone 165 avt(:,:,:) = avt_k(:,:,:) 166 avm(:,:,:) = avm_k(:,:,:) 167 167 168 168 ! Compute surface and bottom friction at T-points … … 183 183 END DO 184 184 185 ! Set surface roughness length 186 SELECT CASE ( nn_z0_met ) 187 ! 188 CASE ( 0 ) ! Constant roughness 185 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 186 CASE ( 0 ) ! Constant roughness 189 187 zhsro(:,:) = rn_hsro 190 188 CASE ( 1 ) ! Standard Charnock formula 191 189 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 192 190 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 191 !!gm zcof = 2._wp * 0.6_wp / 28._wp 192 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustars2(:,:),rsmall) ) ) ! Wave age (eq. 10) 193 193 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 194 194 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 195 195 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 196 !!gm BUG missing a multiplicative coefficient.... 196 197 zhsro(:,:) = hsw(:,:) 197 198 END SELECT 198 199 199 ! Compute shear and dissipation rate200 DO jk = 2, jpkm1201 DO jj = 2, jpjm1202 DO ji = fs_2, fs_jpim1 ! vector opt.203 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) &204 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) &205 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )206 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) &207 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) &208 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )209 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)200 DO jk = 2, jpkm1 !== Compute shear and dissipation rate ==! 201 DO jj = 1, jpjm1 202 DO ji = 1, jpim1 203 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 204 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 205 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 206 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 207 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 208 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 209 ! 210 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / hmxn(ji,jj,jk) 210 211 END DO 211 212 END DO 212 213 END DO 213 !214 ! Lateral boundary conditions (avmu,avmv) (sign unchanged)215 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )216 214 217 215 ! Save tke at before time step 218 216 eb (:,:,:) = en (:,:,:) 219 mxlb(:,:,:) = mxln(:,:,:)217 mxlb(:,:,:) = hmxn(:,:,:) 220 218 221 219 IF( nn_clos == 0 ) THEN ! Mellor-Yamada … … 223 221 DO jj = 2, jpjm1 224 222 DO ji = fs_2, fs_jpim1 ! vector opt. 225 zup = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)223 zup = hmxn(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 226 224 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 227 225 zcoef = ( zup / MAX( zdown, rsmall ) ) … … 247 245 DO jk = 2, jpkm1 248 246 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt.247 DO ji = 2, jpim1 250 248 ! 251 249 ! shear prod. at w-point weightened by mask 252 shear(ji,jj,jk) = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &253 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )250 shear(ji,jj,jk) = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 251 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 254 252 ! 255 253 ! stratif. destruction … … 278 276 ! building the matrix 279 277 zcof = rfact_tke * tmask(ji,jj,jk) 280 ! 281 ! lower diagonal 282 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 283 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 284 ! 285 ! upper diagonal 286 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 287 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 288 ! 289 ! diagonal 290 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 291 & + rdt * zdiss * tmask(ji,jj,jk) 292 ! 293 ! right hand side in en 294 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 278 ! ! lower diagonal 279 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 280 ! ! upper diagonal 281 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 282 ! ! diagonal 283 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 284 ! ! right hand side in en 285 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 295 286 END DO 296 287 END DO … … 300 291 ! 301 292 ! Set surface condition on zwall_psi (1 at the bottom) 302 zwall_psi(:,:, 1) = zwall_psi(:,:,2)303 zwall_psi(:,:,jpk) = 1. 293 zwall_psi(:,:, 1 ) = zwall_psi(:,:,2) 294 zwall_psi(:,:,jpk) = 1._wp 304 295 ! 305 296 ! Surface boundary condition on tke … … 353 344 ! 354 345 CASE ( 0 ) ! Dirichlet 355 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin346 ! ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin 356 347 ! ! Balance between the production and the dissipation terms 357 348 DO jj = 2, jpjm1 … … 503 494 ! building the matrix 504 495 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 505 ! lower diagonal 506 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 507 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 508 ! upper diagonal 509 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 510 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 511 ! diagonal 512 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 513 & + rdt * zdiss * tmask(ji,jj,jk) 514 ! 515 ! right hand side in psi 516 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 496 ! ! lower diagonal 497 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 498 ! ! upper diagonal 499 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 500 ! ! diagonal 501 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 502 ! ! right hand side in psi 503 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 517 504 END DO 518 505 END DO … … 565 552 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 566 553 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 567 568 554 ! 569 555 ! … … 577 563 ! 578 564 CASE ( 0 ) ! Dirichlet 579 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0565 ! ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0 580 566 ! ! Balance between the production and the dissipation terms 581 567 DO jj = 2, jpjm1 … … 700 686 ! Limit dissipation rate under stable stratification 701 687 ! -------------------------------------------------- 702 DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time688 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time 703 689 DO jj = 2, jpjm1 704 690 DO ji = fs_2, fs_jpim1 ! vector opt. 705 691 ! limitation 706 692 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 707 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)693 hmxn(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 708 694 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 709 695 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 710 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )696 IF (ln_length_lim) hmxn(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) ) 711 697 END DO 712 698 END DO … … 733 719 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 734 720 ! 735 ! Store stability function in avmu and avmv736 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)737 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)721 ! Store stability function in z3du and z3dv 722 z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 723 z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 738 724 END DO 739 725 END DO … … 761 747 sh = (rs4 - rs5*gh + rs6*gm) / rcff 762 748 ! 763 ! Store stability function in avmu and avmv764 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)765 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)749 ! Store stability function in z3du and z3dv 750 z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 751 z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 766 752 END DO 767 753 END DO … … 773 759 ! Lines below are useless if GOTM style Dirichlet conditions are used 774 760 775 avmv(:,:,1) = avmv(:,:,2)761 z3dv(:,:,1) = z3dv(:,:,2) 776 762 777 763 DO jj = 2, jpjm1 778 764 DO ji = fs_2, fs_jpim1 ! vector opt. 779 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))765 z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj)) 780 766 END DO 781 767 END DO … … 786 772 DO jj = 2, jpjm1 787 773 DO ji = fs_2, fs_jpim1 ! vector opt. 788 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)789 zav = zsqen * avmu(ji,jj,jk)790 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) *tmask(ji,jj,jk) ! apply mask for zdfmxl routine791 zav = zsqen * avmv(ji,jj,jk)792 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom774 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxn(ji,jj,jk) 775 zav = zsqen * z3du(ji,jj,jk) 776 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 777 zav = zsqen * z3dv(ji,jj,jk) 778 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 793 779 END DO 794 780 END DO 795 781 END DO 796 !797 ! Lateral boundary conditions (sign unchanged)798 782 avt(:,:,1) = 0._wp 783 !!gm I'm sure this is needed to compute the shear term at next time-step 799 784 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 800 801 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 802 DO jj = 2, jpjm1 803 DO ji = fs_2, fs_jpim1 ! vector opt. 804 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 805 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 806 END DO 807 END DO 808 END DO 809 avmu(:,:,1) = 0._wp ; avmv(:,:,1) = 0._wp ! set surface to zero 810 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 811 785 ! 812 786 IF(ln_ctl) THEN 813 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 814 CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls - u: ', mask1=umask, & 815 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 787 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 788 CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk ) 816 789 ENDIF 817 790 ! 818 avt_k (:,:,:) = avt (:,:,:)819 avm_k (:,:,:) = avm(:,:,:)820 avmu_k(:,:,:) = avmu(:,:,:)821 avmv_k(:,:,:) = avmv(:,:,:)791 avt_k(:,:,:) = avt(:,:,:) !== store avt, avm values computed by GLS only ==! 792 avm_k(:,:,:) = avm(:,:,:) 793 ! 794 IF( lrst_oce ) CALL gls_rst( kt, 'WRITE' ) ! write the GLS info in the restart file 822 795 ! 823 796 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 824 797 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 825 !826 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls')827 !798 CALL wrk_dealloc( jpi,jpj,jpk, z3du, z3dv ) 799 ! 800 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') 828 801 ! 829 802 END SUBROUTINE zdf_gls … … 835 808 !! 836 809 !! ** Purpose : Initialization of the vertical eddy diffivity and 837 !! viscosity when using a glsturbulent closure scheme810 !! viscosity computed using a GLS turbulent closure scheme 838 811 !! 839 812 !! ** Method : Read the namzdf_gls namelist and check the parameters 840 !! called at the first timestep (nit000)841 813 !! 842 814 !! ** input : Namlist namzdf_gls … … 892 864 ENDIF 893 865 894 ! !* allocate glsarrays866 ! !* allocate GLS arrays 895 867 IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' ) 896 868 … … 1120 1092 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1121 1093 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1122 1094 ! 1123 1095 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1124 1096 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1125 1097 ! 1126 1098 ! !* Wall proximity function 1127 1099 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1128 1100 1129 ! !* set vertical eddy coef. to the background value 1130 DO jk = 1, jpk 1131 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 1132 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 1133 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 1134 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1135 END DO 1136 ! 1137 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1101 ! !* read or initialize all required files 1102 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxn) 1138 1103 ! 1139 1104 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls_init') … … 1152 1117 !! set to rn_emin or recomputed (nn_igls/=0) 1153 1118 !!---------------------------------------------------------------------- 1154 INTEGER , INTENT(in) :: kt 1155 CHARACTER(len=*), INTENT(in) :: cdrw 1119 INTEGER , INTENT(in) :: kt ! ocean time-step 1120 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1156 1121 ! 1157 1122 INTEGER :: jit, jk ! dummy loop indices 1158 INTEGER :: id1, id2, id3, id4 , id5, id61123 INTEGER :: id1, id2, id3, id4 1159 1124 INTEGER :: ji, jj, ikbu, ikbv 1160 1125 REAL(wp):: cbx, cby … … 1165 1130 IF( ln_rstart ) THEN !* Read the restart file 1166 1131 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 1167 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 1168 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 1169 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 1170 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 1171 id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. ) 1132 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 1133 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 1134 id4 = iom_varid( numror, 'hmxn' , ldstop = .FALSE. ) 1172 1135 ! 1173 IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist 1174 CALL iom_get( numror, jpdom_autoglo, 'en' , en ) 1175 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 1176 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 1177 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 1178 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 1179 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1136 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 1137 CALL iom_get( numror, jpdom_autoglo, 'en' , en ) 1138 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 1139 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 1140 CALL iom_get( numror, jpdom_autoglo, 'hmxn' , hmxn ) 1180 1141 ELSE 1181 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1182 en (:,:,:) = rn_emin 1183 mxln(:,:,:) = 0.05 1184 avt_k (:,:,:) = avt (:,:,:) 1185 avm_k (:,:,:) = avm (:,:,:) 1186 avmu_k(:,:,:) = avmu(:,:,:) 1187 avmv_k(:,:,:) = avmv(:,:,:) 1142 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop' 1143 en (:,:,:) = rn_emin 1144 hmxn (:,:,:) = 0.05_wp 1145 avt_k(:,:,:) = avt(:,:,:) 1146 avm_k(:,:,:) = avm(:,:,:) 1188 1147 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1189 1148 ENDIF 1190 1149 ELSE !* Start from rest 1191 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'1150 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values' 1192 1151 en (:,:,:) = rn_emin 1193 mxln(:,:,:) = 0.05 1152 hmxn(:,:,:) = 0.05_wp 1153 DO jk = 1, jpk 1154 avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 1155 avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk) 1156 END DO 1194 1157 ENDIF 1195 1158 ! … … 1197 1160 ! ! ------------------- 1198 1161 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1199 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1200 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1201 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1202 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1203 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1204 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) 1162 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1163 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 1164 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 1165 CALL iom_rstput( kt, nitrst, numrow, 'hmxn' , hmxn ) 1205 1166 ! 1206 1167 ENDIF -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
r7967 r7990 1 MODULE zdf tmx1 MODULE zdfiwm 2 2 !!======================================================================== 3 !! *** MODULE zdf tmx***3 !! *** MODULE zdfiwm *** 4 4 !! Ocean physics: Internal gravity wave-driven vertical mixing 5 5 !!======================================================================== … … 8 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 9 9 !! 3.6 ! 2016-03 (C. de Lavergne) New param: internal wave-driven mixing 10 !! 4.0 ! 2017-04 (G. Madec) Remove the old tidal mixing param. and key zdftmx(_new)10 !! 4.0 ! 2017-04 (G. Madec) renamed module, remove the old param. and the CPP keys 11 11 !!---------------------------------------------------------------------- 12 12 13 13 !!---------------------------------------------------------------------- 14 !! zdf_ tmx: global momentum & tracer Kz with wave induced Kz15 !! zdf_ tmx_init : global momentum & tracer Kz with wave induced Kz14 !! zdf_iwm : global momentum & tracer Kz with wave induced Kz 15 !! zdf_iwm_init : global momentum & tracer Kz with wave induced Kz 16 16 !!---------------------------------------------------------------------- 17 17 USE oce ! ocean dynamics and tracers variables … … 33 33 PRIVATE 34 34 35 PUBLIC zdf_ tmx! called in step module36 PUBLIC zdf_ tmx_init ! called in nemogcm module37 PUBLIC zdf_ tmx_alloc ! called in nemogcm module38 39 ! !!* Namelist namzdf_ tmx: internal wave-driven mixing *35 PUBLIC zdf_iwm ! called in step module 36 PUBLIC zdf_iwm_init ! called in nemogcm module 37 PUBLIC zdf_iwm_alloc ! called in nemogcm module 38 39 ! !!* Namelist namzdf_iwm : internal wave-driven mixing * 40 40 INTEGER :: nn_zpyc ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 41 41 LOGICAL :: ln_mevar ! variable (=T) or constant (=F) mixing efficiency … … 44 44 REAL(wp) :: r1_6 = 1._wp / 6._wp 45 45 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_ tmx! power available from high-mode wave breaking (W/m2)47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: epyc_ tmx! power available from low-mode, pycnocline-intensified wave breaking (W/m2)48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_ tmx! power available from low-mode, critical slope wave breaking (W/m2)49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_ tmx! WKB decay scale for high-mode energy dissipation (m)50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_ tmx! decay scale for low-mode critical slope dissipation (m)51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emix_ tmx! local energy density available for mixing (W/kg)52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bflx_ tmx! buoyancy flux Kz * N^2 (W/kg)53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pcmap_ tmx! vertically integrated buoyancy flux (W/m2)46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_iwm ! power available from high-mode wave breaking (W/m2) 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: epyc_iwm ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_iwm ! power available from low-mode, critical slope wave breaking (W/m2) 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_iwm ! WKB decay scale for high-mode energy dissipation (m) 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_iwm ! decay scale for low-mode critical slope dissipation (m) 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emix_iwm ! local energy density available for mixing (W/kg) 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bflx_iwm ! buoyancy flux Kz * N^2 (W/kg) 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pcmap_iwm ! vertically integrated buoyancy flux (W/m2) 54 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 55 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_wave ! Internal wave-induced diffusivity … … 64 64 CONTAINS 65 65 66 INTEGER FUNCTION zdf_ tmx_alloc()67 !!---------------------------------------------------------------------- 68 !! *** FUNCTION zdf_ tmx_alloc ***69 !!---------------------------------------------------------------------- 70 ALLOCATE( ebot_ tmx(jpi,jpj), epyc_tmx(jpi,jpj), ecri_tmx(jpi,jpj) , &71 & hbot_ tmx(jpi,jpj), hcri_tmx(jpi,jpj), emix_tmx(jpi,jpj,jpk), &72 & bflx_ tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk), &73 & zav_wave(jpi,jpj,jpk), STAT=zdf_ tmx_alloc )74 ! 75 IF( lk_mpp ) CALL mpp_sum ( zdf_ tmx_alloc )76 IF( zdf_ tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')77 END FUNCTION zdf_ tmx_alloc78 79 80 SUBROUTINE zdf_ tmx( kt )81 !!---------------------------------------------------------------------- 82 !! *** ROUTINE zdf_ tmx***66 INTEGER FUNCTION zdf_iwm_alloc() 67 !!---------------------------------------------------------------------- 68 !! *** FUNCTION zdf_iwm_alloc *** 69 !!---------------------------------------------------------------------- 70 ALLOCATE( ebot_iwm(jpi,jpj), epyc_iwm(jpi,jpj), ecri_iwm(jpi,jpj) , & 71 & hbot_iwm(jpi,jpj), hcri_iwm(jpi,jpj), emix_iwm(jpi,jpj,jpk), & 72 & bflx_iwm(jpi,jpj,jpk), pcmap_iwm(jpi,jpj), zav_ratio(jpi,jpj,jpk), & 73 & zav_wave(jpi,jpj,jpk), STAT=zdf_iwm_alloc ) 74 ! 75 IF( lk_mpp ) CALL mpp_sum ( zdf_iwm_alloc ) 76 IF( zdf_iwm_alloc /= 0 ) CALL ctl_warn('zdf_iwm_alloc: failed to allocate arrays') 77 END FUNCTION zdf_iwm_alloc 78 79 80 SUBROUTINE zdf_iwm( kt ) 81 !!---------------------------------------------------------------------- 82 !! *** ROUTINE zdf_iwm *** 83 83 !! 84 84 !! ** Purpose : add to the vertical mixing coefficients the effect of … … 86 86 !! 87 87 !! ** Method : - internal wave-driven vertical mixing is given by: 88 !! Kz_wave = min( 100 cm2/s, f( Reb = emix_ tmx/( Nu * N^2 ) )89 !! where emix_ tmxis the 3D space distribution of the wave-breaking88 !! Kz_wave = min( 100 cm2/s, f( Reb = emix_iwm /( Nu * N^2 ) ) 89 !! where emix_iwm is the 3D space distribution of the wave-breaking 90 90 !! energy and Nu the molecular kinematic viscosity. 91 91 !! The function f(Reb) is linear (constant mixing efficiency) 92 92 !! if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 93 93 !! 94 !! - Compute emix_ tmx, the 3D power density that allows to compute94 !! - Compute emix_iwm, the 3D power density that allows to compute 95 95 !! Reb and therefrom the wave-induced vertical diffusivity. 96 96 !! This is divided into three components: 97 97 !! 1. Bottom-intensified low-mode dissipation at critical slopes 98 !! emix_ tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx)99 !! / ( 1. - EXP( - H/hcri_ tmx ) ) * hcri_tmx100 !! where hcri_ tmxis the characteristic length scale of the bottom101 !! intensification, ecri_ tmxa map of available power, and H the ocean depth.98 !! emix_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 99 !! / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 100 !! where hcri_iwm is the characteristic length scale of the bottom 101 !! intensification, ecri_iwm a map of available power, and H the ocean depth. 102 102 !! 2. Pycnocline-intensified low-mode dissipation 103 !! emix_ tmx(z) = ( epyc_tmx/ rau0 ) * ( sqrt(rn2(z))^nn_zpyc )103 !! emix_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 104 104 !! / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 105 !! where epyc_ tmxis a map of available power, and nn_zpyc105 !! where epyc_iwm is a map of available power, and nn_zpyc 106 106 !! is the chosen stratification-dependence of the internal wave 107 107 !! energy dissipation. 108 108 !! 3. WKB-height dependent high mode dissipation 109 !! emix_ tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx)110 !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_ tmx) * e3w(z) )111 !! where hbot_ tmxis the characteristic length scale of the WKB bottom112 !! intensification, ebot_ tmxis a map of available power, and z_wkb is the109 !! emix_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 110 !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 111 !! where hbot_iwm is the characteristic length scale of the WKB bottom 112 !! intensification, ebot_iwm is a map of available power, and z_wkb is the 113 113 !! WKB-stretched height above bottom defined as 114 114 !! z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) … … 118 118 !! avt = avt + av_wave 119 119 !! avm = avm + av_wave 120 !! avmu = avmu + mi(av_wave)121 !! avmv = avmv + mj(av_wave)122 120 !! 123 121 !! - if namelist parameter ln_tsdiff = T, account for differential mixing: 124 122 !! avs = avt + av_wave * diffusivity_ratio(Reb) 125 123 !! 126 !! ** Action : - Define emix_ tmxused to compute internal wave-induced mixing127 !! - avt, avs, avm, avmu, avmvincreased by internal wave-driven mixing124 !! ** Action : - Define emix_iwm used to compute internal wave-induced mixing 125 !! - avt, avs, avm, increased by internal wave-driven mixing 128 126 !! 129 127 !! References : de Lavergne et al. 2015, JPO; 2016, in prep. … … 142 140 !!---------------------------------------------------------------------- 143 141 ! 144 IF( nn_timing == 1 ) CALL timing_start('zdf_ tmx')142 IF( nn_timing == 1 ) CALL timing_start('zdf_iwm') 145 143 ! 146 144 CALL wrk_alloc( jpi,jpj, zfact, zhdep ) … … 156 154 DO ji = 1, jpi 157 155 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 158 zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_ tmx(ji,jj) ) )159 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj)156 zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) ) 157 IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 160 158 END DO 161 159 END DO 162 160 163 161 DO jk = 2, jpkm1 ! complete with the level-dependent part 164 emix_ tmx(:,:,jk) = zfact(:,:) * ( EXP( ( gde3w_n(:,:,jk ) - zhdep(:,:) ) / hcri_tmx(:,:) ) &165 & - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_ tmx(:,:) ) ) * wmask(:,:,jk) &162 emix_iwm(:,:,jk) = zfact(:,:) * ( EXP( ( gde3w_n(:,:,jk ) - zhdep(:,:) ) / hcri_iwm(:,:) ) & 163 & - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) ) ) * wmask(:,:,jk) & 166 164 !!gm delta(gde3w_n) = e3t_n !! Please verify the grid-point position w versus t-point 167 165 & / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) … … 170 168 ! !* Pycnocline-intensified mixing: distribute energy over the time-varying 171 169 ! !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 172 170 ! 173 171 SELECT CASE ( nn_zpyc ) 174 172 ! 175 173 CASE ( 1 ) ! Dissipation scales as N (recommended) 176 174 ! 177 175 zfact(:,:) = 0._wp 178 176 DO jk = 2, jpkm1 ! part independent of the level 179 177 zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 180 178 END DO 181 179 ! 182 180 DO jj = 1, jpj 183 181 DO ji = 1, jpi 184 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_ tmx(ji,jj) / ( rau0 * zfact(ji,jj) )185 END DO 186 END DO 187 182 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 183 END DO 184 END DO 185 ! 188 186 DO jk = 2, jpkm1 ! complete with the level-dependent part 189 emix_ tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk)190 END DO 191 187 emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 188 END DO 189 ! 192 190 CASE ( 2 ) ! Dissipation scales as N^2 193 191 ! 194 192 zfact(:,:) = 0._wp 195 193 DO jk = 2, jpkm1 ! part independent of the level 196 194 zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 197 195 END DO 198 196 ! 199 197 DO jj= 1, jpj 200 198 DO ji = 1, jpi 201 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_ tmx(ji,jj) / ( rau0 * zfact(ji,jj) )202 END DO 203 END DO 204 199 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 200 END DO 201 END DO 202 ! 205 203 DO jk = 2, jpkm1 ! complete with the level-dependent part 206 emix_ tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)207 END DO 208 204 emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 205 END DO 206 ! 209 207 END SELECT 210 208 211 209 ! !* WKB-height dependent mixing: distribute energy over the time-varying 212 210 ! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 213 211 ! 214 212 zwkb(:,:,:) = 0._wp 215 213 zfact(:,:) = 0._wp … … 218 216 zwkb(:,:,jk) = zfact(:,:) 219 217 END DO 220 218 ! 221 219 DO jk = 2, jpkm1 222 220 DO jj = 1, jpj 223 221 DO ji = 1, jpi 224 222 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 225 &* tmask(ji,jj,jk) / zfact(ji,jj)223 & * tmask(ji,jj,jk) / zfact(ji,jj) 226 224 END DO 227 225 END DO 228 226 END DO 229 227 zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 230 228 ! 231 229 zweight(:,:,:) = 0._wp 232 230 DO jk = 2, jpkm1 233 zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_ tmx(:,:) * wmask(:,:,jk) &234 & * ( EXP( -zwkb(:,:,jk) / hbot_ tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) ) )235 END DO 236 231 zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk) & 232 & * ( EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) ) ) 233 END DO 234 ! 237 235 zfact(:,:) = 0._wp 238 236 DO jk = 2, jpkm1 ! part independent of the level 239 237 zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 240 238 END DO 241 239 ! 242 240 DO jj = 1, jpj 243 241 DO ji = 1, jpi 244 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_ tmx(ji,jj) / ( rau0 * zfact(ji,jj) )245 END DO 246 END DO 247 242 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 243 END DO 244 END DO 245 ! 248 246 DO jk = 2, jpkm1 ! complete with the level-dependent part 249 emix_ tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk) &247 emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk) & 250 248 & / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 251 249 END DO 252 253 250 ! 254 251 ! Calculate molecular kinematic viscosity 255 252 znu_t(:,:,:) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem) & … … 258 255 znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 259 256 END DO 260 257 ! 261 258 ! Calculate turbulence intensity parameter Reb 262 259 DO jk = 2, jpkm1 263 zReb(:,:,jk) = emix_ tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )264 END DO 265 260 zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 261 END DO 262 ! 266 263 ! Define internal wave-induced diffusivity 267 264 DO jk = 2, jpkm1 268 265 zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6 ! This corresponds to a constant mixing efficiency of 1/6 269 266 END DO 270 267 ! 271 268 IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the 272 269 DO jk = 2, jpkm1 ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes … … 282 279 END DO 283 280 ENDIF 284 281 ! 285 282 DO jk = 2, jpkm1 ! Bound diffusivity by molecular value and 100 cm2/s 286 283 zav_wave(:,:,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp ) * wmask(:,:,jk) 287 284 END DO 288 285 ! 289 286 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave 290 287 ztpc = 0._wp … … 300 297 IF( lk_mpp ) CALL mpp_sum( ztpc ) 301 298 ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 302 299 ! 303 300 IF(lwp) THEN 304 301 WRITE(numout,*) 305 WRITE(numout,*) 'zdf_ tmx : Internal wave-driven mixing (tmx)'302 WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 306 303 WRITE(numout,*) '~~~~~~~ ' 307 304 WRITE(numout,*) … … 339 336 ENDIF 340 337 341 DO jk = 2, jpkm1 !* update momentum diffusivity at wu and wv points342 DO jj = 2, jpjm1343 DO ji = fs_2, fs_jpim1 ! vector opt.344 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj ,jk) ) * wumask(ji,jj,jk)345 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)346 END DO347 END DO348 END DO349 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! lateral boundary condition350 351 338 ! !* output internal wave-driven mixing coefficient 352 339 CALL iom_put( "av_wave", zav_wave ) 353 !* output useful diagnostics: N^2, Kz * N^2 (bflx_ tmx),354 ! vertical integral of rau0 * Kz * N^2 (pcmap_ tmx), energy density (emix_tmx)355 IF( iom_use("bflx_ tmx") .OR. iom_use("pcmap_tmx") ) THEN356 bflx_ tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)357 pcmap_ tmx(:,:) = 0._wp340 !* output useful diagnostics: N^2, Kz * N^2 (bflx_iwm), 341 ! vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm) 342 IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 343 bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 344 pcmap_iwm(:,:) = 0._wp 358 345 DO jk = 2, jpkm1 359 pcmap_ tmx(:,:) = pcmap_tmx(:,:) + e3w_n(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk)360 END DO 361 pcmap_ tmx(:,:) = rau0 * pcmap_tmx(:,:)362 CALL iom_put( "bflx_ tmx", bflx_tmx)363 CALL iom_put( "pcmap_ tmx", pcmap_tmx)346 pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk) 347 END DO 348 pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:) 349 CALL iom_put( "bflx_iwm", bflx_iwm ) 350 CALL iom_put( "pcmap_iwm", pcmap_iwm ) 364 351 ENDIF 365 352 CALL iom_put( "bn2", rn2 ) 366 CALL iom_put( "emix_ tmx", emix_tmx)353 CALL iom_put( "emix_iwm", emix_iwm ) 367 354 368 355 CALL wrk_dealloc( jpi,jpj, zfact, zhdep ) 369 356 CALL wrk_dealloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb ) 370 357 371 IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx- av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)372 ! 373 IF( nn_timing == 1 ) CALL timing_stop('zdf_ tmx')374 ! 375 END SUBROUTINE zdf_ tmx376 377 378 SUBROUTINE zdf_ tmx_init379 !!---------------------------------------------------------------------- 380 !! *** ROUTINE zdf_ tmx_init ***358 IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 359 ! 360 IF( nn_timing == 1 ) CALL timing_stop('zdf_iwm') 361 ! 362 END SUBROUTINE zdf_iwm 363 364 365 SUBROUTINE zdf_iwm_init 366 !!---------------------------------------------------------------------- 367 !! *** ROUTINE zdf_iwm_init *** 381 368 !! 382 369 !! ** Purpose : Initialization of the wave-driven vertical mixing, reading 383 370 !! of input power maps and decay length scales in netcdf files. 384 371 !! 385 !! ** Method : - Read the namzdf_ tmxnamelist and check the parameters372 !! ** Method : - Read the namzdf_iwm namelist and check the parameters 386 373 !! 387 374 !! - Read the input data in NetCDF files : … … 392 379 !! decay scale for critical slope wave-breaking (decay_scale_cri.nc) 393 380 !! 394 !! ** input : - Namlist namzdf_ tmx381 !! ** input : - Namlist namzdf_iwm 395 382 !! - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 396 383 !! decay_scale_bot.nc decay_scale_cri.nc 397 384 !! 398 385 !! ** Action : - Increase by 1 the nstop flag is setting problem encounter 399 !! - Define ebot_ tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx400 !! 401 !! References : de Lavergne et al. 2015, JPO; 2016, in prep.402 !! 386 !! - Define ebot_iwm, epyc_iwm, ecri_iwm, hbot_iwm, hcri_iwm 387 !! 388 !! References : de Lavergne et al. JPO, 2015 ; de Lavergne PhD 2016 389 !! de Lavergne et al. in prep., 2017 403 390 !!---------------------------------------------------------------------- 404 391 INTEGER :: ji, jj, jk ! dummy loop indices … … 407 394 REAL(wp) :: zbot, zpyc, zcri ! local scalars 408 395 !! 409 NAMELIST/namzdf_ tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff410 !!---------------------------------------------------------------------- 411 ! 412 IF( nn_timing == 1 ) CALL timing_start('zdf_ tmx_init')413 ! 414 REWIND( numnam_ref ) ! Namelist namzdf_ tmxin reference namelist : Wave-driven mixing415 READ ( numnam_ref, namzdf_ tmx_new, IOSTAT = ios, ERR = 901)416 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ tmxin reference namelist', lwp )417 ! 418 REWIND( numnam_cfg ) ! Namelist namzdf_ tmxin configuration namelist : Wave-driven mixing419 READ ( numnam_cfg, namzdf_ tmx_new, IOSTAT = ios, ERR = 902 )420 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ tmxin configuration namelist', lwp )421 IF(lwm) WRITE ( numond, namzdf_ tmx_new )396 NAMELIST/namzdf_iwm_new/ nn_zpyc, ln_mevar, ln_tsdiff 397 !!---------------------------------------------------------------------- 398 ! 399 IF( nn_timing == 1 ) CALL timing_start('zdf_iwm_init') 400 ! 401 REWIND( numnam_ref ) ! Namelist namzdf_iwm in reference namelist : Wave-driven mixing 402 READ ( numnam_ref, namzdf_iwm_new, IOSTAT = ios, ERR = 901) 403 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist', lwp ) 404 ! 405 REWIND( numnam_cfg ) ! Namelist namzdf_iwm in configuration namelist : Wave-driven mixing 406 READ ( numnam_cfg, namzdf_iwm_new, IOSTAT = ios, ERR = 902 ) 407 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist', lwp ) 408 IF(lwm) WRITE ( numond, namzdf_iwm_new ) 422 409 ! 423 410 IF(lwp) THEN ! Control print 424 411 WRITE(numout,*) 425 WRITE(numout,*) 'zdf_ tmx_init : internal wave-driven mixing'412 WRITE(numout,*) 'zdf_iwm_init : internal wave-driven mixing' 426 413 WRITE(numout,*) '~~~~~~~~~~~~' 427 WRITE(numout,*) ' Namelist namzdf_ tmx_new : set wave-driven mixing parameters'414 WRITE(numout,*) ' Namelist namzdf_iwm_new : set wave-driven mixing parameters' 428 415 WRITE(numout,*) ' Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 429 416 WRITE(numout,*) ' Variable (T) or constant (F) mixing efficiency = ', ln_mevar … … 435 422 ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 436 423 avmb(:) = 1.4e-6_wp ! viscous molecular value 437 avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_ tmx)424 avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_iwm) 438 425 avtb_2d(:,:) = 1.e0_wp ! uniform 439 426 IF(lwp) THEN ! Control print … … 443 430 ENDIF 444 431 445 ! ! allocate tmxarrays446 IF( zdf_ tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmxarrays' )432 ! ! allocate iwm arrays 433 IF( zdf_iwm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' ) 447 434 ! 448 435 ! ! read necessary fields 449 436 CALL iom_open('mixing_power_bot',inum) ! energy flux for high-mode wave breaking [W/m2] 450 CALL iom_get (inum, jpdom_data, 'field', ebot_ tmx, 1 )437 CALL iom_get (inum, jpdom_data, 'field', ebot_iwm, 1 ) 451 438 CALL iom_close(inum) 452 439 ! 453 440 CALL iom_open('mixing_power_pyc',inum) ! energy flux for pynocline-intensified wave breaking [W/m2] 454 CALL iom_get (inum, jpdom_data, 'field', epyc_ tmx, 1 )441 CALL iom_get (inum, jpdom_data, 'field', epyc_iwm, 1 ) 455 442 CALL iom_close(inum) 456 443 ! 457 444 CALL iom_open('mixing_power_cri',inum) ! energy flux for critical slope wave breaking [W/m2] 458 CALL iom_get (inum, jpdom_data, 'field', ecri_ tmx, 1 )445 CALL iom_get (inum, jpdom_data, 'field', ecri_iwm, 1 ) 459 446 CALL iom_close(inum) 460 447 ! 461 448 CALL iom_open('decay_scale_bot',inum) ! spatially variable decay scale for high-mode wave breaking [m] 462 CALL iom_get (inum, jpdom_data, 'field', hbot_ tmx, 1 )449 CALL iom_get (inum, jpdom_data, 'field', hbot_iwm, 1 ) 463 450 CALL iom_close(inum) 464 451 ! 465 452 CALL iom_open('decay_scale_cri',inum) ! spatially variable decay scale for critical slope wave breaking [m] 466 CALL iom_get (inum, jpdom_data, 'field', hcri_ tmx, 1 )453 CALL iom_get (inum, jpdom_data, 'field', hcri_iwm, 1 ) 467 454 CALL iom_close(inum) 468 455 469 ebot_ tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:)470 epyc_ tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:)471 ecri_ tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:)456 ebot_iwm(:,:) = ebot_iwm(:,:) * ssmask(:,:) 457 epyc_iwm(:,:) = epyc_iwm(:,:) * ssmask(:,:) 458 ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:) 472 459 473 460 ! Set once for all to zero the first and last vertical levels of appropriate variables 474 emix_ tmx(:,:, 1 ) = 0._wp475 emix_ tmx(:,:,jpk) = 0._wp461 emix_iwm (:,:, 1 ) = 0._wp 462 emix_iwm (:,:,jpk) = 0._wp 476 463 zav_ratio(:,:, 1 ) = 0._wp 477 464 zav_ratio(:,:,jpk) = 0._wp … … 479 466 zav_wave (:,:,jpk) = 0._wp 480 467 481 zbot = glob_sum( e1e2t(:,:) * ebot_ tmx(:,:) )482 zpyc = glob_sum( e1e2t(:,:) * epyc_ tmx(:,:) )483 zcri = glob_sum( e1e2t(:,:) * ecri_ tmx(:,:) )468 zbot = glob_sum( e1e2t(:,:) * ebot_iwm(:,:) ) 469 zpyc = glob_sum( e1e2t(:,:) * epyc_iwm(:,:) ) 470 zcri = glob_sum( e1e2t(:,:) * ecri_iwm(:,:) ) 484 471 IF(lwp) THEN 485 472 WRITE(numout,*) ' High-mode wave-breaking energy: ', zbot * 1.e-12_wp, 'TW' … … 488 475 ENDIF 489 476 ! 490 IF( nn_timing == 1 ) CALL timing_stop('zdf_ tmx_init')491 ! 492 END SUBROUTINE zdf_ tmx_init477 IF( nn_timing == 1 ) CALL timing_stop('zdf_iwm_init') 478 ! 479 END SUBROUTINE zdf_iwm_init 493 480 494 481 !!====================================================================== 495 END MODULE zdf tmx482 END MODULE zdfiwm -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
r7953 r7990 2 2 !!====================================================================== 3 3 !! *** MODULE zdfphy *** 4 !! Ocean physics : manager of vertical mixing parametrizations4 !! Vertical ocean physics : manager of all vertical physics packages 5 5 !!====================================================================== 6 6 !! History : 4.0 ! 2017-04 (G. Madec) original code … … 8 8 9 9 !!---------------------------------------------------------------------- 10 !! zdf_phy_init : initialization of all vertical physics pa kages10 !! zdf_phy_init : initialization of all vertical physics packages 11 11 !! zdf_phy : upadate at each time-step the vertical mixing coeff. 12 12 !!---------------------------------------------------------------------- 13 USE par_oce ! mesh and scale factors 14 USE zdf_oce ! TKE vertical mixing 13 USE par_oce ! ocean parameters 14 USE zdf_oce ! vertical physics: shared variables 15 USE zdfbfr ! vertical physics: bottom friction 16 USE zdfric ! vertical physics: Richardson vertical mixing 17 USE zdftke ! vertical physics: TKE vertical mixing 18 USE zdfgls ! vertical physics: GLS vertical mixing 19 USE zdfddm ! vertical physics: double diffusion mixing 20 USE zdfevd ! vertical physics: convection via enhanced vertical diffusion 21 USE zdfiwm ! vertical physics: internal wave-induced mixing 22 USE zdfswm ! vertical physics: surface wave-induced mixing 23 USE zdfmxl ! vertical physics: mixed layer 24 USE tranpc ! convection: non penetrative adjustment 25 USE trc_oce ! variables shared between passive tracer & ocean 15 26 USE sbc_oce ! surface module (only for nn_isf in the option compatibility test) 16 USE zdfbfr ! bottom friction17 USE zdftke ! TKE vertical mixing18 USE zdfgls ! GLS vertical mixing19 USE zdfric ! Richardson vertical mixing20 USE zdfddm ! double diffusion mixing21 USE zdfevd ! enhanced vertical diffusion22 USE zdftmx ! internal tide-induced mixing23 USE zdfqiao !Qiao module wave induced mixing (zdf_qiao routine)24 USE zdfmxl ! Mixed-layer depth (zdf_mxl routine)25 USE tranpc ! convection: non penetrative adjustment26 27 USE sbcrnf ! surface boundary condition: runoff variables 27 28 ! 28 29 USE in_out_manager ! I/O manager 29 30 USE iom ! IOM library 31 USE lbclnk ! lateral boundary conditions 30 32 USE lib_mpp ! distribued memory computing 31 33 … … 33 35 PRIVATE 34 36 35 PUBLIC zdf_phy_init ! routine called by nemogcm.F90 36 PUBLIC zdf_phy ! routine called by step.F90 37 38 39 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 4.0 , NEMO Consortium (2011) 37 PUBLIC zdf_phy_init ! called by nemogcm.F90 38 PUBLIC zdf_phy ! called by step.F90 39 40 INTEGER :: nzdf_phy ! type of vertical closure used 41 ! ! associated indicators 42 INTEGER, PARAMETER :: np_CST = 1 ! Constant Kz 43 INTEGER, PARAMETER :: np_RIC = 2 ! Richardson number dependent Kz 44 INTEGER, PARAMETER :: np_TKE = 3 ! Turbulente Kinetic Eenergy closure scheme for Kz 45 INTEGER, PARAMETER :: np_GLS = 4 ! Generic Length Scale closure scheme for Kz 46 47 !!---------------------------------------------------------------------- 48 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 41 49 !! $Id$ 42 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 51 59 !! 52 60 !! ** Method : Read namelist namzdf, control logicals 61 !! set horizontal shape and vertical profile of background mixing coef. 53 62 !!---------------------------------------------------------------------- 54 63 INTEGER :: ioptio, ios ! local integers … … 58 67 & ln_zdfnpc, nn_npc , nn_npcp, & ! convection : npc 59 68 & ln_zdfddm, rn_avts, rn_hsbfr, & ! double diffusion 60 & ln_zdf tmx, & ! tidalmixing61 & ln_zdf qiao, & ! surface wave-induced mixing69 & ln_zdfswm, & ! surface wave-induced mixing 70 & ln_zdfiwm, & ! internal - - - 62 71 & ln_zdfexp, nn_zdfexp, & ! time-stepping 63 72 & rn_avm0, rn_avt0, nn_avb, nn_havtb ! coefficients 64 65 66 !!org NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp, & 67 !!org & ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp, & 68 !!org & ln_zdfqiao 69 !!---------------------------------------------------------------------- 70 73 !!---------------------------------------------------------------------- 74 ! 75 ! !== Namelist ==! 71 76 REWIND( numnam_ref ) ! Namelist namzdf in reference namelist : Vertical mixing parameters 72 77 READ ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901) 73 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp )74 78 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp ) 79 ! 75 80 REWIND( numnam_cfg ) ! Namelist namzdf in reference namelist : Vertical mixing parameters 76 81 READ ( numnam_cfg, namzdf, IOSTAT = ios, ERR = 902 ) 77 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in configuration namelist', lwp )78 IF(lwm) WRITE ( numond, namzdf )79 80 IF(lwp) THEN !*Parameter print82 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in configuration namelist', lwp ) 83 IF(lwm) WRITE ( numond, namzdf ) 84 ! 85 IF(lwp) THEN ! Parameter print 81 86 WRITE(numout,*) 82 WRITE(numout,*) 'zdf_phy_init 83 WRITE(numout,*) '~~~~~~~~ '87 WRITE(numout,*) 'zdf_phy_init: vertical physics' 88 WRITE(numout,*) '~~~~~~~~~~~~' 84 89 WRITE(numout,*) ' Namelist namzdf : set vertical mixing mixing parameters' 85 90 WRITE(numout,*) ' vertical closure scheme' … … 98 103 WRITE(numout,*) ' maximum avs for dd mixing rn_avts = ', rn_avts 99 104 WRITE(numout,*) ' heat/salt buoyancy flux ratio rn_hsbfr= ', rn_hsbfr 100 WRITE(numout,*) ' surface wave-induced mixing ln_zdfqiao= ', ln_zdfqiao ! surface wave induced mixing 101 WRITE(numout,*) ' tidal mixing ln_zdftmx = ', ln_zdftmx 102 WRITE(numout,*) ' time splitting / backward scheme ln_zdfexp = ', ln_zdfexp 105 WRITE(numout,*) ' gravity wave-induced mixing' 106 WRITE(numout,*) ' surface wave (Qiao et al 2010) ln_zdfswm = ', ln_zdfswm ! surface wave induced mixing 107 WRITE(numout,*) ' internal wave (de Lavergne et al 2017) ln_zdfiwm = ', ln_zdfiwm 108 WRITE(numout,*) ' time-steping scheme' 109 WRITE(numout,*) ' time splitting (T) / implicit (F) ln_zdfexp = ', ln_zdfexp 103 110 WRITE(numout,*) ' number of sub-time step (ln_zdfexp=T) nn_zdfexp = ', nn_zdfexp 104 111 WRITE(numout,*) ' coefficients : ' 105 WRITE(numout,*) ' vertical eddy viscosity rn_avm0 = ', rn_avm0 106 WRITE(numout,*) ' vertical eddy diffusivity rn_avt0 = ', rn_avt0 107 WRITE(numout,*) ' constant background or profile nn_avb = ', nn_avb 108 WRITE(numout,*) ' horizontal variation for avtb nn_havtb = ', nn_havtb 109 ENDIF 110 111 !!gm IF(ln_zdfddm) THEN ! double diffusive mixing' 112 ! avs(:,:,:) = rn_avt0 * wmask(:,:,:) 113 !!gm ENDIF 114 115 ! !* Parameter & logical controls 116 ! ! ---------------------------- 117 ! 118 ! ! ... check of vertical mixing scheme on tracers 119 ! ==> will be done in trazdf module 120 ! 121 ! ! ... check of mixing coefficient 122 IF(lwp) WRITE(numout,*) 123 IF(lwp) WRITE(numout,*) ' vertical mixing option :' 124 ioptio = 0 125 IF( ln_zdfcst ) THEN 126 IF(lwp) WRITE(numout,*) ' constant eddy diffusion coefficients' 127 ioptio = ioptio+1 128 ENDIF 129 IF( ln_zdfric ) THEN 130 IF(lwp) WRITE(numout,*) ' Richardson dependent eddy coefficients' 131 ioptio = ioptio+1 132 ENDIF 133 IF( ln_zdftke ) THEN 134 IF(lwp) WRITE(numout,*) ' TKE dependent eddy coefficients' 135 ioptio = ioptio+1 136 ENDIF 137 IF( ln_zdfgls ) THEN 138 IF(lwp) WRITE(numout,*) ' GLS dependent eddy coefficients' 139 ioptio = ioptio+1 140 ENDIF 141 IF( ioptio == 0 .OR. ioptio > 1 ) & 142 & CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 143 IF( ( ln_zdfric .OR. ln_zdfgls ) .AND. ln_isfcav ) & 144 & CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 145 ! 146 ! ! ... Convection 147 IF(lwp) WRITE(numout,*) 148 IF(lwp) WRITE(numout,*) ' convection :' 149 ! 150 #if defined key_top 151 IF( ln_zdfnpc ) CALL ctl_stop( ' zdf_phy_init: npc scheme is not working with key_top' ) 152 #endif 153 ! 154 ioptio = 0 155 IF( ln_zdfnpc ) THEN 156 IF(lwp) WRITE(numout,*) ' use non penetrative convective scheme' 157 ioptio = ioptio+1 158 ENDIF 159 IF( ln_zdfevd ) THEN 160 IF(lwp) WRITE(numout,*) ' use enhanced vertical dif. scheme' 161 ioptio = ioptio+1 162 ENDIF 163 IF( ln_zdftke ) THEN 164 IF(lwp) WRITE(numout,*) ' use the 1.5 turbulent closure' 165 ENDIF 166 IF( ln_zdfgls ) THEN 167 IF(lwp) WRITE(numout,*) ' use the GLS closure scheme' 168 ENDIF 169 IF ( ioptio > 1 ) CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' ) 170 IF( ioptio == 0 .AND. .NOT.( ln_zdftke .OR. ln_zdfgls ) ) & 171 CALL ctl_stop( ' except for TKE or GLS physics, a convection scheme is', & 172 & ' required: ln_zdfevd or ln_zdfnpc logicals' ) 173 174 ! !* Background eddy viscosity and diffusivity profil 175 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter 112 WRITE(numout,*) ' vertical eddy viscosity rn_avm0 = ', rn_avm0 113 WRITE(numout,*) ' vertical eddy diffusivity rn_avt0 = ', rn_avt0 114 WRITE(numout,*) ' constant background or profile nn_avb = ', nn_avb 115 WRITE(numout,*) ' horizontal variation for avtb nn_havtb = ', nn_havtb 116 ENDIF 117 118 ! !== Background eddy viscosity and diffusivity ==! 119 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter 176 120 avmb(:) = rn_avm0 177 121 avtb(:) = rn_avt0 178 ELSE 122 ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 179 123 avmb(:) = rn_avm0 180 124 avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:) ! m2/s 181 125 IF(ln_sco .AND. lwp) CALL ctl_warn( 'avtb profile not valid in sco' ) 182 126 ENDIF 183 ! 184 IF( ln_rstart ) THEN ! Read avmb, avtb in restart (if exist) 185 ! if ln_traadv_cen, avmb, avtb have been modified in traadv_cen2 module. 186 ! To ensure the restartability, avmb & avtb are written in the restart 187 ! file in traadv_cen2 end read here. 188 IF( iom_varid( numror, 'avmb', ldstop = .FALSE. ) > 0 ) THEN 189 CALL iom_get( numror, jpdom_unknown, 'avmb', avmb ) 190 CALL iom_get( numror, jpdom_unknown, 'avtb', avtb ) 191 ENDIF 192 ENDIF 193 ! ! 2D shape of the avtb 194 avtb_2d(:,:) = 1.e0 ! uniform 195 ! 196 IF( nn_havtb == 1 ) THEN ! decrease avtb in the equatorial band 197 ! -15S -5S : linear decrease from avt0 to avt0/10. 198 ! -5S +5N : cst value avt0/10. 199 ! 5N 15N : linear increase from avt0/10, to avt0 127 ! ! 2D shape of the avtb 128 avtb_2d(:,:) = 1._wp ! uniform 129 ! 130 IF( nn_havtb == 1 ) THEN ! decrease avtb by a factor of ten in the equatorial band 131 ! ! -15S -5S : linear decrease from avt0 to avt0/10. 132 ! ! -5S +5N : cst value avt0/10. 133 ! ! 5N 15N : linear increase from avt0/10, to avt0 200 134 WHERE(-15. <= gphit .AND. gphit < -5 ) avtb_2d = (1. - 0.09 * (gphit + 15.)) 201 135 WHERE( -5. <= gphit .AND. gphit < 5 ) avtb_2d = 0.1 … … 203 137 ENDIF 204 138 ! 205 206 !!gm moved into zdf_phy_init 207 ! 208 CALL zdf_bfr_init ! bottom friction 209 210 ioptio = 0 !== type of vertical turbulent closure ==! (set nzdfphy) 211 ! 212 ! IF( ln_zdfcst ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_CST ; ENDIF 213 ! IF( ln_zdfric ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_RIC ; CALL zdf_ric_init ; ENDIF 214 ! IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init ; ENDIF 215 ! IF( ln_zdfgls ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_GLS ; CALL zdf_gls_init ; ENDIF 216 217 218 ! 219 IF( ln_zdfric ) CALL zdf_ric_init ! Richardson number dependent Kz 220 IF( ln_zdftke ) CALL zdf_tke_init ! TKE closure scheme 221 IF( ln_zdfgls ) CALL zdf_gls_init ! GLS closure scheme 222 IF( ln_zdftmx ) CALL zdf_tmx_init ! tidal vertical mixing 223 !!gm 139 ! ! set 1st/last level Av to zero one for all 140 avt(:,:,1) = 0._wp ; avs(:,:,1) = 0._wp ; avm(:,:,1) = 0._wp 141 142 ! !== Convection ==! 143 ! 144 IF( ln_zdfnpc .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' ) 145 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 146 IF(lwp) THEN 147 WRITE(numout,*) 148 IF ( ln_zdfnpc ) THEN ; WRITE(numout,*) ' convection: use non penetrative convective scheme' 149 ELSEIF( ln_zdfevd ) THEN ; WRITE(numout,*) ' convection: use enhanced vertical diffusion scheme' 150 ELSE ; WRITE(numout,*) ' convection: no specific scheme used' 151 ENDIF 152 ENDIF 153 154 IF(lwp) THEN !== Double Diffusion Mixing parameterization ==! (ddm) 155 WRITE(numout,*) 156 IF( ln_zdfddm ) THEN ; WRITE(numout,*) ' use double diffusive mixing: avs /= avt' 157 ELSE ; WRITE(numout,*) ' No double diffusive mixing: avs = avt' 158 ENDIF 159 ENDIF 160 161 ! !== type of vertical turbulent closure ==! (set nzdf_phy) 162 ioptio = 0 163 IF( ln_zdfcst ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_CST ; ENDIF 164 IF( ln_zdfric ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_RIC ; CALL zdf_ric_init ; ENDIF 165 IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init ; ENDIF 166 IF( ln_zdfgls ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_GLS ; CALL zdf_gls_init ; ENDIF 167 ! 168 IF( ioptio /= 1 ) CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' ) 169 IF( ln_isfcav ) THEN 170 IF( ln_zdfric .OR. ln_zdfgls ) CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' ) 171 ENDIF 172 173 ! !== gravity wave-driven mixing ==! 174 IF( ln_zdfiwm ) CALL zdf_iwm_init ! internal wave-driven mixing 175 IF( ln_zdfswm ) CALL zdf_swm_init ! surface wave-driven mixing 176 177 ! !== bottom friction ==! 178 CALL zdf_bfr_init 179 ! 180 ! !== time-stepping ==! 181 ! Check/update of time stepping done in dynzdf_init/trazdf_init 182 !!gm move it here ? 224 183 ! 225 184 END SUBROUTINE zdf_phy_init 226 185 227 186 228 SUBROUTINE zdf_phy( k stp)187 SUBROUTINE zdf_phy( kt ) 229 188 !!---------------------------------------------------------------------- 230 189 !! *** ROUTINE zdf_phy *** … … 238 197 !! bottom stress..... <<<<====verifier ! 239 198 !!---------------------------------------------------------------------- 240 INTEGER, INTENT(in) :: kstp ! ocean time-step index 241 ! 242 INTEGER :: ji, jj, jk ! dummy loop indice 243 !!---------------------------------------------------------------------- 244 ! 245 CALL zdf_bfr( kstp ) ! bottom friction (if quadratic) 246 ! ! Vertical eddy viscosity and diffusivity coefficients 247 IF( ln_zdfric ) CALL zdf_ric ( kstp ) ! Richardson number dependent Kz 248 IF( ln_zdftke ) CALL zdf_tke ( kstp ) ! TKE closure scheme for Kz 249 IF( ln_zdfgls ) CALL zdf_gls ( kstp ) ! GLS closure scheme for Kz 250 IF( ln_zdfqiao ) CALL zdf_qiao( kstp ) ! Qiao vertical mixing 251 ! 252 IF( ln_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 253 avt (:,:,:) = rn_avt0 * wmask (:,:,:) 254 avm (:,:,:) = rn_avm0 * wmask (:,:,:) 255 avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 256 avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 257 ENDIF 258 ! 259 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 260 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO 261 ENDIF 262 ! 263 IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity 264 ! 265 IF( ln_zdfddm ) THEN ! double diffusive mixing 266 CALL zdf_ddm( kstp ) 267 ELSE ! avs=avt 268 DO jk = 2, jpkm1 ; avs(:,:,jk) = avt(:,:,jk) ; END DO 269 ENDIF 270 ! 271 IF( ln_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing 272 273 CALL zdf_mxl( kstp ) ! mixed layer depth 274 199 INTEGER, INTENT(in) :: kt ! ocean time-step index 200 ! 201 INTEGER :: ji, jj, jk ! dummy loop indice 202 !! --------------------------------------------------------------------- 203 ! 204 CALL zdf_bfr( kt ) !* bottom friction (if quadratic) 205 ! 206 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 207 CASE( np_RIC ) ; CALL zdf_ric( kt ) ! Richardson number dependent Kz 208 CASE( np_TKE ) ; CALL zdf_tke( kt ) ! TKE closure scheme for Kz 209 CASE( np_GLS ) ; CALL zdf_gls( kt ) ! GLS closure scheme for Kz 210 CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 211 avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 212 avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 213 END SELECT 214 ! 215 IF( ln_rnf_mouth ) THEN !* increase diffusivity at rivers mouths 216 DO jk = 2, nkrnf 217 avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk) 218 END DO 219 ENDIF 220 ! 221 IF( ln_zdfevd ) CALL zdf_evd( kt ) !* convection: enhanced vertical eddy diffusivity 222 ! 223 ! !* double diffusive mixing 224 IF( ln_zdfddm ) THEN ! update avt and compute avs 225 CALL zdf_ddm( kt ) 226 ELSE ! same mixing on all tracers 227 avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 228 ENDIF 229 ! 230 ! !* wave-induced mixing 231 IF( ln_zdfswm ) CALL zdf_swm( kt ) ! surface wave (Qiao et al. 2004) 232 IF( ln_zdfiwm ) CALL zdf_iwm( kt ) ! internal wave (de Lavergne et al 2017) 233 234 235 ! !* Lateral boundary conditions (sign unchanged) 236 CALL lbc_lnk( avm, 'W', 1. ) 237 CALL lbc_lnk( avt, 'W', 1. ) 238 CALL lbc_lnk( avs, 'W', 1. ) 239 ! 240 !!gm ===>>> TO BE REMOVED 241 DO jk = 1, jpkm1 !* vertical eddy viscosity at wu- and wv-points 242 DO jj = 2, jpjm1 243 DO ji = 2, jpim1 244 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 245 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 250 !!gm end 251 252 253 CALL zdf_mxl( kt ) !* mixed layer depth, and level 254 255 !!gm !==>> to be moved in zdftke & zdfgls 275 256 ! write TKE or GLS information in the restart file 276 IF( lrst_oce .AND. ln_zdftke ) CALL tke_rst( k stp, 'WRITE' )277 IF( lrst_oce .AND. ln_zdfgls ) CALL gls_rst( k stp, 'WRITE' )278 ! 257 IF( lrst_oce .AND. ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 258 IF( lrst_oce .AND. ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 259 ! 279 260 END SUBROUTINE zdf_phy 280 261 -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r7953 r7990 5 5 !! Richardson number dependent formulation 6 6 !!====================================================================== 7 !! History : OPA ! 1987-09 (P. Andrich) Original code 8 !! 4.0 ! 1991-11 (G. Madec) 9 !! 7.0 ! 1996-01 (G. Madec) complete rewriting of multitasking suppression of common work arrays 10 !! 8.0 ! 1997-06 (G. Madec) complete rewriting of zdfmix 11 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 12 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 13 !! 3.3.1! 2011-09 (P. Oddo) Mixed layer depth parameterization 7 !! History : OPA ! 1987-09 (P. Andrich) Original code 8 !! 4.0 ! 1991-11 (G. Madec) 9 !! 7.0 ! 1996-01 (G. Madec) complete rewriting of multitasking suppression of common work arrays 10 !! 8.0 ! 1997-06 (G. Madec) complete rewriting of zdfmix 11 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 12 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 13 !! 3.3.1! 2011-09 (P. Oddo) Mixed layer depth parameterization 14 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 14 15 !!---------------------------------------------------------------------- 15 16 … … 21 22 USE dom_oce ! ocean space and time domain variables 22 23 USE zdf_oce ! ocean vertical physics 24 USE phycst ! physical constants 25 USE sbc_oce, ONLY : taum 26 USE eosbn2 , ONLY : neos 27 ! 23 28 USE in_out_manager ! I/O manager 24 29 USE lbclnk ! ocean lateral boundary condition (or mpp link) 25 30 USE lib_mpp ! MPP library 26 USE wrk_nemo ! work arrays27 31 USE timing ! Timing 28 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 33 30 USE eosbn2, ONLY : neos31 34 32 35 IMPLICIT NONE … … 36 39 PUBLIC zdf_ric_init ! called by opa.F90 37 40 38 ! !!* Namelist namzdf_ric : Richardson number dependent Kz * 39 INTEGER :: nn_ric ! coefficient of the parameterization 40 REAL(wp) :: rn_avmri ! maximum value of the vertical eddy viscosity 41 REAL(wp) :: rn_alp ! coefficient of the parameterization 42 REAL(wp) :: rn_ekmfc ! Ekman Factor Coeff 43 REAL(wp) :: rn_mldmin ! minimum mixed layer (ML) depth 44 REAL(wp) :: rn_mldmax ! maximum mixed layer depth 45 REAL(wp) :: rn_wtmix ! Vertical eddy Diff. in the ML 46 REAL(wp) :: rn_wvmix ! Vertical eddy Visc. in the ML 47 LOGICAL :: ln_mldw ! Use or not the MLD parameters 48 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmric !: coef. for the horizontal mean at t-point 41 ! !!* Namelist namzdf_ric : Richardson number dependent Kz * 42 INTEGER :: nn_ric ! coefficient of the parameterization 43 REAL(wp) :: rn_avmri ! maximum value of the vertical eddy viscosity 44 REAL(wp) :: rn_alp ! coefficient of the parameterization 45 REAL(wp) :: rn_ekmfc ! Ekman Factor Coeff 46 REAL(wp) :: rn_mldmin ! minimum mixed layer (ML) depth 47 REAL(wp) :: rn_mldmax ! maximum mixed layer depth 48 REAL(wp) :: rn_wtmix ! Vertical eddy Diff. in the ML 49 REAL(wp) :: rn_wvmix ! Vertical eddy Visc. in the ML 50 LOGICAL :: ln_mldw ! Use or not the MLD parameters 50 51 51 52 !! * Substitutions … … 57 58 !!---------------------------------------------------------------------- 58 59 CONTAINS 59 60 INTEGER FUNCTION zdf_ric_alloc()61 !!----------------------------------------------------------------------62 !! *** FUNCTION zdf_ric_alloc ***63 !!----------------------------------------------------------------------64 ALLOCATE( tmric(jpi,jpj,jpk) , STAT= zdf_ric_alloc )65 !66 IF( lk_mpp ) CALL mpp_sum ( zdf_ric_alloc )67 IF( zdf_ric_alloc /= 0 ) CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')68 END FUNCTION zdf_ric_alloc69 70 60 71 61 SUBROUTINE zdf_ric( kt ) … … 105 95 !! and bottom value already set in zdfphy.F90 106 96 !! 97 !! ** Action : avm, avt mixing coeff (inner domain values only) 98 !! 107 99 !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 108 100 !! PFJ Lermusiaux 2001. 109 101 !!---------------------------------------------------------------------- 110 USE phycst, ONLY: rsmall,rau0 111 USE sbc_oce, ONLY: taum 112 !! 113 INTEGER, INTENT( in ) :: kt ! ocean time-step 114 !! 115 INTEGER :: ji, jj, jk ! dummy loop indices 116 REAL(wp) :: zcoef, zdku, zdkv, zri, z05alp, zflageos ! temporary scalars 117 REAL(wp) :: zrhos, zustar 118 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, ekm_dep 119 !!---------------------------------------------------------------------- 120 ! 121 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 122 ! 123 CALL wrk_alloc( jpi,jpj, zwx, ekm_dep ) 124 ! ! =============== 125 DO jk = 2, jpkm1 ! Horizontal slab 126 ! ! =============== 127 ! Richardson number (put in zwx(ji,jj)) 128 ! ----------------- 102 INTEGER, INTENT(in) :: kt ! ocean time-step 103 !! 104 INTEGER :: ji, jj, jk ! dummy loop indices 105 LOGICAL :: ll_av_weight = .TRUE. ! local logical 106 REAL(wp) :: zsh2, zcfRi, zav, zustar ! local scalars 107 REAL(wp), DIMENSION(jpi,jpj) :: zdu2, zdv2, zh_ekm ! 2D workspace 108 !!---------------------------------------------------------------------- 109 ! 110 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 111 ! 112 DO jk = 2, jpkm1 !== avm and avt = F(Richardson number) ==! 113 ! 114 IF( ll_av_weight ) THEN !== viscosity weighted shear & stratification terms 115 ! 116 DO jj = 1, jpjm1 !* viscosity weighted shear term 117 DO ji = 1, jpim1 118 zdu2(ji,jj) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) * wumask(ji,jj,jk) & 119 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 120 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 121 zdv2(ji,jj) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) * wumask(ji,jj,jk) & 122 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 123 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 124 END DO 125 END DO 129 126 DO jj = 2, jpjm1 130 DO ji = fs_2, fs_jpim1 131 zcoef = 0.5 / e3w_n(ji,jj,jk) 132 ! ! shear of horizontal velocity 133 zdku = zcoef * ( ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1) & 134 & -ub(ji-1,jj,jk ) - ub(ji,jj,jk ) ) 135 zdkv = zcoef * ( vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1) & 136 & -vb(ji,jj-1,jk ) - vb(ji,jj,jk ) ) 137 ! ! richardson number (minimum value set to zero) 138 zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 139 zwx(ji,jj) = MAX( zri, 0.e0 ) 140 END DO 141 END DO 142 CALL lbc_lnk( zwx, 'W', 1. ) ! Boundary condition (sign unchanged) 143 144 ! Vertical eddy viscosity and diffusivity coefficients 145 ! ------------------------------------------------------- 146 z05alp = 0.5_wp * rn_alp 147 DO jj = 1, jpjm1 ! Eddy viscosity coefficients (avm) 148 DO ji = 1, fs_jpim1 149 avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric 150 avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric 151 END DO 152 END DO 153 DO jj = 2, jpjm1 ! Eddy diffusivity coefficients (avt) 154 DO ji = fs_2, fs_jpim1 155 avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) ) & 156 & * ( avmu(ji,jj,jk) + avmu(ji-1,jj,jk) & 157 & + avmv(ji,jj,jk) + avmv(ji,jj-1,jk) ) & 158 & + avtb(jk) * tmask(ji,jj,jk) 159 END DO 160 END DO 161 DO jj = 2, jpjm1 ! Add the background coefficient on eddy viscosity 162 DO ji = fs_2, fs_jpim1 163 avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk) 164 avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk) 165 END DO 166 END DO 167 ! ! =============== 168 END DO ! End of slab 169 ! ! =============== 170 ! 171 IF( ln_mldw ) THEN 172 173 ! Compute Ekman depth from wind stress forcing. 174 ! ------------------------------------------------------- 175 zflageos = ( 0.5 + SIGN( 0.5, neos - 1. ) ) * rau0 176 DO jj = 2, jpjm1 177 DO ji = fs_2, fs_jpim1 178 zrhos = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) ) 179 zustar = SQRT( taum(ji,jj) / ( zrhos + rsmall ) ) 180 ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) 181 ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed 182 ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed 183 END DO 127 DO ji = 2, jpim1 !* Richardson number dependent coefficient 128 ! ! shear of horizontal velocity 129 zsh2 = ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 130 & + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 131 ! ! coefficient = F(richardson number) 132 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avt(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) ) ) 133 ! 134 ! !* Vertical eddy viscosity and diffusivity coefficients 135 zav = rn_avmri * zcfRi**nn_ric 136 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 137 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 138 END DO 139 END DO 140 ELSE !== classical Richardson number definition ==! 141 DO jj = 1, jpjm1 !* shear term 142 DO ji = 1, jpim1 143 zdu2(ji,jj) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 144 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 145 zdv2(ji,jj) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 146 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 147 END DO 148 END DO 149 DO jj = 2, jpjm1 150 DO ji = 2, jpim1 !* Richardson number dependent coefficient 151 ! ! shear of horizontal velocity 152 zsh2 = ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 153 & + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 154 ! ! coefficient = F(richardson number) 155 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) ) ) 156 ! 157 ! !* Vertical eddy viscosity and diffusivity coefficients 158 zav = rn_avmri * zcfRi**nn_ric 159 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 160 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 161 END DO 162 END DO 163 ENDIF 164 ! 184 165 END DO 185 186 ! In the first model level vertical diff/visc coeff.s 187 ! are always equal to the namelist values rn_wtmix/rn_wvmix 188 ! ------------------------------------------------------- 189 DO jj = 2, jpjm1 190 DO ji = fs_2, fs_jpim1 191 avmv(ji,jj,1) = MAX( avmv(ji,jj,1), rn_wvmix ) 192 avmu(ji,jj,1) = MAX( avmu(ji,jj,1), rn_wvmix ) 193 avt( ji,jj,1) = MAX( avt(ji,jj,1), rn_wtmix ) 194 END DO 195 END DO 196 197 ! Force the vertical mixing coef within the Ekman depth 198 ! ------------------------------------------------------- 199 DO jk = 2, jpkm1 200 DO jj = 2, jpjm1 201 DO ji = fs_2, fs_jpim1 202 IF( gdept_n(ji,jj,jk) < ekm_dep(ji,jj) ) THEN 203 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix ) 204 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix ) 205 avt( ji,jj,jk) = MAX( avt(ji,jj,jk), rn_wtmix ) 206 ENDIF 207 END DO 208 END DO 209 END DO 210 211 DO jk = 1, jpkm1 212 DO jj = 2, jpjm1 213 DO ji = fs_2, fs_jpim1 214 avmv(ji,jj,jk) = avmv(ji,jj,jk) * vmask(ji,jj,jk) 215 avmu(ji,jj,jk) = avmu(ji,jj,jk) * umask(ji,jj,jk) 216 avt( ji,jj,jk) = avt( ji,jj,jk) * tmask(ji,jj,jk) 217 END DO 218 END DO 219 END DO 220 221 ENDIF 222 223 CALL lbc_lnk( avt , 'W', 1. ) ! Boundary conditions (unchanged sign) 224 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) 225 ! 226 CALL wrk_dealloc( jpi,jpj, zwx, ekm_dep ) 227 ! 228 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') 166 ! 167 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 168 ! 169 DO jj = 2, jpjm1 !* Ekman depth 170 DO ji = 2, jpim1 171 zustar = SQRT( taum(ji,jj) * r1_rau0 ) 172 zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 173 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax ) ) ! set allowed rang 174 END DO 175 END DO 176 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer 177 DO jj = 2, jpjm1 178 DO ji = 2, jpim1 179 IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 180 avm(ji,jj,jk) = MAX( avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) 181 avt(ji,jj,jk) = MAX( avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk) 182 ENDIF 183 END DO 184 END DO 185 END DO 186 END IF 187 ! 188 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') 229 189 ! 230 190 END SUBROUTINE zdf_ric … … 264 224 WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme' 265 225 WRITE(numout,*) '~~~~~~~' 266 WRITE(numout,*) ' Namelist namzdf_ric : set Kz (Ri) parameters'267 WRITE(numout,*) ' maximum vertical viscosity rn_avmri = ', rn_avmri268 WRITE(numout,*) ' coefficient rn_alp = ', rn_alp269 WRITE(numout,*) ' coefficientnn_ric = ', nn_ric270 WRITE(numout,*) ' Ekman Factor Coeff rn_ekmfc = ', rn_ekmfc271 WRITE(numout,*) ' minimum mixed layer depth rn_mldmin = ', rn_mldmin272 WRITE(numout,*) ' maximum mixed layer depth rn_mldmax = ', rn_mldmax273 WRITE(numout,*) ' Vertical eddy Diff. in the ML rn_wtmix = ', rn_wtmix274 WRITE(numout,*) ' Vertical eddy Visc. in the ML rn_wvmix = ', rn_wvmix275 WRITE(numout,*) ' Use the MLD parameterization ln_mldw = ', ln_mldw226 WRITE(numout,*) ' Namelist namzdf_ric : set Kz=F(Ri) parameters' 227 WRITE(numout,*) ' maximum vertical viscosity rn_avmri = ', rn_avmri 228 WRITE(numout,*) ' coefficient rn_alp = ', rn_alp 229 WRITE(numout,*) ' exponent nn_ric = ', nn_ric 230 WRITE(numout,*) ' Ekman layer enhanced mixing ln_mldw = ', ln_mldw 231 WRITE(numout,*) ' Ekman Factor Coeff rn_ekmfc = ', rn_ekmfc 232 WRITE(numout,*) ' minimum mixed layer depth rn_mldmin = ', rn_mldmin 233 WRITE(numout,*) ' maximum mixed layer depth rn_mldmax = ', rn_mldmax 234 WRITE(numout,*) ' Vertical eddy Diff. in the ML rn_wtmix = ', rn_wtmix 235 WRITE(numout,*) ' Vertical eddy Visc. in the ML rn_wvmix = ', rn_wvmix 276 236 ENDIF 277 237 ! 278 ! ! allocate zdfric arrays279 IF( zdf_ric_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )280 !281 DO jk = 1, jpk ! weighting mean array tmric for 4 T-points282 DO jj = 2, jpj ! which accounts for coastal boundary conditions283 DO ji = 2, jpi284 tmric(ji,jj,jk) = tmask(ji,jj,jk) &285 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &286 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )287 END DO288 END DO289 END DO290 tmric(:,1,:) = 0._wp291 !292 238 DO jk = 1, jpk ! Initialization of vertical eddy coef. to the background value 293 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 294 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 295 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 239 avt(:,:,jk) = avtb(jk) * tmask(:,:,jk) 240 avm(:,:,jk) = avmb(jk) * umask(:,:,jk) 296 241 END DO 297 242 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90
r7967 r7990 1 MODULE zdf qiao1 MODULE zdfswm 2 2 !!====================================================================== 3 !! *** MODULE zdf qiao***4 !! Qiao module : vertical mixing enhancement due to surface waves3 !! *** MODULE zdfswm *** 4 !! vertical physics : surface wave-induced mixing 5 5 !!====================================================================== 6 6 !! History : 3.6 ! 2014-10 (E. Clementi) Original code 7 !!---------------------------------------------------------------------- 8 !! zdf_qiao : compute Qiao parameters 7 !! 4.0 ! 2017-04 (G. Madec) debug + simplifications 9 8 !!---------------------------------------------------------------------- 10 9 11 USE in_out_manager ! I/O manager 12 USE lib_mpp ! distribued memory computing library 13 USE sbc_oce ! Surface boundary condition: ocean fields 14 USE zdf_oce 15 USE sbcwave ! wave module 16 USE dom_oce 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 10 !!---------------------------------------------------------------------- 11 !! zdf_swm : update Kz due to surface wave-induced mixing 12 !! zdf_swm_init : initilisation 13 !!---------------------------------------------------------------------- 14 USE dom_oce ! ocean domain variable 15 USE zdf_oce ! vertical physics: mixing coefficients 16 USE sbc_oce ! Surface boundary condition: ocean fields 17 USE sbcwave ! wave module 18 ! 19 USE in_out_manager ! I/O manager 20 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 USE lib_mpp ! distribued memory computing library 18 22 19 23 IMPLICIT NONE 20 24 PRIVATE 21 25 22 PUBLIC zdf_qiao ! routine called in step 26 PUBLIC zdf_swm ! routine called in zdp_phy 27 PUBLIC zdf_swm_init ! routine called in zdf_phy_init 23 28 24 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv25 26 !! * Substitutions27 # include "vectopt_loop_substitute.h90"28 29 !!---------------------------------------------------------------------- 29 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)30 !! $Id: 30 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 31 !! $Id:$ 31 32 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 32 33 !!---------------------------------------------------------------------- 33 34 34 CONTAINS 35 35 36 SUBROUTINE zdf_ qiao( kt )36 SUBROUTINE zdf_swm( kt ) 37 37 !!--------------------------------------------------------------------- 38 !! *** ROUTINE zdf_ qiao***38 !! *** ROUTINE zdf_swm *** 39 39 !! 40 !! ** Purpose :Compute the Qiaoterm (qbv) to be added to40 !! ** Purpose :Compute the swm term (qbv) to be added to 41 41 !! vertical viscosity and diffusivity coeffs. 42 42 !! 43 !! ** Method :qbv = alpha * A * Us(0) * exp (3 * k * z) 43 !! ** Method : Compute the swm term Bv (zqb) and added it to 44 !! vertical viscosity and diffusivity coefficients 45 !! zqb = alpha * A * Us(0) * exp (3 * k * z) 46 !! where alpha is set here to 1 44 47 !! 45 !! ** action : Compute the Qiao wave dependent term46 !! only if ln_wave=.true.48 !! ** action : avt, avs, avm updated by the surface wave-induced mixing 49 !! (inner domain only) 47 50 !! 51 !! reference : Qiao et al. GRL, 2004 48 52 !!--------------------------------------------------------------------- 49 INTEGER, INTENT( in ) ::kt ! ocean time step53 INTEGER, INTENT(in) :: kt ! ocean time step 50 54 ! 51 INTEGER :: jj, ji, jk ! dummy loop indices 55 INTEGER :: ji, jj, jk ! dummy loop indices 56 REAL(wp):: zcoef, zqb ! local scalar 52 57 !!--------------------------------------------------------------------- 53 58 ! 54 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 55 IF( .NOT. ( ln_wave .AND. ln_sdw ) ) & 56 & CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave & 57 & and ln_sdw have to be activated') 58 IF( zdf_qiao_alloc() /= 0 ) & 59 & CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' ) 60 ENDIF 61 62 ! 63 ! Compute the Qiao term Bv (qbv) to be added to 64 ! vertical viscosity and diffusivity 65 ! qbv = alpha * A * Us(0) * exp (3 * k * z) 66 ! alpha here is set to 1 67 !--------------------------------------------------------------------------------- 68 ! 69 !!gm Comment: I don't understand the use of min of 4 gdepw_n to define a quantity at w-point 70 !!gm ==>> this is an error.... 71 DO jk = 1, jpk 72 DO jj = 1, jpjm1 73 DO ji = 1, fs_jpim1 74 qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) * & 75 & EXP(3.0 * wnum(ji,jj) * & 76 & (-MIN( gdepw_n(ji ,jj ,jk), gdepw_n(ji+1,jj ,jk), & 77 & gdepw_n(ji ,jj+1,jk), gdepw_n(ji+1,jj+1,jk)))) & 78 & * wmask(ji,jj,jk) 59 zcoef = 1._wp * 0.353553_wp 60 DO jk = 2, jpkm1 61 DO jj = 2, jpjm1 62 DO ji = 2, jpim1 63 zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk) 64 ! 65 avt(ji,jj,jk) = avt(ji,jj,jk) + zqb 66 avs(ji,jj,jk) = avs(ji,jj,jk) + zqb 67 avm(ji,jj,jk) = avm(ji,jj,jk) + zqb 79 68 END DO 80 69 END DO 81 70 END DO 82 71 ! 83 CALL lbc_lnk( qbv, 'W', 1. ) ! Lateral boundary conditions 84 72 END SUBROUTINE zdf_swm 73 74 75 SUBROUTINE zdf_swm_init 76 !!--------------------------------------------------------------------- 77 !! *** ROUTINE zdf_swm_init *** 78 !! 79 !! ** Purpose : surface wave-induced mixing initialisation 80 !! 81 !! ** Method : check the availability of surface wave fields 82 !!--------------------------------------------------------------------- 85 83 ! 86 ! Interpolate Qiao parameter qbv into the grid_U and grid_V 87 !---------------------------------------------------------- 84 IF(lwp) THEN ! Control print 85 WRITE(numout,*) 86 WRITE(numout,*) 'zdf_swm_init : surface wave-driven mixing' 87 WRITE(numout,*) '~~~~~~~~~~~~' 88 ENDIF 89 IF( .NOT.ln_wave .OR. & 90 & .NOT.ln_sdw ) CALL ctl_stop ( 'zdf_swm_init: ln_zdfswm=T but ln_wave and ln_sdw /= T') 88 91 ! 89 DO jk = 1, jpk 90 DO jj = 1, jpjm1 91 DO ji = 1, fs_jpim1 92 qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk) * & 93 & ( qbv(ji,jj,jk) + qbv(ji+1,jj ,jk) ) 94 qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk) * & 95 & ( qbv(ji,jj,jk) + qbv(ji ,jj+1,jk) ) 96 END DO 97 END DO 98 END DO 99 ! 100 CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. ) ! Lateral boundary conditions 92 END SUBROUTINE zdf_swm_init 101 93 102 ! Enhance vertical mixing coeff.103 !-------------------------------104 !105 !!gm with double diffusion activated, avs is not updated...106 !!gm =====>>> BUG107 DO jk = 1, jpkm1108 DO jj = 1, jpj109 DO ji = 1, jpi110 avmu(ji,jj,jk) = ( avmu(ji,jj,jk) + qbvu(ji,jj,jk) ) * umask(ji,jj,jk)111 avmv(ji,jj,jk) = ( avmv(ji,jj,jk) + qbvv(ji,jj,jk) ) * vmask(ji,jj,jk)112 avt (ji,jj,jk) = ( avt (ji,jj,jk) + qbv (ji,jj,jk) ) * tmask(ji,jj,jk)113 END DO114 END DO115 END DO116 !117 END SUBROUTINE zdf_qiao118 119 120 INTEGER FUNCTION zdf_qiao_alloc()121 !!----------------------------------------------------------------------122 !! *** FUNCTION zdf_qiao_alloc ***123 !!----------------------------------------------------------------------124 ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk), &125 & STAT = zdf_qiao_alloc )126 !127 IF( lk_mpp ) CALL mpp_sum ( zdf_qiao_alloc )128 IF( zdf_qiao_alloc > 0 ) CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')129 !130 END FUNCTION zdf_qiao_alloc131 132 94 !!====================================================================== 133 END MODULE zdf qiao95 END MODULE zdfswm -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7953 r7990 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 29 !! 4.0 ! 2017-04 (G. Madec) Remove CPP keys29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 30 30 !!---------------------------------------------------------------------- 31 31 … … 159 159 !! 160 160 !! ** Action : compute en (now turbulent kinetic energy) 161 !! update avt, avm u, avmv(before vertical eddy coef.)161 !! update avt, avm (before vertical eddy coef.) 162 162 !! 163 163 !! References : Gaspar et al., JGR, 1990, … … 175 175 #endif 176 176 ! 177 IF( kt /= nit000 ) THEN ! restore before value to compute tke 178 avt (:,:,:) = avt_k (:,:,:) 179 avm (:,:,:) = avm_k (:,:,:) 180 avmu(:,:,:) = avmu_k(:,:,:) 181 avmv(:,:,:) = avmv_k(:,:,:) 182 ENDIF 183 ! 184 CALL tke_tke ! now tke (en) 185 ! 186 CALL tke_avn ! now avt, avm, avmu, avmv 187 ! 188 avt_k (:,:,:) = avt (:,:,:) 189 avm_k (:,:,:) = avm (:,:,:) 190 avmu_k(:,:,:) = avmu(:,:,:) 191 avmv_k(:,:,:) = avmv(:,:,:) 177 avt(:,:,:) = avt_k(:,:,:) ! restore before Kz computed with TKE alone 178 avm(:,:,:) = avm_k(:,:,:) 179 ! 180 CALL tke_tke ! now tke (en) 181 ! 182 CALL tke_avn ! now avt, avm, dissl 183 ! 184 avt_k(:,:,:) = avt(:,:,:) ! store Kz computed with TKE alone 185 avm_k(:,:,:) = avm(:,:,:) 192 186 ! 193 187 #if defined key_agrif … … 195 189 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 196 190 #endif 197 ! 191 ! 192 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' ) 193 ! 198 194 END SUBROUTINE zdf_tke 199 195 … … 207 203 !! ** Method : - TKE surface boundary condition 208 204 !! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T) 209 !! - source term due to shear ( saved in avmu, avmv arrays)205 !! - source term due to shear (= Kz dz[Ub] * dz[Un] ) 210 206 !! - Now TKE : resolution of the TKE equation by inverting 211 207 !! a tridiagonal linear system by a "methode de chasse" … … 213 209 !! 214 210 !! ** Action : - en : now turbulent kinetic energy) 215 !! - avmu, avmv : production of TKE by shear at u and v-points216 !! (= Kz dz[Ub] * dz[Un] )217 211 !! --------------------------------------------------------------------- 218 212 INTEGER :: ji, jj, jk ! dummy loop arguments … … 228 222 REAL(wp) :: zzd_up, zzd_lw ! - - 229 223 !!bfr REAL(wp) :: zebot ! - - 230 INTEGER , POINTER, DIMENSION(:,:) :: imlc224 INTEGER , DIMENSION(jpi,jpj) :: imlc 231 225 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 232 226 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw, z3du, z3dv … … 236 230 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 237 231 ! 238 CALL wrk_alloc( jpi,jpj, imlc ) ! integer239 232 CALL wrk_alloc( jpi,jpj, zhlc ) 240 233 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) … … 473 466 END DO 474 467 ENDIF 468 !!gm not sure we need this lbc .... 475 469 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 476 470 ! 477 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer478 471 CALL wrk_dealloc( jpi,jpj, zhlc ) 479 472 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) … … 516 509 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 517 510 !! 518 !! ** Action : - avt : now vertical eddy diffusivity (w-point) 519 !! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points 511 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 520 512 !!---------------------------------------------------------------------- 521 513 INTEGER :: ji, jj, jk ! dummy loop indices 522 REAL(wp) :: zrn2, zraug, zcoef, zav 523 REAL(wp) :: zdku, zri, zsqen! - -524 REAL(wp) :: z dkv, zemxl, zemlm, zemlp! - -514 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 515 REAL(wp) :: zdku, zdkv, zsqen ! - - 516 REAL(wp) :: zemxl, zemlm, zemlp ! - - 525 517 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld 526 518 !!-------------------------------------------------------------------- … … 645 637 646 638 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 647 ! ! Vertical eddy viscosity and diffusivity (avm u, avmv,avt)639 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 648 640 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 649 641 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points … … 658 650 END DO 659 651 END DO 660 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 661 ! 662 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 663 DO jj = 2, jpjm1 664 DO ji = fs_2, fs_jpim1 ! vector opt. 665 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 666 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 667 END DO 668 END DO 669 END DO 670 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 652 ! 671 653 ! 672 654 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt … … 677 659 # if defined key_c1d 678 660 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 679 !!gm bug NO zri here....680 !!gm remove the specific diag for c1d !681 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri682 661 # endif 683 662 END DO … … 685 664 END DO 686 665 ENDIF 687 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 666 !!gm I'm sure this is needed to compute the shear term at next time-step 667 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 668 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 688 669 689 670 IF(ln_ctl) THEN 690 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 691 CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, & 692 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 671 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 672 CALL prt_ctl( tab3d_1=avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 693 673 ENDIF 694 674 ! … … 759 739 ENDIF 760 740 ! 761 IF( ln_zdftmx ) THEN ! Internal wave driven mixing 762 ! ! specific values of rn_emin & rmxl_min are used 763 rn_emin = 1.e-10_wp 764 rmxl_min = 1.e-03_wp 741 IF( ln_zdfiwm ) THEN ! Internal wave-driven mixing 742 rn_emin = 1.e-10_wp ! specific values of rn_emin & rmxl_min are used 743 rmxl_min = 1.e-03_wp ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s) 765 744 IF(lwp) WRITE(numout,*) ' Internal wave-driven mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 766 ELSE 745 ELSE ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s) 767 746 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity 768 747 IF(lwp) WRITE(numout,*) ' minimum mixing length with your parameters rmxl_min = ', rmxl_min … … 798 777 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 799 778 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 800 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)801 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)802 779 END DO 803 780 dissl(:,:,:) = 1.e-12_wp … … 821 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 822 799 ! 823 INTEGER :: jit, jk ! dummy loop indices824 INTEGER :: id1, id2, id3, id4 , id5, id6! local integers800 INTEGER :: jit, jk ! dummy loop indices 801 INTEGER :: id1, id2, id3, id4 ! local integers 825 802 !!---------------------------------------------------------------------- 826 803 ! … … 829 806 IF( ln_rstart ) THEN !* Read the restart file 830 807 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 831 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 832 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 833 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 834 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 835 id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 808 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 809 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 810 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 836 811 ! 837 812 IF( id1 > 0 ) THEN ! 'en' exists 838 813 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 839 IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist 840 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 841 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 842 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 843 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 814 IF( MIN( id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 815 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 816 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 844 817 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 845 818 ELSE ! one at least array is missing 846 CALL tke_avn ! compute avt, avm, a vmu, avmv and dissl (approximation)819 CALL tke_avn ! compute avt, avm, and dissl (approximation) 847 820 ENDIF 848 821 ELSE ! No TKE array found: initialisation 849 822 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 850 en (:,:,:) = rn_emin * tmask(:,:,:) 851 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 852 ! 853 avt_k (:,:,:) = avt (:,:,:) 854 avm_k (:,:,:) = avm (:,:,:) 855 avmu_k(:,:,:) = avmu(:,:,:) 856 avmv_k(:,:,:) = avmv(:,:,:) 823 en (:,:,:) = rn_emin * wmask(:,:,:) 824 CALL tke_avn ! recompute avt, avm, and dissl (approximation) 825 avt_k(:,:,:) = avt(:,:,:) 826 avm_k(:,:,:) = avm(:,:,:) 857 827 ! 858 828 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 829 avt_k(:,:,:) = avt(:,:,:) 830 avm_k(:,:,:) = avm(:,:,:) 859 831 ENDIF 860 832 ELSE !* Start from rest 861 833 en(:,:,:) = rn_emin * tmask(:,:,:) 862 834 DO jk = 1, jpk ! set the Kz to the background value 863 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 864 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 865 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 866 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 835 avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk) 836 avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk) 867 837 END DO 868 838 ENDIF … … 871 841 ! ! ------------------- 872 842 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 873 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 874 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 875 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 876 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 877 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 878 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 843 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 844 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 845 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 846 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 879 847 ! 880 848 ENDIF -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/module_example
r4147 r7990 86 86 !! Give references if exist otherwise suppress these lines 87 87 !!---------------------------------------------------------------------- 88 USE toto_module ! description of the module89 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released90 USE wrk_nemo, ONLY: zztab => wrk_2d_5 ! 2D workspace91 USE wrk_nemo, ONLY: zwx => wrk_3d_12 , zwy => wrk_3d_13 ! 3D workspace92 !!93 88 INTEGER , INTENT(in ) :: kt ! short description 94 89 INTEGER , INTENT(inout) :: pvar1 ! - - … … 100 95 REAL(wp) :: zmlmin, zbbrau ! temporary scalars (DOCTOR : start with z) 101 96 REAL(wp) :: zfact1, zfact2 ! do not use continuation lines in declaration 97 REAL(wp), DIMENSION(jpi,jpj) :: zwrk_2d ! 2D workspace 102 98 !!-------------------------------------------------------------------- 103 104 IF( wrk_in_use(3, 12,13) .OR. wrk_in_use(2, 5 ) THEN 105 CALL ctl_stop('exa_mpl: requested workspace arrays unavailable') ; RETURN 106 ENDIF 107 99 ! 108 100 IF( kt == nit000 ) CALL exa_mpl_init ! Initialization (first time-step only) 109 101 … … 119 111 DO jj = 2, jpjm1 120 112 DO ji = fs_2, fs_jpim1 ! vector opt. 121 avm v(ji,jj,jk) = ....113 avm(ji,jj,jk) = .... 122 114 END DO 123 115 END DO … … 128 120 DO jj = 2, jpjm1 129 121 DO ji = fs_2, fs_jpim1 ! vector opt. 130 avm v(ji,jj,jk) = ...122 avm(ji,jj,jk) = ... 131 123 END DO 132 124 END DO … … 135 127 END SELECT 136 128 ! 137 CALL mpplnk2( avmu, 'U', 1. ) ! Lateral boundary conditions (unchanged sign) 138 ! 139 IF( wrk_not_released(3, 12,13) .OR. wrk_not_released(2, 5 ) THEN 140 CALL ctl_stop('exa_mpl: failed to release workspace arrays') ; RETURN 141 ENDIF 129 CALL lbc_lnk( avm, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 142 130 ! 143 131 END SUBROUTINE exa_mpl -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r7953 r7990 64 64 65 65 USE zdfphy ! vertical physics manager (zdf_phy_init routine) 66 !!gm to be suppressed67 USE zdftmx ! tide-induced vertical mixing (zdf_tmx routine)68 USE zdfbfr ! bottom friction (zdf_bfr routine)69 USE zdftke ! TKE vertical mixing (zdf_tke routine)70 USE zdfgls ! GLS vertical mixing (zdf_gls routine)71 USE zdfddm ! double diffusion mixing (zdf_ddm routine)72 USE zdfevd ! enhanced vertical diffusion (zdf_evd routine)73 USE zdfric ! Richardson vertical mixing (zdf_ric routine)74 USE zdfmxl ! Mixed-layer depth (zdf_mxl routine)75 USE zdfqiao !Qiao module wave induced mixing (zdf_qiao routine)76 !!gm end77 66 78 67 USE step_diu ! Time stepping for diurnal sst -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90
r7681 r7990 119 119 ! ! -------------- 120 120 neln(:,:) = 1 ! euphotic layer level 121 DO jk = 1, jpk 121 DO jk = 1, jpkm1 ! (i.e. 1rst T-level strictly below EL bottom) 122 122 DO jj = 1, jpj 123 123 DO ji = 1, jpi … … 147 147 END SUBROUTINE p2z_opt 148 148 149 149 150 SUBROUTINE p2z_opt_init 150 151 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
