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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
eosbn2.F90 in branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 4619

Last change on this file since 4619 was 4619, checked in by gm, 10 years ago

#1294 : TEOS-10 and Ediag

  • Property svn:keywords set to Id
File size: 63.2 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 (now eos_fzp function)
18   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add alpha/beta used in ldfslp
20   !!            3.7  ! 2012-03  (F. Roquet, G. Madec)  add primitive of alpha and beta used in PE computation
21   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state
22   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module
23   !!             -   ! 2014-02  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80
24   !!----------------------------------------------------------------------
25
26   !!----------------------------------------------------------------------
27   !!   eos            : generic interface of the equation of state
28   !!   eos_insitu     : Compute the in situ density
29   !!   eos_insitu_pot : Compute the insitu and surface referenced potential volumic mass
30   !!   eos_insitu_2d  : Compute the in situ density for 2d fields
31   !!   bn2            : Compute the Brunt-Vaisala frequency
32   !!   eos_rab        : generic interface of in situ thermal/haline expansion ratio
33   !!   eos_rab_3d     : compute in situ thermal/haline expansion ratio
34   !!   eos_rab_2d     : compute in situ thermal/haline expansion ratio for 2d fields
35   !!   eos_fzp        : freezing temperature
36   !!   eos_init       : set eos parameters (namelist)
37   !!----------------------------------------------------------------------
38   USE dom_oce         ! ocean space and time domain
39   USE phycst          ! physical constants
40   !
41   USE in_out_manager  ! I/O manager
42   USE lib_mpp         ! MPP library
43   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
44   USE prtctl          ! Print control
45   USE wrk_nemo        ! Memory Allocation
46   USE timing          ! Timing
47
48   IMPLICIT NONE
49   PRIVATE
50
51   !                   !! * Interface
52   INTERFACE eos
53      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
54   END INTERFACE
55   !
56   INTERFACE eos_rab
57      MODULE PROCEDURE rab_3d, rab_2d
58   END INTERFACE
59   !
60   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules
61   PUBLIC   bn2            ! called by step module
62   PUBLIC   eos_rab        ! called by ldfslp, zdfddm, trabbl
63   PUBLIC   eos_pt_from_ct ! called by sbcssm
64   PUBLIC   eos_fzp        ! called by traadv_cen2 and sbcice_... modules
65   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen module
66   PUBLIC   eos_init       ! called by istate module
67
68   !                                          !!* Namelist (nameos) *
69   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
70   LOGICAL , PUBLIC ::   ln_useCT  = .FALSE.  ! determine if eos_pt_from_ct is used to compute sst_m
71
72   !                                   !!!  simplified eos coefficients
73   ! default value: Vallis 2006
74   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.
75   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.
76   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2       
77   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2       
78   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T 
79   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S 
80   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt 
81   
82   ! EOS parameters
83   REAL(wp) ::   EOS111 , EOS211 , EOS311 , EOS411 , EOS511 , EOS611 , EOS711
84   REAL(wp) ::   EOS121 , EOS221 , EOS321 , EOS421 , EOS521 , EOS621
85   REAL(wp) ::   EOS131 , EOS231 , EOS331 , EOS431 , EOS531
86   REAL(wp) ::   EOS141 , EOS241 , EOS341 , EOS441
87   REAL(wp) ::   EOS151 , EOS251 , EOS351
88   REAL(wp) ::   EOS161 , EOS261
89   REAL(wp) ::   EOS171
90   REAL(wp) ::   EOS112 , EOS212 , EOS312 , EOS412 , EOS512
91   REAL(wp) ::   EOS122 , EOS222 , EOS322 , EOS422
92   REAL(wp) ::   EOS132 , EOS232 , EOS332
93   REAL(wp) ::   EOS142 , EOS242
94   REAL(wp) ::   EOS152
95   REAL(wp) ::   EOS113 , EOS213 , EOS313
96   REAL(wp) ::   EOS123 , EOS223
97   REAL(wp) ::   EOS133
98   
99   ! ALPHA parameters
100   REAL(wp) ::   ALP111 , ALP211 , ALP311 , ALP411 , ALP511 , ALP611
101   REAL(wp) ::   ALP121 , ALP221 , ALP321 , ALP421 , ALP521
102   REAL(wp) ::   ALP131 , ALP231 , ALP331 , ALP431
103   REAL(wp) ::   ALP141 , ALP241 , ALP341
104   REAL(wp) ::   ALP151 , ALP251
105   REAL(wp) ::   ALP161
106   REAL(wp) ::   ALP112 , ALP212 , ALP312 , ALP412
107   REAL(wp) ::   ALP122 , ALP222 , ALP322
108   REAL(wp) ::   ALP132 , ALP232
109   REAL(wp) ::   ALP142
110   REAL(wp) ::   ALP113 , ALP213
111   REAL(wp) ::   ALP123
112   
113   ! BETA parameters
114   REAL(wp) ::   BET111 , BET211 , BET311 , BET411 , BET511 , BET611
115   REAL(wp) ::   BET121 , BET221 , BET321 , BET421 , BET521
116   REAL(wp) ::   BET131 , BET231 , BET331 , BET431
117   REAL(wp) ::   BET141 , BET241 , BET341
118   REAL(wp) ::   BET151 , BET251
119   REAL(wp) ::   BET161
120   REAL(wp) ::   BET112 , BET212 , BET312 , BET412
121   REAL(wp) ::   BET122 , BET222 , BET322
122   REAL(wp) ::   BET132 , BET232
123   REAL(wp) ::   BET142
124   REAL(wp) ::   BET113 , BET213
125   REAL(wp) ::   BET123
126
127   ! PEN parameters
128   REAL(wp) ::   PEN112 , PEN212 , PEN312 , PEN412 , PEN512
129   REAL(wp) ::   PEN122 , PEN222 , PEN322 , PEN422
130   REAL(wp) ::   PEN132 , PEN232 , PEN332
131   REAL(wp) ::   PEN142 , PEN242
132   REAL(wp) ::   PEN152
133   REAL(wp) ::   PEN113 , PEN213 , PEN313
134   REAL(wp) ::   PEN123 , PEN223
135   REAL(wp) ::   PEN133
136   
137   ! ALPHA_PEN parameters
138   REAL(wp) ::   APE112 , APE212 , APE312 , APE412
139   REAL(wp) ::   APE122 , APE222 , APE322
140   REAL(wp) ::   APE132 , APE232
141   REAL(wp) ::   APE142
142   REAL(wp) ::   APE113 , APE213
143   REAL(wp) ::   APE123
144
145   ! BETA_PEN parameters
146   REAL(wp) ::   BPE112 , BPE212 , BPE312 , BPE412
147   REAL(wp) ::   BPE122 , BPE222 , BPE322
148   REAL(wp) ::   BPE132 , BPE232
149   REAL(wp) ::   BPE142
150   REAL(wp) ::   BPE113 , BPE213
151   REAL(wp) ::   BPE123
152
153   !! * Substitutions
154#  include "domzgr_substitute.h90"
155#  include "vectopt_loop_substitute.h90"
156   !!----------------------------------------------------------------------
157   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
158   !! $Id$
159   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
160   !!----------------------------------------------------------------------
161CONTAINS
162
163   SUBROUTINE eos_insitu( pts, prd, pdep )
164      !!----------------------------------------------------------------------
165      !!                   ***  ROUTINE eos_insitu  ***
166      !!
167      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
168      !!       potential temperature and salinity using an equation of state
169      !!       defined through the namelist parameter nn_eos.
170      !!
171      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
172      !!         with   prd    in situ density anomalie     no units
173      !!                t      potential temperature        Celsius
174      !!                s      salinity                     psu
175      !!                z      depth                        meters
176      !!                rho    in situ density              kg/m^3
177      !!                rau0   reference density            kg/m^3
178      !!
179      !!     nn_eos = -1 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
180      !!         Check value: rho = 1028.21583 kg/m^3 for p=3000 dbar, t=3 Celcius, s=35.5 psu
181      !!
182      !!     nn_eos =  0 : Jackett and McDougall (1995) equation of state is used for rho(t,s,z)
183      !!            where the density anomaly of a water parcel at t=4, s=35 moved adiabatically to surface
184      !!            has been removed: rho = rho_JM95 - ( rho_JM95(4,35,z) - rho_JM95(4,35,0) )
185      !!         Check value: rho = 1028.34976 kg/m^3 for p=3000 dbar, t=3 Celcius, s=35.5 psu
186      !!         Check value: rho_JM95 = 1041.83267 kg/m^3 for p=3000 dbar, t=3 Celcius, s=35.5 psu
187      !!
188      !!     nn_eos =  1 : simplified equation of state
189      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
190      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
191      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
192      !!              Vallis like equation: use default values of coefficients
193      !!
194      !! ** Action  :   compute prd , the in situ density (no units)
195      !!
196      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
197      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
198      !!                TEOS-10 Manual, 2010
199      !!                Roquet, Ocean Modelling, in preparation (2013)
200      !!----------------------------------------------------------------------
201      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
202      !                                                               ! 2 : salinity               [psu]
203      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd   ! in situ density            [-]
204      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep  ! depth                      [m]
205      !
206      INTEGER  ::   ji, jj, jk                ! dummy loop indices
207      REAL(wp) ::   zt , zh , zs , zsr, ztm   ! local scalars
208      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
209      !!----------------------------------------------------------------------
210      !
211      IF( nn_timing == 1 )   CALL timing_start('eos')
212      !
213      SELECT CASE( nn_eos )
214      !
215      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
216         !
217         DO jk = 1, jpkm1
218            DO jj = 1, jpj
219               DO ji = 1, jpi
220                  !
221                  zh  = pdep(ji,jj,jk)                                  ! depth
222                  zt  = pts (ji,jj,jk,jp_tem)                           ! temperature
223                  zsr = SQRT( MAX( pts(ji,jj,jk,jp_sal) , 0.1_wp ) )    ! square root salinity
224                  ztm = tmask(ji,jj,jk)                                 ! tmask
225                  !
226                  zn2  = (   EOS133*zsr  + EOS223*zt + EOS123 )*zsr + ( EOS313*zt + EOS213 )*zt + EOS113
227                  !
228                  zn1  = ( ( ( EOS152*zsr   &
229                     &   + EOS242*zt + EOS142 )*zsr   &
230                     &   + ( EOS332*zt + EOS232 )*zt + EOS132 )*zsr   &
231                     &   + ( ( EOS422*zt + EOS322 )*zt + EOS222 )*zt + EOS122 )*zsr   &
232                     &   + ( ( ( EOS512*zt + EOS412 )*zt + EOS312 )*zt + EOS212 )*zt + EOS112
233                  !
234                  zn0  = ( ( ( ( ( EOS171*zsr   &
235                     &   + EOS261*zt + EOS161 )*zsr   &
236                     &   + ( EOS351*zt + EOS251 )*zt + EOS151 )*zsr   &
237                     &   + ( ( EOS441*zt + EOS341 )*zt + EOS241 )*zt + EOS141 )*zsr   &
238                     &   + ( ( ( EOS531*zt + EOS431 )*zt + EOS331 )*zt + EOS231 )*zt + EOS131 )*zsr   &
239                     &   + ( ( ( ( EOS621*zt + EOS521 )*zt + EOS421 )*zt + EOS321 )*zt + EOS221 )*zt + EOS121 )*zsr   &
240                     &   + ( ( ( ( ( EOS711*zt + EOS611 )*zt + EOS511 )*zt + EOS411 )*zt + EOS311 )*zt + EOS211 )*zt + EOS111
241                  !
242                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
243                  !                                               
244                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked)
245               END DO
246            END DO
247         END DO
248         !
249      CASE( 1 )                !==  simplified EOS  ==!
250         !
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  = pdep (ji,jj,jk)
257                  ztm = tmask(ji,jj,jk)
258                  !
259                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
260                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
261                     &  - rn_nu * zt * zs
262                  !                                 
263                  prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
264               END DO
265            END DO
266         END DO
267         !
268      END SELECT
269      !
270      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk )
271      !
272      IF( nn_timing == 1 )   CALL timing_stop('eos')
273      !
274   END SUBROUTINE eos_insitu
275
276
277   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
278      !!----------------------------------------------------------------------
279      !!                  ***  ROUTINE eos_insitu_pot  ***
280      !!
281      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
282      !!      potential volumic mass (Kg/m3) from potential temperature and
283      !!      salinity fields using an equation of state defined through the
284      !!     namelist parameter nn_eos.
285      !!
286      !! ** Action  : - prd  , the in situ density (no units)
287      !!              - prhop, the potential volumic mass (Kg/m3)
288      !!
289      !!----------------------------------------------------------------------
290      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
291      !                                                                ! 2 : salinity               [psu]
292      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
293      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
294      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
295      !
296      INTEGER  ::   ji, jj, jk                ! dummy loop indices
297      REAL(wp) ::   zt , zh , zs , zsr, ztm   ! local scalars
298      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
299      !!----------------------------------------------------------------------
300      !
301      IF( nn_timing == 1 )   CALL timing_start('eos-p')
302      !
303      SELECT CASE ( nn_eos )
304      !
305      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
306         !
307         DO jk = 1, jpkm1
308            DO jj = 1, jpj
309               DO ji = 1, jpi
310                  !
311                  zh  = pdep (ji,jj,jk)                                 ! depth
312                  zt  = pts  (ji,jj,jk,jp_tem)                          ! conservative temperature
313                  zsr = SQRT( MAX( pts(ji,jj,jk,jp_sal) , 0.1_wp ) )    ! square root salinity
314                  ztm = tmask(ji,jj,jk)                                 ! tmask
315                  !
316                  zn2  = (   EOS133*zsr  + EOS223*zt + EOS123 )*zsr + ( EOS313*zt + EOS213 )*zt + EOS113
317                  !
318                  zn1  = ( ( ( EOS152*zsr   &
319                     &   + EOS242*zt + EOS142 )*zsr   &
320                     &   + ( EOS332*zt + EOS232 )*zt + EOS132 )*zsr   &
321                     &   + ( ( EOS422*zt + EOS322 )*zt + EOS222 )*zt + EOS122 )*zsr   &
322                     &   + ( ( ( EOS512*zt + EOS412 )*zt + EOS312 )*zt + EOS212 )*zt + EOS112
323                  !
324                  zn0  = ( ( ( ( ( EOS171*zsr   &
325                     &   + EOS261*zt + EOS161 )*zsr   &
326                     &   + ( EOS351*zt + EOS251 )*zt + EOS151 )*zsr   &
327                     &   + ( ( EOS441*zt + EOS341 )*zt + EOS241 )*zt + EOS141 )*zsr   &
328                     &   + ( ( ( EOS531*zt + EOS431 )*zt + EOS331 )*zt + EOS231 )*zt + EOS131 )*zsr   &
329                     &   + ( ( ( ( EOS621*zt + EOS521 )*zt + EOS421 )*zt + EOS321 )*zt + EOS221 )*zt + EOS121 )*zsr   &
330                     &   + ( ( ( ( ( EOS711*zt + EOS611 )*zt + EOS511 )*zt + EOS411 )*zt + EOS311 )*zt + EOS211 )*zt + EOS111
331                  !
332                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
333                  !
334                  prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
335                  !
336                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
337               END DO
338            END DO
339         END DO
340         !
341      CASE( 1 )                !==  simplified EOS  ==!
342         !
343         DO jk = 1, jpkm1
344            DO jj = 1, jpj
345               DO ji = 1, jpi
346                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
347                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
348                  zh  = pdep (ji,jj,jk)
349                  ztm = tmask(ji,jj,jk)
350                  !                                                     ! potential density referenced at the surface
351                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
352                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
353                     &  - rn_nu * zt * zs
354                  prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
355                  !                                                     ! density anomaly (masked)
356                  zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
357                  prd(ji,jj,jk) = zn * r1_rau0 * ztm
358                  !
359               END DO
360            END DO
361         END DO
362         !
363      END SELECT
364      !
365      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
366      !
367      IF( nn_timing == 1 )   CALL timing_stop('eos-p')
368      !
369   END SUBROUTINE eos_insitu_pot
370
371
372   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
373      !!----------------------------------------------------------------------
374      !!                  ***  ROUTINE eos_insitu_2d  ***
375      !!
376      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
377      !!      potential temperature and salinity using an equation of state
378      !!      defined through the namelist parameter nn_eos. * 2D field case
379      !!
380      !! ** Action  : - prd , the in situ density (no units) (unmasked)
381      !!
382      !!----------------------------------------------------------------------
383      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
384      !                                                           ! 2 : salinity               [psu]
385      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
386      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
387      !
388      INTEGER  ::   ji, jj, jk                ! dummy loop indices
389      REAL(wp) ::   zt , zh , zs , zsr, ztm   ! local scalars
390      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
391      !!----------------------------------------------------------------------
392      !
393      IF( nn_timing == 1 )   CALL timing_start('eos2d')
394      !
395      prd(:,:) = 0._wp
396      !
397      SELECT CASE( nn_eos )
398      !
399      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
400         !
401         DO jj = 1, jpjm1
402            DO ji = 1, fs_jpim1   ! vector opt.
403               !
404               zt    = pts  (ji,jj,jp_tem)            ! interpolated T
405               zsr = SQRT( MAX( pts(ji,jj,jp_sal) , 0.1_wp ) )     ! square root salinity
406               zh    = pdep (ji,jj)            ! depth at the partial step level
407               !
408               !
409               zn2  = (   EOS133*zsr  + EOS223*zt + EOS123 )*zsr + ( EOS313*zt + EOS213 )*zt + EOS113
410               !
411               zn1  = ( ( ( EOS152*zsr   &
412                  &   + EOS242*zt + EOS142 )*zsr   &
413                  &   + ( EOS332*zt + EOS232 )*zt + EOS132 )*zsr   &
414                  &   + ( ( EOS422*zt + EOS322 )*zt + EOS222 )*zt + EOS122 )*zsr   &
415                  &   + ( ( ( EOS512*zt + EOS412 )*zt + EOS312 )*zt + EOS212 )*zt + EOS112
416               !
417               zn0  = ( ( ( ( ( EOS171*zsr   &
418                  &   + EOS261*zt + EOS161 )*zsr   &
419                  &   + ( EOS351*zt + EOS251 )*zt + EOS151 )*zsr   &
420                  &   + ( ( EOS441*zt + EOS341 )*zt + EOS241 )*zt + EOS141 )*zsr   &
421                  &   + ( ( ( EOS531*zt + EOS431 )*zt + EOS331 )*zt + EOS231 )*zt + EOS131 )*zsr   &
422                  &   + ( ( ( ( EOS621*zt + EOS521 )*zt + EOS421 )*zt + EOS321 )*zt + EOS221 )*zt + EOS121 )*zsr   &
423                  &   + ( ( ( ( ( EOS711*zt + EOS611 )*zt + EOS511 )*zt + EOS411 )*zt + EOS311 )*zt + EOS211 )*zt + EOS111
424               !
425               zn   = ( zn2 * zh + zn1 ) * zh + zn0
426               !
427               prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
428               !
429            END DO
430         END DO
431         !
432      CASE( 1 )                !==  simplified EOS  ==!
433         !
434         DO jj = 1, jpjm1
435            DO ji = 1, fs_jpim1   ! vector opt.
436               !
437               zt    = pts  (ji,jj,jp_tem)  - 10._wp
438               zs    = pts  (ji,jj,jp_sal)  - 35._wp
439               zh    = pdep (ji,jj)                         ! depth at the partial step level
440               !
441               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
442                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
443                  &  - rn_nu * zt * zs
444               !
445               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
446               !
447            END DO
448         END DO
449         !
450      END SELECT
451      !
452      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
453      !
454      IF( nn_timing == 1 )   CALL timing_stop('eos2d')
455      !
456   END SUBROUTINE eos_insitu_2d
457
458
459   SUBROUTINE rab_3d( pts, pab )
460      !!----------------------------------------------------------------------
461      !!                 ***  ROUTINE rab_3d  ***
462      !!
463      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points
464      !!
465      !! ** Method  :   calculates alpha / beta at T-points
466      !!       * nn_eos = 0  : polynomial approximation of McDougall (1987)
467      !!                       The alpha/beta is returned as 4-D array pab using the exact expression
468      !!                       based on the JM95 equation of state.
469      !!       * nn_eos = 1  : linear equation of state (temperature only)
470      !!                       We return alpha and beta=0
471      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
472      !!                       We return alpha0 and beta0
473      !!
474      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
475      !!----------------------------------------------------------------------
476      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
477      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
478      !
479      INTEGER  ::   ji, jj, jk                ! dummy loop indices
480      REAL(wp) ::   zt , zh , zs , zsr, ztm   ! local scalars
481      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
482      !!----------------------------------------------------------------------
483      !
484      IF( nn_timing == 1 )   CALL timing_start('rab_3d')
485      !
486      SELECT CASE ( nn_eos )
487      !
488      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
489         !
490         DO jk = 1, jpkm1
491            DO jj = 1, jpj
492               DO ji = 1, jpi
493                  !
494                  zh  = fsdept(ji,jj,jk)                        ! depth
495                  zt  = pts   (ji,jj,jk,jp_tem)                 ! conservative temperature
496                  zsr = SQRT( MAX( pts(ji,jj,jk,jp_sal) , 0.1_wp ) )     ! square root salinity
497                  ztm = tmask(ji,jj,jk)                         ! tmask
498                  !
499                  zn2 = ALP123*zsr + ALP213*zt + ALP113
500                  !
501                  zn1 = ( ( ALP142*zsr   &
502                     &   + ALP232*zt + ALP132 )*zsr   &
503                     &   + ( ALP322*zt + ALP222 )*zt + ALP122 )*zsr   &
504                     &   + ( ( ALP412*zt + ALP312 )*zt + ALP212 )*zt + ALP112
505                  !
506                  zn0 = ( ( ( ( ALP161*zsr   &
507                     &   + ALP251*zt + ALP151 )*zsr   &
508                     &   + ( ALP341*zt + ALP241 )*zt + ALP141 )*zsr   &
509                     &   + ( ( ALP431*zt + ALP331 )*zt + ALP231 )*zt + ALP131 )*zsr   &
510                     &   + ( ( ( ALP521*zt + ALP421 )*zt + ALP321 )*zt + ALP221 )*zt + ALP121 )*zsr   &
511                     &   + ( ( ( ( ALP611*zt + ALP511 )*zt + ALP411 )*zt + ALP311 )*zt + ALP211 )*zt + ALP111
512                     !
513                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
514                  !
515                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
516                  !
517                  !
518                  zn2 = BET123*zsr + BET213*zt + BET113
519                  !
520                  zn1 = ( ( BET142*zsr   &
521                     &   + BET232*zt + BET132 )*zsr   &
522                     &   + ( BET322*zt + BET222 )*zt + BET122 )*zsr   &
523                     &   + ( ( BET412*zt + BET312 )*zt + BET212 )*zt + BET112
524                     !
525                  zn0 = ( ( ( ( BET161*zsr   &
526                     &   + BET251*zt + BET151 )*zsr   &
527                     &   + ( BET341*zt + BET241 )*zt + BET141 )*zsr   &
528                     &   + ( ( BET431*zt + BET331 )*zt + BET231 )*zt + BET131 )*zsr   &
529                     &   + ( ( ( BET521*zt + BET421 )*zt + BET321 )*zt + BET221 )*zt + BET121 )*zsr   &
530                     &   + ( ( ( ( BET611*zt + BET511 )*zt + BET411 )*zt + BET311 )*zt + BET211 )*zt + BET111
531                     !
532                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
533                  !
534                  pab(ji,jj,jk,jp_sal) = zn / zsr * r1_rau0 * ztm
535                  !
536               END DO
537            END DO
538         END DO
539         !
540      CASE( 1 )                  !==  simplified EOS  ==!
541         !
542         DO jk = 1, jpkm1
543            DO jj = 1, jpj
544               DO ji = 1, jpi
545                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
546                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
547                  zh  = fsdept(ji,jj,jk)                 ! depth in meters at t-point
548                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
549                  !                                      ! alpha
550                  pab(ji,jj,jk,jp_tem) = ( rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs ) * r1_rau0 * ztm   ! alpha
551                  pab(ji,jj,jk,jp_sal) = ( rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt ) * r1_rau0 * ztm   ! beta                            ! beta
552               END DO
553            END DO
554         END DO
555         !
556      CASE DEFAULT
557         IF(lwp) WRITE(numout,cform_err)
558         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
559         nstop = nstop + 1
560         !
561      END SELECT
562      !
563      IF( nn_timing == 1 )   CALL timing_stop('rab_3d')
564      !
565   END SUBROUTINE rab_3d
566
567
568   SUBROUTINE rab_2d( pts, pdep, pab )
569      !!----------------------------------------------------------------------
570      !!                 ***  ROUTINE rab_2d  ***
571      !!
572      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
573      !!
574      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
575      !!----------------------------------------------------------------------
576      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
577      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m]
578      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
579      !
580      INTEGER  ::   ji, jj, jk                ! dummy loop indices
581      REAL(wp) ::   zt , zh , zs , zsr, ztm   ! local scalars
582      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
583      !!----------------------------------------------------------------------
584      !
585      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
586      !
587      pab(:,:,:) = 0._wp
588      !
589      SELECT CASE ( nn_eos )
590      !
591      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
592         !
593         DO jj = 1, jpjm1
594            DO ji = 1, fs_jpim1   ! vector opt.
595               !
596               zt    = pts  (ji,jj,jp_tem)            ! interpolated T
597               zsr = SQRT( MAX( pts(ji,jj,jp_sal) , 0.1_wp ) )     ! square root salinity
598               zh    = pdep (ji,jj)                   ! depth at the partial step level
599               !
600               !
601               zn2 = ALP123*zsr + ALP213*zt + ALP113
602               !
603               zn1 = ( ( ALP142*zsr   &
604                  &   + ALP232*zt + ALP132 )*zsr   &
605                  &   + ( ALP322*zt + ALP222 )*zt + ALP122 )*zsr   &
606                  &   + ( ( ALP412*zt + ALP312 )*zt + ALP212 )*zt + ALP112
607               !
608               zn0 = ( ( ( ( ALP161*zsr   &
609                  &   + ALP251*zt + ALP151 )*zsr   &
610                  &   + ( ALP341*zt + ALP241 )*zt + ALP141 )*zsr   &
611                  &   + ( ( ALP431*zt + ALP331 )*zt + ALP231 )*zt + ALP131 )*zsr   &
612                  &   + ( ( ( ALP521*zt + ALP421 )*zt + ALP321 )*zt + ALP221 )*zt + ALP121 )*zsr   &
613                  &   + ( ( ( ( ALP611*zt + ALP511 )*zt + ALP411 )*zt + ALP311 )*zt + ALP211 )*zt + ALP111
614               !
615               zn   = ( zn2 * zh + zn1 ) * zh + zn0
616               !
617               pab(ji,jj,jp_tem) = zn * r1_rau0
618               !
619               !
620               zn2 = BET123*zsr + BET213*zt + BET113
621               !
622               zn1 = ( ( BET142*zsr   &
623                  &   + BET232*zt + BET132 )*zsr   &
624                  &   + ( BET322*zt + BET222 )*zt + BET122 )*zsr   &
625                  &   + ( ( BET412*zt + BET312 )*zt + BET212 )*zt + BET112
626               !
627               zn0 = ( ( ( ( BET161*zsr   &
628                  &   + BET251*zt + BET151 )*zsr   &
629                  &   + ( BET341*zt + BET241 )*zt + BET141 )*zsr   &
630                  &   + ( ( BET431*zt + BET331 )*zt + BET231 )*zt + BET131 )*zsr   &
631                  &   + ( ( ( BET521*zt + BET421 )*zt + BET321 )*zt + BET221 )*zt + BET121 )*zsr   &
632                  &   + ( ( ( ( BET611*zt + BET511 )*zt + BET411 )*zt + BET311 )*zt + BET211 )*zt + BET111
633               !
634               zn   = ( zn2 * zh + zn1 ) * zh + zn0
635               !
636               pab(ji,jj,jp_sal) = zn / zsr * r1_rau0
637               !
638               !
639            END DO
640         END DO
641         !
642      CASE( 1 )                  !==  simplified EOS  ==!
643         !
644         DO jj = 1, jpjm1
645            DO ji = 1, fs_jpim1   ! vector opt.
646               !
647               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
648               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
649               zh    = pdep (ji,jj)                   ! depth at the partial step level
650               !
651               !                                    ! alpha
652               pab(ji,jj,jp_tem) = ( rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs ) * r1_rau0   ! alpha
653               pab(ji,jj,jp_sal) = ( rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt ) * r1_rau0   ! beta                            ! beta
654               !
655            END DO
656         END DO
657         !
658      CASE DEFAULT
659         IF(lwp) WRITE(numout,cform_err)
660         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
661         nstop = nstop + 1
662         !
663      END SELECT
664      !
665      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
666      !
667   END SUBROUTINE rab_2d
668
669
670   SUBROUTINE bn2( pts, pab, pn2 )
671      !!----------------------------------------------------------------------
672      !!                  ***  ROUTINE bn2  ***
673      !!
674      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the
675      !!                time-step of the input arguments
676      !!
677      !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
678      !!      where alpha and beta are given in pab, and computed on T-points.
679      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module.
680      !!
681      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point
682      !!
683      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
684      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
685      !!                McDougall, J. Phys. Oceano., 1987
686      !!----------------------------------------------------------------------
687      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu]
688      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1]
689      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
690      !
691      INTEGER  ::   ji, jj, jk      ! dummy loop indices
692      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
693      !!----------------------------------------------------------------------
694      !
695      IF( nn_timing == 1 ) CALL timing_start('bn2')
696      !
697      DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
698         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
699            DO ji = 1, jpi
700               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
701                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
702               !
703               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
704               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
705               !
706               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
707                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
708                  &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
709            END DO
710         END DO
711      END DO
712      !
713      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
714      !
715      IF( nn_timing == 1 )   CALL timing_stop('bn2')
716      !
717   END SUBROUTINE bn2
718
719
720   FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
721      !!----------------------------------------------------------------------
722      !!                 ***  ROUTINE eos_pt_from_ct  ***
723      !!
724      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celcius]
725      !!
726      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
727      !!       checkvalue: pt=20.0239165517474 Celsius for s=35.7psu, t=20degC
728      !!
729      !! Reference  :   TEOS-10, UNESCO
730      !!                Rational fit, see Roquet (2013)
731      !!----------------------------------------------------------------------
732      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp [Celcius]
733      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity   [psu]
734      ! Leave result array automatic rather than making explicitly allocated
735      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celcius]
736      !
737      INTEGER  ::   ji, jj               ! dummy loop indices
738      REAL(wp) ::   zt , zsr, ztm        ! local scalars
739      REAL(wp) ::   zn , zd              ! local scalars
740      !!----------------------------------------------------------------------
741      !
742      IF ( nn_timing == 1 )   CALL timing_start('eos_pt_from_ct')
743      !
744      DO jj = 1, jpj
745          DO ji = 1, jpi
746              zt  = ctmp   (ji,jj)
747              zsr = SQRT( MAX( psal(ji,jj) , 0.1_wp ) )              ! square root salinity
748              ztm = tmask(ji,jj,1)                          ! tmask
749              !
750              zn = ( ( ( ( -1.960202285944569e-04_wp*zsr    &
751                  &    - 6.698743070819174e-05_wp*zt + 4.456449354084589e-03_wp )*zsr    &
752                  &    + ( 1.424140773841990e-05_wp*zt - 1.961306826357777e-04_wp )*zt &
753                  &        - 1.535458707717669e-02_wp )*zsr    &
754                  &    + ( ( 2.814957335530553e-06_wp*zt + 1.534287192323960e-04_wp )*zt &
755                  &        + 2.701363554476190e-02_wp )*zt - 1.869037589289360e-02_wp )*zsr    &
756                  &    + ( ( ( -1.646758381788031e-08_wp*zt + 4.568111226988701e-07_wp )*zt &
757                  &        - 7.853794523448087e-04_wp )*zt - 8.000001510387515e-03_wp )*zt &
758                  &        - 1.132820958753811e-03_wp )*zsr    &
759                  &    + ( ( ( ( -3.885399051933490e-09_wp*zt + 6.608853649947305e-07_wp )*zt &
760                  &        - 1.507746029483144e-04_wp )*zt - 3.457378558412451e-03_wp )*zt &
761                  &        - 7.616362942562558e-01_wp )*zt - 2.067472332529282e-01_wp
762              !
763              zd = ( ( 8.421053090759674e-04_wp*zsr    &
764                  &    - 8.267668312138185e-04_wp*zt - 6.920238162292441e-02_wp )*zsr    &
765                  &    + ( 5.063023899703211e-05_wp*zt + 1.488626985392304e-02_wp )*zt &
766                  &        + 1.541411683853927e-01_wp )*zsr    &
767                  &    + ( ( 6.537389208598346e-07_wp*zt + 1.856814792031252e-03_wp )*zt &
768                  &        + 1.592477436776734e-01_wp )*zt + 1.406742653394891e+01_wp
769              !
770              ptmp(ji,jj) = ( zt + zn / zd ) * ztm
771              !
772          END DO
773      END DO
774      !
775      IF( nn_timing == 1 )   CALL timing_stop('eos_pt_from_ct')
776      !
777   END FUNCTION eos_pt_from_ct
778
779
780   FUNCTION eos_fzp( psal, pdep ) RESULT( ptf )
781      !!----------------------------------------------------------------------
782      !!                 ***  ROUTINE eos_fzp  ***
783      !!
784      !! ** Purpose :   Compute the freezing point temperature [Celcius]
785      !!
786      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by
787      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
788      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
789      !!
790      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
791      !!----------------------------------------------------------------------
792      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
793      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
794      ! Leave result array automatic rather than making explicitly allocated
795      REAL(wp), DIMENSION(jpi,jpj) ::   ptf   ! freezing temperature [Celcius]
796      !
797      INTEGER  ::   ji, jj               ! dummy loop indices
798      REAL(wp) ::   zt , zsr             ! local scalars
799      !!----------------------------------------------------------------------
800      !
801      SELECT CASE ( nn_eos )
802      !
803      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
804         !
805         DO jj = 1, jpj
806            DO ji = 1, jpi
807               zsr= SQRT( ABS( psal(ji,jj) ) ) * 0.1_wp   ! square root salinity x 0.1
808               ptf(ji,jj) = 0.017947064327968_wp &
809                    &    + zsr*zsr*(-6.07609909992982_wp + zsr*(4.883198653547851_wp &
810                    &    + zsr*(-11.8808160123054_wp + zsr*(13.34658511480257_wp &
811                    &    + zsr*(-8.722761043208607_wp + 2.082038908808201_wp*zsr)))))
812            END DO
813         END DO
814         !
815         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.69e-4 * pdep(:,:)
816         !
817      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
818         !
819         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
820            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
821         !
822         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
823         !
824      CASE DEFAULT
825         IF(lwp) WRITE(numout,cform_err)
826         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
827         nstop = nstop + 1
828         !
829      END SELECT
830      !
831   END FUNCTION eos_fzp
832
833
834   SUBROUTINE eos_pen( pts, pab_pe, ppen )
835      !!----------------------------------------------------------------------
836      !!                 ***  ROUTINE eos_pen  ***
837      !!
838      !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
839      !!
840      !! ** Method  :   PE is defined analytically as the vertical
841      !!                   primitive of EOS times -g integrated between 0 and z>0.
842      !!                pen is the nonlinear PE anomaly: ppen = ( PE - rau0 gz ) / rau0 gz - rd
843      !!                                                      = 1/z * /int_0^z rd dz - rd
844      !!                                where rd is the density anomaly (see eos_rhd function)
845      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
846      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
847      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS
848      !!
849      !!       * nn_eos = -1 : polynomial TEOS-10
850      !!       * nn_eos =  0 : Jackett and McDougall (1995) (should not be used, using polynomial TEOS-10 formulation)
851      !!       * nn_eos =  1 : Vallis equation of state (Vallis 2006)
852      !!
853      !! ** Action  : - pen         : PE anomaly given at T-points
854      !!            : - pab_pe  : given at T-points
855      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe
856      !!                    pab_pe(:,:,:,jp_sal) is beta_pe
857      !!----------------------------------------------------------------------
858      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
859      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
860      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   ppen     ! potential energy anomaly
861      !
862      INTEGER  ::   ji, jj, jk               ! dummy loop indices
863      REAL(wp) ::   zt , zh , zsr, ztm , zs  ! local scalars
864      REAL(wp) ::   zn
865      !!----------------------------------------------------------------------
866      !
867      IF ( nn_timing == 1 )   CALL timing_start('eos_pen')
868      !
869      SELECT CASE ( nn_eos )
870      !
871      CASE ( -1 , 0 )               ! polyTEOS10 (used also for JM95 case as an approximation)
872         !
873         DO jk = 1, jpkm1
874            DO jj = 1, jpj
875               DO ji = 1, jpi
876                  zt  = pts (ji,jj,jk,jp_tem)
877                  zh  = fsdept(ji,jj,jk)                              ! depth
878                  zsr = SQRT( MAX( pts(ji,jj,jk,jp_sal) , 0.1_wp ) )  ! square root salinity
879                  ztm = tmask(ji,jj,jk)                               ! tmask
880                  !
881                  !
882                  zn  = ( ( ( PEN133*zsr   &
883                     &   + PEN223*zt + PEN123 )*zsr   &
884                     &   + ( PEN313*zt + PEN213 )*zt + PEN113 )*zh   &
885                     &   + ( ( ( PEN152*zsr   &
886                     &   + PEN242*zt + PEN142 )*zsr   &
887                     &   + ( PEN332*zt + PEN232 )*zt + PEN132 )*zsr   &
888                     &   + ( ( PEN422*zt + PEN322 )*zt + PEN222 )*zt + PEN122 )*zsr   &
889                     &   + ( ( ( PEN512*zt + PEN412 )*zt + PEN312 )*zt + PEN212 )*zt + PEN112 )*zh
890                  !
891                  !                              ! potential energy non-linear anomaly
892                  ppen(ji,jj,jk)  = zn * r1_rau0 * ztm
893                  !
894                  !
895                  zn  = ( ( ( APE123 )*zsr   &
896                     &   + APE213*zt + APE113 )*zh   &
897                     &   + ( ( ( APE142 )*zsr   &
898                     &   + APE232*zt + APE132 )*zsr   &
899                     &   + ( APE322*zt + APE222 )*zt + APE122 )*zsr   &
900                     &   + ( ( APE412*zt + APE312 )*zt + APE212 )*zt + APE112 )*zh
901                  !
902                  !                              ! alphaPE non-linear anomaly
903                  pab_pe(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
904                  !
905                  !
906                  zn  = ( ( ( BPE123 )*zsr   &
907                     &   + BPE213*zt + BPE113 )*zh   &
908                     &   + ( ( ( BPE142 )*zsr   &
909                     &   + BPE232*zt + BPE132 )*zsr   &
910                     &   + ( BPE322*zt + BPE222 )*zt + BPE122 )*zsr   &
911                     &   + ( ( BPE412*zt + BPE312 )*zt + BPE212 )*zt + BPE112 )*zh
912                  !
913                  !                              ! betaPE non-linear anomaly
914                  pab_pe(ji,jj,jk,jp_sal) = zn / zsr * r1_rau0 * ztm
915                  !
916               END DO
917            END DO
918         END DO
919         !
920      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
921         !
922         DO jk = 1, jpkm1
923            DO jj = 1, jpj
924               DO ji = 1, jpi
925                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
926                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
927                  zh  = fsdept(ji,jj,jk)               ! depth in meters  at t-point
928                  ztm = tmask(ji,jj,jk)                ! tmask
929                  zn  = 0.5_wp * zh * r1_rau0 * ztm
930                  !                                    ! Potential Energy
931                  ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
932                  !                                    ! alphaPE
933                  pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
934                  pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
935                  !
936               END DO
937            END DO
938         END DO
939         !
940      CASE DEFAULT
941         IF(lwp) WRITE(numout,cform_err)
942         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
943         nstop = nstop + 1
944         !
945      END SELECT
946      !
947      IF( nn_timing == 1 )   CALL timing_stop('eos_pen')
948      !
949   END SUBROUTINE eos_pen
950
951
952   SUBROUTINE eos_init
953      !!----------------------------------------------------------------------
954      !!                 ***  ROUTINE eos_init  ***
955      !!
956      !! ** Purpose :   initializations for the equation of state
957      !!
958      !! ** Method  :   Read the namelist nameos and control the parameters
959      !!----------------------------------------------------------------------
960      INTEGER  ::   ios   ! local integer
961      !!
962      NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
963         &                                             rn_lambda2, rn_mu2, rn_nu
964      !!----------------------------------------------------------------------
965      !
966      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
967      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
968901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
969      !
970      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
971      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
972902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
973      WRITE( numond, nameos )
974      !
975      rau0        = 1026._wp      !: volumic mass of reference     [kg/m3]
976      rcp         = 3992._wp      !: heat capacity                 [J/K]
977      !
978      IF(lwp) THEN                ! Control print
979         WRITE(numout,*)
980         WRITE(numout,*) 'eos_init : equation of state'
981         WRITE(numout,*) '~~~~~~~~'
982         WRITE(numout,*) '          Namelist nameos : set eos parameters'
983         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
984         IF( ln_useCT )   THEN
985            WRITE(numout,*) '             model uses Conservative Temperature'
986            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
987         ENDIF
988      ENDIF
989      !
990      SELECT CASE( nn_eos )         ! check option
991      !
992      CASE( -1 )                       !==  polynomial TEOS-10  ==!
993         IF(lwp) WRITE(numout,*)
994         IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
995         !
996         EOS111 =  9.9984112007e+02_wp
997         EOS211 =  6.4904564681e-02_wp
998         EOS311 = -8.1866708294e-03_wp
999         EOS411 =  8.4327846386e-05_wp
1000         EOS511 = -9.9853850905e-07_wp
1001         EOS611 =  9.5342506671e-09_wp
1002         EOS711 = -5.5942725501e-11_wp
1003         EOS121 =  4.8745609554e-04_wp
1004         EOS221 = -7.9137653036e-04_wp
1005         EOS321 =  1.1241300058e-04_wp
1006         EOS421 = -4.2658318975e-06_wp
1007         EOS521 =  2.7589924936e-08_wp
1008         EOS621 =  5.6573234128e-10_wp
1009         EOS131 =  8.2409990767e-01_wp
1010         EOS231 = -3.6746935996e-03_wp
1011         EOS331 =  2.8088034997e-05_wp
1012         EOS431 =  3.9432934562e-07_wp
1013         EOS531 = -8.6081969467e-09_wp
1014         EOS141 = -8.5860758561e-03_wp
1015         EOS241 =  1.0048170774e-04_wp
1016         EOS341 = -1.1219591564e-06_wp
1017         EOS441 = -1.4550322467e-09_wp
1018         EOS151 =  1.2633903654e-03_wp
1019         EOS251 = -7.2707200365e-06_wp
1020         EOS351 =  2.6241351559e-08_wp
1021         EOS161 = -7.4242043203e-05_wp
1022         EOS261 =  3.8374457797e-07_wp
1023         EOS171 =  1.2722088050e-06_wp
1024         EOS112 =  4.3845092294e-04_wp
1025         EOS212 = -3.5237528758e-05_wp
1026         EOS312 =  5.7586626030e-07_wp
1027         EOS412 = -6.3139238271e-09_wp
1028         EOS512 =  3.9434377793e-11_wp
1029         EOS122 = -7.6672754907e-06_wp
1030         EOS222 =  9.0194877799e-08_wp
1031         EOS322 = -1.6940551903e-09_wp
1032         EOS422 =  3.1783762285e-11_wp
1033         EOS132 = -5.9870305749e-06_wp
1034         EOS232 =  1.0325176608e-07_wp
1035         EOS332 = -8.0185331501e-10_wp
1036         EOS142 = -5.6324524012e-07_wp
1037         EOS242 =  6.6221210305e-10_wp
1038         EOS152 =  3.9558159454e-08_wp
1039         EOS113 = -1.6847693220e-09_wp
1040         EOS213 =  5.3653075114e-10_wp
1041         EOS313 = -1.0864160933e-11_wp
1042         EOS123 = -1.6883150438e-09_wp
1043         EOS223 =  6.5982407189e-12_wp
1044         EOS133 =  2.7370939782e-10_wp
1045         !
1046         ALP111 = -6.4904564681e-02_wp
1047         ALP211 =  1.6373341659e-02_wp
1048         ALP311 = -2.5298353916e-04_wp
1049         ALP411 =  3.9941540362e-06_wp
1050         ALP511 = -4.7671253336e-08_wp
1051         ALP611 =  3.3565635301e-10_wp
1052         ALP121 =  7.9137653036e-04_wp
1053         ALP221 = -2.2482600117e-04_wp
1054         ALP321 =  1.2797495693e-05_wp
1055         ALP421 = -1.1035969974e-07_wp
1056         ALP521 = -2.8286617064e-09_wp
1057         ALP131 =  3.6746935996e-03_wp
1058         ALP231 = -5.6176069994e-05_wp
1059         ALP331 = -1.1829880369e-06_wp
1060         ALP431 =  3.4432787787e-08_wp
1061         ALP141 = -1.0048170774e-04_wp
1062         ALP241 =  2.2439183129e-06_wp
1063         ALP341 =  4.3650967401e-09_wp
1064         ALP151 =  7.2707200365e-06_wp
1065         ALP251 = -5.2482703118e-08_wp
1066         ALP161 = -3.8374457797e-07_wp
1067         ALP112 =  3.5237528758e-05_wp
1068         ALP212 = -1.1517325206e-06_wp
1069         ALP312 =  1.8941771481e-08_wp
1070         ALP412 = -1.5773751117e-10_wp
1071         ALP122 = -9.0194877799e-08_wp
1072         ALP222 =  3.3881103805e-09_wp
1073         ALP322 = -9.5351286854e-11_wp
1074         ALP132 = -1.0325176608e-07_wp
1075         ALP232 =  1.6037066300e-09_wp
1076         ALP142 = -6.6221210305e-10_wp
1077         ALP113 = -5.3653075114e-10_wp
1078         ALP213 =  2.1728321866e-11_wp
1079         ALP123 = -6.5982407189e-12_wp
1080         !
1081         BET111 =  2.4372804777e-04_wp
1082         BET211 = -3.9568826518e-04_wp
1083         BET311 =  5.6206500292e-05_wp
1084         BET411 = -2.1329159488e-06_wp
1085         BET511 =  1.3794962468e-08_wp
1086         BET611 =  2.8286617064e-10_wp
1087         BET121 =  8.2409990767e-01_wp
1088         BET221 = -3.6746935996e-03_wp
1089         BET321 =  2.8088034997e-05_wp
1090         BET421 =  3.9432934562e-07_wp
1091         BET521 = -8.6081969467e-09_wp
1092         BET131 = -1.2879113784e-02_wp
1093         BET231 =  1.5072256162e-04_wp
1094         BET331 = -1.6829387347e-06_wp
1095         BET431 = -2.1825483700e-09_wp
1096         BET141 =  2.5267807308e-03_wp
1097         BET241 = -1.4541440073e-05_wp
1098         BET341 =  5.2482703118e-08_wp
1099         BET151 = -1.8560510801e-04_wp
1100         BET251 =  9.5936144494e-07_wp
1101         BET161 =  3.8166264151e-06_wp
1102         BET112 = -3.8336377454e-06_wp
1103         BET212 =  4.5097438900e-08_wp
1104         BET312 = -8.4702759513e-10_wp
1105         BET412 =  1.5891881142e-11_wp
1106         BET122 = -5.9870305749e-06_wp
1107         BET222 =  1.0325176608e-07_wp
1108         BET322 = -8.0185331501e-10_wp
1109         BET132 = -8.4486786018e-07_wp
1110         BET232 =  9.9331815458e-10_wp
1111         BET142 =  7.9116318907e-08_wp
1112         BET113 = -8.4415752190e-10_wp
1113         BET213 =  3.2991203594e-12_wp
1114         BET123 =  2.7370939782e-10_wp
1115         !
1116         PEN112 = -2.1922546147e-04_wp
1117         PEN212 =  1.7618764379e-05_wp
1118         PEN312 = -2.8793313015e-07_wp
1119         PEN412 =  3.1569619136e-09_wp
1120         PEN512 = -1.9717188896e-11_wp
1121         PEN122 =  3.8336377454e-06_wp
1122         PEN222 = -4.5097438900e-08_wp
1123         PEN322 =  8.4702759513e-10_wp
1124         PEN422 = -1.5891881142e-11_wp
1125         PEN132 =  2.9935152875e-06_wp
1126         PEN232 = -5.1625883042e-08_wp
1127         PEN332 =  4.0092665751e-10_wp
1128         PEN142 =  2.8162262006e-07_wp
1129         PEN242 = -3.3110605153e-10_wp
1130         PEN152 = -1.9779079727e-08_wp
1131         PEN113 =  1.1231795480e-09_wp
1132         PEN213 = -3.5768716743e-10_wp
1133         PEN313 =  7.2427739554e-12_wp
1134         PEN123 =  1.1255433625e-09_wp
1135         PEN223 = -4.3988271459e-12_wp
1136         PEN133 = -1.8247293188e-10_wp
1137         !
1138         APE112 = -1.7618764379e-05_wp
1139         APE212 =  5.7586626030e-07_wp
1140         APE312 = -9.4708857407e-09_wp
1141         APE412 =  7.8868755586e-11_wp
1142         APE122 =  4.5097438900e-08_wp
1143         APE222 = -1.6940551903e-09_wp
1144         APE322 =  4.7675643427e-11_wp
1145         APE132 =  5.1625883042e-08_wp
1146         APE232 = -8.0185331501e-10_wp
1147         APE142 =  3.3110605153e-10_wp
1148         APE113 =  3.5768716743e-10_wp
1149         APE213 = -1.4485547911e-11_wp
1150         APE123 =  4.3988271459e-12_wp
1151         !
1152         BPE112 =  1.9168188727e-06_wp
1153         BPE212 = -2.2548719450e-08_wp
1154         BPE312 =  4.2351379756e-10_wp
1155         BPE412 = -7.9459405711e-12_wp
1156         BPE122 =  2.9935152875e-06_wp
1157         BPE222 = -5.1625883042e-08_wp
1158         BPE322 =  4.0092665751e-10_wp
1159         BPE132 =  4.2243393009e-07_wp
1160         BPE232 = -4.9665907729e-10_wp
1161         BPE142 = -3.9558159454e-08_wp
1162         BPE113 =  5.6277168127e-10_wp
1163         BPE213 = -2.1994135730e-12_wp
1164         BPE123 = -1.8247293188e-10_wp
1165         !
1166      CASE( 0 )                        !==  polynomial EOS-80 formulation  ==!
1167         !
1168         IF(lwp) WRITE(numout,*)
1169         IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)'
1170         !
1171         EOS111 =  9.9984240021e+02_wp
1172         EOS211 =  6.8051135248e-02_wp
1173         EOS311 = -9.0991456041e-03_wp
1174         EOS411 =  9.9238754931e-05_wp
1175         EOS511 = -1.0602210389e-06_wp
1176         EOS611 =  5.1296071368e-09_wp
1177         EOS711 =  1.1887242136e-11_wp
1178         EOS121 = -1.9799119661e-04_wp
1179         EOS221 = -1.6676932531e-05_wp
1180         EOS321 =  3.2084304200e-06_wp
1181         EOS421 = -4.8983568948e-08_wp
1182         EOS521 = -1.5610111448e-09_wp
1183         EOS621 =  2.8366680541e-11_wp
1184         EOS131 =  8.2448582544e-01_wp
1185         EOS231 = -4.0764319966e-03_wp
1186         EOS331 =  7.4898579574e-05_wp
1187         EOS431 = -7.7884546946e-07_wp
1188         EOS531 =  5.0405569756e-09_wp
1189         EOS141 = -5.7150204345e-03_wp
1190         EOS241 =  1.0157923923e-04_wp
1191         EOS341 = -1.6325795501e-06_wp
1192         EOS441 = -1.6144820180e-09_wp
1193         EOS151 =  4.8005813592e-04_wp
1194         EOS251 =  1.6199868493e-07_wp
1195         EOS351 =  7.8805435082e-09_wp
1196         EOS161 =  1.9978872493e-07_wp
1197         EOS261 = -3.1437087783e-08_wp
1198         EOS171 =  2.5463113215e-08_wp
1199         EOS112 =  4.3756680685e-04_wp
1200         EOS212 = -3.7298355935e-05_wp
1201         EOS312 =  6.5348660941e-07_wp
1202         EOS412 = -7.2558835910e-09_wp
1203         EOS512 =  3.8273323221e-11_wp
1204         EOS122 =  1.4369434305e-06_wp
1205         EOS222 =  5.6521871685e-08_wp
1206         EOS322 = -1.2136936896e-08_wp
1207         EOS422 =  2.1030996727e-10_wp
1208         EOS132 = -1.0262750144e-05_wp
1209         EOS232 =  2.0858571138e-07_wp
1210         EOS332 = -1.4826357415e-09_wp
1211         EOS142 =  8.5860757162e-08_wp
1212         EOS242 = -5.2140945761e-09_wp
1213         EOS152 =  8.2162682979e-09_wp
1214         EOS113 = -5.5011801985e-09_wp
1215         EOS213 =  5.7020541017e-10_wp
1216         EOS313 = -1.0193603551e-11_wp
1217         EOS123 = -1.0893703284e-10_wp
1218         EOS223 = -4.1944954508e-12_wp
1219         EOS133 =  1.1857534169e-10_wp
1220         !
1221         ALP111 = -6.8051135248e-02_wp
1222         ALP211 =  1.8198291208e-02_wp
1223         ALP311 = -2.9771626479e-04_wp
1224         ALP411 =  4.2408841556e-06_wp
1225         ALP511 = -2.5648035684e-08_wp
1226         ALP611 = -7.1323452819e-11_wp
1227         ALP121 =  1.6676932531e-05_wp
1228         ALP221 = -6.4168608399e-06_wp
1229         ALP321 =  1.4695070684e-07_wp
1230         ALP421 =  6.2440445791e-09_wp
1231         ALP521 = -1.4183340270e-10_wp
1232         ALP131 =  4.0764319966e-03_wp
1233         ALP231 = -1.4979715915e-04_wp
1234         ALP331 =  2.3365364084e-06_wp
1235         ALP431 = -2.0162227903e-08_wp
1236         ALP141 = -1.0157923923e-04_wp
1237         ALP241 =  3.2651591003e-06_wp
1238         ALP341 =  4.8434460540e-09_wp
1239         ALP151 = -1.6199868493e-07_wp
1240         ALP251 = -1.5761087016e-08_wp
1241         ALP161 =  3.1437087783e-08_wp
1242         ALP112 =  3.7298355935e-05_wp
1243         ALP212 = -1.3069732188e-06_wp
1244         ALP312 =  2.1767650773e-08_wp
1245         ALP412 = -1.5309329288e-10_wp
1246         ALP122 = -5.6521871685e-08_wp
1247         ALP222 =  2.4273873792e-08_wp
1248         ALP322 = -6.3092990180e-10_wp
1249         ALP132 = -2.0858571138e-07_wp
1250         ALP232 =  2.9652714830e-09_wp
1251         ALP142 =  5.2140945761e-09_wp
1252         ALP113 = -5.7020541017e-10_wp
1253         ALP213 =  2.0387207103e-11_wp
1254         ALP123 =  4.1944954508e-12_wp
1255         !
1256         BET111 = -9.8995598303e-05_wp
1257         BET211 = -8.3384662657e-06_wp
1258         BET311 =  1.6042152100e-06_wp
1259         BET411 = -2.4491784474e-08_wp
1260         BET511 = -7.8050557238e-10_wp
1261         BET611 =  1.4183340270e-11_wp
1262         BET121 =  8.2448582544e-01_wp
1263         BET221 = -4.0764319966e-03_wp
1264         BET321 =  7.4898579574e-05_wp
1265         BET421 = -7.7884546946e-07_wp
1266         BET521 =  5.0405569756e-09_wp
1267         BET131 = -8.5725306517e-03_wp
1268         BET231 =  1.5236885884e-04_wp
1269         BET331 = -2.4488693252e-06_wp
1270         BET431 = -2.4217230270e-09_wp
1271         BET141 =  9.6011627184e-04_wp
1272         BET241 =  3.2399736986e-07_wp
1273         BET341 =  1.5761087016e-08_wp
1274         BET151 =  4.9947181231e-07_wp
1275         BET251 = -7.8592719459e-08_wp
1276         BET161 =  7.6389339644e-08_wp
1277         BET112 =  7.1847171523e-07_wp
1278         BET212 =  2.8260935842e-08_wp
1279         BET312 = -6.0684684481e-09_wp
1280         BET412 =  1.0515498363e-10_wp
1281         BET122 = -1.0262750144e-05_wp
1282         BET222 =  2.0858571138e-07_wp
1283         BET322 = -1.4826357415e-09_wp
1284         BET132 =  1.2879113574e-07_wp
1285         BET232 = -7.8211418642e-09_wp
1286         BET142 =  1.6432536596e-08_wp
1287         BET113 = -5.4468516419e-11_wp
1288         BET213 = -2.0972477254e-12_wp
1289         BET123 =  1.1857534169e-10_wp
1290         !
1291         PEN112 = -2.1878340343e-04_wp
1292         PEN212 =  1.8649177967e-05_wp
1293         PEN312 = -3.2674330470e-07_wp
1294         PEN412 =  3.6279417955e-09_wp
1295         PEN512 = -1.9136661610e-11_wp
1296         PEN122 = -7.1847171523e-07_wp
1297         PEN222 = -2.8260935842e-08_wp
1298         PEN322 =  6.0684684481e-09_wp
1299         PEN422 = -1.0515498363e-10_wp
1300         PEN132 =  5.1313750720e-06_wp
1301         PEN232 = -1.0429285569e-07_wp
1302         PEN332 =  7.4131787075e-10_wp
1303         PEN142 = -4.2930378581e-08_wp
1304         PEN242 =  2.6070472881e-09_wp
1305         PEN152 = -4.1081341490e-09_wp
1306         PEN113 =  3.6674534657e-09_wp
1307         PEN213 = -3.8013694011e-10_wp
1308         PEN313 =  6.7957357010e-12_wp
1309         PEN123 =  7.2624688558e-11_wp
1310         PEN223 =  2.7963303006e-12_wp
1311         PEN133 = -7.9050227794e-11_wp
1312         !
1313         APE112 = -1.8649177967e-05_wp
1314         APE212 =  6.5348660941e-07_wp
1315         APE312 = -1.0883825386e-08_wp
1316         APE412 =  7.6546646441e-11_wp
1317         APE122 =  2.8260935842e-08_wp
1318         APE222 = -1.2136936896e-08_wp
1319         APE322 =  3.1546495090e-10_wp
1320         APE132 =  1.0429285569e-07_wp
1321         APE232 = -1.4826357415e-09_wp
1322         APE142 = -2.6070472881e-09_wp
1323         APE113 =  3.8013694011e-10_wp
1324         APE213 = -1.3591471402e-11_wp
1325         APE123 = -2.7963303006e-12_wp
1326         !
1327         BPE112 = -3.5923585761e-07_wp
1328         BPE212 = -1.4130467921e-08_wp
1329         BPE312 =  3.0342342240e-09_wp
1330         BPE412 = -5.2577491816e-11_wp
1331         BPE122 =  5.1313750720e-06_wp
1332         BPE222 = -1.0429285569e-07_wp
1333         BPE322 =  7.4131787075e-10_wp
1334         BPE132 = -6.4395567871e-08_wp
1335         BPE232 =  3.9105709321e-09_wp
1336         BPE142 = -8.2162682979e-09_wp
1337         BPE113 =  3.6312344279e-11_wp
1338         BPE213 =  1.3981651503e-12_wp
1339         BPE123 = -7.9050227794e-11_wp
1340         !
1341      CASE( 1 )                        !==  Simplified EOS     ==!
1342         IF(lwp) THEN
1343            WRITE(numout,*)
1344            WRITE(numout,*) '          use of simplified eos:    rhd(dT=T-10,dS=S-35,Z) = '
1345            WRITE(numout,*) '             [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0'
1346            WRITE(numout,*)
1347            WRITE(numout,*) '             thermal exp. coef.    rn_a0      = ', rn_a0
1348            WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0
1349            WRITE(numout,*) '             cabbeling coef.       rn_lambda1 = ', rn_lambda1
1350            WRITE(numout,*) '             cabbeling coef.       rn_lambda2 = ', rn_lambda2
1351            WRITE(numout,*) '             thermobar. coef.      rn_mu1     = ', rn_mu1
1352            WRITE(numout,*) '             thermobar. coef.      rn_mu2     = ', rn_mu2
1353            WRITE(numout,*) '             2nd cabbel. coef.     rn_nu      = ', rn_nu
1354            WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization '
1355         ENDIF
1356         !
1357      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
1358         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
1359         CALL ctl_stop( ctmp1 )
1360         !
1361      END SELECT
1362      !
1363      r1_rau0     = 1._wp / rau0
1364      r1_rcp      = 1._wp / rcp
1365      r1_rau0_rcp = 1._wp / ( rau0 * rcp )
1366      !
1367      IF(lwp) WRITE(numout,*)
1368      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
1369      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
1370      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
1371      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
1372      !
1373   END SUBROUTINE eos_init
1374
1375   !!======================================================================
1376END MODULE eosbn2
Note: See TracBrowser for help on using the repository browser.