Ignore:
Timestamp:
2014-12-15T17:42:49+01:00 (6 years ago)
Author:
timgraham
Message:

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r4147 r4990  
    4040   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile 
    4141   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: Bottom friction coefficients set in zdfbfr 
     42   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top friction coefficients set in zdfbfr 
    4243   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s] 
    4344   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm  , avt     !: vertical viscosity & diffusivity coef at w-pt [m2/s] 
     
    5758      ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) ,                         & 
    5859         &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      & 
     60         &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      & 
    5961         &     avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk)           ,      & 
    6062         &     avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk)           , STAT = zdf_oce_alloc ) 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4624 r4990  
    4141   REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM) 
    4242   LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM) 
     43   REAL(wp), PUBLIC ::   rn_tfri1     ! top drag coefficient (linear case)  (PUBLIC for TAM) 
     44   REAL(wp), PUBLIC ::   rn_tfri2     ! top drag coefficient (non linear case) (PUBLIC for TAM) 
     45   REAL(wp), PUBLIC ::   rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM) 
     46   REAL(wp), PUBLIC ::   rn_tfeb2     ! background top turbulent kinetic energy  [m2/s2] (PUBLIC for TAM) 
     47   REAL(wp), PUBLIC ::   rn_tfrien    ! local factor to enhance coefficient tfri (PUBLIC for TAM) 
     48   LOGICAL , PUBLIC ::   ln_tfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM) 
     49 
    4350   LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM) 
    4451   REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
     52   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
    4553   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction 
    46    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM) 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
    4755 
    4856   !! * Substitutions 
     
    6068      !!                ***  FUNCTION zdf_bfr_alloc  *** 
    6169      !!---------------------------------------------------------------------- 
    62       ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
     70      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
    6371      ! 
    6472      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc ) 
     
    8896      INTEGER  ::   ikbt, ikbu, ikbv             ! local integers 
    8997      REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars 
    90       REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt 
     98      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt, ztfrt 
    9199      !!---------------------------------------------------------------------- 
    92100      ! 
     
    101109      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only 
    102110         ! 
    103          CALL wrk_alloc( jpi, jpj, zbfrt ) 
     111         CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt ) 
    104112 
    105113         IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient 
    106114 
    107 #  if defined key_vectopt_loop 
    108             DO jj = 1, 1 
    109 !CDIR NOVERRCHK 
    110                DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    111 #  else 
    112 !CDIR NOVERRCHK 
    113115            DO jj = 1, jpj 
    114 !CDIR NOVERRCHK 
    115116               DO ji = 1, jpi 
    116 #  endif 
    117117                  ikbt = mbkt(ji,jj) 
    118 ! JC: possible WAD implementation should modify line below if layers vanish 
     118!! JC: possible WAD implementation should modify line below if layers vanish 
    119119                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
    120120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    121121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
     122! (ISF) 
     123                  ikbt = mikt(ji,jj) 
     124! JC: possible WAD implementation should modify line below if layers vanish 
     125                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     126                  ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     127                  ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
     128 
    122129               END DO 
    123130            END DO 
     
    125132         ELSE 
    126133            zbfrt(:,:) = bfrcoef2d(:,:) 
    127          ENDIF 
    128  
    129 # if defined key_vectopt_loop 
    130          DO jj = 1, 1 
    131 !CDIR NOVERRCHK 
    132             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    133 # else 
    134 !CDIR NOVERRCHK 
     134            ztfrt(:,:) = tfrcoef2d(:,:) 
     135         ENDIF 
     136 
    135137         DO jj = 2, jpjm1 
    136 !CDIR NOVERRCHK 
    137138            DO ji = 2, jpim1 
    138 # endif 
    139139               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
    140140               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     
    150150               bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu 
    151151               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv 
     152               ! 
     153               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     154               IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
     155                  bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
     156                               &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
     157                               &          * zecu * (1._wp - umask(ji,jj,1)) 
     158               END IF 
     159               IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
     160                  bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
     161                               &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
     162                               &          * zecv * (1._wp - vmask(ji,jj,1)) 
     163               END IF 
     164               ! (ISF) ======================================================================== 
     165               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
     166               ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
     167               ! 
     168               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     169                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     170               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     171                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
     172               ! 
     173               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
     174               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
     175               ! 
     176               tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
     177               tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
     178               ! (ISF) END ==================================================================== 
     179               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     180               IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
     181                  tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
     182                               &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
     183                               &          * zecu * (1._wp - umask(ji,jj,1)) 
     184               END IF 
     185               IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
     186                  tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
     187                               &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
     188                               &          * zecv * (1._wp - vmask(ji,jj,1)) 
     189               END IF 
    152190            END DO 
    153191         END DO 
    154  
    155192         ! 
    156193         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
     
    158195         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    159196            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    160          CALL wrk_dealloc( jpi,jpj,zbfrt ) 
     197         CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt ) 
    161198      ENDIF 
    162199      ! 
     
    183220      INTEGER   ::   ios 
    184221      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars 
     222      REAL(wp)  ::   zmintfr, zmaxtfr   ! temporary scalars 
    185223      REAL(wp)  ::   ztmp, zfru, zfrv   !    -         - 
    186224      !! 
    187225      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 
    188                     &  rn_bfrien, ln_bfrimp, ln_loglayer 
     226                    &          rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, & 
     227                    &  rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer 
    189228      !!---------------------------------------------------------------------- 
    190229      ! 
     
    215254         bfrua(:,:) = 0.e0 
    216255         bfrva(:,:) = 0.e0 
     256         tfrua(:,:) = 0.e0 
     257         tfrva(:,:) = 0.e0 
    217258         ! 
    218259      CASE( 1 ) 
    219260         IF(lwp) WRITE(numout,*) '      linear botton friction' 
    220          IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1 
     261         IF(lwp) WRITE(numout,*) '      bottom friction coef.   rn_bfri1  = ', rn_bfri1 
    221262         IF( ln_bfr2d ) THEN 
    222263            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    223264            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
     265         ENDIF 
     266         IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_bfri1 
     267         IF( ln_tfr2d ) THEN 
     268            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     269            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
    224270         ENDIF 
    225271         ! 
     
    252298            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    253299         ENDIF 
     300         IF(lwp) WRITE(numout,*) '      quadratic top    friction' 
     301         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_tfri2 
     302         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max 
     303         IF(lwp) WRITE(numout,*) '      background tke   rn_tfeb2  = ', rn_tfeb2 
     304         IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d = ', ln_loglayer 
     305         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_tfrz0 [m] = ', rn_tfrz0 
     306         IF( rn_tfrz0<=0.e0 ) THEN 
     307            WRITE(ctmp1,*) '      bottom roughness must be strictly positive' 
     308            CALL ctl_stop( ctmp1 ) 
     309         ENDIF 
     310         IF( ln_tfr2d ) THEN 
     311            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     312            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
     313         ENDIF 
    254314         ! 
    255315         IF(ln_bfr2d) THEN 
     
    265325         ! 
    266326         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 
    267 #  if defined key_vectopt_loop 
    268             DO jj = 1, 1 
    269 !CDIR NOVERRCHK 
    270                DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    271 #  else 
    272 !CDIR NOVERRCHK 
    273327            DO jj = 1, jpj 
    274 !CDIR NOVERRCHK 
    275328               DO ji = 1, jpi 
    276 #  endif 
    277329                  ikbt = mbkt(ji,jj) 
    278330                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     
    308360      zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
    309361      zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
    310       ! 
    311 #  if defined key_vectopt_loop 
    312       DO jj = 1, 1 
    313 !CDIR NOVERRCHK 
    314          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    315 #  else 
    316 !CDIR NOVERRCHK 
     362      zmintfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
     363      zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
     364      ! 
    317365      DO jj = 2, jpjm1 
    318 !CDIR NOVERRCHK 
    319366         DO ji = 2, jpim1 
    320 #  endif 
    321367             ikbu = mbku(ji,jj)       ! deepest ocean level at u- and v-points 
    322368             ikbv = mbkv(ji,jj) 
     
    352398         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
    353399         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 
     400         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
    354401         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
    355402      ENDIF 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r4624 r4990  
    66   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing 
    77   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    8    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_zdfddm   ||   defined key_esopa 
     
    1819   USE dom_oce         ! ocean space and time domain variables  
    1920   USE zdf_oce         ! ocean vertical physics variables 
     21   USE eosbn2         ! equation of state 
     22   ! 
    2023   USE in_out_manager  ! I/O manager 
    2124   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    3437   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag 
    3538 
    36    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point 
    37    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio 
    38  
    39    !                      !!* Namelist namzdf_ddm : double diffusive mixing * 
    40    REAL(wp) ::   rn_avts   ! maximum value of avs for salt fingering 
    41    REAL(wp) ::   rn_hsbfr  ! heat/salt buoyancy flux ratio 
     39   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point 
     40 
     41   !                       !!* Namelist namzdf_ddm : double diffusive mixing * 
     42   REAL(wp) ::   rn_avts    ! maximum value of avs for salt fingering 
     43   REAL(wp) ::   rn_hsbfr   ! heat/salt buoyancy flux ratio 
    4244 
    4345   !! * Substitutions 
     46#  include "domzgr_substitute.h90" 
    4447#  include "vectopt_loop_substitute.h90" 
    4548   !!---------------------------------------------------------------------- 
    46    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     49   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    4750   !! $Id$ 
    4851   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5457      !!                ***  ROUTINE zdf_ddm_alloc  *** 
    5558      !!---------------------------------------------------------------------- 
    56       ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc ) 
    57       ! 
     59      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc ) 
    5860      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc ) 
    5961      IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays') 
     
    7173      !!      diffusive mixing (i.e. salt fingering and diffusive layering) 
    7274      !!      following Merryfield et al. (1999). The rate of double diffusive  
    73       !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S] 
    74       !!      which is computed in rn2.F 
     75      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]): 
    7576      !!         * salt fingering (Schmitt 1981): 
    76       !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 ) 
    77       !!      for Rrau > 1 and rn2 > 0 : zavfs = O 
    78       !!      otherwise                : zavft = 0.7 zavs / Rrau 
     77      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 ) 
     78      !!      for R > 1 and rn2 > 0 : zavfs = O 
     79      !!      otherwise                : zavft = 0.7 zavs / R 
    7980      !!         * diffusive layering (Federov 1988): 
    80       !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6   
    81       !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) ) 
     81      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) ) 
    8282      !!      otherwise                   : zavdt = 0  
    83       !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85) 
    84       !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau       
     83      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85) 
     84      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R       
    8585      !!      otherwise                     : zavds = 0  
    8686      !!         * update the eddy diffusivity: 
     
    9696      ! 
    9797      INTEGER  ::   ji, jj , jk     ! dummy loop indices 
    98       REAL(wp) ::   zinr, zrr       ! temporary scalars 
    99       REAL(wp) ::   zavft, zavfs    !    -         - 
    100       REAL(wp) ::   zavdt, zavds    !    -         - 
    101       REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3 
     98      REAL(wp) ::   zaw, zbw, zrw   ! local scalars 
     99      REAL(wp) ::   zdt, zds 
     100      REAL(wp) ::   zinr, zrr       !   -      - 
     101      REAL(wp) ::   zavft, zavfs    !   -      - 
     102      REAL(wp) ::   zavdt, zavds    !   -      - 
     103      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
    102104      !!---------------------------------------------------------------------- 
    103105      ! 
    104106      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm') 
    105107      ! 
    106       CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
    107  
     108      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
     109      ! 
    108110      !                                                ! =============== 
    109111      DO jk = 2, jpkm1                                 ! Horizontal slab 
     
    111113         ! Define the mask  
    112114         ! --------------- 
    113          rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau 
     115         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
     116            DO ji = 1, jpi 
     117               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     118                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     119               ! 
     120               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  & 
     121                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     122               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  & 
     123                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     124               ! 
     125               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 
     126               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) )  
     127               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp 
     128               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau 
     129            END DO 
     130         END DO 
    114131 
    115132         DO jj = 1, jpj                                     ! indicators: 
     
    119136               ELSE                                       ;   zmsks(ji,jj) = 1._wp 
    120137               ENDIF 
    121                ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere             
    122                IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp 
     138               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere             
     139               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp 
    123140               ELSE                                       ;   zmskf(ji,jj) = 1._wp 
    124141               ENDIF 
    125142               ! diffusive layering indicators:  
    126                !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere 
    127                IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp 
     143               !     ! mskdl1=1 if 0< R <1; 0 elsewhere 
     144               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp 
    128145               ELSE                                       ;   zmskd1(ji,jj) = 1._wp 
    129146               ENDIF 
    130                !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere 
    131                IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp 
     147               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere 
     148               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp 
    132149               ELSE                                       ;   zmskd2(ji,jj) = 1._wp 
    133150               ENDIF 
    134                !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere 
    135                IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
    136                ELSE                                                         ;   zmskd3(ji,jj) = 1._wp 
     151               !   mskdl3=1 if 0.5< R <1; 0 elsewhere 
     152               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
     153               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp 
    137154               ENDIF 
    138155            END DO 
     
    149166!CDIR NOVERRCHK 
    150167            DO ji = 1, jpi 
    151                zinr = 1./rrau(ji,jj,jk) 
     168               zinr = 1._wp / zrau(ji,jj) 
    152169               ! salt fingering 
    153                zrr = rrau(ji,jj,jk)/rn_hsbfr 
     170               zrr = zrau(ji,jj) / rn_hsbfr 
    154171               zrr = zrr * zrr 
    155172               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj) 
     
    157174               ! diffusive layering 
    158175               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj) 
    159                zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   & 
    160                   &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  ) 
     176               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   & 
     177                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    161178               ! add to the eddy viscosity coef. previously computed 
    162179               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
     
    196213      ENDIF 
    197214      ! 
    198       CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
     215      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
    199216      ! 
    200217      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm') 
     
    212229      !!              called by zdf_ddm at the first timestep (nit000) 
    213230      !!---------------------------------------------------------------------- 
     231      INTEGER ::   ios   ! local integer 
     232      !! 
    214233      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr 
    215       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    216234      !!---------------------------------------------------------------------- 
    217235      ! 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r3294 r4990  
    7878         ! 
    7979         DO jk = 1, jpkm1  
    80 #if defined key_vectopt_loop 
    81             DO jj = 1, 1                     ! big loop forced 
    82                DO ji = jpi+2, jpij    
    83 #else 
    8480            DO jj = 2, jpj             ! no vector opt. 
    8581               DO ji = 2, jpi 
    86 #endif 
    8782#if defined key_zdfkpp 
    8883                  ! no evd mixing in the boundary layer with KPP 
     
    110105         DO jk = 1, jpkm1 
    111106!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!  
    112 #if defined key_vectopt_loop 
    113             DO jj = 1, 1                     ! big loop forced 
    114                DO ji = 1, jpij    
    115 #else 
    116107            DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call) 
    117108               DO ji = 1, jpi 
    118 #endif 
    119109#if defined key_zdfkpp 
    120110                  ! no evd mixing in the boundary layer with KPP 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r4677 r4990  
    1414   !!---------------------------------------------------------------------- 
    1515   USE par_oce         ! mesh and scale factors 
     16   USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test) 
    1617   USE ldftra_oce      ! ocean active tracers: lateral physics 
    1718   USE ldfdyn_oce      ! ocean dynamics lateral physics 
     
    117118      IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   & 
    118119         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 
     120      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 )   & 
     121         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 
    119122      ! 
    120123      !                               ! ... Convection 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90

    r4624 r4990  
    2626   USE phycst         ! physical constants 
    2727   USE eosbn2         ! equation of state 
    28    USE zdfddm         ! double diffusion mixing 
     28   USE zdfddm         ! double diffusion mixing (avs array) 
     29   USE lib_mpp        ! MPP library 
     30   USE trd_oce        ! ocean trends definition 
     31   USE trdtra         ! tracers trends 
     32   ! 
    2933   USE in_out_manager ! I/O manager 
    30    USE lib_mpp        ! MPP library 
    31    USE wrk_nemo       ! work arrays 
    3234   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3335   USE prtctl         ! Print control 
    34    USE trdmod_oce     ! ocean trends definition 
    35    USE trdtra         ! tracers trends 
     36   USE wrk_nemo       ! work arrays 
    3637   USE timing         ! Timing 
    37    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3839 
    3940   IMPLICIT NONE 
     
    246247      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t 
    247248#if defined key_zdfddm 
    248       REAL(wp) ::   zrrau, zds, zavdds, zavddt,zinr   ! double diffusion mixing 
    249       REAL(wp), POINTER, DIMENSION(:,:)   ::     zdifs 
     249      REAL(wp) ::   zrw, zkm1s                    ! local scalars 
     250      REAL(wp) ::   zrrau, zdt, zds, zavdds, zavddt, zinr   ! double diffusion mixing 
     251      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdifs 
    250252      REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps 
    251       REAL(wp) ::                            zkm1s 
    252253      REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs 
    253254      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus 
     
    274275#endif 
    275276 
    276       zviscos(:,:,:) = 0. 
    277       zblcm  (:,:  ) = 0.  
    278       zdiffut(:,:,:) = 0. 
    279       zblct  (:,:  ) = 0.  
     277      zviscos(:,:,:) = 0._wp 
     278      zblcm  (:,:  ) = 0._wp 
     279      zdiffut(:,:,:) = 0._wp 
     280      zblct  (:,:  ) = 0._wp  
    280281#if defined key_zdfddm 
    281       zdiffus(:,:,:) = 0. 
    282       zblcs  (:,:  ) = 0.  
    283 #endif 
    284       ghats(:,:,:) = 0. 
    285       
    286       zBo   (:,:) = 0. 
    287       zBosol(:,:) = 0. 
    288       zustar(:,:) = 0. 
    289  
    290  
     282      zdiffus(:,:,:) = 0._wp 
     283      zblcs  (:,:  ) = 0._wp  
     284#endif 
     285      ghats  (:,:,:) = 0._wp 
     286      zBo    (:,:  ) = 0._wp 
     287      zBosol (:,:  ) = 0._wp 
     288      zustar (:,:  ) = 0._wp 
     289      ! 
    291290      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    292291      ! I. Interior diffusivity and viscosity at w points ( T interfaces) 
     
    332331                  avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri     
    333332               ENDIF 
     333               ! 
    334334#if defined key_zdfddm  
    335                avs (ji,jj,jk) =  avt (ji,jj,jk)               
     335               ! 
    336336               !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90 
    337                ! ------------------------------------------------------------------ 
    338                ! only retains positive value of rrau 
    339                zrrau = MAX( rrau(ji,jj,jk), epsln ) 
    340                zds   = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) 
    341                IF( zrrau > 1. .AND. zds > 0.) THEN 
    342                   ! 
    343                   ! Salt fingering case. 
    344                   !--------------------- 
    345                   ! Compute interior diffusivity for double diffusive mixing of 
    346                   ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 
    347                   ! After that set interior diffusivity for double diffusive mixing 
    348                   ! of temperature 
     337               ! ------------------------- 
     338               avs (ji,jj,jk) =  avt (ji,jj,jk)    
     339 
     340               ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
     341               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     342                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     343               ! 
     344               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk) 
     345               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk) 
     346               ! 
     347               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 
     348               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) )  
     349               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp 
     350               zrrau = MAX(  epsln , zdt / zds  )    ! only retains positive value of zrau 
     351               ! 
     352               IF( zrrau > 1. .AND. zds > 0.) THEN                        ! Salt fingering case. 
     353                  !                                                       !--------------------- 
     354                  ! Compute interior diffusivity for double diffusive mixing of salinity.  
     355                  ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 
     356                  ! After that set interior diffusivity for double diffusive mixing of temperature 
    349357                  zavdds = MIN( zrrau, Rrho0 ) 
    350358                  zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 ) 
     
    353361                  zavdds = difssf * zavdds  
    354362                  zavddt = 0.7 * zavdds 
    355                ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN 
    356363                  ! 
    357                   ! Diffusive convection case. 
    358                   !--------------------------- 
    359                   ! Compute interior diffusivity for double diffusive mixing of 
    360                   ! temperature (Marmorino and Caldwell, 1976);  
     364               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN   ! Diffusive convection case. 
     365                  !                                                       !--------------------------- 
     366                  ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976);  
    361367                  ! Compute interior diffusivity for double diffusive mixing of salinity  
    362368                  zinr   = 1. / zrrau 
    363369                  zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) )  
    364370                  zavddt = difsdc * zavddt 
    365                   IF( zrrau < 0.5) THEN 
    366                      zavdds = zavddt * 0.15 * zrrau 
    367                   ELSE 
    368                      zavdds = zavddt * (1.85 * zrrau - 0.85 )  
     371                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau 
     372                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 )  
    369373                  ENDIF 
    370374               ELSE 
     
    385389      !--------------------------------------------------------------------- 
    386390      DO jj = 2, jpjm1 
    387          DO ji = fs_2, fs_jpim1      
    388             IF( nn_eos < 1) THEN    
    389                zt     = tsn(ji,jj,1,jp_tem) 
    390                zs     = tsn(ji,jj,1,jp_sal) - 35.0 
    391                zh     = fsdept(ji,jj,1) 
    392                !  potential volumic mass 
    393                zrhos  = rhop(ji,jj,1) 
    394                zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta 
    395                   &                               - 0.203814e-03 ) * zt   & 
    396                   &                               + 0.170907e-01 ) * zt   & 
    397                   &   + 0.665157e-01                                      & 
    398                   &   +     ( - 0.678662e-05 * zs                         & 
    399                   &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    400                   &   +   ( ( - 0.302285e-13 * zh                         & 
    401                   &           - 0.251520e-11 * zs                         & 
    402                   &           + 0.512857e-12 * zt * zt           ) * zh   & 
    403                   &           - 0.164759e-06 * zs                         & 
    404                   &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    405                   &                               + 0.380374e-04 ) * zh 
    406  
    407                zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta 
    408                   &                            - 0.301985e-05 ) * zt      & 
    409                   &   + 0.785567e-03                                      & 
    410                   &   + (     0.515032e-08 * zs                           & 
    411                   &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    412                   &   +(  (   0.121551e-17 * zh                           & 
    413                   &         - 0.602281e-15 * zs                           & 
    414                   &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    415                   &                             + 0.408195e-10   * zs     & 
    416                   &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    417                   &                             - 0.121555e-07 ) * zh 
    418  
    419                zthermal = zbeta * zalbet / ( rcp * zrhos + epsln ) 
    420                zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs 
    421             ELSE 
    422                zrhos    = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) ) 
    423                zthermal = rn_alpha / ( rcp * zrhos + epsln ) 
    424                zhalin   = rn_beta * tsn(ji,jj,1,jp_sal) * rcs 
    425                zbeta    = rn_beta 
    426             ENDIF 
     391         DO ji = fs_2, fs_jpim1            
     392            zrhos    = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1) 
     393            zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln ) 
     394            zbeta    = rab_n(ji,jj,1,jp_sal) 
     395            zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs 
     396            ! 
    427397            ! Radiative surface buoyancy force 
    428398            zBosol(ji,jj) = grav * zthermal * qsr(ji,jj) 
     
    435405            ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          & 
    436406               &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1)  
    437          ENDDO 
    438       ENDDO 
     407         END DO 
     408      END DO 
    439409 
    440410      zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. )  
     
    447417            ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    448418            zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) ) 
    449          ENDDO 
    450       ENDDO 
     419         END DO 
     420      END DO 
    451421 
    452422!CDIR NOVERRCHK   
     
    12701240         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    12711241!!bug gm jpttdzdf ==> jpttkpp 
    1272          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt ) 
    1273          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds ) 
     1242         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     1243         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    12741244         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    12751245      ENDIF 
     
    13401310         IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic 
    13411311            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:) 
    1342             CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:) ) 
     1312            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) ) 
    13431313         ENDIF 
    13441314         ! 
     
    13751345      !!---------------------------------------------------------------------- 
    13761346      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     1347      INTEGER  ::   ios            ! local integer 
    13771348#if ! defined key_kppcustom 
    13781349      INTEGER  ::   jm             ! dummy loop indices      
     
    13821353      REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars 
    13831354#endif 
    1384       INTEGER  ::   ios            ! Local integer output status for namelist read 
    13851355      REAL(wp) ::   zhbf           ! tempory scalars 
    13861356      LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r4245 r4990  
    66   !! History :  1.0  ! 2003-08  (G. Madec)  original code 
    77   !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop 
     8   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl  
     9   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop  
    810   !!---------------------------------------------------------------------- 
    911   !!   zdf_mxl      : Compute the turbocline and mixed layer depths. 
     
    1416   USE in_out_manager  ! I/O manager 
    1517   USE prtctl          ! Print control 
     18   USE phycst          ! physical constants 
    1619   USE iom             ! I/O library 
    1720   USE lib_mpp         ! MPP library 
     
    2528   PUBLIC   zdf_mxl       ! called by step.F90 
    2629 
    27    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    ! density criterion for mixed layer depth 
    28    REAL(wp), PUBLIC ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
    29  
    3030   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     34 
     35   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
     36   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
    3437 
    3538   !! * Substitutions 
     
    7073      !!      eddy diffusivity coefficient (resulting from the vertical physics 
    7174      !!      alone, not the isopycnal part, see trazdf.F) fall below a given 
    72       !!      value defined locally (avt_c here taken equal to 5 cm/s2) 
     75      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default) 
    7376      !! 
    7477      !! ** Action  :   nmln, hmld, hmlp, hmlpt 
    7578      !!---------------------------------------------------------------------- 
    7679      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    77       !! 
    78       INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    79       INTEGER  ::   iikn, iiki          ! temporary integer within a do loop 
    80       INTEGER, POINTER, DIMENSION(:,:) ::   imld                ! temporary workspace 
     80      ! 
     81      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     82      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer 
     83      REAL(wp) ::   zN2_c        ! local scalar 
     84      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace 
    8185      !!---------------------------------------------------------------------- 
    8286      ! 
     
    9498 
    9599      ! w-level of the mixing and mixed layers 
    96       nmln(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point 
    97       imld(:,:) = mbkt(:,:) + 1 
    98       DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
     100      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point 
     101      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
     102      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     103      DO jk = nlb10, jpkm1 
     104         DO jj = 1, jpj                ! Mixed layer level: w-level  
     105            DO ji = 1, jpi 
     106               ikt = mbkt(ji,jj) 
     107               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk) 
     108               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     109            END DO 
     110         END DO 
     111      END DO 
     112      ! 
     113      ! w-level of the turbocline 
     114      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point 
     115      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10  
    99116         DO jj = 1, jpj 
    100117            DO ji = 1, jpi 
    101                IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c )   nmln(ji,jj) = jk      ! Mixed layer 
    102                IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = jk      ! Turbocline  
     118               imkt = mikt(ji,jj) 
     119               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline  
    103120            END DO 
    104121         END DO 
     
    109126            iiki = imld(ji,jj) 
    110127            iikn = nmln(ji,jj) 
    111             hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * tmask(ji,jj,1)    ! Turbocline depth  
    112             hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * tmask(ji,jj,1)    ! Mixed layer depth 
    113             hmlpt(ji,jj) = fsdept(ji,jj,iikn-1)                     ! depth of the last T-point inside the mixed layer 
     128            imkt = mikt(ji,jj) 
     129            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth  
     130            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
     131            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    114132         END DO 
    115133      END DO 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r4624 r4990  
    241241      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242242         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     243            IF (mikt(ji,jj) .GT. 1) THEN 
     244               en(ji,jj,mikt(ji,jj))=rn_emin 
     245            ELSE 
     246               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     247            END IF 
    244248         END DO 
    245249      END DO 
     
    291295         END DO 
    292296         !                               ! finite LC depth 
    293 # if defined key_vectopt_loop 
    294          DO jj = 1, 1 
    295             DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    296 # else 
    297297         DO jj = 1, jpj  
    298298            DO ji = 1, jpi 
    299 # endif 
    300299               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 
    301300            END DO 
    302301         END DO 
    303302         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    304 !CDIR NOVERRCHK 
    305303         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    306 !CDIR NOVERRCHK 
    307304            DO jj = 2, jpjm1 
    308 !CDIR NOVERRCHK 
    309305               DO ji = fs_2, fs_jpim1   ! vector opt. 
    310306                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     
    342338      END DO 
    343339      ! 
    344       DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    345          DO jj = 2, jpjm1 
    346             DO ji = fs_2, fs_jpim1   ! vector opt. 
     340      DO jj = 2, jpjm1 
     341         DO ji = fs_2, fs_jpim1   ! vector opt. 
     342            DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
    347343               zcof   = zfact1 * tmask(ji,jj,jk) 
    348344               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
     
    360356               !                                   ! right hand side in en 
    361357               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    362                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk) 
    363             END DO 
    364          END DO 
    365       END DO 
    366       !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    367       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    368          DO jj = 2, jpjm1 
    369             DO ji = fs_2, fs_jpim1    ! vector opt. 
     358                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     359                  &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     360            END DO 
     361            !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     362            DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    370363               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    371364            END DO 
    372          END DO 
    373       END DO 
    374       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    375          DO ji = fs_2, fs_jpim1    ! vector opt. 
    376             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    377          END DO 
    378       END DO 
    379       DO jk = 3, jpkm1 
    380          DO jj = 2, jpjm1 
    381             DO ji = fs_2, fs_jpim1    ! vector opt. 
     365            ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     366            zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
     367            ! 
     368            DO jk = mikt(ji,jj)+2, jpkm1 
    382369               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    383370            END DO 
    384          END DO 
    385       END DO 
    386       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    387          DO ji = fs_2, fs_jpim1    ! vector opt. 
     371            ! 
     372            ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    388373            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    389          END DO 
    390       END DO 
    391       DO jk = jpk-2, 2, -1 
    392          DO jj = 2, jpjm1 
    393             DO ji = fs_2, fs_jpim1    ! vector opt. 
     374            ! 
     375            DO jk = jpk-2, mikt(ji,jj)+1, -1 
    394376               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    395377            END DO 
    396          END DO 
    397       END DO 
    398       DO jk = 2, jpkm1                             ! set the minimum value of tke 
    399          DO jj = 2, jpjm1 
    400             DO ji = fs_2, fs_jpim1   ! vector opt. 
     378            ! 
     379            DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    401380               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
    402381            END DO 
     
    412391               DO ji = fs_2, fs_jpim1   ! vector opt. 
    413392                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    414                      &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     393                     &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    415394               END DO 
    416395            END DO 
     
    421400               jk = nmln(ji,jj) 
    422401               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                   &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     402                  &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    424403            END DO 
    425404         END DO 
     
    433412                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    434413                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    435                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress  
     414                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    436415                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    437416                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    438417                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    439                      &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 
     418                     &                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
    440419               END DO 
    441420            END DO 
     
    506485      ! 
    507486      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    508          zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    509          zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  ) 
    510       ELSE                          ! surface set to the minimum value 
    511          zmxlm(:,:,1) = rn_mxl0 
     487         DO jj = 2, jpjm1 
     488            DO ji = fs_2, fs_jpim1 
     489               IF (mikt(ji,jj) .GT. 1) THEN 
     490                  zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
     491               ELSE 
     492                  zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     493                  zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
     494               END IF 
     495            END DO 
     496         END DO 
     497      ELSE  
     498         DO jj = 2, jpjm1 
     499            DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
     500               zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
     501            END DO 
     502         END DO 
    512503      ENDIF 
    513504      zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    514505      ! 
    515506!CDIR NOVERRCHK 
    516       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    517 !CDIR NOVERRCHK 
    518          DO jj = 2, jpjm1 
    519 !CDIR NOVERRCHK 
    520             DO ji = fs_2, fs_jpim1   ! vector opt. 
     507      DO jj = 2, jpjm1 
     508!CDIR NOVERRCHK 
     509         DO ji = fs_2, fs_jpim1   ! vector opt. 
     510            !CDIR NOVERRCHK 
     511            DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    521512               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    522513               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    523514            END DO 
     515            zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
    524516         END DO 
    525517      END DO 
     
    533525      ! 
    534526      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    535          DO jk = 2, jpkm1 
    536             DO jj = 2, jpjm1 
    537                DO ji = fs_2, fs_jpim1   ! vector opt. 
    538                   zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   & 
     527         DO jj = 2, jpjm1 
     528            DO ji = fs_2, fs_jpim1   ! vector opt. 
     529               DO jk = mikt(ji,jj)+1, jpkm1 
     530                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    539531                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    540532                  zmxlm(ji,jj,jk) = zemxl 
     
    545537         ! 
    546538      CASE ( 1 )           ! bounded by the vertical scale factor 
    547          DO jk = 2, jpkm1 
    548             DO jj = 2, jpjm1 
    549                DO ji = fs_2, fs_jpim1   ! vector opt. 
     539         DO jj = 2, jpjm1 
     540            DO ji = fs_2, fs_jpim1   ! vector opt. 
     541               DO jk = mikt(ji,jj)+1, jpkm1 
    550542                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    551543                  zmxlm(ji,jj,jk) = zemxl 
     
    556548         ! 
    557549      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    558          DO jk = 2, jpkm1         ! from the surface to the bottom : 
    559             DO jj = 2, jpjm1 
    560                DO ji = fs_2, fs_jpim1   ! vector opt. 
     550         DO jj = 2, jpjm1 
     551            DO ji = fs_2, fs_jpim1   ! vector opt. 
     552               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
    561553                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    562554               END DO 
    563             END DO 
    564          END DO 
    565          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
    566             DO jj = 2, jpjm1 
    567                DO ji = fs_2, fs_jpim1   ! vector opt. 
     555               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
    568556                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    569557                  zmxlm(ji,jj,jk) = zemxl 
     
    574562         ! 
    575563      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    576          DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
    577             DO jj = 2, jpjm1 
    578                DO ji = fs_2, fs_jpim1   ! vector opt. 
     564         DO jj = 2, jpjm1 
     565            DO ji = fs_2, fs_jpim1   ! vector opt. 
     566               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
    579567                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    580568               END DO 
    581             END DO 
    582          END DO 
    583          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
    584             DO jj = 2, jpjm1 
    585                DO ji = fs_2, fs_jpim1   ! vector opt. 
     569               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
    586570                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    587571               END DO 
     
    628612      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    629613      ! 
    630       DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points 
     614      DO jj = 2, jpjm1 
     615         DO ji = fs_2, fs_jpim1   ! vector opt. 
     616            DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
     617               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
     618            END DO 
     619            DO jk = mikv(ji,jj)+1, jpkm1 
     620               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     621            END DO 
     622         END DO 
     623      END DO 
     624      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     625      ! 
     626      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    631627         DO jj = 2, jpjm1 
    632628            DO ji = fs_2, fs_jpim1   ! vector opt. 
    633                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    634                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
    635             END DO 
    636          END DO 
    637       END DO 
    638       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    639       ! 
    640       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    641          DO jk = 2, jpkm1 
    642             DO jj = 2, jpjm1 
    643                DO ji = fs_2, fs_jpim1   ! vector opt. 
     629               DO jk = mikt(ji,jj)+1, jpkm1 
    644630                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    645631                  !                                          ! shear 
     
    655641                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    656642# if defined key_c1d 
    657                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number 
    658                   e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
     643                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     644                  e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
    659645# endif 
    660646              END DO 
     
    740726      ! 
    741727      !                               !* Check of some namelist values 
    742       IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    743       IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    744       IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
    745 #if ! key_coupled 
    746       IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    747 #endif 
     728      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
     729      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
     730      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     731      IF( nn_etau == 3 .AND. .NOT. lk_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    748732 
    749733      IF( ln_mxl0 ) THEN 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r4624 r4990  
    126126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column 
    127127      DO jk = 2, jpkm1 
    128          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) 
     128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 
    129129      END DO 
    130130 
     
    135135      END DO 
    136136 
    137       DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
    138          zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s 
     137      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     138         DO ji = 1, jpi 
     139            DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
     140               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     141            END DO 
     142         END DO 
    139143      END DO 
    140144 
     
    162166      !                          !   Update  mixing coefs  !                           
    163167      !                          ! ----------------------- ! 
    164       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    165          avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk) 
    166          avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) 
    167          DO jj = 2, jpjm1 
    168             DO ji = fs_2, fs_jpim1  ! vector opt. 
     168      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     169         DO ji = 1, jpi 
     170            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     171               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 
     172               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 
     173            END DO 
     174         END DO 
     175      END DO 
     176       
     177      DO jj = 2, jpjm1 
     178         DO ji = fs_2, fs_jpim1  ! vector opt. 
     179            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    169180               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    170181               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     
    237248      zsum2(:,:) = 0.e0 
    238249      DO jk= 2, jpk 
    239          zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 
    240          zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                 
     250         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
     251         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)                
    241252      END DO 
    242253      DO jj = 1, jpj 
     
    274285      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column 
    275286      DO jk = 2, jpkm1 
    276          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) 
     287         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
    277288      END DO 
    278289 
     
    284295 
    285296      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s 
    286          zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s 
     297         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1)   ! kz max = 120 cm2/s 
    287298      END DO 
    288299 
     
    414425      ! only the energy available for mixing is taken into account, 
    415426      ! (mixing efficiency tidal dissipation efficiency) 
    416       en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1) 
     427      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:) 
    417428 
    418429      ! Vertical structure (az_tmx) 
     
    443454         ztpc = 0.e0 
    444455         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
    445          DO jk= 2, jpkm1 
    446             DO jj = 1, jpj 
    447                DO ji = 1, jpi 
     456         DO jj = 1, jpj 
     457            DO ji = 1, jpi 
     458               DO jk= mikt(ji,jj)+1, jpkm1 
    448459                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
    449460               END DO 
     
    459470         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )    
    460471         zkz(:,:) = 0.e0 
    461          DO jk = 2, jpkm1 
    462          DO jj = 1, jpj 
    463             DO ji = 1, jpi 
     472         DO jj = 1, jpj 
     473            DO ji = 1, jpi 
     474               DO jk = mikt(ji,jj)+1, jpkm1 
    464475               zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk) 
    465             END DO 
    466          END DO 
     476               END DO 
     477            END DO 
    467478         END DO 
    468479         ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz 
     
    484495         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 
    485496 
    486          DO jk = 2, jpkm1 
    487             zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s 
     497         DO jj = 1, jpj 
     498            DO ji = 1, jpi 
     499               DO jk = mikt(ji,jj)+1, jpkm1 
     500                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     501               END DO 
     502            END DO 
    488503         END DO 
    489504         ztpc = 0.e0 
Note: See TracChangeset for help on using the changeset viewer.