New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 2371 for branches/nemo_v3_3_beta/NEMOGCM/NEMO – NEMO

Ignore:
Timestamp:
2010-11-10T16:38:27+01:00 (13 years ago)
Author:
acc
Message:

nemo_v3_3_beta. Improvement of the Griffies operator code (#680). Some aspects are still undergoing scientific assessment but the code structure is ready for release. Dynamical allocation of the large triad arrays has been introduced to reduce memory when the new operator is not in use.

Location:
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r2287 r2371  
    8787      ! 1. Compute lateral diffusive coefficient  
    8888      ! ---------------------------------------- 
    89  
    90       DO jk = 1, jpk 
     89      IF (ln_traldf_grif) THEN 
     90         DO jk = 1, jpk 
    9191#  if defined key_vectopt_loop   
    92 !CDIR NOVERRCHK  
    93          DO ji = 1, jpij   ! vector opt. 
    94             ! Take the max of N^2 and zero then take the vertical sum 
    95             ! of the square root of the resulting N^2 ( required to compute 
    96             ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
    97             zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
    98             zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
    99             ! Compute elements required for the inverse time scale of baroclinic 
    100             ! eddies using the isopycnal slopes calculated in ldfslp.F : 
    101             ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    102             ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk) 
     92            !CDIR NOVERRCHK  
     93            DO ji = 1, jpij   ! vector opt. 
     94               ! Take the max of N^2 and zero then take the vertical sum 
     95               ! of the square root of the resulting N^2 ( required to compute 
     96               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
     97               zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
     98               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
     99               ! Compute elements required for the inverse time scale of baroclinic 
     100               ! eddies using the isopycnal slopes calculated in ldfslp.F : 
     101               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     102               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk) 
     103               zah(ji,1) = zah(ji,1) + zn2 * wslp2(ji,1,jk) * ze3w 
     104               zhw(ji,1) = zhw(ji,1) + ze3w 
     105            END DO 
     106#  else 
     107            DO jj = 2, jpjm1 
     108               !CDIR NOVERRCHK  
     109               DO ji = 2, jpim1 
     110                  ! Take the max of N^2 and zero then take the vertical sum  
     111                  ! of the square root of the resulting N^2 ( required to compute  
     112                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
     113                  zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
     114                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
     115                  ! Compute elements required for the inverse time scale of baroclinic 
     116                  ! eddies using the isopycnal slopes calculated in ldfslp.F :  
     117                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     118                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     119                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 
     120                  zhw(ji,jj) = zhw(ji,jj) + ze3w 
     121               END DO 
     122            END DO 
     123#  endif 
     124         END DO 
     125      ELSE 
     126         DO jk = 1, jpk 
     127#  if defined key_vectopt_loop   
     128            !CDIR NOVERRCHK  
     129            DO ji = 1, jpij   ! vector opt. 
     130               ! Take the max of N^2 and zero then take the vertical sum 
     131               ! of the square root of the resulting N^2 ( required to compute 
     132               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
     133               zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
     134               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
     135               ! Compute elements required for the inverse time scale of baroclinic 
     136               ! eddies using the isopycnal slopes calculated in ldfslp.F : 
     137               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     138               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk) 
    103139               zah(ji,1) = zah(ji,1) + zn2   & 
    104                               * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    & 
    105                                 + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   & 
    106                               * ze3w 
    107             zhw(ji,1) = zhw(ji,1) + ze3w 
    108          END DO 
     140                    * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    & 
     141                    + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   & 
     142                    * ze3w 
     143               zhw(ji,1) = zhw(ji,1) + ze3w 
     144            END DO 
    109145#  else 
    110          DO jj = 2, jpjm1 
    111 !CDIR NOVERRCHK  
    112             DO ji = 2, jpim1 
    113                ! Take the max of N^2 and zero then take the vertical sum  
    114                ! of the square root of the resulting N^2 ( required to compute  
    115                ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    116                zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
    117                zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
    118                ! Compute elements required for the inverse time scale of baroclinic 
    119                ! eddies using the isopycnal slopes calculated in ldfslp.F :  
    120                ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    121                ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
    122                zah(ji,jj) = zah(ji,jj) + zn2   & 
    123                               * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    & 
    124                                 + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  & 
    125                               * ze3w 
    126                zhw(ji,jj) = zhw(ji,jj) + ze3w 
    127             END DO  
    128          END DO  
     146            DO jj = 2, jpjm1 
     147               !CDIR NOVERRCHK  
     148               DO ji = 2, jpim1 
     149                  ! Take the max of N^2 and zero then take the vertical sum  
     150                  ! of the square root of the resulting N^2 ( required to compute  
     151                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
     152                  zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
     153                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
     154                  ! Compute elements required for the inverse time scale of baroclinic 
     155                  ! eddies using the isopycnal slopes calculated in ldfslp.F :  
     156                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     157                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     158                  zah(ji,jj) = zah(ji,jj) + zn2   & 
     159                       * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    & 
     160                       + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  & 
     161                       * ze3w 
     162                  zhw(ji,jj) = zhw(ji,jj) + ze3w 
     163               END DO 
     164            END DO 
    129165#  endif 
    130       END DO  
     166         END DO 
     167      END IF 
    131168 
    132169      DO jj = 2, jpjm1 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2344 r2371  
    99   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90 
    1010   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates 
    11    !!            3.3  ! 2006-10  (C. Harris, G. Nurser)  add ldf_slp_grif (Griffies operator) 
     11   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator 
    1212   !!---------------------------------------------------------------------- 
    1313#if   defined key_ldfslp   ||   defined key_esopa 
    … …  
    1515   !!   'key_ldfslp'                      Rotation of lateral mixing tensor 
    1616   !!---------------------------------------------------------------------- 
    17    !!   ldf_slp      : compute the slopes of neutral surface 
    18    !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface 
     17   !!   ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator) 
     18   !!   ldf_slp      : calculates the slopes of neutral surface   (Madec operator) 
     19   !!   ldf_slp_mxl  : calculates the slopes at the base of the mixed layer (Madec operator) 
    1920   !!   ldf_slp_init : initialization of the slopes computation 
    2021   !!---------------------------------------------------------------------- 
    2122   USE oce             ! ocean dynamics and tracers 
    2223   USE dom_oce         ! ocean space and time domain 
    23    USE ldftra_oce 
    24    USE ldfdyn_oce 
     24   USE ldftra_oce      ! lateral diffusion: traceur 
     25   USE ldfdyn_oce      ! lateral diffusion: dynamics 
    2526   USE phycst          ! physical constants 
    2627   USE zdfmxl          ! mixed layer depth 
    27    USE eosbn2 
     28   USE eosbn2          ! equation of states 
    2829   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2930   USE in_out_manager  ! I/O manager 
    … …  
    3334   PRIVATE 
    3435 
    35    PUBLIC   ldf_slp         ! routine called by step.F90 
    36    PUBLIC   ldf_slp_init    ! routine called by opa.F90 
    37    PUBLIC   ldf_slp_grif    ! routine called by step.F90 
    38  
    39    LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag 
    40    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points 
    41    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points 
    42    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   wslp2                !: wslp**2 from Griffies quarter cells 
    43    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   alpha, beta          !: alpha,beta at T points 
    44    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd 
    45    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psix_eiv 
    46    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psiy_eiv 
    47     
    48    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt 
    49    REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer 
    50    REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer 
     36   PUBLIC   ldf_slp        ! routine called by step.F90 
     37   PUBLIC   ldf_slp_grif   ! routine called by step.F90 
     38   PUBLIC   ldf_slp_init   ! routine called by opa.F90 
     39 
     40   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag 
     41   !                                                                             !! Madec operator 
     42   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   uslp, wslpi          !: i_slope at U- and W-points 
     43   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   vslp, wslpj          !: j-slope at V- and W-points 
     44   !                                                                             !! Griffies operator 
     45   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   wslp2                !: wslp**2 from Griffies quarter cells 
     46   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials  
     47   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate 
     48 
     49   !                                                              !! Madec operator 
     50   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   omlmask           ! mask of the surface mixed layer at T-pt 
     51   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer 
     52   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer 
     53 
     54   REAL(wp) ::   repsln = 1.e-25_wp       ! tiny value used as minium of di(rho), dj(rho) and dk(rho) 
    5155 
    5256   !! * Substitutions 
    … …  
    6771      !!  
    6872      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal 
    69       !!              surfaces referenced locally) ('key_traldfiso'). 
     73      !!              surfaces referenced locally) (ln_traldf_iso=T). 
    7074      !!  
    7175      !! ** Method  :   The slope in the i-direction is computed at U- and  
    … …  
    347351      ! 
    348352   END SUBROUTINE ldf_slp 
    349  
     353    
    350354 
    351355   SUBROUTINE ldf_slp_grif ( kt ) 
    352      !!---------------------------------------------------------------------- 
    353      !!                 ***  ROUTINE ldf_slp_grif  *** 
    354      !! 
    355      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope 
    356      !!      of iso-pycnal surfaces referenced locally) ('key_traldfiso') 
    357      !!      at W-points using the Griffies quarter-cells.  Also calculates 
    358      !!      alpha and beta at T-points for use in the Griffies isopycnal 
    359      !!      scheme. 
    360      !! 
    361      !! ** Method  :   The slope in the i-direction is computed at U- and 
    362      !!      W-points (uslp, wslpi) and the slope in the j-direction is 
    363      !!      computed at V- and W-points (vslp, wslpj). 
    364      !! 
    365      !! ** Action : - alpha, beta 
    366      !!               wslp2 squared slope of neutral surfaces at w-points. 
    367      !!---------------------------------------------------------------------- 
    368      USE oce, zdit  => ua   ! use ua as workspace 
    369      USE oce, zdis  => va   ! use va as workspace 
    370      USE oce, zdjt  => ta   ! use ta as workspace 
    371      USE oce, zdjs  => sa   ! use sa as workspace 
    372      !! 
    373      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
    374      !! 
    375      INTEGER  ::   ji, jj, jk, ip, jp, kp  ! dummy loop indices 
    376      INTEGER  ::   iku, ikv                ! local integer 
    377      REAL(wp) ::   zt, zs, zh, zt2, zsp5, zp1t1           ! local scalars 
    378      REAL(wp) ::     zdenr, zrhotmp, zdndt, zdddt         !   -      - 
    379      REAL(wp) ::     zdnds, zddds, znum, zden             !   -      - 
    380      REAL(wp) ::     zslope, za_sxe, zslopec, zdsloper    !   -      - 
    381      REAL(wp) ::     zfact, zepsln, zatempw,zatempu,zatempv    !   -      - 
    382      REAL(wp) ::     ze1ur, ze2vr, ze3wr, zdxt, zdxs, zdyt, zdys, zdzt, zdzs, zvolf  
    383      REAL(wp) ::     zr_slpmax, zdxrho, zdyrho, zabs_dzrho 
    384      REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) ::   zsx,zsy 
    385      REAL(wp), DIMENSION(jpi,jpj    ,0:1,0:1) ::   zsx_ml_base, zsy_ml_base 
    386      REAL(wp), DIMENSION(jpi,jpj,jpk)         ::   zdkt, zdks 
    387      REAL(wp), DIMENSION(jpi,jpj) ::   zr_ml_basew 
    388      !!---------------------------------------------------------------------- 
    389  
    390      ! 0. Local constant initialization 
    391      ! -------------------------------- 
    392      zr_slpmax = 1.0_wp/slpmax 
    393  
    394      ! zslopec=0.004 
    395      ! zdsloper=1000.0 
    396      zepsln=1e-25 
    397  
    398      SELECT CASE ( nn_eos ) 
    399  
    400      CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
    401  
    402         !                                                ! =============== 
    403         DO jk = 1, jpk                                   ! Horizontal slab 
    404            !                                             ! =============== 
    405            DO jj = 1, jpjm1 
    406               DO ji = 1, fs_jpim1 
    407                  zt = tb(ji,jj,jk)          ! potential temperature 
    408                  zs = sb(ji,jj,jk) - 35.0   ! salinity anomaly (s-35) 
    409                  zh = fsdept(ji,jj,jk)      ! depth in meters 
    410  
    411                  beta(ji,jj,jk) = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      & 
    412                       &                            - 0.301985e-05 ) * zt      & 
    413                       &   + 0.785567e-03                                      & 
    414                       &   + (     0.515032e-08 * zs                           & 
    415                       &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    416                       &   +(  (   0.121551e-17 * zh                           & 
    417                       &         - 0.602281e-15 * zs                           & 
    418                       &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    419                       &                             + 0.408195e-10   * zs     & 
    420                       &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    421                       &                             - 0.121555e-07 ) * zh 
    422  
    423                  alpha(ji,jj,jk) = - beta(ji,jj,jk) *                       & 
    424                       &     (((( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   & 
    425                       &                               - 0.203814e-03 ) * zt   & 
    426                       &                               + 0.170907e-01 ) * zt   & 
    427                       &   + 0.665157e-01                                      & 
    428                       &   +     ( - 0.678662e-05 * zs                         & 
    429                       &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    430                       &   +   ( ( - 0.302285e-13 * zh                         & 
    431                       &           - 0.251520e-11 * zs                         & 
    432                       &           + 0.512857e-12 * zt * zt           ) * zh   & 
    433                       &           - 0.164759e-06 * zs                         & 
    434                       &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    435                       &                               + 0.380374e-04 ) * zh) 
    436  
    437               END DO 
    438            END DO 
    439         END DO 
    440  
    441      CASE ( 1 ) 
    442  
    443         alpha(:,:,:)=-rn_alpha 
    444         beta(:,:,:)=0.0 
    445  
    446      CASE ( 2 ) 
    447  
    448         alpha(:,:,:)=-rn_alpha 
    449         beta (:,:,:)=rn_beta 
    450  
    451      CASE ( 3 ) 
    452  
    453         DO jk =1, jpk 
    454            DO jj = 1, jpjm1 
    455               DO ji = 1, fs_jpim1 
    456                  zt = tb(ji,jj,jk) 
    457                  zs = sb(ji,jj,jk) 
    458                  zh = fsdept(ji,jj,jk) 
    459                  zt2 = zt**2 
    460                  zsp5 = sqrt(ABS(zs)) 
    461                  zp1t1=zh*zt 
    462                  znum=9.99843699e+02+zt*(7.35212840e+00+zt*(-5.45928211e-02+3.98476704e-04*zt)) & 
    463                       +zs*(2.96938239e+00-7.23268813e-03*zt+2.12382341e-03*zs)  & 
    464                       +zh*(1.04004591e-02+1.03970529e-07*zt2+5.18761880e-06*zs+ & 
    465                       zh*(-3.24041825e-08-1.23869360e-11*zt2)) 
    466                  zden=1.00000000e+00+zt*(7.28606739e-03+zt*(-4.60835542e-05+zt*(3.68390573e-07+zt*1.80809186e-10))) & 
    467                       +zs*(2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(4.76534122e-06+1.63410736e-09*zt2)) & 
    468                       + zh*(5.30848875e-06+zh*zt*(-3.03175128e-16*zt2-1.27934137e-17*zh)) 
    469                  zdenr=1.0/zden 
    470                  zrhotmp=znum*zdenr 
    471                  zdndt=7.35212840e+00+zt*(-1.091856422e-01+1.195430112e-03*zt)-7.23268813e-03*zs & 
    472                       +zp1t1*(2.07941058e-07-2.4773872e-11*zh) 
    473                  zdddt=7.28606739e-03+zt*(-9.21671084e-05+zt*(1.105171719e-06+7.23236744e-10*zt))  & 
    474                       +zs*(-9.27062484e-06+zt*(-5.35030929e-10*zt+3.26821472e-09*zsp5))  & 
    475                       +zh*zh*(-9.09525384e-16*zt2-1.27934137e-17*zh) 
    476                  zdnds=2.96938239e+00-7.23268813e-03*zt+2*2.12382341e-03*zs+5.18761880e-06*zh 
    477                  zddds=2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(7.14801183e-06+2.45116104e-09*zt2) 
    478                  alpha(ji,jj,jk)=(zdndt-zrhotmp*zdddt)*zdenr 
    479                  beta(ji,jj,jk)=zdenr*(zdnds-zrhotmp*zddds) 
    480  
    481               END DO 
    482            END DO 
    483         END DO 
    484  
    485      CASE DEFAULT 
    486  
    487         IF(lwp) WRITE(numout,cform_err) 
    488         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
    489         nstop = nstop + 1 
    490  
    491      END SELECT 
    492  
    493      CALL lbc_lnk( alpha, 'T', 1. ) 
    494      CALL lbc_lnk( beta, 'T', 1. ) 
    495  
    496  
    497      ! --------------------------------------- 
    498      ! 1. Horizontal tracer gradients at T-level jk 
    499      ! --------------------------------------- 
    500      DO jk = 1, jpkm1 
    501         DO jj = 1, jpjm1 
    502            DO ji = 1, fs_jpim1   ! vector opt. 
    503               ! i-gradient of T and S at jj 
    504               zdit (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk) 
    505               zdis (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk) 
    506               ! j-gradient of T and S at jj 
    507               zdjt (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk) 
    508               zdjs (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
    509            END DO 
    510         END DO 
    511      END DO 
    512  
    513      IF( ln_zps ) THEN      ! partial steps correction at the last level 
     356      !!---------------------------------------------------------------------- 
     357      !!                 ***  ROUTINE ldf_slp_grif  *** 
     358      !! 
     359      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope 
     360      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 
     361      !!      at W-points using the Griffies quarter-cells.   
     362      !! 
     363      !! ** Method  :   calculates alpha and beta at T-points  
     364      !! 
     365      !! ** Action : - triadi_g, triadj_g   T-pts i- and j-slope triads relative to geopot. (used for eiv) 
     366      !!             - triadi , triadj    T-pts i- and j-slope triads relative to model-coordinate 
     367      !!             - wslp2              squared slope of neutral surfaces at w-points. 
     368      !!---------------------------------------------------------------------- 
     369      USE oce,   zdit  => ua   ! use ua as workspace 
     370      USE oce,   zdis  => va   ! use va as workspace 
     371      USE oce,   zdjt  => ta   ! use ta as workspace 
     372      USE oce,   zdjs  => sa   ! use sa as workspace 
     373      !! 
     374      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
     375      !! 
     376      INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices 
     377      INTEGER  ::   iku, ikv                ! temporary integer 
     378      REAL(wp) ::   zfacti, zfactj, zatempw,zatempu,zatempv   ! local scalars 
     379      REAL(wp) ::   zbu, zbv, zbti, zbtj 
     380      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 
     381      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 
     382      REAL(wp) ::   zdzrho_raw 
     383      REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) ::   zdzrho, zdyrho, zdxrho     ! Horizontal and vertical density gradients 
     384      REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   zti_mlb, ztj_mlb 
     385      REAL(wp), DIMENSION(jpi,jpj,jpk)     ::   zdkt, zdks 
     386      REAL(wp), DIMENSION(jpi,jpj,jpk)     ::   zalpha, zbeta       ! alpha, beta at T points, at depth fsgdept 
     387      REAL(wp), DIMENSION(jpi,jpj)         ::   z1_mlbw 
     388      !!---------------------------------------------------------------------- 
     389 
     390      !--------------------------------! 
     391      !  Some preliminary calculation  ! 
     392      !--------------------------------! 
     393      ! 
     394      CALL eos_alpbet( tsb, zalpha, zbeta )     !==  before thermal and haline expension coeff. at T-points  ==! 
     395      ! 
     396      DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==! 
     397         DO jj = 1, jpjm1 
     398            DO ji = 1, fs_jpim1   ! vector opt. 
     399               zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk)    ! i-gradient of T and S at jj 
     400               zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 
     401               zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk)    ! j-gradient of T and S at jj 
     402               zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
     403            END DO 
     404         END DO 
     405      END DO 
     406      IF( ln_zps ) THEN                               ! partial steps: correction at the last level 
    514407# if defined key_vectopt_loop 
    515      DO jj = 1, 1 
    516         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     408         DO jj = 1, 1 
     409            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    517410# else 
    518            DO jj = 1, jpjm1 
    519               DO ji = 1, jpim1 
     411         DO jj = 1, jpjm1 
     412            DO ji = 1, jpim1 
    520413# endif 
    521                  ! last ocean level 
    522                  iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
    523                  ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
    524                  ! i-gradient of T and S 
    525                  zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem) 
    526                  zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 
    527                  ! j-gradient of T and S 
    528                  zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem) 
    529                  zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 
    530               END DO 
    531            END DO 
    532         ENDIF 
    533  
    534         ! --------------------------------------- 
    535         ! 1. Vertical tracer gradient at w-level jk 
    536         ! --------------------------------------- 
    537         DO jk = 2, jpk 
    538            zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
    539            zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
    540         END DO 
    541  
    542         zdkt(:,:,1) = 0.0 
    543         zdks(:,:,1) = 0.0 
    544  
    545         ! --------------------------------------- 
    546         ! Depth of the w-point below ML base 
    547         ! --------------------------------------- 
    548         DO jj = 1, jpj 
    549            DO ji = 1, jpi 
    550               jk = nmln(ji,jj) 
    551               zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1) 
    552            END DO 
    553         END DO 
    554  
    555  
    556         wslp2(:,:,:)=0.0 
    557         tfw(:,:,:) = 0.0 
    558         sfw(:,:,:) = 0.0 
    559         ftu(:,:,:) = 0.0 
    560         fsu(:,:,:) = 0.0 
    561         ftv(:,:,:) = 0.0 
    562         fsv(:,:,:) = 0.0 
    563         ftud(:,:,:) = 0.0 
    564         fsud(:,:,:) = 0.0 
    565         ftvd(:,:,:) = 0.0 
    566         fsvd(:,:,:) = 0.0 
    567         psix_eiv(:,:,:) = 0.0 
    568         psiy_eiv(:,:,:) = 0.0 
    569  
    570         ! ---------------------------------------------------------------------- 
    571         ! x-z plane 
    572         ! ---------------------------------------------------------------------- 
    573  
    574         ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1) 
    575         DO ip=0,1 
    576            DO kp=0,1 
    577  
    578               DO jk = 1, jpkm1 
    579                  DO jj = 1, jpjm1 
    580                     DO ji = 1, fs_jpim1   ! vector opt. 
    581  
    582                        ze1ur=1.0/e1u(ji,jj) 
    583                        zdxt=zdit(ji,jj,jk)*ze1ur 
    584                        zdxs=zdis(ji,jj,jk)*ze1ur 
    585  
    586                        ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 
    587                        zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
    588                        zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
    589                        ! Calculate the horizontal and vertical density gradient 
    590                        zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs 
    591                        zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln 
    592                        ! Limit by slpmax, and mask by psi-point 
    593                        zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) & 
    594                             & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho)) 
    595                     END DO 
    596                  END DO 
    597               END DO 
    598  
    599            END DO 
    600         END DO 
    601  
    602         ! calculate slope of triad immediately below mixed-layer base 
    603         DO ip=0,1 
    604            DO kp=0,1 
    605               DO jj = 1, jpjm1 
    606                  DO ji = 1, fs_jpim1 
    607                     jk = nmln(ji+ip,jj) 
    608                     zsx_ml_base(ji+ip,jj,1-ip,kp)=zsx(ji+ip,jj,jk+1-kp,1-ip,kp) 
    609                  END DO 
    610               END DO 
    611            END DO 
    612         END DO 
    613  
    614         ! Below ML use limited zsx as is 
    615         ! Inside ML replace by linearly reducing zsx_ml_base towards surface 
    616         DO ip=0,1 
    617            DO kp=0,1 
    618  
    619               DO jk = 1, jpkm1 
    620  
    621                  DO jj = 1, jpjm1 
    622  
    623                     DO ji = 1, fs_jpim1   ! vector opt. 
    624                        ! k index of uppermost point(s) of triad is jk+kp-1 
    625                        ! must be .ge. nmln(ji,jj) for zfact=1. 
    626                        !                   otherwise  zfact=0. 
    627                        zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)) 
    628                        zsx(ji+ip,jj,jk,1-ip,kp) = zfact*zsx(ji+ip,jj,jk,1-ip,kp) + & 
    629                             & (1.0_wp-zfact)*(fsdepw(ji+ip,jj,jk+kp)*zr_ml_basew(ji+ip,jj))*zsx_ml_base(ji+ip,jj,1-ip,kp)  
    630                     END DO 
    631  
    632                  END DO 
    633  
    634               END DO 
    635            END DO 
    636         END DO 
    637  
    638         ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points 
    639         DO ip=0,1 
    640            DO kp=0,1 
    641  
    642               DO jk = 1, jpkm1 
    643  
    644                  DO jj = 1, jpjm1 
    645  
    646                     DO ji = 1, fs_jpim1   ! vector opt. 
    647  
    648                        ze1ur=1.0/e1u(ji,jj) 
    649                        zdxt=zdit(ji,jj,jk)*ze1ur 
    650                        zdxs=zdis(ji,jj,jk)*ze1ur 
    651  
    652                        ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 
    653                        zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
    654                        zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
    655                        zslope=zsx(ji+ip,jj,jk,1-ip,kp) 
    656  
    657                        zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
    658  
    659                        ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur 
    660                        fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur 
    661                        ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur 
    662                        fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur 
    663                        tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt 
    664                        sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs 
    665                        wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ & 
    666                             & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2 
    667                        psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 
    668  
    669                     END DO 
    670                  END DO 
    671  
    672               END DO 
    673            END DO 
    674         END DO 
    675  
    676         ! ---------------------------------------------------------------------- 
    677         ! y-z plane 
    678         ! ---------------------------------------------------------------------- 
    679  
    680         ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1) 
    681         DO jp=0,1 
    682            DO kp=0,1 
    683  
    684               DO jk = 1, jpkm1 
    685                  DO jj = 1, jpjm1 
    686                     DO ji = 1, fs_jpim1   ! vector opt. 
    687  
    688                        ze2vr=1.0/e2v(ji,jj) 
    689                        zdyt=zdjt(ji,jj,jk)*ze2vr 
    690                        zdys=zdjs(ji,jj,jk)*ze2vr 
    691  
    692                        ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 
    693                        zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
    694                        zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
    695                        ! Calculate the horizontal and vertical density gradient 
    696                        zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys 
    697                        zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln 
    698                        ! Limit by slpmax, and mask by psi-point 
    699                        zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) & 
    700                             & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho)) 
    701                     END DO 
    702                  END DO 
    703               END DO 
    704  
    705            END DO 
    706         END DO 
    707  
    708         ! calculate slope of triad immediately below mixed-layer base 
    709         DO jp=0,1 
    710            DO kp=0,1 
    711               DO jj = 1, jpjm1 
    712                  DO ji = 1, fs_jpim1 
    713                     jk = nmln(ji,jj+jp) 
    714                     zsy_ml_base(ji,jj+jp,1-jp,kp)=zsy(ji,jj+jp,jk+1-kp,1-jp,kp) 
    715                  END DO 
    716               END DO 
    717            END DO 
    718         END DO 
    719  
    720         ! Below ML use limited zsy as is 
    721         ! Inside ML replace by linearly reducing zsy_ml_base towards surface 
    722         DO jp=0,1 
    723            DO kp=0,1 
    724               DO jk = 1, jpkm1 
    725                  DO jj = 1, jpjm1 
    726                     DO ji = 1, fs_jpim1   ! vector opt. 
    727                        ! k index of uppermost point(s) of triad is jk+kp-1 
    728                        ! must be .ge. nmln(ji,jj) for zfact=1. 
    729                        !                   otherwise  zfact=0. 
    730                        zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)) 
    731                        zsy(ji,jj+jp,jk,1-jp,kp) = zfact*zsy(ji,jj+jp,jk,1-jp,kp) + & 
    732                            & (1.0_wp-zfact)*(fsdepw(ji,jj+jp,jk+kp)*zr_ml_basew(ji,jj+jp))*zsy_ml_base(ji,jj+jp,1-jp,kp)  
    733                     END DO 
    734  
    735                  END DO 
    736  
    737               END DO 
    738            END DO 
    739         END DO 
    740  
    741         ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points 
    742         DO jp=0,1 
    743            DO kp=0,1 
    744               DO jk = 1, jpkm1 
    745                  DO jj = 1, jpjm1 
    746                     DO ji = 1, fs_jpim1   ! vector opt. 
    747  
    748                        ze2vr=1.0/e2v(ji,jj) 
    749                        zdyt=zdjt(ji,jj,jk)*ze2vr 
    750                        zdys=zdjs(ji,jj,jk)*ze2vr 
    751  
    752                        ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 
    753                        zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
    754                        zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
    755                        zslope=zsy(ji,jj+jp,jk,1-jp,kp) 
    756  
    757                        zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
    758  
    759                        ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr 
    760                        fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr 
    761                        ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr 
    762                        fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr 
    763                        tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt 
    764                        sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys 
    765                        wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ & 
    766                             & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2 
    767                        psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 
    768  
    769                     END DO 
    770                  END DO 
    771               END DO 
    772            END DO 
    773         END DO 
    774  
    775         tfw(:,:,1)=0.e0 
    776         sfw(:,:,1)=0.e0 
    777         wslp2(:,:,1)=0.e0 
    778  
    779         CALL lbc_lnk( wslp2, 'W', 1. ) 
    780         CALL lbc_lnk( tfw, 'W', 1. ) 
    781         CALL lbc_lnk( sfw, 'W', 1. ) 
    782         CALL lbc_lnk( ftu, 'U', -1. ) 
    783         CALL lbc_lnk( fsu, 'U', -1. ) 
    784         CALL lbc_lnk( ftv, 'V', -1. ) 
    785         CALL lbc_lnk( fsv, 'V', -1. ) 
    786         CALL lbc_lnk( ftud, 'U', -1. ) 
    787         CALL lbc_lnk( fsud, 'U', -1. ) 
    788         CALL lbc_lnk( ftvd, 'V', -1. ) 
    789         CALL lbc_lnk( fsvd, 'V', -1. ) 
    790         CALL lbc_lnk( psix_eiv, 'U', -1. ) 
    791         CALL lbc_lnk( psiy_eiv, 'V', -1. ) 
     414               iku = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1, 2 )    ! last ocean level 
     415               ikv = MAX( MIN( mbathy(ji,jj), mbathy(ji,jj+1  ) ) - 1, 2 ) 
     416               zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem)                          ! i-gradient of T and S 
     417               zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 
     418               zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem)                          ! j-gradient of T and S 
     419               zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 
     420            END DO 
     421         END DO 
     422      ENDIF 
     423      ! 
     424      zdkt(:,:,1) = 0._wp                    !==  before vertical T & S gradient at w-level  ==! 
     425      zdks(:,:,1) = 0._wp 
     426      DO jk = 2, jpk 
     427         zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
     428         zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
     429      END DO 
     430      ! 
     431      ! 
     432      DO jl = 0, 1                           !==  density i-, j-, and k-gradients  ==! 
     433         ip = jl   ;   jp = jl         ! guaranteed nonzero gradients ( absolute value larger than repsln) 
     434         DO jk = 1, jpkm1                          ! done each pair of triad 
     435            DO jj = 1, jpjm1                       ! NB: not masked due to the minimum value set 
     436               DO ji = 1, fs_jpim1   ! vector opt.  
     437                  zdxrho_raw = ( zalpha(ji+ip,jj   ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj   ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 
     438                  zdyrho_raw = ( zalpha(ji   ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji   ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 
     439                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )    ! keep the sign 
     440                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX(   repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     441               END DO 
     442            END DO 
     443         END DO 
     444      END DO 
     445     DO kp = 0, 1                           !==  density i-, j-, and k-gradients  ==! 
     446         DO jk = 1, jpkm1                          ! done each pair of triad 
     447            DO jj = 1, jpj                       ! NB: not masked due to the minimum value set 
     448               DO ji = 1, jpi   ! vector opt.  
     449                   zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   & 
     450                     &       / fse3w(ji,jj,jk+kp) 
     451                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )                    ! force zdzrho >= repsln 
     452               END DO 
     453            END DO 
     454         END DO 
     455      END DO 
     456      ! 
     457      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==! 
     458         DO ji = 1, jpi                            ! i.e. 1 / (hmld+e3t(nmln))  where hmld=depw(nmln) 
     459            jk = MIN( nmln(ji,jj), mbathy(ji,jj) - 1 ) + 1     ! MIN in case ML depth is the ocean depth 
     460            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk) 
     461         END DO 
     462      END DO 
     463      ! 
     464      !                                      !==  intialisations to zero  ==! 
     465      ! 
     466      wslp2  (:,:,:)     = 0._wp                                           ! wslp2 will be cumulated 3D field set to zero 
     467      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero 
     468      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp 
     469!!gm _iso set to zero missing 
     470      triadi (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero 
     471      triadj (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp 
     472       
     473      !-------------------------------------! 
     474      !  Triads just below the Mixed Layer  ! 
     475      !-------------------------------------! 
     476      ! 
     477      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 
     478         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED  
     479            DO jj = 1, jpjm1 
     480               DO ji = 1, fs_jpim1 
     481                  ip = jl   ;   jp = jl 
     482                  jk = MIN( nmln(ji+ip,jj), mbathy(ji+ip,jj) - 1 ) + 1     ! ML level+1 (MIN in case ML depth is the ocean depth) 
     483                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      & 
     484                     &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk) 
     485                  jk = MIN( nmln(ji,jj+jp), mbathy(ji,jj+jp) - 1 ) + 1 
     486                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      & 
     487                     &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
     488                  zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
     489                  ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     490               END DO 
     491            END DO 
     492         END DO 
     493      END DO 
     494 
     495      !-------------------------------------! 
     496      !  Triads with surface limits         ! 
     497      !-------------------------------------! 
     498      ! 
     499      DO kp = 0, 1                        ! k-index of triads 
     500         DO jl = 0, 1 
     501            ip = jl   ;   jp = jl         ! i- and j-indices of triads (i-k and j-k planes) 
     502            DO jk = 1, jpkm1 
     503               DO jj = 1, jpjm1 
     504                  DO ji = 1, fs_jpim1   ! vector opt. 
     505                     ! 
     506                     ! Calculate slope relative to geopotentials used for GM skew fluxes 
     507                     ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth) 
     508                     ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 
     509                     ! masked by umask taken at the level of dz(rho) 
     510                     ! 
     511                     ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 
     512                     ! 
     513                     zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked 
     514                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp) 
     515                     zti_coord = ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
     516                     ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
     517                  ! unmasked 
     518                     zti_g_raw = zti_raw + zti_coord      ! ref to geopot surfaces 
     519                     ztj_g_raw = ztj_raw + ztj_coord 
     520                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
     521                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     522                     ! 
     523                     ! Below  ML use limited zti_g as is 
     524                     ! Inside ML replace by linearly reducing sx_mlb towards surface 
     525                     ! 
     526                     zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp )  ! k index of uppermost point(s) of triad is jk+kp-1 
     527                     zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp )  ! must be .ge. nmln(ji,jj) for zfact=1 
     528                     !                                                          !                   otherwise  zfact=0 
     529                     zti_g_lim =           zfacti   * zti_g_lim                       & 
     530                        &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   & 
     531                        &                           * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 
     532                     ztj_g_lim =           zfactj   * ztj_g_lim                       & 
     533                        &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   & 
     534                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp)                   ! masked 
     535                     ! 
     536                     triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) 
     537                     triadj_g(ji   ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp) 
     538                     ! 
     539                     ! Get coefficients of isoneutral diffusion tensor 
     540                     ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 
     541                     ! 2. We require that isoneutral diffusion  gives no vertical buoyancy flux 
     542                     !     i.e. 33 term = (real slope* 31, 13 terms) 
     543                     ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms 
     544                     ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 
     545                     ! 
     546                     zti_lim  = zti_g_lim - zti_coord    ! remove the coordinate slope  ==> relative to coordinate surfaces 
     547                     ztj_lim  = ztj_g_lim - ztj_coord                                  
     548                     zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp)      ! square of limited slopes            ! masked <<== 
     549                     ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 
     550                     ! 
     551                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   ) 
     552                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   ) 
     553                     zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 
     554                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 
     555                     ! 
     556                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim2 / zti_raw                                          ! masked 
     557                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 
     558                     ! 
     559!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked 
     560                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2             ! masked 
     561                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 
     562                  END DO 
     563               END DO 
     564            END DO 
     565         END DO 
     566      END DO 
     567      ! 
     568      wslp2(:,:,1) = 0._wp                ! force the surface wslp to zero 
     569       
     570      CALL lbc_lnk( wslp2, 'W', 1. )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked 
    792571      ! 
    793572   END SUBROUTINE ldf_slp_grif 
    … …  
    1008787      !!---------------------------------------------------------------------- 
    1009788      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     789      INTEGER ::   ierr         ! local integer 
    1010790      !!---------------------------------------------------------------------- 
    1011791       
    1012792      IF(lwp) THEN     
    1013793         WRITE(numout,*) 
    1014          WRITE(numout,*) 'ldf_slp : direction of lateral mixing' 
    1015          WRITE(numout,*) '~~~~~~~' 
     794         WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 
     795         WRITE(numout,*) '~~~~~~~~~~~~' 
    1016796      ENDIF 
    1017  
    1018       ! Direction of lateral diffusion (tracers and/or momentum) 
    1019       ! ------------------------------ 
    1020       ! set the slope to zero (even in s-coordinates) 
    1021  
    1022       uslp (:,:,:) = 0.e0 
    1023       vslp (:,:,:) = 0.e0 
    1024       wslpi(:,:,:) = 0.e0 
    1025       wslpj(:,:,:) = 0.e0 
    1026  
    1027       uslpml (:,:) = 0.e0 
    1028       vslpml (:,:) = 0.e0 
    1029       wslpiml(:,:) = 0.e0 
    1030       wslpjml(:,:) = 0.e0 
    1031  
    1032       IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 
    1033          IF(lwp) THEN 
    1034             WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
     797       
     798      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes 
     799         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 
     800         ALLOCATE( triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1)                      , STAT=ierr ) 
     801         IF( ierr > 0 ) THEN 
     802            CALL ctl_stop( 'ldf_slp_init : unable to allocate Griffies operator slope ' )   ;   RETURN 
    1035803         ENDIF 
    1036  
    1037          ! geopotential diffusion in s-coordinates on tracers and/or momentum 
    1038          ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 
    1039          ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 
    1040  
    1041          ! set the slope of diffusion to the slope of s-surfaces 
    1042          !      ( c a u t i o n : minus sign as fsdep has positive value ) 
    1043          DO jk = 1, jpk 
    1044             DO jj = 2, jpjm1 
    1045                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1046                   uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 
    1047                   vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
    1048                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
    1049                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
     804         ! 
     805         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 
     806         ! 
     807         IF( ( ln_traldf_hor .AND. ln_dynldf_hor ) .AND. ln_sco )   & 
     808            &     CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator ',   & 
     809            &                                              'in s-coordinate not supported' ) 
     810         ! 
     811      ELSE                             ! Madec operator : slopes at u-, v-, and w-points 
     812         ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                & 
     813            &   omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj)   , vslpml(jpi,jpj)    , wslpiml(jpi,jpj)   , wslpjml(jpi,jpj) , STAT=ierr ) 
     814         IF( ierr > 0 ) THEN 
     815            CALL ctl_stop( 'ldf_slp_init : unable to allocate Madec operator slope ' )   ;   RETURN 
     816         ENDIF 
     817 
     818         ! Direction of lateral diffusion (tracers and/or momentum) 
     819         ! ------------------------------ 
     820         uslp (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp      ! set the slope to zero (even in s-coordinates) 
     821         vslp (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp 
     822         wslpi(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp 
     823         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp 
     824 
     825!!gm I no longer understand this..... 
     826         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 
     827            IF(lwp) THEN 
     828               WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
     829            ENDIF 
     830 
     831            ! geopotential diffusion in s-coordinates on tracers and/or momentum 
     832            ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 
     833            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 
     834 
     835            ! set the slope of diffusion to the slope of s-surfaces 
     836            !      ( c a u t i o n : minus sign as fsdep has positive value ) 
     837            DO jk = 1, jpk 
     838               DO jj = 2, jpjm1 
     839                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     840                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 
     841                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
     842                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
     843                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
     844                  END DO 
    1050845               END DO 
    1051846            END DO 
    1052          END DO 
    1053          ! Lateral boundary conditions on the slopes 
    1054          CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
    1055          CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
    1056       ENDIF 
    1057       ! 
     847            ! Lateral boundary conditions on the slopes 
     848            CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
     849            CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
     850         ENDIF 
     851      ENDIF      ! 
    1058852   END SUBROUTINE ldf_slp_init 
    1059853 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90

    r2287 r2371  
    6767      NAMELIST/namtra_ldf/ ln_traldf_lap  , ln_traldf_bilap,                  & 
    6868         &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   & 
    69          &                 ln_traldf_grif ,                                   & 
     69         &                 ln_traldf_grif , ln_traldf_gdia,                   & 
    7070         &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       & 
    7171         &                 rn_slpmax 
    … …  
    8282         WRITE(numout,*) 'ldf_tra_init : lateral tracer physics' 
    8383         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    84          WRITE(numout,*) '   Namelist namtra_ldf : lateral mixing coefficients' 
     84         WRITE(numout,*) '   Namelist namtra_ldf : lateral mixing parameters (type, direction, coefficients)' 
    8585         WRITE(numout,*) '      laplacian operator            ln_traldf_lap   = ', ln_traldf_lap 
    8686         WRITE(numout,*) '      bilaplacian operator          ln_traldf_bilap = ', ln_traldf_bilap 
    87          WRITE(numout,*) '      griffies    operator          ln_traldf_grif  = ', ln_traldf_grif  
     87         WRITE(numout,*) '      iso-level                     ln_traldf_level = ', ln_traldf_level 
     88         WRITE(numout,*) '      horizontal (geopotential)     ln_traldf_hor   = ', ln_traldf_hor 
     89         WRITE(numout,*) '      iso-neutral                   ln_traldf_iso   = ', ln_traldf_iso 
     90         WRITE(numout,*) '      iso-neutral (Griffies)        ln_traldf_grif  = ', ln_traldf_grif 
     91         WRITE(numout,*) '      Griffies strmfn diagnostics   ln_traldf_gdia  = ', ln_traldf_gdia 
    8892         WRITE(numout,*) '      lateral eddy diffusivity      rn_aht_0        = ', rn_aht_0 
    8993         WRITE(numout,*) '      background hor. diffusivity   rn_ahtb_0       = ', rn_ahtb_0 
    9094         WRITE(numout,*) '      eddy induced velocity coef.   rn_aeiv_0       = ', rn_aeiv_0 
     95         WRITE(numout,*) '      maximum isoppycnal slope      rn_slpmax       = ', rn_slpmax 
     96         WRITE(numout,*) '   + griffies operator internal controls not set via the namelist (experimental): ' 
     97         WRITE(numout,*) '      calculate triads twice        l_triad_iso     = ', l_triad_iso 
     98         WRITE(numout,*) '      no Shapiro filter             l_no_smooth     = ', l_no_smooth 
    9199         WRITE(numout,*) 
    92100      ENDIF 
    93101 
    94       slpmax = rn_slpmax 
    95102      !                                ! convert DOCTOR namelist names into OLD names 
    96103      aht0  = rn_aht_0 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90

    r2287 r2371  
    2121   LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.    !: iso-neutral direction 
    2222   LOGICAL , PUBLIC ::   ln_traldf_grif  = .FALSE.   !: griffies skew flux 
     23   LOGICAL , PUBLIC ::   ln_traldf_gdia  = .FALSE.   !: griffies skew flux streamfunction diagnostics 
    2324   REAL(wp), PUBLIC ::   rn_aht_0        = 2000._wp  !: lateral eddy diffusivity (m2/s) 
    2425   REAL(wp), PUBLIC ::   rn_ahtb_0       =    0._wp  !: lateral background eddy diffusivity (m2/s) 
    … …  
    2728 
    2829   REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names 
    29    REAL(wp), PUBLIC ::   slpmax                     !: slope limit  
     30   LOGICAL , PUBLIC ::   l_triad_iso     = .FALSE.   !: calculate triads twice 
     31   LOGICAL , PUBLIC ::   l_no_smooth     = .FALSE.   !: no Shapiro smoothing 
    3032 
    3133#if defined key_traldf_c3d 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r2287 r2371  
    1717   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function 
    1818   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
     19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp 
    1920   !!---------------------------------------------------------------------- 
    2021 
    … …  
    2627   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    2728   !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
     29   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff. 
    2830   !!   tfreez         : Compute the surface freezing temperature 
    2931   !!   eos_init       : set eos parameters (namelist) 
    … …  
    3840   PRIVATE 
    3941 
    40    !! * Interface  
     42   !                   !! * Interface  
    4143   INTERFACE eos 
    4244      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
    … …  
    4951   PUBLIC   eos_init   ! called by istate module 
    5052   PUBLIC   bn2        ! called by step module 
     53   PUBLIC   eos_alpbet ! called by ldfslp module 
    5154   PUBLIC   tfreez     ! called by sbcice_... modules 
    5255 
    53    !                                        !!* Namelist (nameos) * 
    54    INTEGER , PUBLIC ::   nn_eos   = 0        !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    55    REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4  !: thermal expension coeff. (linear equation of state) 
    56    REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4  !: saline  expension coeff. (linear equation of state) 
    57  
    58    REAL(wp), PUBLIC ::   ralpbet           !: alpha / beta ratio 
     56   !                                          !!* Namelist (nameos) * 
     57   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     58   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 
     59   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state) 
     60 
     61   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
    5962    
    6063   !! * Substitutions 
    … …  
    6467   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    6568   !! $Id$ 
    66    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6770   !!---------------------------------------------------------------------- 
    68  
    6971CONTAINS 
    7072 
    … …  
    135137                  ! 
    136138                  ! compute volumic mass pure water at atm pressure 
    137                   zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt   & 
    138                      &      -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     139                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     140                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    139141                  ! seawater volumic mass atm pressure 
    140                   zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt        & 
    141                      &                   -4.0899e-3 ) *zt+0.824493 
    142                   zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    143                   zr4= 4.8314e-4 
     142                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     143                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     144                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     145                  zr4= 4.8314e-4_wp 
    144146                  ! 
    145147                  ! potential volumic mass (reference to the surface) 
    … …  
    147149                  ! 
    148150                  ! add the compression terms 
    149                   ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    150                   zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     151                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     152                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    151153                  zb = zbw + ze * zs 
    152154                  ! 
    153                   zd = -2.042967e-2 
    154                   zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    155                   zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 
     155                  zd = -2.042967e-2_wp 
     156                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     157                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    156158                  za = ( zd*zsr + zc ) *zs + zaw 
    157159                  ! 
    158                   zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    159                   za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    160                   zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 
     160                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     161                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     162                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    161163                  zk0= ( zb1*zsr + za1 )*zs + zkw 
    162164                  ! 
    163165                  ! masked in situ density anomaly 
    164                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     166                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    165167                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
    166168               END DO 
    … …  
    170172      CASE( 1 )                !==  Linear formulation function of temperature only  ==! 
    171173         DO jk = 1, jpkm1 
    172             prd(:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
     174            prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    173175         END DO 
    174176         ! 
    … …  
    258260                  ! 
    259261                  ! compute volumic mass pure water at atm pressure 
    260                   zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   & 
    261                      &                       -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     262                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     263                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    262264                  ! seawater volumic mass atm pressure 
    263                   zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   & 
    264                      &                                   -4.0899e-3 ) *zt+0.824493 
    265                   zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    266                   zr4= 4.8314e-4 
     265                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
     266                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
     267                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     268                  zr4= 4.8314e-4_wp 
    267269                  ! 
    268270                  ! potential volumic mass (reference to the surface) 
    … …  
    273275                  ! 
    274276                  ! add the compression terms 
    275                   ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    276                   zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     277                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     278                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    277279                  zb = zbw + ze * zs 
    278280                  ! 
    279                   zd = -2.042967e-2 
    280                   zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    281                   zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 
     281                  zd = -2.042967e-2_wp 
     282                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     283                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    282284                  za = ( zd*zsr + zc ) *zs + zaw 
    283285                  ! 
    284                   zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    285                   za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    286                   zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 
     286                  zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp 
     287                  za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp 
     288                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    287289                  zk0= ( zb1*zsr + za1 )*zs + zkw 
    288290                  ! 
    289291                  ! masked in situ density anomaly 
    290                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     292                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    291293                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
    292294               END DO 
    … …  
    296298      CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    297299         DO jk = 1, jpkm1 
    298             prd  (:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk) 
    299             prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
     300            prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk) 
     301            prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
    300302         END DO 
    301303         ! 
    … …  
    303305         DO jk = 1, jpkm1 
    304306            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    305             prhop(:,:,jk) = ( 1.e0     + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
     307            prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
    306308         END DO 
    307309         ! 
    … …  
    382384               ! 
    383385               ! compute volumic mass pure water at atm pressure 
    384                zr1 = ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   & 
    385                   &                        -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     386               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     387                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    386388               ! seawater volumic mass atm pressure 
    387                zr2 = ( ( ( 5.3875e-9*zt-8.2467e-7 )*zt+7.6438e-5 ) *zt   & 
    388                   &                                   -4.0899e-3 ) *zt+0.824493 
    389                zr3 = ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    390                zr4 = 4.8314e-4 
     389               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
     390                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
     391               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
     392               zr4 = 4.8314e-4_wp 
    391393               ! 
    392394               ! potential volumic mass (reference to the surface) 
    … …  
    394396               ! 
    395397               ! add the compression terms 
    396                ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    397                zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     398               ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     399               zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    398400               zb = zbw + ze * zs 
    399401               ! 
    400                zd = -2.042967e-2 
    401                zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    402                zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt -4.721788 
     402               zd =    -2.042967e-2_wp 
     403               zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     404               zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 
    403405               za = ( zd*zsr + zc ) *zs + zaw 
    404406               ! 
    405                zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    406                za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    407                zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt   & 
    408                   &                                       +2098.925 ) *zt+190925.6 
     407               zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp 
     408               za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp 
     409               zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   & 
     410                  &                             +2098.925_wp ) *zt+190925.6_wp 
    409411               zk0= ( zb1*zsr + za1 )*zs + zkw 
    410412               ! 
    411413               ! masked in situ density anomaly 
    412                prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     414               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
    413415            END DO 
    414416         END DO 
    … …  
    417419         DO jj = 1, jpjm1 
    418420            DO ji = 1, fs_jpim1   ! vector opt. 
    419                prd(ji,jj) = ( 0.0285 - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
     421               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    420422            END DO 
    421423         END DO 
    … …  
    490492                  zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point 
    491493                  ! 
    492                   zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta 
    493                      &                               - 0.203814e-03 ) * zt   & 
    494                      &                               + 0.170907e-01 ) * zt   & 
    495                      &   + 0.665157e-01                                      & 
    496                      &   +     ( - 0.678662e-05 * zs                         & 
    497                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    498                      &   +   ( ( - 0.302285e-13 * zh                         & 
    499                      &           - 0.251520e-11 * zs                         & 
    500                      &           + 0.512857e-12 * zt * zt           ) * zh   & 
    501                      &           - 0.164759e-06 * zs                         & 
    502                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    503                      &                               + 0.380374e-04 ) * zh 
     494                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta 
     495                     &                                  - 0.203814e-03_wp ) * zt   & 
     496                     &                                  + 0.170907e-01_wp ) * zt   & 
     497                     &   +         0.665157e-01_wp                                 & 
     498                     &   +     ( - 0.678662e-05_wp * zs                            & 
     499                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     500                     &   +   ( ( - 0.302285e-13_wp * zh                            & 
     501                     &           - 0.251520e-11_wp * zs                            & 
     502                     &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     503                     &           - 0.164759e-06_wp * zs                            & 
     504                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     505                     &                                  + 0.380374e-04_wp ) * zh 
    504506                     ! 
    505                   zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta 
    506                      &                            - 0.301985e-05 ) * zt      & 
    507                      &   + 0.785567e-03                                      & 
    508                      &   + (     0.515032e-08 * zs                           & 
    509                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    510                      &   +(  (   0.121551e-17 * zh                           & 
    511                      &         - 0.602281e-15 * zs                           & 
    512                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    513                      &                             + 0.408195e-10   * zs     & 
    514                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    515                      &                             - 0.121555e-07 ) * zh 
     507                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta 
     508                     &                               - 0.301985e-05_wp ) * zt      & 
     509                     &   +       0.785567e-03_wp                                   & 
     510                     &   + (     0.515032e-08_wp * zs                              & 
     511                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     & 
     512                     &   + ( (   0.121551e-17_wp * zh                              & 
     513                     &         - 0.602281e-15_wp * zs                              & 
     514                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     & 
     515                     &                                + 0.408195e-10_wp   * zs     & 
     516                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     & 
     517                     &                                - 0.121555e-07_wp ) * zh 
    516518                     ! 
    517519                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2  
    … …  
    521523                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    522524                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    523                   IF ( ABS( zds) <= 1.e-20 ) zds = 1.e-20 
     525                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    524526                  rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    525527#endif 
    … …  
    544546               DO ji = 1, jpi 
    545547                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )   
    546                   IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20 
     548                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    547549                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    548550               END DO 
    … …  
    560562 
    561563 
     564   SUBROUTINE eos_alpbet( pts, palph, pbeta ) 
     565      !!---------------------------------------------------------------------- 
     566      !!                 ***  ROUTINE ldf_slp_grif  *** 
     567      !! 
     568      !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points.  
     569      !! 
     570      !! ** Method  :   calculates alpha and beta at T-points  
     571      !!       * nn_eos = 0  : UNESCO sea water properties 
     572      !!         The brunt-vaisala frequency is computed using the polynomial 
     573      !!      polynomial expression of McDougall (1987): 
     574      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
     575      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
     576      !!      computed and used in zdfddm module : 
     577      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
     578      !!       * nn_eos = 1  : linear equation of state (temperature only) 
     579      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
     580      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
     581      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
     582      !!       * nn_eos = 3  : Jackett JAOT 2003 ??? 
     583      !! 
     584      !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point 
     585      !!---------------------------------------------------------------------- 
     586      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity 
     587      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff. 
     588      !! 
     589      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     590      REAL(wp) ::   zt, zs, zh   ! local scalars  
     591      !!---------------------------------------------------------------------- 
     592      ! 
     593      SELECT CASE ( nn_eos ) 
     594      ! 
     595      CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
     596         DO jk = 1, jpk 
     597            DO jj = 1, jpj 
     598               DO ji = 1, jpi 
     599                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
     600                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
     601                  zh = fsdept(ji,jj,jk)              ! depth in meters  
     602                  ! 
     603                  pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   & 
     604                     &                                        - 0.301985e-05_wp ) * zt   & 
     605                     &           + 0.785567e-03_wp                                       & 
     606                     &           + (     0.515032e-08_wp * zs                            & 
     607                     &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   & 
     608                     &           + ( (   0.121551e-17_wp * zh                            & 
     609                     &                 - 0.602281e-15_wp * zs                            & 
     610                     &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   & 
     611                     &                                        + 0.408195e-10_wp   * zs   & 
     612                     &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   & 
     613                     &                                        - 0.121555e-07_wp ) * zh 
     614                     ! 
     615                  palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             & 
     616                      &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
     617                      &                                  - 0.203814e-03_wp ) * zt   & 
     618                      &                                  + 0.170907e-01_wp ) * zt   & 
     619                      &   + 0.665157e-01_wp                                         & 
     620                      &   +     ( - 0.678662e-05_wp * zs                            & 
     621                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     622                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
     623                      &           - 0.251520e-11_wp * zs                            & 
     624                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     625                      &           - 0.164759e-06_wp * zs                            & 
     626                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     627                      &                                  + 0.380374e-04_wp ) * zh) 
     628               END DO 
     629            END DO 
     630         END DO 
     631         ! 
     632      CASE ( 1 ) 
     633         palph(:,:,:) = - rn_alpha 
     634         pbeta(:,:,:) =   0._wp 
     635         ! 
     636      CASE ( 2 ) 
     637         palph(:,:,:) = - rn_alpha 
     638         pbeta(:,:,:) =   rn_beta 
     639         ! 
     640      CASE DEFAULT 
     641         IF(lwp) WRITE(numout,cform_err) 
     642         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     643         nstop = nstop + 1 
     644         ! 
     645      END SELECT 
     646      ! 
     647   END SUBROUTINE eos_alpbet 
     648 
     649 
    562650   FUNCTION tfreez( psal ) RESULT( ptf ) 
    563651      !!---------------------------------------------------------------------- 
    … …  
    576664      !!---------------------------------------------------------------------- 
    577665      ! 
    578       ptf(:,:) = ( - 0.0575 + 1.710523e-3 * SQRT( psal(:,:) )   & 
    579          &                  - 2.154996e-4 *       psal(:,:)   ) * psal(:,:) 
     666      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
     667         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
    580668      ! 
    581669   END FUNCTION tfreez 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r2287 r2371  
    7373      CASE ( 1 )     
    7474         IF ( ln_traldf_grif ) THEN 
    75             CALL tra_ldf_iso_grif    ( kt )           ! Griffies quarter-cell formulation 
     75 
     76            CALL tra_ldf_iso_grif    ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0  )           ! Griffies quarter-cell formulation 
    7677         ELSE 
    7778            CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian 
    … …  
    8586         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    8687         IF ( ln_traldf_grif ) THEN 
    87             CALL tra_ldf_iso_grif    ( kt ) 
     88            CALL tra_ldf_iso_grif    ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 ) 
    8889         ELSE 
    8990            CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )   
    … …  
    133134      INTEGER ::   ioptio, ierr         ! temporary integers  
    134135!       
    135 !     NAMELIST/namtra_ldf/ ln_traldf_lap  , ln_traldf_bilap,                  & 
    136 !        &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   & 
    137 !        &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0 
    138136      !!---------------------------------------------------------------------- 
    139137 
    … …  
    141139      ! =============================================== 
    142140     
    143 !     REWIND( numnam )                ! Namelist namtra_ldf already read in ldftra module 
    144 !     READ  ( numnam, namtra_ldf )     
    145  
    146141      IF(lwp) THEN                    ! Namelist print 
    147142         WRITE(numout,*) 
    148143         WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator' 
    149144         WRITE(numout,*) '~~~~~~~~~~~' 
    150          WRITE(numout,*) '   Namelist namtra_ldf : set lateral mixing parameters (type, direction, coefficients)' 
    151          WRITE(numout,*) '      laplacian operator          ln_traldf_lap   = ', ln_traldf_lap 
    152          WRITE(numout,*) '      bilaplacian operator        ln_traldf_bilap = ', ln_traldf_bilap 
    153          WRITE(numout,*) '      iso-level                   ln_traldf_level = ', ln_traldf_level 
    154          WRITE(numout,*) '      horizontal (geopotential)   ln_traldf_hor   = ', ln_traldf_hor 
    155          WRITE(numout,*) '      iso-neutral                 ln_traldf_iso   = ', ln_traldf_iso 
    156          WRITE(numout,*) '      iso-neutral (Griffies)      ln_traldf_grif  = ', ln_traldf_grif 
     145         WRITE(numout,*) '   Namelist namtra_ldf already read in ldftra module' 
     146         WRITE(numout,*) '   see ldf_tra_init report for lateral mixing parameters' 
     147         WRITE(numout,*) 
    157148      ENDIF 
    158149 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r2287 r2371  
    11MODULE traldf_iso_grif 
    2    !!============================================================================== 
     2   !!====================================================================== 
    33   !!                   ***  MODULE  traldf_iso_grif  *** 
    4    !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend 
    5    !!============================================================================== 
    6    !! History :   9.0  !  06-10  (C. Harris) 
     4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend 
     5   !!====================================================================== 
     6   !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)   
     7   !!                !          Griffies operator version adapted from traldf_iso.F90 
    78   !!---------------------------------------------------------------------- 
    89#if   defined key_ldfslp   ||   defined key_esopa 
    … …  
    1415   !!                       and with the vertical part of 
    1516   !!                       the isopycnal or geopotential s-coord. operator  
    16    !!                       vector optimization, use k-j-i loops. 
    17    !!---------------------------------------------------------------------- 
    18  
     17   !!---------------------------------------------------------------------- 
    1918   USE oce             ! ocean dynamics and active tracers 
    2019   USE dom_oce         ! ocean space and time domain 
    2120   USE ldftra_oce      ! ocean active tracers: lateral physics 
    22    USE trdmod          ! ocean active tracers trends  
    23    USE trdmod_oce      ! ocean variables trends 
    2421   USE zdf_oce         ! ocean vertical physics 
    2522   USE in_out_manager  ! I/O manager 
     23   USE iom             ! 
    2624   USE ldfslp          ! iso-neutral slopes 
    2725   USE diaptr          ! poleward transport diagnostics 
    28    USE prtctl          ! Print control 
     26   USE trc_oce         ! share passive tracers/Ocean variables 
     27#if defined key_diaar5 
     28   USE phycst          ! physical constants 
     29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     30#endif 
    2931 
    3032   IMPLICIT NONE 
    … …  
    3234 
    3335   PUBLIC tra_ldf_iso_grif   ! routine called by traldf.F90 
     36 
     37   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   psix_eiv 
     38   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   psiy_eiv 
     39   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   ah_wslp2 
    3440 
    3541   !! * Substitutions 
    3642#  include "domzgr_substitute.h90" 
    3743#  include "ldftra_substitute.h90" 
     44#  include "vectopt_loop_substitute.h90" 
    3845#  include "ldfeiv_substitute.h90" 
    39 #  include "vectopt_loop_substitute.h90" 
    40  
    4146   !!---------------------------------------------------------------------- 
    4247   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    4752CONTAINS 
    4853 
    49    SUBROUTINE tra_ldf_iso_grif( kt ) 
     54  SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv,  & 
     55       &                                ptb, pta, kjpt, pahtb0 ) 
     56    !!---------------------------------------------------------------------- 
     57    !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
     58    !! 
     59    !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
     60    !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and  
     61    !!      add it to the general trend of tracer equation. 
     62    !! 
     63    !! ** Method  :   The horizontal component of the lateral diffusive trends  
     64    !!      is provided by a 2nd order operator rotated along neural or geopo- 
     65    !!      tential surfaces to which an eddy induced advection can be added 
     66    !!      It is computed using before fields (forward in time) and isopyc- 
     67    !!      nal or geopotential slopes computed in routine ldfslp. 
     68    !! 
     69    !!      1st part :  masked horizontal derivative of T  ( di[ t ] ) 
     70    !!      ========    with partial cell update if ln_zps=T. 
     71    !! 
     72    !!      2nd part :  horizontal fluxes of the lateral mixing operator 
     73    !!      ========     
     74    !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
     75    !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ] 
     76    !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 
     77    !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ] 
     78    !!      take the horizontal divergence of the fluxes: 
     79    !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
     80    !!      Add this trend to the general trend (ta,sa): 
     81    !!         ta = ta + difft 
     82    !! 
     83    !!      3rd part: vertical trends of the lateral mixing operator 
     84    !!      ========  (excluding the vertical flux proportional to dk[t] ) 
     85    !!      vertical fluxes associated with the rotated lateral mixing: 
     86    !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ] 
     87    !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  } 
     88    !!      take the horizontal divergence of the fluxes: 
     89    !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
     90    !!      Add this trend to the general trend (ta,sa): 
     91    !!         pta = pta + difft 
     92    !! 
     93    !! ** Action :   Update pta arrays with the before rotated diffusion 
     94    !!---------------------------------------------------------------------- 
     95    USE oce,   zftu => ua   ! use ua as workspace 
     96    USE oce,   zftv => va   ! use va as workspace 
     97    !! 
     98    INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
     99    CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     100    INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     101    REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     102    REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
     103    REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     104    REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
     105    !! 
     106    INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices 
     107    INTEGER  ::  ip,jp,kp        ! dummy loop indices 
     108    INTEGER  ::  iku, ikv        ! temporary integer 
     109    INTEGER  ::  ierr            ! temporary integer 
     110    REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
     111    REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
     112    REAL(wp) ::  zcoef0, zbtr                  !   -      - 
     113    REAL(wp), DIMENSION(jpi,jpj,0:1) ::   zdkt               ! 2D+1 workspace 
     114    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw   ! 3D workspace 
     115      ! 
     116    REAL(wp) ::   zslope_skew,zslope_iso,zslope2, zbu, zbv 
     117    REAL(wp) ::   ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt 
     118    REAL(wp) ::   ah,zah_slp,zaei_slp 
     119#if defined key_diaar5 
     120    REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                ! 2D workspace 
     121    REAL(wp)                         ::   zztmp              ! local scalar 
     122#endif 
    50123      !!---------------------------------------------------------------------- 
    51       !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
    52       !! 
    53       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    54       !!      trend for a laplacian tensor (except the dz[ dz[.] ] term) and  
    55       !!      add it to the general trend of tracer equation. 
    56       !! 
    57       !! ** Method  :   The horizontal component of the lateral diffusive trends  
    58       !!      is provided by a 2nd order operator rotated along neutral or geopo- 
    59       !!      tential surfaces to which an eddy induced advection can be added 
    60       !!      It is computed using before fields (forward in time) 
    61       !! 
    62       !!      1st part :  masked horizontal derivative of T & S ( di[ t ] ) 
    63       !!      ========    with partial cell update if ln_zps=T. 
    64       !! 
    65       !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    66       !!      ========     
    67       !!      take the horizontal divergence of the fluxes: 
    68       !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
    69       !!      Add this trend to the general trend (ta,sa): 
    70       !!         ta = ta + difft 
    71       !! 
    72       !!      3rd part: vertical trends of the lateral mixing operator 
    73       !!      ========  (excluding the vertical flux proportional to dk[t] ) 
    74       !!      vertical fluxes associated with the rotated lateral mixing: 
    75  
    76       !!      take the horizontal divergence of the fluxes: 
    77       !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
    78       !!      Add this trend to the general trend (ta,sa): 
    79       !!         ta = ta + difft 
    80       !! 
    81       !! ** Action : 
    82       !!         Update (ta,sa) arrays with the before rotated diffusion trend 
    83       !!            (except the dk[ dk[.] ] term) 
    84       !! 
     124 
     125    IF( kt == nit000 )  THEN 
     126       IF(lwp) WRITE(numout,*) 
     127       IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 
     128       IF(lwp) WRITE(numout,*) '                   WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 
     129       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     130       ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 
     131       IF( ierr > 0 ) THEN 
     132          CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' )   ;   RETURN 
     133       ENDIF 
     134       IF( ln_traldf_gdia ) THEN 
     135          ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 
     136          IF( ierr > 0 ) THEN 
     137             CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' )   ;   RETURN 
     138          ENDIF 
     139       ENDIF 
     140    ENDIF 
     141 
     142      !                                                
    85143      !!---------------------------------------------------------------------- 
    86       USE oce           , zftv => ua   ! use ua as workspace 
    87       USE oce           , zfsv => va   ! use va as workspace 
    88       !! 
    89       INTEGER, INTENT( in ) ::   kt        ! ocean time-step index 
    90       !! 
    91       INTEGER  ::   ji, jj, jk, ip, jp, kp    ! dummy loop indices 
    92       INTEGER  ::   iku, ikv                  ! temporary integer 
    93       REAL(wp) ::   zatempu, zdx, zta         ! temporary scalars 
    94       REAL(wp) ::   zatempv, zdy, zsa         !    "         " 
    95       REAL(wp) ::   zslopec, zdsloper, zepsln !    "         " 
    96       REAL(wp) ::   zsxe, za_sxe, zfact       !    "         " 
    97       REAL(wp) ::   zbtr                      !    "         " 
    98       REAL(wp), DIMENSION(2) ::   zdelta   ! 1D workspace 
    99       REAL(wp), DIMENSION(jpi,jpj)     ::   zftu   ! 2D workspace 
    100       REAL(wp), DIMENSION(jpi,jpj)     ::   zfsu   !    "           " 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zdkt     ! 3D workspace 
    102       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdis, zdjs, zdks     !  "      " 
    103  
    104  
    105  
     144      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv  
    106145      !!---------------------------------------------------------------------- 
    107146 
    108       IF( kt == nit000 ) THEN 
    109          IF(lwp) WRITE(numout,*) 
    110          IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator' 
    111          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    112       ENDIF 
    113  
    114       IF ( .NOT. lk_traldf_eiv ) THEN 
    115         fsaeiu(:,:,:)=0.0 
    116         fsaeiv(:,:,:)=0.0 
    117         fsaeiw(:,:,:)=0.0 
    118       ENDIF 
    119  
    120       DO jk=1,jpkm1 
    121       DO jj=1,jpjm1 
    122       DO ji=1,fs_jpim1 
    123                ftu(ji,jj,jk)=ftud(ji,jj,jk)+ftu(ji,jj,jk)*(fsahtu(ji,jj,jk)-fsaeiu(ji,jj,jk)) 
    124                fsu(ji,jj,jk)=fsud(ji,jj,jk)+fsu(ji,jj,jk)*(fsahtu(ji,jj,jk)-fsaeiu(ji,jj,jk)) 
    125                ftv(ji,jj,jk)=ftvd(ji,jj,jk)+ftv(ji,jj,jk)*(fsahtv(ji,jj,jk)-fsaeiv(ji,jj,jk)) 
    126                fsv(ji,jj,jk)=fsvd(ji,jj,jk)+fsv(ji,jj,jk)*(fsahtv(ji,jj,jk)-fsaeiv(ji,jj,jk)) 
     147!!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 
     148 
     149    ah_wslp2(:,:,:) = 0.0 
     150    IF( ln_traldf_gdia ) THEN 
     151       psix_eiv(:,:,:) = 0.0 
     152       psiy_eiv(:,:,:) = 0.0 
     153    ENDIF 
     154 
     155    DO ip=0,1 
     156       DO kp=0,1 
     157          DO jk=1,jpkm1 
     158             DO jj=1,jpjm1 
     159                DO ji=1,fs_jpim1 
     160                   ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 
     161                   zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
     162                   ah = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk) 
     163                   zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     164                   zslope2=(zslope_skew & 
     165                        & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur & 
     166                        & )**2 
     167                   ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 
     168                        & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 
     169                   IF( ln_traldf_gdia ) THEN 
     170                      zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 
     171                      psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
     172                   ENDIF 
     173                END DO 
     174             END DO 
     175          END DO 
     176       END DO 
     177    END DO 
     178 
     179    DO jp=0,1 
     180       DO kp=0,1 
     181          DO jk=1,jpkm1 
     182             DO jj=1,jpjm1 
     183                DO ji=1,fs_jpim1 
     184                   ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 
     185                   zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
     186                   ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 
     187                   zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     188                   zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 
     189                        & )**2 
     190                   ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 
     191                        & ah*((zbv*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*zslope2 
     192                   IF( ln_traldf_gdia ) THEN 
     193                      zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew     !fsaeit(ji,jj+jp,jk)*zslope_skew 
     194                      psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
     195                   ENDIF 
     196                END DO 
     197             END DO 
     198          END DO 
     199       END DO 
    127200      END DO 
    128       END DO 
    129       END DO 
    130  
    131       DO jk = 1, jpkm1 
    132  
    133          ! II.4 Second derivative (divergence) and add to the general trend 
    134          ! ---------------------------------------------------------------- 
    135          DO jj = 2 , jpjm1 
    136             DO ji = fs_2, fs_jpim1   ! vector opt. 
    137                zbtr= 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
    138                zta = zbtr * ( ftu(ji,jj,jk) - ftu(ji-1,jj,jk) + ftv(ji,jj,jk) - ftv(ji,jj-1,jk)  ) 
    139                zsa = zbtr * ( fsu(ji,jj,jk) - fsu(ji-1,jj,jk) + fsv(ji,jj,jk) - fsv(ji,jj-1,jk)  ) 
    140                ta (ji,jj,jk) = ta (ji,jj,jk) + zta 
    141                sa (ji,jj,jk) = sa (ji,jj,jk) + zsa 
     201 
     202      ! 
     203      !                                                          ! =========== 
     204      DO jn = 1, kjpt                                            ! tracer loop 
     205         !                                                       ! =========== 
     206         ! Zero fluxes for each tracer 
     207         ztfw(:,:,:) = 0._wp 
     208         zftu(:,:,:) = 0._wp 
     209         zftv(:,:,:) = 0._wp 
     210         !                                                
     211         DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==! 
     212            DO jj = 1, jpjm1 
     213               DO ji = 1, fs_jpim1   ! vector opt. 
     214                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     215                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     216               END DO 
    142217            END DO 
    143218         END DO 
    144          !                                          ! =============== 
    145       END DO                                        !   End of slab   
    146       !                                             ! =============== 
    147  
    148       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN   ! Poleward diffusive heat and salt transports 
    149          pht_ldf(:) = ptr_vj( zftv(:,:,:) ) 
    150          pst_ldf(:) = ptr_vj( zfsv(:,:,:) ) 
    151       ENDIF 
    152  
    153       !!---------------------------------------------------------------------- 
    154       !!   III - vertical trend of T & S (extra diagonal terms only) 
    155       !!---------------------------------------------------------------------- 
    156  
    157       ! I.5 Divergence of vertical fluxes added to the general tracer trend 
    158       ! ------------------------------------------------------------------- 
    159  
    160       DO jk=1,jpk 
    161       DO jj=2,jpjm1 
    162       DO ji=fs_2,fs_jpim1 
    163                tfw(ji,jj,jk)=tfw(ji,jj,jk)*(fsahtw(ji,jj,jk)+fsaeiw(ji,jj,jk)) 
    164                sfw(ji,jj,jk)=sfw(ji,jj,jk)*(fsahtw(ji,jj,jk)+fsaeiw(ji,jj,jk)) 
    165       END DO 
    166       END DO 
    167       END DO 
    168  
    169       DO jk = 1, jpkm1 
    170          DO jj = 2, jpjm1 
    171             DO ji = fs_2, fs_jpim1   ! vector opt. 
    172                zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
    173                zta  = (  tfw(ji,jj,jk) - tfw(ji,jj,jk+1)  ) * zbtr 
    174                zsa  = (  sfw(ji,jj,jk) - sfw(ji,jj,jk+1)  ) * zbtr 
    175                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    176                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
     219         IF( ln_zps ) THEN                               ! partial steps: correction at the last level 
     220# if defined key_vectopt_loop 
     221            DO jj = 1, 1 
     222               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     223# else 
     224            DO jj = 1, jpjm1 
     225               DO ji = 1, jpim1 
     226# endif 
     227                  iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1                ! last level 
     228                  ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
     229                  zdit(ji,jj,iku) = pgu(ji,jj,jn)           
     230                  zdjt(ji,jj,ikv) = pgv(ji,jj,jn)       
     231               END DO 
    177232            END DO 
    178          END DO 
    179       END DO 
    180       ! 
    181       DO jk=1,jpk 
    182       DO jj=2,jpjm1 
    183       DO ji=fs_2,fs_jpim1 
    184        psix_eiv(ji,jj,jk) = psix_eiv(ji,jj,jk)*fsaeiu(ji,jj,jk) 
    185        psiy_eiv(ji,jj,jk) = psiy_eiv(ji,jj,jk)*fsaeiv(ji,jj,jk) 
    186       END DO 
    187       END DO 
    188       END DO 
    189  
    190    END SUBROUTINE tra_ldf_iso_grif 
     233         ENDIF 
     234 
     235         !!---------------------------------------------------------------------- 
     236         !!   II - horizontal trend  (full) 
     237         !!---------------------------------------------------------------------- 
     238         ! 
     239         DO jk = 1, jpkm1 
     240            ! 
     241            !                    !==  Vertical tracer gradient at level jk and jk+1 
     242            zdkt(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
     243            ! 
     244            !                          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     245            IF( jk == 1 ) THEN   ;   zdkt(:,:,0) = zdkt(:,:,1) 
     246            ELSE                 ;   zdkt(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
     247            ENDIF 
     248 
     249            IF( .NOT. l_triad_iso ) THEN 
     250               triadi = triadi_g 
     251               triadj = triadj_g 
     252            ENDIF 
     253 
     254            DO ip=0,1              !==  Horizontal & vertical fluxes 
     255               DO kp=0,1 
     256                  DO jj=1,jpjm1 
     257                     DO ji=1,fs_jpim1 
     258                        ze1ur = 1._wp / e1u(ji,jj) 
     259                        zdxt = zdit(ji,jj,jk) * ze1ur 
     260                        ze3wr = 1._wp/fse3w(ji+ip,jj,jk+kp) 
     261                        zdzt  = zdkt(ji+ip,jj,kp) * ze3wr 
     262                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     263                        zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
     264 
     265                        zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     266                        ah = fsahtu(ji,jj,jk)!*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)             ===>>  ???? 
     267                        zah_slp = ah*zslope_iso 
     268                        zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew 
     269                        zftu(ji,jj,jk)  = zftu(ji,jj,jk) - (ah*zdxt + (zah_slp - zaei_slp)*zdzt)*zbu*ze1ur 
     270                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp)*zdxt *zbu*ze3wr 
     271                     END DO 
     272                  END DO 
     273               END DO 
     274            END DO 
     275 
     276            DO jp=0,1 
     277               DO kp=0,1 
     278                  DO jj=1,jpjm1 
     279                     DO ji=1,fs_jpim1 
     280                        ze2vr = 1._wp/e2v(ji,jj) 
     281                        zdyt = zdjt(ji,jj,jk)*ze2vr 
     282                        ze3wr = 1._wp/fse3w(ji,jj+jp,jk+kp) 
     283                        zdzt = zdkt(ji,jj+jp,kp) * ze3wr 
     284                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     285                        zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 
     286                        zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     287                        ah    = fsahtv(ji,jj,jk)!*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk) 
     288                        zah_slp = ah * zslope_iso 
     289                        zaei_slp = fsaeiw(ji,jj+jp,jk)*zslope_skew!fsaeit(ji,jj+jp,jk)*zslope_skew 
     290                        zftv(ji,jj,jk) = zftv(ji,jj,jk) - (ah*zdyt + (zah_slp - zaei_slp)*zdzt)*zbv*ze2vr 
     291                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp)*zdyt*zbv*ze3wr 
     292                     END DO 
     293                  END DO 
     294               END DO 
     295            END DO 
     296 
     297          !                        !==  divergence and add to the general trend  ==! 
     298          DO jj = 2 , jpjm1 
     299             DO ji = fs_2, fs_jpim1   ! vector opt. 
     300                zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     301                pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (   zftu(ji-1,jj,jk) - zftu(ji,jj,jk)   & 
     302                  &                                            + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   ) 
     303             END DO 
     304          END DO 
     305          ! 
     306       END DO 
     307       ! 
     308       DO jk = 1, jpkm1            !== Divergence of vertical fluxes added to the general tracer trend 
     309          DO jj = 2, jpjm1 
     310             DO ji = fs_2, fs_jpim1   ! vector opt. 
     311                pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     312                  &                                 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     313             END DO 
     314          END DO 
     315       END DO 
     316       ! 
     317       !                            ! "Poleward" diffusive heat or salt transports (T-S case only) 
     318       IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     319         IF( jn == jp_tem)   pht_ldf(:) = ptr_vj( zftv(:,:,:) )        ! 3.3  names 
     320         IF( jn == jp_sal)   pst_ldf(:) = ptr_vj( zftv(:,:,:) ) 
     321       ENDIF 
     322 
     323#if defined key_diaar5 
     324       IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     325          zztmp = 0.5 * rau0 * rcp  
     326          z2d(:,:) = 0.e0  
     327          DO jk = 1, jpkm1 
     328             DO jj = 2, jpjm1 
     329                DO ji = fs_2, fs_jpim1   ! vector opt. 
     330                   z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk)       & 
     331                        &       * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk)  
     332                END DO 
     333             END DO 
     334          END DO 
     335          CALL lbc_lnk( z2d, 'U', -1. ) 
     336          CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     337          z2d(:,:) = 0.e0  
     338          DO jk = 1, jpkm1 
     339             DO jj = 2, jpjm1 
     340                DO ji = fs_2, fs_jpim1   ! vector opt. 
     341                   z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk)       & 
     342                        &       * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk)  
     343                END DO 
     344             END DO 
     345          END DO 
     346          CALL lbc_lnk( z2d, 'V', -1. ) 
     347          CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     348       END IF 
     349#endif 
     350       ! 
     351    END DO 
     352    ! 
     353  END SUBROUTINE tra_ldf_iso_grif 
    191354 
    192355#else 
    … …  
    196359CONTAINS 
    197360   SUBROUTINE tra_ldf_iso_grif( kt )               ! Empty routine 
    198       WRITE(*,*) 'tra_ldf_iso_grif: You should not have seen this print! error?', kt 
     361      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt 
    199362   END SUBROUTINE tra_ldf_iso_grif 
    200363#endif 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r2287 r2371  
    9191      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace 
    9292      USE oce    , ONLY :   zws   => va   ! va  -          - 
     93      USE traldf_iso_grif    , ONLY :   ah_wslp2 
    9394      !!  
    9495      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    … …  
    145146                  DO jj = 2, jpjm1 
    146147                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    147                         zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
     148                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
     149                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing 
    148150                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
    149151                     END DO 
    … …  
    201203                  DO jj = 2, jpjm1 
    202204                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    203                         zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
     205                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing 
     206                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
    204207                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
    205208                     END DO 
Note: See TracChangeset for help on using the changeset viewer.