source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

  • Property svn:keywords set to Id
File size: 89.0 KB
Line 
1MODULE eosbn2
2   !!==============================================================================
3   !!                       ***  MODULE  eosbn2  ***
4   !! Ocean diagnostic variable : equation of state - in situ and potential density
5   !!                                               - Brunt-Vaisala frequency
6   !!==============================================================================
7   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
8   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
9   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
10   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
11   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
12   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
13   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
14   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
15   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
16   !!             -   ! 2003-08  (G. Madec)  F90, free form
17   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function
18   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
20   !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add pen_ddt_dds used in trdpen
21   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state
22   !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta
23   !!----------------------------------------------------------------------
24
25   !!----------------------------------------------------------------------
26   !!   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
35   !!   eos_init       : set eos parameters (namelist)
36   !!----------------------------------------------------------------------
37   USE dom_oce         ! ocean space and time domain
38   USE phycst          ! physical constants
39   USE zdfddm          ! vertical physics: double diffusion
40   USE in_out_manager  ! I/O manager
41   USE lib_mpp         ! MPP library
42   USE prtctl          ! Print control
43   USE wrk_nemo        ! Memory Allocation
44   USE timing          ! Timing
45
46   IMPLICIT NONE
47   PRIVATE
48
49   !                   !! * Interface
50   INTERFACE eos
51      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
52   END INTERFACE
53   INTERFACE bn2
54      MODULE PROCEDURE eos_bn2
55   END INTERFACE
56
57   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules
58   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
62   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen
63   PUBLIC   tfreez         ! called by sbcice_... modules
64
65   !                                          !!* 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
72
73   !! * Substitutions
74#  include "domzgr_substitute.h90"
75#  include "vectopt_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
78   !! $Id$
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE eos_insitu( pts, prd )
84      !!----------------------------------------------------------------------
85      !!                   ***  ROUTINE eos_insitu  ***
86      !!
87      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
88      !!       potential temperature and salinity using an equation of state
89      !!       defined through the namelist parameter nn_eos.
90      !!
91      !! ** Method  :   5 cases:
92      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state.
93      !!         Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar,
94      !!          t = 3 deg celcius, s=35.5 psu
95      !!      nn_eos = 0 : standard NEMO equation of state, based on
96      !!         Jackett and McDougall (1995) equation of state.
97      !!         the in situ density is computed directly as a function of
98      !!         potential temperature relative to the surface (the opa t
99      !!         variable), salt and pressure (assuming no pressure variation
100      !!         along geopotential surfaces, i.e. the pressure p in decibars
101      !!         is approximated by the depth in meters.
102      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
103      !!         with pressure                      p        decibars
104      !!              potential temperature         t        deg celsius
105      !!              salinity                      s        psu
106      !!              reference volumic mass        rau0     kg/m**3
107      !!              in situ volumic mass          rho      kg/m**3
108      !!              in situ density anomalie      prd      no units
109      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
110      !!          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.
119      !!      Note that no boundary condition problem occurs in this routine
120      !!      as pts are defined over the whole domain.
121      !!
122      !! ** Action  :   compute prd , the in situ density (no units)
123      !!
124      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
125      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
126      !!----------------------------------------------------------------------
127      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
128      !                                                      ! 2 : salinity               [psu]
129      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-]
130      !
131      INTEGER  ::   ji, jj, jk           ! dummy loop indices
132      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   !   -      -
137      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
138      !!----------------------------------------------------------------------
139      !
140      IF( nn_timing == 1 ) CALL timing_start('eos')
141      !
142      CALL wrk_alloc( jpi, jpj, jpk, zws )
143      !
144      SELECT CASE( nn_eos )
145      !
146      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==!
147         !
148!CDIR NOVERRCHK
149         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
150         !
151         DO jk = 1, jpkm1
152            DO jj = 1, jpj
153               DO ji = 1, jpi
154                  zt = pts   (ji,jj,jk,jp_tem)
155                  zs = pts   (ji,jj,jk,jp_sal)
156                  zh = fsdept(ji,jj,jk)        ! depth
157                  zsr= zws   (ji,jj,jk)        ! square root salinity
158                  !
159                  ! 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   &
161                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
162                  ! seawater volumic mass atm pressure
163                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
164                     &                      -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
170                  !
171                  ! 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
185                  !
186                  ! 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)
189               END DO
190            END DO
191         END DO
192         !
193      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==!
194         !
195!CDIR NOVERRCHK
196         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
197         !
198         DO jk = 1, jpkm1
199            DO jj = 1, jpj
200               DO ji = 1, jpi
201                  zt = pts   (ji,jj,jk,jp_tem)
202                  zs = pts   (ji,jj,jk,jp_sal)
203                  zh = fsdept(ji,jj,jk)        ! depth
204                  zsr= zws   (ji,jj,jk)        ! square root salinity
205                  !
206                  ! 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                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
209                  ! seawater volumic mass atm pressure
210                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
211                     &                      -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
217                  !
218                  ! 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
232                  !
233                  ! 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)
236               END DO
237            END DO
238         END DO
239         !
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  ==!
251         DO jk = 1, jpkm1
252            DO jj = 1, jpj
253               DO ji = 1, jpi
254                  zt = pts   (ji,jj,jk,jp_tem) - 10._wp
255                  zs = pts   (ji,jj,jk,jp_sal) - 35._wp
256                  zh = fsdept(ji,jj,jk)        ! depth
257                  ! 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)
261               END DO
262            END DO
263         END DO
264         !
265      END SELECT
266      !
267      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk )
268      !
269      CALL wrk_dealloc( jpi, jpj, jpk, zws )
270      !
271      IF( nn_timing == 1 ) CALL timing_stop('eos')
272      !
273   END SUBROUTINE eos_insitu
274
275
276   SUBROUTINE eos_insitu_pot( pts, prd, prhop )
277      !!----------------------------------------------------------------------
278      !!                  ***  ROUTINE eos_insitu_pot  ***
279      !!
280      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
281      !!      potential volumic mass (Kg/m3) from potential temperature and
282      !!      salinity fields using an equation of state defined through the
283      !!     namelist parameter nn_eos.
284      !!
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.
320      !!
321      !! ** Action  : - prd  , the in situ density (no units)
322      !!              - prhop, the potential volumic mass (Kg/m3)
323      !!
324      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
325      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
326      !!                Brown and Campana, Mon. Weather Rev., 1978
327      !!----------------------------------------------------------------------
328      !!
329      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
330      !                                                                ! 2 : salinity               [psu]
331      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
332      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
333      !
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               !   -      -
337      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
338      !!----------------------------------------------------------------------
339      !
340      IF( nn_timing == 1 ) CALL timing_start('eos-p')
341      !
342      CALL wrk_alloc( jpi, jpj, jpk, zws )
343      !
344      SELECT CASE ( nn_eos )
345      !
346      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==!
347!CDIR NOVERRCHK
348         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
349         !
350         DO jk = 1, jpkm1
351            DO jj = 1, jpj
352               DO ji = 1, jpi
353                  zt = pts   (ji,jj,jk,jp_tem)
354                  zs = pts   (ji,jj,jk,jp_sal)
355                  zh = fsdept(ji,jj,jk)        ! depth
356                  zsr= zws   (ji,jj,jk)        ! square root salinity
357                  !
358                  ! 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
361                  ! 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)
372                  !
373                  ! 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)
391               END DO
392            END DO
393         END DO
394         !
395      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==!
396!CDIR NOVERRCHK
397         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
398         !
399         DO jk = 1, jpkm1
400            DO jj = 1, jpj
401               DO ji = 1, jpi
402                  zt = pts   (ji,jj,jk,jp_tem)
403                  zs = pts   (ji,jj,jk,jp_sal)
404                  zh = fsdept(ji,jj,jk)        ! depth
405                  zsr= zws   (ji,jj,jk)        ! square root salinity
406                  !
407                  ! 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
410                  ! 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)
421                  !
422                  ! 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
436                  !
437                  ! 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)
440               END DO
441            END DO
442         END DO
443         !
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  ==!
457         DO jk = 1, jpkm1
458            DO jj = 1, jpj
459               DO ji = 1, jpi
460                  zt = pts   (ji,jj,jk,jp_tem) - 10._wp
461                  zs = pts   (ji,jj,jk,jp_sal) - 35._wp
462                  zh = fsdept(ji,jj,jk)        ! depth
463                  ! 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)
467                  ! 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)
470                 !
471               END DO
472            END DO
473         END DO
474         !
475      END SELECT
476      !
477      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
478      !
479      CALL wrk_dealloc( jpi, jpj, jpk, zws )
480      !
481      IF( nn_timing == 1 ) CALL timing_stop('eos-p')
482      !
483   END SUBROUTINE eos_insitu_pot
484
485
486   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
487      !!----------------------------------------------------------------------
488      !!                  ***  ROUTINE eos_insitu_2d  ***
489      !!
490      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
491      !!      potential temperature and salinity using an equation of state
492      !!      defined through the namelist parameter nn_eos. * 2D field case
493      !!
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.
522      !!
523      !! ** Action  : - prd , the in situ density (no units)
524      !!
525      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
526      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
527      !!----------------------------------------------------------------------
528      !!
529      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
530      !                                                           ! 2 : salinity               [psu]
531      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m]
532      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
533      !!
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        !    -         -
537      REAL(wp), POINTER, DIMENSION(:,:) :: zws
538      !!----------------------------------------------------------------------
539      !
540      IF( nn_timing == 1 ) CALL timing_start('eos2d')
541      !
542      CALL wrk_alloc( jpi, jpj, zws )
543      !
544
545      prd(:,:) = 0._wp
546
547      SELECT CASE( nn_eos )
548      !
549      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==!
550      !
551!CDIR NOVERRCHK
552         DO jj = 1, jpjm1
553!CDIR NOVERRCHK
554            DO ji = 1, fs_jpim1   ! vector opt.
555               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
556            END DO
557         END DO
558         DO jj = 1, jpjm1
559            DO ji = 1, fs_jpim1   ! vector opt.
560               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
561               zt    = pts  (ji,jj,jp_tem)     ! interpolated T
562               zs    = pts  (ji,jj,jp_sal)     ! interpolated S
563               zsr   = zws  (ji,jj)            ! square root of interpolated S
564               zh    = pdep (ji,jj)            ! depth at the partial step level
565               !
566               ! 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
569               ! 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
577               !
578               ! 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
582               !
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
587               !
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
592               !
593               ! masked in situ density anomaly
594               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
595            END DO
596         END DO
597         !
598      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==!
599      !
600!CDIR NOVERRCHK
601         DO jj = 1, jpjm1
602!CDIR NOVERRCHK
603            DO ji = 1, fs_jpim1   ! vector opt.
604               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
605            END DO
606         END DO
607         DO jj = 1, jpjm1
608            DO ji = 1, fs_jpim1   ! vector opt.
609               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
610               zt    = pts  (ji,jj,jp_tem)     ! interpolated T
611               zs    = pts  (ji,jj,jp_sal)     ! interpolated S
612               zsr   = zws  (ji,jj)            ! square root of interpolated S
613               zh    = pdep (ji,jj)            ! depth at the partial step level
614               !
615               ! 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
618               ! 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
626               !
627               ! 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
631               !
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
636               !
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
642               !
643               ! 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 )  ==!
649         DO jj = 1, jpjm1
650            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
666               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! interpolated T
667               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! interpolated S
668               zh    = pdep (ji,jj)            ! depth at the partial step level
669               ! 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
673               !
674            END DO
675         END DO
676         !
677      END SELECT
678
679      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
680      !
681      CALL wrk_dealloc( jpi, jpj, zws )
682      !
683      IF( nn_timing == 1 ) CALL timing_stop('eos2d')
684      !
685   END SUBROUTINE eos_insitu_2d
686
687
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
736      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
737      !!----------------------------------------------------------------------
738
739      !
740      IF( nn_timing == 1 ) CALL timing_start('bn2')
741      !
742      CALL wrk_alloc( jpi, jpj, jpk, zws )
743      !
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      !
977      SELECT CASE ( nn_eos )
978      !
979      CASE ( -1 )               ! Jackett and McDougall (1995) formulation
980         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
981         DO jk = 1, jpkm1
982            DO jj = 1, jpj
983               DO ji = 1, jpi
984                  zt = pts   (ji,jj,jk,jp_tem)
985                  zs = pts   (ji,jj,jk,jp_sal)
986                  zh = fsdept(ji,jj,jk)        ! depth
987                  zsr= zws   (ji,jj,jk)        ! square root salinity
988                  !
989                  ! 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   &
991                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
992                  ! seawater volumic mass atm pressure
993                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
994                     &                      -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
1000                  !
1001                  ! 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
1020              ! 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
1031                  ! 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)
1041               END DO
1042            END DO
1043         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
1150         !
1151         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1152         DO jk = 1, jpkm1
1153            DO jj = 1, jpj
1154               DO ji = 1, jpi
1155                  zt = pts   (ji,jj,jk,jp_tem)
1156                  zs = pts   (ji,jj,jk,jp_sal)
1157                  zh = fsdept(ji,jj,jk)        ! depth
1158                  zsr= zws   (ji,jj,jk)        ! square root salinity
1159                  !
1160                  ! 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   &
1162                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
1163                  ! seawater volumic mass atm pressure
1164                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
1165                     &                      -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
1171                  !
1172                  ! 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
1191              ! 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)
1202                  ! 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)
1210               END DO
1211            END DO
1212         END DO
1213         !
1214      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
1215         !
1216         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1217         DO jk = 1, jpkm1
1218            DO jj = 1, jpj
1219               DO ji = 1, jpi
1220                  zt = pts   (ji,jj,jk,jp_tem)
1221                  zs = pts   (ji,jj,jk,jp_sal)
1222                  zh = fsdept(ji,jj,jk)        ! depth
1223                  zsr= zws   (ji,jj,jk)        ! square root salinity
1224                  !
1225                  ! 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   &
1227                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
1228                  ! seawater volumic mass atm pressure
1229                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
1230                     &                      -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
1236                  !
1237                  ! 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
1256              ! 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)
1267                  ! 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)
1275               END DO
1276            END DO
1277         END DO
1278         !
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
1365         !
1366         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1367         DO jk = 1, jpkm1
1368            DO jj = 1, jpj
1369               DO ji = 1, jpi
1370                  zt = pts   (ji,jj,jk,jp_tem)
1371                  zs = pts   (ji,jj,jk,jp_sal)
1372                  zh = fsdept(ji,jj,jk)        ! depth
1373                  zsr= zws   (ji,jj,jk)        ! square root salinity
1374                  !
1375                  ! 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   &
1377                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
1378                  ! seawater volumic mass atm pressure
1379                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
1380                     &                      -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
1386                  !
1387                  ! 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                  !
1454               END DO
1455            END DO
1456         END DO
1457         !
1458      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
1459         !
1460         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1461         DO jk = 1, jpkm1
1462            DO jj = 1, jpj
1463               DO ji = 1, jpi
1464                  zt = pts   (ji,jj,jk,jp_tem)
1465                  zs = pts   (ji,jj,jk,jp_sal)
1466                  zh = fsdept(ji,jj,jk)        ! depth
1467                  zsr= zws   (ji,jj,jk)        ! square root salinity
1468                  !
1469                  ! 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   &
1471                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
1472                  ! seawater volumic mass atm pressure
1473                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
1474                     &                      -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
1480                  !
1481                  ! 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
1510                  ! 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)
1547                  !
1548               END DO
1549            END DO
1550         END DO
1551         !
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  ==!
1567         DO jk = 1, jpkm1
1568            DO jj = 1, jpj
1569               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)
1579                  !
1580               END DO
1581            END DO
1582         END DO
1583         pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! betaPE=beta
1584         !
1585      CASE DEFAULT
1586         IF(lwp) WRITE(numout,cform_err)
1587         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
1588         nstop = nstop + 1
1589         !
1590      END SELECT
1591      !
1592      CALL wrk_dealloc( jpi, jpj, jpk, zws )
1593      !
1594      IF( nn_timing == 1 ) CALL timing_stop('eos_pen')
1595      !
1596   END SUBROUTINE eos_pen
1597
1598
1599
1600   FUNCTION tfreez( psal ) RESULT( ptf )
1601      !!----------------------------------------------------------------------
1602      !!                 ***  ROUTINE tfreez  ***
1603      !!
1604      !! ** Purpose :   Compute the sea surface freezing temperature [Celcius]
1605      !!
1606      !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???)
1607      !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
1608      !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
1609      !!
1610      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
1611      !!----------------------------------------------------------------------
1612      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
1613      ! Leave result array automatic rather than making explicitly allocated
1614      REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius]
1615      !!----------------------------------------------------------------------
1616      !
1617      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
1618         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
1619      !
1620   END FUNCTION tfreez
1621
1622
1623   SUBROUTINE eos_init
1624      !!----------------------------------------------------------------------
1625      !!                 ***  ROUTINE eos_init  ***
1626      !!
1627      !! ** Purpose :   initializations for the equation of state
1628      !!
1629      !! ** Method  :   Read the namelist nameos and control the parameters
1630      !!----------------------------------------------------------------------
1631      NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta
1632      !!----------------------------------------------------------------------
1633      !
1634      REWIND( numnam )            ! Read Namelist nameos : equation of state
1635      READ  ( numnam, nameos )
1636      !
1637      IF(lwp) THEN                ! Control print
1638         WRITE(numout,*)
1639         WRITE(numout,*) 'eos_init : equation of state'
1640         WRITE(numout,*) '~~~~~~~~'
1641         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
1645      ENDIF
1646      !
1647      SELECT CASE( nn_eos )         ! check option
1648      !
1649      CASE( -1 )                        !==  Jackett and McDougall (1995) formulation  ==!
1650         IF(lwp) WRITE(numout,*)
1651         IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1995) equation of state'
1652         !
1653      CASE( 0 )                        !==  modified Jackett and McDougall (1995) formulation  ==!
1654         IF(lwp) WRITE(numout,*)
1655         IF(lwp) WRITE(numout,*) '          use of modified Jackett & McDougall (1995) equation of state'
1656         IF(lwp) WRITE(numout,*) '                and McDougall (1987) Brunt-Vaisala frequency'
1657         !
1658      CASE( 1 )                        !==  Linear formulation = F( temperature )  ==!
1659         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
1673         !
1674      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
1675         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
1676         CALL ctl_stop( ctmp1 )
1677         !
1678      END SELECT
1679      !
1680   END SUBROUTINE eos_init
1681
1682   !!======================================================================
1683END MODULE eosbn2
Note: See TracBrowser for help on using the repository browser.