source: NEMO/trunk/src/OCE/TRA/eosbn2.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 10 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 71.8 KB
Line 
1MODULE eosbn2
2   !!==============================================================================
3   !!                       ***  MODULE  eosbn2  ***
4   !! Equation Of Seawater : in situ density - Brunt-Vaisala frequency
5   !!==============================================================================
6   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
7   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
8   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
9   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
10   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
11   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
12   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
13   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
14   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
15   !!             -   ! 2003-08  (G. Madec)  F90, free form
16   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function (now eos_fzp function)
17   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
18   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add alpha/beta used in ldfslp
19   !!            3.7  ! 2012-03  (F. Roquet, G. Madec)  add primitive of alpha and beta used in PE computation
20   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state
21   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module
22   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80
23   !!             -   ! 2015-06  (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF
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_pt_from_ct: compute the potential temperature from the Conservative Temperature
33   !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio
34   !!   eos_rab_3d    : compute in situ thermal/haline expansion ratio
35   !!   eos_rab_2d    : compute in situ thermal/haline expansion ratio for 2d fields
36   !!   eos_fzp_2d    : freezing temperature for 2d fields
37   !!   eos_fzp_0d    : freezing temperature for scalar
38   !!   eos_init      : set eos parameters (namelist)
39   !!----------------------------------------------------------------------
40   USE dom_oce        ! ocean space and time domain
41   USE phycst         ! physical constants
42   USE stopar         ! Stochastic T/S fluctuations
43   USE stopts         ! Stochastic T/S fluctuations
44   !
45   USE in_out_manager ! I/O manager
46   USE lib_mpp        ! MPP library
47   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
48   USE prtctl         ! Print control
49   USE lbclnk         ! ocean lateral boundary conditions
50   USE timing         ! Timing
51
52   IMPLICIT NONE
53   PRIVATE
54
55   !                  !! * Interface
56   INTERFACE eos
57      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
58   END INTERFACE
59   !
60   INTERFACE eos_rab
61      MODULE PROCEDURE rab_3d, rab_2d, rab_0d
62   END INTERFACE
63   !
64   INTERFACE eos_fzp 
65      MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
66   END INTERFACE
67   !
68   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules
69   PUBLIC   bn2            ! called by step module
70   PUBLIC   eos_rab        ! called by ldfslp, zdfddm, trabbl
71   PUBLIC   eos_pt_from_ct ! called by sbcssm
72   PUBLIC   eos_fzp        ! called by traadv_cen2 and sbcice_... modules
73   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen module
74   PUBLIC   eos_init       ! called by istate module
75
76   !                               !!** Namelist nameos **
77   LOGICAL , PUBLIC ::   ln_TEOS10
78   LOGICAL , PUBLIC ::   ln_EOS80
79   LOGICAL , PUBLIC ::   ln_SEOS
80
81   ! Parameters
82   LOGICAL , PUBLIC    ::   l_useCT         ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise
83   INTEGER , PUBLIC    ::   neos            ! Identifier for equation of state used
84
85   INTEGER , PARAMETER ::   np_teos10 = -1  ! parameter for using TEOS10
86   INTEGER , PARAMETER ::   np_eos80  =  0  ! parameter for using EOS80
87   INTEGER , PARAMETER ::   np_seos   = 1   ! parameter for using Simplified Equation of state
88
89   !                               !!!  simplified eos coefficients (default value: Vallis 2006)
90   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.
91   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.
92   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2       
93   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2       
94   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T 
95   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S 
96   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt 
97   
98   ! TEOS10/EOS80 parameters
99   REAL(wp) ::   r1_S0, r1_T0, r1_Z0, rdeltaS
100   
101   ! EOS parameters
102   REAL(wp) ::   EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
103   REAL(wp) ::   EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510
104   REAL(wp) ::   EOS020 , EOS120 , EOS220 , EOS320 , EOS420
105   REAL(wp) ::   EOS030 , EOS130 , EOS230 , EOS330
106   REAL(wp) ::   EOS040 , EOS140 , EOS240
107   REAL(wp) ::   EOS050 , EOS150
108   REAL(wp) ::   EOS060
109   REAL(wp) ::   EOS001 , EOS101 , EOS201 , EOS301 , EOS401
110   REAL(wp) ::   EOS011 , EOS111 , EOS211 , EOS311
111   REAL(wp) ::   EOS021 , EOS121 , EOS221
112   REAL(wp) ::   EOS031 , EOS131
113   REAL(wp) ::   EOS041
114   REAL(wp) ::   EOS002 , EOS102 , EOS202
115   REAL(wp) ::   EOS012 , EOS112
116   REAL(wp) ::   EOS022
117   REAL(wp) ::   EOS003 , EOS103
118   REAL(wp) ::   EOS013 
119   
120   ! ALPHA parameters
121   REAL(wp) ::   ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
122   REAL(wp) ::   ALP010 , ALP110 , ALP210 , ALP310 , ALP410
123   REAL(wp) ::   ALP020 , ALP120 , ALP220 , ALP320
124   REAL(wp) ::   ALP030 , ALP130 , ALP230
125   REAL(wp) ::   ALP040 , ALP140
126   REAL(wp) ::   ALP050
127   REAL(wp) ::   ALP001 , ALP101 , ALP201 , ALP301
128   REAL(wp) ::   ALP011 , ALP111 , ALP211
129   REAL(wp) ::   ALP021 , ALP121
130   REAL(wp) ::   ALP031
131   REAL(wp) ::   ALP002 , ALP102
132   REAL(wp) ::   ALP012
133   REAL(wp) ::   ALP003
134   
135   ! BETA parameters
136   REAL(wp) ::   BET000 , BET100 , BET200 , BET300 , BET400 , BET500
137   REAL(wp) ::   BET010 , BET110 , BET210 , BET310 , BET410
138   REAL(wp) ::   BET020 , BET120 , BET220 , BET320
139   REAL(wp) ::   BET030 , BET130 , BET230
140   REAL(wp) ::   BET040 , BET140
141   REAL(wp) ::   BET050
142   REAL(wp) ::   BET001 , BET101 , BET201 , BET301
143   REAL(wp) ::   BET011 , BET111 , BET211
144   REAL(wp) ::   BET021 , BET121
145   REAL(wp) ::   BET031
146   REAL(wp) ::   BET002 , BET102
147   REAL(wp) ::   BET012
148   REAL(wp) ::   BET003
149
150   ! PEN parameters
151   REAL(wp) ::   PEN000 , PEN100 , PEN200 , PEN300 , PEN400
152   REAL(wp) ::   PEN010 , PEN110 , PEN210 , PEN310
153   REAL(wp) ::   PEN020 , PEN120 , PEN220
154   REAL(wp) ::   PEN030 , PEN130
155   REAL(wp) ::   PEN040
156   REAL(wp) ::   PEN001 , PEN101 , PEN201
157   REAL(wp) ::   PEN011 , PEN111
158   REAL(wp) ::   PEN021
159   REAL(wp) ::   PEN002 , PEN102
160   REAL(wp) ::   PEN012
161   
162   ! ALPHA_PEN parameters
163   REAL(wp) ::   APE000 , APE100 , APE200 , APE300
164   REAL(wp) ::   APE010 , APE110 , APE210
165   REAL(wp) ::   APE020 , APE120
166   REAL(wp) ::   APE030
167   REAL(wp) ::   APE001 , APE101
168   REAL(wp) ::   APE011
169   REAL(wp) ::   APE002
170
171   ! BETA_PEN parameters
172   REAL(wp) ::   BPE000 , BPE100 , BPE200 , BPE300
173   REAL(wp) ::   BPE010 , BPE110 , BPE210
174   REAL(wp) ::   BPE020 , BPE120
175   REAL(wp) ::   BPE030
176   REAL(wp) ::   BPE001 , BPE101
177   REAL(wp) ::   BPE011
178   REAL(wp) ::   BPE002
179
180   !! * Substitutions
181#  include "do_loop_substitute.h90"
182   !!----------------------------------------------------------------------
183   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
184   !! $Id$
185   !! Software governed by the CeCILL license (see ./LICENSE)
186   !!----------------------------------------------------------------------
187CONTAINS
188
189   SUBROUTINE eos_insitu( pts, prd, pdep )
190      !!----------------------------------------------------------------------
191      !!                   ***  ROUTINE eos_insitu  ***
192      !!
193      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
194      !!       potential temperature and salinity using an equation of state
195      !!       selected in the nameos namelist
196      !!
197      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
198      !!         with   prd    in situ density anomaly      no units
199      !!                t      TEOS10: CT or EOS80: PT      Celsius
200      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu
201      !!                z      depth                        meters
202      !!                rho    in situ density              kg/m^3
203      !!                rau0   reference density            kg/m^3
204      !!
205      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
206      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
207      !!
208      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
209      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
210      !!
211      !!     ln_seos : simplified equation of state
212      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
213      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
214      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
215      !!              Vallis like equation: use default values of coefficients
216      !!
217      !! ** Action  :   compute prd , the in situ density (no units)
218      !!
219      !! References :   Roquet et al, Ocean Modelling, in preparation (2014)
220      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
221      !!                TEOS-10 Manual, 2010
222      !!----------------------------------------------------------------------
223      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
224      !                                                               ! 2 : salinity               [psu]
225      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd   ! in situ density            [-]
226      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep  ! depth                      [m]
227      !
228      INTEGER  ::   ji, jj, jk                ! dummy loop indices
229      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
230      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
231      !!----------------------------------------------------------------------
232      !
233      IF( ln_timing )   CALL timing_start('eos-insitu')
234      !
235      SELECT CASE( neos )
236      !
237      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
238         !
239         DO_3D_11_11( 1, jpkm1 )
240            !
241            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
242            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
243            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
244            ztm = tmask(ji,jj,jk)                                         ! tmask
245            !
246            zn3 = EOS013*zt   &
247               &   + EOS103*zs+EOS003
248               !
249            zn2 = (EOS022*zt   &
250               &   + EOS112*zs+EOS012)*zt   &
251               &   + (EOS202*zs+EOS102)*zs+EOS002
252               !
253            zn1 = (((EOS041*zt   &
254               &   + EOS131*zs+EOS031)*zt   &
255               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
256               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
257               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
258               !
259            zn0 = (((((EOS060*zt   &
260               &   + EOS150*zs+EOS050)*zt   &
261               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
262               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
263               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
264               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
265               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
266               !
267            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
268            !
269            prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked)
270            !
271         END_3D
272         !
273      CASE( np_seos )                !==  simplified EOS  ==!
274         !
275         DO_3D_11_11( 1, jpkm1 )
276            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
277            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
278            zh  = pdep (ji,jj,jk)
279            ztm = tmask(ji,jj,jk)
280            !
281            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
282               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
283               &  - rn_nu * zt * zs
284               !                                 
285            prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
286         END_3D
287         !
288      END SELECT
289      !
290      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
291      !
292      IF( ln_timing )   CALL timing_stop('eos-insitu')
293      !
294   END SUBROUTINE eos_insitu
295
296
297   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
298      !!----------------------------------------------------------------------
299      !!                  ***  ROUTINE eos_insitu_pot  ***
300      !!
301      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
302      !!      potential volumic mass (Kg/m3) from potential temperature and
303      !!      salinity fields using an equation of state selected in the
304      !!     namelist.
305      !!
306      !! ** Action  : - prd  , the in situ density (no units)
307      !!              - prhop, the potential volumic mass (Kg/m3)
308      !!
309      !!----------------------------------------------------------------------
310      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
311      !                                                                ! 2 : salinity               [psu]
312      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
313      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
314      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
315      !
316      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
317      INTEGER  ::   jdof
318      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
319      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
320      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
321      !!----------------------------------------------------------------------
322      !
323      IF( ln_timing )   CALL timing_start('eos-pot')
324      !
325      SELECT CASE ( neos )
326      !
327      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
328         !
329         ! Stochastic equation of state
330         IF ( ln_sto_eos ) THEN
331            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
332            ALLOCATE(zn_sto(1:2*nn_sto_eos))
333            ALLOCATE(zsign(1:2*nn_sto_eos))
334            DO jsmp = 1, 2*nn_sto_eos, 2
335              zsign(jsmp)   = 1._wp
336              zsign(jsmp+1) = -1._wp
337            END DO
338            !
339            DO_3D_11_11( 1, jpkm1 )
340               !
341               ! compute density (2*nn_sto_eos) times:
342               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
343               ! (2) for t-dt, s-ds (with the opposite fluctuation)
344               DO jsmp = 1, nn_sto_eos*2
345                  jdof   = (jsmp + 1) / 2
346                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
347                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
348                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
349                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
350                  ztm    = tmask(ji,jj,jk)                                         ! tmask
351                  !
352                  zn3 = EOS013*zt   &
353                     &   + EOS103*zs+EOS003
354                     !
355                  zn2 = (EOS022*zt   &
356                     &   + EOS112*zs+EOS012)*zt   &
357                     &   + (EOS202*zs+EOS102)*zs+EOS002
358                     !
359                  zn1 = (((EOS041*zt   &
360                     &   + EOS131*zs+EOS031)*zt   &
361                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
362                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
363                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
364                     !
365                  zn0_sto(jsmp) = (((((EOS060*zt   &
366                     &   + EOS150*zs+EOS050)*zt   &
367                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
368                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
369                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
370                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
371                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
372                     !
373                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
374               END DO
375               !
376               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
377               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
378               DO jsmp = 1, nn_sto_eos*2
379                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
380                  !
381                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked)
382               END DO
383               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
384               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
385            END_3D
386            DEALLOCATE(zn0_sto,zn_sto,zsign)
387         ! Non-stochastic equation of state
388         ELSE
389            DO_3D_11_11( 1, jpkm1 )
390               !
391               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
392               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
393               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
394               ztm = tmask(ji,jj,jk)                                         ! tmask
395               !
396               zn3 = EOS013*zt   &
397                  &   + EOS103*zs+EOS003
398                  !
399               zn2 = (EOS022*zt   &
400                  &   + EOS112*zs+EOS012)*zt   &
401                  &   + (EOS202*zs+EOS102)*zs+EOS002
402                  !
403               zn1 = (((EOS041*zt   &
404                  &   + EOS131*zs+EOS031)*zt   &
405                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
406                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
407                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
408                  !
409               zn0 = (((((EOS060*zt   &
410                  &   + EOS150*zs+EOS050)*zt   &
411                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
412                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
413                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
414                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
415                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
416                  !
417               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
418               !
419               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
420               !
421               prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
422            END_3D
423         ENDIF
424         
425      CASE( np_seos )                !==  simplified EOS  ==!
426         !
427         DO_3D_11_11( 1, jpkm1 )
428            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
429            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
430            zh  = pdep (ji,jj,jk)
431            ztm = tmask(ji,jj,jk)
432            !                                                     ! potential density referenced at the surface
433            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
434               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
435               &  - rn_nu * zt * zs
436            prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
437            !                                                     ! density anomaly (masked)
438            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
439            prd(ji,jj,jk) = zn * r1_rau0 * ztm
440            !
441         END_3D
442         !
443      END SELECT
444      !
445      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
446      !
447      IF( ln_timing )   CALL timing_stop('eos-pot')
448      !
449   END SUBROUTINE eos_insitu_pot
450
451
452   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
453      !!----------------------------------------------------------------------
454      !!                  ***  ROUTINE eos_insitu_2d  ***
455      !!
456      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
457      !!      potential temperature and salinity using an equation of state
458      !!      selected in the nameos namelist. * 2D field case
459      !!
460      !! ** Action  : - prd , the in situ density (no units) (unmasked)
461      !!
462      !!----------------------------------------------------------------------
463      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
464      !                                                           ! 2 : salinity               [psu]
465      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
466      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
467      !
468      INTEGER  ::   ji, jj, jk                ! dummy loop indices
469      REAL(wp) ::   zt , zh , zs              ! local scalars
470      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
471      !!----------------------------------------------------------------------
472      !
473      IF( ln_timing )   CALL timing_start('eos2d')
474      !
475      prd(:,:) = 0._wp
476      !
477      SELECT CASE( neos )
478      !
479      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
480         !
481         DO_2D_11_11
482            !
483            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
484            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
485            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
486            !
487            zn3 = EOS013*zt   &
488               &   + EOS103*zs+EOS003
489               !
490            zn2 = (EOS022*zt   &
491               &   + EOS112*zs+EOS012)*zt   &
492               &   + (EOS202*zs+EOS102)*zs+EOS002
493               !
494            zn1 = (((EOS041*zt   &
495               &   + EOS131*zs+EOS031)*zt   &
496               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
497               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
498               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
499               !
500            zn0 = (((((EOS060*zt   &
501               &   + EOS150*zs+EOS050)*zt   &
502               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
503               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
504               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
505               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
506               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
507               !
508            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
509            !
510            prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
511            !
512         END_2D
513         !
514      CASE( np_seos )                !==  simplified EOS  ==!
515         !
516         DO_2D_11_11
517            !
518            zt    = pts  (ji,jj,jp_tem)  - 10._wp
519            zs    = pts  (ji,jj,jp_sal)  - 35._wp
520            zh    = pdep (ji,jj)                         ! depth at the partial step level
521            !
522            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
523               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
524               &  - rn_nu * zt * zs
525               !
526            prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
527            !
528         END_2D
529         !
530      END SELECT
531      !
532      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
533      !
534      IF( ln_timing )   CALL timing_stop('eos2d')
535      !
536   END SUBROUTINE eos_insitu_2d
537
538
539   SUBROUTINE rab_3d( pts, pab, Kmm )
540      !!----------------------------------------------------------------------
541      !!                 ***  ROUTINE rab_3d  ***
542      !!
543      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points
544      !!
545      !! ** Method  :   calculates alpha / beta at T-points
546      !!
547      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
548      !!----------------------------------------------------------------------
549      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
550      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
551      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
552      !
553      INTEGER  ::   ji, jj, jk                ! dummy loop indices
554      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
555      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
556      !!----------------------------------------------------------------------
557      !
558      IF( ln_timing )   CALL timing_start('rab_3d')
559      !
560      SELECT CASE ( neos )
561      !
562      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
563         !
564         DO_3D_11_11( 1, jpkm1 )
565            !
566            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
567            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
568            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
569            ztm = tmask(ji,jj,jk)                                         ! tmask
570            !
571            ! alpha
572            zn3 = ALP003
573            !
574            zn2 = ALP012*zt + ALP102*zs+ALP002
575            !
576            zn1 = ((ALP031*zt   &
577               &   + ALP121*zs+ALP021)*zt   &
578               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
579               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
580               !
581            zn0 = ((((ALP050*zt   &
582               &   + ALP140*zs+ALP040)*zt   &
583               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
584               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
585               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
586               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
587               !
588            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
589            !
590            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
591            !
592            ! beta
593            zn3 = BET003
594            !
595            zn2 = BET012*zt + BET102*zs+BET002
596            !
597            zn1 = ((BET031*zt   &
598               &   + BET121*zs+BET021)*zt   &
599               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
600               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
601               !
602            zn0 = ((((BET050*zt   &
603               &   + BET140*zs+BET040)*zt   &
604               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
605               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
606               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
607               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
608               !
609            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
610            !
611            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
612            !
613         END_3D
614         !
615      CASE( np_seos )                  !==  simplified EOS  ==!
616         !
617         DO_3D_11_11( 1, jpkm1 )
618            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
619            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
620            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
621            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
622            !
623            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
624            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
625            !
626            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
627            pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
628            !
629         END_3D
630         !
631      CASE DEFAULT
632         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
633         CALL ctl_stop( 'rab_3d:', ctmp1 )
634         !
635      END SELECT
636      !
637      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
638         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
639      !
640      IF( ln_timing )   CALL timing_stop('rab_3d')
641      !
642   END SUBROUTINE rab_3d
643
644
645   SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
646      !!----------------------------------------------------------------------
647      !!                 ***  ROUTINE rab_2d  ***
648      !!
649      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
650      !!
651      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
652      !!----------------------------------------------------------------------
653      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
654      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
655      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m]
656      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
657      !
658      INTEGER  ::   ji, jj, jk                ! dummy loop indices
659      REAL(wp) ::   zt , zh , zs              ! local scalars
660      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
661      !!----------------------------------------------------------------------
662      !
663      IF( ln_timing )   CALL timing_start('rab_2d')
664      !
665      pab(:,:,:) = 0._wp
666      !
667      SELECT CASE ( neos )
668      !
669      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
670         !
671         DO_2D_11_11
672            !
673            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
674            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
675            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
676            !
677            ! alpha
678            zn3 = ALP003
679            !
680            zn2 = ALP012*zt + ALP102*zs+ALP002
681            !
682            zn1 = ((ALP031*zt   &
683               &   + ALP121*zs+ALP021)*zt   &
684               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
685               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
686               !
687            zn0 = ((((ALP050*zt   &
688               &   + ALP140*zs+ALP040)*zt   &
689               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
690               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
691               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
692               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
693               !
694            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
695            !
696            pab(ji,jj,jp_tem) = zn * r1_rau0
697            !
698            ! beta
699            zn3 = BET003
700            !
701            zn2 = BET012*zt + BET102*zs+BET002
702            !
703            zn1 = ((BET031*zt   &
704               &   + BET121*zs+BET021)*zt   &
705               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
706               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
707               !
708            zn0 = ((((BET050*zt   &
709               &   + BET140*zs+BET040)*zt   &
710               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
711               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
712               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
713               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
714               !
715            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
716            !
717            pab(ji,jj,jp_sal) = zn / zs * r1_rau0
718            !
719            !
720         END_2D
721         !
722      CASE( np_seos )                  !==  simplified EOS  ==!
723         !
724         DO_2D_11_11
725            !
726            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
727            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
728            zh    = pdep (ji,jj)                   ! depth at the partial step level
729            !
730            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
731            pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
732            !
733            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
734            pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
735            !
736         END_2D
737         !
738      CASE DEFAULT
739         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
740         CALL ctl_stop( 'rab_2d:', ctmp1 )
741         !
742      END SELECT
743      !
744      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
745         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
746      !
747      IF( ln_timing )   CALL timing_stop('rab_2d')
748      !
749   END SUBROUTINE rab_2d
750
751
752   SUBROUTINE rab_0d( pts, pdep, pab, Kmm )
753      !!----------------------------------------------------------------------
754      !!                 ***  ROUTINE rab_0d  ***
755      !!
756      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
757      !!
758      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
759      !!----------------------------------------------------------------------
760      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
761      REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
762      REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m]
763      REAL(wp), DIMENSION(jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
764      !
765      REAL(wp) ::   zt , zh , zs              ! local scalars
766      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
767      !!----------------------------------------------------------------------
768      !
769      IF( ln_timing )   CALL timing_start('rab_0d')
770      !
771      pab(:) = 0._wp
772      !
773      SELECT CASE ( neos )
774      !
775      CASE( np_teos10, np_eos80 )      !==  polynomial TEOS-10 / EOS-80 ==!
776         !
777         !
778         zh  = pdep * r1_Z0                                  ! depth
779         zt  = pts (jp_tem) * r1_T0                           ! temperature
780         zs  = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
781         !
782         ! alpha
783         zn3 = ALP003
784         !
785         zn2 = ALP012*zt + ALP102*zs+ALP002
786         !
787         zn1 = ((ALP031*zt   &
788            &   + ALP121*zs+ALP021)*zt   &
789            &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
790            &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
791            !
792         zn0 = ((((ALP050*zt   &
793            &   + ALP140*zs+ALP040)*zt   &
794            &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
795            &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
796            &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
797            &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
798            !
799         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
800         !
801         pab(jp_tem) = zn * r1_rau0
802         !
803         ! beta
804         zn3 = BET003
805         !
806         zn2 = BET012*zt + BET102*zs+BET002
807         !
808         zn1 = ((BET031*zt   &
809            &   + BET121*zs+BET021)*zt   &
810            &   + (BET211*zs+BET111)*zs+BET011)*zt   &
811            &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
812            !
813         zn0 = ((((BET050*zt   &
814            &   + BET140*zs+BET040)*zt   &
815            &   + (BET230*zs+BET130)*zs+BET030)*zt   &
816            &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
817            &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
818            &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
819            !
820         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
821         !
822         pab(jp_sal) = zn / zs * r1_rau0
823         !
824         !
825         !
826      CASE( np_seos )                  !==  simplified EOS  ==!
827         !
828         zt    = pts(jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
829         zs    = pts(jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
830         zh    = pdep                   ! depth at the partial step level
831         !
832         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
833         pab(jp_tem) = zn * r1_rau0   ! alpha
834         !
835         zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
836         pab(jp_sal) = zn * r1_rau0   ! beta
837         !
838      CASE DEFAULT
839         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
840         CALL ctl_stop( 'rab_0d:', ctmp1 )
841         !
842      END SELECT
843      !
844      IF( ln_timing )   CALL timing_stop('rab_0d')
845      !
846   END SUBROUTINE rab_0d
847
848
849   SUBROUTINE bn2( pts, pab, pn2, Kmm )
850      !!----------------------------------------------------------------------
851      !!                  ***  ROUTINE bn2  ***
852      !!
853      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the
854      !!                time-step of the input arguments
855      !!
856      !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
857      !!      where alpha and beta are given in pab, and computed on T-points.
858      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module.
859      !!
860      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point
861      !!
862      !!----------------------------------------------------------------------
863      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
864      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
865      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
866      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
867      !
868      INTEGER  ::   ji, jj, jk      ! dummy loop indices
869      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
870      !!----------------------------------------------------------------------
871      !
872      IF( ln_timing )   CALL timing_start('bn2')
873      !
874      DO_3D_11_11( 2, jpkm1 )
875         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
876            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
877            !
878         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
879         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
880         !
881         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
882            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
883            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
884      END_3D
885      !
886      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk )
887      !
888      IF( ln_timing )   CALL timing_stop('bn2')
889      !
890   END SUBROUTINE bn2
891
892
893   FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
894      !!----------------------------------------------------------------------
895      !!                 ***  ROUTINE eos_pt_from_ct  ***
896      !!
897      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celsius]
898      !!
899      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
900      !!       checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
901      !!
902      !! Reference  :   TEOS-10, UNESCO
903      !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
904      !!----------------------------------------------------------------------
905      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp   [Celsius]
906      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity     [psu]
907      ! Leave result array automatic rather than making explicitly allocated
908      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celsius]
909      !
910      INTEGER  ::   ji, jj               ! dummy loop indices
911      REAL(wp) ::   zt , zs , ztm        ! local scalars
912      REAL(wp) ::   zn , zd              ! local scalars
913      REAL(wp) ::   zdeltaS , z1_S0 , z1_T0
914      !!----------------------------------------------------------------------
915      !
916      IF( ln_timing )   CALL timing_start('eos_pt_from_ct')
917      !
918      zdeltaS = 5._wp
919      z1_S0   = 0.875_wp/35.16504_wp
920      z1_T0   = 1._wp/40._wp
921      !
922      DO_2D_11_11
923         !
924         zt  = ctmp   (ji,jj) * z1_T0
925         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
926         ztm = tmask(ji,jj,1)
927         !
928         zn = ((((-2.1385727895e-01_wp*zt   &
929            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
930            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
931            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
932            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
933            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
934            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
935            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
936            !
937         zd = (2.0035003456_wp*zt   &
938            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
939            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
940            !
941         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
942            !
943      END_2D
944      !
945      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct')
946      !
947   END FUNCTION eos_pt_from_ct
948
949
950   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep )
951      !!----------------------------------------------------------------------
952      !!                 ***  ROUTINE eos_fzp  ***
953      !!
954      !! ** Purpose :   Compute the freezing point temperature [Celsius]
955      !!
956      !! ** Method  :   UNESCO freezing point (ptf) in Celsius is given by
957      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
958      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
959      !!
960      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
961      !!----------------------------------------------------------------------
962      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
963      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
964      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celsius]
965      !
966      INTEGER  ::   ji, jj          ! dummy loop indices
967      REAL(wp) ::   zt, zs, z1_S0   ! local scalars
968      !!----------------------------------------------------------------------
969      !
970      SELECT CASE ( neos )
971      !
972      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==!
973         !
974         z1_S0 = 1._wp / 35.16504_wp
975         DO_2D_11_11
976            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
977            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
978               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
979         END_2D
980         ptf(:,:) = ptf(:,:) * psal(:,:)
981         !
982         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
983         !
984      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
985         !
986         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
987            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
988            !
989         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
990         !
991      CASE DEFAULT
992         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
993         CALL ctl_stop( 'eos_fzp_2d:', ctmp1 )
994         !
995      END SELECT     
996      !
997  END SUBROUTINE eos_fzp_2d
998
999
1000  SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
1001      !!----------------------------------------------------------------------
1002      !!                 ***  ROUTINE eos_fzp  ***
1003      !!
1004      !! ** Purpose :   Compute the freezing point temperature [Celsius]
1005      !!
1006      !! ** Method  :   UNESCO freezing point (ptf) in Celsius is given by
1007      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
1008      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
1009      !!
1010      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
1011      !!----------------------------------------------------------------------
1012      REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu]
1013      REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m]
1014      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celsius]
1015      !
1016      REAL(wp) :: zs   ! local scalars
1017      !!----------------------------------------------------------------------
1018      !
1019      SELECT CASE ( neos )
1020      !
1021      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==!
1022         !
1023         zs  = SQRT( ABS( psal ) / 35.16504_wp )           ! square root salinity
1024         ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
1025                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
1026         ptf = ptf * psal
1027         !
1028         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
1029         !
1030      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
1031         !
1032         ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   &
1033            &                - 2.154996e-4_wp *       psal   ) * psal
1034            !
1035         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
1036         !
1037      CASE DEFAULT
1038         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
1039         CALL ctl_stop( 'eos_fzp_0d:', ctmp1 )
1040         !
1041      END SELECT
1042      !
1043   END SUBROUTINE eos_fzp_0d
1044
1045
1046   SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm )
1047      !!----------------------------------------------------------------------
1048      !!                 ***  ROUTINE eos_pen  ***
1049      !!
1050      !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
1051      !!
1052      !! ** Method  :   PE is defined analytically as the vertical
1053      !!                   primitive of EOS times -g integrated between 0 and z>0.
1054      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
1055      !!                                                      = 1/z * /int_0^z rd dz - rd
1056      !!                                where rd is the density anomaly (see eos_rhd function)
1057      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
1058      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
1059      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS
1060      !!
1061      !! ** Action  : - pen         : PE anomaly given at T-points
1062      !!            : - pab_pe  : given at T-points
1063      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe
1064      !!                    pab_pe(:,:,:,jp_sal) is beta_pe
1065      !!----------------------------------------------------------------------
1066      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
1067      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
1068      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
1069      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   ppen     ! potential energy anomaly
1070      !
1071      INTEGER  ::   ji, jj, jk                ! dummy loop indices
1072      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
1073      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
1074      !!----------------------------------------------------------------------
1075      !
1076      IF( ln_timing )   CALL timing_start('eos_pen')
1077      !
1078      SELECT CASE ( neos )
1079      !
1080      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
1081         !
1082         DO_3D_11_11( 1, jpkm1 )
1083            !
1084            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
1085            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
1086            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
1087            ztm = tmask(ji,jj,jk)                                         ! tmask
1088            !
1089            ! potential energy non-linear anomaly
1090            zn2 = (PEN012)*zt   &
1091               &   + PEN102*zs+PEN002
1092               !
1093            zn1 = ((PEN021)*zt   &
1094               &   + PEN111*zs+PEN011)*zt   &
1095               &   + (PEN201*zs+PEN101)*zs+PEN001
1096               !
1097            zn0 = ((((PEN040)*zt   &
1098               &   + PEN130*zs+PEN030)*zt   &
1099               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
1100               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
1101               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
1102               !
1103            zn  = ( zn2 * zh + zn1 ) * zh + zn0
1104            !
1105            ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm
1106            !
1107            ! alphaPE non-linear anomaly
1108            zn2 = APE002
1109            !
1110            zn1 = (APE011)*zt   &
1111               &   + APE101*zs+APE001
1112               !
1113            zn0 = (((APE030)*zt   &
1114               &   + APE120*zs+APE020)*zt   &
1115               &   + (APE210*zs+APE110)*zs+APE010)*zt   &
1116               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
1117               !
1118            zn  = ( zn2 * zh + zn1 ) * zh + zn0
1119            !                             
1120            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
1121            !
1122            ! betaPE non-linear anomaly
1123            zn2 = BPE002
1124            !
1125            zn1 = (BPE011)*zt   &
1126               &   + BPE101*zs+BPE001
1127               !
1128            zn0 = (((BPE030)*zt   &
1129               &   + BPE120*zs+BPE020)*zt   &
1130               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
1131               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
1132               !
1133            zn  = ( zn2 * zh + zn1 ) * zh + zn0
1134            !                             
1135            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
1136            !
1137         END_3D
1138         !
1139      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
1140         !
1141         DO_3D_11_11( 1, jpkm1 )
1142            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
1143            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
1144            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
1145            ztm = tmask(ji,jj,jk)                ! tmask
1146            zn  = 0.5_wp * zh * r1_rau0 * ztm
1147            !                                    ! Potential Energy
1148            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
1149            !                                    ! alphaPE
1150            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
1151            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
1152            !
1153         END_3D
1154         !
1155      CASE DEFAULT
1156         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
1157         CALL ctl_stop( 'eos_pen:', ctmp1 )
1158         !
1159      END SELECT
1160      !
1161      IF( ln_timing )   CALL timing_stop('eos_pen')
1162      !
1163   END SUBROUTINE eos_pen
1164
1165
1166   SUBROUTINE eos_init
1167      !!----------------------------------------------------------------------
1168      !!                 ***  ROUTINE eos_init  ***
1169      !!
1170      !! ** Purpose :   initializations for the equation of state
1171      !!
1172      !! ** Method  :   Read the namelist nameos and control the parameters
1173      !!----------------------------------------------------------------------
1174      INTEGER  ::   ios   ! local integer
1175      INTEGER  ::   ioptio   ! local integer
1176      !!
1177      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
1178         &                                             rn_lambda2, rn_mu2, rn_nu
1179      !!----------------------------------------------------------------------
1180      !
1181      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
1182901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' )
1183      !
1184      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
1185902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' )
1186      IF(lwm) WRITE( numond, nameos )
1187      !
1188      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
1189      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
1190      !
1191      IF(lwp) THEN                ! Control print
1192         WRITE(numout,*)
1193         WRITE(numout,*) 'eos_init : equation of state'
1194         WRITE(numout,*) '~~~~~~~~'
1195         WRITE(numout,*) '   Namelist nameos : Chosen the Equation Of Seawater (EOS)'
1196         WRITE(numout,*) '      TEOS-10 : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_TEOS10 = ', ln_TEOS10
1197         WRITE(numout,*) '      EOS-80  : rho=F(Potential    Temperature, Practical Salinity, depth)   ln_EOS80  = ', ln_EOS80
1198         WRITE(numout,*) '      S-EOS   : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_SEOS   = ', ln_SEOS
1199      ENDIF
1200
1201      ! Check options for equation of state & set neos based on logical flags
1202      ioptio = 0
1203      IF( ln_TEOS10 ) THEN   ;   ioptio = ioptio+1   ;   neos = np_teos10   ;   ENDIF
1204      IF( ln_EOS80  ) THEN   ;   ioptio = ioptio+1   ;   neos = np_eos80    ;   ENDIF
1205      IF( ln_SEOS   ) THEN   ;   ioptio = ioptio+1   ;   neos = np_seos     ;   ENDIF
1206      IF( ioptio /= 1 )   CALL ctl_stop("Exactly one equation of state option must be selected")
1207      !
1208      SELECT CASE( neos )         ! check option
1209      !
1210      CASE( np_teos10 )                       !==  polynomial TEOS-10  ==!
1211         IF(lwp) WRITE(numout,*)
1212         IF(lwp) WRITE(numout,*) '   ==>>>   use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
1213         !
1214         l_useCT = .TRUE.                          ! model temperature is Conservative temperature
1215         !
1216         rdeltaS = 32._wp
1217         r1_S0  = 0.875_wp/35.16504_wp
1218         r1_T0  = 1._wp/40._wp
1219         r1_Z0  = 1.e-4_wp
1220         !
1221         EOS000 = 8.0189615746e+02_wp
1222         EOS100 = 8.6672408165e+02_wp
1223         EOS200 = -1.7864682637e+03_wp
1224         EOS300 = 2.0375295546e+03_wp
1225         EOS400 = -1.2849161071e+03_wp
1226         EOS500 = 4.3227585684e+02_wp
1227         EOS600 = -6.0579916612e+01_wp
1228         EOS010 = 2.6010145068e+01_wp
1229         EOS110 = -6.5281885265e+01_wp
1230         EOS210 = 8.1770425108e+01_wp
1231         EOS310 = -5.6888046321e+01_wp
1232         EOS410 = 1.7681814114e+01_wp
1233         EOS510 = -1.9193502195_wp
1234         EOS020 = -3.7074170417e+01_wp
1235         EOS120 = 6.1548258127e+01_wp
1236         EOS220 = -6.0362551501e+01_wp
1237         EOS320 = 2.9130021253e+01_wp
1238         EOS420 = -5.4723692739_wp
1239         EOS030 = 2.1661789529e+01_wp
1240         EOS130 = -3.3449108469e+01_wp
1241         EOS230 = 1.9717078466e+01_wp
1242         EOS330 = -3.1742946532_wp
1243         EOS040 = -8.3627885467_wp
1244         EOS140 = 1.1311538584e+01_wp
1245         EOS240 = -5.3563304045_wp
1246         EOS050 = 5.4048723791e-01_wp
1247         EOS150 = 4.8169980163e-01_wp
1248         EOS060 = -1.9083568888e-01_wp
1249         EOS001 = 1.9681925209e+01_wp
1250         EOS101 = -4.2549998214e+01_wp
1251         EOS201 = 5.0774768218e+01_wp
1252         EOS301 = -3.0938076334e+01_wp
1253         EOS401 = 6.6051753097_wp
1254         EOS011 = -1.3336301113e+01_wp
1255         EOS111 = -4.4870114575_wp
1256         EOS211 = 5.0042598061_wp
1257         EOS311 = -6.5399043664e-01_wp
1258         EOS021 = 6.7080479603_wp
1259         EOS121 = 3.5063081279_wp
1260         EOS221 = -1.8795372996_wp
1261         EOS031 = -2.4649669534_wp
1262         EOS131 = -5.5077101279e-01_wp
1263         EOS041 = 5.5927935970e-01_wp
1264         EOS002 = 2.0660924175_wp
1265         EOS102 = -4.9527603989_wp
1266         EOS202 = 2.5019633244_wp
1267         EOS012 = 2.0564311499_wp
1268         EOS112 = -2.1311365518e-01_wp
1269         EOS022 = -1.2419983026_wp
1270         EOS003 = -2.3342758797e-02_wp
1271         EOS103 = -1.8507636718e-02_wp
1272         EOS013 = 3.7969820455e-01_wp
1273         !
1274         ALP000 = -6.5025362670e-01_wp
1275         ALP100 = 1.6320471316_wp
1276         ALP200 = -2.0442606277_wp
1277         ALP300 = 1.4222011580_wp
1278         ALP400 = -4.4204535284e-01_wp
1279         ALP500 = 4.7983755487e-02_wp
1280         ALP010 = 1.8537085209_wp
1281         ALP110 = -3.0774129064_wp
1282         ALP210 = 3.0181275751_wp
1283         ALP310 = -1.4565010626_wp
1284         ALP410 = 2.7361846370e-01_wp
1285         ALP020 = -1.6246342147_wp
1286         ALP120 = 2.5086831352_wp
1287         ALP220 = -1.4787808849_wp
1288         ALP320 = 2.3807209899e-01_wp
1289         ALP030 = 8.3627885467e-01_wp
1290         ALP130 = -1.1311538584_wp
1291         ALP230 = 5.3563304045e-01_wp
1292         ALP040 = -6.7560904739e-02_wp
1293         ALP140 = -6.0212475204e-02_wp
1294         ALP050 = 2.8625353333e-02_wp
1295         ALP001 = 3.3340752782e-01_wp
1296         ALP101 = 1.1217528644e-01_wp
1297         ALP201 = -1.2510649515e-01_wp
1298         ALP301 = 1.6349760916e-02_wp
1299         ALP011 = -3.3540239802e-01_wp
1300         ALP111 = -1.7531540640e-01_wp
1301         ALP211 = 9.3976864981e-02_wp
1302         ALP021 = 1.8487252150e-01_wp
1303         ALP121 = 4.1307825959e-02_wp
1304         ALP031 = -5.5927935970e-02_wp
1305         ALP002 = -5.1410778748e-02_wp
1306         ALP102 = 5.3278413794e-03_wp
1307         ALP012 = 6.2099915132e-02_wp
1308         ALP003 = -9.4924551138e-03_wp
1309         !
1310         BET000 = 1.0783203594e+01_wp
1311         BET100 = -4.4452095908e+01_wp
1312         BET200 = 7.6048755820e+01_wp
1313         BET300 = -6.3944280668e+01_wp
1314         BET400 = 2.6890441098e+01_wp
1315         BET500 = -4.5221697773_wp
1316         BET010 = -8.1219372432e-01_wp
1317         BET110 = 2.0346663041_wp
1318         BET210 = -2.1232895170_wp
1319         BET310 = 8.7994140485e-01_wp
1320         BET410 = -1.1939638360e-01_wp
1321         BET020 = 7.6574242289e-01_wp
1322         BET120 = -1.5019813020_wp
1323         BET220 = 1.0872489522_wp
1324         BET320 = -2.7233429080e-01_wp
1325         BET030 = -4.1615152308e-01_wp
1326         BET130 = 4.9061350869e-01_wp
1327         BET230 = -1.1847737788e-01_wp
1328         BET040 = 1.4073062708e-01_wp
1329         BET140 = -1.3327978879e-01_wp
1330         BET050 = 5.9929880134e-03_wp
1331         BET001 = -5.2937873009e-01_wp
1332         BET101 = 1.2634116779_wp
1333         BET201 = -1.1547328025_wp
1334         BET301 = 3.2870876279e-01_wp
1335         BET011 = -5.5824407214e-02_wp
1336         BET111 = 1.2451933313e-01_wp
1337         BET211 = -2.4409539932e-02_wp
1338         BET021 = 4.3623149752e-02_wp
1339         BET121 = -4.6767901790e-02_wp
1340         BET031 = -6.8523260060e-03_wp
1341         BET002 = -6.1618945251e-02_wp
1342         BET102 = 6.2255521644e-02_wp
1343         BET012 = -2.6514181169e-03_wp
1344         BET003 = -2.3025968587e-04_wp
1345         !
1346         PEN000 = -9.8409626043_wp
1347         PEN100 = 2.1274999107e+01_wp
1348         PEN200 = -2.5387384109e+01_wp
1349         PEN300 = 1.5469038167e+01_wp
1350         PEN400 = -3.3025876549_wp
1351         PEN010 = 6.6681505563_wp
1352         PEN110 = 2.2435057288_wp
1353         PEN210 = -2.5021299030_wp
1354         PEN310 = 3.2699521832e-01_wp
1355         PEN020 = -3.3540239802_wp
1356         PEN120 = -1.7531540640_wp
1357         PEN220 = 9.3976864981e-01_wp
1358         PEN030 = 1.2324834767_wp
1359         PEN130 = 2.7538550639e-01_wp
1360         PEN040 = -2.7963967985e-01_wp
1361         PEN001 = -1.3773949450_wp
1362         PEN101 = 3.3018402659_wp
1363         PEN201 = -1.6679755496_wp
1364         PEN011 = -1.3709540999_wp
1365         PEN111 = 1.4207577012e-01_wp
1366         PEN021 = 8.2799886843e-01_wp
1367         PEN002 = 1.7507069098e-02_wp
1368         PEN102 = 1.3880727538e-02_wp
1369         PEN012 = -2.8477365341e-01_wp
1370         !
1371         APE000 = -1.6670376391e-01_wp
1372         APE100 = -5.6087643219e-02_wp
1373         APE200 = 6.2553247576e-02_wp
1374         APE300 = -8.1748804580e-03_wp
1375         APE010 = 1.6770119901e-01_wp
1376         APE110 = 8.7657703198e-02_wp
1377         APE210 = -4.6988432490e-02_wp
1378         APE020 = -9.2436260751e-02_wp
1379         APE120 = -2.0653912979e-02_wp
1380         APE030 = 2.7963967985e-02_wp
1381         APE001 = 3.4273852498e-02_wp
1382         APE101 = -3.5518942529e-03_wp
1383         APE011 = -4.1399943421e-02_wp
1384         APE002 = 7.1193413354e-03_wp
1385         !
1386         BPE000 = 2.6468936504e-01_wp
1387         BPE100 = -6.3170583896e-01_wp
1388         BPE200 = 5.7736640125e-01_wp
1389         BPE300 = -1.6435438140e-01_wp
1390         BPE010 = 2.7912203607e-02_wp
1391         BPE110 = -6.2259666565e-02_wp
1392         BPE210 = 1.2204769966e-02_wp
1393         BPE020 = -2.1811574876e-02_wp
1394         BPE120 = 2.3383950895e-02_wp
1395         BPE030 = 3.4261630030e-03_wp
1396         BPE001 = 4.1079296834e-02_wp
1397         BPE101 = -4.1503681096e-02_wp
1398         BPE011 = 1.7676120780e-03_wp
1399         BPE002 = 1.7269476440e-04_wp
1400         !
1401      CASE( np_eos80 )                        !==  polynomial EOS-80 formulation  ==!
1402         !
1403         IF(lwp) WRITE(numout,*)
1404         IF(lwp) WRITE(numout,*) '   ==>>>   use of EOS-80 equation of state (pot. temp. and pract. salinity)'
1405         !
1406         l_useCT = .FALSE.                         ! model temperature is Potential temperature
1407         rdeltaS = 20._wp
1408         r1_S0  = 1._wp/40._wp
1409         r1_T0  = 1._wp/40._wp
1410         r1_Z0  = 1.e-4_wp
1411         !
1412         EOS000 = 9.5356891948e+02_wp
1413         EOS100 = 1.7136499189e+02_wp
1414         EOS200 = -3.7501039454e+02_wp
1415         EOS300 = 5.1856810420e+02_wp
1416         EOS400 = -3.7264470465e+02_wp
1417         EOS500 = 1.4302533998e+02_wp
1418         EOS600 = -2.2856621162e+01_wp
1419         EOS010 = 1.0087518651e+01_wp
1420         EOS110 = -1.3647741861e+01_wp
1421         EOS210 = 8.8478359933_wp
1422         EOS310 = -7.2329388377_wp
1423         EOS410 = 1.4774410611_wp
1424         EOS510 = 2.0036720553e-01_wp
1425         EOS020 = -2.5579830599e+01_wp
1426         EOS120 = 2.4043512327e+01_wp
1427         EOS220 = -1.6807503990e+01_wp
1428         EOS320 = 8.3811577084_wp
1429         EOS420 = -1.9771060192_wp
1430         EOS030 = 1.6846451198e+01_wp
1431         EOS130 = -2.1482926901e+01_wp
1432         EOS230 = 1.0108954054e+01_wp
1433         EOS330 = -6.2675951440e-01_wp
1434         EOS040 = -8.0812310102_wp
1435         EOS140 = 1.0102374985e+01_wp
1436         EOS240 = -4.8340368631_wp
1437         EOS050 = 1.2079167803_wp
1438         EOS150 = 1.1515380987e-01_wp
1439         EOS060 = -2.4520288837e-01_wp
1440         EOS001 = 1.0748601068e+01_wp
1441         EOS101 = -1.7817043500e+01_wp
1442         EOS201 = 2.2181366768e+01_wp
1443         EOS301 = -1.6750916338e+01_wp
1444         EOS401 = 4.1202230403_wp
1445         EOS011 = -1.5852644587e+01_wp
1446         EOS111 = -7.6639383522e-01_wp
1447         EOS211 = 4.1144627302_wp
1448         EOS311 = -6.6955877448e-01_wp
1449         EOS021 = 9.9994861860_wp
1450         EOS121 = -1.9467067787e-01_wp
1451         EOS221 = -1.2177554330_wp
1452         EOS031 = -3.4866102017_wp
1453         EOS131 = 2.2229155620e-01_wp
1454         EOS041 = 5.9503008642e-01_wp
1455         EOS002 = 1.0375676547_wp
1456         EOS102 = -3.4249470629_wp
1457         EOS202 = 2.0542026429_wp
1458         EOS012 = 2.1836324814_wp
1459         EOS112 = -3.4453674320e-01_wp
1460         EOS022 = -1.2548163097_wp
1461         EOS003 = 1.8729078427e-02_wp
1462         EOS103 = -5.7238495240e-02_wp
1463         EOS013 = 3.8306136687e-01_wp
1464         !
1465         ALP000 = -2.5218796628e-01_wp
1466         ALP100 = 3.4119354654e-01_wp
1467         ALP200 = -2.2119589983e-01_wp
1468         ALP300 = 1.8082347094e-01_wp
1469         ALP400 = -3.6936026529e-02_wp
1470         ALP500 = -5.0091801383e-03_wp
1471         ALP010 = 1.2789915300_wp
1472         ALP110 = -1.2021756164_wp
1473         ALP210 = 8.4037519952e-01_wp
1474         ALP310 = -4.1905788542e-01_wp
1475         ALP410 = 9.8855300959e-02_wp
1476         ALP020 = -1.2634838399_wp
1477         ALP120 = 1.6112195176_wp
1478         ALP220 = -7.5817155402e-01_wp
1479         ALP320 = 4.7006963580e-02_wp
1480         ALP030 = 8.0812310102e-01_wp
1481         ALP130 = -1.0102374985_wp
1482         ALP230 = 4.8340368631e-01_wp
1483         ALP040 = -1.5098959754e-01_wp
1484         ALP140 = -1.4394226233e-02_wp
1485         ALP050 = 3.6780433255e-02_wp
1486         ALP001 = 3.9631611467e-01_wp
1487         ALP101 = 1.9159845880e-02_wp
1488         ALP201 = -1.0286156825e-01_wp
1489         ALP301 = 1.6738969362e-02_wp
1490         ALP011 = -4.9997430930e-01_wp
1491         ALP111 = 9.7335338937e-03_wp
1492         ALP211 = 6.0887771651e-02_wp
1493         ALP021 = 2.6149576513e-01_wp
1494         ALP121 = -1.6671866715e-02_wp
1495         ALP031 = -5.9503008642e-02_wp
1496         ALP002 = -5.4590812035e-02_wp
1497         ALP102 = 8.6134185799e-03_wp
1498         ALP012 = 6.2740815484e-02_wp
1499         ALP003 = -9.5765341718e-03_wp
1500         !
1501         BET000 = 2.1420623987_wp
1502         BET100 = -9.3752598635_wp
1503         BET200 = 1.9446303907e+01_wp
1504         BET300 = -1.8632235232e+01_wp
1505         BET400 = 8.9390837485_wp
1506         BET500 = -1.7142465871_wp
1507         BET010 = -1.7059677327e-01_wp
1508         BET110 = 2.2119589983e-01_wp
1509         BET210 = -2.7123520642e-01_wp
1510         BET310 = 7.3872053057e-02_wp
1511         BET410 = 1.2522950346e-02_wp
1512         BET020 = 3.0054390409e-01_wp
1513         BET120 = -4.2018759976e-01_wp
1514         BET220 = 3.1429341406e-01_wp
1515         BET320 = -9.8855300959e-02_wp
1516         BET030 = -2.6853658626e-01_wp
1517         BET130 = 2.5272385134e-01_wp
1518         BET230 = -2.3503481790e-02_wp
1519         BET040 = 1.2627968731e-01_wp
1520         BET140 = -1.2085092158e-01_wp
1521         BET050 = 1.4394226233e-03_wp
1522         BET001 = -2.2271304375e-01_wp
1523         BET101 = 5.5453416919e-01_wp
1524         BET201 = -6.2815936268e-01_wp
1525         BET301 = 2.0601115202e-01_wp
1526         BET011 = -9.5799229402e-03_wp
1527         BET111 = 1.0286156825e-01_wp
1528         BET211 = -2.5108454043e-02_wp
1529         BET021 = -2.4333834734e-03_wp
1530         BET121 = -3.0443885826e-02_wp
1531         BET031 = 2.7786444526e-03_wp
1532         BET002 = -4.2811838287e-02_wp
1533         BET102 = 5.1355066072e-02_wp
1534         BET012 = -4.3067092900e-03_wp
1535         BET003 = -7.1548119050e-04_wp
1536         !
1537         PEN000 = -5.3743005340_wp
1538         PEN100 = 8.9085217499_wp
1539         PEN200 = -1.1090683384e+01_wp
1540         PEN300 = 8.3754581690_wp
1541         PEN400 = -2.0601115202_wp
1542         PEN010 = 7.9263222935_wp
1543         PEN110 = 3.8319691761e-01_wp
1544         PEN210 = -2.0572313651_wp
1545         PEN310 = 3.3477938724e-01_wp
1546         PEN020 = -4.9997430930_wp
1547         PEN120 = 9.7335338937e-02_wp
1548         PEN220 = 6.0887771651e-01_wp
1549         PEN030 = 1.7433051009_wp
1550         PEN130 = -1.1114577810e-01_wp
1551         PEN040 = -2.9751504321e-01_wp
1552         PEN001 = -6.9171176978e-01_wp
1553         PEN101 = 2.2832980419_wp
1554         PEN201 = -1.3694684286_wp
1555         PEN011 = -1.4557549876_wp
1556         PEN111 = 2.2969116213e-01_wp
1557         PEN021 = 8.3654420645e-01_wp
1558         PEN002 = -1.4046808820e-02_wp
1559         PEN102 = 4.2928871430e-02_wp
1560         PEN012 = -2.8729602515e-01_wp
1561         !
1562         APE000 = -1.9815805734e-01_wp
1563         APE100 = -9.5799229402e-03_wp
1564         APE200 = 5.1430784127e-02_wp
1565         APE300 = -8.3694846809e-03_wp
1566         APE010 = 2.4998715465e-01_wp
1567         APE110 = -4.8667669469e-03_wp
1568         APE210 = -3.0443885826e-02_wp
1569         APE020 = -1.3074788257e-01_wp
1570         APE120 = 8.3359333577e-03_wp
1571         APE030 = 2.9751504321e-02_wp
1572         APE001 = 3.6393874690e-02_wp
1573         APE101 = -5.7422790533e-03_wp
1574         APE011 = -4.1827210323e-02_wp
1575         APE002 = 7.1824006288e-03_wp
1576         !
1577         BPE000 = 1.1135652187e-01_wp
1578         BPE100 = -2.7726708459e-01_wp
1579         BPE200 = 3.1407968134e-01_wp
1580         BPE300 = -1.0300557601e-01_wp
1581         BPE010 = 4.7899614701e-03_wp
1582         BPE110 = -5.1430784127e-02_wp
1583         BPE210 = 1.2554227021e-02_wp
1584         BPE020 = 1.2166917367e-03_wp
1585         BPE120 = 1.5221942913e-02_wp
1586         BPE030 = -1.3893222263e-03_wp
1587         BPE001 = 2.8541225524e-02_wp
1588         BPE101 = -3.4236710714e-02_wp
1589         BPE011 = 2.8711395266e-03_wp
1590         BPE002 = 5.3661089288e-04_wp
1591         !
1592      CASE( np_seos )                        !==  Simplified EOS     ==!
1593
1594         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
1595         
1596         IF(lwp) THEN
1597            WRITE(numout,*)
1598            WRITE(numout,*) '   ==>>>   use of simplified eos:    '
1599            WRITE(numout,*) '              rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
1600            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rau0'
1601            WRITE(numout,*) '              with the following coefficients :'
1602            WRITE(numout,*) '                 thermal exp. coef.    rn_a0      = ', rn_a0
1603            WRITE(numout,*) '                 saline  cont. coef.   rn_b0      = ', rn_b0
1604            WRITE(numout,*) '                 cabbeling coef.       rn_lambda1 = ', rn_lambda1
1605            WRITE(numout,*) '                 cabbeling coef.       rn_lambda2 = ', rn_lambda2
1606            WRITE(numout,*) '                 thermobar. coef.      rn_mu1     = ', rn_mu1
1607            WRITE(numout,*) '                 thermobar. coef.      rn_mu2     = ', rn_mu2
1608            WRITE(numout,*) '                 2nd cabbel. coef.     rn_nu      = ', rn_nu
1609            WRITE(numout,*) '              Caution: rn_beta0=0 incompatible with ddm parameterization '
1610         ENDIF
1611         l_useCT = .TRUE.          ! Use conservative temperature
1612         !
1613      CASE DEFAULT                     !==  ERROR in neos  ==!
1614         WRITE(ctmp1,*) '          bad flag value for neos = ', neos, '. You should never see this error'
1615         CALL ctl_stop( ctmp1 )
1616         !
1617      END SELECT
1618      !
1619      rau0_rcp    = rau0 * rcp 
1620      r1_rau0     = 1._wp / rau0
1621      r1_rcp      = 1._wp / rcp
1622      r1_rau0_rcp = 1._wp / rau0_rcp 
1623      !
1624      IF(lwp) THEN
1625         IF( l_useCT )   THEN
1626            WRITE(numout,*)
1627            WRITE(numout,*) '   ==>>>   model uses Conservative Temperature'
1628            WRITE(numout,*) '           Important: model must be initialized with CT and SA fields'
1629         ELSE
1630            WRITE(numout,*)
1631            WRITE(numout,*) '   ==>>>   model does not use Conservative Temperature'
1632         ENDIF
1633      ENDIF
1634      !
1635      IF(lwp) WRITE(numout,*)
1636      IF(lwp) WRITE(numout,*) '   Associated physical constant'
1637      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
1638      IF(lwp) WRITE(numout,*) '      1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
1639      IF(lwp) WRITE(numout,*) '      ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
1640      IF(lwp) WRITE(numout,*) '      rau0 * rcp                       rau0_rcp = ', rau0_rcp
1641      IF(lwp) WRITE(numout,*) '      1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
1642      !
1643   END SUBROUTINE eos_init
1644
1645   !!======================================================================
1646END MODULE eosbn2
Note: See TracBrowser for help on using the repository browser.