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 4933 for branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2014-12-01T11:11:43+01:00 (10 years ago)
Author:
cetlod
Message:

dev_CNRS_CICE : merging CNRS and CICE branche

Location:
branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r4624 r4933  
    1515   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d 
    1616   !!             -   ! 2003-08  (G. Madec)  F90, free form 
    17    !!            3.0  ! 2006-08  (G. Madec)  add tfreez function 
     17   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function (now eos_fzp 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 
     19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add alpha/beta used in ldfslp 
     20   !!            3.7  ! 2012-03  (F. Roquet, G. Madec)  add primitive of alpha and beta used in PE computation 
     21   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state 
     22   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module 
     23   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80 
    2024   !!---------------------------------------------------------------------- 
    2125 
    … …  
    2327   !!   eos            : generic interface of the equation of state 
    2428   !!   eos_insitu     : Compute the in situ density 
    25    !!   eos_insitu_pot : Compute the insitu and surface referenced potential 
    26    !!                    volumic mass 
     29   !!   eos_insitu_pot : Compute the insitu and surface referenced potential volumic mass 
    2730   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    28    !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
    29    !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio 
    30    !!   tfreez         : Compute the surface freezing temperature 
     31   !!   bn2            : Compute the Brunt-Vaisala frequency 
     32   !!   eos_rab        : generic interface of in situ thermal/haline expansion ratio  
     33   !!   eos_rab_3d     : compute in situ thermal/haline expansion ratio 
     34   !!   eos_rab_2d     : compute in situ thermal/haline expansion ratio for 2d fields 
     35   !!   eos_fzp        : freezing temperature 
    3136   !!   eos_init       : set eos parameters (namelist) 
    3237   !!---------------------------------------------------------------------- 
    3338   USE dom_oce         ! ocean space and time domain 
    3439   USE phycst          ! physical constants 
    35    USE zdfddm          ! vertical physics: double diffusion 
     40   ! 
    3641   USE in_out_manager  ! I/O manager 
    3742   USE lib_mpp         ! MPP library 
     43   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3844   USE prtctl          ! Print control 
    3945   USE wrk_nemo        ! Memory Allocation 
     46   USE lbclnk         ! ocean lateral boundary conditions 
    4047   USE timing          ! Timing 
    4148 
    … …  
    4754      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
    4855   END INTERFACE 
    49    INTERFACE bn2 
    50       MODULE PROCEDURE eos_bn2 
     56   ! 
     57   INTERFACE eos_rab 
     58      MODULE PROCEDURE rab_3d, rab_2d 
    5159   END INTERFACE 
    52  
    53    PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules 
    54    PUBLIC   eos_init   ! called by istate module 
    55    PUBLIC   bn2        ! called by step module 
    56    PUBLIC   eos_alpbet ! called by ldfslp module 
    57    PUBLIC   tfreez     ! called by sbcice_... modules 
    58  
    59    !                                  !!* Namelist (nameos) * 
    60    INTEGER , PUBLIC ::   nn_eos       !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    61    REAL(wp), PUBLIC ::   rn_alpha     !: thermal expension coeff. (linear equation of state) 
    62    REAL(wp), PUBLIC ::   rn_beta      !: saline  expension coeff. (linear equation of state) 
    63  
    64    REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
     60   ! 
     61   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules 
     62   PUBLIC   bn2            ! called by step module 
     63   PUBLIC   eos_rab        ! called by ldfslp, zdfddm, trabbl 
     64   PUBLIC   eos_pt_from_ct ! called by sbcssm 
     65   PUBLIC   eos_fzp        ! called by traadv_cen2 and sbcice_... modules 
     66   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen module 
     67   PUBLIC   eos_init       ! called by istate module 
     68 
     69   !                                          !!* Namelist (nameos) * 
     70   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     71   LOGICAL , PUBLIC ::   ln_useCT  = .FALSE.  ! determine if eos_pt_from_ct is used to compute sst_m 
     72 
     73   !                                   !!!  simplified eos coefficients 
     74   ! default value: Vallis 2006 
     75   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.  
     76   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.  
     77   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2         
     78   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2         
     79   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T   
     80   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S   
     81   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt   
     82    
     83   ! TEOS10/EOS80 parameters 
     84   REAL(wp) ::   r1_S0, r1_T0, r1_Z0, rdeltaS 
     85    
     86   ! EOS parameters 
     87   REAL(wp) ::   EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 
     88   REAL(wp) ::   EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510 
     89   REAL(wp) ::   EOS020 , EOS120 , EOS220 , EOS320 , EOS420 
     90   REAL(wp) ::   EOS030 , EOS130 , EOS230 , EOS330 
     91   REAL(wp) ::   EOS040 , EOS140 , EOS240 
     92   REAL(wp) ::   EOS050 , EOS150 
     93   REAL(wp) ::   EOS060 
     94   REAL(wp) ::   EOS001 , EOS101 , EOS201 , EOS301 , EOS401 
     95   REAL(wp) ::   EOS011 , EOS111 , EOS211 , EOS311 
     96   REAL(wp) ::   EOS021 , EOS121 , EOS221 
     97   REAL(wp) ::   EOS031 , EOS131 
     98   REAL(wp) ::   EOS041 
     99   REAL(wp) ::   EOS002 , EOS102 , EOS202 
     100   REAL(wp) ::   EOS012 , EOS112 
     101   REAL(wp) ::   EOS022 
     102   REAL(wp) ::   EOS003 , EOS103 
     103   REAL(wp) ::   EOS013  
     104    
     105   ! ALPHA parameters 
     106   REAL(wp) ::   ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 
     107   REAL(wp) ::   ALP010 , ALP110 , ALP210 , ALP310 , ALP410 
     108   REAL(wp) ::   ALP020 , ALP120 , ALP220 , ALP320 
     109   REAL(wp) ::   ALP030 , ALP130 , ALP230 
     110   REAL(wp) ::   ALP040 , ALP140 
     111   REAL(wp) ::   ALP050 
     112   REAL(wp) ::   ALP001 , ALP101 , ALP201 , ALP301 
     113   REAL(wp) ::   ALP011 , ALP111 , ALP211 
     114   REAL(wp) ::   ALP021 , ALP121 
     115   REAL(wp) ::   ALP031 
     116   REAL(wp) ::   ALP002 , ALP102 
     117   REAL(wp) ::   ALP012 
     118   REAL(wp) ::   ALP003 
     119    
     120   ! BETA parameters 
     121   REAL(wp) ::   BET000 , BET100 , BET200 , BET300 , BET400 , BET500 
     122   REAL(wp) ::   BET010 , BET110 , BET210 , BET310 , BET410 
     123   REAL(wp) ::   BET020 , BET120 , BET220 , BET320 
     124   REAL(wp) ::   BET030 , BET130 , BET230 
     125   REAL(wp) ::   BET040 , BET140 
     126   REAL(wp) ::   BET050 
     127   REAL(wp) ::   BET001 , BET101 , BET201 , BET301 
     128   REAL(wp) ::   BET011 , BET111 , BET211 
     129   REAL(wp) ::   BET021 , BET121 
     130   REAL(wp) ::   BET031 
     131   REAL(wp) ::   BET002 , BET102 
     132   REAL(wp) ::   BET012 
     133   REAL(wp) ::   BET003 
     134 
     135   ! PEN parameters 
     136   REAL(wp) ::   PEN000 , PEN100 , PEN200 , PEN300 , PEN400 
     137   REAL(wp) ::   PEN010 , PEN110 , PEN210 , PEN310 
     138   REAL(wp) ::   PEN020 , PEN120 , PEN220 
     139   REAL(wp) ::   PEN030 , PEN130 
     140   REAL(wp) ::   PEN040 
     141   REAL(wp) ::   PEN001 , PEN101 , PEN201 
     142   REAL(wp) ::   PEN011 , PEN111 
     143   REAL(wp) ::   PEN021 
     144   REAL(wp) ::   PEN002 , PEN102 
     145   REAL(wp) ::   PEN012 
     146    
     147   ! ALPHA_PEN parameters 
     148   REAL(wp) ::   APE000 , APE100 , APE200 , APE300 
     149   REAL(wp) ::   APE010 , APE110 , APE210 
     150   REAL(wp) ::   APE020 , APE120 
     151   REAL(wp) ::   APE030 
     152   REAL(wp) ::   APE001 , APE101 
     153   REAL(wp) ::   APE011 
     154   REAL(wp) ::   APE002 
     155 
     156   ! BETA_PEN parameters 
     157   REAL(wp) ::   BPE000 , BPE100 , BPE200 , BPE300 
     158   REAL(wp) ::   BPE010 , BPE110 , BPE210 
     159   REAL(wp) ::   BPE020 , BPE120 
     160   REAL(wp) ::   BPE030 
     161   REAL(wp) ::   BPE001 , BPE101 
     162   REAL(wp) ::   BPE011 
     163   REAL(wp) ::   BPE002 
    65164 
    66165   !! * Substitutions 
    … …  
    68167#  include "vectopt_loop_substitute.h90" 
    69168   !!---------------------------------------------------------------------- 
    70    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     169   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    71170   !! $Id$ 
    72171   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    82181      !!       defined through the namelist parameter nn_eos. 
    83182      !! 
    84       !! ** Method  :   3 cases: 
    85       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
    86       !!         the in situ density is computed directly as a function of 
    87       !!         potential temperature relative to the surface (the opa t 
    88       !!         variable), salt and pressure (assuming no pressure variation 
    89       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    90       !!         is approximated by the depth in meters. 
    91       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    92       !!         with pressure                      p        decibars 
    93       !!              potential temperature         t        deg celsius 
    94       !!              salinity                      s        psu 
    95       !!              reference volumic mass        rau0     kg/m**3 
    96       !!              in situ volumic mass          rho      kg/m**3 
    97       !!              in situ density anomalie      prd      no units 
    98       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
    99       !!          t = 40 deg celcius, s=40 psu 
    100       !!      nn_eos = 1 : linear equation of state function of temperature only 
    101       !!              prd(t) = 0.0285 - rn_alpha * t 
    102       !!      nn_eos = 2 : linear equation of state function of temperature and 
    103       !!               salinity 
    104       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
    105       !!      Note that no boundary condition problem occurs in this routine 
    106       !!      as pts are defined over the whole domain. 
     183      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
     184      !!         with   prd    in situ density anomaly      no units 
     185      !!                t      TEOS10: CT or EOS80: PT      Celsius 
     186      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu 
     187      !!                z      depth                        meters 
     188      !!                rho    in situ density              kg/m^3 
     189      !!                rau0   reference density            kg/m^3 
     190      !! 
     191      !!     nn_eos = -1 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     192      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celcius, sa=35.5 g/kg 
     193      !! 
     194      !!     nn_eos =  0 : polynomial EOS-80 equation of state is used for rho(t,s,z). 
     195      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celcius, sp=35.5 psu 
     196      !! 
     197      !!     nn_eos =  1 : simplified equation of state 
     198      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0 
     199      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
     200      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     201      !!              Vallis like equation: use default values of coefficients 
    107202      !! 
    108203      !! ** Action  :   compute prd , the in situ density (no units) 
    109204      !! 
    110       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
    111       !!---------------------------------------------------------------------- 
    112       !! 
    113       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    114       !                                                      ! 2 : salinity               [psu] 
    115       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-] 
    116       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m] 
    117       !! 
    118       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    119       REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    120       REAL(wp) ::   zr1, zr2, zr3, zr4   !   -      - 
    121       REAL(wp) ::   zrhop, ze, zbw, zb   !   -      - 
    122       REAL(wp) ::   zd , zc , zaw, za    !   -      - 
    123       REAL(wp) ::   zb1, za1, zkw, zk0   !   -      - 
    124       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    125       !!---------------------------------------------------------------------- 
    126  
    127       ! 
    128       IF( nn_timing == 1 ) CALL timing_start('eos') 
    129       ! 
    130       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     205      !! References :   Roquet et al, Ocean Modelling, in preparation (2014) 
     206      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     207      !!                TEOS-10 Manual, 2010 
     208      !!---------------------------------------------------------------------- 
     209      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
     210      !                                                               ! 2 : salinity               [psu] 
     211      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd   ! in situ density            [-] 
     212      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep  ! depth                      [m] 
     213      ! 
     214      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     215      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
     216      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     217      !!---------------------------------------------------------------------- 
     218      ! 
     219      IF( nn_timing == 1 )   CALL timing_start('eos-insitu') 
    131220      ! 
    132221      SELECT CASE( nn_eos ) 
    133222      ! 
    134       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    135 !CDIR NOVERRCHK 
    136          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     223      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    137224         ! 
    138225         DO jk = 1, jpkm1 
    139226            DO jj = 1, jpj 
    140227               DO ji = 1, jpi 
    141                   zt = pts   (ji,jj,jk,jp_tem) 
    142                   zs = pts   (ji,jj,jk,jp_sal) 
    143                   zh = pdep(ji,jj,jk)        ! depth 
    144                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    145                   ! 
    146                   ! compute volumic mass pure water at atm pressure 
    147                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    148                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    149                   ! seawater volumic mass atm pressure 
    150                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    151                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    152                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    153                   zr4= 4.8314e-4_wp 
    154                   ! 
    155                   ! potential volumic mass (reference to the surface) 
    156                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    157                   ! 
    158                   ! add the compression terms 
    159                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    160                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    161                   zb = zbw + ze * zs 
    162                   ! 
    163                   zd = -2.042967e-2_wp 
    164                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    165                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    166                   za = ( zd*zsr + zc ) *zs + zaw 
    167                   ! 
    168                   zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    169                   za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    170                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    171                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    172                   ! 
    173                   ! masked in situ density anomaly 
    174                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    175                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     228                  ! 
     229                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     230                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     231                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     232                  ztm = tmask(ji,jj,jk)                                         ! tmask 
     233                  ! 
     234                  zn3 = EOS013*zt   & 
     235                     &   + EOS103*zs+EOS003 
     236                     ! 
     237                  zn2 = (EOS022*zt   & 
     238                     &   + EOS112*zs+EOS012)*zt   & 
     239                     &   + (EOS202*zs+EOS102)*zs+EOS002 
     240                     ! 
     241                  zn1 = (((EOS041*zt   & 
     242                     &   + EOS131*zs+EOS031)*zt   & 
     243                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     244                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     245                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     246                     ! 
     247                  zn0 = (((((EOS060*zt   & 
     248                     &   + EOS150*zs+EOS050)*zt   & 
     249                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     250                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     251                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     252                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     253                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     254                     ! 
     255                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     256                  ! 
     257                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     258                  ! 
    176259               END DO 
    177260            END DO 
    178261         END DO 
    179262         ! 
    180       CASE( 1 )                !==  Linear formulation function of temperature only  ==! 
     263      CASE( 1 )                !==  simplified EOS  ==! 
     264         ! 
    181265         DO jk = 1, jpkm1 
    182             prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
     266            DO jj = 1, jpj 
     267               DO ji = 1, jpi 
     268                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     269                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     270                  zh  = pdep (ji,jj,jk) 
     271                  ztm = tmask(ji,jj,jk) 
     272                  ! 
     273                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     274                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     275                     &  - rn_nu * zt * zs 
     276                     !                                  
     277                  prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked) 
     278               END DO 
     279            END DO 
    183280         END DO 
    184281         ! 
    185       CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==! 
    186          DO jk = 1, jpkm1 
    187             prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    188          END DO 
    189          ! 
    190282      END SELECT 
    191283      ! 
    192       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk ) 
    193       ! 
    194       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    195       ! 
    196       IF( nn_timing == 1 ) CALL timing_stop('eos') 
     284      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', ovlap=1, kdim=jpk ) 
     285      ! 
     286      IF( nn_timing == 1 )   CALL timing_stop('eos-insitu') 
    197287      ! 
    198288   END SUBROUTINE eos_insitu 
    … …  
    208298      !!     namelist parameter nn_eos. 
    209299      !! 
    210       !! ** Method  : 
    211       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
    212       !!         the in situ density is computed directly as a function of 
    213       !!         potential temperature relative to the surface (the opa t 
    214       !!         variable), salt and pressure (assuming no pressure variation 
    215       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    216       !!         is approximated by the depth in meters. 
    217       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    218       !!              rhop(t,s)  = rho(t,s,0) 
    219       !!         with pressure                      p        decibars 
    220       !!              potential temperature         t        deg celsius 
    221       !!              salinity                      s        psu 
    222       !!              reference volumic mass        rau0     kg/m**3 
    223       !!              in situ volumic mass          rho      kg/m**3 
    224       !!              in situ density anomalie      prd      no units 
    225       !! 
    226       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
    227       !!          t = 40 deg celcius, s=40 psu 
    228       !! 
    229       !!      nn_eos = 1 : linear equation of state function of temperature only 
    230       !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t 
    231       !!              rhop(t,s)  = rho(t,s) 
    232       !! 
    233       !!      nn_eos = 2 : linear equation of state function of temperature and 
    234       !!               salinity 
    235       !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
    236       !!                       = rn_beta * s - rn_alpha * tn - 1. 
    237       !!              rhop(t,s)  = rho(t,s) 
    238       !!      Note that no boundary condition problem occurs in this routine 
    239       !!      as (tn,sn) or (ta,sa) are defined over the whole domain. 
    240       !! 
    241300      !! ** Action  : - prd  , the in situ density (no units) 
    242301      !!              - prhop, the potential volumic mass (Kg/m3) 
    243302      !! 
    244       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
    245       !!                Brown and Campana, Mon. Weather Rev., 1978 
    246       !!---------------------------------------------------------------------- 
    247       !! 
     303      !!---------------------------------------------------------------------- 
    248304      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius] 
    249305      !                                                                ! 2 : salinity               [psu] 
    … …  
    252308      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    253309      ! 
    254       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    255       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars 
    256       REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0               !   -      - 
    257       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    258       !!---------------------------------------------------------------------- 
    259       ! 
    260       IF( nn_timing == 1 ) CALL timing_start('eos-p') 
    261       ! 
    262       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     310      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     311      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
     312      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     313      !!---------------------------------------------------------------------- 
     314      ! 
     315      IF( nn_timing == 1 )   CALL timing_start('eos-pot') 
    263316      ! 
    264317      SELECT CASE ( nn_eos ) 
    265318      ! 
    266       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    267 !CDIR NOVERRCHK 
    268          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     319      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    269320         ! 
    270321         DO jk = 1, jpkm1 
    271322            DO jj = 1, jpj 
    272323               DO ji = 1, jpi 
    273                   zt = pts   (ji,jj,jk,jp_tem) 
    274                   zs = pts   (ji,jj,jk,jp_sal) 
    275                   zh = pdep(ji,jj,jk)        ! depth 
    276                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    277                   ! 
    278                   ! compute volumic mass pure water at atm pressure 
    279                   zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    280                      &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    281                   ! seawater volumic mass atm pressure 
    282                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
    283                      &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
    284                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    285                   zr4= 4.8314e-4_wp 
    286                   ! 
    287                   ! potential volumic mass (reference to the surface) 
    288                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    289                   ! 
    290                   ! save potential volumic mass 
    291                   prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 
    292                   ! 
    293                   ! add the compression terms 
    294                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    295                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    296                   zb = zbw + ze * zs 
    297                   ! 
    298                   zd = -2.042967e-2_wp 
    299                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    300                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    301                   za = ( zd*zsr + zc ) *zs + zaw 
    302                   ! 
    303                   zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp 
    304                   za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp 
    305                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    306                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    307                   ! 
    308                   ! masked in situ density anomaly 
    309                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    310                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     324                  ! 
     325                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     326                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     327                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     328                  ztm = tmask(ji,jj,jk)                                         ! tmask 
     329                  ! 
     330                  zn3 = EOS013*zt   & 
     331                     &   + EOS103*zs+EOS003 
     332                     ! 
     333                  zn2 = (EOS022*zt   & 
     334                     &   + EOS112*zs+EOS012)*zt   & 
     335                     &   + (EOS202*zs+EOS102)*zs+EOS002 
     336                     ! 
     337                  zn1 = (((EOS041*zt   & 
     338                     &   + EOS131*zs+EOS031)*zt   & 
     339                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     340                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     341                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     342                     ! 
     343                  zn0 = (((((EOS060*zt   & 
     344                     &   + EOS150*zs+EOS050)*zt   & 
     345                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     346                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     347                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     348                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     349                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     350                     ! 
     351                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     352                  ! 
     353                  prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     354                  ! 
     355                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
    311356               END DO 
    312357            END DO 
    313358         END DO 
    314359         ! 
    315       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
     360      CASE( 1 )                !==  simplified EOS  ==! 
     361         ! 
    316362         DO jk = 1, jpkm1 
    317             prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    318             prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
     363            DO jj = 1, jpj 
     364               DO ji = 1, jpi 
     365                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     366                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     367                  zh  = pdep (ji,jj,jk) 
     368                  ztm = tmask(ji,jj,jk) 
     369                  !                                                     ! potential density referenced at the surface 
     370                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
     371                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
     372                     &  - rn_nu * zt * zs 
     373                  prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 
     374                  !                                                     ! density anomaly (masked) 
     375                  zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
     376                  prd(ji,jj,jk) = zn * r1_rau0 * ztm 
     377                  ! 
     378               END DO 
     379            END DO 
    319380         END DO 
    320381         ! 
    321       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    322          DO jk = 1, jpkm1 
    323             prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    324             prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
    325          END DO 
    326          ! 
    327382      END SELECT 
    328383      ! 
    329       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 
    330       ! 
    331       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    332       ! 
    333       IF( nn_timing == 1 ) CALL timing_stop('eos-p') 
     384      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 
     385      ! 
     386      IF( nn_timing == 1 )   CALL timing_stop('eos-pot') 
    334387      ! 
    335388   END SUBROUTINE eos_insitu_pot 
    … …  
    344397      !!      defined through the namelist parameter nn_eos. * 2D field case 
    345398      !! 
    346       !! ** Method : 
    347       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
    348       !!         the in situ density is computed directly as a function of 
    349       !!         potential temperature relative to the surface (the opa t 
    350       !!         variable), salt and pressure (assuming no pressure variation 
    351       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    352       !!         is approximated by the depth in meters. 
    353       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    354       !!         with pressure                      p        decibars 
    355       !!              potential temperature         t        deg celsius 
    356       !!              salinity                      s        psu 
    357       !!              reference volumic mass        rau0     kg/m**3 
    358       !!              in situ volumic mass          rho      kg/m**3 
    359       !!              in situ density anomalie      prd      no units 
    360       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
    361       !!          t = 40 deg celcius, s=40 psu 
    362       !!      nn_eos = 1 : linear equation of state function of temperature only 
    363       !!              prd(t) = 0.0285 - rn_alpha * t 
    364       !!      nn_eos = 2 : linear equation of state function of temperature and 
    365       !!               salinity 
    366       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
    367       !!      Note that no boundary condition problem occurs in this routine 
    368       !!      as pts are defined over the whole domain. 
    369       !! 
    370       !! ** Action  : - prd , the in situ density (no units) 
    371       !! 
    372       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
    373       !!---------------------------------------------------------------------- 
    374       !! 
     399      !! ** Action  : - prd , the in situ density (no units) (unmasked) 
     400      !! 
     401      !!---------------------------------------------------------------------- 
    375402      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    376403      !                                                           ! 2 : salinity               [psu] 
    377       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m] 
     404      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
    378405      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    379       !! 
    380       INTEGER  ::   ji, jj                    ! dummy loop indices 
    381       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary scalars 
    382       REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask        !    -         - 
    383       REAL(wp), POINTER, DIMENSION(:,:) :: zws 
    384       !!---------------------------------------------------------------------- 
    385       ! 
    386       IF( nn_timing == 1 ) CALL timing_start('eos2d') 
    387       ! 
    388       CALL wrk_alloc( jpi, jpj, zws ) 
    389       ! 
    390  
     406      ! 
     407      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     408      REAL(wp) ::   zt , zh , zs              ! local scalars 
     409      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     410      !!---------------------------------------------------------------------- 
     411      ! 
     412      IF( nn_timing == 1 )   CALL timing_start('eos2d') 
     413      ! 
    391414      prd(:,:) = 0._wp 
    392  
     415      ! 
    393416      SELECT CASE( nn_eos ) 
    394417      ! 
    395       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    396       ! 
    397 !CDIR NOVERRCHK 
     418      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     419         ! 
    398420         DO jj = 1, jpjm1 
    399 !CDIR NOVERRCHK 
    400421            DO ji = 1, fs_jpim1   ! vector opt. 
    401                zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 
     422               ! 
     423               zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     424               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     425               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     426               ! 
     427               zn3 = EOS013*zt   & 
     428                  &   + EOS103*zs+EOS003 
     429                  ! 
     430               zn2 = (EOS022*zt   & 
     431                  &   + EOS112*zs+EOS012)*zt   & 
     432                  &   + (EOS202*zs+EOS102)*zs+EOS002 
     433                  ! 
     434               zn1 = (((EOS041*zt   & 
     435                  &   + EOS131*zs+EOS031)*zt   & 
     436                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     437                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     438                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     439                  ! 
     440               zn0 = (((((EOS060*zt   & 
     441                  &   + EOS150*zs+EOS050)*zt   & 
     442                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     443                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     444                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     445                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     446                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     447                  ! 
     448               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     449               ! 
     450               prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly 
     451               ! 
    402452            END DO 
    403453         END DO 
     454         ! 
     455         CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions 
     456         ! 
     457      CASE( 1 )                !==  simplified EOS  ==! 
     458         ! 
    404459         DO jj = 1, jpjm1 
    405460            DO ji = 1, fs_jpim1   ! vector opt. 
    406                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
    407                zt    = pts  (ji,jj,jp_tem)            ! interpolated T 
    408                zs    = pts  (ji,jj,jp_sal)            ! interpolated S 
    409                zsr   = zws  (ji,jj)            ! square root of interpolated S 
    410                zh    = pdep (ji,jj)            ! depth at the partial step level 
    411                ! 
    412                ! compute volumic mass pure water at atm pressure 
    413                zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    414                   &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    415                ! seawater volumic mass atm pressure 
    416                zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
    417                   &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
    418                zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
    419                zr4 = 4.8314e-4_wp 
    420                ! 
    421                ! potential volumic mass (reference to the surface) 
    422                zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    423                ! 
    424                ! add the compression terms 
    425                ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    426                zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    427                zb = zbw + ze * zs 
    428                ! 
    429                zd =    -2.042967e-2_wp 
    430                zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    431                zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 
    432                za = ( zd*zsr + zc ) *zs + zaw 
    433                ! 
    434                zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp 
    435                za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp 
    436                zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   & 
    437                   &                             +2098.925_wp ) *zt+190925.6_wp 
    438                zk0= ( zb1*zsr + za1 )*zs + zkw 
    439                ! 
    440                ! masked in situ density anomaly 
    441                prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     461               ! 
     462               zt    = pts  (ji,jj,jp_tem)  - 10._wp 
     463               zs    = pts  (ji,jj,jp_sal)  - 35._wp 
     464               zh    = pdep (ji,jj)                         ! depth at the partial step level 
     465               ! 
     466               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     467                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     468                  &  - rn_nu * zt * zs 
     469                  ! 
     470               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly 
     471               ! 
    442472            END DO 
    443473         END DO 
    444474         ! 
    445       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    446          DO jj = 1, jpjm1 
    447             DO ji = 1, fs_jpim1   ! vector opt. 
    448                prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    449             END DO 
    450          END DO 
    451          ! 
    452       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    453          DO jj = 1, jpjm1 
    454             DO ji = 1, fs_jpim1   ! vector opt. 
    455                prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    456             END DO 
    457          END DO 
     475         CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions 
    458476         ! 
    459477      END SELECT 
    460  
     478      ! 
    461479      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
    462480      ! 
    463       CALL wrk_dealloc( jpi, jpj, zws ) 
    464       ! 
    465       IF( nn_timing == 1 ) CALL timing_stop('eos2d') 
     481      IF( nn_timing == 1 )   CALL timing_stop('eos2d') 
    466482      ! 
    467483   END SUBROUTINE eos_insitu_2d 
    468484 
    469485 
    470    SUBROUTINE eos_bn2( pts, pn2 ) 
    471       !!---------------------------------------------------------------------- 
    472       !!                  ***  ROUTINE eos_bn2  *** 
    473       !! 
    474       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time- 
    475       !!      step of the input arguments 
    476       !! 
    477       !! ** Method : 
    478       !!       * nn_eos = 0  : UNESCO sea water properties 
    479       !!         The brunt-vaisala frequency is computed using the polynomial 
    480       !!      polynomial expression of McDougall (1987): 
    481       !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
    482       !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    483       !!      computed and used in zdfddm module : 
    484       !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
    485       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    486       !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
    487       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    488       !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    489       !!      The use of potential density to compute N^2 introduces e r r o r 
    490       !!      in the sign of N^2 at great depths. We recommand the use of 
    491       !!      nn_eos = 0, except for academical studies. 
    492       !!        Macro-tasked on horizontal slab (jk-loop) 
    493       !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
    494       !!      and is never used at this level. 
    495       !! 
    496       !! ** Action  : - pn2 : the brunt-vaisala frequency 
    497       !! 
    498       !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 
    499       !!---------------------------------------------------------------------- 
    500       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    501       !                                                               ! 2 : salinity               [psu] 
    502       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1] 
    503       !! 
    504       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    505       REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars 
    506 #if defined key_zdfddm 
    507       REAL(wp) ::   zds   ! local scalars 
    508 #endif 
    509       !!---------------------------------------------------------------------- 
    510  
    511       ! 
    512       IF( nn_timing == 1 ) CALL timing_start('bn2') 
    513       ! 
    514       ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    515       ! -------------------------- 
    516       ! 
    517       SELECT CASE( nn_eos ) 
    518       ! 
    519       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    520          DO jk = 2, jpkm1 
     486   SUBROUTINE rab_3d( pts, pab ) 
     487      !!---------------------------------------------------------------------- 
     488      !!                 ***  ROUTINE rab_3d  *** 
     489      !! 
     490      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points 
     491      !! 
     492      !! ** Method  :   calculates alpha / beta at T-points 
     493      !! 
     494      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
     495      !!---------------------------------------------------------------------- 
     496      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
     497      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
     498      ! 
     499      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     500      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
     501      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     502      !!---------------------------------------------------------------------- 
     503      ! 
     504      IF( nn_timing == 1 )   CALL timing_start('rab_3d') 
     505      ! 
     506      SELECT CASE ( nn_eos ) 
     507      ! 
     508      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     509         ! 
     510         DO jk = 1, jpkm1 
    521511            DO jj = 1, jpj 
    522512               DO ji = 1, jpi 
    523                   zgde3w = grav / fse3w(ji,jj,jk) 
    524                   zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt 
    525                   zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0  ! salinity anomaly (s-35) at w-pt 
    526                   zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    527                   ! 
    528                   zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta 
    529                      &                                  - 0.203814e-03_wp ) * zt   & 
    530                      &                                  + 0.170907e-01_wp ) * zt   & 
    531                      &   +         0.665157e-01_wp                                 & 
    532                      &   +     ( - 0.678662e-05_wp * zs                            & 
    533                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    534                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
    535                      &           - 0.251520e-11_wp * zs                            & 
    536                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    537                      &           - 0.164759e-06_wp * zs                            & 
    538                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    539                      &                                  + 0.380374e-04_wp ) * zh 
    540                      ! 
    541                   zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta 
    542                      &                               - 0.301985e-05_wp ) * zt      & 
    543                      &   +       0.785567e-03_wp                                   & 
    544                      &   + (     0.515032e-08_wp * zs                              & 
    545                      &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     & 
    546                      &   + ( (   0.121551e-17_wp * zh                              & 
    547                      &         - 0.602281e-15_wp * zs                              & 
    548                      &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     & 
    549                      &                                + 0.408195e-10_wp   * zs     & 
    550                      &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     & 
    551                      &                                - 0.121555e-07_wp ) * zh 
    552                      ! 
    553                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    554                      &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    555                      &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    556 #if defined key_zdfddm 
    557                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    558                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    559                   IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    560                   rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    561 #endif 
     513                  ! 
     514                  zh  = fsdept(ji,jj,jk) * r1_Z0                                ! depth 
     515                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     516                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     517                  ztm = tmask(ji,jj,jk)                                         ! tmask 
     518                  ! 
     519                  ! alpha 
     520                  zn3 = ALP003 
     521                  ! 
     522                  zn2 = ALP012*zt + ALP102*zs+ALP002 
     523                  ! 
     524                  zn1 = ((ALP031*zt   & 
     525                     &   + ALP121*zs+ALP021)*zt   & 
     526                     &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     527                     &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     528                     ! 
     529                  zn0 = ((((ALP050*zt   & 
     530                     &   + ALP140*zs+ALP040)*zt   & 
     531                     &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     532                     &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     533                     &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     534                     &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     535                     ! 
     536                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     537                  ! 
     538                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 
     539                  ! 
     540                  ! beta 
     541                  zn3 = BET003 
     542                  ! 
     543                  zn2 = BET012*zt + BET102*zs+BET002 
     544                  ! 
     545                  zn1 = ((BET031*zt   & 
     546                     &   + BET121*zs+BET021)*zt   & 
     547                     &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     548                     &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     549                     ! 
     550                  zn0 = ((((BET050*zt   & 
     551                     &   + BET140*zs+BET040)*zt   & 
     552                     &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     553                     &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     554                     &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     555                     &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     556                     ! 
     557                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     558                  ! 
     559                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 
     560                  ! 
    562561               END DO 
    563562            END DO 
    564563         END DO 
    565564         ! 
    566       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    567          DO jk = 2, jpkm1 
    568             pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
    569          END DO 
    570          ! 
    571       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    572          DO jk = 2, jpkm1 
    573             pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      & 
    574                &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   & 
    575                &               / fse3w(:,:,jk) * tmask(:,:,jk) 
    576          END DO 
    577 #if defined key_zdfddm 
    578          DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
     565      CASE( 1 )                  !==  simplified EOS  ==! 
     566         ! 
     567         DO jk = 1, jpkm1 
    579568            DO jj = 1, jpj 
    580569               DO ji = 1, jpi 
    581                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
    582                   IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    583                   rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     570                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     571                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     572                  zh  = fsdept(ji,jj,jk)                 ! depth in meters at t-point 
     573                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
     574                  ! 
     575                  zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     576                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha 
     577                  ! 
     578                  zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     579                  pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta 
     580                  ! 
    584581               END DO 
    585582            END DO 
    586583         END DO 
    587 #endif 
    588       END SELECT 
    589  
    590       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
    591 #if defined key_zdfddm 
    592       IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
    593 #endif 
    594       ! 
    595       IF( nn_timing == 1 ) CALL timing_stop('bn2') 
    596       ! 
    597    END SUBROUTINE eos_bn2 
    598  
    599  
    600    SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 
    601       !!---------------------------------------------------------------------- 
    602       !!                 ***  ROUTINE eos_alpbet  *** 
    603       !! 
    604       !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points 
    605       !! 
    606       !! ** Method  :   calculates alpha / beta ratio at T-points 
    607       !!       * nn_eos = 0  : UNESCO sea water properties 
    608       !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 
    609       !!                       polynomial expression of McDougall (1987). 
    610       !!                       Scalar beta0 is returned = 1. 
    611       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    612       !!                       The ratio is undefined, so we return alpha as palpbet 
    613       !!                       Scalar beta0 is returned = 0. 
    614       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    615       !!                       The alpha/beta ratio is returned as ralpbet 
    616       !!                       Scalar beta0 is returned = 1. 
    617       !! 
    618       !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points 
    619       !!            :   beta0   : 1. or 0. 
    620       !!---------------------------------------------------------------------- 
    621       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
    622       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio 
    623       REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T) 
    624       !! 
    625       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    626       REAL(wp) ::   zt, zs, zh   ! local scalars 
    627       !!---------------------------------------------------------------------- 
    628       ! 
    629       IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 
    630       ! 
    631       SELECT CASE ( nn_eos ) 
    632       ! 
    633       CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
    634          DO jk = 1, jpk 
    635             DO jj = 1, jpj 
    636                DO ji = 1, jpi 
    637                   zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
    638                   zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    639                   zh = fsdept(ji,jj,jk)               ! depth in meters 
    640                   ! 
    641                   palpbet(ji,jj,jk) =                                              & 
    642                      &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
    643                      &                                  - 0.203814e-03_wp ) * zt   & 
    644                      &                                  + 0.170907e-01_wp ) * zt   & 
    645                      &   + 0.665157e-01_wp                                         & 
    646                      &   +     ( - 0.678662e-05_wp * zs                            & 
    647                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    648                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
    649                      &           - 0.251520e-11_wp * zs                            & 
    650                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    651                      &           - 0.164759e-06_wp * zs                            & 
    652                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    653                      &                                  + 0.380374e-04_wp ) * zh 
    654                END DO 
    655             END DO 
    656          END DO 
    657          beta0 = 1._wp 
    658          ! 
    659       CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
    660          palpbet(:,:,:) = rn_alpha 
    661          beta0 = 0._wp 
    662          ! 
    663       CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
    664          palpbet(:,:,:) = ralpbet 
    665          beta0 = 1._wp 
    666584         ! 
    667585      CASE DEFAULT 
    … …  
    672590      END SELECT 
    673591      ! 
    674       IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 
    675       ! 
    676    END SUBROUTINE eos_alpbet 
    677  
    678  
    679    FUNCTION tfreez( psal, pdep ) RESULT( ptf ) 
     592      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
     593         &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk ) 
     594      ! 
     595      IF( nn_timing == 1 )   CALL timing_stop('rab_3d') 
     596      ! 
     597   END SUBROUTINE rab_3d 
     598 
     599 
     600   SUBROUTINE rab_2d( pts, pdep, pab ) 
     601      !!---------------------------------------------------------------------- 
     602      !!                 ***  ROUTINE rab_2d  *** 
     603      !! 
     604      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked) 
     605      !! 
     606      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
     607      !!---------------------------------------------------------------------- 
     608      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity 
     609      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m] 
     610      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio 
     611      ! 
     612      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     613      REAL(wp) ::   zt , zh , zs              ! local scalars 
     614      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     615      !!---------------------------------------------------------------------- 
     616      ! 
     617      IF( nn_timing == 1 ) CALL timing_start('rab_2d') 
     618      ! 
     619      pab(:,:,:) = 0._wp 
     620      ! 
     621      SELECT CASE ( nn_eos ) 
     622      ! 
     623      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     624         ! 
     625         DO jj = 1, jpjm1 
     626            DO ji = 1, fs_jpim1   ! vector opt. 
     627               ! 
     628               zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     629               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     630               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     631               ! 
     632               ! alpha 
     633               zn3 = ALP003 
     634               ! 
     635               zn2 = ALP012*zt + ALP102*zs+ALP002 
     636               ! 
     637               zn1 = ((ALP031*zt   & 
     638                  &   + ALP121*zs+ALP021)*zt   & 
     639                  &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     640                  &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     641                  ! 
     642               zn0 = ((((ALP050*zt   & 
     643                  &   + ALP140*zs+ALP040)*zt   & 
     644                  &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     645                  &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     646                  &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     647                  &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     648                  ! 
     649               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     650               ! 
     651               pab(ji,jj,jp_tem) = zn * r1_rau0 
     652               ! 
     653               ! beta 
     654               zn3 = BET003 
     655               ! 
     656               zn2 = BET012*zt + BET102*zs+BET002 
     657               ! 
     658               zn1 = ((BET031*zt   & 
     659                  &   + BET121*zs+BET021)*zt   & 
     660                  &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     661                  &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     662                  ! 
     663               zn0 = ((((BET050*zt   & 
     664                  &   + BET140*zs+BET040)*zt   & 
     665                  &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     666                  &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     667                  &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     668                  &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     669                  ! 
     670               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     671               ! 
     672               pab(ji,jj,jp_sal) = zn / zs * r1_rau0 
     673               ! 
     674               ! 
     675            END DO 
     676         END DO 
     677         ! 
     678         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions 
     679         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                     
     680         ! 
     681      CASE( 1 )                  !==  simplified EOS  ==! 
     682         ! 
     683         DO jj = 1, jpjm1 
     684            DO ji = 1, fs_jpim1   ! vector opt. 
     685               ! 
     686               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     687               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     688               zh    = pdep (ji,jj)                   ! depth at the partial step level 
     689               ! 
     690               zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     691               pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha 
     692               ! 
     693               zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     694               pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta 
     695               ! 
     696            END DO 
     697         END DO 
     698         ! 
     699         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions 
     700         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                     
     701         ! 
     702      CASE DEFAULT 
     703         IF(lwp) WRITE(numout,cform_err) 
     704         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     705         nstop = nstop + 1 
     706         ! 
     707      END SELECT 
     708      ! 
     709      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
     710         &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
     711      ! 
     712      IF( nn_timing == 1 )   CALL timing_stop('rab_2d') 
     713      ! 
     714   END SUBROUTINE rab_2d 
     715 
     716 
     717   SUBROUTINE bn2( pts, pab, pn2 ) 
     718      !!---------------------------------------------------------------------- 
     719      !!                  ***  ROUTINE bn2  *** 
     720      !! 
     721      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
     722      !!                time-step of the input arguments 
     723      !! 
     724      !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 
     725      !!      where alpha and beta are given in pab, and computed on T-points. 
     726      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module. 
     727      !! 
     728      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point  
     729      !! 
     730      !!---------------------------------------------------------------------- 
     731      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
     732      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1] 
     733      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
     734      ! 
     735      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     736      REAL(wp) ::   zaw, zbw, zrw   ! local scalars 
     737      !!---------------------------------------------------------------------- 
     738      ! 
     739      IF( nn_timing == 1 ) CALL timing_start('bn2') 
     740      ! 
     741      DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 ) 
     742         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90 
     743            DO ji = 1, jpi 
     744               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     745                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     746                  ! 
     747               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     748               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
     749               ! 
     750               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
     751                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
     752                  &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     753            END DO 
     754         END DO 
     755      END DO 
     756      ! 
     757      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
     758      ! 
     759      IF( nn_timing == 1 )   CALL timing_stop('bn2') 
     760      ! 
     761   END SUBROUTINE bn2 
     762 
     763 
     764   FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp ) 
     765      !!---------------------------------------------------------------------- 
     766      !!                 ***  ROUTINE eos_pt_from_ct  *** 
     767      !! 
     768      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celcius] 
     769      !! 
     770      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm 
     771      !!       checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC 
     772      !! 
     773      !! Reference  :   TEOS-10, UNESCO 
     774      !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC) 
     775      !!---------------------------------------------------------------------- 
     776      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp [Celcius] 
     777      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity   [psu] 
     778      ! Leave result array automatic rather than making explicitly allocated 
     779      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celcius] 
     780      ! 
     781      INTEGER  ::   ji, jj               ! dummy loop indices 
     782      REAL(wp) ::   zt , zs , ztm        ! local scalars 
     783      REAL(wp) ::   zn , zd              ! local scalars 
     784      REAL(wp) ::   zdeltaS , z1_S0 , z1_T0 
     785      !!---------------------------------------------------------------------- 
     786      ! 
     787      IF ( nn_timing == 1 )   CALL timing_start('eos_pt_from_ct') 
     788      ! 
     789      zdeltaS = 5._wp 
     790      z1_S0   = 0.875_wp/35.16504_wp 
     791      z1_T0   = 1._wp/40._wp 
     792      ! 
     793      DO jj = 1, jpj 
     794         DO ji = 1, jpi 
     795            ! 
     796            zt  = ctmp   (ji,jj) * z1_T0 
     797            zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
     798            ztm = tmask(ji,jj,1) 
     799            ! 
     800            zn = ((((-2.1385727895e-01_wp*zt   & 
     801               &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
     802               &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
     803               &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
     804               &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
     805               &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
     806               &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
     807               &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
     808               ! 
     809            zd = (2.0035003456_wp*zt   & 
     810               &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
     811               &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
     812               ! 
     813            ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
     814               ! 
     815         END DO 
     816      END DO 
     817      ! 
     818      IF( nn_timing == 1 )   CALL timing_stop('eos_pt_from_ct') 
     819      ! 
     820   END FUNCTION eos_pt_from_ct 
     821 
     822 
     823   FUNCTION eos_fzp( psal, pdep ) RESULT( ptf ) 
     824      !!---------------------------------------------------------------------- 
     825      !!                 ***  ROUTINE eos_fzp  *** 
     826      !! 
     827      !! ** Purpose :   Compute the freezing point temperature [Celcius] 
     828      !! 
     829      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by 
     830      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z 
     831      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m 
     832      !! 
     833      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
     834      !!---------------------------------------------------------------------- 
     835      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu] 
     836      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m] 
     837      REAL(wp), DIMENSION(jpi,jpj)                          ::   ptf   ! freezing temperature [Celcius] 
     838      ! 
     839      INTEGER  ::   ji, jj   ! dummy loop indices 
     840      REAL(wp) ::   zt, zs   ! local scalars 
     841      !!---------------------------------------------------------------------- 
     842      ! 
     843      SELECT CASE ( nn_eos ) 
     844      ! 
     845      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==! 
     846         ! 
     847         DO jj = 1, jpj 
     848            DO ji = 1, jpi 
     849               zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 )           ! square root salinity 
     850               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
     851                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
     852            END DO 
     853         END DO 
     854         ptf(:,:) = ptf(:,:) * psal(:,:) 
     855         ! 
     856         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
     857         ! 
     858      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==! 
     859         ! 
     860         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
     861            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
     862            ! 
     863         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
     864         ! 
     865      CASE DEFAULT 
     866         IF(lwp) WRITE(numout,cform_err) 
     867         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     868         nstop = nstop + 1 
     869         ! 
     870      END SELECT 
     871      ! 
     872   END FUNCTION eos_fzp 
     873 
     874 
     875   SUBROUTINE eos_pen( pts, pab_pe, ppen ) 
     876      !!---------------------------------------------------------------------- 
     877      !!                 ***  ROUTINE eos_pen  *** 
     878      !! 
     879      !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points 
     880      !! 
     881      !! ** Method  :   PE is defined analytically as the vertical  
     882      !!                   primitive of EOS times -g integrated between 0 and z>0. 
     883      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd 
     884      !!                                                      = 1/z * /int_0^z rd dz - rd  
     885      !!                                where rd is the density anomaly (see eos_rhd function) 
     886      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S: 
     887      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT 
     888      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS 
     889      !! 
     890      !! ** Action  : - pen         : PE anomaly given at T-points 
     891      !!            : - pab_pe  : given at T-points 
     892      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe 
     893      !!                    pab_pe(:,:,:,jp_sal) is beta_pe 
     894      !!---------------------------------------------------------------------- 
     895      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity 
     896      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe 
     897      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   ppen     ! potential energy anomaly 
     898      ! 
     899      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     900      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
     901      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      - 
     902      !!---------------------------------------------------------------------- 
     903      ! 
     904      IF( nn_timing == 1 )   CALL timing_start('eos_pen') 
     905      ! 
     906      SELECT CASE ( nn_eos ) 
     907      ! 
     908      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     909         ! 
     910         DO jk = 1, jpkm1 
     911            DO jj = 1, jpj 
     912               DO ji = 1, jpi 
     913                  ! 
     914                  zh  = fsdept(ji,jj,jk) * r1_Z0                                ! depth 
     915                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     916                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     917                  ztm = tmask(ji,jj,jk)                                         ! tmask 
     918                  ! 
     919                  ! potential energy non-linear anomaly 
     920                  zn2 = (PEN012)*zt   & 
     921                     &   + PEN102*zs+PEN002 
     922                     ! 
     923                  zn1 = ((PEN021)*zt   & 
     924                     &   + PEN111*zs+PEN011)*zt   & 
     925                     &   + (PEN201*zs+PEN101)*zs+PEN001 
     926                     ! 
     927                  zn0 = ((((PEN040)*zt   & 
     928                     &   + PEN130*zs+PEN030)*zt   & 
     929                     &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
     930                     &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
     931                     &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
     932                     ! 
     933                  zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     934                  ! 
     935                  ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm 
     936                  ! 
     937                  ! alphaPE non-linear anomaly 
     938                  zn2 = APE002 
     939                  ! 
     940                  zn1 = (APE011)*zt   & 
     941                     &   + APE101*zs+APE001 
     942                     ! 
     943                  zn0 = (((APE030)*zt   & 
     944                     &   + APE120*zs+APE020)*zt   & 
     945                     &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
     946                     &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
     947                     ! 
     948                  zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     949                  !                               
     950                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 
     951                  ! 
     952                  ! betaPE non-linear anomaly 
     953                  zn2 = BPE002 
     954                  ! 
     955                  zn1 = (BPE011)*zt   & 
     956                     &   + BPE101*zs+BPE001 
     957                     ! 
     958                  zn0 = (((BPE030)*zt   & 
     959                     &   + BPE120*zs+BPE020)*zt   & 
     960                     &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
     961                     &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
     962                     ! 
     963                  zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     964                  !                               
     965                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 
     966                  ! 
     967               END DO 
     968            END DO 
     969         END DO 
     970         ! 
     971      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     972         ! 
     973         DO jk = 1, jpkm1 
     974            DO jj = 1, jpj 
     975               DO ji = 1, jpi 
     976                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
     977                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
     978                  zh  = fsdept(ji,jj,jk)               ! depth in meters  at t-point 
     979                  ztm = tmask(ji,jj,jk)                ! tmask 
     980                  zn  = 0.5_wp * zh * r1_rau0 * ztm 
     981                  !                                    ! Potential Energy 
     982                  ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
     983                  !                                    ! alphaPE 
     984                  pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
     985                  pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
     986                  ! 
     987               END DO 
     988            END DO 
     989         END DO 
     990         ! 
     991      CASE DEFAULT 
     992         IF(lwp) WRITE(numout,cform_err) 
     993         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     994         nstop = nstop + 1 
     995         ! 
     996      END SELECT 
     997      ! 
     998      IF( nn_timing == 1 )   CALL timing_stop('eos_pen') 
     999      ! 
     1000   END SUBROUTINE eos_pen 
     1001 
     1002 
     1003   SUBROUTINE eos_init 
    6801004      !!---------------------------------------------------------------------- 
    6811005      !!                 ***  ROUTINE eos_init  *** 
    6821006      !! 
    683       !! ** Purpose :   Compute the sea surface freezing temperature [Celcius] 
    684       !! 
    685       !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???) 
    686       !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p 
    687       !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars 
    688       !! 
    689       !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
    690       !!---------------------------------------------------------------------- 
    691       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu] 
    692       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [decibars] 
    693       ! Leave result array automatic rather than making explicitly allocated 
    694       REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius] 
    695       !!---------------------------------------------------------------------- 
    696       ! 
    697       ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
    698          &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
    699       IF ( PRESENT( pdep ) ) THEN    
    700          ptf(:,:) = ptf(:,:) - 7.53e-4_wp * pdep(:,:) 
    701       ENDIF 
    702       ! 
    703    END FUNCTION tfreez 
    704  
    705  
    706    SUBROUTINE eos_init 
    707       !!---------------------------------------------------------------------- 
    708       !!                 ***  ROUTINE eos_init  *** 
    709       !! 
    7101007      !! ** Purpose :   initializations for the equation of state 
    7111008      !! 
    7121009      !! ** Method  :   Read the namelist nameos and control the parameters 
    7131010      !!---------------------------------------------------------------------- 
    714       NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta 
    715       !!---------------------------------------------------------------------- 
    716       INTEGER  ::   ios 
     1011      INTEGER  ::   ios   ! local integer 
     1012      !! 
     1013      NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   & 
     1014         &                                             rn_lambda2, rn_mu2, rn_nu 
     1015      !!---------------------------------------------------------------------- 
    7171016      ! 
    7181017      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state 
    7191018      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 
    7201019901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp ) 
    721  
     1020      ! 
    7221021      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state 
    7231022      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 
    7241023902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp ) 
    725       IF(lwm) WRITE( numond, nameos ) 
     1024      WRITE( numond, nameos ) 
     1025      ! 
     1026      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     1027      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
    7261028      ! 
    7271029      IF(lwp) THEN                ! Control print 
    … …  
    7311033         WRITE(numout,*) '          Namelist nameos : set eos parameters' 
    7321034         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos 
    733          WRITE(numout,*) '             thermal exp. coef. (linear)    rn_alpha = ', rn_alpha 
    734          WRITE(numout,*) '             saline  exp. coef. (linear)    rn_beta  = ', rn_beta 
     1035         IF( ln_useCT )   THEN 
     1036            WRITE(numout,*) '             model uses Conservative Temperature' 
     1037            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
     1038         ENDIF 
    7351039      ENDIF 
    7361040      ! 
    7371041      SELECT CASE( nn_eos )         ! check option 
    7381042      ! 
    739       CASE( 0 )                        !==  Jackett and McDougall (1994) formulation  ==! 
     1043      CASE( -1 )                       !==  polynomial TEOS-10  ==! 
    7401044         IF(lwp) WRITE(numout,*) 
    741          IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1994) equation of state and' 
    742          IF(lwp) WRITE(numout,*) '                 McDougall (1987) Brunt-Vaisala frequency' 
    743          ! 
    744       CASE( 1 )                        !==  Linear formulation = F( temperature )  ==! 
     1045         IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)' 
     1046         ! 
     1047         rdeltaS = 32._wp 
     1048         r1_S0  = 0.875_wp/35.16504_wp 
     1049         r1_T0  = 1._wp/40._wp 
     1050         r1_Z0  = 1.e-4_wp 
     1051         ! 
     1052         EOS000 = 8.0189615746e+02_wp 
     1053         EOS100 = 8.6672408165e+02_wp 
     1054         EOS200 = -1.7864682637e+03_wp 
     1055         EOS300 = 2.0375295546e+03_wp 
     1056         EOS400 = -1.2849161071e+03_wp 
     1057         EOS500 = 4.3227585684e+02_wp 
     1058         EOS600 = -6.0579916612e+01_wp 
     1059         EOS010 = 2.6010145068e+01_wp 
     1060         EOS110 = -6.5281885265e+01_wp 
     1061         EOS210 = 8.1770425108e+01_wp 
     1062         EOS310 = -5.6888046321e+01_wp 
     1063         EOS410 = 1.7681814114e+01_wp 
     1064         EOS510 = -1.9193502195_wp 
     1065         EOS020 = -3.7074170417e+01_wp 
     1066         EOS120 = 6.1548258127e+01_wp 
     1067         EOS220 = -6.0362551501e+01_wp 
     1068         EOS320 = 2.9130021253e+01_wp 
     1069         EOS420 = -5.4723692739_wp 
     1070         EOS030 = 2.1661789529e+01_wp 
     1071         EOS130 = -3.3449108469e+01_wp 
     1072         EOS230 = 1.9717078466e+01_wp 
     1073         EOS330 = -3.1742946532_wp 
     1074         EOS040 = -8.3627885467_wp 
     1075         EOS140 = 1.1311538584e+01_wp 
     1076         EOS240 = -5.3563304045_wp 
     1077         EOS050 = 5.4048723791e-01_wp 
     1078         EOS150 = 4.8169980163e-01_wp 
     1079         EOS060 = -1.9083568888e-01_wp 
     1080         EOS001 = 1.9681925209e+01_wp 
     1081         EOS101 = -4.2549998214e+01_wp 
     1082         EOS201 = 5.0774768218e+01_wp 
     1083         EOS301 = -3.0938076334e+01_wp 
     1084         EOS401 = 6.6051753097_wp 
     1085         EOS011 = -1.3336301113e+01_wp 
     1086         EOS111 = -4.4870114575_wp 
     1087         EOS211 = 5.0042598061_wp 
     1088         EOS311 = -6.5399043664e-01_wp 
     1089         EOS021 = 6.7080479603_wp 
     1090         EOS121 = 3.5063081279_wp 
     1091         EOS221 = -1.8795372996_wp 
     1092         EOS031 = -2.4649669534_wp 
     1093         EOS131 = -5.5077101279e-01_wp 
     1094         EOS041 = 5.5927935970e-01_wp 
     1095         EOS002 = 2.0660924175_wp 
     1096         EOS102 = -4.9527603989_wp 
     1097         EOS202 = 2.5019633244_wp 
     1098         EOS012 = 2.0564311499_wp 
     1099         EOS112 = -2.1311365518e-01_wp 
     1100         EOS022 = -1.2419983026_wp 
     1101         EOS003 = -2.3342758797e-02_wp 
     1102         EOS103 = -1.8507636718e-02_wp 
     1103         EOS013 = 3.7969820455e-01_wp 
     1104         ! 
     1105         ALP000 = -6.5025362670e-01_wp 
     1106         ALP100 = 1.6320471316_wp 
     1107         ALP200 = -2.0442606277_wp 
     1108         ALP300 = 1.4222011580_wp 
     1109         ALP400 = -4.4204535284e-01_wp 
     1110         ALP500 = 4.7983755487e-02_wp 
     1111         ALP010 = 1.8537085209_wp 
     1112         ALP110 = -3.0774129064_wp 
     1113         ALP210 = 3.0181275751_wp 
     1114         ALP310 = -1.4565010626_wp 
     1115         ALP410 = 2.7361846370e-01_wp 
     1116         ALP020 = -1.6246342147_wp 
     1117         ALP120 = 2.5086831352_wp 
     1118         ALP220 = -1.4787808849_wp 
     1119         ALP320 = 2.3807209899e-01_wp 
     1120         ALP030 = 8.3627885467e-01_wp 
     1121         ALP130 = -1.1311538584_wp 
     1122         ALP230 = 5.3563304045e-01_wp 
     1123         ALP040 = -6.7560904739e-02_wp 
     1124         ALP140 = -6.0212475204e-02_wp 
     1125         ALP050 = 2.8625353333e-02_wp 
     1126         ALP001 = 3.3340752782e-01_wp 
     1127         ALP101 = 1.1217528644e-01_wp 
     1128         ALP201 = -1.2510649515e-01_wp 
     1129         ALP301 = 1.6349760916e-02_wp 
     1130         ALP011 = -3.3540239802e-01_wp 
     1131         ALP111 = -1.7531540640e-01_wp 
     1132         ALP211 = 9.3976864981e-02_wp 
     1133         ALP021 = 1.8487252150e-01_wp 
     1134         ALP121 = 4.1307825959e-02_wp 
     1135         ALP031 = -5.5927935970e-02_wp 
     1136         ALP002 = -5.1410778748e-02_wp 
     1137         ALP102 = 5.3278413794e-03_wp 
     1138         ALP012 = 6.2099915132e-02_wp 
     1139         ALP003 = -9.4924551138e-03_wp 
     1140         ! 
     1141         BET000 = 1.0783203594e+01_wp 
     1142         BET100 = -4.4452095908e+01_wp 
     1143         BET200 = 7.6048755820e+01_wp 
     1144         BET300 = -6.3944280668e+01_wp 
     1145         BET400 = 2.6890441098e+01_wp 
     1146         BET500 = -4.5221697773_wp 
     1147         BET010 = -8.1219372432e-01_wp 
     1148         BET110 = 2.0346663041_wp 
     1149         BET210 = -2.1232895170_wp 
     1150         BET310 = 8.7994140485e-01_wp 
     1151         BET410 = -1.1939638360e-01_wp 
     1152         BET020 = 7.6574242289e-01_wp 
     1153         BET120 = -1.5019813020_wp 
     1154         BET220 = 1.0872489522_wp 
     1155         BET320 = -2.7233429080e-01_wp 
     1156         BET030 = -4.1615152308e-01_wp 
     1157         BET130 = 4.9061350869e-01_wp 
     1158         BET230 = -1.1847737788e-01_wp 
     1159         BET040 = 1.4073062708e-01_wp 
     1160         BET140 = -1.3327978879e-01_wp 
     1161         BET050 = 5.9929880134e-03_wp 
     1162         BET001 = -5.2937873009e-01_wp 
     1163         BET101 = 1.2634116779_wp 
     1164         BET201 = -1.1547328025_wp 
     1165         BET301 = 3.2870876279e-01_wp 
     1166         BET011 = -5.5824407214e-02_wp 
     1167         BET111 = 1.2451933313e-01_wp 
     1168         BET211 = -2.4409539932e-02_wp 
     1169         BET021 = 4.3623149752e-02_wp 
     1170         BET121 = -4.6767901790e-02_wp 
     1171         BET031 = -6.8523260060e-03_wp 
     1172         BET002 = -6.1618945251e-02_wp 
     1173         BET102 = 6.2255521644e-02_wp 
     1174         BET012 = -2.6514181169e-03_wp 
     1175         BET003 = -2.3025968587e-04_wp 
     1176         ! 
     1177         PEN000 = -9.8409626043_wp 
     1178         PEN100 = 2.1274999107e+01_wp 
     1179         PEN200 = -2.5387384109e+01_wp 
     1180         PEN300 = 1.5469038167e+01_wp 
     1181         PEN400 = -3.3025876549_wp 
     1182         PEN010 = 6.6681505563_wp 
     1183         PEN110 = 2.2435057288_wp 
     1184         PEN210 = -2.5021299030_wp 
     1185         PEN310 = 3.2699521832e-01_wp 
     1186         PEN020 = -3.3540239802_wp 
     1187         PEN120 = -1.7531540640_wp 
     1188         PEN220 = 9.3976864981e-01_wp 
     1189         PEN030 = 1.2324834767_wp 
     1190         PEN130 = 2.7538550639e-01_wp 
     1191         PEN040 = -2.7963967985e-01_wp 
     1192         PEN001 = -1.3773949450_wp 
     1193         PEN101 = 3.3018402659_wp 
     1194         PEN201 = -1.6679755496_wp 
     1195         PEN011 = -1.3709540999_wp 
     1196         PEN111 = 1.4207577012e-01_wp 
     1197         PEN021 = 8.2799886843e-01_wp 
     1198         PEN002 = 1.7507069098e-02_wp 
     1199         PEN102 = 1.3880727538e-02_wp 
     1200         PEN012 = -2.8477365341e-01_wp 
     1201         ! 
     1202         APE000 = -1.6670376391e-01_wp 
     1203         APE100 = -5.6087643219e-02_wp 
     1204         APE200 = 6.2553247576e-02_wp 
     1205         APE300 = -8.1748804580e-03_wp 
     1206         APE010 = 1.6770119901e-01_wp 
     1207         APE110 = 8.7657703198e-02_wp 
     1208         APE210 = -4.6988432490e-02_wp 
     1209         APE020 = -9.2436260751e-02_wp 
     1210         APE120 = -2.0653912979e-02_wp 
     1211         APE030 = 2.7963967985e-02_wp 
     1212         APE001 = 3.4273852498e-02_wp 
     1213         APE101 = -3.5518942529e-03_wp 
     1214         APE011 = -4.1399943421e-02_wp 
     1215         APE002 = 7.1193413354e-03_wp 
     1216         ! 
     1217         BPE000 = 2.6468936504e-01_wp 
     1218         BPE100 = -6.3170583896e-01_wp 
     1219         BPE200 = 5.7736640125e-01_wp 
     1220         BPE300 = -1.6435438140e-01_wp 
     1221         BPE010 = 2.7912203607e-02_wp 
     1222         BPE110 = -6.2259666565e-02_wp 
     1223         BPE210 = 1.2204769966e-02_wp 
     1224         BPE020 = -2.1811574876e-02_wp 
     1225         BPE120 = 2.3383950895e-02_wp 
     1226         BPE030 = 3.4261630030e-03_wp 
     1227         BPE001 = 4.1079296834e-02_wp 
     1228         BPE101 = -4.1503681096e-02_wp 
     1229         BPE011 = 1.7676120780e-03_wp 
     1230         BPE002 = 1.7269476440e-04_wp 
     1231         ! 
     1232      CASE( 0 )                        !==  polynomial EOS-80 formulation  ==! 
     1233         ! 
    7451234         IF(lwp) WRITE(numout,*) 
    746          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )' 
    747          IF( lk_zdfddm ) CALL ctl_stop( '          double diffusive mixing parameterization requires',   & 
    748               &                         ' that T and S are used as state variables' ) 
    749          ! 
    750       CASE( 2 )                        !==  Linear formulation = F( temperature , salinity )  ==! 
    751          ralpbet = rn_alpha / rn_beta 
    752          IF(lwp) WRITE(numout,*) 
    753          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 
     1235         IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)' 
     1236         ! 
     1237         rdeltaS = 20._wp 
     1238         r1_S0  = 1._wp/40._wp 
     1239         r1_T0  = 1._wp/40._wp 
     1240         r1_Z0  = 1.e-4_wp 
     1241         ! 
     1242         EOS000 = 9.5356891948e+02_wp 
     1243         EOS100 = 1.7136499189e+02_wp 
     1244         EOS200 = -3.7501039454e+02_wp 
     1245         EOS300 = 5.1856810420e+02_wp 
     1246         EOS400 = -3.7264470465e+02_wp 
     1247         EOS500 = 1.4302533998e+02_wp 
     1248         EOS600 = -2.2856621162e+01_wp 
     1249         EOS010 = 1.0087518651e+01_wp 
     1250         EOS110 = -1.3647741861e+01_wp 
     1251         EOS210 = 8.8478359933_wp 
     1252         EOS310 = -7.2329388377_wp 
     1253         EOS410 = 1.4774410611_wp 
     1254         EOS510 = 2.0036720553e-01_wp 
     1255         EOS020 = -2.5579830599e+01_wp 
     1256         EOS120 = 2.4043512327e+01_wp 
     1257         EOS220 = -1.6807503990e+01_wp 
     1258         EOS320 = 8.3811577084_wp 
     1259         EOS420 = -1.9771060192_wp 
     1260         EOS030 = 1.6846451198e+01_wp 
     1261         EOS130 = -2.1482926901e+01_wp 
     1262         EOS230 = 1.0108954054e+01_wp 
     1263         EOS330 = -6.2675951440e-01_wp 
     1264         EOS040 = -8.0812310102_wp 
     1265         EOS140 = 1.0102374985e+01_wp 
     1266         EOS240 = -4.8340368631_wp 
     1267         EOS050 = 1.2079167803_wp 
     1268         EOS150 = 1.1515380987e-01_wp 
     1269         EOS060 = -2.4520288837e-01_wp 
     1270         EOS001 = 1.0748601068e+01_wp 
     1271         EOS101 = -1.7817043500e+01_wp 
     1272         EOS201 = 2.2181366768e+01_wp 
     1273         EOS301 = -1.6750916338e+01_wp 
     1274         EOS401 = 4.1202230403_wp 
     1275         EOS011 = -1.5852644587e+01_wp 
     1276         EOS111 = -7.6639383522e-01_wp 
     1277         EOS211 = 4.1144627302_wp 
     1278         EOS311 = -6.6955877448e-01_wp 
     1279         EOS021 = 9.9994861860_wp 
     1280         EOS121 = -1.9467067787e-01_wp 
     1281         EOS221 = -1.2177554330_wp 
     1282         EOS031 = -3.4866102017_wp 
     1283         EOS131 = 2.2229155620e-01_wp 
     1284         EOS041 = 5.9503008642e-01_wp 
     1285         EOS002 = 1.0375676547_wp 
     1286         EOS102 = -3.4249470629_wp 
     1287         EOS202 = 2.0542026429_wp 
     1288         EOS012 = 2.1836324814_wp 
     1289         EOS112 = -3.4453674320e-01_wp 
     1290         EOS022 = -1.2548163097_wp 
     1291         EOS003 = 1.8729078427e-02_wp 
     1292         EOS103 = -5.7238495240e-02_wp 
     1293         EOS013 = 3.8306136687e-01_wp 
     1294         ! 
     1295         ALP000 = -2.5218796628e-01_wp 
     1296         ALP100 = 3.4119354654e-01_wp 
     1297         ALP200 = -2.2119589983e-01_wp 
     1298         ALP300 = 1.8082347094e-01_wp 
     1299         ALP400 = -3.6936026529e-02_wp 
     1300         ALP500 = -5.0091801383e-03_wp 
     1301         ALP010 = 1.2789915300_wp 
     1302         ALP110 = -1.2021756164_wp 
     1303         ALP210 = 8.4037519952e-01_wp 
     1304         ALP310 = -4.1905788542e-01_wp 
     1305         ALP410 = 9.8855300959e-02_wp 
     1306         ALP020 = -1.2634838399_wp 
     1307         ALP120 = 1.6112195176_wp 
     1308         ALP220 = -7.5817155402e-01_wp 
     1309         ALP320 = 4.7006963580e-02_wp 
     1310         ALP030 = 8.0812310102e-01_wp 
     1311         ALP130 = -1.0102374985_wp 
     1312         ALP230 = 4.8340368631e-01_wp 
     1313         ALP040 = -1.5098959754e-01_wp 
     1314         ALP140 = -1.4394226233e-02_wp 
     1315         ALP050 = 3.6780433255e-02_wp 
     1316         ALP001 = 3.9631611467e-01_wp 
     1317         ALP101 = 1.9159845880e-02_wp 
     1318         ALP201 = -1.0286156825e-01_wp 
     1319         ALP301 = 1.6738969362e-02_wp 
     1320         ALP011 = -4.9997430930e-01_wp 
     1321         ALP111 = 9.7335338937e-03_wp 
     1322         ALP211 = 6.0887771651e-02_wp 
     1323         ALP021 = 2.6149576513e-01_wp 
     1324         ALP121 = -1.6671866715e-02_wp 
     1325         ALP031 = -5.9503008642e-02_wp 
     1326         ALP002 = -5.4590812035e-02_wp 
     1327         ALP102 = 8.6134185799e-03_wp 
     1328         ALP012 = 6.2740815484e-02_wp 
     1329         ALP003 = -9.5765341718e-03_wp 
     1330         ! 
     1331         BET000 = 2.1420623987_wp 
     1332         BET100 = -9.3752598635_wp 
     1333         BET200 = 1.9446303907e+01_wp 
     1334         BET300 = -1.8632235232e+01_wp 
     1335         BET400 = 8.9390837485_wp 
     1336         BET500 = -1.7142465871_wp 
     1337         BET010 = -1.7059677327e-01_wp 
     1338         BET110 = 2.2119589983e-01_wp 
     1339         BET210 = -2.7123520642e-01_wp 
     1340         BET310 = 7.3872053057e-02_wp 
     1341         BET410 = 1.2522950346e-02_wp 
     1342         BET020 = 3.0054390409e-01_wp 
     1343         BET120 = -4.2018759976e-01_wp 
     1344         BET220 = 3.1429341406e-01_wp 
     1345         BET320 = -9.8855300959e-02_wp 
     1346         BET030 = -2.6853658626e-01_wp 
     1347         BET130 = 2.5272385134e-01_wp 
     1348         BET230 = -2.3503481790e-02_wp 
     1349         BET040 = 1.2627968731e-01_wp 
     1350         BET140 = -1.2085092158e-01_wp 
     1351         BET050 = 1.4394226233e-03_wp 
     1352         BET001 = -2.2271304375e-01_wp 
     1353         BET101 = 5.5453416919e-01_wp 
     1354         BET201 = -6.2815936268e-01_wp 
     1355         BET301 = 2.0601115202e-01_wp 
     1356         BET011 = -9.5799229402e-03_wp 
     1357         BET111 = 1.0286156825e-01_wp 
     1358         BET211 = -2.5108454043e-02_wp 
     1359         BET021 = -2.4333834734e-03_wp 
     1360         BET121 = -3.0443885826e-02_wp 
     1361         BET031 = 2.7786444526e-03_wp 
     1362         BET002 = -4.2811838287e-02_wp 
     1363         BET102 = 5.1355066072e-02_wp 
     1364         BET012 = -4.3067092900e-03_wp 
     1365         BET003 = -7.1548119050e-04_wp 
     1366         ! 
     1367         PEN000 = -5.3743005340_wp 
     1368         PEN100 = 8.9085217499_wp 
     1369         PEN200 = -1.1090683384e+01_wp 
     1370         PEN300 = 8.3754581690_wp 
     1371         PEN400 = -2.0601115202_wp 
     1372         PEN010 = 7.9263222935_wp 
     1373         PEN110 = 3.8319691761e-01_wp 
     1374         PEN210 = -2.0572313651_wp 
     1375         PEN310 = 3.3477938724e-01_wp 
     1376         PEN020 = -4.9997430930_wp 
     1377         PEN120 = 9.7335338937e-02_wp 
     1378         PEN220 = 6.0887771651e-01_wp 
     1379         PEN030 = 1.7433051009_wp 
     1380         PEN130 = -1.1114577810e-01_wp 
     1381         PEN040 = -2.9751504321e-01_wp 
     1382         PEN001 = -6.9171176978e-01_wp 
     1383         PEN101 = 2.2832980419_wp 
     1384         PEN201 = -1.3694684286_wp 
     1385         PEN011 = -1.4557549876_wp 
     1386         PEN111 = 2.2969116213e-01_wp 
     1387         PEN021 = 8.3654420645e-01_wp 
     1388         PEN002 = -1.4046808820e-02_wp 
     1389         PEN102 = 4.2928871430e-02_wp 
     1390         PEN012 = -2.8729602515e-01_wp 
     1391         ! 
     1392         APE000 = -1.9815805734e-01_wp 
     1393         APE100 = -9.5799229402e-03_wp 
     1394         APE200 = 5.1430784127e-02_wp 
     1395         APE300 = -8.3694846809e-03_wp 
     1396         APE010 = 2.4998715465e-01_wp 
     1397         APE110 = -4.8667669469e-03_wp 
     1398         APE210 = -3.0443885826e-02_wp 
     1399         APE020 = -1.3074788257e-01_wp 
     1400         APE120 = 8.3359333577e-03_wp 
     1401         APE030 = 2.9751504321e-02_wp 
     1402         APE001 = 3.6393874690e-02_wp 
     1403         APE101 = -5.7422790533e-03_wp 
     1404         APE011 = -4.1827210323e-02_wp 
     1405         APE002 = 7.1824006288e-03_wp 
     1406         ! 
     1407         BPE000 = 1.1135652187e-01_wp 
     1408         BPE100 = -2.7726708459e-01_wp 
     1409         BPE200 = 3.1407968134e-01_wp 
     1410         BPE300 = -1.0300557601e-01_wp 
     1411         BPE010 = 4.7899614701e-03_wp 
     1412         BPE110 = -5.1430784127e-02_wp 
     1413         BPE210 = 1.2554227021e-02_wp 
     1414         BPE020 = 1.2166917367e-03_wp 
     1415         BPE120 = 1.5221942913e-02_wp 
     1416         BPE030 = -1.3893222263e-03_wp 
     1417         BPE001 = 2.8541225524e-02_wp 
     1418         BPE101 = -3.4236710714e-02_wp 
     1419         BPE011 = 2.8711395266e-03_wp 
     1420         BPE002 = 5.3661089288e-04_wp 
     1421         ! 
     1422      CASE( 1 )                        !==  Simplified EOS     ==! 
     1423         IF(lwp) THEN 
     1424            WRITE(numout,*) 
     1425            WRITE(numout,*) '          use of simplified eos:    rhd(dT=T-10,dS=S-35,Z) = ' 
     1426            WRITE(numout,*) '             [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0' 
     1427            WRITE(numout,*) 
     1428            WRITE(numout,*) '             thermal exp. coef.    rn_a0      = ', rn_a0 
     1429            WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0 
     1430            WRITE(numout,*) '             cabbeling coef.       rn_lambda1 = ', rn_lambda1 
     1431            WRITE(numout,*) '             cabbeling coef.       rn_lambda2 = ', rn_lambda2 
     1432            WRITE(numout,*) '             thermobar. coef.      rn_mu1     = ', rn_mu1 
     1433            WRITE(numout,*) '             thermobar. coef.      rn_mu2     = ', rn_mu2 
     1434            WRITE(numout,*) '             2nd cabbel. coef.     rn_nu      = ', rn_nu 
     1435            WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization ' 
     1436         ENDIF 
    7541437         ! 
    7551438      CASE DEFAULT                     !==  ERROR in nn_eos  ==! 
    … …  
    7591442      END SELECT 
    7601443      ! 
     1444      r1_rau0     = 1._wp / rau0 
     1445      r1_rcp      = 1._wp / rcp 
     1446      r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 
     1447      ! 
     1448      IF(lwp) WRITE(numout,*) 
     1449      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m^3' 
     1450      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg' 
     1451      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
     1452      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
     1453      ! 
    7611454   END SUBROUTINE eos_init 
    7621455 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r4499 r4933  
    44   !! Ocean  tracers:  horizontal & vertical advective trend 
    55   !!====================================================================== 
    6    !! History :  8.2  ! 2001-08  (G. Madec, E. Durand) trahad+trazad=traadv  
    7    !!            1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  ! 2004-08  (C. Talandier) New trends organization 
     6   !! History :  OPA  ! 2001-08  (G. Madec, E. Durand) v8.2 trahad+trazad=traadv  
     7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
     8   !!             -   ! 2004-08  (C. Talandier) New trends organization 
    99   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization 
    1010   !!            2.0  ! 2006-04  (R. Benshila, G. Madec) Step reorganization 
    … …  
    2121   USE dom_oce         ! ocean space and time domain 
    2222   USE eosbn2          ! equation of state 
    23    USE trdmod_oce      ! tracers trends 
    24    USE trdtra          ! tracers trends 
     23   USE trd_oce         ! trends: ocean variables 
     24   USE trdtra          ! trends manager: tracers  
    2525   USE closea          ! closed sea 
    2626   USE sbcrnf          ! river runoffs 
    … …  
    3737   PRIVATE 
    3838 
    39    PUBLIC   tra_adv_cen2       ! routine called by step.F90 
    40    PUBLIC   ups_orca_set       ! routine used by traadv_cen2_jki.F90 
    41  
    42    LOGICAL  :: l_trd       ! flag to compute trends 
     39   PUBLIC   tra_adv_cen2   ! routine called by traadv.F90 
    4340 
    4441   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits  
    … …  
    5552 
    5653   SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn,     & 
    57       &                                 ptb, ptn, pta, kjpt   )  
     54      &                                         ptb, ptn, pta, kjpt   )  
    5855      !!---------------------------------------------------------------------- 
    5956      !!                  ***  ROUTINE tra_adv_cen2  *** 
    … …  
    8582      !!       * Add this trend now to the general trend of tracer (ta,sa): 
    8683      !!               pta = pta + ztra 
    87       !!       * trend diagnostic ('key_trdtra' defined): the trend is 
     84      !!       * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 
    8885      !!      saved for diagnostics. The trends saved is expressed as 
    89       !!      Uh.gradh(T), i.e. 
    90       !!                     save trend = ztra + ptn divn 
     86      !!      Uh.gradh(T), i.e.  save trend = ztra + ptn divn 
    9187      !! 
    9288      !!         Part II : vertical advection 
    … …  
    104100      !!         Add this trend now to the general trend of tracer (ta,sa): 
    105101      !!             pta = pta + ztra 
    106       !!         Trend diagnostic ('key_trdtra' defined): the trend is 
     102      !!         Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 
    107103      !!      saved for diagnostics. The trends saved is expressed as : 
    108104      !!             save trend =  w.gradz(T) = ztra - ptn divn. 
    … …  
    111107      !!              - save trends if needed 
    112108      !!---------------------------------------------------------------------- 
    113       USE oce     , ONLY:   zwx => ua        , zwy  => va          ! (ua,va) used as 3D workspace 
    114       ! 
    115109      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    116110      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    … …  
    128122      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      - 
    129123      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      - 
    130       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez  
    131       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind 
     124      REAL(wp), POINTER, DIMENSION(:,:)   :: zfzp         ! 2D workspace 
     125      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy     ! 3D     - 
     126      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind    !  -     - 
    132127      !!---------------------------------------------------------------------- 
    133128      ! 
    134129      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2') 
    135130      ! 
    136       CALL wrk_alloc( jpi, jpj, ztfreez ) 
    137       CALL wrk_alloc( jpi, jpj, jpk, zwz, zind ) 
     131      CALL wrk_alloc( jpi, jpj, zfzp ) 
     132      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwz, zind ) 
    138133      ! 
    139134 
    … …  
    144139         IF(lwp) WRITE(numout,*) 
    145140         ! 
    146          IF ( .NOT. ALLOCATED( upsmsk ) )  THEN 
     141         IF( .NOT. ALLOCATED( upsmsk ) )  THEN 
    147142             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 
    148143             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array') 
    … …  
    162157      ENDIF 
    163158      ! 
    164       l_trd = .FALSE. 
    165       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    166       ! 
    167159      ! Upstream / centered scheme indicator 
    168160      ! ------------------------------------ 
    169161!!gm  not strickly exact : the freezing point should be computed at each ocean levels... 
    170162!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations 
    171       ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) ) 
     163      zfzp(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) 
    172164      DO jk = 1, jpk 
    173165         DO jj = 1, jpj 
    174166            DO ji = 1, jpi 
    175167               !                                        ! below ice covered area (if tn < "freezing"+0.1 ) 
    176                IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0 
    177                ELSE                                                      ;   zice = 0.e0 
     168               IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN   ;   zice = 1._wp 
     169               ELSE                                                   ;   zice = 0._wp 
    178170               ENDIF 
    179171               zind(ji,jj,jk) = MAX (   & 
    … …  
    260252         END DO 
    261253 
    262          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    263          IF( l_trd ) THEN 
    264             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
    265             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    266             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     254         !                                 ! trend diagnostics 
     255         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.    & 
     256            &( cdtype == 'TRC' .AND. l_trdtrc ) )   THEN 
     257            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     258            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     259            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    267260         END IF 
    268261         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    269262         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    270            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    271            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     263           IF( jn == jp_tem )   htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
     264           IF( jn == jp_sal )   str_adv(:) = ptr_vj( zwy(:,:,:) ) 
    272265         ENDIF 
    273266         ! 
    274       ENDDO 
     267      END DO 
    275268 
    276269      ! ---------------------------  required in restart file to ensure restartability) 
    … …  
    281274      ENDIF 
    282275      ! 
    283       CALL wrk_dealloc( jpi, jpj, ztfreez ) 
    284       CALL wrk_dealloc( jpi, jpj, jpk, zwz, zind ) 
     276      CALL wrk_dealloc( jpi, jpj, zfzp ) 
     277      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwz, zind ) 
    285278      ! 
    286279      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2') 
    … …  
    303296      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers 
    304297      !!---------------------------------------------------------------------- 
    305        
    306298      ! 
    307299      IF( nn_timing == 1 )  CALL timing_start('ups_orca_set') 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

    r4499 r4933  
    1616   !!---------------------------------------------------------------------- 
    1717   USE oce            ! ocean dynamics and active tracers 
     18   USE trc_oce        ! share passive tracers/Ocean variables 
    1819   USE dom_oce        ! ocean space and time domain 
    19    USE trdmod_oce     ! tracers trends  
    20    USE trdtra         ! tracers trends  
    21    USE in_out_manager ! I/O manager 
     20   USE trd_oce        ! trends: ocean variables 
     21   USE trdtra         ! tracers trends manager 
    2222   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    23    USE trabbl         ! tracers: bottom boundary layer 
    24    USE lib_mpp        ! distribued memory computing 
    25    USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
     23   USE sbcrnf          ! river runoffs 
    2624   USE diaptr         ! poleward transport diagnostics 
    27    USE trc_oce        ! share passive tracers/Ocean variables 
     25   ! 
    2826   USE wrk_nemo       ! Memory Allocation 
    2927   USE timing         ! Timing 
    3028   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    31    USE eosbn2          ! equation of state 
    32    USE sbcrnf          ! river runoffs 
     29   USE in_out_manager ! I/O manager 
     30   USE lib_mpp        ! distribued memory computing 
     31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    3332 
    3433   IMPLICIT NONE 
    3534   PRIVATE 
    3635 
    37    PUBLIC   tra_adv_muscl  ! routine called by step.F90 
    38  
    39    LOGICAL  :: l_trd                        ! flag to compute trends 
    40    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits 
    41    !                                                             !  and in closed seas (orca 2 and 4 configurations) 
    42    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind         !: mixed upstream/centered index 
     36   PUBLIC   tra_adv_muscl   ! routine called by traadv.F90 
     37    
     38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits 
     39   !                                                           !  and in closed seas (orca 2 and 4 configurations) 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index 
     41    
    4342   !! * Substitutions 
    4443#  include "domzgr_substitute.h90" 
    … …  
    5150CONTAINS 
    5251 
    53    SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 
    54       &                                        ptb, pta, kjpt, ld_msc_ups ) 
     52   SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   & 
     53      &                                                ptb, pta, kjpt, ld_msc_ups ) 
    5554      !!---------------------------------------------------------------------- 
    5655      !!                    ***  ROUTINE tra_adv_muscl  *** 
    … …  
    6867      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 
    6968      !!---------------------------------------------------------------------- 
    70       USE oce     , ONLY:   zwx   => ua    , zwy   => va          ! (ua,va) used as workspace 
    71       ! 
    7269      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    7370      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    … …  
    7976      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field 
    8077      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    81  
    82       ! 
    83       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     78      ! 
     79      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices 
     80      INTEGER  ::   ierr                      ! local integer 
    8481      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars 
    8582      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      - 
    8683      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      - 
    87       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy 
    88       INTEGER  ::   ierr 
     84      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace 
     85      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -  
    8986      !!---------------------------------------------------------------------- 
    9087      ! 
    9188      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl') 
    9289      ! 
    93       CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy ) 
    94       ! 
    95  
     90      CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
     91      ! 
    9692      IF( kt == kit000 )  THEN 
    9793         IF(lwp) WRITE(numout,*) 
    … …  
    117113 
    118114         ! 
    119          ! Upstream / centered scheme indicator 
     115         ! Upstream / MUSCL scheme indicator 
    120116         ! ------------------------------------ 
     117!!gm  useless 
    121118         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed 
     119!!gm 
    122120         ! 
    123121         IF( ld_msc_ups )  THEN 
    124             DO jk = 1, jpk 
    125                DO jj = 1, jpj 
    126                   DO ji = 1, jpi 
    127                      xind(ji,jj,jk) = 1  - MAX (           & 
    128                         rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows) 
    129                         upsmsk(ji,jj) ) * tmask(ji,jj,jk)     ! some of some straits 
    130                   END DO 
    131                END DO 
     122            DO jk = 1, jpkm1 
     123               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed 
     124                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows) 
     125                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits 
    132126            END DO 
    133127         ENDIF  
    134128         ! 
    135129      ENDIF  
    136       ! 
    137       l_trd = .FALSE. 
    138       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    139        
     130      !       
    140131      !                                                     ! =========== 
    141132      DO jn = 1, kjpt                                       ! tracer loop 
    … …  
    192183                  zalpha = 0.5 - z0u 
    193184                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 
    194                   zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk)) 
    195                   zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk)) 
     185                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
     186                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
    196187                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    197188                  ! 
    … …  
    199190                  zalpha = 0.5 - z0v 
    200191                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    201                   zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk)) 
    202                   zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk)) 
     192                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
     193                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
    203194                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    204195               END DO 
    … …  
    222213         END DO         
    223214         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    224          IF( l_trd )  THEN 
    225             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) ) 
    226             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) ) 
     215         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
     216            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN 
     217            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
     218            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    227219         END IF 
    228220         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    … …  
    274266                  zalpha = 0.5 + z0w 
    275267                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr  
    276                   zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1)) 
    277                   zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  )) 
     268                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
     269                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    278270                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    279271               END DO  
    … …  
    281273         END DO 
    282274 
    283          ! Compute & add the vertical advective trend 
    284          DO jk = 1, jpkm1 
     275         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend 
    285276            DO jj = 2, jpjm1       
    286277               DO ji = fs_2, fs_jpim1   ! vector opt. 
    287                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     278                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    288279                  ! vertical advective trends  
    289280                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    … …  
    294285         END DO 
    295286         !                                 ! Save the vertical advective trends for diagnostic 
    296          IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    297          ! 
    298       ENDDO 
    299       ! 
    300       CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy ) 
     287         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
     288            &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
     289            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     290         ! 
     291      END DO 
     292      ! 
     293      CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
    301294      ! 
    302295      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl') 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

    r4499 r4933  
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce             ! ocean dynamics and active tracers 
     15   USE trc_oce         ! share passive tracers/Ocean variables 
    1516   USE dom_oce         ! ocean space and time domain 
    16    USE trdmod_oce      ! tracers trends  
    17    USE trdtra          ! tracers trends  
     17   USE trd_oce         ! trends: ocean variables 
     18   USE trdtra          ! trends manager: tracers  
    1819   USE in_out_manager  ! I/O manager 
    1920   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
    20    USE trabbl          ! tracers: bottom boundary layer 
     21   USE diaptr          ! poleward transport diagnostics 
     22   ! 
    2123   USE lib_mpp         ! distribued memory computing 
    2224   USE lbclnk          ! ocean lateral boundary condition (or mpp link)  
    23    USE diaptr          ! poleward transport diagnostics 
    24    USE trc_oce         ! share passive tracers/Ocean variables 
    2525   USE wrk_nemo        ! Memory Allocation 
    2626   USE timing          ! Timing 
    … …  
    3131 
    3232   PUBLIC   tra_adv_muscl2        ! routine called by step.F90 
    33  
    34    LOGICAL  :: l_trd       ! flag to compute trends 
    3533 
    3634   !! * Substitutions 
    … …  
    6159      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 
    6260      !!---------------------------------------------------------------------- 
    63       USE oce     , ONLY:   zwx   => ua    , zwy   => va         ! (ua,va) used as 3D workspace 
    64       !! 
    6561      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    6662      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    … …  
    7672      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      - 
    7773      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      - 
    78       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy 
     74      REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy , zwx, zwy 
    7975      !!---------------------------------------------------------------------- 
    8076      ! 
    8177      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl2') 
    8278      ! 
    83       CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy ) 
     79      CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
    8480      ! 
    8581 
    … …  
    9086      ENDIF 
    9187      ! 
    92       l_trd = .FALSE. 
    93       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    94  
    9588      !                                                          ! =========== 
    9689      DO jn = 1, kjpt                                            ! tracer loop 
    … …  
    200193         END DO 
    201194         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    202          IF( l_trd ) THEN 
    203             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) ) 
    204             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) ) 
     195         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
     196            &( cdtype == 'TRC' .AND. l_trdtrc )      ) THEN 
     197            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
     198            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    205199         END IF 
    206200 
    … …  
    284278         END DO 
    285279         !                       ! trend diagnostics (contribution of upstream fluxes) 
    286          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     280         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
     281            &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
     282            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    287283         ! 
    288284      END DO 
    289285      ! 
    290       CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy ) 
     286      CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
    291287      ! 
    292288      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl2') 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r4499 r4933  
    1717   USE oce             ! ocean dynamics and active tracers 
    1818   USE dom_oce         ! ocean space and time domain 
    19    USE trdmod_oce      ! ocean space and time domain 
    20    USE trdtra          ! ocean tracers trends  
    21    USE trabbl          ! advective term in the BBL 
     19   USE trc_oce         ! share passive tracers/Ocean variables 
     20   USE trd_oce         ! trends: ocean variables 
     21   USE trdtra          ! trends manager: tracers  
     22   USE dynspg_oce      ! surface pressure gradient variables 
     23   USE diaptr          ! poleward transport diagnostics 
     24   ! 
    2225   USE lib_mpp         ! distribued memory computing 
    2326   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    24    USE dynspg_oce      ! surface pressure gradient variables 
    2527   USE in_out_manager  ! I/O manager 
    26    USE diaptr          ! poleward transport diagnostics 
    27    USE trc_oce         ! share passive tracers/Ocean variables 
    2828   USE wrk_nemo        ! Memory Allocation 
    2929   USE timing          ! Timing 
    … …  
    9393      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    9494      !!---------------------------------------------------------------------- 
    95  
    9695      ! 
    9796      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck') 
    … …  
    103102         IF(lwp) WRITE(numout,*) 
    104103      ENDIF 
    105       ! 
    106104      l_trd = .FALSE. 
    107       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    108  
     105      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
     106      ! 
    109107      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
    110108      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )  
    … …  
    124122      !! 
    125123      !!---------------------------------------------------------------------- 
    126       USE oce     , ONLY:   zwx => ua       ! ua used as workspace 
    127       ! 
    128124      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    129125      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    … …  
    136132      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    137133      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
    138       REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd 
     134      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 
    139135      !---------------------------------------------------------------------- 
    140136      ! 
    141       CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd ) 
     137      CALL wrk_alloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd ) 
    142138      !                                                          ! =========== 
    143139      DO jn = 1, kjpt                                            ! tracer loop 
    … …  
    233229         END DO 
    234230         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    235          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
     231         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
    236232         ! 
    237233      END DO 
    238234      ! 
    239       CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd ) 
     235      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd ) 
    240236      ! 
    241237   END SUBROUTINE tra_adv_qck_i 
    … …  
    247243      !! 
    248244      !!---------------------------------------------------------------------- 
    249       USE oce     , ONLY:   zwy => ua       ! ua used as workspace 
    250       ! 
    251245      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    252246      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    … …  
    259253      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    260254      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
    261       REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd 
     255      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd 
    262256      !---------------------------------------------------------------------- 
    263257      ! 
    264       CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd ) 
     258      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd ) 
    265259      ! 
    266260      !                                                          ! =========== 
    … …  
    359353         END DO 
    360354         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    361          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     355         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    362356         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    363357         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    … …  
    368362      END DO 
    369363      ! 
    370       CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd ) 
     364      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd ) 
    371365      ! 
    372366   END SUBROUTINE tra_adv_qck_j 
    … …  
    378372      !! 
    379373      !!---------------------------------------------------------------------- 
    380       USE oce, ONLY:   zwz => ua   ! ua used as workspace 
    381       ! 
    382374      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    383375      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    … …  
    389381      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    390382      REAL(wp) ::   zbtr , ztra      ! local scalars 
    391       !!---------------------------------------------------------------------- 
    392  
     383      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz 
     384      !!---------------------------------------------------------------------- 
     385      ! 
     386      CALL wrk_alloc( jpi, jpj, jpk, zwz ) 
    393387      !                                                          ! =========== 
    394388      DO jn = 1, kjpt                                            ! tracer loop 
    … …  
    422416         END DO 
    423417         !                                 ! Save the vertical advective trends for diagnostic 
    424          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     418         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    425419         ! 
    426420      END DO 
     421      ! 
     422      CALL wrk_dealloc( jpi, jpj, jpk, zwz ) 
    427423      ! 
    428424   END SUBROUTINE tra_adv_cen2_k 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4499 r4933  
    2222   USE oce            ! ocean dynamics and active tracers 
    2323   USE dom_oce        ! ocean space and time domain 
    24    USE trdmod_oce     ! tracers trends 
     24   USE trc_oce        ! share passive tracers/Ocean variables 
     25   USE trd_oce        ! trends: ocean variables 
    2526   USE trdtra         ! tracers trends 
    26    USE in_out_manager ! I/O manager 
    2727   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
     28   USE diaptr         ! poleward transport diagnostics 
     29   ! 
    2830   USE lib_mpp        ! MPP library 
    2931   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    30    USE diaptr         ! poleward transport diagnostics 
    31    USE trc_oce        ! share passive tracers/Ocean variables 
     32   USE in_out_manager ! I/O manager 
    3233   USE wrk_nemo       ! Memory Allocation 
    3334   USE timing         ! Timing 
    … …  
    9394         IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype 
    9495         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     96         ! 
     97         l_trd = .FALSE. 
     98         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    9599      ENDIF 
    96       ! 
    97       l_trd = .FALSE. 
    98       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    99100      ! 
    100101      IF( l_trd )  THEN 
    … …  
    228229            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    229230             
    230             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )    
    231             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
    232             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
     231            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )    
     232            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
     233            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
    233234         END IF 
    234235         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    … …  
    261262      !!       in-space based differencing for fluid 
    262263      !!---------------------------------------------------------------------- 
    263       ! 
    264       !!---------------------------------------------------------------------- 
    265264      REAL(wp), DIMENSION(jpk)         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    266265      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
    267266      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions 
    268267      ! 
    269       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    270       INTEGER ::   ikm1         ! local integer 
     268      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     269      INTEGER  ::   ikm1         ! local integer 
    271270      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars 
    272271      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    … …  
    278277      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 
    279278      ! 
    280  
    281279      zbig  = 1.e+40_wp 
    282280      zrtrn = 1.e-15_wp 
    283281      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp 
    284282 
    285  
    286283      ! Search local extrema 
    287284      ! -------------------- 
    288285      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    289       zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   & 
    290          &        paft * tmask - zbig * ( 1.e0 - tmask )  ) 
    291       zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   & 
    292          &        paft * tmask + zbig * ( 1.e0 - tmask )  ) 
     286      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   & 
     287         &        paft * tmask - zbig * ( 1._wp - tmask )  ) 
     288      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   & 
     289         &        paft * tmask + zbig * ( 1._wp - tmask )  ) 
    293290 
    294291      DO jk = 1, jpkm1 
    … …  
    334331         DO jj = 2, jpjm1 
    335332            DO ji = fs_2, fs_jpim1   ! vector opt. 
    336                zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    337                zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     333               zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     334               zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
    338335               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
    339                paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 
    340  
    341                zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
    342                zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     336               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
     337 
     338               zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     339               zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
    343340               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
    344                pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 
     341               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    345342 
    346343      ! monotonic flux in the k direction, i.e. pcc 
    … …  
    349346               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
    350347               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 
    351                pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb ) 
     348               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 
    352349            END DO 
    353350         END DO 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r4499 r4933  
    1414   USE oce            ! ocean dynamics and active tracers 
    1515   USE dom_oce        ! ocean space and time domain 
    16    USE trdmod_oce     ! ocean space and time domain 
    17    USE trdtra 
    18    USE lib_mpp 
     16   USE trc_oce        ! share passive tracers/Ocean variables 
     17   USE trd_oce        ! trends: ocean variables 
     18   USE trdtra         ! trends manager: tracers  
     19   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
     20   USE diaptr         ! poleward transport diagnostics 
     21   ! 
     22   USE lib_mpp        ! I/O library 
    1923   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    2024   USE in_out_manager ! I/O manager 
    21    USE diaptr         ! poleward transport diagnostics 
    22    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    23    USE trc_oce        ! share passive tracers/Ocean variables 
    2425   USE wrk_nemo       ! Memory Allocation 
    2526   USE timing         ! Timing 
    … …  
    5152      !!      and add it to the general trend of passive tracer equations. 
    5253      !! 
    53       !! ** Method  :   The upstream biased 3rd order scheme (UBS) is based on an 
     54      !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order 
    5455      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 
    5556      !!      It is only used in the horizontal direction. 
    5657      !!      For example the i-component of the advective fluxes are given by : 
    5758      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0 
    58       !!          zwx = !  or  
     59      !!          ztu = !  or  
    5960      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0 
    6061      !!      where zltu is the second derivative of the before temperature field: 
    … …  
    7677      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.  
    7778      !!---------------------------------------------------------------------- 
    78       USE oce     , ONLY:   zwx  => ua       , zwy  => va         ! (ua,va) used as workspace 
    79       ! 
    8079      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    8180      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    … …  
    9897      CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
    9998      ! 
    100  
    10199      IF( kt == kit000 )  THEN 
    102100         IF(lwp) WRITE(numout,*) 
    … …  
    151149                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) ) 
    152150                  ! UBS advective fluxes 
    153                   zwx(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
    154                   zwy(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
     151                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
     152                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
    155153               END DO 
    156154            END DO 
    … …  
    159157         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends 
    160158 
    161          ! Horizontal advective trends 
    162          DO jk = 1, jpkm1 
    163             !  Tracer flux divergence at t-point added to the general trend 
     159         DO jk = 1, jpkm1                 ! Horizontal advective trends 
    164160            DO jj = 2, jpjm1 
    165161               DO ji = fs_2, fs_jpim1   ! vector opt. 
    166                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    167                   ! horizontal advective 
    168                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    169                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    170                   ! add it to the general tracer trends 
    171                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     162                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                        & 
     163                     &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    & 
     164                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    172165               END DO 
    173166            END DO 
    … …  
    178171         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 
    179172 
    180          ! 3. Save the horizontal advective trends for diagnostic 
    181          ! ------------------------------------------------------ 
    182          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    183          IF( l_trd ) THEN 
    184              CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
    185              CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     173         !                 
     174         IF( l_trd ) THEN                  ! trend diagnostics 
     175             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) ) 
     176             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) ) 
    186177         END IF 
    187178         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    188179         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    189             IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    190             IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     180            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( ztv(:,:,:) ) 
     181            IF( jn == jp_sal )  str_adv(:) = ptr_vj( ztv(:,:,:) ) 
    191182         ENDIF 
    192183          
    … …  
    265256               END DO 
    266257            END DO 
    267             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv ) 
     258            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 
    268259         ENDIF 
    269260         ! 
    270       ENDDO 
     261      END DO 
    271262      ! 
    272263      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
    … …  
    290281      !!       in-space based differencing for fluid 
    291282      !!---------------------------------------------------------------------- 
    292       ! 
    293283      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt   ! vertical profile of tracer time-step 
    294284      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field 
    … …  
    306296      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo ) 
    307297      ! 
    308  
    309298      zbig  = 1.e+40_wp 
    310299      zrtrn = 1.e-15_wp 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r4624 r4933  
    1818   USE dom_oce         ! domain: ocean 
    1919   USE phycst          ! physical constants 
    20    USE trdmod_oce      ! trends: ocean variables  
    21    USE trdtra          ! trends: active tracers  
     20   USE trd_oce         ! trends: ocean variables 
     21   USE trdtra          ! trends manager: tracers  
    2222   USE in_out_manager  ! I/O manager 
    2323   USE prtctl          ! Print control 
    … …  
    8484      ! 
    8585      !                             !  Add the geothermal heat flux trend on temperature 
    86 #if defined key_vectopt_loop 
    87       DO jj = 1, 1 
    88          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    89 #else 
    9086      DO jj = 2, jpjm1 
    9187         DO ji = 2, jpim1 
    92 #endif 
    9388            ik = mbkt(ji,jj) 
    9489            zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik) 
    … …  
    9994      IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics 
    10095         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    101          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbc, ztrdt ) 
     96         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 
    10297         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
    10398      ENDIF 
    … …  
    130125      INTEGER  ::   inum                ! temporary logical unit 
    131126      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    132       !! 
     127      ! 
    133128      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst  
    134129      !!---------------------------------------------------------------------- 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4624 r4933  
    1212   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC 
    1313   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
     14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta 
    1415   !!---------------------------------------------------------------------- 
    1516#if   defined key_trabbl   ||   defined key_esopa 
    … …  
    2829   USE phycst         ! physical constant 
    2930   USE eosbn2         ! equation of state 
    30    USE trdmod_oce     ! trends: ocean variables 
     31   USE trd_oce     ! trends: ocean variables 
    3132   USE trdtra         ! trends: active tracers 
    32    USE iom            ! IOM server 
     33   ! 
     34   USE iom            ! IOM library                
    3335   USE in_out_manager ! I/O manager 
    3436   USE lbclnk         ! ocean lateral boundary conditions 
    … …  
    3638   USE wrk_nemo       ! Memory Allocation 
    3739   USE timing         ! Timing 
    38  
     40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3941 
    4042   IMPLICIT NONE 
    … …  
    5759   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s] 
    5860 
    59    LOGICAL , PUBLIC ::   l_bbl                  !: flag to compute bbl diffu. flux coef and transport 
     61   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport 
    6062 
    6163   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer 
    … …  
    8486         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     & 
    8587         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          & 
    86          &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc                ) 
     88         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc ) 
    8789         ! 
    8890      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc ) 
    … …  
    104106      !!---------------------------------------------------------------------- 
    105107      INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
    106       !! 
     108      ! 
    107109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    108110      !!---------------------------------------------------------------------- 
    … …  
    110112      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl') 
    111113      ! 
    112       IF( l_trdtra )   THEN                        !* Save ta and sa trends 
     114      IF( l_trdtra )   THEN                         !* Save ta and sa trends 
    113115         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    114116         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    … …  
    116118      ENDIF 
    117119 
    118       IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
    119  
    120       IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl 
     120      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
     121 
     122      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl 
    121123         ! 
    122124         CALL tra_bbl_dif( tsb, tsa, jpts ) 
    123125         IF( ln_ctl )  & 
    124126         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
    125          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     127            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    126128         ! lateral boundary conditions ; just need for outputs 
    127129         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. ) 
    … …  
    131133      END IF 
    132134 
    133       IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl 
     135      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    134136         ! 
    135137         CALL tra_bbl_adv( tsb, tsa, jpts ) 
    136138         IF(ln_ctl)   & 
    137139         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   & 
    138          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     140            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    139141         ! lateral boundary conditions ; just need for outputs 
    140142         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. ) 
    … …  
    147149         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    148150         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    149          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt ) 
    150          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds ) 
     151         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 
     152         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 
    151153         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    152154      ENDIF 
    … …  
    164166      !!                advection terms. 
    165167      !! 
    166       !! ** Method  : 
    167       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     168      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) : 
    168169      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    169170      !!      along bottom slope gradient) an additional lateral 2nd order 
    … …  
    179180      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    180181      !!---------------------------------------------------------------------- 
    181       ! 
    182182      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    183183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    … …  
    196196      DO jn = 1, kjpt                                     ! tracer loop 
    197197         !                                                ! =========== 
    198 #  if defined key_vectopt_loop 
    199          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    200             DO ji = 1, jpij 
    201 #else 
    202198         DO jj = 1, jpj 
    203199            DO ji = 1, jpi 
    204 #endif 
    205                ik = mbkt(ji,jj)                        ! bottom T-level index 
    206                zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S 
     200               ik = mbkt(ji,jj)                              ! bottom T-level index 
     201               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S 
    207202            END DO 
    208203         END DO 
    209          !                                                ! Compute the trend 
    210 #  if defined key_vectopt_loop 
    211          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    212             DO ji = jpi+1, jpij-jpi-1 
    213 #  else 
    214          DO jj = 2, jpjm1 
     204         !                
     205         DO jj = 2, jpjm1                                    ! Compute the trend 
    215206            DO ji = 2, jpim1 
    216 #  endif 
    217                ik = mbkt(ji,jj)                            ! bottom T-level index 
     207               ik = mbkt(ji,jj)                              ! bottom T-level index 
    218208               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik) 
    219209               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
    … …  
    264254      DO jn = 1, kjpt                                            ! tracer loop 
    265255         !                                                       ! =========== 
    266 # if defined key_vectopt_loop 
    267          DO jj = 1, 1 
    268             DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling) 
    269 # else 
    270256         DO jj = 1, jpjm1 
    271257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west 
    272 # endif 
    273258               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection 
    274259                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf) 
    … …  
    333318      !!                advection terms. 
    334319      !! 
    335       !! ** Method  : 
    336       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     320      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) : 
    337321      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    338322      !!      along bottom slope gradient) an additional lateral 2nd order 
    … …  
    342326      !!      a downslope velocity of 20 cm/s if the condition for slope 
    343327      !!      convection is satified) 
    344       !!        * advective bbl (nn_bbl_adv=1 or 2) : 
     328      !!              * advective bbl (nn_bbl_adv=1 or 2) : 
    345329      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity 
    346330      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation 
    … …  
    353337      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    354338      !!---------------------------------------------------------------------- 
    355       ! 
    356339      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index 
    357       INTEGER         , INTENT(in   ) ::   kit000          ! first time step index 
     340      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
    358341      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    359342      !! 
    360343      INTEGER  ::   ji, jj                    ! dummy loop indices 
    361344      INTEGER  ::   ik                        ! local integers 
    362       INTEGER  ::   iis , iid , ijs , ijd     !   -       - 
    363       INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       - 
    364       REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars 
    365       REAL(wp) ::   zgdrho, zt, zs, zh        !   -      - 
    366       !! 
    367       REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function 
    368       REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep 
    369       !!----------------------- zv_bbl----------------------------------------------- 
    370       ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients 
    371       ! ================            pft :  potential temperature in degrees celcius 
    372       !                             pfs :  salinity anomaly (s-35) in psu 
    373       !                             pfh :  depth in meters 
    374       ! nn_eos = 0  (Jackett and McDougall 1994 formulation) 
    375       fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta 
    376          ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    & 
    377                                    - 0.203814e-03 ) * pft                    & 
    378                                    + 0.170907e-01 ) * pft                    & 
    379                                    + 0.665157e-01                            & 
    380          +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   & 
    381          +  ( ( - 0.302285e-13 * pfh                                         & 
    382                 - 0.251520e-11 * pfs                                         & 
    383                 + 0.512857e-12 * pft * pft          ) * pfh                  & 
    384                                      - 0.164759e-06   * pfs                  & 
    385              +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  & 
    386                                      + 0.380374e-04 ) * pfh 
    387       fsbeta( pft, pfs, pfh ) =                                              &   ! beta 
    388          ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      & 
    389                                  - 0.301985e-05 ) * pft                      & 
    390                                  + 0.785567e-03                              & 
    391          + (     0.515032e-08 * pfs                                          & 
    392                + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   & 
    393                +(  (   0.121551e-17 * pfh                                    & 
    394                      - 0.602281e-15 * pfs                                    & 
    395                      - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             & 
    396                                           + 0.408195e-10   * pfs             & 
    397                  + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             & 
    398                                           - 0.121555e-07 ) * pfh 
    399       !!---------------------------------------------------------------------- 
    400  
     345      INTEGER  ::   iis, iid, ikus, ikud      !   -       - 
     346      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       - 
     347      REAL(wp) ::   za, zb, zgdrho            ! local scalars 
     348      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      - 
     349      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace 
     350      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace 
     351      !!---------------------------------------------------------------------- 
    401352      ! 
    402353      IF( nn_timing == 1 )  CALL timing_start( 'bbl') 
    403354      ! 
    404       CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    405       ! 
    406  
    407355      IF( kt == kit000 )  THEN 
    408356         IF(lwp)  WRITE(numout,*) 
    … …  
    410358         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~' 
    411359      ENDIF 
    412  
    413       !                                        !* bottom temperature, salinity, velocity and depth 
    414 #if defined key_vectopt_loop 
    415       DO jj = 1, 1   ! vector opt. (forced unrolling) 
    416          DO ji = 1, jpij 
    417 #else 
     360      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity) 
    418361      DO jj = 1, jpj 
    419362         DO ji = 1, jpi 
    420 #endif 
    421             ik = mbkt(ji,jj)                        ! bottom T-level index 
    422             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S 
    423             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 
    424             zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth 
     363            ik = mbkt(ji,jj)                             ! bottom T-level index 
     364            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S 
     365            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 
    425366            ! 
    426             zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity 
    427             zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
     367            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth 
     368            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity 
     369            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
    428370         END DO 
    429371      END DO 
    430  
     372      ! 
     373      CALL eos_rab( zts, zdep, zab ) 
     374      ! 
    431375      !                                   !-------------------! 
    432376      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   ! 
    433377         !                                !-------------------! 
    434378         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 
    435             DO ji = 1, jpim1 
    436                !                                                ! i-direction 
    437                zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth 
    438                zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    439                zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    440                !                                                         ! masked bbl i-gradient of density 
    441                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    442                   &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
     379            DO ji = 1, fs_jpim1   ! vector opt. 
     380               !                                                   ! i-direction 
     381               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point 
     382               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     383               !                                                         ! 2*masked bottom density gradient 
     384               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     385                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
    443386               ! 
    444                zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
    445                ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff. 
     387               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
     388               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff. 
    446389               ! 
    447                !                                                ! j-direction 
    448                zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth 
    449                zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    450                zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    451                !                                                         ! masked bbl j-gradient of density 
    452                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    453                   &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
     390               !                                                   ! j-direction 
     391               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point 
     392               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     393               !                                                         ! 2*masked bottom density gradient 
     394               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     395                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
    454396               ! 
    455                zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
     397               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
    456398               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) 
    457                ! 
    458399            END DO 
    459400         END DO 
    … …  
    469410            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0 
    470411               DO ji = 1, fs_jpim1   ! vector opt. 
    471                   !                                               ! i-direction 
    472                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth 
    473                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    474                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    475                   !                                                           ! masked bbl i-gradient of density 
    476                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    477                      &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
    478                   ! 
    479                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope 
    480                   zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope 
    481                   ! 
    482                   !                                                           ! bbl velocity 
     412                  !                                                  ! i-direction 
     413                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     414                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     415                  !                                                          ! 2*masked bottom density gradient  
     416                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     417                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
     418                  ! 
     419                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope 
     420                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope 
     421                  ! 
     422                  !                                                          ! bbl velocity 
    483423                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) 
    484424                  ! 
    485                   !                                               ! j-direction 
    486                   zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth 
    487                   zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    488                   zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    489                   !                                                           ! masked bbl j-gradient of density 
    490                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    491                      &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
    492                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope 
    493                   zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope 
    494                   ! 
    495                   !                                                           ! bbl velocity 
     425                  !                                                  ! j-direction 
     426                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     427                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     428                  !                                                          ! 2*masked bottom density gradient 
     429                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     430                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
     431                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope 
     432                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope 
     433                  ! 
     434                  !                                                          ! bbl transport 
    496435                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) 
    497436               END DO 
    … …  
    502441            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down 
    503442               DO ji = 1, fs_jpim1   ! vector opt. 
    504                   !                                         ! i-direction 
     443                  !                                                  ! i-direction 
    505444                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf) 
    506                   iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 
    507                   ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj) 
    508                   ! 
    509                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    510                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth 
    511                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    512                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    513                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    514                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
    515                      &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1) 
    516                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    517                   ! 
    518                   !                                             ! bbl transport (down-slope direction) 
     445                  iid  = ji + MAX( 0, mgrhu(ji,jj) ) 
     446                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 
     447                  ! 
     448                  ikud = mbku_d(ji,jj) 
     449                  ikus = mbku(ji,jj) 
     450                  ! 
     451                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     452                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     453                  !                                                          !   masked bottom density gradient 
     454                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    & 
     455                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1) 
     456                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     457                  ! 
     458                  !                                                          ! bbl transport (down-slope direction) 
    519459                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) ) 
    520460                  ! 
    521                   !                                         ! j-direction 
     461                  !                                                  ! j-direction 
    522462                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf) 
    523                   ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 
    524                   ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj) 
    525                   ! 
    526                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    527                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth 
    528                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0 
    529                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 
    530                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    531                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
    532                      &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1) 
    533                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    534                   ! 
    535                   !                                             ! bbl transport (down-slope direction) 
     463                  ijd  = jj + MAX( 0, mgrhv(ji,jj) ) 
     464                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 
     465                  ! 
     466                  ikvd = mbkv_d(ji,jj) 
     467                  ikvs = mbkv(ji,jj) 
     468                  ! 
     469                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     470                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     471                  !                                                          !   masked bottom density gradient 
     472                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    & 
     473                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1) 
     474                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     475                  ! 
     476                  !                                                          ! bbl transport (down-slope direction) 
    536477                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) ) 
    537478               END DO 
    … …  
    541482      ENDIF 
    542483      ! 
    543       CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    544       ! 
    545484      IF( nn_timing == 1 )  CALL timing_stop( 'bbl') 
    546485      ! 
    … …  
    558497      !!---------------------------------------------------------------------- 
    559498      INTEGER ::   ji, jj               ! dummy loop indices 
    560       INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer 
    561       INTEGER  ::   ios                 ! Local integer output status for namelist read 
     499      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer 
     500      INTEGER ::   ios                  !   -      - 
    562501      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
    563502      !! 
    … …  
    598537      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)' 
    599538 
    600       IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' ) 
    601  
    602539      !                             !* vertical index of  "deep" bottom u- and v-points 
    603540      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv) 
    … …  
    607544         END DO 
    608545      END DO 
    609       ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
     546      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
    610547      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    611548      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    612549 
    613                                      !* sign of grad(H) at u- and v-points 
    614       mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0. 
     550                                        !* sign of grad(H) at u- and v-points 
     551      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
    615552      DO jj = 1, jpjm1 
    616553         DO ji = 1, jpim1 
    … …  
    621558 
    622559      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point 
    623          DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0) 
     560         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0) 
    624561            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) ) 
    625562            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r4624 r4933  
    2828   USE dom_oce        ! ocean: domain variables 
    2929   USE c1d            ! 1D vertical configuration 
    30    USE trdmod_oce     ! ocean: trend variables 
    31    USE trdtra         ! active tracers: trends 
     30   USE trd_oce        ! trends: ocean variables 
     31   USE trdtra         ! trends manager: tracers  
    3232   USE zdf_oce        ! ocean: vertical physics 
    3333   USE phycst         ! physical constants 
    … …  
    4848   PUBLIC   dtacof_zoom  ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90 
    4949 
     50!!gm  why all namelist variable public????   only ln_tradmp should be sufficient 
     51 
    5052   !                               !!* Namelist namtra_dmp : T & S newtonian damping * 
    5153   LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag 
    … …  
    112114      ! 
    113115      CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta ) 
     116      ! 
    114117      !                           !==   input T-S data at kt   ==! 
    115118      CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt 
    … …  
    172175      ! 
    173176      IF( l_trdtra )   THEN       ! trend diagnostic 
    174          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ttrdmp ) 
    175          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, strdmp ) 
     177         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 
     178         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 
    176179      ENDIF 
    177180      !                           ! Control print 
    … …  
    194197      !! ** Method  :   read the namtra_dmp namelist and check the parameters 
    195198      !!---------------------------------------------------------------------- 
     199      INTEGER  ::   ios   ! Local integer output status for namelist read 
     200      !! 
    196201      NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file 
    197       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    198       !!---------------------------------------------------------------------- 
    199  
     202      !!---------------------------------------------------------------------- 
     203      ! 
    200204      REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 
    201205      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    202206901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 
    203  
     207      ! 
    204208      REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term 
    205209      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    … …  
    228232         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 
    229233         ! 
     234!!gm  I don't understand the specificities of c1d case......    
     235!!gm  to be check with the autor of these lines 
     236          
    230237#if ! defined key_c1d 
    231238         SELECT CASE ( nn_hdmp ) 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4488 r4933  
    2323   USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine) 
    2424   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine) 
    25    USE trdmod_oce      ! ocean space and time domain 
    26    USE trdtra          ! ocean active tracers trends 
     25   USE trd_oce         ! trends: ocean variables 
     26   USE trdtra          ! trends manager: tracers  
     27   ! 
    2728   USE prtctl          ! Print control 
    2829   USE in_out_manager  ! I/O manager 
    … …  
    3536   PRIVATE 
    3637 
    37    PUBLIC   tra_ldf         ! called by step.F90  
    38    PUBLIC   tra_ldf_init    ! called by opa.F90  
     38   PUBLIC   tra_ldf        ! called by step.F90  
     39   PUBLIC   tra_ldf_init   ! called by opa.F90  
    3940   ! 
    4041   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... namlist logicals) 
    … …  
    112113         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    113114         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    114          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_ldf, ztrdt ) 
    115          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_ldf, ztrds ) 
     115         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
     116         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 
    116117         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    117118      ENDIF 
    … …  
    174175            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    175176         ENDIF 
    176          IF ( ln_zps ) THEN             ! z-coordinate 
     177         IF ( ln_zps ) THEN             ! zps-coordinate 
    177178            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed 
    178179            IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation) 
    179180            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    180181         ENDIF 
    181          IF ( ln_sco ) THEN             ! z-coordinate 
     182         IF ( ln_sco ) THEN             ! s-coordinate 
    182183            IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation) 
    183184            IF ( ln_traldf_hor   )   nldf = 1      ! horizontal (   rotation) 
    … …  
    192193            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    193194         ENDIF 
    194          IF ( ln_zps ) THEN             ! z-coordinate 
     195         IF ( ln_zps ) THEN             ! zps-coordinate 
    195196            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed  
    196197            IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation) 
    197198            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    198199         ENDIF 
    199          IF ( ln_sco ) THEN             ! z-coordinate 
     200         IF ( ln_sco ) THEN             ! s-coordinate 
    200201            IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation) 
    201202            IF ( ln_traldf_hor   )   nldf = 3      ! horizontal (   rotation) 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r3632 r4933  
    252252         END DO 
    253253         IF( ln_zps.and.l_grad_zps ) THEN              ! partial steps: correction at the last level 
    254 # if defined key_vectopt_loop 
    255             DO jj = 1, 1 
    256                DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    257 # else 
    258254            DO jj = 1, jpjm1 
    259255               DO ji = 1, jpim1 
    260 # endif 
    261256                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    262257                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r4313 r4933  
    22   !!============================================================================== 
    33   !!                       ***  MODULE  tranpc  *** 
    4    !! Ocean active tracers:  non penetrative convection scheme 
     4   !! Ocean active tracers:  non penetrative convective adjustment scheme 
    55   !!============================================================================== 
    66   !! History :  1.0  ! 1990-09  (G. Madec)  Original code 
    … …  
    99   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90 
    1010   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
     11   !!            3.7  ! 2014-06  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 
    1112   !!---------------------------------------------------------------------- 
    1213 
    … …  
    1415   !!   tra_npc : apply the non penetrative convection scheme 
    1516   !!---------------------------------------------------------------------- 
    16    USE oce             ! ocean dynamics and active tracers  
     17   USE oce             ! ocean dynamics and active tracers 
    1718   USE dom_oce         ! ocean space and time domain 
     19   USE phycst          ! physical constants 
    1820   USE zdf_oce         ! ocean vertical physics 
    19    USE trdmod_oce      ! ocean active tracer trends 
     21   USE trd_oce         ! ocean active tracer trends 
    2022   USE trdtra          ! ocean active tracer trends 
    21    USE eosbn2          ! equation of state (eos routine)  
     23   USE eosbn2          ! equation of state (eos routine) 
     24   ! 
    2225   USE lbclnk          ! lateral boundary conditions (or mpp link) 
    2326   USE in_out_manager  ! I/O manager 
    … …  
    2932   PRIVATE 
    3033 
    31    PUBLIC   tra_npc       ! routine called by step.F90 
     34   PUBLIC   tra_npc    ! routine called by step.F90 
    3235 
    3336   !! * Substitutions 
    3437#  include "domzgr_substitute.h90" 
    35    !!---------------------------------------------------------------------- 
    36    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    37    !! $Id$  
     38#  include "vectopt_loop_substitute.h90" 
     39   !!---------------------------------------------------------------------- 
     40   !! NEMO/OPA 3.6 , NEMO Consortium (2014) 
     41   !! $Id$ 
    3842   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3943   !!---------------------------------------------------------------------- 
    … …  
    4448      !!                  ***  ROUTINE tranpc  *** 
    4549      !! 
    46       !! ** Purpose :   Non penetrative convective adjustment scheme. solve  
     50      !! ** Purpose : Non-penetrative convective adjustment scheme. solve 
    4751      !!      the static instability of the water column on after fields 
    4852      !!      while conserving heat and salt contents. 
    4953      !! 
    50       !! ** Method  :   The algorithm used converges in a maximium of jpk  
    51       !!      iterations. instabilities are treated when the vertical density 
    52       !!      gradient is less than 1.e-5. 
    53       !!      l_trdtra=T: the trend associated with this algorithm is saved. 
     54      !! ** Method  : updated algorithm able to deal with non-linear equation of state 
     55      !!              (i.e. static stability computed locally) 
    5456      !! 
    5557      !! ** Action  : - (ta,sa) after the application od the npc scheme 
    56       !!              - save the associated trends (ttrd,strd) ('key_trdtra') 
     58      !!              - send the associated trends for on-line diagnostics (l_trdtra=T) 
    5759      !! 
    58       !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
     60      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
    5961      !!---------------------------------------------------------------------- 
    60       ! 
    6162      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    6263      ! 
    6364      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6465      INTEGER  ::   inpcc        ! number of statically instable water column 
    65       INTEGER  ::   inpci        ! number of iteration for npc scheme 
    66       INTEGER  ::   jiter, jkdown, jkp        ! ??? 
    67       INTEGER  ::   ikbot, ik, ikup, ikdown   ! ??? 
    68       REAL(wp) ::   ze3tot, zta, zsa, zraua, ze3dwn 
    69       REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx, zwy, zwz 
    70       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds, zrhop 
     66      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers 
     67      LOGICAL  ::   l_bottom_reached, l_column_treated 
     68      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
     69      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     70      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
     71      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
     72      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta 
     73      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2  
     74      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta 
     75      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
     76      ! 
     77      !!LB debug: 
     78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 
     79      INTEGER :: ilc1, jlc1, klc1, nncpu 
     80      LOGICAL :: lp_monitor_point = .FALSE. 
     81      !!LB debug. 
    7182      !!---------------------------------------------------------------------- 
    7283      ! 
    7384      IF( nn_timing == 1 )  CALL timing_start('tra_npc') 
    7485      ! 
    75       CALL wrk_alloc(jpi, jpj, jpk, zrhop ) 
    76       CALL wrk_alloc(jpi, jpk, zwx, zwy, zwz ) 
    77       ! 
    7886      IF( MOD( kt, nn_npc ) == 0 ) THEN 
    79  
    80          inpcc = 0 
    81          inpci = 0 
    82  
    83          CALL eos( tsa, rhd, zrhop, fsdept_n(:,:,:) )         ! Potential density 
    84  
    85          IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     87         ! 
     88         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2 
     89         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta 
     90         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj 
     91         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj 
     92 
     93         IF( l_trdtra )   THEN                    !* Save initial after fields 
    8694            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    8795            ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
    … …  
    8997         ENDIF 
    9098 
    91          !                                                ! =============== 
    92          DO jj = 1, jpj                                   !  Vertical slab 
    93             !                                             ! =============== 
    94             !  Static instability pointer  
    95             ! ---------------------------- 
    96             DO jk = 1, jpkm1 
    97                DO ji = 1, jpi 
    98                   zwx(ji,jk) = ( zrhop(ji,jj,jk) - zrhop(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 
    99                END DO 
    100             END DO 
    101  
    102             ! 1.1 do not consider the boundary points 
    103  
    104             ! even if east-west cyclic b. c. do not considere ji=1 or jpi 
    105             DO jk = 1, jpkm1 
    106                zwx( 1 ,jk) = 0.e0 
    107                zwx(jpi,jk) = 0.e0 
    108             END DO 
    109             ! even if south-symmetric b. c. used, do not considere jj=1 
    110             IF( jj == 1 )   zwx(:,:) = 0.e0 
    111  
    112             DO jk = 1, jpkm1 
    113                DO ji = 1, jpi 
    114                   zwx(ji,jk) = 1. 
    115                   IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0 
    116                END DO 
    117             END DO 
    118  
    119             zwy(:,1) = 0.e0 
    120             DO ji = 1, jpi 
    121                DO jk = 1, jpkm1 
    122                   zwy(ji,1) = zwy(ji,1) + zwx(ji,jk) 
    123                END DO 
    124             END DO 
    125  
    126             zwz(1,1) = 0.e0 
    127             DO ji = 1, jpi 
    128                zwz(1,1) = zwz(1,1) + zwy(ji,1) 
    129             END DO 
    130  
    131             inpcc = inpcc + NINT( zwz(1,1) ) 
    132  
    133  
    134             ! 2. Vertical mixing for each instable portion of the density profil 
    135             ! ------------------------------------------------------------------ 
    136  
    137             IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable : 
    138                DO ji = 1, jpi 
    139                   IF( zwy(ji,1) /= 0.e0 ) THEN 
     99         !LB debug: 
     100         IF( lwp .AND. l_LB_debug ) THEN 
     101            WRITE(numout,*) 
     102            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea 
     103         ENDIF 
     104         !LBdebug: Monitoring of 1 column subject to convection... 
     105         IF( l_LB_debug ) THEN 
     106            ! Location of 1 known convection spot to follow what's happening in the water column 
     107            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 
     108            nncpu = 15  ; ! the CPU domain contains the convection spot 
     109            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 
     110            !nncpu = 54  ; ! the CPU domain contains the convection spot 
     111            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
     112         ENDIF 
     113         !LBdebug. 
     114 
     115         CALL eos_rab( tsa, zab )         ! after alpha and beta 
     116         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala 
     117         
     118         inpcc = 0 
     119 
     120         DO jj = 2, jpjm1                 ! interior column only 
     121            DO ji = fs_2, fs_jpim1 
     122               ! 
     123               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
     124                  !                                     ! consider one ocean column  
     125                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature 
     126                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity 
     127 
     128                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha  
     129                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta   
     130                  zvn2(:)         = zn2(ji,jj,:)            ! N^2  
     131                  
     132                  IF( l_LB_debug ) THEN                  !LB debug: 
     133                     lp_monitor_point = .FALSE. 
     134                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
     135                     ! writing only if on CPU domain where conv region is: 
     136                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point  
     137                      
     138                     IF(lp_monitor_point) THEN 
     139                        WRITE(numout,*) '' ;WRITE(numout,*) '' ; 
     140                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt 
     141                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') , ji,jj 
     142                        DO jk = 1, klc1 
     143                           WRITE(numout,*) jk, zvn2(jk) 
     144                        END DO 
     145                        WRITE(numout,*) ' ' 
     146                     ENDIF 
     147                  ENDIF                                  !LB debug  end 
     148 
     149                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
     150                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2) 
     151                  ilayer = 0 
     152                  jiter  = 0 
     153                  l_column_treated = .FALSE. 
     154                  
     155                  DO WHILE ( .NOT. l_column_treated ) 
    140156                     ! 
    141                      ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level 
     157                     jiter = jiter + 1 
     158                     
     159                     IF( jiter >= 400 ) EXIT 
     160                     
     161                     l_bottom_reached = .FALSE. 
     162 
     163                     DO WHILE ( .NOT. l_bottom_reached ) 
     164 
     165                        ik = ik + 1 
     166                        
     167                        !! Checking level ik for instability 
     168                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     169 
     170                        IF( zvn2(ik) < 0. ) THEN ! Instability found! 
     171 
     172                           ikm  = ik              ! first level whith negative N2 
     173                           ilayer = ilayer + 1    ! yet another layer found.... 
     174                           IF(jiter == 1) inpcc = inpcc + 1 
     175 
     176                           IF(l_LB_debug .AND. lp_monitor_point) & 
     177                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 
     178                              & ' inpcc =', inpcc 
     179 
     180                           !! Case we mix with upper regions where N2==0: 
     181                           !! All the points above ikup where N2 == 0 must also be mixed => we go 
     182                           !! upward to find a new ikup, where the layer doesn't have N2==0 
     183                           ikup = ikm 
     184                           DO jk = ikm, 2, -1 
     185                              ikup = ikup - 1 
     186                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 
     187                           END DO 
     188                           
     189                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 
     190                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 
     191 
     192                           
     193                           IF( lp_monitor_point )   WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 
     194                           
     195                           zsum_temp = 0._wp 
     196                           zsum_sali = 0._wp 
     197                           zsum_alfa = 0._wp 
     198                           zsum_beta = 0._wp 
     199                           zsum_z    = 0._wp 
     200                                                     
     201                           DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column 
     202                              ! 
     203                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk 
     204                              ! 
     205                              zdz       = fse3t(ji,jj,jk) 
     206                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz 
     207                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz 
     208                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz 
     209                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
     210                              zsum_z    = zsum_z    + zdz 
     211                              ! 
     212                              !! EXIT if we found the bottom of the unstable portion of the water column     
     213                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT 
     214                           END DO 
     215                           
     216                           !ik     = jk !LB remove? 
     217                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 
     218                           
     219                           IF(l_LB_debug .AND. lp_monitor_point) & 
     220                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer 
     221                           
     222                           ! Mixing Temperature and salinity between ikup and ikdown: 
     223                           zta   = zsum_temp/zsum_z 
     224                           zsa   = zsum_sali/zsum_z 
     225                           zalfa = zsum_alfa/zsum_z 
     226                           zbeta = zsum_beta/zsum_z 
     227 
     228                           IF(l_LB_debug .AND. lp_monitor_point) THEN 
     229                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
     230                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
     231                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa 
     232                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
     233                           ENDIF 
     234 
     235                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column 
     236                           DO jk = ikup, ikdown 
     237                              zvts(jk,jp_tem) = zta 
     238                              zvts(jk,jp_sal) = zsa 
     239                              zvab(jk,jp_tem) = zalfa 
     240                              zvab(jk,jp_sal) = zbeta 
     241                           END DO 
     242                           ! 
     243                           !! Before updating N2, it is possible that another unstable 
     244                           !! layer exists underneath the one we just homogeneized! 
     245                           ik = ikdown 
     246                           !  
     247                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN 
     248                        ! 
     249                        IF( ik == ikbot ) l_bottom_reached = .TRUE. 
     250                        ! 
     251                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
     252 
     253                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1' 
     254                     
     255                     ! ******* At this stage ik == ikbot ! ******* 
     256                     
     257                     IF( ilayer > 0 ) THEN 
     258                        !! least an unstable layer has been found 
     259                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     260                        !! => Need to re-compute N2! will use Alpha and Beta! 
     261                        ! 
     262                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!      
     263                           !! Doing exactly as in eosbn2.F90: 
     264                           !! * Except that we only are interested in the sign of N2 !!! 
     265                           !!   => just considering the vertical gradient of density 
     266                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     267                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     268                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     269                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     270                           
     271                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     272                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     273                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     274                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     275                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )   
     276                        END DO 
     277 
     278                        IF(l_LB_debug .AND. lp_monitor_point) THEN 
     279                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 
     280                              & jiter, ji,jj 
     281                           DO jk = 1, klc1 
     282                              WRITE(numout,*) jk, zvn2(jk) 
     283                           END DO 
     284                           WRITE(numout,*) ' ' 
     285                        ENDIF 
     286 
     287                        ik     = 1  ! starting again at the surface for the next iteration 
     288                        ilayer = 0 
     289                     ENDIF 
    142290                     ! 
    143                      DO jiter = 1, jpk          ! vertical iteration 
    144                         ! 
    145                         ! search of ikup : the first static instability from the sea surface 
    146                         ! 
    147                         ik = 0 
    148 220                     CONTINUE 
    149                         ik = ik + 1 
    150                         IF( ik >= ikbot ) GO TO 200 
    151                         zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1) 
    152                         IF( zwx(ji,ik) <= 0.e0 ) GO TO 220 
    153                         ikup = ik 
    154                         ! the density profil is instable below ikup 
    155                         ! ikdown : bottom of the instable portion of the density profil 
    156                         ! search of ikdown and vertical mixing from ikup to ikdown 
    157                         ! 
    158                         ze3tot= fse3t(ji,jj,ikup) 
    159                         zta   = tsa  (ji,jj,ikup,jp_tem) 
    160                         zsa   = tsa  (ji,jj,ikup,jp_sal) 
    161                         zraua = zrhop(ji,jj,ikup) 
    162                         ! 
    163                         DO jkdown = ikup+1, ikbot-1 
    164                            IF( zraua <= zrhop(ji,jj,jkdown) ) THEN 
    165                               ikdown = jkdown 
    166                               GO TO 240 
    167                            ENDIF 
    168                            ze3dwn =  fse3t(ji,jj,jkdown) 
    169                            ze3tot =  ze3tot + ze3dwn 
    170                            zta   = ( zta*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_tem)*ze3dwn )/ze3tot 
    171                            zsa   = ( zsa*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_sal)*ze3dwn )/ze3tot 
    172                            zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot 
    173                            inpci = inpci+1 
    174                         END DO 
    175                         ikdown = ikbot-1 
    176 240                     CONTINUE 
    177                         ! 
    178                         DO jkp = ikup, ikdown-1 
    179                            tsa  (ji,jj,jkp,jp_tem) = zta 
    180                            tsa  (ji,jj,jkp,jp_sal) = zsa 
    181                            zrhop(ji,jj,jkp       ) = zraua 
    182                         END DO 
    183                         IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN 
    184                            tsa  (ji,jj,jkp,jp_tem) = zta 
    185                            tsa  (ji,jj,jkp,jp_sal) = zsa 
    186                            zrhop(ji,jj,ikdown    ) = zraua 
    187                         ENDIF 
    188                      END DO 
    189                   ENDIF 
    190 200               CONTINUE 
    191                END DO 
    192                ! <<-- no more static instability on slab jj 
    193             ENDIF 
    194             !                                             ! =============== 
    195          END DO                                           !   End of slab 
    196          !                                                ! =============== 
    197          !  
    198          IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic 
    199             ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    200             ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    201             CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_npc, ztrdt ) 
    202             CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_npc, ztrds ) 
     291                     IF( ik >= ikbot ) THEN 
     292                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---' 
     293                        l_column_treated = .TRUE. 
     294                     ENDIF 
     295                     ! 
     296                  END DO ! DO WHILE ( .NOT. l_column_treated ) 
     297 
     298                  !! Updating tsa: 
     299                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem) 
     300                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 
     301 
     302                  !! lolo:  Should we update something else???? 
     303                  !! => like alpha and beta? 
     304 
     305                  IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '' 
     306 
     307               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN 
     308 
     309            END DO ! ji 
     310         END DO ! jj 
     311         ! 
     312         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
     313            z1_r2dt = 1._wp / (2._wp * rdt) 
     314            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt 
     315            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt 
     316            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
     317            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 
    203318            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    204319         ENDIF 
    205        
    206          ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign) 
    207          ! ------------------------------============ 
     320         ! 
    208321         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    209        
    210  
    211          !  2. non penetrative convective scheme statistics 
    212          !  ----------------------------------------------- 
    213          IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN 
    214             IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   & 
    215                &                   ' water column : ',inpcc, ' number of iteration : ',inpci 
    216          ENDIF 
    217          ! 
    218       ENDIF 
    219       ! 
    220       CALL wrk_dealloc(jpi, jpj, jpk, zrhop ) 
    221       CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz ) 
     322         ! 
     323         IF(lwp) THEN 
     324            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 
     325            WRITE(numout,*)' => number of statically instable water column : ',inpcc 
     326            WRITE(numout,*) '' ; WRITE(numout,*) '' 
     327         ENDIF 
     328         ! 
     329         CALL wrk_dealloc(jpi, jpj, jpk, zn2 ) 
     330         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab ) 
     331         CALL wrk_dealloc(jpk, zvn2 ) 
     332         CALL wrk_dealloc(jpk, 2, zvts, zvab ) 
     333         ! 
     334      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN 
    222335      ! 
    223336      IF( nn_timing == 1 )  CALL timing_stop('tra_npc') 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r4328 r4933  
    2727   USE dom_oce         ! ocean space and time domain variables  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    29    USE zdf_oce         ! ??? 
     29   USE zdf_oce         ! ocean vertical mixing 
    3030   USE domvvl          ! variable volume 
    3131   USE dynspg_oce      ! surface     pressure gradient variables 
    3232   USE dynhpg          ! hydrostatic pressure gradient  
    33    USE trdmod_oce      ! ocean space and time domain variables  
    34    USE trdtra          ! ocean active tracers trends  
    35    USE phycst 
    36    USE bdy_oce 
     33   USE trd_oce         ! trends: ocean variables 
     34   USE trdtra          ! trends manager: tracers  
     35   USE traqsr          ! penetrative solar radiation (needed for nksr) 
     36   USE phycst          ! physical constant 
     37   USE ldftra_oce      ! lateral physics on tracers 
     38   USE bdy_oce         ! BDY open boundary condition variables 
    3739   USE bdytra          ! open boundary condition (bdy_tra routine) 
     40   ! 
    3841   USE in_out_manager  ! I/O manager 
    3942   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4043   USE prtctl          ! Print control 
    41    USE traqsr          ! penetrative solar radiation (needed for nksr) 
     44   USE wrk_nemo        ! Memory allocation 
     45   USE timing          ! Timing 
    4246#if defined key_agrif 
    4347   USE agrif_opa_update 
    4448   USE agrif_opa_interp 
    4549#endif 
    46    USE wrk_nemo        ! Memory allocation 
    47    USE timing          ! Timing 
    4850 
    4951   IMPLICIT NONE 
    … …  
    8082      !!             at the local domain   boundaries through lbc_lnk call,  
    8183      !!             at the one-way open boundaries (lk_bdy=T),  
    82       !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
     84      !!             at the AGRIF zoom   boundaries (lk_agrif=T) 
    8385      !! 
    8486      !!              - Update lateral boundary conditions on AGRIF children 
    … …  
    127129         ztrdt(:,:,:) = tsn(:,:,:,jp_tem)  
    128130         ztrds(:,:,:) = tsn(:,:,:,jp_sal) 
     131         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
     132            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
     133            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
     134         ENDIF 
    129135      ENDIF 
    130136 
    … …  
    150156      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt      
    151157         DO jk = 1, jpkm1 
    152             zfact = 1.e0_wp / r2dtra(jk)              
     158            zfact = 1._wp / r2dtra(jk)              
    153159            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 
    154160            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 
    155161         END DO 
    156          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt ) 
    157          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds ) 
     162         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 
     163         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 
    158164         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    159165      END IF 
    … …  
    163169         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask ) 
    164170      ! 
    165       ! 
    166       IF( nn_timing == 1 )  CALL timing_stop('tra_nxt') 
     171      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt') 
    167172      ! 
    168173   END SUBROUTINE tra_nxt 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4834 r4933  
    2121   USE sbc_oce         ! surface boundary condition: ocean 
    2222   USE trc_oce         ! share SMS/Ocean variables 
    23    USE trdmod_oce      ! ocean variables trends 
    24    USE trdtra          ! ocean active tracers trends  
     23   USE trd_oce        ! trends: ocean variables 
     24   USE trdtra         ! trends manager: tracers 
    2525   USE in_out_manager  ! I/O manager 
    2626   USE phycst          ! physical constants 
    … …  
    169169               DO ji = 1, jpi 
    170170                  IF ( qsr(ji,jj) /= 0._wp ) THEN 
    171                      oatte(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 
    172                      iatte(ji,jj) = oatte(ji,jj) 
     171                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 
    173172                  ENDIF 
    174173               END DO 
    … …  
    241240                        zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
    242241                        zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
    243                         oatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2)  
    244                         iatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zcoef + zcoef ) * tmask(ji,jj,2) 
     242                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2)  
    245243                     END DO 
    246244                  END DO 
    … …  
    259257               ! clem: store attenuation coefficient of the first ocean level 
    260258               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    261                   oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    262                   iatte(:,:) = oatte(:,:) 
     259                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    263260               ENDIF 
    264261           ENDIF 
    … …  
    287284                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
    288285                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
    289                         oatte(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 
    290                         iatte(ji,jj) = oatte(ji,jj) 
     286                        fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 
    291287                     END DO 
    292288                  END DO 
    … …  
    302298               ! clem: store attenuation coefficient of the first ocean level 
    303299               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    304                   oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    305                   iatte(:,:) = oatte(:,:) 
     300                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    306301               ENDIF 
    307302               ! 
    … …  
    334329      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    335330         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    336          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt ) 
     331         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    337332         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )  
    338333      ENDIF 
    … …  
    384379      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init') 
    385380      ! 
    386       ! clem init for oatte and iatte 
     381      ! Default value for fraqsr_1lev 
    387382      IF( .NOT. ln_rstart ) THEN 
    388          oatte(:,:) = 1._wp 
    389          iatte(:,:) = 1._wp 
     383         fraqsr_1lev(:,:) = 1._wp 
    390384      ENDIF 
    391385      ! 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r3764 r4933  
    1818   USE dom_oce         ! ocean space domain variables 
    1919   USE phycst          ! physical constant 
     20   USE sbcmod          ! ln_rnf   
     21   USE sbcrnf          ! River runoff   
    2022   USE traqsr          ! solar radiation penetration 
    21    USE trdmod_oce      ! ocean trends  
    22    USE trdtra          ! ocean trends 
     23   USE trd_oce         ! trends: ocean variables 
     24   USE trdtra          ! trends manager: tracers  
     25   ! 
    2326   USE in_out_manager  ! I/O manager 
    2427   USE prtctl          ! Print control 
    25    USE sbcrnf          ! River runoff   
    26    USE sbcmod          ! ln_rnf   
    27    USE iom 
     28   USE iom             ! I/O library 
    2829   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2930   USE wrk_nemo        ! Memory Allocation 
    … …  
    3940#  include "vectopt_loop_substitute.h90" 
    4041   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     42   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    4243   !! $Id$ 
    4344   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    9192      !!         where emp, the surface freshwater budget (evaporation minus 
    9293      !!         precipitation minus runoff) given in kg/m2/s is divided 
    93       !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.     
     94      !!         by rau0 (density of sea water) to obtain m/s.     
    9495      !!         Note: even though Fwe does not appear explicitly for  
    9596      !!         temperature in this routine, the heat carried by the water 
    … …  
    107108      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated 
    108109      !!                with the tracer surface boundary condition  
    109       !!              - save the trend it in ttrd ('key_trdtra') 
     110      !!              - send trends to trdtra module (l_trdtra=T) 
    110111      !!---------------------------------------------------------------------- 
    111112      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    … …  
    124125      ENDIF 
    125126 
    126       IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     127      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    127128         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )  
    128129         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    … …  
    137138 
    138139      !---------------------------------------- 
    139       !        EMP, EMPS and QNS effects 
     140      !        EMP, SFX and QNS effects 
    140141      !---------------------------------------- 
    141142      !                                          Set before sbc tracer content fields 
    … …  
    146147              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
    147148            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file' 
    148             zfact = 0.5e0 
     149            zfact = 0.5_wp 
    149150            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend 
    150151            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend 
    151152         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
    152             zfact = 1.e0 
    153             sbc_tsc_b(:,:,:) = 0.e0 
     153            zfact = 1._wp 
     154            sbc_tsc_b(:,:,:) = 0._wp 
    154155         ENDIF 
    155156      ELSE                                         ! Swap of forcing fields 
    156157         !                                         ! ---------------------- 
    157          zfact = 0.5e0 
     158         zfact = 0.5_wp 
    158159         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 
    159160      ENDIF 
    … …  
    226227      ENDIF 
    227228  
    228       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
     229      IF( l_trdtra )   THEN                      ! send trends for further diagnostics 
    229230         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    230231         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    231          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_nsr, ztrdt ) 
    232          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_nsr, ztrds ) 
     232         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
     233         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 
    233234         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    234235      ENDIF 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r3294 r4933  
    1919   USE sbc_oce         ! surface boundary condition: ocean 
    2020   USE dynspg_oce 
    21  
    2221   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine) 
    2322   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine) 
    24  
    2523   USE ldftra_oce      ! ocean active tracers: lateral physics 
    26    USE trdmod_oce      ! ocean active tracers: lateral physics 
    27    USE trdtra      ! ocean tracers trends  
     24   USE trd_oce         ! trends: ocean variables 
     25   USE trdtra          ! trends manager: tracers  
     26   ! 
    2827   USE in_out_manager  ! I/O manager 
    2928   USE prtctl          ! Print control 
    … …  
    3231   USE wrk_nemo        ! Memory allocation 
    3332   USE timing          ! Timing 
    34  
    3533 
    3634   IMPLICIT NONE 
    … …  
    4745#  include "vectopt_loop_substitute.h90" 
    4846   !!---------------------------------------------------------------------- 
    49    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     47   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    5048   !! $Id$ 
    5149   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    9694            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 
    9795         END DO 
    98          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt ) 
    99          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds ) 
     96         CALL lbc_lnk( ztrdt, 'T', 1. ) 
     97         CALL lbc_lnk( ztrds, 'T', 1. ) 
     98         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     99         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    100100         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    101101      ENDIF 

  • branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r3294 r4933  
    7474      !!          Idem for di(s) and dj(s)           
    7575      !! 
    76       !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at  
    77       !!      the good depth zh from interpolated T and S for the different 
    78       !!      formulation of the equation of state (eos). 
     76      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     77      !!      depth zh from interpolated T and S for the different formulations 
     78      !!      of the equation of state (eos). 
    7979      !!      Gradient formulation for rho : 
    80       !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 
     80      !!          di(rho) = rd~ - rd(i,j,k)   or  rd(i+1,j,k) - rd~ 
    8181      !! 
    8282      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
    8383      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points  
    8484      !!---------------------------------------------------------------------- 
    85       ! 
    8685      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    8786      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
    8887      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    8988      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    90       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    91       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
     89      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
     90      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
    9291      ! 
    9392      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    9493      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    9594      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    96       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    97       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj    ! interpolated value of tracer 
     95      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     96      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
    9897      !!---------------------------------------------------------------------- 
    9998      ! 
    10099      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
    101100      ! 
    102       CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    103       CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           )  
    104       ! 
    105101      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    106102         ! 
    107 # if defined key_vectopt_loop 
    108          jj = 1 
    109          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    110 # else 
    111103         DO jj = 1, jpjm1 
    112104            DO ji = 1, jpim1 
    113 # endif 
    114105               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    115106               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    … …  
    121112                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    122113                  ! interpolated values of tracers 
    123                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     114                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    124115                  ! gradient of  tracers 
    125116                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    … …  
    127118                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    128119                  ! interpolated values of tracers 
    129                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     120                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    130121                  ! gradient of tracers 
    131122                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    … …  
    136127                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    137128                  ! interpolated values of tracers 
    138                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     129                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    139130                  ! gradient of tracers 
    140131                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    … …  
    142133                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    143134                  ! interpolated values of tracers 
    144                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     135                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    145136                  ! gradient of tracers 
    146137                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    147138               ENDIF 
    148 # if ! defined key_vectopt_loop 
    149139            END DO 
    150 # endif 
    151140         END DO 
    152141         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    … …  
    156145      ! horizontal derivative of density anomalies (rd) 
    157146      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    158 # if defined key_vectopt_loop 
    159          jj = 1 
    160          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    161 # else 
    162147         DO jj = 1, jpjm1 
    163148            DO ji = 1, jpim1 
    164 # endif 
    165149               iku = mbku(ji,jj) 
    166150               ikv = mbkv(ji,jj) 
    … …  
    173157               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
    174158               ENDIF 
    175 # if ! defined key_vectopt_loop 
    176159            END DO 
    177 # endif 
    178160         END DO 
    179161 
    … …  
    184166 
    185167         ! Gradient of density at the last level  
    186 # if defined key_vectopt_loop 
    187          jj = 1 
    188          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    189 # else 
    190168         DO jj = 1, jpjm1 
    191169            DO ji = 1, jpim1 
    192 # endif 
    193170               iku = mbku(ji,jj) 
    194171               ikv = mbkv(ji,jj) 
    195172               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    196173               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    197                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    198                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
     174               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     175               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
    199176               ENDIF 
    200                IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    201                ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
     177               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     178               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
    202179               ENDIF 
    203 # if ! defined key_vectopt_loop 
    204180            END DO 
    205 # endif 
    206181         END DO 
    207182         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
    208183         ! 
    209184      END IF 
    210       ! 
    211       CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    212       CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj           )  
    213185      ! 
    214186      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
Note: See TracChangeset for help on using the changeset viewer.