Changeset 8160


Ignore:
Timestamp:
2017-06-10T09:31:34+02:00 (3 years ago)
Author:
gm
Message:

#1883 (HPC-09) - step-8: remove split-explicit calculation of (tra/dyn)zdf + associated changes in namelist

Location:
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM
Files:
4 deleted
19 edited

Legend:

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

    r8143 r8160  
    340340   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    341341   ! 
    342    !                       ! time-stepping 
    343    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    344       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    345    ! 
    346342   !                       ! coefficients 
    347343   rn_avm0     =   0.1e-6     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg

    r8143 r8160  
    279279   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    280280   ! 
    281    !                       ! time-stepping 
    282    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    283       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    284    ! 
    285281   !                       ! coefficients 
    286282   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg

    r8143 r8160  
    268268   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    269269   ! 
    270    !                       ! time-stepping 
    271    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    272       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    273    ! 
    274270   !                       ! coefficients 
    275271   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg

    r8143 r8160  
    222222   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    223223   ! 
    224    !                       ! time-stepping 
    225    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    226       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    227    ! 
    228224   !                       !  Coefficients 
    229225   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg

    r8143 r8160  
    259259   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    260260   ! 
    261    !                       ! time-stepping 
    262    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    263       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    264    ! 
    265261   !                       !  Coefficients 
    266262   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8143 r8160  
    914914   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    915915   ! 
    916    !                       ! time-stepping 
    917    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    918       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    919    ! 
    920916   !                       ! coefficients 
    921917   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg

    r8143 r8160  
    339339   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    340340   ! 
    341    !                       ! time-stepping 
    342    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    343    ! 
    344341   !                       ! coefficients 
    345342   rn_avm0     =   1.e-3     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg

    r8143 r8160  
    201201   ln_zdfnpc   = .false.      !  Non-Penetrative Convective algorithm 
    202202   ! 
    203    !                       ! time-stepping 
    204    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    205    ! 
    206203   !                       ! coefficients 
    207204   rn_avm0     =   1.e-4      !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg

    r8143 r8160  
    198198   ln_zdfnpc   = .false.      !  Non-Penetrative Convective algorithm 
    199199   ! 
    200    !                       ! time-stepping 
    201    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    202    ! 
    203200   !                       ! coefficients 
    204201   rn_avm0     =   1.e-4      !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg

    r8143 r8160  
    359359   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T ) 
    360360   ! 
    361    !                       ! time-stepping 
    362    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme 
    363       nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T 
    364    ! 
    365361   !                       ! coefficients 
    366362   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r7990 r8160  
    66   !! History :  1.0  !  2005-11  (G. Madec)  Original code 
    77   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     8   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option 
    89   !!---------------------------------------------------------------------- 
    910 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   dyn_zdf       : Update the momentum trend with the vertical diffusion 
    12    !!   dyn_zdf_init  : initializations of the vertical diffusion scheme 
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and tracers variables 
     15   USE phycst         ! physical constants 
    1516   USE dom_oce        ! ocean space and time domain variables  
     17   USE sbc_oce        ! surface boundary condition: ocean 
    1618   USE zdf_oce        ! ocean vertical physics variables 
    17    USE dynzdf_exp     ! vertical diffusion: explicit (dyn_zdf_exp     routine) 
    18    USE dynzdf_imp     ! vertical diffusion: implicit (dyn_zdf_imp     routine) 
     19   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     20   USE dynadv   , ONLY: ln_dynadv_vec ! dynamics: advection form 
    1921   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    2022   USE trd_oce        ! trends: ocean variables 
     
    3032   PRIVATE 
    3133 
    32    PUBLIC   dyn_zdf       !  routine called by step.F90 
    33    PUBLIC   dyn_zdf_init  !  routine called by opa.F90 
    34  
    35    INTEGER  ::   nzdf = 0   ! type vertical diffusion algorithm used, defined from ln_zdf... namlist logicals 
     34   PUBLIC   dyn_zdf   !  routine called by step.F90 
     35 
     36   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise  
    3637 
    3738   !! * Substitutions 
    3839#  include "vectopt_loop_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     41   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    4142   !! $Id$ 
    4243   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     
    5758      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf') 
    5859      ! 
    59       !                                          ! set time step 
     60      !                             !* set time step 
    6061      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping) 
    6162      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog) 
    6263      ENDIF 
    6364 
    64       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     65      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    6566         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    6667         ztrdu(:,:,:) = ua(:,:,:) 
    6768         ztrdv(:,:,:) = va(:,:,:) 
    6869      ENDIF 
    69  
    70       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 
    71       ! 
    72       CASE ( 0 )   ;   CALL dyn_zdf_exp( kt, r2dt )      ! explicit scheme 
    73       CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme 
    74       ! 
    75       END SELECT 
    76  
     70      !                             !* compute lateral mixing trend and add it to the general trend 
     71      ! 
     72      CALL dyn_zdf_imp( kt, r2dt ) 
     73      ! 
    7774      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    7875         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
     
    9087 
    9188 
    92    SUBROUTINE dyn_zdf_init 
     89   SUBROUTINE dyn_zdf_imp( kt, p2dt ) 
    9390      !!---------------------------------------------------------------------- 
    94       !!                 ***  ROUTINE dyn_zdf_init  *** 
     91      !!                  ***  ROUTINE dyn_zdf_imp  *** 
     92      !!                    
     93      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion 
     94      !!              together with the Leap-Frog time stepping using an  
     95      !!              implicit scheme. 
    9596      !! 
    96       !! ** Purpose :   initialization of the vertical diffusion scheme 
     97      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing 
     98      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf. 
     99      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise 
     100      !!               - update the after velocity with the implicit vertical mixing. 
     101      !!      This requires to solver the following system:  
     102      !!         ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 
     103      !!      with the following surface/top/bottom boundary condition: 
     104      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
     105      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
    97106      !! 
    98       !! ** Method  :   implicit (euler backward) scheme (default) 
    99       !!                explicit (time-splitting) scheme if ln_zdfexp=T 
     107      !! ** Action :   (ua,va) after velocity  
     108      !!--------------------------------------------------------------------- 
     109      INTEGER , INTENT(in) ::  kt     ! ocean time-step index 
     110      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
     111      ! 
     112      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     113      INTEGER  ::   ikbu, ikbv    ! local integers 
     114      REAL(wp) ::   zzwi, ze3ua   ! local scalars 
     115      REAL(wp) ::   zzws, ze3va   !   -      - 
     116      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwd, zws 
    100117      !!---------------------------------------------------------------------- 
    101       USE zdftke 
    102       USE zdfgls 
    103       !!---------------------------------------------------------------------- 
    104       ! 
    105       ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module) 
    106       IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme 
    107       ELSE                   ;   nzdf = 1           ! use implicit scheme 
    108       ENDIF 
    109       ! 
    110       ! Force implicit schemes 
    111       IF( ln_zdftke .OR. ln_zdfgls   )   nzdf = 1   ! TKE or GLS physics 
    112       IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics 
    113       IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    114       ! 
    115       IF(lwp) THEN                                  ! Print the choice 
    116          WRITE(numout,*) 
    117          WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 
    118          WRITE(numout,*) '~~~~~~~~~~~' 
    119          IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme' 
    120          IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme' 
    121       ENDIF 
    122       ! 
    123    END SUBROUTINE dyn_zdf_init 
     118      ! 
     119      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp') 
     120      ! 
     121      IF( kt == nit000 ) THEN 
     122         IF(lwp) WRITE(numout,*) 
     123         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
     124         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     125         ! 
     126         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator 
     127         ELSE                   ;    r_vvl = 1._wp 
     128         ENDIF 
     129      ENDIF 
     130      ! 
     131      !              !==  Time step dynamics  ==! 
     132      ! 
     133      IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity 
     134         DO jk = 1, jpkm1 
     135            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     136            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     137         END DO 
     138      ELSE                                         ! applied on thickness weighted velocity 
     139         DO jk = 1, jpkm1 
     140            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
     141               &          + p2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
     142            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
     143               &          + p2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
     144         END DO 
     145      ENDIF 
     146      ! 
     147      !              !==  Apply semi-implicit bottom friction  ==! 
     148      ! 
     149      ! Only needed for semi-implicit bottom friction setup. The explicit 
     150      ! bottom friction has been included in "u(v)a" which act as the R.H.S 
     151      ! column vector of the tri-diagonal matrix equation 
     152      ! 
     153      IF( ln_drgimp ) THEN 
     154         DO jj = 2, jpjm1 
     155            DO ji = 2, jpim1 
     156               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
     157               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     158               avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 
     159               avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 
     160            END DO 
     161         END DO 
     162         IF ( ln_isfcav ) THEN 
     163            DO jj = 2, jpjm1 
     164               DO ji = 2, jpim1 
     165                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
     166                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     167                  ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 
     168                  avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 
     169                  avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 
     170               END DO 
     171            END DO 
     172         END IF 
     173      ENDIF 
     174      ! 
     175      ! With split-explicit free surface, barotropic stress is treated explicitly 
     176      ! Update velocities at the bottom. 
     177      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
     178      !            not lead to the effective stress seen over the whole barotropic loop.  
     179      ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
     180      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
     181         DO jk = 1, jpkm1        ! remove barotropic velocities 
     182            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 
     183            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
     184         END DO 
     185         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
     186            DO ji = fs_2, fs_jpim1   ! vector opt. 
     187               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     188               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     189               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
     190               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
     191               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     192               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     193            END DO 
     194         END DO 
     195         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
     196            DO jj = 2, jpjm1         
     197               DO ji = fs_2, fs_jpim1   ! vector opt. 
     198                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
     199                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     200                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
     201                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
     202                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     203                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     204               END DO 
     205            END DO 
     206         END IF 
     207      ENDIF 
     208      ! 
     209      !              !==  Vertical diffusion on u  ==! 
     210      ! 
     211      ! Matrix and second member construction 
     212      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     213      ! non zero value at the ocean bottom depending on the bottom friction used. 
     214      ! 
     215      DO jk = 1, jpkm1        ! Matrix 
     216         DO jj = 2, jpjm1  
     217            DO ji = fs_2, fs_jpim1   ! vector opt. 
     218               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
     219               zzwi = - p2dt * avmu(ji,jj,jk  ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) 
     220               zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 
     221               zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk  ) 
     222               zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 
     223               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     224            END DO 
     225         END DO 
     226      END DO 
     227      DO jj = 2, jpjm1        ! Surface boundary conditions 
     228         DO ji = fs_2, fs_jpim1   ! vector opt. 
     229            zwi(ji,jj,1) = 0._wp 
     230            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     231         END DO 
     232      END DO 
     233 
     234      ! Matrix inversion starting from the first level 
     235      !----------------------------------------------------------------------- 
     236      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     237      ! 
     238      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     239      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     240      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     241      !        (        ...               )( ...  ) ( ...  ) 
     242      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     243      ! 
     244      !   m is decomposed in the product of an upper and a lower triangular matrix 
     245      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
     246      !   The solution (the after velocity) is in ua 
     247      !----------------------------------------------------------------------- 
     248      ! 
     249      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     250         DO jj = 2, jpjm1    
     251            DO ji = fs_2, fs_jpim1   ! vector opt. 
     252               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     253            END DO 
     254         END DO 
     255      END DO 
     256      ! 
     257      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
     258         DO ji = fs_2, fs_jpim1   ! vector opt. 
     259            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
     260            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     261               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
     262         END DO 
     263      END DO 
     264      DO jk = 2, jpkm1 
     265         DO jj = 2, jpjm1 
     266            DO ji = fs_2, fs_jpim1 
     267               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
     268            END DO 
     269         END DO 
     270      END DO 
     271      ! 
     272      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
     273         DO ji = fs_2, fs_jpim1   ! vector opt. 
     274            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     275         END DO 
     276      END DO 
     277      DO jk = jpk-2, 1, -1 
     278         DO jj = 2, jpjm1 
     279            DO ji = fs_2, fs_jpim1 
     280               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     281            END DO 
     282         END DO 
     283      END DO 
     284      ! 
     285      !              !==  Vertical diffusion on v  ==! 
     286      ! 
     287      ! Matrix and second member construction 
     288      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
     289      ! non zero value at the ocean bottom depending on the bottom friction used 
     290      ! 
     291      DO jk = 1, jpkm1        ! Matrix 
     292         DO jj = 2, jpjm1    
     293            DO ji = fs_2, fs_jpim1   ! vector opt. 
     294               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     295               zzwi = - p2dt * avmv (ji,jj,jk  ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) 
     296               zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 
     297               zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     298               zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     299               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     300            END DO 
     301         END DO 
     302      END DO 
     303      DO jj = 2, jpjm1        ! Surface boundary conditions 
     304         DO ji = fs_2, fs_jpim1   ! vector opt. 
     305            zwi(ji,jj,1) = 0._wp 
     306            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     307         END DO 
     308      END DO 
     309 
     310      ! Matrix inversion 
     311      !----------------------------------------------------------------------- 
     312      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     313      ! 
     314      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     315      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     316      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     317      !        (        ...               )( ...  ) ( ...  ) 
     318      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     319      ! 
     320      !   m is decomposed in the product of an upper and lower triangular matrix 
     321      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
     322      !   The solution (after velocity) is in 2d array va 
     323      !----------------------------------------------------------------------- 
     324      ! 
     325      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     326         DO jj = 2, jpjm1    
     327            DO ji = fs_2, fs_jpim1   ! vector opt. 
     328               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     329            END DO 
     330         END DO 
     331      END DO 
     332      ! 
     333      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
     334         DO ji = fs_2, fs_jpim1   ! vector opt.           
     335            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
     336            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     337               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
     338         END DO 
     339      END DO 
     340      DO jk = 2, jpkm1 
     341         DO jj = 2, jpjm1 
     342            DO ji = fs_2, fs_jpim1   ! vector opt. 
     343               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
     344            END DO 
     345         END DO 
     346      END DO 
     347      ! 
     348      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
     349         DO ji = fs_2, fs_jpim1   ! vector opt. 
     350            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     351         END DO 
     352      END DO 
     353      DO jk = jpk-2, 1, -1 
     354         DO jj = 2, jpjm1 
     355            DO ji = fs_2, fs_jpim1 
     356               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     357            END DO 
     358         END DO 
     359      END DO 
     360      ! 
     361      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp') 
     362      ! 
     363   END SUBROUTINE dyn_zdf_imp 
    124364 
    125365   !!============================================================================== 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r7954 r8160  
    44   !! Ocean active tracers:  vertical component of the tracer mixing trend 
    55   !!============================================================================== 
    6    !! History :  1.0  ! 2005-11  (G. Madec)  Original code 
    7    !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
     6   !! History :  1.0  !  2005-11  (G. Madec)  Original code 
     7   !!            3.0  !  2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
     8   !!            4.0  !  2017-06  (G. Madec)  remove explict time-stepping option 
    89   !!---------------------------------------------------------------------- 
    910 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   tra_zdf       : Update the tracer trend with the vertical diffusion 
    12    !!   tra_zdf_init  : initialisation of the computation 
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and tracers variables 
     
    2020   USE ldftra         ! lateral diffusion: eddy diffusivity 
    2121   USE ldfslp         ! lateral diffusion: iso-neutral slope  
    22    USE trazdf_exp     ! vertical diffusion: explicit (tra_zdf_exp routine) 
    23    USE trazdf_imp     ! vertical diffusion: implicit (tra_zdf_imp routine) 
    2422   USE trd_oce        ! trends: ocean variables 
    2523   USE trdtra         ! trends: tracer trend manager 
     
    2927   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3028   USE lib_mpp        ! MPP library 
    31    USE wrk_nemo       ! Memory allocation 
    3229   USE timing         ! Timing 
    3330 
     
    3532   PRIVATE 
    3633 
    37    PUBLIC   tra_zdf        ! routine called by step.F90 
    38    PUBLIC   tra_zdf_init   ! routine called by nemogcm.F90 
    39  
    40    INTEGER ::   nzdf = 0   ! type vertical diffusion algorithm used (defined from ln_zdf...  namlist logicals) 
     34   PUBLIC   tra_zdf       ! called by step.F90 
     35   PUBLIC   tra_zdf_imp   ! called by trczdf.F90 
    4136 
    4237   !! * Substitutions 
     
    5853      ! 
    5954      INTEGER  ::   jk                   ! Dummy loop indices 
    60       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace 
     55      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace 
    6156      !!--------------------------------------------------------------------- 
    6257      ! 
     
    6964      ENDIF 
    7065      ! 
    71       IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    72          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     66      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
     67         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    7368         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    7469         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    7570      ENDIF 
    7671      ! 
    77       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 
    78       CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme  
    79       CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme  
    80       END SELECT 
     72      !                                      !* compute lateral mixing trend and add it to the general trend 
     73      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )  
     74 
    8175!!gm WHY here !   and I don't like that ! 
    8276      ! DRAKKAR SSS control { 
     
    9791         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    9892         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    99          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     93         DEALLOCATE( ztrdt , ztrds ) 
    10094      ENDIF 
    10195      !                                          ! print mean trends (used for debugging) 
     
    107101   END SUBROUTINE tra_zdf 
    108102 
    109  
    110    SUBROUTINE tra_zdf_init 
     103  
     104   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )  
    111105      !!---------------------------------------------------------------------- 
    112       !!                 ***  ROUTINE tra_zdf_init  *** 
    113       !! 
    114       !! ** Purpose :   Choose the vertical mixing scheme 
    115       !! 
    116       !! ** Method  :   Set nzdf from ln_zdfexp 
    117       !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T) 
    118       !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F) 
    119       !!      NB: rotation of lateral mixing operator or TKE & GLS schemes, 
    120       !!          an implicit scheme is required. 
    121       !!---------------------------------------------------------------------- 
    122       ! 
    123       ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module) 
    124       IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme 
    125       ELSE                   ;   nzdf = 1           ! use implicit scheme 
    126       ENDIF 
    127       ! 
    128       ! Force implicit schemes when absolutely needed 
    129       IF( .NOT.ln_zdfcst             )   nzdf = 1   ! turbulent closure schemes 
    130       IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics 
    131       IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    132       ! 
    133       IF(lwp) THEN 
    134          WRITE(numout,*) 
    135          WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 
    136          WRITE(numout,*) '~~~~~~~~~~~' 
    137          IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme' 
    138          IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme' 
    139       ENDIF 
    140       ! 
    141    END SUBROUTINE tra_zdf_init 
     106      !!                  ***  ROUTINE tra_zdf_imp  *** 
     107      !! 
     108      !! ** Purpose :   Compute the after tracer through a implicit computation 
     109      !!     of the vertical tracer diffusion (including the vertical component  
     110      !!     of lateral mixing (only for 2nd order operator, for fourth order  
     111      !!     it is already computed and add to the general trend in traldf)  
     112      !! 
     113      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by: 
     114      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 
     115      !!      It is computed using a backward time scheme (t=after field) 
     116      !!      which provide directly the after tracer field. 
     117      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers 
     118      !!      Surface and bottom boundary conditions: no diffusive flux on 
     119      !!      both tracers (bottom, applied through the masked field avt). 
     120      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
     121      !! 
     122      !! ** Action  : - pta  becomes the after tracer 
     123      !!--------------------------------------------------------------------- 
     124      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
     125      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index 
     126      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     127      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
     128      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields 
     130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field 
     131      ! 
     132      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     133      REAL(wp) ::  zrhs             ! local scalars 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws 
     135      !!--------------------------------------------------------------------- 
     136      ! 
     137      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp') 
     138      ! 
     139      IF( kt == kit000 )  THEN 
     140         IF(lwp)WRITE(numout,*) 
     141         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 
     142         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 
     143      ENDIF 
     144      !                                               ! ============= ! 
     145      DO jn = 1, kjpt                                 !  tracer loop  ! 
     146         !                                            ! ============= ! 
     147         !  Matrix construction 
     148         ! -------------------- 
     149         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 
     150         ! 
     151         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   & 
     152            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN 
     153            ! 
     154            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 
     155            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk) 
     156            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk) 
     157            ENDIF 
     158            zwt(:,:,1) = 0._wp 
     159            ! 
     160            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
     161               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
     162                  DO jk = 2, jpkm1 
     163                     DO jj = 2, jpjm1 
     164                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     165                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     166                        END DO 
     167                     END DO 
     168                  END DO 
     169               ELSE                          ! standard or triad iso-neutral operator 
     170                  DO jk = 2, jpkm1 
     171                     DO jj = 2, jpjm1 
     172                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     173                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     174                        END DO 
     175                     END DO 
     176                  END DO 
     177               ENDIF 
     178            ENDIF 
     179            ! 
     180            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
     181            DO jk = 1, jpkm1 
     182               DO jj = 2, jpjm1 
     183                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     184!!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that 
     185                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
     186                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
     187                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     188                 END DO 
     189               END DO 
     190            END DO 
     191            ! 
     192            !! Matrix inversion from the first level 
     193            !!---------------------------------------------------------------------- 
     194            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     195            ! 
     196            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     197            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     198            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     199            !        (        ...               )( ...  ) ( ...  ) 
     200            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     201            ! 
     202            !   m is decomposed in the product of an upper and lower triangular matrix. 
     203            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 
     204            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal 
     205            !   and "superior" (above diagonal) components of the tridiagonal system. 
     206            !   The solution will be in the 4d array pta. 
     207            !   The 3d array zwt is used as a work space array. 
     208            !   En route to the solution pta is used a to evaluate the rhs and then  
     209            !   used as a work space array: its value is modified. 
     210            ! 
     211            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
     212               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction) 
     213                  zwt(ji,jj,1) = zwd(ji,jj,1) 
     214               END DO 
     215            END DO 
     216            DO jk = 2, jpkm1 
     217               DO jj = 2, jpjm1 
     218                  DO ji = fs_2, fs_jpim1 
     219                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     220                  END DO 
     221               END DO 
     222            END DO 
     223            ! 
     224         ENDIF  
     225         !          
     226         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     227            DO ji = fs_2, fs_jpim1 
     228               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 
     229            END DO 
     230         END DO 
     231         DO jk = 2, jpkm1 
     232            DO jj = 2, jpjm1 
     233               DO ji = fs_2, fs_jpim1 
     234                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side 
     235                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
     236               END DO 
     237            END DO 
     238         END DO 
     239         ! 
     240         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
     241            DO ji = fs_2, fs_jpim1 
     242               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     243            END DO 
     244         END DO 
     245         DO jk = jpk-2, 1, -1 
     246            DO jj = 2, jpjm1 
     247               DO ji = fs_2, fs_jpim1 
     248                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
     249                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     250               END DO 
     251            END DO 
     252         END DO 
     253         !                                            ! ================= ! 
     254      END DO                                          !  end tracer loop  ! 
     255      !                                               ! ================= ! 
     256      ! 
     257      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp') 
     258      ! 
     259   END SUBROUTINE tra_zdf_imp 
    142260 
    143261   !!============================================================================== 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r8143 r8160  
    3737   LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag 
    3838   LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag 
    39    !                             ! time-stepping 
    40    LOGICAL , PUBLIC ::   ln_zdfexp   !: explicit vertical diffusion scheme flag 
    41    INTEGER , PUBLIC ::      nn_zdfexp   !: number of sub-time step (explicit time stepping) 
    4239   !                             ! coefficients  
    4340   REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90

    r8143 r8160  
    190190      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' ) 
    191191      ! 
    192       IF( ln_drgimp .AND. ln_zdfexp )   CALL ctl_stop( 'zdf_drg_init: ln_drgimp=T requires ln_zdfexp' ) 
    193  
    194192      ! 
    195193      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r8143 r8160  
    327327      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag 
    328328      zd_lw(:,:,2) = 0._wp 
    329       zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
     329      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
    330330      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    331331          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
     
    556556         ! 
    557557         ! One level below 
    558          zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
     558         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
    559559         zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 
    560560         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     
    577577         ! 
    578578         ! Set psi vertical flux at the surface: 
    579          zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     579         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    580580         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    581581         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     
    841841      !! 
    842842      !!---------------------------------------------------------------------- 
    843       USE dynzdf_exp 
    844       USE trazdf_exp 
    845       ! 
    846843      INTEGER ::   jk    ! dummy loop indices 
    847844      INTEGER ::   ios   ! Local integer output status for namelist read 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

    r8143 r8160  
    7373         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing 
    7474         &             ln_zdfiwm,                                    &     ! internal  -      -      - 
    75          &             ln_zdfexp, nn_zdfexp,                         &     ! time-stepping 
    7675         &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients 
    7776      !!---------------------------------------------------------------------- 
     
    110109         WRITE(numout,*) '         surface  wave (Qiao et al 2010)         ln_zdfswm = ', ln_zdfswm                                          ! surface wave induced mixing 
    111110         WRITE(numout,*) '         internal wave (de Lavergne et al 2017)  ln_zdfiwm = ', ln_zdfiwm 
    112          WRITE(numout,*) '      time-steping scheme' 
    113          WRITE(numout,*) '         time splitting (T) / implicit (F)       ln_zdfexp = ', ln_zdfexp 
    114          WRITE(numout,*) '         number of sub-time step (ln_zdfexp=T)   nn_zdfexp = ', nn_zdfexp 
    115111         WRITE(numout,*) '      coefficients : ' 
    116112         WRITE(numout,*) '         vertical eddy viscosity                 rn_avm0   = ', rn_avm0 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r7954 r8160  
    471471                            CALL tra_adv_init      ! horizontal & vertical advection 
    472472                            CALL tra_ldf_init      ! lateral mixing 
    473                             CALL tra_zdf_init      ! vertical mixing and after tracer fields 
    474473 
    475474      !                                      ! Dynamics 
     
    479478                            CALL dyn_ldf_init      ! lateral mixing 
    480479                            CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure 
    481                             CALL dyn_zdf_init      ! vertical diffusion 
    482480                            CALL dyn_spg_init      ! surface pressure gradient 
    483481 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8143 r8160  
    344344   END SUBROUTINE stp 
    345345    
     346   !!====================================================================== 
    346347END MODULE step 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90

    r7953 r8160  
    1717   USE oce_trc       ! ocean dynamics and active tracers 
    1818   USE trd_oce       ! trends: ocean variables 
    19    USE trazdf_exp    ! vertical diffusion: explicit (tra_zdf_exp     routine) 
    20    USE trazdf_imp    ! vertical diffusion: implicit (tra_zdf_imp     routine) 
     19   USE trazdf        ! tracer: vertical diffusion 
     20!!gm do we really need this ? 
    2121   USE trcldf        ! passive tracers: lateral diffusion 
     22!!gm 
    2223   USE trdtra        ! trends manager: tracers  
    2324   USE prtctl_trc    ! Print control 
Note: See TracChangeset for help on using the changeset viewer.