Changeset 7990


Ignore:
Timestamp:
2017-04-29T17:24:54+02:00 (3 years ago)
Author:
gm
Message:

#1880 (HPC-09): OPA remove avmu, avmv from zdf modules + move CALL tke(gls)_rst & gls_rst in zdftke(gls) + rename zdftmx and zdfqiao

Location:
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM
Files:
27 edited
2 moved

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg

    r7954 r7990  
    324324      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    325325   ! 
    326    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    327    ! 
    328    ln_zdfqiao  = .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 ) 
    329329   ! 
    330330   !                       ! time-stepping 
     
    353353/ 
    354354!----------------------------------------------------------------------- 
    355 &namzdf_tmx    !   internal wave-driven mixing parameterization         (ln_zdftmx =T) 
     355&namzdf_iwm    !   internal wave-driven mixing parameterization         (ln_zdfiwm =T) 
    356356!----------------------------------------------------------------------- 
    357357/ 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg

    r7954 r7990  
    275275      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    276276   ! 
    277    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    278    ! 
    279    ln_zdfqiao  = .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 ) 
    280280   ! 
    281281   !                       ! time-stepping 
    282282   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
     283      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    283284   ! 
    284285   !                       ! coefficients 
     
    301302/ 
    302303!----------------------------------------------------------------------- 
    303 &namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T) 
     304&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
    304305!----------------------------------------------------------------------- 
    305306/ 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg

    r7954 r7990  
    264264      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    265265   ! 
    266    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    267    ! 
    268    ln_zdfqiao  = .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 ) 
    269269   ! 
    270270   !                       ! time-stepping 
     
    292292/ 
    293293!----------------------------------------------------------------------- 
    294 &namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T) 
     294&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
    295295!----------------------------------------------------------------------- 
    296296/ 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg

    r7954 r7990  
    218218      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    219219   ! 
    220    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    221    ! 
    222    ln_zdfqiao  = .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 ) 
    223223   ! 
    224224   !                       ! time-stepping 
    225    ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme 
    226       nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T 
     225   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
     226      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    227227   ! 
    228228   !                       !  Coefficients 
     
    251251/ 
    252252!----------------------------------------------------------------------- 
    253 &namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx=T) 
    254 !----------------------------------------------------------------------- 
    255    ln_tmx_itf  = .false.   !  ITF specific parameterisation 
    256 / 
    257 !----------------------------------------------------------------------- 
    258253&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi) 
    259254!----------------------------------------------------------------------- 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg

    r7954 r7990  
    226226      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    227227   ! 
    228    ln_zdftmx   = .true.    ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    229    ! 
    230    ln_zdfqiao  = .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 ) 
    231231   ! 
    232232   !                       ! time-stepping 
    233    ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme 
    234       nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T 
     233   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
     234      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    235235   ! 
    236236   !                       !  Coefficients 
     
    245245/ 
    246246!----------------------------------------------------------------------- 
    247 &namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T) 
     247&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
    248248!----------------------------------------------------------------------- 
    249249   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  
    2727/ 
    2828!----------------------------------------------------------------------- 
    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) 
    3130!----------------------------------------------------------------------- 
    3231/ 
     
    255254      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    256255   ! 
    257    ln_zdftmx   = .true.    ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    258    ! 
    259    ln_zdfqiao  = .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 ) 
    260259   ! 
    261260   !                       ! time-stepping 
    262    ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme 
    263       nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T 
     261   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
     262      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    264263   ! 
    265264   !                       !  Coefficients 
     
    274273/ 
    275274!----------------------------------------------------------------------- 
    276 &namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx =T) 
     275&namzdf_iwm    !   tidal mixing parameterization                        (ln_zdfiwm =T) 
    277276!----------------------------------------------------------------------- 
    278277   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  
    398398         <field id="avm_evd"      long_name="convective enhancement of vertical viscosity"     standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" /> 
    399399 
    400          <!-- variables available with ln_zdftmx =T --> 
     400         <!-- variables available with ln_zdfiwm =T --> 
    401401         <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    /> 
    402402         <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" /> 
    406406 
    407407         <!-- variables available with diaar5 -->    
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7954 r7990  
    1111!!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp) 
    1212!!              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) 
    1414!!              9 - miscellaneous    (nammpp, namctl) 
    1515!!             10 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto) 
     
    875875!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T) 
    876876!!    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) 
    878878!!====================================================================== 
    879879! 
     
    899899      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    900900   ! 
    901    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    902    ! 
    903    ln_zdfqiao  = .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 ) 
    904904   ! 
    905905   !                       ! time-stepping 
     
    919919   rn_alp      =    5.     !  coefficient of the parameterization 
    920920   nn_ric      =    2      !  coefficient of the parameterization 
    921    rn_ekmfc    =    0.7    !  Factor in the Ekman depth Equation 
    922    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-layer 
    925    rn_wvmix    =   10.0    !  vertical eddy diffusion coeff [m2/s] in the mixed-layer 
    926    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 
    927927/ 
    928928!----------------------------------------------------------------------- 
     
    974974/ 
    975975!----------------------------------------------------------------------- 
    976 &namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T) 
     976&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T) 
    977977!----------------------------------------------------------------------- 
    978978   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  
    344344   ln_zdfddm   = .false.   ! double diffusive mixing 
    345345   ! 
    346    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    347    ! 
    348    ln_zdfqiao  = .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 ) 
    349349   ! 
    350350   !                       ! time-stepping 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg

    r7954 r7990  
    223223/ 
    224224!----------------------------------------------------------------------- 
    225 &namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx =T) 
     225&namzdf_iwm    !   tidal mixing parameterization                        (ln_zdfiwm =T) 
    226226!----------------------------------------------------------------------- 
    227227/ 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg

    r7954 r7990  
    328328!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T) 
    329329!!    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) 
    331331!!====================================================================== 
    332332!----------------------------------------------------------------------- 
     
    351351      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio 
    352352   ! 
    353    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx) 
    354    ! 
    355    ln_zdfqiao  = .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 ) 
    356356   ! 
    357357   !                       ! 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 
    359360   ! 
    360361   !                       ! coefficients 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt

    r7954 r7990  
    66ORCA2_OFF_TRC OPA_SRC OFF_SRC TOP_SRC 
    77ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC 
     8GYRE_PISCES OPA_SRC TOP_SRC 
    89GYRE_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  
    99   USE dom_oce         ! ocean space and time domain 
    1010   USE zdf_oce         ! ocean vertical physics     
    11    USE zdfgls, ONLY: mxln 
     11   USE zdfgls   , ONLY : hmxn 
    1212   USE in_out_manager  ! I/O units 
    1313   USE iom             ! I/0 library 
     
    107107      IF( ln_zdfgls ) THEN 
    108108         en_25h(:,:,:) = en(:,:,:) 
    109          rmxln_25h(:,:,:) = mxln(:,:,:) 
     109         rmxln_25h(:,:,:) = hmxn(:,:,:) 
    110110      ENDIF 
    111111#if defined key_lim3 || defined key_lim2 
     
    169169         IF( ln_zdfgls ) THEN 
    170170            en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
    171             rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:) 
     171            rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + hmxn(:,:,:) 
    172172         ENDIF 
    173173         cnt_25h = cnt_25h + 1 
     
    217217         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
    218218         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(:,:,:)) 
    220220         CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
    221221         ! Write vertical physics 
    222          zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     222         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    223223         CALL iom_put("avt25h", zw3d )   ! diffusivity 
    224          zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     224         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    225225         CALL iom_put("avm25h", zw3d)   ! viscosity 
    226226         IF( ln_zdftke ) THEN 
    227             zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     227            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    228228            CALL iom_put("tke25h", zw3d)   ! tke 
    229229         ENDIF 
    230230         IF( ln_zdfgls ) THEN 
    231             zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     231            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    232232            CALL iom_put("tke25h", zw3d)   ! tke 
    233             zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     233            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    234234            CALL iom_put( "mxln25h",zw3d) 
    235235         ENDIF 
     
    249249         IF( ln_zdfgls ) THEN 
    250250            en_25h(:,:,:) = en(:,:,:) 
    251             rmxln_25h(:,:,:) = mxln(:,:,:) 
     251            rmxln_25h(:,:,:) = hmxn(:,:,:) 
    252252         ENDIF 
    253253         cnt_25h = 1 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r7953 r7990  
    4242   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
    44  
    4544CONTAINS 
    4645    
     
    5150      !! ** Purpose :   compute the vertical ocean dynamics physics. 
    5251      !!--------------------------------------------------------------------- 
    53       !! 
    5452      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    5553      ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r7753 r7990  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit scheme 
     14   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit time-stepping 
    1515   !!                   together with the Leap-Frog time integration. 
    1616   !!---------------------------------------------------------------------- 
     
    3939#  include "vectopt_loop_substitute.h90" 
    4040   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     41   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    4242   !! $Id$ 
    4343   !! 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  
    99   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields 
    1010   !!---------------------------------------------------------------------- 
     11 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   namsbc_cpl      : coupled formulation namlist 
     
    974975      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    975976      !!---------------------------------------------------------------------- 
    976       USE zdf_oce,  ONLY : ln_zdfqiao 
     977      USE zdf_oce,  ONLY : ln_zdfswm 
    977978 
    978979      IMPLICIT NONE 
     
    11591160      !                                                      !      Wave mean period     ! 
    11601161      !                                                      ! ========================= !  
    1161          IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1) 
     1162         IF( srcv(jpr_wper)%laction )   wmp(:,:) = frcv(jpr_wper)%z3(:,:,1) 
    11621163      ! 
    11631164      !                                                      ! ========================= !  
    11641165      !                                                      !  Significant wave height  ! 
    11651166      !                                                      ! ========================= !  
    1166          IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
     1167         IF( srcv(jpr_hsig)%laction )   hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
    11671168      ! 
    11681169      !                                                      ! ========================= !  
    1169       !                                                      !    Vertical mixing Qiao   ! 
     1170      !                                                      !    surface wave mixing    ! 
    11701171      !                                                      ! ========================= !  
    1171          IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1) 
     1172         IF( srcv(jpr_wnum)%laction .AND. ln_zdfswm )  wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1) 
    11721173 
    11731174         ! 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  
    1919   USE oce            ! ocean variables 
    2020   USE sbc_oce        ! Surface boundary condition: ocean fields 
    21    USE zdf_oce,  ONLY : ln_zdfqiao 
     21   USE zdf_oce,  ONLY : ln_zdfswm 
    2222   USE bdy_oce        ! open boundary condition variables 
    2323   USE domvvl         ! domain: variable volume layers 
     
    227227         ! 
    228228         ! Read also wave number if needed, so that it is available in coupling routines 
    229          IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     229         IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN 
    230230            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing 
    231231            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     
    345345         vsd(:,:,:) = 0._wp 
    346346         wsd(:,:,:) = 0._wp 
    347          ! Wave number needed only if ln_zdfqiao=T 
     347         ! Wave number needed only if ln_zdfswm=T 
    348348         IF( .NOT. cpl_wnum ) THEN 
    349349            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  
    1616   PUBLIC  zdf_oce_alloc    ! Called in nemogcm.F90 
    1717 
    18    !                                 !!* namelist namzdf: vertical diffusion * 
    19 !!gm 
     18   !                            !!* namelist namzdf: vertical physics * 
    2019   !                             ! vertical closure scheme flags 
    2120   LOGICAL , PUBLIC ::   ln_zdfcst   !: constant coefficients 
     
    2322   LOGICAL , PUBLIC ::   ln_zdftke   !: Turbulent Kinetic Energy closure 
    2423   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 
    3425   LOGICAL , PUBLIC ::   ln_zdfevd   !: convection: enhanced vertical diffusion flag 
    3526   INTEGER , PUBLIC ::      nn_evdm     !: =0/1 flag to apply enhanced avm or not 
     
    3829   INTEGER , PUBLIC ::      nn_npc      !: non penetrative convective scheme call  frequency 
    3930   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) 
    4241   !                             ! coefficients  
    4342   REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s) 
     
    4746 
    4847 
    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] 
    5153   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt 
    5254   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 
    6057 
    6158   !!---------------------------------------------------------------------- 
     
    7168      !!---------------------------------------------------------------------- 
    7269      ! 
    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 ) 
    8175         ! 
    8276      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  
    88   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    99   !!            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  
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2223   USE prtctl         ! Print control 
    2324   USE lib_mpp        ! MPP library 
    24    USE wrk_nemo       ! work arrays 
    2525   USE timing         ! Timing 
    2626 
     
    6363      !!      avt = avt + zavft + zavdt 
    6464      !!      avs = avs + zavfs + zavds 
    65       !!      avmu, avmv are required to remain at least above avt and avs. 
     65      !!      avm is required to remain at least above avt and avs. 
    6666      !!       
    6767      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S 
     
    7777      REAL(wp) ::   zavft, zavfs    !   -      - 
    7878      REAL(wp) ::   zavdt, zavds    !   -      - 
    79       REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
     79      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
    8080      !!---------------------------------------------------------------------- 
    8181      ! 
    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') 
    8583      ! 
    8684      !                                                ! =============== 
     
    8987         ! Define the mask  
    9088         ! --------------- 
     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          
    9193         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
    9294            DO ji = 1, jpi 
    9395               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
     96!!gm please, use e3w_n below  
    9497                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )  
    9598               ! 
     
    156159            END DO 
    157160         END DO 
    158  
    159  
    160          ! Increase avmu, avmv if necessary 
    161          ! -------------------------------- 
    162 !!gm to be changed following the definition of avm. 
    163          DO jj = 1, jpjm1 
    164             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 DO 
    172          END DO 
    173161         !                                                ! =============== 
    174162      END DO                                              !   End of slab 
    175163      !                                                   ! =============== 
    176164      ! 
    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  
    183165      IF(ln_ctl) THEN 
    184166         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) 
    187167      ENDIF 
    188       ! 
    189       CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
    190168      ! 
    191169      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm') 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r7931 r7990  
    88   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    99   !!            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)  
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2324   USE iom             ! for iom_put 
    2425   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    25    USE wrk_nemo        ! work arrays 
    2626   USE timing          ! Timing 
    2727 
     
    4545      !!      sivity coefficients when a static instability is encountered. 
    4646      !! 
    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. 
    5054      !! 
    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 
    5456      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! ocean time-step indexocean time step 
     57      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step 
    5658      ! 
    5759      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    58       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zavt_evd, zavm_evd 
     60      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd 
    5961      !!---------------------------------------------------------------------- 
    6062      ! 
     
    6870      ENDIF 
    6971      ! 
    70       CALL wrk_alloc( jpi,jpj,jpk,   zavt_evd, zavm_evd )  
    7172      ! 
    7273      zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
     
    7475      SELECT CASE ( nn_evdm ) 
    7576      ! 
    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) 
    7778         ! 
    7879         zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application 
    7980         ! 
     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         ! 
    8087         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 
    8390                  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) 
    9093                  ENDIF 
    9194               END DO 
    9295            END DO 
    9396         END DO  
    94          CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions 
    95          CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. ) 
    9697         ! 
    9798         zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd 
    9899         CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
    99100         ! 
    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 
    101107         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 
    105110                  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) 
    107112               END DO 
    108113            END DO 
     
    110115         ! 
    111116      END SELECT  
    112  
     117      ! 
    113118      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
    114119      CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
    115120      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 )  
    118121      ! 
    119122      IF( nn_timing == 1 )  CALL timing_stop('zdf_evd') 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7953 r7990  
    55   !!                 turbulent closure parameterization 
    66   !!====================================================================== 
    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  
    910   !!---------------------------------------------------------------------- 
    1011 
     
    3637   PRIVATE 
    3738 
    38    PUBLIC   zdf_gls        ! routine called in step module 
    39    PUBLIC   zdf_gls_init   ! routine called in zdfphy module 
    40    PUBLIC   gls_rst        ! routine called in zdfphy module 
     39   PUBLIC   zdf_gls        ! called in zdfphy 
     40   PUBLIC   zdf_gls_init   ! called in zdfphy 
     41   PUBLIC   gls_rst        ! called in zdfphy 
    4142 
    4243   ! 
    43    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length 
     44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxn    !: now mixing length 
    4445   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function 
    4546   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
     
    113114      !!                ***  FUNCTION zdf_gls_alloc  *** 
    114115      !!---------------------------------------------------------------------- 
    115       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
     116      ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    116117         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
    117118         ! 
     
    148149      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal 
    149150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal 
     151      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3du, z3dv  ! u- and v-shears 
    150152      !!-------------------------------------------------------------------- 
    151153      ! 
    152       IF( nn_timing == 1 )  CALL timing_start('zdf_gls') 
     154      IF( nn_timing == 1 )   CALL timing_start('zdf_gls') 
    153155      ! 
    154156      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 
    157160      ! Preliminary computing 
    158161 
    159162      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    160163 
    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(:,:,:) 
    167167 
    168168      ! Compute surface and bottom friction at T-points 
     
    183183      END DO     
    184184 
    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           
    189187         zhsro(:,:) = rn_hsro 
    190188      CASE ( 1 )             ! Standard Charnock formula 
    191189         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 
    192190      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) 
    193193         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    194194         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    195195      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
     196!!gm  BUG missing a multiplicative coefficient.... 
    196197         zhsro(:,:) = hsw(:,:) 
    197198      END SELECT 
    198199 
    199       ! Compute shear and dissipation rate 
    200       DO jk = 2, jpkm1 
    201          DO jj = 2, jpjm1 
    202             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) 
    210211            END DO 
    211212         END DO 
    212213      END DO 
    213       ! 
    214       ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 
    215       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. ) 
    216214 
    217215      ! Save tke at before time step 
    218216      eb  (:,:,:) = en  (:,:,:) 
    219       mxlb(:,:,:) = mxln(:,:,:) 
     217      mxlb(:,:,:) = hmxn(:,:,:) 
    220218 
    221219      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
     
    223221            DO jj = 2, jpjm1  
    224222               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) 
    226224                  zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 
    227225                  zcoef = ( zup / MAX( zdown, rsmall ) ) 
     
    247245      DO jk = 2, jpkm1 
    248246         DO jj = 2, jpjm1 
    249             DO ji = fs_2, fs_jpim1   ! vector opt. 
     247            DO ji = 2, jpim1 
    250248               ! 
    251249               ! 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) ) 
    254252               ! 
    255253               ! stratif. destruction 
     
    278276               ! building the matrix 
    279277               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) 
    295286            END DO 
    296287         END DO 
     
    300291      ! 
    301292      ! 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 
    304295      ! 
    305296      ! Surface boundary condition on tke 
     
    353344      ! 
    354345      CASE ( 0 )             ! Dirichlet  
    355          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 
     346         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin 
    356347         !                      ! Balance between the production and the dissipation terms 
    357348         DO jj = 2, jpjm1 
     
    503494               ! building the matrix 
    504495               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) 
    517504            END DO 
    518505         END DO 
     
    565552      zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    566553      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
    567  
    568554      !    
    569555      ! 
     
    577563      ! 
    578564      CASE ( 0 )             ! Dirichlet  
    579          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 
     565         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0 
    580566         !                      ! Balance between the production and the dissipation terms 
    581567         DO jj = 2, jpjm1 
     
    700686      ! Limit dissipation rate under stable stratification 
    701687      ! -------------------------------------------------- 
    702       DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time 
     688      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time 
    703689         DO jj = 2, jpjm1 
    704690            DO ji = fs_2, fs_jpim1    ! vector opt. 
    705691               ! limitation 
    706692               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) 
    708694               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
    709695               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) ) 
    711697            END DO 
    712698         END DO 
     
    733719                  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) 
    734720                  ! 
    735                   ! Store stability function in avmu and avmv 
    736                   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) 
    738724               END DO 
    739725            END DO 
     
    761747                  sh = (rs4 - rs5*gh + rs6*gm) / rcff 
    762748                  ! 
    763                   ! Store stability function in avmu and avmv 
    764                   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) 
    766752               END DO 
    767753            END DO 
     
    773759      ! Lines below are useless if GOTM style Dirichlet conditions are used 
    774760 
    775       avmv(:,:,1) = avmv(:,:,2) 
     761      z3dv(:,:,1) = z3dv(:,:,2) 
    776762 
    777763      DO jj = 2, jpjm1 
    778764         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)) 
    780766         END DO 
    781767      END DO 
     
    786772         DO jj = 2, jpjm1 
    787773            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 routine 
    791                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 bottom 
     774               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 
    793779            END DO 
    794780         END DO 
    795781      END DO 
    796       ! 
    797       ! Lateral boundary conditions (sign unchanged) 
    798782      avt(:,:,1)  = 0._wp 
     783!!gm I'm sure this is needed to compute the shear term at next time-step 
    799784      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      ! 
    812786      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 ) 
    816789      ENDIF 
    817790      ! 
    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 
    822795      ! 
    823796      CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro ) 
    824797      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') 
    828801      ! 
    829802   END SUBROUTINE zdf_gls 
     
    835808      !!                      
    836809      !! ** Purpose :   Initialization of the vertical eddy diffivity and  
    837       !!      viscosity when using a gls turbulent closure scheme 
     810      !!              viscosity computed using a GLS turbulent closure scheme 
    838811      !! 
    839812      !! ** Method  :   Read the namzdf_gls namelist and check the parameters 
    840       !!      called at the first timestep (nit000) 
    841813      !! 
    842814      !! ** input   :   Namlist namzdf_gls 
     
    892864      ENDIF 
    893865 
    894       !                                !* allocate gls arrays 
     866      !                                !* allocate GLS arrays 
    895867      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' ) 
    896868 
     
    11201092      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
    11211093      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
    1122  
     1094      ! 
    11231095      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
    11241096      rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke 
    1125  
     1097      ! 
    11261098      !                                !* Wall proximity function 
    11271099      zwall (:,:,:) = 1._wp * tmask(:,:,:) 
    11281100 
    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) 
    11381103      ! 
    11391104      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init') 
     
    11521117      !!                set to rn_emin or recomputed (nn_igls/=0) 
    11531118      !!---------------------------------------------------------------------- 
    1154       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1155       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
     1119      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     1120      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    11561121      ! 
    11571122      INTEGER ::   jit, jk   ! dummy loop indices 
    1158       INTEGER ::   id1, id2, id3, id4, id5, id6 
     1123      INTEGER ::   id1, id2, id3, id4 
    11591124      INTEGER ::   ji, jj, ikbu, ikbv 
    11601125      REAL(wp)::   cbx, cby 
     
    11651130         IF( ln_rstart ) THEN                   !* Read the restart file 
    11661131            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. ) 
    11721135            ! 
    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  ) 
    11801141            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(:,:,:) 
    11881147               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO 
    11891148            ENDIF 
    11901149         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' 
    11921151            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 
    11941157         ENDIF 
    11951158         ! 
     
    11971160         !                                   ! ------------------- 
    11981161         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  ) 
    12051166         ! 
    12061167      ENDIF 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

    r7967 r7990  
    1 MODULE zdftmx 
     1MODULE zdfiwm 
    22   !!======================================================================== 
    3    !!                       ***  MODULE  zdftmx  *** 
     3   !!                       ***  MODULE  zdfiwm  *** 
    44   !! Ocean physics: Internal gravity wave-driven vertical mixing 
    55   !!======================================================================== 
     
    88   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase 
    99   !!            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 
    1111   !!---------------------------------------------------------------------- 
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz 
    15    !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz 
     14   !!   zdf_iwm       : global     momentum & tracer Kz with wave induced Kz 
     15   !!   zdf_iwm_init  : global     momentum & tracer Kz with wave induced Kz 
    1616   !!---------------------------------------------------------------------- 
    1717   USE oce            ! ocean dynamics and tracers variables 
     
    3333   PRIVATE 
    3434 
    35    PUBLIC   zdf_tmx         ! called in step module  
    36    PUBLIC   zdf_tmx_init    ! called in nemogcm module  
    37    PUBLIC   zdf_tmx_alloc   ! called in nemogcm module 
    38  
    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 * 
    4040   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
    4141   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
     
    4444   REAL(wp) ::  r1_6 = 1._wp / 6._wp 
    4545 
    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) 
    5454   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T) 
    5555   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity 
     
    6464CONTAINS 
    6565 
    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_alloc 
    78  
    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  *** 
    8383      !!                    
    8484      !! ** Purpose :   add to the vertical mixing coefficients the effect of 
     
    8686      !! 
    8787      !! ** 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_tmx is the 3D space distribution of the wave-breaking  
     88      !!                  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  
    9090      !!              energy and Nu the molecular kinematic viscosity. 
    9191      !!              The function f(Reb) is linear (constant mixing efficiency) 
    9292      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 
    9393      !! 
    94       !!              - Compute emix_tmx, the 3D power density that allows to compute 
     94      !!              - Compute emix_iwm, the 3D power density that allows to compute 
    9595      !!              Reb and therefrom the wave-induced vertical diffusivity. 
    9696      !!              This is divided into three components: 
    9797      !!                 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_tmx 
    100       !!              where hcri_tmx is the characteristic length scale of the bottom  
    101       !!              intensification, ecri_tmx a 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. 
    102102      !!                 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 ) 
    104104      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
    105       !!              where epyc_tmx is a map of available power, and nn_zpyc 
     105      !!              where epyc_iwm is a map of available power, and nn_zpyc 
    106106      !!              is the chosen stratification-dependence of the internal wave 
    107107      !!              energy dissipation. 
    108108      !!                 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_tmx is the characteristic length scale of the WKB bottom  
    112       !!              intensification, ebot_tmx is a map of available power, and z_wkb is the 
     109      !!                     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 
    113113      !!              WKB-stretched height above bottom defined as 
    114114      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 
     
    118118      !!                     avt  = avt  +    av_wave 
    119119      !!                     avm  = avm  +    av_wave 
    120       !!                     avmu = avmu + mi(av_wave) 
    121       !!                     avmv = avmv + mj(av_wave) 
    122120      !! 
    123121      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing: 
    124122      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb) 
    125123      !! 
    126       !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing 
    127       !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing     
     124      !! ** Action  : - Define emix_iwm used to compute internal wave-induced mixing 
     125      !!              - avt, avs, avm, increased by internal wave-driven mixing     
    128126      !! 
    129127      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
     
    142140      !!---------------------------------------------------------------------- 
    143141      ! 
    144       IF( nn_timing == 1 )   CALL timing_start('zdf_tmx') 
     142      IF( nn_timing == 1 )   CALL timing_start('zdf_iwm') 
    145143      ! 
    146144      CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
     
    156154         DO ji = 1, jpi 
    157155            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) 
    160158         END DO 
    161159      END DO 
    162160 
    163161      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)   & 
    166164!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point 
    167165            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     
    170168      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
    171169      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
    172  
     170      ! 
    173171      SELECT CASE ( nn_zpyc ) 
    174  
     172      ! 
    175173      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    176  
     174         ! 
    177175         zfact(:,:) = 0._wp 
    178176         DO jk = 2, jpkm1              ! part independent of the level 
    179177            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    180178         END DO 
    181  
     179         ! 
    182180         DO jj = 1, jpj 
    183181            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         ! 
    188186         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         ! 
    192190      CASE ( 2 )               ! Dissipation scales as N^2 
    193  
     191         ! 
    194192         zfact(:,:) = 0._wp 
    195193         DO jk = 2, jpkm1              ! part independent of the level 
    196194            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    197195         END DO 
    198  
     196         ! 
    199197         DO jj= 1, jpj 
    200198            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         ! 
    205203         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         ! 
    209207      END SELECT 
    210208 
    211209      !                        !* WKB-height dependent mixing: distribute energy over the time-varying  
    212210      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    213        
     211      ! 
    214212      zwkb(:,:,:) = 0._wp 
    215213      zfact(:,:) = 0._wp 
     
    218216         zwkb(:,:,jk) = zfact(:,:) 
    219217      END DO 
    220  
     218      ! 
    221219      DO jk = 2, jpkm1 
    222220         DO jj = 1, jpj 
    223221            DO ji = 1, jpi 
    224222               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) 
    226224            END DO 
    227225         END DO 
    228226      END DO 
    229227      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
    230  
     228      ! 
    231229      zweight(:,:,:) = 0._wp 
    232230      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      ! 
    237235      zfact(:,:) = 0._wp 
    238236      DO jk = 2, jpkm1              ! part independent of the level 
    239237         zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    240238      END DO 
    241  
     239      ! 
    242240      DO jj = 1, jpj 
    243241         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      ! 
    248246      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)   & 
    250248            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    251249      END DO 
    252  
    253  
     250      ! 
    254251      ! Calculate molecular kinematic viscosity 
    255252      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
     
    258255         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    259256      END DO 
    260  
     257      ! 
    261258      ! Calculate turbulence intensity parameter Reb 
    262259      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      ! 
    266263      ! Define internal wave-induced diffusivity 
    267264      DO jk = 2, jpkm1 
    268265         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    269266      END DO 
    270  
     267      ! 
    271268      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    272269         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
     
    282279         END DO 
    283280      ENDIF 
    284  
     281      ! 
    285282      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    286283         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
    287284      END DO 
    288  
     285      ! 
    289286      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    290287         ztpc = 0._wp 
     
    300297         IF( lk_mpp )   CALL mpp_sum( ztpc ) 
    301298         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    302   
     299         ! 
    303300         IF(lwp) THEN 
    304301            WRITE(numout,*) 
    305             WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)' 
     302            WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)' 
    306303            WRITE(numout,*) '~~~~~~~ ' 
    307304            WRITE(numout,*) 
     
    339336      ENDIF 
    340337 
    341       DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points 
    342          DO jj = 2, jpjm1 
    343             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 DO 
    347          END DO 
    348       END DO 
    349       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition 
    350  
    351338      !                             !* output internal wave-driven mixing coefficient 
    352339      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") ) THEN 
    356          bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    357          pcmap_tmx(:,:) = 0._wp 
     340                                    !* 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 
    358345         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 ) 
    364351      ENDIF 
    365352      CALL iom_put( "bn2", rn2 ) 
    366       CALL iom_put( "emix_tmx", emix_tmx ) 
     353      CALL iom_put( "emix_iwm", emix_iwm ) 
    367354       
    368355      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
    369356      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    370357 
    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_tmx 
    376  
    377  
    378    SUBROUTINE zdf_tmx_init 
    379       !!---------------------------------------------------------------------- 
    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  *** 
    381368      !!                      
    382369      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading 
    383370      !!              of input power maps and decay length scales in netcdf files. 
    384371      !! 
    385       !! ** Method  : - Read the namzdf_tmx namelist and check the parameters 
     372      !! ** Method  : - Read the namzdf_iwm namelist and check the parameters 
    386373      !! 
    387374      !!              - Read the input data in NetCDF files : 
     
    392379      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc) 
    393380      !! 
    394       !! ** input   : - Namlist namzdf_tmx 
     381      !! ** input   : - Namlist namzdf_iwm 
    395382      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 
    396383      !!              decay_scale_bot.nc decay_scale_cri.nc 
    397384      !! 
    398385      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter 
    399       !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx 
    400       !! 
    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 
    403390      !!---------------------------------------------------------------------- 
    404391      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    407394      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
    408395      !! 
    409       NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff 
    410       !!---------------------------------------------------------------------- 
    411       ! 
    412       IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init') 
    413       ! 
    414       REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing 
    415       READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901) 
    416 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp ) 
    417       ! 
    418       REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing 
    419       READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 ) 
    420 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in 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) 
     403901   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 ) 
     407902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist', lwp ) 
     408      IF(lwm) WRITE ( numond, namzdf_iwm_new ) 
    422409      ! 
    423410      IF(lwp) THEN                  ! Control print 
    424411         WRITE(numout,*) 
    425          WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing' 
     412         WRITE(numout,*) 'zdf_iwm_init : internal wave-driven mixing' 
    426413         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' 
    428415         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 
    429416         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar 
     
    435422      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 
    436423      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)     
    438425      avtb_2d(:,:) = 1.e0_wp     ! uniform  
    439426      IF(lwp) THEN                  ! Control print 
     
    443430      ENDIF 
    444431             
    445       !                             ! allocate tmx arrays 
    446       IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' ) 
     432      !                             ! allocate iwm arrays 
     433      IF( zdf_iwm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' ) 
    447434      ! 
    448435      !                             ! read necessary fields 
    449436      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 )  
    451438      CALL iom_close(inum) 
    452439      ! 
    453440      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 ) 
    455442      CALL iom_close(inum) 
    456443      ! 
    457444      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 ) 
    459446      CALL iom_close(inum) 
    460447      ! 
    461448      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 ) 
    463450      CALL iom_close(inum) 
    464451      ! 
    465452      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 ) 
    467454      CALL iom_close(inum) 
    468455 
    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(:,:) 
    472459 
    473460      ! Set once for all to zero the first and last vertical levels of appropriate variables 
    474       emix_tmx (:,:, 1 ) = 0._wp 
    475       emix_tmx (:,:,jpk) = 0._wp 
     461      emix_iwm (:,:, 1 ) = 0._wp 
     462      emix_iwm (:,:,jpk) = 0._wp 
    476463      zav_ratio(:,:, 1 ) = 0._wp 
    477464      zav_ratio(:,:,jpk) = 0._wp 
     
    479466      zav_wave (:,:,jpk) = 0._wp 
    480467 
    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(:,:) ) 
    484471      IF(lwp) THEN 
    485472         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
     
    488475      ENDIF 
    489476      ! 
    490       IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init') 
    491       ! 
    492    END SUBROUTINE zdf_tmx_init 
     477      IF( nn_timing == 1 )  CALL timing_stop('zdf_iwm_init') 
     478      ! 
     479   END SUBROUTINE zdf_iwm_init 
    493480 
    494481   !!====================================================================== 
    495 END MODULE zdftmx 
     482END MODULE zdfiwm 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

    r7953 r7990  
    22   !!====================================================================== 
    33   !!                      ***  MODULE  zdfphy  *** 
    4    !! Ocean physics :   manager of vertical mixing parametrizations 
     4   !! Vertical ocean physics :   manager of all vertical physics packages 
    55   !!====================================================================== 
    66   !! History :  4.0  !  2017-04  (G. Madec)  original code 
     
    88 
    99   !!---------------------------------------------------------------------- 
    10    !!   zdf_phy_init  : initialization of all vertical physics pakages 
     10   !!   zdf_phy_init  : initialization of all vertical physics packages 
    1111   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff.  
    1212   !!---------------------------------------------------------------------- 
    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            
    1526   USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test) 
    16    USE zdfbfr         ! bottom friction 
    17    USE zdftke         ! TKE vertical mixing 
    18    USE zdfgls         ! GLS vertical mixing 
    19    USE zdfric         ! Richardson vertical mixing    
    20    USE zdfddm         ! double diffusion mixing       
    21    USE zdfevd         ! enhanced vertical diffusion 
    22    USE zdftmx         ! internal tide-induced mixing 
    23    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 adjustment 
    2627   USE sbcrnf         ! surface boundary condition: runoff variables 
    2728   ! 
    2829   USE in_out_manager ! I/O manager 
    2930   USE iom            ! IOM library 
     31   USE lbclnk         ! lateral boundary conditions 
    3032   USE lib_mpp        ! distribued memory computing 
    3133 
     
    3335   PRIVATE 
    3436 
    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) 
    4149   !! $Id$ 
    4250   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5159      !! 
    5260      !! ** Method  :   Read namelist namzdf, control logicals  
     61      !!                set horizontal shape and vertical profile of background mixing coef. 
    5362      !!---------------------------------------------------------------------- 
    5463      INTEGER ::   ioptio, ios   ! local integers 
     
    5867         &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc 
    5968         &             ln_zdfddm, rn_avts, rn_hsbfr,                 &     ! double diffusion 
    60          &             ln_zdftmx,                                    &     ! tidal mixing 
    61          &             ln_zdfqiao,                                   &     ! surface wave-induced mixing 
     69         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing 
     70         &             ln_zdfiwm,                                    &     ! internal  -      -      - 
    6271         &             ln_zdfexp, nn_zdfexp,                         &     ! time-stepping 
    6372         &             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  ==! 
    7176      REWIND( numnam_ref )              ! Namelist namzdf in reference namelist : Vertical mixing parameters 
    7277      READ  ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901) 
    73 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp ) 
    74  
     78901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp ) 
     79      ! 
    7580      REWIND( numnam_cfg )              ! Namelist namzdf in reference namelist : Vertical mixing parameters 
    7681      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 print 
     82902   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 
    8186         WRITE(numout,*) 
    82          WRITE(numout,*) 'zdf_phy_init : vertical physics' 
    83          WRITE(numout,*) '~~~~~~~~' 
     87         WRITE(numout,*) 'zdf_phy_init: vertical physics' 
     88         WRITE(numout,*) '~~~~~~~~~~~~' 
    8489         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters' 
    8590         WRITE(numout,*) '      vertical closure scheme' 
     
    98103         WRITE(numout,*) '         maximum avs for dd mixing                 rn_avts = ', rn_avts 
    99104         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 
    103110         WRITE(numout,*) '         number of sub-time step (ln_zdfexp=T)   nn_zdfexp = ', nn_zdfexp 
    104111         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 
    176120         avmb(:) = rn_avm0 
    177121         avtb(:) = rn_avt0                      
    178       ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
     122      ELSE                               ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
    179123         avmb(:) = rn_avm0 
    180124         avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s 
    181125         IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' ) 
    182126      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 
    200134           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.)) 
    201135           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1 
     
    203137      ENDIF 
    204138      ! 
    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 ? 
    224183      ! 
    225184   END SUBROUTINE zdf_phy_init 
    226185 
    227186 
    228    SUBROUTINE zdf_phy( kstp ) 
     187   SUBROUTINE zdf_phy( kt ) 
    229188      !!---------------------------------------------------------------------- 
    230189      !!                     ***  ROUTINE zdf_phy  *** 
     
    238197      !!                bottom stress.....                               <<<<====verifier ! 
    239198      !!---------------------------------------------------------------------- 
    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     
    275256                                                      ! write TKE or GLS information in the restart file 
    276       IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kstp, 'WRITE' ) 
    277       IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kstp, '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      !       
    279260   END SUBROUTINE zdf_phy 
    280261 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r7953 r7990  
    55   !!                 Richardson number dependent formulation 
    66   !!====================================================================== 
    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  
    1415   !!---------------------------------------------------------------------- 
    1516 
     
    2122   USE dom_oce        ! ocean space and time domain variables 
    2223   USE zdf_oce        ! ocean vertical physics 
     24   USE phycst         ! physical constants 
     25   USE sbc_oce,  ONLY :   taum 
     26   USE eosbn2 ,  ONLY :   neos 
     27   ! 
    2328   USE in_out_manager ! I/O manager 
    2429   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    2530   USE lib_mpp        ! MPP library 
    26    USE wrk_nemo       ! work arrays 
    2731   USE timing         ! Timing 
    2832   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2933 
    30    USE eosbn2, ONLY : neos 
    3134 
    3235   IMPLICIT NONE 
     
    3639   PUBLIC   zdf_ric_init    ! called by opa.F90 
    3740 
    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 
    5051 
    5152   !! * Substitutions 
     
    5758   !!---------------------------------------------------------------------- 
    5859CONTAINS 
    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_alloc 
    69  
    7060 
    7161   SUBROUTINE zdf_ric( kt ) 
     
    10595      !!      and bottom value already set in zdfphy.F90 
    10696      !! 
     97      !! ** Action  :   avm, avt  mixing coeff (inner domain values only) 
     98      !! 
    10799      !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 
    108100      !!              PFJ Lermusiaux 2001. 
    109101      !!---------------------------------------------------------------------- 
    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 
    129126         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         ! 
    184165      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') 
    229189      ! 
    230190   END SUBROUTINE zdf_ric 
     
    264224         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme' 
    265225         WRITE(numout,*) '~~~~~~~' 
    266          WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters' 
    267          WRITE(numout,*) '      maximum vertical viscosity     rn_avmri  = ', rn_avmri 
    268          WRITE(numout,*) '      coefficient                    rn_alp    = ', rn_alp 
    269          WRITE(numout,*) '      coefficient                    nn_ric    = ', nn_ric 
    270          WRITE(numout,*) '      Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc 
    271          WRITE(numout,*) '      minimum mixed layer depth      rn_mldmin = ', rn_mldmin 
    272          WRITE(numout,*) '      maximum mixed layer depth      rn_mldmax = ', rn_mldmax 
    273          WRITE(numout,*) '      Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix 
    274          WRITE(numout,*) '      Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix 
    275          WRITE(numout,*) '      Use the MLD parameterization   ln_mldw   = ', ln_mldw 
     226         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 
    276236      ENDIF 
    277237      ! 
    278       !                              ! allocate zdfric arrays 
    279       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-points 
    282          DO jj = 2, jpj              ! which accounts for coastal boundary conditions             
    283             DO ji = 2, jpi 
    284                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 DO 
    288          END DO 
    289       END DO 
    290       tmric(:,1,:) = 0._wp 
    291       ! 
    292238      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) 
    296241      END DO 
    297242      ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90

    r7967 r7990  
    1 MODULE zdfqiao 
     1MODULE zdfswm 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  zdfqiao  *** 
    4    !! Qiao module      : vertical mixing enhancement due to surface waves  
     3   !!                       ***  MODULE  zdfswm  *** 
     4   !! vertical physics :   surface wave-induced mixing  
    55   !!====================================================================== 
    66   !! 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 
    98   !!---------------------------------------------------------------------- 
    109 
    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 
    1822    
    1923   IMPLICIT NONE 
    2024   PRIVATE 
    2125 
    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 
    2328 
    24    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv 
    25  
    26    !! * Substitutions 
    27 #  include "vectopt_loop_substitute.h90" 
    2829   !!---------------------------------------------------------------------- 
    29    !! NEMO/OPA 4.0 , NEMO Consortium (2011)  
    30    !! $Id: $ 
     30   !! NEMO/OPA 4.0 , NEMO Consortium (2017)  
     31   !! $Id:$ 
    3132   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3233   !!---------------------------------------------------------------------- 
    33  
    3434CONTAINS 
    3535 
    36    SUBROUTINE zdf_qiao( kt ) 
     36   SUBROUTINE zdf_swm( kt ) 
    3737      !!--------------------------------------------------------------------- 
    38       !!                     ***  ROUTINE zdf_qiao *** 
     38      !!                     ***  ROUTINE zdf_swm *** 
    3939      !! 
    40       !! ** Purpose :Compute the Qiao term (qbv) to be added to 
     40      !! ** Purpose :Compute the swm term (qbv) to be added to 
    4141      !!             vertical viscosity and diffusivity coeffs.   
    4242      !! 
    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 
    4447      !!              
    45       !! ** action  :Compute the Qiao wave dependent term  
    46       !!             only if ln_wave=.true. 
     48      !! ** action  :    avt, avs, avm updated by the surface wave-induced mixing 
     49      !!                               (inner domain only) 
    4750      !!                
     51      !! reference : Qiao et al. GRL, 2004 
    4852      !!--------------------------------------------------------------------- 
    49       INTEGER, INTENT( in  ) ::  kt   ! ocean time step 
     53      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    5054      ! 
    51       INTEGER :: jj, ji, jk   ! dummy loop indices 
     55      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     56      REAL(wp)::   zcoef, zqb   ! local scalar 
    5257      !!--------------------------------------------------------------------- 
    5358      ! 
    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 
    7968            END DO 
    8069         END DO 
    8170      END DO 
    8271      ! 
    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      !!--------------------------------------------------------------------- 
    8583      ! 
    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') 
    8891      ! 
    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 
    10193 
    102       ! Enhance vertical mixing coeff.          
    103       !------------------------------- 
    104       ! 
    105 !!gm  with double diffusion activated, avs is not updated... 
    106 !!gm                    =====>>> BUG 
    107       DO jk = 1, jpkm1 
    108          DO jj = 1, jpj 
    109             DO ji = 1, jpi 
    110                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 DO 
    114          END DO 
    115       END DO 
    116       ! 
    117    END SUBROUTINE zdf_qiao 
    118  
    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_alloc 
    131        
    13294   !!====================================================================== 
    133 END MODULE zdfqiao 
     95END MODULE zdfswm 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r7953 r7990  
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    2828   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    29    !!            4.0  !  2017-04  (G. Madec)  Remove CPP keys 
     29   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only  
    3030   !!---------------------------------------------------------------------- 
    3131 
     
    159159      !! 
    160160      !! ** Action  :   compute en (now turbulent kinetic energy) 
    161       !!                update avt, avmu, avmv (before vertical eddy coef.) 
     161      !!                update avt, avm (before vertical eddy coef.) 
    162162      !! 
    163163      !! References : Gaspar et al., JGR, 1990, 
     
    175175#endif 
    176176      ! 
    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(:,:,:)  
    192186      ! 
    193187#if defined key_agrif 
     
    195189      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    196190#endif       
    197      !  
     191      !  
     192      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' ) 
     193      ! 
    198194  END SUBROUTINE zdf_tke 
    199195 
     
    207203      !! ** Method  : - TKE surface boundary condition 
    208204      !!              - 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] ) 
    210206      !!              - Now TKE : resolution of the TKE equation by inverting 
    211207      !!                a tridiagonal linear system by a "methode de chasse" 
     
    213209      !! 
    214210      !! ** Action  : - en : now turbulent kinetic energy) 
    215       !!              - avmu, avmv : production of TKE by shear at u and v-points 
    216       !!                (= Kz dz[Ub] * dz[Un] ) 
    217211      !! --------------------------------------------------------------------- 
    218212      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
     
    228222      REAL(wp) ::   zzd_up, zzd_lw                  !    -         - 
    229223!!bfr      REAL(wp) ::   zebot                           !    -         - 
    230       INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc 
     224      INTEGER , DIMENSION(jpi,jpj) ::   imlc 
    231225      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc 
    232226      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 
     
    236230      IF( nn_timing == 1 )  CALL timing_start('tke_tke') 
    237231      ! 
    238       CALL wrk_alloc( jpi,jpj,       imlc )    ! integer 
    239232      CALL wrk_alloc( jpi,jpj,       zhlc )  
    240233      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
     
    473466         END DO 
    474467      ENDIF 
     468!!gm not sure we need this lbc .... 
    475469      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    476470      ! 
    477       CALL wrk_dealloc( jpi,jpj,       imlc )    ! integer 
    478471      CALL wrk_dealloc( jpi,jpj,       zhlc )  
    479472      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
     
    516509      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 
    517510      !! 
    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) 
    520512      !!---------------------------------------------------------------------- 
    521513      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    522       REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars 
    523       REAL(wp) ::   zdku, zri, zsqen            !   -      - 
    524       REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      - 
     514      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
     515      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
     516      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
    525517      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld 
    526518      !!-------------------------------------------------------------------- 
     
    645637 
    646638      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    647       !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt) 
     639      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt) 
    648640      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    649641      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points 
     
    658650         END DO 
    659651      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      ! 
    671653      ! 
    672654      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
     
    677659# if defined key_c1d 
    678660                  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 Ri 
    682661# endif 
    683662              END DO 
     
    685664         END DO 
    686665      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) 
    688669 
    689670      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 ) 
    693673      ENDIF 
    694674      ! 
     
    759739      ENDIF 
    760740      ! 
    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) 
    765744         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) 
    767746         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity 
    768747         IF(lwp) WRITE(numout,*) '      minimum mixing length with your parameters rmxl_min = ', rmxl_min 
     
    798777         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    799778         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    800          avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
    801          avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    802779      END DO 
    803780      dissl(:,:,:) = 1.e-12_wp 
     
    821798     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    822799     ! 
    823      INTEGER ::   jit, jk   ! dummy loop indices 
    824      INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers 
     800     INTEGER ::   jit, jk              ! dummy loop indices 
     801     INTEGER ::   id1, id2, id3, id4   ! local integers 
    825802     !!---------------------------------------------------------------------- 
    826803     ! 
     
    829806        IF( ln_rstart ) THEN                   !* Read the restart file 
    830807           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. ) 
    836811           ! 
    837812           IF( id1 > 0 ) THEN                       ! 'en' exists 
    838813              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 ) 
    844817                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 
    845818              ELSE                                                 ! one at least array is missing 
    846                  CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation) 
     819                 CALL tke_avn                                          ! compute avt, avm, and dissl (approximation) 
    847820              ENDIF 
    848821           ELSE                                     ! No TKE array found: initialisation 
    849822              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(:,:,:) 
    857827              ! 
    858828              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
     829              avt_k(:,:,:) = avt(:,:,:) 
     830              avm_k(:,:,:) = avm(:,:,:) 
    859831           ENDIF 
    860832        ELSE                                   !* Start from rest 
    861833           en(:,:,:) = rn_emin * tmask(:,:,:) 
    862834           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) 
    867837           END DO 
    868838        ENDIF 
     
    871841        !                                   ! ------------------- 
    872842        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 ) 
    879847        ! 
    880848     ENDIF 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/module_example

    r4147 r7990  
    8686      !!                Give references if exist otherwise suppress these lines 
    8787      !!---------------------------------------------------------------------- 
    88       USE toto_module      ! description of the module 
    89       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    90       USE wrk_nemo, ONLY:   zztab => wrk_2d_5                     ! 2D workspace 
    91       USE wrk_nemo, ONLY:   zwx => wrk_3d_12 , zwy => wrk_3d_13   ! 3D workspace 
    92       !! 
    9388      INTEGER , INTENT(in   )                     ::   kt      ! short description  
    9489      INTEGER , INTENT(inout)                     ::   pvar1   !   -         - 
     
    10095      REAL(wp) ::   zmlmin, zbbrau   ! temporary scalars     (DOCTOR : start with z) 
    10196      REAL(wp) ::   zfact1, zfact2   ! do not use continuation lines in declaration 
     97      REAL(wp), DIMENSION(jpi,jpj) ::   zwrk_2d   ! 2D workspace 
    10298      !!-------------------------------------------------------------------- 
    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      ! 
    108100      IF( kt == nit000  )   CALL exa_mpl_init    ! Initialization (first time-step only) 
    109101 
     
    119111            DO jj = 2, jpjm1 
    120112               DO ji = fs_2, fs_jpim1   ! vector opt. 
    121                   avmv(ji,jj,jk) = .... 
     113                  avm(ji,jj,jk) = .... 
    122114               END DO 
    123115            END DO 
     
    128120            DO jj = 2, jpjm1 
    129121               DO ji = fs_2, fs_jpim1   ! vector opt. 
    130                   avmv(ji,jj,jk) = ... 
     122                  avm(ji,jj,jk) = ... 
    131123               END DO 
    132124            END DO 
     
    135127      END SELECT 
    136128      ! 
    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) 
    142130      ! 
    143131   END SUBROUTINE exa_mpl 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7953 r7990  
    6464 
    6565   USE zdfphy         ! vertical physics manager      (zdf_phy_init routine) 
    66 !!gm to be suppressed 
    67    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 end 
    7766 
    7867   USE step_diu        ! Time stepping for diurnal sst 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90

    r7681 r7990  
    119119      !                                          ! -------------- 
    120120      neln(:,:) = 1                                   ! euphotic layer level 
    121       DO jk = 1, jpk                                  ! (i.e. 1rst T-level strictly below EL bottom) 
     121      DO jk = 1, jpkm1                                ! (i.e. 1rst T-level strictly below EL bottom) 
    122122         DO jj = 1, jpj 
    123123           DO ji = 1, jpi 
     
    147147   END SUBROUTINE p2z_opt 
    148148 
     149 
    149150   SUBROUTINE p2z_opt_init 
    150151      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.