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/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 8485

Last change on this file since 8485 was 7993, checked in by frrh, 7 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.

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