Changeset 3893


Ignore:
Timestamp:
2013-04-24T16:00:59+02:00 (7 years ago)
Author:
gm
Message:

dev_r3858_CNRS3_Ediag: #927 keep only Vallis and UNESCO eos introduce alpha beta and their primitive

Location:
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3876 r3893  
    564564&nameos        !   ocean physical parameters 
    565565!----------------------------------------------------------------------- 
    566    nn_eos      =   0       !  type of equation of state and Brunt-Vaisala frequency 
    567                            !     = 0, UNESCO (formulation of Jackett and McDougall (1994) and of McDougall (1987) ) 
    568                            !     = 1, linear: rho(T)   = rau0 * ( 1.028 - ralpha * T ) 
    569                            !     = 2, linear: rho(T,S) = rau0 * ( rbeta * S - ralpha * T ) 
    570    rn_alpha    =   2.0e-4  !  thermal expension coefficient (nn_eos= 1 or 2) 
    571    rn_beta     =   7.7e-4  !  saline  expension coefficient (nn_eos= 2) 
     566   nn_eos      =  0        !  type of equation of state and Brunt-Vaisala frequency 
     567                           !     =-1, UNESCO (formulation of Jackett and McDougall (1995) ) 
     568                           !     = 0, UNESCO (modified formulation of Jackett and McDougall (1994) and of McDougall (1987) ) 
     569                           !     = 1, simplified eos (based on Vallis, 2006):  
     570                           !             prd(T,S,z)  =  - rn_alpha * (1+rn_cabbel*(T-rn_temp0)+rn_thermo*z) * (T-rn_temp0)  
     571                           !                            + rn_beta  * (S-rn_salt0)  
     572                           !                            + rn_gamma * z 
     573   rn_alph0    =   2.0e-4  !  thermal expension coefficient (nn_eos= 1) 
     574   rn_beta0    =   7.7e-4  !  saline  expension coefficient (nn_eos= 1) 
     575   rn_gamm0    =   4.39e-6 !  depth expansion coeff. (=0 for linear eos) 
     576   rn_cabbel   =   2.5e-2  !  cabbeling coeff. (=0 for linear eos) 
     577   rn_thermo   =   1.1e-4  !  thermobaric coeff. (=0 for linear eos) 
    572578/ 
    573579!----------------------------------------------------------------------- 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r3764 r3893  
    3232   USE phycst          ! physical constants 
    3333   USE dtatsd          ! data temperature and salinity   (dta_tsd routine) 
    34    USE in_out_manager  ! I/O manager 
    35    USE iom             ! I/O library 
    3634   USE zpshde          ! partial step: hor. derivative (zps_hde routine) 
    3735   USE eosbn2          ! equation of state            (eos bn2 routine) 
     
    4240   USE dynspg_ts       ! pressure gradient schemes 
    4341   USE sol_oce         ! ocean solver variables 
     42   ! 
     43   USE in_out_manager  ! I/O manager 
     44   USE iom             ! I/O library 
    4445   USE lib_mpp         ! MPP library 
    4546   USE restart         ! restart 
     
    5657#  include "vectopt_loop_substitute.h90" 
    5758   !!---------------------------------------------------------------------- 
    58    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     59   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
    5960   !! $Id$ 
    6061   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    8182      CALL dta_tsd_init                       ! Initialisation of T & S input data 
    8283 
    83       rhd  (:,:,:  ) = 0.e0 
    84       rhop (:,:,:  ) = 0.e0 
    85       rn2  (:,:,:  ) = 0.e0  
    86       tsa  (:,:,:,:) = 0.e0     
     84      rhd  (:,:,:  ) = 0._wp 
     85      rhop (:,:,:  ) = 0._wp 
     86      rn2  (:,:,:  ) = 0._wp 
     87      tsa  (:,:,:,:) = 0._wp    
     88      rab_b(:,:,:,:) = 0._wp 
     89      rab_n(:,:,:,:) = 0._wp 
    8790 
    8891      IF( ln_rstart ) THEN                    ! Restart from a file 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r3876 r3893  
    404404      IF( lrst_oce )   CALL flt_rst( kt, 'WRITE' )      ! write filtered free surface arrays in restart file 
    405405      ! 
    406       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt') 
     406      IF( nn_timing == 1 )   CALL timing_stop('dyn_spg_flt') 
    407407      ! 
    408408   END SUBROUTINE dyn_spg_flt 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r3848 r3893  
    2828   USE zdfmxl         ! mixed layer depth 
    2929   USE eosbn2         ! equation of states 
     30   ! 
     31   USE in_out_manager ! I/O manager 
    3032   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    31    USE in_out_manager ! I/O manager 
    3233   USE prtctl         ! Print control 
    3334   USE wrk_nemo       ! work arrays 
     
    393394      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 
    394395      REAL(wp) ::   zdzrho_raw 
    395       REAL(wp) ::   zbeta0 
    396396      REAL(wp), POINTER, DIMENSION(:,:)     ::   z1_mlbw 
    397       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalbet 
    398397      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients 
    399398      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only 
     
    403402      ! 
    404403      CALL wrk_alloc( jpi,jpj, z1_mlbw ) 
    405       CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 
    406404      CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    407405      CALL wrk_alloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
     
    410408      !  Some preliminary calculation  ! 
    411409      !--------------------------------! 
    412       ! 
    413       CALL eos_alpbet( tsb, zalbet, zbeta0 )  !==  before local thermal/haline expension ratio at T-points  ==! 
    414410      ! 
    415411      DO jl = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==! 
     
    419415            DO jj = 1, jpjm1                  ! NB: not masked ==>  a minimum value is set 
    420416               DO ji = 1, fs_jpim1            ! vector opt. 
    421                   zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point 
    422                   zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    423                   zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point 
    424                   zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    425                   zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
    426                   zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
    427                   zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     417                  zdit = ( tsb(ji+1,jj  ,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point 
     418                  zdis = ( tsb(ji+1,jj  ,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     419                  zdjt = ( tsb(ji  ,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point 
     420                  zdjs = ( tsb(ji  ,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     421                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) / e1u(ji,jj) 
     422                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj) 
     423                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    428424                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
    429425               END DO 
     
    431427         END DO 
    432428         ! 
    433          IF( ln_zps.and.l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom 
     429         IF( ln_zps .AND. l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom 
    434430# if defined key_vectopt_loop 
    435431            DO jj = 1, 1 
     
    442438                  zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature 
    443439                  zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity 
    444                   zdxrho_raw = ( - zalbet(ji+ip,jj   ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
    445                   zdyrho_raw = ( - zalbet(ji   ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
     440                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) / e1u(ji,jj) 
     441                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj) 
    446442                  zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    447443                  zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     
    463459                     zdks = 0._wp 
    464460                  ENDIF 
    465                   zdzrho_raw = ( - zalbet(ji   ,jj   ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 
    466                   zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )    ! force zdzrho >= repsln 
     461                  zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp) 
     462                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw )    ! force zdzrho >= repsln 
    467463                 END DO 
    468464            END DO 
     
    608604      ! 
    609605      CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 
    610       CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 
    611606      CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    612607      CALL wrk_dealloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r3876 r3893  
    1717   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function 
    1818   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    19    !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp 
     19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_ab used in ldfslp 
    2020   !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add pen_ddt_dds used in trdpen 
    2121   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state 
    22    !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta 
     22   !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta to be used in dynhpg 
    2323   !!---------------------------------------------------------------------- 
    2424 
    2525   !!---------------------------------------------------------------------- 
    2626   !!   eos            : generic interface of the equation of state 
    27    !!   eos_insitu     : Compute the in situ density 
    28    !!   eos_insitu_pot : Compute the insitu and surface referenced potential 
    29    !!                    volumic mass 
    30    !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    31    !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
    32    !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio 
    33    !!   eos_expansion_ratio : calculates the partial derivative of in situ density with respect to T and S 
    34    !!   tfreez         : Compute the surface freezing temperature 
     27   !!   eos_insitu     : density anomaly (rho/rho0 - 1) 
     28   !!   eos_insitu_pot : density and surface referenced potential density anomalies 
     29   !!   eos_insitu_2d  : density anomaly for 2d fields 
     30   !!   eos_bn2        : Brunt-Vaisala frequency 
     31   !!   eos_rab        : partial derivative of density anomaly with respect to T and S 
     32   !!   tfreez         : surface freezing temperature 
    3533   !!   eos_init       : set eos parameters (namelist) 
    3634   !!---------------------------------------------------------------------- 
     
    4947   !                   !! * Interface 
    5048   INTERFACE eos 
    51       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
     49      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, & 
     50         &                eos_ab, eos_insitu_ab 
    5251   END INTERFACE 
    5352   INTERFACE bn2 
    54       MODULE PROCEDURE eos_bn2 
     53      MODULE PROCEDURE eos_bn2, eos_bn2_ab 
    5554   END INTERFACE 
    5655 
    5756   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules 
    5857   PUBLIC   eos_init       ! called by istate module 
    59    PUBLIC   bn2            ! called by step module 
    60    PUBLIC   eos_alpbet     ! called by ldfslp module 
    61    PUBLIC   eos_alpha_beta ! used for density diagnostics in dyntra 
     58   PUBLIC   bn2            ! called by step module, update alpha/beta and compute bn2 
    6259   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen 
    6360   PUBLIC   tfreez         ! called by sbcice_... modules 
    6461 
    6562   !                                          !!* Namelist (nameos) * 
    66    INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1/2/3 type of eq. of state and associated Brunt-Vaisala frequ. 
    67    REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 
    68    REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state) 
    69  
    70    REAL(wp) ::   vallis_ratio = 0   ! 1027/rau0 
    71    REAL(wp) ::   vallis_diff  = 0   ! (1027-rau0)/rau0 
     63   INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1 type of eq. of state and Brunt-Vaisala frequ. 
     64 
     65   REAL(wp) ::   rn_alph0  = 1.67e-4 !: thermal expansion coeff. (simplified equation of state) 
     66   REAL(wp) ::   rn_beta0  = 7.8e-4 !: saline  expansion coeff. (simplified equation of state) 
     67   REAL(wp) ::   rn_gamm0  = 4.4e-6 !: depth expansion coeff. (simplified equation of state) 
     68   REAL(wp) ::   rn_cabbel = 3.0e-2 !: cabbeling coeff. (simplified equation of state) 
     69   REAL(wp) ::   rn_thermo = 1.1e-4 !: thermobaric coeff. (simplified equation of state) 
     70   REAL(wp) ::   offset    = 0      !: offset applied on Vallis density anomaly, so that rho(10,35,0)=1027 
    7271 
    7372   !! * Substitutions 
     
    8988      !!       defined through the namelist parameter nn_eos. 
    9089      !! 
    91       !! ** Method  :   5 cases: 
     90      !! ** Method  :   3 cases: 
    9291      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
    9392      !!         Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 
     
    10099      !!         along geopotential surfaces, i.e. the pressure p in decibars 
    101100      !!         is approximated by the depth in meters. 
    102       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    103       !!         with pressure                      p        decibars 
     101      !!              prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
     102      !!         with depth                         z        meters 
    104103      !!              potential temperature         t        deg celsius 
    105104      !!              salinity                      s        psu 
     
    107106      !!              in situ volumic mass          rho      kg/m**3 
    108107      !!              in situ density anomalie      prd      no units 
    109       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
     108      !!         Check value: rho = 1060.93298 kg/m**3 for z=10000 dbar, 
    110109      !!          t = 40 deg celcius, s=40 psu 
    111       !!      nn_eos = 1 : linear equation of state function of temperature only 
    112       !!              prd(t) = 0.0285 - rn_alpha * t 
    113       !!      nn_eos = 2 : linear equation of state function of temperature and 
    114       !!               salinity 
    115       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
    116       !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
    117       !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0 
    118       !!              where the pressure p in decibars is approximated by the depth in meters. 
     110      !!      nn_eos = 1 : Vallis simplified equation of state (Vallis book, 2006) 
     111      !!              prd(t,s,z) =  ( rho(t,s,p) - rau0 ) / rau0 
     112      !!                         = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z 
     113      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
     114      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     115      !!              Vallis (2006) equation: use default values of coefficients  
    119116      !!      Note that no boundary condition problem occurs in this routine 
    120117      !!      as pts are defined over the whole domain. 
     
    122119      !! ** Action  :   compute prd , the in situ density (no units) 
    123120      !! 
    124       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    125       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    126       !!---------------------------------------------------------------------- 
     121      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     122      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     123      !!---------------------------------------------------------------------- 
     124      !! 
    127125      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    128126      !                                                      ! 2 : salinity               [psu] 
    129127      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-] 
    130       ! 
     128      !! 
    131129      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    132130      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    133       REAL(wp) ::   zr1, zr2, zr3, zr4   !   -      - 
    134       REAL(wp) ::   zrhop, ze, zbw, zb   !   -      - 
    135       REAL(wp) ::   zd , zc , zaw, za    !   -      - 
    136       REAL(wp) ::   zb1, za1, zkw, zk0   !   -      - 
     131      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     132      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     133      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     134      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    137135      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    138136      !!---------------------------------------------------------------------- 
     137 
    139138      ! 
    140139      IF( nn_timing == 1 ) CALL timing_start('eos') 
    141140      ! 
    142       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    143       ! 
    144141      SELECT CASE( nn_eos ) 
    145142      ! 
    146143      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
     144         ! 
     145         CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    147146         ! 
    148147!CDIR NOVERRCHK 
     
    158157                  ! 
    159158                  ! compute volumic mass pure water at atm pressure 
    160                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     159                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    161160                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    162161                  ! seawater volumic mass atm pressure 
    163                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     162                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    164163                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    165                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    166                   zr4= 4.8314e-4_wp 
    167                   ! 
    168                   ! potential volumic mass (reference to the surface) 
    169                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     164                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     165                  zn4= 4.8314e-4_wp 
     166                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    170167                  ! 
    171168                  ! add the compression terms 
    172                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    173                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    174                   zb = zbw + ze * zs 
    175                   ! 
    176                   zd = -1.480266e-4_wp 
    177                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    178                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    179                   za = ( zd*zsr + zc ) *zs + zaw 
    180                   ! 
    181                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    182                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    183                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    184                   zk0= ( zb1*zsr + za1 )*zs + zkw 
     169                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     170                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     171                  za = za1 + za2 * zs 
     172                  ! 
     173                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     174                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     175                  zb3 = 1.480266e-4_wp 
     176                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
     177                  ! 
     178                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     179                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     180                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     181                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
    185182                  ! 
    186183                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
    187                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    188                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     184                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     185                     &             - 1._wp  ) * tmask(ji,jj,jk) 
    189186               END DO 
    190187            END DO 
    191188         END DO 
    192189         ! 
     190         CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     191         ! 
    193192      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
     193         ! 
     194         CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    194195         ! 
    195196!CDIR NOVERRCHK 
     
    205206                  ! 
    206207                  ! compute volumic mass pure water at atm pressure 
    207                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     208                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    208209                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    209210                  ! seawater volumic mass atm pressure 
    210                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     211                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    211212                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    212                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    213                   zr4= 4.8314e-4_wp 
    214                   ! 
    215                   ! potential volumic mass (reference to the surface) 
    216                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     213                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     214                  zn4= 4.8314e-4_wp 
     215                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    217216                  ! 
    218217                  ! add the compression terms 
    219                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    220                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    221                   zb = zbw + ze * zs 
    222                   ! 
    223                   zd = -2.042967e-2_wp 
    224                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    225                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    226                   za = ( zd*zsr + zc ) *zs + zaw 
    227                   ! 
    228                   zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    229                   za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    230                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    231                   zk0= ( zb1*zsr + za1 )*zs + zkw 
     218                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     219                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     220                  za = za1 + za2 * zs 
     221                  ! 
     222                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     223                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     224                  zb3= 2.042967e-2_wp 
     225                  zb = ( zb3*zsr + zb2 ) *zs + zb1 
     226                  ! 
     227                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     228                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     229                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     230                  zc = ( zc3*zsr + zc2 )*zs + zc1 
    232231                  ! 
    233232                  ! masked in situ density anomaly 
    234                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    235                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     233                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     234                     &             - 1._wp  ) * tmask(ji,jj,jk) 
    236235               END DO 
    237236            END DO 
    238237         END DO 
    239238         ! 
    240       CASE( 1 )                !==  Linear formulation function of temperature only  ==! 
    241          DO jk = 1, jpkm1 
    242             prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    243          END DO 
    244          ! 
    245       CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==! 
    246          DO jk = 1, jpkm1 
    247             prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    248          END DO 
    249          ! 
    250       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     239         CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     240         ! 
     241      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    251242         DO jk = 1, jpkm1 
    252243            DO jj = 1, jpj 
     
    256247                  zh = fsdept(ji,jj,jk)        ! depth 
    257248                  ! masked in situ density anomaly 
    258                   prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             & 
    259                      &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
    260                      &   ) * tmask(ji,jj,jk) 
     249                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
     250                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
    261251               END DO 
    262252            END DO 
     
    267257      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk ) 
    268258      ! 
    269       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    270       ! 
    271259      IF( nn_timing == 1 ) CALL timing_stop('eos') 
    272260      ! 
     
    278266      !!                  ***  ROUTINE eos_insitu_pot  *** 
    279267      !! 
    280       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
     268      !! ** Purpose :   Compute the in situ density (ratio (rho-rau0)/rau0) and the 
    281269      !!      potential volumic mass (Kg/m3) from potential temperature and 
    282270      !!      salinity fields using an equation of state defined through the 
    283271      !!     namelist parameter nn_eos. 
    284272      !! 
    285       !! ** Method  :   5 cases: 
    286       !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
    287       !!      nn_eos = 0 : standard NEMO equation of state, based on  
    288       !!         Jackett and McDougall (1995) equation of state. 
    289       !!         the in situ density is computed directly as a function of 
    290       !!         potential temperature relative to the surface (the opa t 
    291       !!         variable), salt and pressure (assuming no pressure variation 
    292       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    293       !!         is approximated by the depth in meters. 
    294       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    295       !!              rhop(t,s)  = rho(t,s,0) 
    296       !!         with pressure                      p        decibars 
    297       !!              potential temperature         t        deg celsius 
    298       !!              salinity                      s        psu 
    299       !!              reference volumic mass        rau0     kg/m**3 
    300       !!              in situ volumic mass          rho      kg/m**3 
    301       !!              in situ density anomalie      prd      no units 
    302       !! 
    303       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
    304       !!          t = 40 deg celcius, s=40 psu 
    305       !! 
    306       !!      nn_eos = 1 : linear equation of state function of temperature only 
    307       !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t 
    308       !!              rhop(t,s)  = rho(t,s) 
    309       !! 
    310       !!      nn_eos = 2 : linear equation of state function of temperature and 
    311       !!               salinity 
    312       !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
    313       !!                       = rn_beta * s - rn_alpha * tn - 1. 
    314       !!              rhop(t,s)  = rho(t,s) 
    315       !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
    316       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    317       !!              rhop(t,s)  = rho(t,s,0) 
    318       !!      Note that no boundary condition problem occurs in this routine 
    319       !!      as (tn,sn) or (ta,sa) are defined over the whole domain. 
     273      !! ** Method  :   Same as for eos_insitu 
    320274      !! 
    321275      !! ** Action  : - prd  , the in situ density (no units) 
     
    332286      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
    333287      ! 
    334       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    335       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars 
    336       REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0               !   -      - 
     288      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     289      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
     290      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     291      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     292      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     293      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    337294      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    338295      !!---------------------------------------------------------------------- 
     
    357314                  ! 
    358315                  ! compute volumic mass pure water at atm pressure 
    359                   zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    360                      &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     316                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     317                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    361318                  ! seawater volumic mass atm pressure 
    362                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
    363                      &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
    364                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    365                   zr4= 4.8314e-4_wp 
    366                   ! 
    367                   ! potential volumic mass (reference to the surface) 
    368                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    369                   ! 
    370                   ! save potential volumic mass 
    371                   prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 
     319                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     320                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     321                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     322                  zn4= 4.8314e-4_wp 
     323                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    372324                  ! 
    373325                  ! add the compression terms 
    374                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    375                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    376                   zb = zbw + ze * zs 
    377                   ! 
    378                   zd = -1.480266e-4_wp 
    379                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    380                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    381                   za = ( zd*zsr + zc ) *zs + zaw 
    382                   ! 
    383                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    384                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    385                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    386                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    387                   ! 
    388                   ! masked in situ density anomaly 
    389                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    390                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     326                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     327                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     328                  za = za1 + za2 * zs 
     329                  ! 
     330                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     331                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     332                  zb3 = 1.480266e-4_wp 
     333                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
     334                  ! 
     335                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     336                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     337                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     338                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
     339                  ! 
     340                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
     341                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     342                     &             - 1._wp  ) * tmask(ji,jj,jk) 
     343              prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
    391344               END DO 
    392345            END DO 
     
    406359                  ! 
    407360                  ! compute volumic mass pure water at atm pressure 
    408                   zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    409                      &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     361                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     362                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    410363                  ! seawater volumic mass atm pressure 
    411                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
    412                      &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
    413                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    414                   zr4= 4.8314e-4_wp 
    415                   ! 
    416                   ! potential volumic mass (reference to the surface) 
    417                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    418                   ! 
    419                   ! save potential volumic mass 
    420                   prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 
     364                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     365                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     366                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     367                  zn4= 4.8314e-4_wp 
     368                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    421369                  ! 
    422370                  ! add the compression terms 
    423                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    424                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    425                   zb = zbw + ze * zs 
    426                   ! 
    427                   zd = -2.042967e-2_wp 
    428                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    429                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    430                   za = ( zd*zsr + zc ) *zs + zaw 
    431                   ! 
    432                   zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp 
    433                   za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp 
    434                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    435                   zk0= ( zb1*zsr + za1 )*zs + zkw 
     371                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     372                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     373                  za = za1 + za2 * zs 
     374                  ! 
     375                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     376                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     377                  zb3= 2.042967e-2_wp 
     378                  zb = ( zb3*zsr + zb2 ) *zs + zb1 
     379                  ! 
     380                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     381                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     382                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     383                  zc = ( zc3*zsr + zc2 )*zs + zc1 
    436384                  ! 
    437385                  ! masked in situ density anomaly 
    438                   prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    439                      &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     386                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     387                     &             - 1._wp  ) * tmask(ji,jj,jk) 
     388              prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
    440389               END DO 
    441390            END DO 
    442391         END DO 
    443392         ! 
    444       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    445          DO jk = 1, jpkm1 
    446             prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    447             prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
    448          END DO 
    449          ! 
    450       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    451          DO jk = 1, jpkm1 
    452             prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    453             prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
    454          END DO 
    455          ! 
    456       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     393      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    457394         DO jk = 1, jpkm1 
    458395            DO jj = 1, jpj 
     
    462399                  zh = fsdept(ji,jj,jk)        ! depth 
    463400                  ! masked in situ density anomaly 
    464                   prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             & 
    465                      &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
    466                      &   ) * tmask(ji,jj,jk) 
     401                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
     402                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
    467403                  ! masked in situ density anomaly 
    468                  prhop(ji,jj,jk) = ( 1.0_wp - ( 1.67e-4_wp + 0.5e-5_wp*zt ) *zt + 0.78e-3_wp *zs )   & 
    469                      &              * 1027._wp * tmask(ji,jj,jk) 
     404                 prhop(ji,jj,jk) = ( 1.0_wp + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt ) *zt + rn_beta0 *zs )  & 
     405                     &              * rau0 * tmask(ji,jj,jk) 
    470406                 ! 
    471407               END DO 
     
    492428      !!      defined through the namelist parameter nn_eos. * 2D field case 
    493429      !! 
    494       !! ** Method  :   5 cases: 
    495       !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
    496       !!      nn_eos = 0 : standard NEMO equation of state, based on  
    497       !!         Jackett and McDougall (1995) equation of state. 
    498       !!         the in situ density is computed directly as a function of 
    499       !!         potential temperature relative to the surface (the opa t 
    500       !!         variable), salt and pressure (assuming no pressure variation 
    501       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    502       !!         is approximated by the depth in meters. 
    503       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
    504       !!         with pressure                      p        decibars 
    505       !!              potential temperature         t        deg celsius 
    506       !!              salinity                      s        psu 
    507       !!              reference volumic mass        rau0     kg/m**3 
    508       !!              in situ volumic mass          rho      kg/m**3 
    509       !!              in situ density anomalie      prd      no units 
    510       !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 
    511       !!          t = 40 deg celcius, s=40 psu 
    512       !!      nn_eos = 1 : linear equation of state function of temperature only 
    513       !!              prd(t) = 0.0285 - rn_alpha * t 
    514       !!      nn_eos = 2 : linear equation of state function of temperature and 
    515       !!               salinity 
    516       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
    517       !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
    518       !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0 
    519       !!              where the pressure p in decibars is approximated by the depth in meters. 
    520       !!      Note that no boundary condition problem occurs in this routine 
    521       !!      as pts are defined over the whole domain. 
     430      !! ** Method  :   Same as for eos_insitu 
    522431      !! 
    523432      !! ** Action  : - prd , the in situ density (no units) 
     
    532441      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    533442      !! 
    534       INTEGER  ::   ji, jj                    ! dummy loop indices 
    535       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary scalars 
    536       REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask        !    -         - 
     443      INTEGER  ::   ji, jj              ! dummy loop indices 
     444      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
     445      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     446      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     447      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     448      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    537449      REAL(wp), POINTER, DIMENSION(:,:) :: zws 
    538450      !!---------------------------------------------------------------------- 
     
    558470         DO jj = 1, jpjm1 
    559471            DO ji = 1, fs_jpim1   ! vector opt. 
    560                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
    561472               zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
    562473               zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
     
    565476               ! 
    566477               ! compute volumic mass pure water at atm pressure 
    567                zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    568                   &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     478               zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     479                  &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    569480               ! seawater volumic mass atm pressure 
    570                zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
    571                   &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
    572                zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
    573                zr4 = 4.8314e-4_wp 
    574                ! 
    575                ! potential volumic mass (reference to the surface) 
    576                zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     481               zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     482                  &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     483               zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     484               zn4= 4.8314e-4_wp 
     485               zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    577486               ! 
    578487               ! add the compression terms 
    579                ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    580                zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    581                zb = zbw + ze * zs 
     488               za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     489               za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     490               za = za1 + za2 * zs 
    582491               ! 
    583                zd = -1.480266e-4_wp 
    584                zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    585                zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    586                za = ( zd*zsr + zc ) *zs + zaw 
     492               zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     493               zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     494               zb3 = 1.480266e-4_wp 
     495               zb  = ( zb3*zsr + zb2 ) *zs + zb1 
    587496               ! 
    588                zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    589                za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    590                zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    591                zk0= ( zb1*zsr + za1 )*zs + zkw 
     497               zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     498               zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     499               zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     500               zc  = ( zc3*zsr + zc2 )*zs + zc1 
    592501               ! 
    593502               ! masked in situ density anomaly 
    594                prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     503               prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     504                  &             - 1._wp  ) * tmask(ji,jj,1) 
    595505            END DO 
    596506         END DO 
     
    607517         DO jj = 1, jpjm1 
    608518            DO ji = 1, fs_jpim1   ! vector opt. 
    609                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
    610519               zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
    611520               zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
     
    614523               ! 
    615524               ! compute volumic mass pure water at atm pressure 
    616                zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
    617                   &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     525               zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     526                  &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    618527               ! seawater volumic mass atm pressure 
    619                zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
    620                   &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
    621                zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
    622                zr4 = 4.8314e-4_wp 
    623                ! 
    624                ! potential volumic mass (reference to the surface) 
    625                zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     528               zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     529                  &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     530               zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     531               zn4= 4.8314e-4_wp 
     532               zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    626533               ! 
    627534               ! add the compression terms 
    628                ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    629                zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    630                zb = zbw + ze * zs 
     535               za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     536               za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     537               za = za1 + za2 * zs 
    631538               ! 
    632                zd =    -2.042967e-2_wp 
    633                zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    634                zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 
    635                za = ( zd*zsr + zc ) *zs + zaw 
     539               zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     540               zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     541               zb3= 2.042967e-2_wp 
     542               zb = ( zb3*zsr + zb2 ) *zs + zb1 
    636543               ! 
    637                zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp 
    638                za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp 
    639                zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   & 
    640                   &                             +2098.925_wp ) *zt+190925.6_wp 
    641                zk0= ( zb1*zsr + za1 )*zs + zkw 
     544               zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     545               zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     546               zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     547               zc = ( zc3*zsr + zc2 )*zs + zc1 
    642548               ! 
    643549               ! masked in situ density anomaly 
    644                prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
    645             END DO 
    646          END DO 
    647          ! 
    648       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
     550               prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
     551                  &             - 1._wp  ) * tmask(ji,jj,1) 
     552            END DO 
     553         END DO 
     554         ! 
     555      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    649556         DO jj = 1, jpjm1 
    650557            DO ji = 1, fs_jpim1   ! vector opt. 
    651                prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    652             END DO 
    653          END DO 
    654          ! 
    655       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    656          DO jj = 1, jpjm1 
    657             DO ji = 1, fs_jpim1   ! vector opt. 
    658                prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    659             END DO 
    660          END DO 
    661          ! 
    662       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
    663          DO jj = 1, jpjm1 
    664             DO ji = 1, fs_jpim1   ! vector opt. 
    665                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
    666558               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! interpolated T 
    667559               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! interpolated S 
    668560               zh    = pdep (ji,jj)            ! depth at the partial step level 
    669561               ! masked in situ density anomaly 
    670                prd(ji,jj) = vallis_diff + vallis_ratio * (                                                & 
    671                   &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
    672                   &   ) * zmask 
     562               prd(ji,jj) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt  &                                           & 
     563                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,1) 
    673564               ! 
    674565            END DO 
     
    686577 
    687578 
    688    SUBROUTINE eos_bn2( pts, pn2 ) 
    689       !!---------------------------------------------------------------------- 
    690       !!                  ***  ROUTINE eos_bn2  *** 
    691       !! 
    692       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time- 
    693       !!      step of the input arguments 
    694       !! 
    695       !! ** Method  :   5 cases: 
    696       !!       * nn_eos = -1  : The brunt-vaisala frequency is computed for 
    697       !!       the Jackett and McDougall (1995) equation of state: 
    698       !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
    699       !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    700       !!      computed and used in zdfddm module : 
    701       !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
    702       !!       * nn_eos = 0  : The brunt-vaisala frequency is based on the 
    703       !!            formulation of McDougall (1987). 
    704       !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    705       !!      computed and used in zdfddm module : 
    706       !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
    707       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    708       !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
    709       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    710       !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    711       !!      The use of potential density to compute N^2 introduces e r r o r 
    712       !!      in the sign of N^2 at great depths. We recommand the use of 
    713       !!      nn_eos = 0, except for academical studies. 
    714       !!        Macro-tasked on horizontal slab (jk-loop) 
    715       !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
    716       !!      and is never used at this level. 
    717       !! 
    718       !! ** Action  : - pn2 : the brunt-vaisala frequency 
    719       !! 
    720       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    721       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    722       !!                McDougall, J. Phys. Oceano., 1987 
    723       !!---------------------------------------------------------------------- 
    724       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    725       !                                                               ! 2 : salinity               [psu] 
    726       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1] 
    727       !! 
    728       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    729       REAL(wp) ::   zgde3w, zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop      ! local scalars 
    730       REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM    !   -      - 
    731       REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom   !   -      - 
    732       REAL(wp) ::   zalpbet, zalpha, zbeta                                  !   -      - 
    733 #if defined key_zdfddm 
    734       REAL(wp) ::   zds   ! local scalars 
    735 #endif 
     579 
     580   SUBROUTINE eos_ab( pts, pab ) 
     581      !!---------------------------------------------------------------------- 
     582      !!                 ***  ROUTINE eos_ab  *** 
     583      !! 
     584      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points 
     585      !! 
     586      !! ** Method  :   calculates alpha and beta at T-points 
     587      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
     588      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
     589      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006) 
     590      !!                       alpha and beta ratios are returned as 3-D arrays. 
     591      !! 
     592      !! ** Action  : - pab : partial derivative of in situ density with respect to T & S at T-points 
     593      !!---------------------------------------------------------------------- 
     594      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
     595      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! alpha and beta 
     596      ! 
     597      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     598      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
     599      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     600      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     601      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     602      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
     603      REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
     604      REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
    736605      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    737606      !!---------------------------------------------------------------------- 
    738  
    739       ! 
    740       IF( nn_timing == 1 ) CALL timing_start('bn2') 
     607      ! 
     608      IF ( nn_timing == 1 )   CALL timing_start('eos_ab') 
    741609      ! 
    742610      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    743611      ! 
    744       ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    745       ! -------------------------- 
    746       ! 
    747       SELECT CASE( nn_eos ) 
    748       ! 
    749       CASE( -1 )                !==  Jackett and McDougall (1995)  ==! 
    750          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    751          DO jk = 2, jpkm1 
    752             DO jj = 1, jpj 
    753                DO ji = 1, jpi 
    754                   zgde3w = grav / fse3w(ji,jj,jk) 
    755                   zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt 
    756                   zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) )         ! salinity at w-pt 
    757                   zsr = 0.5 * ( zws(ji,jj,jk) + zws(ji,jj,jk-1) )        ! square root of S at w-pt 
    758                   zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    759                   ! 
    760                   ! compute volumic mass pure water at atm pressure 
    761                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    762                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    763                   ! seawater volumic mass atm pressure 
    764                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    765                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    766                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    767                   zr4= 4.8314e-4_wp 
    768                   ! 
    769                   ! potential volumic mass (reference to the surface) 
    770                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
    771                   ! 
    772                   ! add the compression terms 
    773                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    774                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    775                   zb = zbw + ze * zs 
    776                   ! 
    777                   zd = -1.480266e-4_wp 
    778                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    779                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    780                   za = ( zd*zsr + zc ) *zs + zaw 
    781                   ! 
    782                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    783                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    784                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    785                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    786                   ! 
    787                   zM= ( zb*zh - za )*zh + zk0 
    788                   zDenom= 1._wp - zh / zM 
    789                   ! compute in-situ density 
    790                   zrho = zrhop/zDenom 
    791               ! alpha 
    792                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    793                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    794                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    795                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    796                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    797                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    798                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    799                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    800                   zdDenom = zh * zdM / (zM*zM) 
    801                   zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
    802                   ! beta 
    803                   zdzb  = ze 
    804                   zdza  = -3.0644505e-2_wp*zsr+zc 
    805                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    806                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    807                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    808                   zdDenom = zh * zdM / (zM*zM) 
    809                   zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
    810                   ! alpha/beta 
    811                   zalpbet = zalpha / zbeta 
    812                   ! N2 
    813                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    814                      &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    815                      &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    816 #if defined key_zdfddm 
    817                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    818                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    819                   IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    820                   rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    821 #endif 
    822                END DO 
    823             END DO 
    824          END DO 
    825          ! 
    826       CASE( 0 )                !==  McDougall (1987) formulation  ==! 
    827          DO jk = 2, jpkm1 
    828             DO jj = 1, jpj 
    829                DO ji = 1, jpi 
    830                   zgde3w = grav / fse3w(ji,jj,jk) 
    831                   zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt 
    832                   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 
    833                   zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    834                   ! 
    835                   zalpbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt  &   ! ratio alpha/beta 
    836                      &                                  - 0.203814e-03_wp ) * zt   & 
    837                      &                                  + 0.170907e-01_wp ) * zt   & 
    838                      &   +         0.665157e-01_wp                                 & 
    839                      &   +     ( - 0.678662e-05_wp * zs                            & 
    840                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    841                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
    842                      &           - 0.251520e-11_wp * zs                            & 
    843                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    844                      &           - 0.164759e-06_wp * zs                            & 
    845                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    846                      &                                  + 0.380374e-04_wp ) * zh 
    847                      ! 
    848                   zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta 
    849                      &                               - 0.301985e-05_wp ) * zt      & 
    850                      &   +       0.785567e-03_wp                                   & 
    851                      &   + (     0.515032e-08_wp * zs                              & 
    852                      &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     & 
    853                      &   + ( (   0.121551e-17_wp * zh                              & 
    854                      &         - 0.602281e-15_wp * zs                              & 
    855                      &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     & 
    856                      &                                + 0.408195e-10_wp   * zs     & 
    857                      &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     & 
    858                      &                                - 0.121555e-07_wp ) * zh 
    859                      ! 
    860                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    861                      &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    862                      &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    863 #if defined key_zdfddm 
    864                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    865                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    866                   IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    867                   rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    868 #endif 
    869                END DO 
    870             END DO 
    871          END DO 
    872          ! 
    873       CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    874          DO jk = 2, jpkm1 
    875             pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )        & 
    876               &                / fse3w(:,:,jk) * tmask(:,:,jk) 
    877          END DO 
    878          ! 
    879       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    880          DO jk = 2, jpkm1 
    881             pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      & 
    882                &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   & 
    883                &               / fse3w(:,:,jk) * tmask(:,:,jk) 
    884          END DO 
    885 #if defined key_zdfddm 
    886          DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    887             DO jj = 1, jpj 
    888                DO ji = 1, jpi 
    889                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
    890                   IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    891                   rrau(ji,jj,jk) = rn_alpha / rn_beta * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    892                END DO 
    893             END DO 
    894          END DO 
    895 #endif 
    896          ! 
    897       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
    898          DO jk = 2, jpkm1 
    899             DO jj = 1, jpj 
    900                DO ji = 1, jpi 
    901                   zgde3w = grav / fse3w(ji,jj,jk) 
    902                   zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) - 10._wp ! potential temperature at w-pt 
    903                   zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35._wp ! salinity anomaly (s-35) at w-pt 
    904                   zh = fsdepw(ji,jj,jk)                                                 ! depth in meters  at w-point 
    905                   ! 
    906                   zalpha = 1027._wp * ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) / rau0 ! alpha 
    907                   zbeta  = 1027._wp * 0.78e-3_wp / rau0  ! beta 
    908                   zalpbet = zalpha / zbeta  ! ratio alpha/beta 
    909                   ! 
    910                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    911                      &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    912                      &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    913 #if defined key_zdfddm 
    914                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    915                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    916                   IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    917                   rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    918 #endif 
    919                END DO 
    920             END DO 
    921          END DO 
    922          ! 
    923       END SELECT 
    924  
    925       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
    926 #if defined key_zdfddm 
    927       IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
    928 #endif 
    929       ! 
    930       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    931       ! 
    932       IF( nn_timing == 1 ) CALL timing_stop('bn2') 
    933       ! 
    934    END SUBROUTINE eos_bn2 
    935  
    936  
    937    SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 
    938       !!---------------------------------------------------------------------- 
    939       !!                 ***  ROUTINE eos_alpbet  *** 
    940       !! 
    941       !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points 
    942       !! 
    943       !! ** Method  :   calculates alpha / beta ratio at T-points 
    944       !!       * nn_eos = -1  : alpha / beta ratio is computed for 
    945       !!       the Jackett and McDougall (1995) equation of state: 
    946       !!                       Scalar beta0 is returned = 1. 
    947       !!       * nn_eos = 0  : alpha / beta ratio using the formulation of McDougall (1987). 
    948       !!                       Scalar beta0 is returned = 1. 
    949       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    950       !!                       The ratio is undefined, so we return alpha as palpbet 
    951       !!                       Scalar beta0 is returned = 0. 
    952       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    953       !!                       The alpha/beta ratio is returned as palpbet 
    954       !!                       Scalar beta0 is returned = 1. 
    955       !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
    956       !!                       Scalar beta0 is returned = 1. 
    957       !! 
    958       !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points 
    959       !!            :   beta0   : 1. or 0. 
    960       !!---------------------------------------------------------------------- 
    961       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
    962       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio 
    963       REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T) 
    964       !! 
    965       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    966       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
    967       REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
    968       REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
    969       REAL(wp) ::   zalpha, zbeta                                ! local scalars 
    970       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    971       !!---------------------------------------------------------------------- 
    972       ! 
    973       IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 
    974       ! 
    975       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    976       ! 
    977612      SELECT CASE ( nn_eos ) 
    978613      ! 
    979614      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     615         ! 
    980616         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    981617         DO jk = 1, jpkm1 
     
    988624                  ! 
    989625                  ! compute volumic mass pure water at atm pressure 
    990                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     626                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    991627                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    992628                  ! seawater volumic mass atm pressure 
    993                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     629                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    994630                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    995                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    996                   zr4= 4.8314e-4_wp 
    997                   ! 
    998                   ! potential volumic mass (reference to the surface) 
    999                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     631                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     632                  zn4= 4.8314e-4_wp 
     633                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    1000634                  ! 
    1001635                  ! add the compression terms 
    1002                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    1003                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    1004                   zb = zbw + ze * zs 
    1005                   ! 
    1006                   zd = -1.480266e-4_wp 
    1007                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    1008                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    1009                   za = ( zd*zsr + zc ) *zs + zaw 
    1010                   ! 
    1011                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    1012                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    1013                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    1014                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    1015                   ! 
    1016                   zM= ( zb*zh - za )*zh + zk0 
    1017                   zDenom= 1._wp - zh / zM 
    1018                   ! compute in-situ density 
    1019                   zrho = zrhop/zDenom 
     636                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     637                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     638                  za = za1 + za2 * zs 
     639                  ! 
     640                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     641                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     642                  zb3 = 1.480266e-4_wp 
     643                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
     644                  ! 
     645                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     646                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
     647                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
     648                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
     649                  ! 
     650                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    1020651              ! alpha 
    1021                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1022                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    1023                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    1024                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    1025                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    1026                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    1027                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1028                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1029                   zdDenom = zh * zdM / (zM*zM) 
    1030                   zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
     652                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
     653                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
     654                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
     655                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
     656                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 
     657                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
     658                      
     659                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     660                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     661                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
     662                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
     663                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     664                  ! 
     665                  pab(ji,jj,jk,jp_tem) =  - ( & 
     666                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     667                     &   / rau0 * tmask(ji,jj,jk) 
     668                  ! 
    1031669                  ! beta 
    1032                   zdzb  = ze 
    1033                   zdza  = -3.0644505e-2_wp*zsr+zc 
    1034                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    1035                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    1036                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1037                   zdDenom = zh * zdM / (zM*zM) 
    1038                   zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
    1039                   ! alpha/beta 
    1040                   palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 
     670                  zdza  = za2 
     671                  zdzb  = 2.220399e-4_wp*zsr+zb2 
     672                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     673                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     674                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     675                  pab(ji,jj,jk,jp_sal) =    ( & 
     676                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     677                     &   / rau0 * tmask(ji,jj,jk) 
     678                  ! 
    1041679               END DO 
    1042680            END DO 
    1043681         END DO 
    1044          beta0 = 1._wp 
    1045          ! 
    1046       CASE ( 0 )               ! McDougall (1987) formulation 
    1047          DO jk = 1, jpk 
    1048             DO jj = 1, jpj 
    1049                DO ji = 1, jpi 
    1050                   zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
    1051                   zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    1052                   zh = fsdept(ji,jj,jk)               ! depth in meters 
    1053                   ! 
    1054                   palpbet(ji,jj,jk) =                                              & 
    1055                      &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
    1056                      &                                  - 0.203814e-03_wp ) * zt   & 
    1057                      &                                  + 0.170907e-01_wp ) * zt   & 
    1058                      &   + 0.665157e-01_wp                                         & 
    1059                      &   +     ( - 0.678662e-05_wp * zs                            & 
    1060                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    1061                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
    1062                      &           - 0.251520e-11_wp * zs                            & 
    1063                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    1064                      &           - 0.164759e-06_wp * zs                            & 
    1065                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    1066                      &                                  + 0.380374e-04_wp ) * zh 
    1067                END DO 
    1068             END DO 
    1069          END DO 
    1070          beta0 = 1._wp 
    1071          ! 
    1072       CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
    1073          palpbet(:,:,:) = rn_alpha 
    1074          beta0 = 0._wp 
    1075          ! 
    1076       CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
    1077          palpbet(:,:,:) = rn_alpha / rn_beta 
    1078          beta0 = 1._wp 
    1079          ! 
    1080       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
    1081          DO jk = 1, jpkm1 
    1082             DO jj = 1, jpj 
    1083                DO ji = 1, jpi 
    1084                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature 
    1085                   zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    1086                   zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    1087                   ! 
    1088                   zalpha = ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) ! alpha/vallis_ratio 
    1089                   zbeta  = 0.78e-3_wp  ! beta/vallis_ratio 
    1090                   ! 
    1091                   palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 
    1092                END DO 
    1093             END DO 
    1094          END DO 
    1095          beta0 = 1._wp 
    1096          ! 
    1097       CASE DEFAULT 
    1098          IF(lwp) WRITE(numout,cform_err) 
    1099          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
    1100          nstop = nstop + 1 
    1101          ! 
    1102       END SELECT 
    1103       ! 
    1104       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    1105       ! 
    1106       IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 
    1107       ! 
    1108    END SUBROUTINE eos_alpbet 
    1109  
    1110  
    1111    SUBROUTINE eos_alpha_beta( pts, palpha, pbeta ) 
    1112       !!---------------------------------------------------------------------- 
    1113       !!                 ***  ROUTINE eos_alpha_beta  *** 
    1114       !! 
    1115       !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points 
    1116       !! 
    1117       !! ** Method  :   calculates alpha and beta at T-points 
    1118       !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
    1119       !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
    1120       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    1121       !!                       palpha = rn_alpha 
    1122       !!                       pbeta = 0 
    1123       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    1124       !!                       palpha = rn_alpha 
    1125       !!                       pbeta = rn_beta 
    1126       !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
    1127       !!                       alpha and beta ratios are returned as 3-D arrays. 
    1128       !! 
    1129       !! ** Action  : - palpha : thermal expansion ratio at T-points 
    1130       !!            : - pbeta  : haline expansion ratio at T-points 
    1131       !!---------------------------------------------------------------------- 
    1132       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts      ! pot. temperature & salinity 
    1133       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpha   ! alpha 
    1134       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbeta    ! beta 
    1135       ! 
    1136       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    1137       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
    1138       REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
    1139       REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
    1140       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    1141       !!---------------------------------------------------------------------- 
    1142       ! 
    1143       IF ( nn_timing == 1 )   CALL timing_start('eos_alpha_beta') 
    1144       ! 
    1145       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    1146       ! 
    1147       SELECT CASE ( nn_eos ) 
    1148       ! 
    1149       CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     682         ! 
     683      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
    1150684         ! 
    1151685         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     
    1159693                  ! 
    1160694                  ! compute volumic mass pure water at atm pressure 
    1161                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     695                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    1162696                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    1163697                  ! seawater volumic mass atm pressure 
    1164                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     698                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    1165699                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    1166                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    1167                   zr4= 4.8314e-4_wp 
    1168                   ! 
    1169                   ! potential volumic mass (reference to the surface) 
    1170                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     700                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     701                  zn4= 4.8314e-4_wp 
     702                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    1171703                  ! 
    1172704                  ! add the compression terms 
    1173                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    1174                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    1175                   zb = zbw + ze * zs 
    1176                   ! 
    1177                   zd = -1.480266e-4_wp 
    1178                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    1179                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    1180                   za = ( zd*zsr + zc ) *zs + zaw 
    1181                   ! 
    1182                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    1183                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    1184                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    1185                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    1186                   ! 
    1187                   zM= ( zb*zh - za )*zh + zk0 
    1188                   zDenom= 1._wp - zh / zM 
    1189                   ! compute in-situ density 
    1190                   zrho = zrhop/zDenom 
     705                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     706                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     707                  za = za1 + za2 * zs 
     708                  ! 
     709                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     710                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     711                  zb3= 2.042967e-2_wp 
     712                  zb = ( zb3*zsr + zb2 ) *zs + zb1 
     713                  ! 
     714                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     715                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     716                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     717                  zc = ( zc3*zsr + zc2 )*zs + zc1 
     718                  ! 
     719                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    1191720              ! alpha 
    1192                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1193                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    1194                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    1195                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    1196                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    1197                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    1198                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1199                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1200                   zdDenom = zh * zdM / (zM*zM) 
    1201                   palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     721                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     722                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
     723                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     724                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     725                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     726                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     727                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     728                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     729                  ! 
     730                  pab(ji,jj,jk,jp_tem) =  - ( & 
     731                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     732                     &   / rau0 * tmask(ji,jj,jk) 
    1202733                  ! beta 
    1203                   zdzb  = ze 
    1204                   zdza  = -3.0644505e-2_wp*zsr+zc 
    1205                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    1206                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    1207                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1208                   zdDenom = zh * zdM / (zM*zM) 
    1209                   pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     734                  zdza  = za2 
     735                  zdzb  = 3.0644505e-2_wp*zsr+zb2 
     736                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     737                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     738                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     739                  pab(ji,jj,jk,jp_sal) =    ( & 
     740                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     741                     &   / rau0 * tmask(ji,jj,jk) 
     742                  ! 
    1210743               END DO 
    1211744            END DO 
    1212745         END DO 
    1213746         ! 
    1214       CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
     747      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     748         DO jk = 1, jpkm1 
     749            DO jj = 1, jpj 
     750               DO ji = 1, jpi 
     751                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature anomaly (t-T0) 
     752                  zh = fsdept(ji,jj,jk)               ! depth in meters at t-point 
     753                  ! 
     754                  pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh ) * tmask(ji,jj,jk)   ! alpha 
     755               END DO 
     756            END DO 
     757         END DO 
     758         pab(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta 
     759         ! 
     760      CASE DEFAULT 
     761         IF(lwp) WRITE(numout,cform_err) 
     762         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     763         nstop = nstop + 1 
     764         ! 
     765      END SELECT 
     766      ! 
     767      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     768      ! 
     769      IF( nn_timing == 1 )   CALL timing_stop('eos_ab') 
     770      ! 
     771   END SUBROUTINE eos_ab 
     772 
     773 
     774 
     775   SUBROUTINE eos_insitu_ab( pts, prd, pab ) 
     776      !!---------------------------------------------------------------------- 
     777      !!                 ***  ROUTINE eos_insitu_ab  *** 
     778      !! 
     779      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms 
     780      !!                  and the insitu density anomaly at T-points 
     781      !! 
     782      !! ** Method  :   calculates rhd, alpha, beta at T-points 
     783      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
     784      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
     785      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006) 
     786      !!                       alpha and beta ratios are returned as 3-D arrays. 
     787      !! 
     788      !! ** Action  : - pab : thermal and haline expansion ratios at T-points 
     789      !!              - prd , the density anomaly (no units) 
     790      !! 
     791      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     792      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     793      !!---------------------------------------------------------------------- 
     794      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity    [Celcius,psu] 
     795      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   prd   ! density anomaly                [-] 
     796      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! partial derivative of density anomaly 
     797      !                                                               ! with respect to T and S        [1/Celcius,1/psu] 
     798      ! 
     799      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     800      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
     801      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     802      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     803      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     804      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
     805      REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
     806      REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
     807      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
     808      !!---------------------------------------------------------------------- 
     809      ! 
     810      IF ( nn_timing == 1 )   CALL timing_start('eos_insitu_ab') 
     811      ! 
     812      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     813      ! 
     814      SELECT CASE ( nn_eos ) 
     815      ! 
     816      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
    1215817         ! 
    1216818         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     
    1224826                  ! 
    1225827                  ! compute volumic mass pure water at atm pressure 
    1226                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     828                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    1227829                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    1228830                  ! seawater volumic mass atm pressure 
    1229                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     831                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    1230832                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    1231                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    1232                   zr4= 4.8314e-4_wp 
    1233                   ! 
    1234                   ! potential volumic mass (reference to the surface) 
    1235                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     833                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     834                  zn4= 4.8314e-4_wp 
     835                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    1236836                  ! 
    1237837                  ! add the compression terms 
    1238                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    1239                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    1240                   zb = zbw + ze * zs 
    1241                   ! 
    1242                   zd = -2.042967e-2_wp 
    1243                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    1244                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    1245                   za = ( zd*zsr + zc ) *zs + zaw 
    1246                   ! 
    1247                   zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    1248                   za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    1249                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    1250                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    1251                   ! 
    1252                   zM= ( zb*zh - za )*zh + zk0 
    1253                   zDenom= 1._wp - zh / zM 
    1254                   ! compute in-situ density 
    1255                   zrho = zrhop/zDenom 
     838                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     839                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     840                  za = za1 + za2 * zs 
     841                  ! 
     842                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     843                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     844                  zb3 = 1.480266e-4_wp 
     845                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
     846                  ! 
     847                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     848                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
     849                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
     850                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
     851                  ! 
     852                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
     853                  ! 
     854                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
     855                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) / rau0  - 1._wp  ) * tmask(ji,jj,jk) 
     856                  ! 
    1256857              ! alpha 
    1257                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1258                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    1259                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    1260                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    1261                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    1262                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    1263                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1264                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1265                   zdDenom = zh * zdM / (zM*zM) 
    1266                   palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     858                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
     859                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
     860                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
     861                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
     862                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 
     863                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
     864                      
     865                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     866                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     867                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
     868                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
     869                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     870                  ! 
     871                  pab(ji,jj,jk,jp_tem) =  - ( & 
     872                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     873                     &   / rau0 * tmask(ji,jj,jk) 
     874                  ! 
    1267875                  ! beta 
    1268                   zdzb  = ze 
    1269                   zdza  = -3.0644505e-2_wp*zsr+zc 
    1270                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    1271                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    1272                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1273                   zdDenom = zh * zdM / (zM*zM) 
    1274                   pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     876                  zdza  = za2 
     877                  zdzb  = 2.220399e-4_wp*zsr+zb2 
     878                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     879                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     880                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     881                  pab(ji,jj,jk,jp_sal) =    ( & 
     882                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
     883                     &   / rau0 * tmask(ji,jj,jk) 
     884                  ! 
    1275885               END DO 
    1276886            END DO 
    1277887         END DO 
    1278888         ! 
    1279       CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
    1280          palpha(:,:,:) = rn_alpha 
    1281          pbeta (:,:,:) = 0._wp 
    1282          ! 
    1283       CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
    1284          palpha(:,:,:) = rn_alpha 
    1285          pbeta (:,:,:) = rn_beta 
    1286          ! 
    1287       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
    1288          DO jk = 1, jpkm1 
    1289             DO jj = 1, jpj 
    1290                DO ji = 1, jpi 
    1291                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10) 
    1292                   zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point 
    1293                   palpha(ji,jj,jk) = vallis_ratio * & 
    1294                       &  ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alpha 
    1295                   ! 
    1296                END DO 
    1297             END DO 
    1298          END DO 
    1299          pbeta (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! beta 
    1300          ! 
    1301       CASE DEFAULT 
    1302          IF(lwp) WRITE(numout,cform_err) 
    1303          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
    1304          nstop = nstop + 1 
    1305          ! 
    1306       END SELECT 
    1307       ! 
    1308       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    1309       ! 
    1310       IF( nn_timing == 1 ) CALL timing_stop('eos_alpha_beta') 
    1311       ! 
    1312    END SUBROUTINE eos_alpha_beta 
    1313  
    1314  
    1315  
    1316    SUBROUTINE eos_pen( pts, palphaPE, pbetaPE, pPEanom ) 
    1317       !!---------------------------------------------------------------------- 
    1318       !!                 ***  ROUTINE eos_pen  *** 
    1319       !! 
    1320       !! ** Purpose :   Calculates alpha_PE, beta_PE and PE at T-points 
    1321       !! 
    1322       !! ** Method  :   PE is defined analytically as the vertical  
    1323       !!                   primitive of EOS times -g integrated between 0 and z>0. 
    1324       !!                PEanom is the PE anomaly: PEanom = ( PE - rau0 gz ) / rau0 gz 
    1325       !!                                           = -1/z * /int_z^0 rd , where rd density anomaly 
    1326       !!                dalphaPE and dbetaPE are partial derivatives of PE anomaly with respect to T and S: 
    1327       !!                    alphaPE = - 1/(rau0 gz) * dPE/dT = - dPEanom/dT 
    1328       !!                    betaPE  =   1/(rau0 gz) * dPE/dS = - dPEanom/dS 
    1329       !! 
    1330       !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
    1331       !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state. 
    1332       !!       * nn_eos = 1  : linear equation of state (temperature only) 
    1333       !!                       palpha = rau0*g*rn_alpha*z 
    1334       !!                       pbeta = 0 
    1335       !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    1336       !!                       palpha = rau0*g*rn_alpha*z 
    1337       !!                       pbeta = -rau0*g*rn_beta*z 
    1338       !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
    1339       !! 
    1340       !! ** Action  : - pPEanom     : PE anomaly given at T-points 
    1341       !!            : - palphaPE    : given at T-points 
    1342       !!            : - pbetaPE     : given at T-points 
    1343       !!---------------------------------------------------------------------- 
    1344       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
    1345       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palphaPE  ! partial derivative of PE anom 
    1346       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbetaPE   ! with respect to T and S, resp. 
    1347       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pPEanom   ! potential energy anomaly 
    1348       ! 
    1349       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    1350       REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
    1351       REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
    1352       REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
    1353       REAL(wp) ::   zap1, zdelta, zsdelta, zAp, zBp, zlogBp, zCp, zatanCp ! 
    1354       REAL(wp) ::   zddelta, zdAp, zdBp, zdlogBp, zdCp, zP, zdP, zEp      ! 
    1355       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    1356       !!---------------------------------------------------------------------- 
    1357       ! 
    1358       IF ( nn_timing == 1 )   CALL timing_start('eos_pen') 
    1359       ! 
    1360       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    1361       ! 
    1362       SELECT CASE ( nn_eos ) 
    1363       ! 
    1364       CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     889      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
    1365890         ! 
    1366891         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     
    1374899                  ! 
    1375900                  ! compute volumic mass pure water at atm pressure 
    1376                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     901                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    1377902                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    1378903                  ! seawater volumic mass atm pressure 
    1379                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     904                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    1380905                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    1381                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    1382                   zr4= 4.8314e-4_wp 
    1383                   ! 
    1384                   ! potential volumic mass (reference to the surface) 
    1385                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     906                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     907                  zn4= 4.8314e-4_wp 
     908                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    1386909                  ! 
    1387910                  ! add the compression terms 
    1388                   ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    1389                   zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    1390                   zb = zbw + ze * zs 
    1391                   ! 
    1392                   zd = -1.480266e-4_wp 
    1393                   zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
    1394                   zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
    1395                   za = ( zd*zsr + zc ) *zs + zaw 
    1396                   ! 
    1397                   zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    1398                   za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    1399                   zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
    1400                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    1401                   ! 
    1402                   zM= ( zb*zh - za )*zh + zk0 
    1403                   zDenom= 1._wp - zh / zM 
    1404                   ! compute in-situ density 
    1405                   zrho = zrhop/zDenom 
    1406                   ! 
    1407                   zap1    = 1._wp + za 
    1408                   zdelta  = 4._wp*zb*zk0 - zap1*zap1 
    1409                   zsdelta = sqrt( abs( zdelta ) ) 
    1410                   zAp     = zap1 / zb / zsdelta 
    1411                   zBp     = ( zM - zh ) / zk0 
    1412                   zlogBp  = log(abs(zBp)) 
    1413                   zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1) 
    1414                   zatanCp = atan(zCp) 
    1415                   zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 
    1416                   ! compute potential energy 
    1417                   pPEanom(ji,jj,jk)     = ( ( zrhop * zP ) / rau0 / zh - 1._wp ) * tmask(ji,jj,jk) !!wrong 
    1418                   ! 
    1419               ! dPEdt 
    1420                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1421                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    1422                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    1423                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    1424                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    1425                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    1426                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1427                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1428                   zdDenom = zh * zdM / (zM*zM) 
    1429                   ! 
    1430                   zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
    1431                   zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
    1432                   zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
    1433                   zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
    1434                   zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
    1435                   ! 
    1436                   palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk)  
    1437                   ! 
    1438               ! dPEds 
    1439                   zdzb  = ze 
    1440                   zdza  = -3.0644505e-2_wp*zsr+zc 
    1441                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    1442                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    1443                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1444                   zdDenom = zh * zdM / (zM*zM) 
    1445                   ! 
    1446                   zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
    1447                   zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
    1448                   zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
    1449                   zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
    1450                   zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
    1451                   ! 
    1452                   pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
    1453                   ! 
     911                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     912                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     913                  za = za1 + za2 * zs 
     914                  ! 
     915                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     916                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     917                  zb3= 2.042967e-2_wp 
     918                  zb = ( zb3*zsr + zb2 ) *zs + zb1 
     919                  ! 
     920                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     921                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     922                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     923                  zc = ( zc3*zsr + zc2 )*zs + zc1 
     924                  ! 
     925                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
     926                  ! 
     927                  ! masked in situ density anomaly (rho/rho0 - 1). Important: rau0=1035, should be 1027! 
     928                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) * r1_rau0  - 1._wp  ) * tmask(ji,jj,jk) 
     929                  ! 
     930                  !                                 ! alpha 
     931                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     932                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
     933                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     934                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     935                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     936                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     937                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     938                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     939                  ! 
     940                  pab(ji,jj,jk,jp_tem) =  - (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     941                     &                    * r1_rau0 * tmask(ji,jj,jk) 
     942                  ! 
     943                  !                                 ! beta 
     944                  zdza  = za2 
     945                  zdzb  = 3.0644505e-2_wp*zsr+zb2 
     946                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     947                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     948                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     949                  pab(ji,jj,jk,jp_sal) = (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     950                     &                 * r1_rau0 * tmask(ji,jj,jk) 
    1454951               END DO 
    1455952            END DO 
    1456953         END DO 
    1457954         ! 
    1458       CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
     955      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     956         DO jk = 1, jpkm1 
     957            DO jj = 1, jpj 
     958               DO ji = 1, jpi 
     959                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0) 
     960                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point 
     961                  ! masked in situ density anomaly 
     962                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
     963                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
     964                  ! alpha 
     965                  pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh )  & 
     966                      &                       * tmask(ji,jj,jk)   ! alpha 
     967                  ! 
     968               END DO 
     969            END DO 
     970         END DO 
     971         ! beta 
     972         pab (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta 
     973         ! 
     974      CASE DEFAULT 
     975         IF(lwp) WRITE(numout,cform_err) 
     976         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     977         nstop = nstop + 1 
     978         ! 
     979      END SELECT 
     980      ! 
     981      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     982      ! 
     983      IF( nn_timing == 1 ) CALL timing_stop('eos_insitu_ab') 
     984      ! 
     985   END SUBROUTINE eos_insitu_ab 
     986 
     987 
     988 
     989   SUBROUTINE eos_pen( pts, pab_pe, pen ) 
     990      !!---------------------------------------------------------------------- 
     991      !!                 ***  ROUTINE eos_pen  *** 
     992      !! 
     993      !! ** Purpose :   Calculates alpha_PE, beta_PE and PE at T-points 
     994      !! 
     995      !! ** Method  :   PE is defined analytically as the vertical  
     996      !!                   primitive of EOS times -g integrated between 0 and z>0. 
     997      !!                pen is the PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz 
     998      !!                                           = -1/z * /int_z^0 rd , where rd density anomaly 
     999      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S: 
     1000      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT = - d(pen)/dT 
     1001      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS = - d(pen)/dS 
     1002      !! 
     1003      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
     1004      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state. 
     1005      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006) 
     1006      !! 
     1007      !! ** Action  : - pen         : PE anomaly given at T-points 
     1008      !!            : - pab_pe  : given at T-points 
     1009      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe 
     1010      !!                    pab_pe(:,:,:,jp_sal) is beta_pe 
     1011      !!---------------------------------------------------------------------- 
     1012      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity 
     1013      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe 
     1014      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pen     ! potential energy anomaly 
     1015      ! 
     1016      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     1017      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
     1018      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     1019      REAL(wp) ::   zn,  za1, za2, za    !   -      - 
     1020      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
     1021      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
     1022      REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
     1023      REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
     1024      REAL(wp) ::   zdelta, zsdelta      !   -      - 
     1025      REAL(wp) ::   zsgn, zAp, zlogBp    !   -      - 
     1026      REAL(wp) ::   zCp, zatCp, zP       !   -      - 
     1027      REAL(wp) ::   zddelta, zdAp        !   -      - 
     1028      REAL(wp) ::   zdlogBp, zdCp, zdP   !   -      - 
     1029      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
     1030      !!---------------------------------------------------------------------- 
     1031      ! 
     1032      IF ( nn_timing == 1 )   CALL timing_start('eos_pen') 
     1033      ! 
     1034      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     1035      ! 
     1036      SELECT CASE ( nn_eos ) 
     1037      ! 
     1038      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
    14591039         ! 
    14601040         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     
    14641044                  zt = pts   (ji,jj,jk,jp_tem) 
    14651045                  zs = pts   (ji,jj,jk,jp_sal) 
    1466                   zh = fsdept(ji,jj,jk)        ! depth 
     1046                  zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0) 
    14671047                  zsr= zws   (ji,jj,jk)        ! square root salinity 
    14681048                  ! 
    14691049                  ! compute volumic mass pure water at atm pressure 
    1470                   zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1050                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    14711051                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    14721052                  ! seawater volumic mass atm pressure 
    1473                   zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1053                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    14741054                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    1475                   zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    1476                   zr4= 4.8314e-4_wp 
    1477                   ! 
    1478                   ! potential volumic mass (reference to the surface) 
    1479                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1055                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1056                  zn4= 4.8314e-4_wp 
     1057                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    14801058                  ! 
    14811059                  ! add the compression terms 
    1482                   ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    1483                   zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    1484                   zb = zbw + ze * zs 
    1485                   ! 
    1486                   zd = -2.042967e-2_wp 
    1487                   zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
    1488                   zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    1489                   za = ( zd*zsr + zc ) *zs + zaw 
    1490                   ! 
    1491                   zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    1492                   za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    1493                   zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    1494                   zk0= ( zb1*zsr + za1 )*zs + zkw 
    1495                   ! 
    1496                   zM= ( zb*zh - za )*zh + zk0 
    1497                   zDenom= 1._wp - zh / zM 
    1498                   ! compute in-situ density 
    1499                   zrho = zrhop/zDenom 
    1500                   ! 
    1501                   zap1    = 1._wp + za 
    1502                   zdelta  = 4._wp*zb*zk0 - zap1*zap1 
    1503                   zsdelta = sqrt( abs( zdelta ) ) 
    1504                   zAp     = zap1 / zb / zsdelta 
    1505                   zBp     = ( zM - zh ) / zk0 
    1506                   zlogBp  = log(abs(zBp)) 
    1507                   zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1) 
    1508                   zatanCp = atan(zCp) 
    1509                   zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 
     1060                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     1061                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     1062                  za = za1 + za2 * zs 
     1063                  ! 
     1064                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
     1065                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
     1066                  zb3 = 1.480266e-4_wp 
     1067                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
     1068                  ! 
     1069                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     1070                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
     1071                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
     1072                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
     1073                  ! 
     1074                  zdelta = 4._wp * za * zc - zb * zb 
     1075                  ! ddt 
     1076                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
     1077                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
     1078                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
     1079                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
     1080                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 
     1081                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
     1082                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1083                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1084                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
     1085                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1086                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     1087                  ! 
     1088                  ! stability condition for the calculation (add a small perturbation on T if needed) 
     1089                  IF( ABS(zdelta) < 1.e-5 .OR. ABS(za)<1.e-9 ) THEN ; 
     1090                     zt = zt + 1.e-3_wp 
     1091                     ! 
     1092                     zn = zn + 1.e-3_wp * zdn 
     1093                     za = za + 1.e-3_wp * zdza 
     1094                     zb = zb + 1.e-3_wp * zdzb 
     1095                     zc = zc + 1.e-3_wp * zdzc 
     1096                     ! 
     1097                     zdelta = 4._wp * za * zc - zb * zb 
     1098                     ! 
     1099                  ENDIF 
     1100                  ! 
     1101                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
     1102                  zsgn   = SIGN( 1. , zdelta ) 
     1103                  zsdelta = SQRT( zsgn * zdelta ) 
     1104                  zAp     = - zb / za / zsdelta 
     1105                  zlogBp  = - LOG( zc * zm_inv ) 
     1106                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb) 
     1107                  IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp); 
     1108                  ELSE  ;    zatCp = ATANH(zCp) 
     1109                  ENDIF 
     1110                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 
     1111                  ! 
    15101112                  ! compute potential energy 
    1511                   pPEanom(ji,jj,jk)     = - grav * zrhop * zP * tmask(ji,jj,jk) 
    1512                   ! 
    1513               ! dPEdt 
    1514                   zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1515                   zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
    1516                   zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    1517                      &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    1518                   zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    1519                      &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    1520                      &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1521                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1522                   zdDenom = zh * zdM / (zM*zM) 
    1523                   ! 
    1524                   zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
    1525                   zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
    1526                   zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
    1527                   zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
    1528                   zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
    1529                   ! 
    1530                   palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
    1531                   ! 
    1532               ! dPEds 
    1533                   zdzb  = ze 
    1534                   zdza  = -3.0644505e-2_wp*zsr+zc 
    1535                   zdzk0 = 1.5_wp*zb1*zsr+za1 
    1536                   zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
    1537                   zdM = (zdzb*zh-zdza)*zh+zdzk0 
    1538                   zdDenom = zh * zdM / (zM*zM) 
    1539                   ! 
    1540                   zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
    1541                   zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
    1542                   zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
    1543                   zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
    1544                   zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
    1545                   ! 
    1546                   pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
     1113                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk) 
     1114                  ! 
     1115              ! alphaPE 
     1116                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
     1117                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     1118                  zdlogBp = zdm*zm_inv - zdzc/zc 
     1119                  zdCp    = zCp * (  zddelta - (2._wp*zdzc+zh*zdzb) / (2._wp*zc+zh*zb)  ) 
     1120                  zdP     = (  0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp   & 
     1121                     &       + zAp*zdCp / ( 1._wp + zsgn*zCp*zCp )  ) / zh 
     1122                  ! 
     1123                  pab_pe(ji,jj,jk,jp_tem)    = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk) 
     1124                  ! dds 
     1125                  zdza  = za2 
     1126                  zdzb  = 2.220399e-4_wp*zsr+zb2 
     1127                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     1128                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     1129                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     1130                  ! betaPE 
     1131                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
     1132                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     1133                  zdlogBp = zdm*zm_inv - zdzc/zc 
     1134                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb)) 
     1135                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + & 
     1136                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
     1137                  ! 
     1138                  pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    & 
     1139                                                    &           / rau0 * tmask(ji,jj,jk) 
    15471140                  ! 
    15481141               END DO 
     
    15501143         END DO 
    15511144         ! 
    1552       CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
    1553          palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 
    1554          pbetaPE (:,:,:) = 0._wp 
    1555          DO jk = 1, jpkm1 
    1556             pPEanom(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    1557          END DO 
    1558          ! 
    1559       CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
    1560          palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 
    1561          pbetaPE (:,:,:) = rn_beta  * tmask(:,:,:) 
    1562          DO jk = 1, jpkm1 
    1563             pPEanom(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    1564          END DO 
    1565          ! 
    1566       CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     1145      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
     1146         ! 
     1147         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    15671148         DO jk = 1, jpkm1 
    15681149            DO jj = 1, jpj 
    15691150               DO ji = 1, jpi 
    1570                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10) 
    1571                   zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    1572                   zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point 
    1573                   ! 
    1574                   palphaPE(ji,jj,jk) = vallis_ratio * & 
    1575                       &  ( 1.67e-4_wp * ( 1._wp + 0.55e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alphaPE 
    1576                   pPEanom(ji,jj,jk)     = vallis_diff + vallis_ratio *   & 
    1577                       &   ( - ( 1.67e-4_wp*(1._wp+0.55e-4_wp*zh) + 0.5e-5_wp*zt )*zt + 0.78e-3_wp*zs + 2.195e-6_wp*zh ) & 
    1578                       &   * tmask(ji,jj,jk) 
     1151                  ! 
     1152                  zt = pts   (ji,jj,jk,jp_tem) 
     1153                  zs = pts   (ji,jj,jk,jp_sal) 
     1154                  zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0) 
     1155                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     1156                  ! 
     1157                  ! compute volumic mass pure water at atm pressure 
     1158                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1159                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     1160                  ! seawater volumic mass atm pressure 
     1161                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1162                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     1163                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1164                  zn4= 4.8314e-4_wp 
     1165                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
     1166                  ! 
     1167                  ! add the compression terms 
     1168                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     1169                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     1170                  za = za1 + za2 * zs 
     1171                  ! 
     1172                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
     1173                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
     1174                  zb3= 2.042967e-2_wp 
     1175                  zb = ( zb3*zsr + zb2 ) *zs + zb1 
     1176                  ! 
     1177                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     1178                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     1179                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     1180                  zc = ( zc3*zsr + zc2 )*zs + zc1 
     1181                  ! 
     1182                  zdelta = 4._wp * za * zc - zb * zb 
     1183                  ! 
     1184                  ! dds 
     1185                  zdza  = za2 
     1186                  zdzb  = 3.0644505e-2_wp*zsr+zb2 
     1187                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     1188                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
     1189                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
     1190                  ! stability condition for the calculation (add a small perturbation on S if needed) 
     1191                  IF( ABS(zdelta) < 1.e-5 .OR. ABS(za)<1.e-9 ) THEN ; !stability condition for the calculation 
     1192                     zs = zs + 1.e-3_wp 
     1193                     zsr= zsr + 0.5e-3_wp/zsr 
     1194                     ! 
     1195                     zn = zn + 1.e-3_wp * zdn 
     1196                     za = za + 1.e-3_wp * zdza 
     1197                     zb = zb + 1.e-3_wp * zdzb 
     1198                     zc = zc + 1.e-3_wp * zdzc 
     1199                     ! 
     1200                     zdelta = 4._wp * za * zc - zb * zb 
     1201                     ! 
     1202                  ENDIF 
     1203                  ! 
     1204                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
     1205                  zsgn   = SIGN( 1. , zdelta ) 
     1206                  zsdelta = SQRT( zsgn * zdelta ) 
     1207                  zAp     = - zb / za / zsdelta 
     1208                  zlogBp  = - LOG( zc * zm_inv ) 
     1209                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb) 
     1210                  IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp); 
     1211                  ELSE  ;    zatCp = ATANH(zCp) 
     1212                  ENDIF 
     1213                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 
     1214                  ! 
     1215                  ! compute potential energy 
     1216                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk) 
     1217                  ! 
     1218                  ! betaPE 
     1219                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
     1220                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     1221                  zdlogBp = zdm*zm_inv - zdzc/zc 
     1222                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb)) 
     1223                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + & 
     1224                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
     1225                  ! 
     1226                  pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    & 
     1227                                                    &           / rau0 * tmask(ji,jj,jk) 
     1228                  ! 
     1229              ! ddt 
     1230                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1231                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
     1232                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1233                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1234                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1235                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1236                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1237                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     1238              ! alphaPE 
     1239                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
     1240                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     1241                  zdlogBp = zdm*zm_inv - zdzc/zc 
     1242                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb)) 
     1243                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + & 
     1244                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
     1245                  ! 
     1246                  pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk) 
    15791247                  ! 
    15801248               END DO 
    15811249            END DO 
    15821250         END DO 
    1583          pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! betaPE=beta 
     1251         ! 
     1252      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     1253         DO jk = 1, jpkm1 
     1254            DO jj = 1, jpj 
     1255               DO ji = 1, jpi 
     1256                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0) 
     1257                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point 
     1258                  pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) )  & 
     1259                      &          * tmask(ji,jj,jk)                                     ! alphaPE 
     1260                  pen(ji,jj,jk) = ( - rn_alph0 * ( 1._wp + rn_cabbel*zt + 0.5_wp*rn_thermo*zh ) *zt   &                                          & 
     1261                      &          + rn_beta0 * zs + 0.5_wp*rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
     1262                  ! 
     1263               END DO 
     1264            END DO 
     1265         END DO 
     1266         pab_pe (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! betaPE 
    15841267         ! 
    15851268      CASE DEFAULT 
     
    15981281 
    15991282 
     1283   SUBROUTINE eos_bn2_ab( pts, pn2 ) 
     1284      !!---------------------------------------------------------------------- 
     1285      !!                  ***  ROUTINE eos_bn2_ab  *** 
     1286      !! 
     1287      !! ** Purpose :   Update alpha/beta then compute the local Brunt-Vaisala  
     1288      !!      frequency at the time-step of the input arguments 
     1289      !! 
     1290      !! ** Method  :   call eos_ab to obtain alpha and beta coefficients 
     1291      !!                then interpolate them vertically, and compute pn2 
     1292      !!        Macro-tasked on horizontal slab (jk-loop) 
     1293      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
     1294      !!      and is never used at this level. 
     1295      !! 
     1296      !! ** Action  : - pn2 : the brunt-vaisala frequency 
     1297      !! 
     1298      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     1299      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     1300      !!                McDougall, J. Phys. Oceano., 1987 
     1301      !!---------------------------------------------------------------------- 
     1302      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
     1303      !                                                               ! 2 : salinity               [psu] 
     1304      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1] 
     1305      !! 
     1306      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     1307      REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars 
     1308#if defined key_zdfddm 
     1309      REAL(wp) ::   zds   ! local scalars 
     1310#endif 
     1311      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zab 
     1312      !!---------------------------------------------------------------------- 
     1313      ! 
     1314      ! 
     1315      IF( nn_timing == 1 ) CALL timing_start('eos_bn2_ab') 
     1316      ! 
     1317      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
     1318      ! -------------------------- 
     1319      ! 
     1320      CALL wrk_alloc( jpi, jpj, jpk, jpts, zab ) 
     1321      ! 
     1322      CALL eos_ab( pts, zab ) 
     1323      ! 
     1324      DO jk = 2, jpkm1 
     1325         DO jj = 1, jpj 
     1326            DO ji = 1, jpi 
     1327               ! 
     1328               zgde3w = 1._wp / fse3w(ji,jj,jk) 
     1329               zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w 
     1330               ! 
     1331               za = ( zab(ji,jj,jk,jp_tem) +   & 
     1332                  &   ( zab(ji,jj,jk-1,jp_tem) - zab(ji,jj,jk,jp_tem) ) * zr1 )   & 
     1333                  &   * tmask(ji,jj,jk) 
     1334               zb = ( zab(ji,jj,jk,jp_sal) +   & 
     1335                  &   ( zab(ji,jj,jk-1,jp_sal) - zab(ji,jj,jk,jp_sal) ) * zr1 )   & 
     1336                  &   * tmask(ji,jj,jk) 
     1337               ! N2 
     1338               pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2 
     1339                  &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     1340                  &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     1341#if defined key_zdfddm 
     1342               zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
     1343               IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
     1344               rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     1345#endif 
     1346            END DO 
     1347         END DO 
     1348      END DO 
     1349      ! 
     1350 
     1351      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
     1352#if defined key_zdfddm 
     1353      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
     1354#endif 
     1355      ! 
     1356      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zab ) 
     1357      ! 
     1358      IF( nn_timing == 1 ) CALL timing_stop('eos_bn2_ab') 
     1359      ! 
     1360   END SUBROUTINE eos_bn2_ab 
     1361 
     1362 
     1363   SUBROUTINE eos_bn2( pts, pab, pn2 ) 
     1364      !!---------------------------------------------------------------------- 
     1365      !!                  ***  ROUTINE eos_bn2  *** 
     1366      !! 
     1367      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
     1368      !!                time-step of the input arguments 
     1369      !! 
     1370      !! ** Method  :   interpolate alpha/beta vertically, and compute pn2 
     1371      !!        Macro-tasked on horizontal slab (jk-loop) 
     1372      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
     1373      !!      and is never used at this level. 
     1374      !! 
     1375      !! ** Action  : - pn2 : the brunt-vaisala frequency 
     1376      !! 
     1377      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     1378      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     1379      !!                McDougall, J. Phys. Oceano., 1987 
     1380      !!---------------------------------------------------------------------- 
     1381      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
     1382      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! partial derivative of in situ  
     1383      !                                                              ! density with respect to T and S [1/Celcius,1/psu] 
     1384      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency    [1/s^2] 
     1385      !! 
     1386      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     1387      REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars 
     1388#if defined key_zdfddm 
     1389      REAL(wp) ::   zds   ! local scalars 
     1390#endif 
     1391      !!---------------------------------------------------------------------- 
     1392      ! 
     1393      ! 
     1394      IF( nn_timing == 1 ) CALL timing_start('eos_bn2') 
     1395      ! 
     1396      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
     1397      ! -------------------------- 
     1398      ! 
     1399      DO jk = 2, jpkm1 
     1400         DO jj = 1, jpj 
     1401            DO ji = 1, jpi 
     1402               zgde3w = 1._wp / fse3w(ji,jj,jk) 
     1403               zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w 
     1404               ! 
     1405               za = ( pab(ji,jj,jk,jp_tem) + ( pab(ji,jj,jk-1,jp_tem) - pab(ji,jj,jk,jp_tem) ) * zr1 ) * tmask(ji,jj,jk) 
     1406               zb = ( pab(ji,jj,jk,jp_sal) + ( pab(ji,jj,jk-1,jp_sal) - pab(ji,jj,jk,jp_sal) ) * zr1 ) * tmask(ji,jj,jk) 
     1407               ! 
     1408               pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2 
     1409                  &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     1410                  &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     1411#if defined key_zdfddm 
     1412               zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
     1413               IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
     1414               rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     1415#endif 
     1416            END DO 
     1417         END DO 
     1418      END DO 
     1419      ! 
     1420 
     1421      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
     1422#if defined key_zdfddm 
     1423      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
     1424#endif 
     1425      ! 
     1426      IF( nn_timing == 1 )   CALL timing_stop('eos_bn2') 
     1427      ! 
     1428   END SUBROUTINE eos_bn2 
     1429 
     1430 
     1431 
    16001432   FUNCTION tfreez( psal ) RESULT( ptf ) 
    16011433      !!---------------------------------------------------------------------- 
     
    16291461      !! ** Method  :   Read the namelist nameos and control the parameters 
    16301462      !!---------------------------------------------------------------------- 
    1631       NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta 
     1463      NAMELIST/nameos/ nn_eos, rn_alph0, rn_beta0, rn_gamm0, rn_cabbel, rn_thermo 
    16321464      !!---------------------------------------------------------------------- 
    16331465      ! 
     
    16401472         WRITE(numout,*) '~~~~~~~~' 
    16411473         WRITE(numout,*) '          Namelist nameos : set eos parameters' 
    1642          WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos 
    1643          WRITE(numout,*) '             thermal exp. coef. (linear)    rn_alpha = ', rn_alpha 
    1644          WRITE(numout,*) '             saline  exp. coef. (linear)    rn_beta  = ', rn_beta 
     1474         WRITE(numout,*) '             flag for eq. of state and N^2    nn_eos    = ', nn_eos 
     1475         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_alph0  = ', rn_alph0 
     1476         WRITE(numout,*) '             saline  exp. coef. (nn_eos=1)    rn_beta0  = ', rn_beta0 
     1477         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_gamm0  = ', rn_gamm0 
     1478         WRITE(numout,*) '             saline  exp. coef. (nn_eos=1)    rn_cabbel = ', rn_cabbel 
     1479         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_thermo = ', rn_thermo 
    16451480      ENDIF 
    16461481      ! 
     
    16561491         IF(lwp) WRITE(numout,*) '                and McDougall (1987) Brunt-Vaisala frequency' 
    16571492         ! 
    1658       CASE( 1 )                        !==  Linear formulation = F( temperature )  ==! 
     1493      CASE( 1 )                        !==  Vallis (2006) formulation  ==! 
    16591494         IF(lwp) WRITE(numout,*) 
    1660          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )' 
    1661          IF( lk_zdfddm ) CALL ctl_stop( '          double diffusive mixing parameterization requires',   & 
    1662               &                         ' that T and S are used as state variables' ) 
    1663          ! 
    1664       CASE( 2 )                        !==  Linear formulation = F( temperature , salinity )  ==! 
    1665          IF(lwp) WRITE(numout,*) 
    1666          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 
    1667          ! 
    1668       CASE( 3 )                        !==  Vallis (2006) formulation  ==! 
    1669          IF(lwp) WRITE(numout,*) 
    1670          IF(lwp) WRITE(numout,*) '          use of Vallis (2006) equation of state' 
    1671          vallis_ratio = 1027._wp / rau0 
    1672          vallis_diff = (1027._wp-rau0) / rau0 
     1495         IF(lwp) WRITE(numout,*) '          use of simplified eos:    rhd(T,S,z) = ' 
     1496         IF(lwp) WRITE(numout,*) '             -alph0*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta0*(S-S0) + gamm0*z' 
     1497         IF( lk_zdfddm .AND. rn_beta0==0. )   & 
     1498              &   CALL ctl_stop( '          double diffusive mixing parameterization requires that rn_beta<>0') 
     1499         offset = 1027._wp / rau0 - 1._wp 
    16731500         ! 
    16741501      CASE DEFAULT                     !==  ERROR in nn_eos  ==! 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90

    r3876 r3893  
    7979      ! 
    8080 
     81!!gm end 
    8182      IF( ln_tra_mxl .OR. ln_vor_trd )   CALL ctl_stop( 'ML tracer and Barotropic vorticity diags are still using old IOIPSL' ) 
    8283!!gm end 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

    r3876 r3893  
    276276      ENDIF 
    277277      ! 
    278       nkstp     = nit000 - 1 
    279 !!gm check the value of nkstp  here nit000 should be OK 
    280       ! 
    281278   END SUBROUTINE trd_ken_init 
    282279 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdpen.F90

    r3876 r3893  
    3333   INTEGER ::   nkstp   ! current time step  
    3434 
    35    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alphaPE, betaPE   ! partial derivative of PEanom with respect to T and S 
     35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   ab_pe   ! partial derivative of PE anomaly with respect to T and S 
    3636 
    3737   !! * Substitutions 
     
    5050      !!                  ***  FUNCTION trd_tra_alloc  *** 
    5151      !!--------------------------------------------------------------------- 
    52       ALLOCATE( alphaPE(jpi,jpj,jpk) , betaPE(jpi,jpj,jpk) , STAT= trd_pen_alloc ) 
     52      ALLOCATE( ab_pe(jpi,jpj,jpk,jpts) , STAT= trd_pen_alloc ) 
    5353      ! 
    5454      IF( lk_mpp             )   CALL mpp_sum ( trd_pen_alloc ) 
    55       IF( trd_pen_alloc /= 0 )   CALL ctl_warn('trd_pen_alloc: failed to allocate arrays') 
     55      IF( trd_pen_alloc /= 0 )   CALL ctl_warn( 'trd_pen_alloc: failed to allocate arrays' ) 
    5656   END FUNCTION trd_pen_alloc 
    5757 
     
    7979      IF ( kt /= nkstp ) THEN   ! full eos: set partial derivatives at the 1st call of kt time step 
    8080         nkstp = kt 
    81          CALL eos_pen( tsn, alphaPE, betaPE, zpe ) 
    82          CALL iom_put( "alphaPE", alphaPE ) 
    83          CALL iom_put( "betaPE", betaPE ) 
    84          CALL iom_put( "PEanom", zpe ) 
     81         CALL eos_pen( tsn, ab_PE, zpe ) 
     82         CALL iom_put( "alphaPE", ab_pe(:,:,:,jp_tem) ) 
     83         CALL iom_put( "betaPE" , ab_pe(:,:,:,jp_sal) ) 
     84         CALL iom_put( "PEanom" , zpe ) 
    8585      ENDIF 
    8686      ! 
    8787      DO jk = 1, jpkm1 
    88          zpe(:,:,jk) = ( - alphaPE(:,:,jk) * ptrdx(:,:,jk)   & 
    89             &            + betaPE (:,:,jk) * ptrdy(:,:,jk)  ) 
     88         zpe(:,:,jk) = ( - ab_pe(:,:,jk,jp_tem) * ptrdx(:,:,jk)   & 
     89            &            + ab_pe(:,:,jk,jp_sal) * ptrdy(:,:,jk)  ) 
    9090      END DO 
    9191 
     
    9696                                IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface 
    9797                                   CALL wrk_alloc( jpi, jpj, z2d ) 
    98                                    z2d(:,:) = wn(:,:,1) * ( - alphaPE(:,:,1) * tsn(:,:,1,jp_tem)    & 
    99                                       &                     + betaPE (:,:,1) * tsn(:,:,1,jp_sal)  ) / fse3t(:,:,1) 
     98                                   z2d(:,:) = wn(:,:,1) * ( - ab_pe(:,:,1,jp_tem) * tsn(:,:,1,jp_tem)    & 
     99                                      &                     + ab_pe(:,:,1,jp_sal) * tsn(:,:,1,jp_sal)  ) / fse3t(:,:,1) 
    100100                                   CALL iom_put( "petrd_sad" , z2d ) 
    101101                                   CALL wrk_dealloc( jpi, jpj, z2d ) 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r3625 r3893  
    2020   !! dynamics and tracer fields                            ! before ! now    ! after  ! the after trends becomes the fields 
    2121   !! --------------------------                            ! fields ! fields ! trends ! only after tra_zdf and dyn_spg 
    22    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s] 
    23    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
    24    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s] 
    25    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity           [s-1] 
    26    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdivb,  hdivn          !: horizontal divergence        [s-1] 
    27    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn            !: 4D T-S fields        [Celcius,psu]  
    28    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2   [s-2] 
    29    ! 
    30    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsa             !: 4D T-S trends fields & work array  
     22   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity          [m/s] 
     23   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity          [m/s] 
     24   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity              [m/s] 
     25   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity             [s-1] 
     26   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdivb,  hdivn          !: horizontal divergence          [s-1] 
     27   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                  [Celcius,psu]  
     28   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celcius-1,psu-1] 
     29   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2] 
    3130   ! 
    3231   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units] 
    3332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3] 
    3433 
    35    !! free surface                                      !  before  ! now    ! after  ! 
    36    !! ------------                                      !  fields  ! fields ! trends ! 
     34   !! free surface                                              !  before  ! now    ! after  ! 
     35   !! ------------                                              !  fields  ! fields ! trends ! 
    3736   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m] 
    3837   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 
     
    7574         &      hdivb(jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             & 
    7675         &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     & 
     76         &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             & 
    7777         &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
    7878         ! 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3878 r3893  
    104104      ! Ocean physics update                (ua, va, tsa used as workspace) 
    105105      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    106                          CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency 
    107                          CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency 
     106                         CALL eos( tsb, rab_b )              ! before thermal/haline expansion coef.  
     107                         CALL bn2( tsb, rab_b, rn2b )        ! before Brunt-Vaisala frequency 
     108                         CALL eos( tsn, rab_n )              ! now thermal/haline expansion coef.  
     109                         CALL bn2( tsn, rab_n, rn2  )        ! now    Brunt-Vaisala frequency 
     110                          
    108111      ! 
    109112      !  VERTICAL PHYSICS 
Note: See TracChangeset for help on using the changeset viewer.