source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 74.4 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(inout) ::   prd    ! in situ density            [-]
315      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::   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      CALL lbc_lnk( prd, 'T', 1.0_wp )
460      !
461      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
462      !
463      IF( nn_timing == 1 )   CALL timing_stop('eos-pot')
464      !
465   END SUBROUTINE eos_insitu_pot
466
467
468   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
469      !!----------------------------------------------------------------------
470      !!                  ***  ROUTINE eos_insitu_2d  ***
471      !!
472      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
473      !!      potential temperature and salinity using an equation of state
474      !!      defined through the namelist parameter nn_eos. * 2D field case
475      !!
476      !! ** Action  : - prd , the in situ density (no units) (unmasked)
477      !!
478      !!----------------------------------------------------------------------
479      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
480      !                                                           ! 2 : salinity               [psu]
481      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
482      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
483      !
484      INTEGER  ::   ji, jj, jk                ! dummy loop indices
485      REAL(wp) ::   zt , zh , zs              ! local scalars
486      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
487      !!----------------------------------------------------------------------
488      !
489      IF( nn_timing == 1 )   CALL timing_start('eos2d')
490      !
491      prd(:,:) = 0._wp
492      !
493      SELECT CASE( nn_eos )
494      !
495      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
496         !
497         DO jj = 1, jpjm1
498            DO ji = 1, fs_jpim1   ! vector opt.
499               !
500               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
501               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
502               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
503               !
504               zn3 = EOS013*zt   &
505                  &   + EOS103*zs+EOS003
506                  !
507               zn2 = (EOS022*zt   &
508                  &   + EOS112*zs+EOS012)*zt   &
509                  &   + (EOS202*zs+EOS102)*zs+EOS002
510                  !
511               zn1 = (((EOS041*zt   &
512                  &   + EOS131*zs+EOS031)*zt   &
513                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
514                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
515                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
516                  !
517               zn0 = (((((EOS060*zt   &
518                  &   + EOS150*zs+EOS050)*zt   &
519                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
520                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
521                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
522                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
523                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
524                  !
525               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
526               !
527               prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
528               !
529            END DO
530         END DO
531         !
532         CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions
533         !
534      CASE( 1 )                !==  simplified EOS  ==!
535         !
536         DO jj = 1, jpjm1
537            DO ji = 1, fs_jpim1   ! vector opt.
538               !
539               zt    = pts  (ji,jj,jp_tem)  - 10._wp
540               zs    = pts  (ji,jj,jp_sal)  - 35._wp
541               zh    = pdep (ji,jj)                         ! depth at the partial step level
542               !
543               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
544                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
545                  &  - rn_nu * zt * zs
546                  !
547               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
548               !
549            END DO
550         END DO
551         !
552         CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions
553         !
554      END SELECT
555      !
556      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
557      !
558      IF( nn_timing == 1 )   CALL timing_stop('eos2d')
559      !
560   END SUBROUTINE eos_insitu_2d
561
562
563   SUBROUTINE rab_3d( pts, pab )
564      !!----------------------------------------------------------------------
565      !!                 ***  ROUTINE rab_3d  ***
566      !!
567      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points
568      !!
569      !! ** Method  :   calculates alpha / beta at T-points
570      !!
571      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
572      !!----------------------------------------------------------------------
573      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
574      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
575      !
576      INTEGER  ::   ji, jj, jk                ! dummy loop indices
577      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
578      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
579      !!----------------------------------------------------------------------
580      !
581      IF( nn_timing == 1 )   CALL timing_start('rab_3d')
582      !
583      SELECT CASE ( nn_eos )
584      !
585      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
586         !
587         DO jk = 1, jpkm1
588            DO jj = 1, jpj
589               DO ji = 1, jpi
590                  !
591                  zh  = fsdept(ji,jj,jk) * r1_Z0                                ! depth
592                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
593                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
594                  ztm = tmask(ji,jj,jk)                                         ! tmask
595                  !
596                  ! alpha
597                  zn3 = ALP003
598                  !
599                  zn2 = ALP012*zt + ALP102*zs+ALP002
600                  !
601                  zn1 = ((ALP031*zt   &
602                     &   + ALP121*zs+ALP021)*zt   &
603                     &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
604                     &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
605                     !
606                  zn0 = ((((ALP050*zt   &
607                     &   + ALP140*zs+ALP040)*zt   &
608                     &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
609                     &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
610                     &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
611                     &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
612                     !
613                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
614                  !
615                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
616                  !
617                  ! beta
618                  zn3 = BET003
619                  !
620                  zn2 = BET012*zt + BET102*zs+BET002
621                  !
622                  zn1 = ((BET031*zt   &
623                     &   + BET121*zs+BET021)*zt   &
624                     &   + (BET211*zs+BET111)*zs+BET011)*zt   &
625                     &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
626                     !
627                  zn0 = ((((BET050*zt   &
628                     &   + BET140*zs+BET040)*zt   &
629                     &   + (BET230*zs+BET130)*zs+BET030)*zt   &
630                     &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
631                     &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
632                     &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
633                     !
634                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
635                  !
636                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
637                  !
638               END DO
639            END DO
640         END DO
641         !
642      CASE( 1 )                  !==  simplified EOS  ==!
643         !
644         DO jk = 1, jpkm1
645            DO jj = 1, jpj
646               DO ji = 1, jpi
647                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
648                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
649                  zh  = fsdept(ji,jj,jk)                 ! depth in meters at t-point
650                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
651                  !
652                  zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
653                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
654                  !
655                  zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
656                  pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
657                  !
658               END DO
659            END DO
660         END DO
661         !
662      CASE DEFAULT
663         IF(lwp) WRITE(numout,cform_err)
664         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
665         nstop = nstop + 1
666         !
667      END SELECT
668      !
669      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
670         &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk )
671      !
672      IF( nn_timing == 1 )   CALL timing_stop('rab_3d')
673      !
674   END SUBROUTINE rab_3d
675
676   SUBROUTINE rab_2d( pts, pdep, pab )
677      !!----------------------------------------------------------------------
678      !!                 ***  ROUTINE rab_2d  ***
679      !!
680      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
681      !!
682      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
683      !!----------------------------------------------------------------------
684      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
685      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m]
686      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
687      !
688      INTEGER  ::   ji, jj, jk                ! dummy loop indices
689      REAL(wp) ::   zt , zh , zs              ! local scalars
690      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
691      !!----------------------------------------------------------------------
692      !
693      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
694      !
695      pab(:,:,:) = 0._wp
696      !
697      SELECT CASE ( nn_eos )
698      !
699      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
700         !
701         DO jj = 1, jpjm1
702            DO ji = 1, fs_jpim1   ! vector opt.
703               !
704               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
705               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
706               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
707               !
708               ! alpha
709               zn3 = ALP003
710               !
711               zn2 = ALP012*zt + ALP102*zs+ALP002
712               !
713               zn1 = ((ALP031*zt   &
714                  &   + ALP121*zs+ALP021)*zt   &
715                  &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
716                  &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
717                  !
718               zn0 = ((((ALP050*zt   &
719                  &   + ALP140*zs+ALP040)*zt   &
720                  &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
721                  &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
722                  &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
723                  &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
724                  !
725               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
726               !
727               pab(ji,jj,jp_tem) = zn * r1_rau0
728               !
729               ! beta
730               zn3 = BET003
731               !
732               zn2 = BET012*zt + BET102*zs+BET002
733               !
734               zn1 = ((BET031*zt   &
735                  &   + BET121*zs+BET021)*zt   &
736                  &   + (BET211*zs+BET111)*zs+BET011)*zt   &
737                  &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
738                  !
739               zn0 = ((((BET050*zt   &
740                  &   + BET140*zs+BET040)*zt   &
741                  &   + (BET230*zs+BET130)*zs+BET030)*zt   &
742                  &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
743                  &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
744                  &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
745                  !
746               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
747               !
748               pab(ji,jj,jp_sal) = zn / zs * r1_rau0
749               !
750               !
751            END DO
752         END DO
753         !
754         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions
755         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                   
756         !
757      CASE( 1 )                  !==  simplified EOS  ==!
758         !
759         DO jj = 1, jpjm1
760            DO ji = 1, fs_jpim1   ! vector opt.
761               !
762               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
763               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
764               zh    = pdep (ji,jj)                   ! depth at the partial step level
765               !
766               zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
767               pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
768               !
769               zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
770               pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
771               !
772            END DO
773         END DO
774         !
775         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions
776         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                   
777         !
778      CASE DEFAULT
779         IF(lwp) WRITE(numout,cform_err)
780         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
781         nstop = nstop + 1
782         !
783      END SELECT
784      !
785      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
786         &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
787      !
788      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
789      !
790   END SUBROUTINE rab_2d
791
792
793   SUBROUTINE rab_0d( pts, pdep, pab )
794      !!----------------------------------------------------------------------
795      !!                 ***  ROUTINE rab_0d  ***
796      !!
797      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
798      !!
799      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
800      !!----------------------------------------------------------------------
801      REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
802      REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m]
803      REAL(wp), DIMENSION(jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
804      !
805      REAL(wp) ::   zt , zh , zs              ! local scalars
806      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
807      !!----------------------------------------------------------------------
808      !
809      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
810      !
811      pab(:) = 0._wp
812      !
813      SELECT CASE ( nn_eos )
814      !
815      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
816         !
817         !
818         zh  = pdep * r1_Z0                                  ! depth
819         zt  = pts (jp_tem) * r1_T0                           ! temperature
820         zs  = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
821         !
822         ! alpha
823         zn3 = ALP003
824         !
825         zn2 = ALP012*zt + ALP102*zs+ALP002
826         !
827         zn1 = ((ALP031*zt   &
828            &   + ALP121*zs+ALP021)*zt   &
829            &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
830            &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
831            !
832         zn0 = ((((ALP050*zt   &
833            &   + ALP140*zs+ALP040)*zt   &
834            &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
835            &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
836            &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
837            &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
838            !
839         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
840         !
841         pab(jp_tem) = zn * r1_rau0
842         !
843         ! beta
844         zn3 = BET003
845         !
846         zn2 = BET012*zt + BET102*zs+BET002
847         !
848         zn1 = ((BET031*zt   &
849            &   + BET121*zs+BET021)*zt   &
850            &   + (BET211*zs+BET111)*zs+BET011)*zt   &
851            &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
852            !
853         zn0 = ((((BET050*zt   &
854            &   + BET140*zs+BET040)*zt   &
855            &   + (BET230*zs+BET130)*zs+BET030)*zt   &
856            &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
857            &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
858            &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
859            !
860         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
861         !
862         pab(jp_sal) = zn / zs * r1_rau0
863         !
864         !
865         !
866      CASE( 1 )                  !==  simplified EOS  ==!
867         !
868         zt    = pts(jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
869         zs    = pts(jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
870         zh    = pdep                    ! depth at the partial step level
871         !
872         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
873         pab(jp_tem) = zn * r1_rau0   ! alpha
874         !
875         zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
876         pab(jp_sal) = zn * r1_rau0   ! beta
877         !
878      CASE DEFAULT
879         IF(lwp) WRITE(numout,cform_err)
880         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
881         nstop = nstop + 1
882         !
883      END SELECT
884      !
885      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
886      !
887   END SUBROUTINE rab_0d
888
889
890   SUBROUTINE bn2( pts, pab, pn2 )
891      !!----------------------------------------------------------------------
892      !!                  ***  ROUTINE bn2  ***
893      !!
894      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the
895      !!                time-step of the input arguments
896      !!
897      !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
898      !!      where alpha and beta are given in pab, and computed on T-points.
899      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module.
900      !!
901      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point
902      !!
903      !!----------------------------------------------------------------------
904      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu]
905      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1]
906      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
907      !
908      INTEGER  ::   ji, jj, jk      ! dummy loop indices
909      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
910      !!----------------------------------------------------------------------
911      !
912      IF( nn_timing == 1 ) CALL timing_start('bn2')
913      !
914      DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
915         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
916            DO ji = 1, jpi
917               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
918                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
919                  !
920               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
921               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
922               !
923               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
924                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
925                  &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
926            END DO
927         END DO
928      END DO
929      !
930      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
931      !
932      IF( nn_timing == 1 )   CALL timing_stop('bn2')
933      !
934   END SUBROUTINE bn2
935
936
937   FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
938      !!----------------------------------------------------------------------
939      !!                 ***  ROUTINE eos_pt_from_ct  ***
940      !!
941      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celcius]
942      !!
943      !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
944      !!       checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
945      !!
946      !! Reference  :   TEOS-10, UNESCO
947      !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
948      !!----------------------------------------------------------------------
949      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp [Celcius]
950      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity   [psu]
951      ! Leave result array automatic rather than making explicitly allocated
952      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celcius]
953      !
954      INTEGER  ::   ji, jj               ! dummy loop indices
955      REAL(wp) ::   zt , zs , ztm        ! local scalars
956      REAL(wp) ::   zn , zd              ! local scalars
957      REAL(wp) ::   zdeltaS , z1_S0 , z1_T0
958      !!----------------------------------------------------------------------
959      !
960      IF ( nn_timing == 1 )   CALL timing_start('eos_pt_from_ct')
961      !
962      zdeltaS = 5._wp
963      z1_S0   = 0.875_wp/35.16504_wp
964      z1_T0   = 1._wp/40._wp
965      !
966      DO jj = 1, jpj
967         DO ji = 1, jpi
968            !
969            zt  = ctmp   (ji,jj) * z1_T0
970            zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
971            ztm = tmask(ji,jj,1)
972            !
973            zn = ((((-2.1385727895e-01_wp*zt   &
974               &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
975               &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
976               &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
977               &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
978               &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
979               &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
980               &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
981               !
982            zd = (2.0035003456_wp*zt   &
983               &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
984               &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
985               !
986            ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
987               !
988         END DO
989      END DO
990      !
991      IF( nn_timing == 1 )   CALL timing_stop('eos_pt_from_ct')
992      !
993   END FUNCTION eos_pt_from_ct
994
995
996   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep )
997      !!----------------------------------------------------------------------
998      !!                 ***  ROUTINE eos_fzp  ***
999      !!
1000      !! ** Purpose :   Compute the freezing point temperature [Celcius]
1001      !!
1002      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by
1003      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
1004      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
1005      !!
1006      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
1007      !!----------------------------------------------------------------------
1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
1009      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
1010      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celcius]
1011      !
1012      INTEGER  ::   ji, jj   ! dummy loop indices
1013      REAL(wp) ::   zt, zs   ! local scalars
1014      !!----------------------------------------------------------------------
1015      !
1016      SELECT CASE ( nn_eos )
1017      !
1018      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
1019         !
1020         DO jj = 1, jpj
1021            DO ji = 1, jpi
1022               zs= SQRT( ABS( psal(ji,jj) ) / 35.16504_wp )           ! square root salinity
1023               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
1024                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
1025            END DO
1026         END DO
1027         ptf(:,:) = ptf(:,:) * psal(:,:)
1028         !
1029         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
1030         !
1031      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
1032         !
1033         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
1034            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
1035            !
1036         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
1037         !
1038      CASE DEFAULT
1039         IF(lwp) WRITE(numout,cform_err)
1040         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
1041         nstop = nstop + 1
1042         !
1043      END SELECT     
1044      !
1045  END SUBROUTINE eos_fzp_2d
1046
1047  SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
1048      !!----------------------------------------------------------------------
1049      !!                 ***  ROUTINE eos_fzp  ***
1050      !!
1051      !! ** Purpose :   Compute the freezing point temperature [Celcius]
1052      !!
1053      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by
1054      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
1055      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
1056      !!
1057      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
1058      !!----------------------------------------------------------------------
1059      REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu]
1060      REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m]
1061      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celcius]
1062      !
1063      REAL(wp) :: zs   ! local scalars
1064      !!----------------------------------------------------------------------
1065      !
1066      SELECT CASE ( nn_eos )
1067      !
1068      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
1069         !
1070         zs  = SQRT( ABS( psal ) / 35.16504_wp )           ! square root salinity
1071         ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
1072                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
1073         ptf = ptf * psal
1074         !
1075         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
1076         !
1077      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
1078         !
1079         ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   &
1080            &                - 2.154996e-4_wp *       psal   ) * psal
1081            !
1082         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
1083         !
1084      CASE DEFAULT
1085         IF(lwp) WRITE(numout,cform_err)
1086         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
1087         nstop = nstop + 1
1088         !
1089      END SELECT
1090      !
1091   END SUBROUTINE eos_fzp_0d
1092
1093
1094   SUBROUTINE eos_pen( pts, pab_pe, ppen )
1095      !!----------------------------------------------------------------------
1096      !!                 ***  ROUTINE eos_pen  ***
1097      !!
1098      !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
1099      !!
1100      !! ** Method  :   PE is defined analytically as the vertical
1101      !!                   primitive of EOS times -g integrated between 0 and z>0.
1102      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
1103      !!                                                      = 1/z * /int_0^z rd dz - rd
1104      !!                                where rd is the density anomaly (see eos_rhd function)
1105      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
1106      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
1107      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS
1108      !!
1109      !! ** Action  : - pen         : PE anomaly given at T-points
1110      !!            : - pab_pe  : given at T-points
1111      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe
1112      !!                    pab_pe(:,:,:,jp_sal) is beta_pe
1113      !!----------------------------------------------------------------------
1114      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
1115      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
1116      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   ppen     ! potential energy anomaly
1117      !
1118      INTEGER  ::   ji, jj, jk                ! dummy loop indices
1119      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
1120      REAL(wp) ::   zn , zn0, zn1, zn2        !   -      -
1121      !!----------------------------------------------------------------------
1122      !
1123      IF( nn_timing == 1 )   CALL timing_start('eos_pen')
1124      !
1125      SELECT CASE ( nn_eos )
1126      !
1127      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
1128         !
1129         DO jk = 1, jpkm1
1130            DO jj = 1, jpj
1131               DO ji = 1, jpi
1132                  !
1133                  zh  = fsdept(ji,jj,jk) * r1_Z0                                ! depth
1134                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
1135                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
1136                  ztm = tmask(ji,jj,jk)                                         ! tmask
1137                  !
1138                  ! potential energy non-linear anomaly
1139                  zn2 = (PEN012)*zt   &
1140                     &   + PEN102*zs+PEN002
1141                     !
1142                  zn1 = ((PEN021)*zt   &
1143                     &   + PEN111*zs+PEN011)*zt   &
1144                     &   + (PEN201*zs+PEN101)*zs+PEN001
1145                     !
1146                  zn0 = ((((PEN040)*zt   &
1147                     &   + PEN130*zs+PEN030)*zt   &
1148                     &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
1149                     &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
1150                     &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
1151                     !
1152                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
1153                  !
1154                  ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm
1155                  !
1156                  ! alphaPE non-linear anomaly
1157                  zn2 = APE002
1158                  !
1159                  zn1 = (APE011)*zt   &
1160                     &   + APE101*zs+APE001
1161                     !
1162                  zn0 = (((APE030)*zt   &
1163                     &   + APE120*zs+APE020)*zt   &
1164                     &   + (APE210*zs+APE110)*zs+APE010)*zt   &
1165                     &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
1166                     !
1167                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
1168                  !                             
1169                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
1170                  !
1171                  ! betaPE non-linear anomaly
1172                  zn2 = BPE002
1173                  !
1174                  zn1 = (BPE011)*zt   &
1175                     &   + BPE101*zs+BPE001
1176                     !
1177                  zn0 = (((BPE030)*zt   &
1178                     &   + BPE120*zs+BPE020)*zt   &
1179                     &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
1180                     &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
1181                     !
1182                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
1183                  !                             
1184                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
1185                  !
1186               END DO
1187            END DO
1188         END DO
1189         !
1190      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
1191         !
1192         DO jk = 1, jpkm1
1193            DO jj = 1, jpj
1194               DO ji = 1, jpi
1195                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
1196                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
1197                  zh  = fsdept(ji,jj,jk)               ! depth in meters  at t-point
1198                  ztm = tmask(ji,jj,jk)                ! tmask
1199                  zn  = 0.5_wp * zh * r1_rau0 * ztm
1200                  !                                    ! Potential Energy
1201                  ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
1202                  !                                    ! alphaPE
1203                  pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
1204                  pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
1205                  !
1206               END DO
1207            END DO
1208         END DO
1209         !
1210      CASE DEFAULT
1211         IF(lwp) WRITE(numout,cform_err)
1212         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
1213         nstop = nstop + 1
1214         !
1215      END SELECT
1216      !
1217      IF( nn_timing == 1 )   CALL timing_stop('eos_pen')
1218      !
1219   END SUBROUTINE eos_pen
1220
1221
1222   SUBROUTINE eos_init
1223      !!----------------------------------------------------------------------
1224      !!                 ***  ROUTINE eos_init  ***
1225      !!
1226      !! ** Purpose :   initializations for the equation of state
1227      !!
1228      !! ** Method  :   Read the namelist nameos and control the parameters
1229      !!----------------------------------------------------------------------
1230      INTEGER  ::   ios   ! local integer
1231      !!
1232      NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
1233         &                                             rn_lambda2, rn_mu2, rn_nu
1234      !!----------------------------------------------------------------------
1235      !
1236      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
1237      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
1238901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
1239      !
1240      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
1241      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
1242902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
1243      IF(lwm .AND. nprint > 2) WRITE( numond, nameos )
1244      !
1245      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
1246      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
1247      !
1248      IF(lwp) THEN                ! Control print
1249         WRITE(numout,*)
1250         WRITE(numout,*) 'eos_init : equation of state'
1251         WRITE(numout,*) '~~~~~~~~'
1252         WRITE(numout,*) '          Namelist nameos : set eos parameters'
1253         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
1254         IF( ln_useCT )   THEN
1255            WRITE(numout,*) '             model uses Conservative Temperature'
1256            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
1257         ELSE
1258            WRITE(numout,*) '             model does not use Conservative Temperature'
1259         ENDIF
1260         IF(lflush) CALL flush(numout)
1261      ENDIF
1262      !
1263      ! Consistency check on ln_useCT and nn_eos
1264      IF ((nn_eos .EQ. -1) .AND. (.NOT. ln_useCT)) THEN
1265         CALL ctl_stop("ln_useCT should be set to True if using TEOS-10 (nn_eos=-1)")
1266      ELSE IF ((nn_eos .NE. -1) .AND. (ln_useCT)) THEN
1267         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)")
1268      ENDIF
1269      !
1270      SELECT CASE( nn_eos )         ! check option
1271      !
1272      CASE( -1 )                       !==  polynomial TEOS-10  ==!
1273         IF(lwp) WRITE(numout,*)
1274         IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
1275         IF(lwp .AND. lflush) CALL flush(numout)
1276         !
1277         rdeltaS = 32._wp
1278         r1_S0  = 0.875_wp/35.16504_wp
1279         r1_T0  = 1._wp/40._wp
1280         r1_Z0  = 1.e-4_wp
1281         !
1282         EOS000 = 8.0189615746e+02_wp
1283         EOS100 = 8.6672408165e+02_wp
1284         EOS200 = -1.7864682637e+03_wp
1285         EOS300 = 2.0375295546e+03_wp
1286         EOS400 = -1.2849161071e+03_wp
1287         EOS500 = 4.3227585684e+02_wp
1288         EOS600 = -6.0579916612e+01_wp
1289         EOS010 = 2.6010145068e+01_wp
1290         EOS110 = -6.5281885265e+01_wp
1291         EOS210 = 8.1770425108e+01_wp
1292         EOS310 = -5.6888046321e+01_wp
1293         EOS410 = 1.7681814114e+01_wp
1294         EOS510 = -1.9193502195_wp
1295         EOS020 = -3.7074170417e+01_wp
1296         EOS120 = 6.1548258127e+01_wp
1297         EOS220 = -6.0362551501e+01_wp
1298         EOS320 = 2.9130021253e+01_wp
1299         EOS420 = -5.4723692739_wp
1300         EOS030 = 2.1661789529e+01_wp
1301         EOS130 = -3.3449108469e+01_wp
1302         EOS230 = 1.9717078466e+01_wp
1303         EOS330 = -3.1742946532_wp
1304         EOS040 = -8.3627885467_wp
1305         EOS140 = 1.1311538584e+01_wp
1306         EOS240 = -5.3563304045_wp
1307         EOS050 = 5.4048723791e-01_wp
1308         EOS150 = 4.8169980163e-01_wp
1309         EOS060 = -1.9083568888e-01_wp
1310         EOS001 = 1.9681925209e+01_wp
1311         EOS101 = -4.2549998214e+01_wp
1312         EOS201 = 5.0774768218e+01_wp
1313         EOS301 = -3.0938076334e+01_wp
1314         EOS401 = 6.6051753097_wp
1315         EOS011 = -1.3336301113e+01_wp
1316         EOS111 = -4.4870114575_wp
1317         EOS211 = 5.0042598061_wp
1318         EOS311 = -6.5399043664e-01_wp
1319         EOS021 = 6.7080479603_wp
1320         EOS121 = 3.5063081279_wp
1321         EOS221 = -1.8795372996_wp
1322         EOS031 = -2.4649669534_wp
1323         EOS131 = -5.5077101279e-01_wp
1324         EOS041 = 5.5927935970e-01_wp
1325         EOS002 = 2.0660924175_wp
1326         EOS102 = -4.9527603989_wp
1327         EOS202 = 2.5019633244_wp
1328         EOS012 = 2.0564311499_wp
1329         EOS112 = -2.1311365518e-01_wp
1330         EOS022 = -1.2419983026_wp
1331         EOS003 = -2.3342758797e-02_wp
1332         EOS103 = -1.8507636718e-02_wp
1333         EOS013 = 3.7969820455e-01_wp
1334         !
1335         ALP000 = -6.5025362670e-01_wp
1336         ALP100 = 1.6320471316_wp
1337         ALP200 = -2.0442606277_wp
1338         ALP300 = 1.4222011580_wp
1339         ALP400 = -4.4204535284e-01_wp
1340         ALP500 = 4.7983755487e-02_wp
1341         ALP010 = 1.8537085209_wp
1342         ALP110 = -3.0774129064_wp
1343         ALP210 = 3.0181275751_wp
1344         ALP310 = -1.4565010626_wp
1345         ALP410 = 2.7361846370e-01_wp
1346         ALP020 = -1.6246342147_wp
1347         ALP120 = 2.5086831352_wp
1348         ALP220 = -1.4787808849_wp
1349         ALP320 = 2.3807209899e-01_wp
1350         ALP030 = 8.3627885467e-01_wp
1351         ALP130 = -1.1311538584_wp
1352         ALP230 = 5.3563304045e-01_wp
1353         ALP040 = -6.7560904739e-02_wp
1354         ALP140 = -6.0212475204e-02_wp
1355         ALP050 = 2.8625353333e-02_wp
1356         ALP001 = 3.3340752782e-01_wp
1357         ALP101 = 1.1217528644e-01_wp
1358         ALP201 = -1.2510649515e-01_wp
1359         ALP301 = 1.6349760916e-02_wp
1360         ALP011 = -3.3540239802e-01_wp
1361         ALP111 = -1.7531540640e-01_wp
1362         ALP211 = 9.3976864981e-02_wp
1363         ALP021 = 1.8487252150e-01_wp
1364         ALP121 = 4.1307825959e-02_wp
1365         ALP031 = -5.5927935970e-02_wp
1366         ALP002 = -5.1410778748e-02_wp
1367         ALP102 = 5.3278413794e-03_wp
1368         ALP012 = 6.2099915132e-02_wp
1369         ALP003 = -9.4924551138e-03_wp
1370         !
1371         BET000 = 1.0783203594e+01_wp
1372         BET100 = -4.4452095908e+01_wp
1373         BET200 = 7.6048755820e+01_wp
1374         BET300 = -6.3944280668e+01_wp
1375         BET400 = 2.6890441098e+01_wp
1376         BET500 = -4.5221697773_wp
1377         BET010 = -8.1219372432e-01_wp
1378         BET110 = 2.0346663041_wp
1379         BET210 = -2.1232895170_wp
1380         BET310 = 8.7994140485e-01_wp
1381         BET410 = -1.1939638360e-01_wp
1382         BET020 = 7.6574242289e-01_wp
1383         BET120 = -1.5019813020_wp
1384         BET220 = 1.0872489522_wp
1385         BET320 = -2.7233429080e-01_wp
1386         BET030 = -4.1615152308e-01_wp
1387         BET130 = 4.9061350869e-01_wp
1388         BET230 = -1.1847737788e-01_wp
1389         BET040 = 1.4073062708e-01_wp
1390         BET140 = -1.3327978879e-01_wp
1391         BET050 = 5.9929880134e-03_wp
1392         BET001 = -5.2937873009e-01_wp
1393         BET101 = 1.2634116779_wp
1394         BET201 = -1.1547328025_wp
1395         BET301 = 3.2870876279e-01_wp
1396         BET011 = -5.5824407214e-02_wp
1397         BET111 = 1.2451933313e-01_wp
1398         BET211 = -2.4409539932e-02_wp
1399         BET021 = 4.3623149752e-02_wp
1400         BET121 = -4.6767901790e-02_wp
1401         BET031 = -6.8523260060e-03_wp
1402         BET002 = -6.1618945251e-02_wp
1403         BET102 = 6.2255521644e-02_wp
1404         BET012 = -2.6514181169e-03_wp
1405         BET003 = -2.3025968587e-04_wp
1406         !
1407         PEN000 = -9.8409626043_wp
1408         PEN100 = 2.1274999107e+01_wp
1409         PEN200 = -2.5387384109e+01_wp
1410         PEN300 = 1.5469038167e+01_wp
1411         PEN400 = -3.3025876549_wp
1412         PEN010 = 6.6681505563_wp
1413         PEN110 = 2.2435057288_wp
1414         PEN210 = -2.5021299030_wp
1415         PEN310 = 3.2699521832e-01_wp
1416         PEN020 = -3.3540239802_wp
1417         PEN120 = -1.7531540640_wp
1418         PEN220 = 9.3976864981e-01_wp
1419         PEN030 = 1.2324834767_wp
1420         PEN130 = 2.7538550639e-01_wp
1421         PEN040 = -2.7963967985e-01_wp
1422         PEN001 = -1.3773949450_wp
1423         PEN101 = 3.3018402659_wp
1424         PEN201 = -1.6679755496_wp
1425         PEN011 = -1.3709540999_wp
1426         PEN111 = 1.4207577012e-01_wp
1427         PEN021 = 8.2799886843e-01_wp
1428         PEN002 = 1.7507069098e-02_wp
1429         PEN102 = 1.3880727538e-02_wp
1430         PEN012 = -2.8477365341e-01_wp
1431         !
1432         APE000 = -1.6670376391e-01_wp
1433         APE100 = -5.6087643219e-02_wp
1434         APE200 = 6.2553247576e-02_wp
1435         APE300 = -8.1748804580e-03_wp
1436         APE010 = 1.6770119901e-01_wp
1437         APE110 = 8.7657703198e-02_wp
1438         APE210 = -4.6988432490e-02_wp
1439         APE020 = -9.2436260751e-02_wp
1440         APE120 = -2.0653912979e-02_wp
1441         APE030 = 2.7963967985e-02_wp
1442         APE001 = 3.4273852498e-02_wp
1443         APE101 = -3.5518942529e-03_wp
1444         APE011 = -4.1399943421e-02_wp
1445         APE002 = 7.1193413354e-03_wp
1446         !
1447         BPE000 = 2.6468936504e-01_wp
1448         BPE100 = -6.3170583896e-01_wp
1449         BPE200 = 5.7736640125e-01_wp
1450         BPE300 = -1.6435438140e-01_wp
1451         BPE010 = 2.7912203607e-02_wp
1452         BPE110 = -6.2259666565e-02_wp
1453         BPE210 = 1.2204769966e-02_wp
1454         BPE020 = -2.1811574876e-02_wp
1455         BPE120 = 2.3383950895e-02_wp
1456         BPE030 = 3.4261630030e-03_wp
1457         BPE001 = 4.1079296834e-02_wp
1458         BPE101 = -4.1503681096e-02_wp
1459         BPE011 = 1.7676120780e-03_wp
1460         BPE002 = 1.7269476440e-04_wp
1461         !
1462      CASE( 0 )                        !==  polynomial EOS-80 formulation  ==!
1463         !
1464         IF(lwp) WRITE(numout,*)
1465         IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)'
1466         IF(lwp .AND. lflush) CALL flush(numout)
1467         !
1468         rdeltaS = 20._wp
1469         r1_S0  = 1._wp/40._wp
1470         r1_T0  = 1._wp/40._wp
1471         r1_Z0  = 1.e-4_wp
1472         !
1473         EOS000 = 9.5356891948e+02_wp
1474         EOS100 = 1.7136499189e+02_wp
1475         EOS200 = -3.7501039454e+02_wp
1476         EOS300 = 5.1856810420e+02_wp
1477         EOS400 = -3.7264470465e+02_wp
1478         EOS500 = 1.4302533998e+02_wp
1479         EOS600 = -2.2856621162e+01_wp
1480         EOS010 = 1.0087518651e+01_wp
1481         EOS110 = -1.3647741861e+01_wp
1482         EOS210 = 8.8478359933_wp
1483         EOS310 = -7.2329388377_wp
1484         EOS410 = 1.4774410611_wp
1485         EOS510 = 2.0036720553e-01_wp
1486         EOS020 = -2.5579830599e+01_wp
1487         EOS120 = 2.4043512327e+01_wp
1488         EOS220 = -1.6807503990e+01_wp
1489         EOS320 = 8.3811577084_wp
1490         EOS420 = -1.9771060192_wp
1491         EOS030 = 1.6846451198e+01_wp
1492         EOS130 = -2.1482926901e+01_wp
1493         EOS230 = 1.0108954054e+01_wp
1494         EOS330 = -6.2675951440e-01_wp
1495         EOS040 = -8.0812310102_wp
1496         EOS140 = 1.0102374985e+01_wp
1497         EOS240 = -4.8340368631_wp
1498         EOS050 = 1.2079167803_wp
1499         EOS150 = 1.1515380987e-01_wp
1500         EOS060 = -2.4520288837e-01_wp
1501         EOS001 = 1.0748601068e+01_wp
1502         EOS101 = -1.7817043500e+01_wp
1503         EOS201 = 2.2181366768e+01_wp
1504         EOS301 = -1.6750916338e+01_wp
1505         EOS401 = 4.1202230403_wp
1506         EOS011 = -1.5852644587e+01_wp
1507         EOS111 = -7.6639383522e-01_wp
1508         EOS211 = 4.1144627302_wp
1509         EOS311 = -6.6955877448e-01_wp
1510         EOS021 = 9.9994861860_wp
1511         EOS121 = -1.9467067787e-01_wp
1512         EOS221 = -1.2177554330_wp
1513         EOS031 = -3.4866102017_wp
1514         EOS131 = 2.2229155620e-01_wp
1515         EOS041 = 5.9503008642e-01_wp
1516         EOS002 = 1.0375676547_wp
1517         EOS102 = -3.4249470629_wp
1518         EOS202 = 2.0542026429_wp
1519         EOS012 = 2.1836324814_wp
1520         EOS112 = -3.4453674320e-01_wp
1521         EOS022 = -1.2548163097_wp
1522         EOS003 = 1.8729078427e-02_wp
1523         EOS103 = -5.7238495240e-02_wp
1524         EOS013 = 3.8306136687e-01_wp
1525         !
1526         ALP000 = -2.5218796628e-01_wp
1527         ALP100 = 3.4119354654e-01_wp
1528         ALP200 = -2.2119589983e-01_wp
1529         ALP300 = 1.8082347094e-01_wp
1530         ALP400 = -3.6936026529e-02_wp
1531         ALP500 = -5.0091801383e-03_wp
1532         ALP010 = 1.2789915300_wp
1533         ALP110 = -1.2021756164_wp
1534         ALP210 = 8.4037519952e-01_wp
1535         ALP310 = -4.1905788542e-01_wp
1536         ALP410 = 9.8855300959e-02_wp
1537         ALP020 = -1.2634838399_wp
1538         ALP120 = 1.6112195176_wp
1539         ALP220 = -7.5817155402e-01_wp
1540         ALP320 = 4.7006963580e-02_wp
1541         ALP030 = 8.0812310102e-01_wp
1542         ALP130 = -1.0102374985_wp
1543         ALP230 = 4.8340368631e-01_wp
1544         ALP040 = -1.5098959754e-01_wp
1545         ALP140 = -1.4394226233e-02_wp
1546         ALP050 = 3.6780433255e-02_wp
1547         ALP001 = 3.9631611467e-01_wp
1548         ALP101 = 1.9159845880e-02_wp
1549         ALP201 = -1.0286156825e-01_wp
1550         ALP301 = 1.6738969362e-02_wp
1551         ALP011 = -4.9997430930e-01_wp
1552         ALP111 = 9.7335338937e-03_wp
1553         ALP211 = 6.0887771651e-02_wp
1554         ALP021 = 2.6149576513e-01_wp
1555         ALP121 = -1.6671866715e-02_wp
1556         ALP031 = -5.9503008642e-02_wp
1557         ALP002 = -5.4590812035e-02_wp
1558         ALP102 = 8.6134185799e-03_wp
1559         ALP012 = 6.2740815484e-02_wp
1560         ALP003 = -9.5765341718e-03_wp
1561         !
1562         BET000 = 2.1420623987_wp
1563         BET100 = -9.3752598635_wp
1564         BET200 = 1.9446303907e+01_wp
1565         BET300 = -1.8632235232e+01_wp
1566         BET400 = 8.9390837485_wp
1567         BET500 = -1.7142465871_wp
1568         BET010 = -1.7059677327e-01_wp
1569         BET110 = 2.2119589983e-01_wp
1570         BET210 = -2.7123520642e-01_wp
1571         BET310 = 7.3872053057e-02_wp
1572         BET410 = 1.2522950346e-02_wp
1573         BET020 = 3.0054390409e-01_wp
1574         BET120 = -4.2018759976e-01_wp
1575         BET220 = 3.1429341406e-01_wp
1576         BET320 = -9.8855300959e-02_wp
1577         BET030 = -2.6853658626e-01_wp
1578         BET130 = 2.5272385134e-01_wp
1579         BET230 = -2.3503481790e-02_wp
1580         BET040 = 1.2627968731e-01_wp
1581         BET140 = -1.2085092158e-01_wp
1582         BET050 = 1.4394226233e-03_wp
1583         BET001 = -2.2271304375e-01_wp
1584         BET101 = 5.5453416919e-01_wp
1585         BET201 = -6.2815936268e-01_wp
1586         BET301 = 2.0601115202e-01_wp
1587         BET011 = -9.5799229402e-03_wp
1588         BET111 = 1.0286156825e-01_wp
1589         BET211 = -2.5108454043e-02_wp
1590         BET021 = -2.4333834734e-03_wp
1591         BET121 = -3.0443885826e-02_wp
1592         BET031 = 2.7786444526e-03_wp
1593         BET002 = -4.2811838287e-02_wp
1594         BET102 = 5.1355066072e-02_wp
1595         BET012 = -4.3067092900e-03_wp
1596         BET003 = -7.1548119050e-04_wp
1597         !
1598         PEN000 = -5.3743005340_wp
1599         PEN100 = 8.9085217499_wp
1600         PEN200 = -1.1090683384e+01_wp
1601         PEN300 = 8.3754581690_wp
1602         PEN400 = -2.0601115202_wp
1603         PEN010 = 7.9263222935_wp
1604         PEN110 = 3.8319691761e-01_wp
1605         PEN210 = -2.0572313651_wp
1606         PEN310 = 3.3477938724e-01_wp
1607         PEN020 = -4.9997430930_wp
1608         PEN120 = 9.7335338937e-02_wp
1609         PEN220 = 6.0887771651e-01_wp
1610         PEN030 = 1.7433051009_wp
1611         PEN130 = -1.1114577810e-01_wp
1612         PEN040 = -2.9751504321e-01_wp
1613         PEN001 = -6.9171176978e-01_wp
1614         PEN101 = 2.2832980419_wp
1615         PEN201 = -1.3694684286_wp
1616         PEN011 = -1.4557549876_wp
1617         PEN111 = 2.2969116213e-01_wp
1618         PEN021 = 8.3654420645e-01_wp
1619         PEN002 = -1.4046808820e-02_wp
1620         PEN102 = 4.2928871430e-02_wp
1621         PEN012 = -2.8729602515e-01_wp
1622         !
1623         APE000 = -1.9815805734e-01_wp
1624         APE100 = -9.5799229402e-03_wp
1625         APE200 = 5.1430784127e-02_wp
1626         APE300 = -8.3694846809e-03_wp
1627         APE010 = 2.4998715465e-01_wp
1628         APE110 = -4.8667669469e-03_wp
1629         APE210 = -3.0443885826e-02_wp
1630         APE020 = -1.3074788257e-01_wp
1631         APE120 = 8.3359333577e-03_wp
1632         APE030 = 2.9751504321e-02_wp
1633         APE001 = 3.6393874690e-02_wp
1634         APE101 = -5.7422790533e-03_wp
1635         APE011 = -4.1827210323e-02_wp
1636         APE002 = 7.1824006288e-03_wp
1637         !
1638         BPE000 = 1.1135652187e-01_wp
1639         BPE100 = -2.7726708459e-01_wp
1640         BPE200 = 3.1407968134e-01_wp
1641         BPE300 = -1.0300557601e-01_wp
1642         BPE010 = 4.7899614701e-03_wp
1643         BPE110 = -5.1430784127e-02_wp
1644         BPE210 = 1.2554227021e-02_wp
1645         BPE020 = 1.2166917367e-03_wp
1646         BPE120 = 1.5221942913e-02_wp
1647         BPE030 = -1.3893222263e-03_wp
1648         BPE001 = 2.8541225524e-02_wp
1649         BPE101 = -3.4236710714e-02_wp
1650         BPE011 = 2.8711395266e-03_wp
1651         BPE002 = 5.3661089288e-04_wp
1652         !
1653      CASE( 1 )                        !==  Simplified EOS     ==!
1654         IF(lwp) THEN
1655            WRITE(numout,*)
1656            WRITE(numout,*) '          use of simplified eos:    rhd(dT=T-10,dS=S-35,Z) = '
1657            WRITE(numout,*) '             [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0'
1658            WRITE(numout,*)
1659            WRITE(numout,*) '             thermal exp. coef.    rn_a0      = ', rn_a0
1660            WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0
1661            WRITE(numout,*) '             cabbeling coef.       rn_lambda1 = ', rn_lambda1
1662            WRITE(numout,*) '             cabbeling coef.       rn_lambda2 = ', rn_lambda2
1663            WRITE(numout,*) '             thermobar. coef.      rn_mu1     = ', rn_mu1
1664            WRITE(numout,*) '             thermobar. coef.      rn_mu2     = ', rn_mu2
1665            WRITE(numout,*) '             2nd cabbel. coef.     rn_nu      = ', rn_nu
1666            WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization '
1667            IF(lflush) CALL flush(numout)
1668         ENDIF
1669         !
1670      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
1671         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
1672         CALL ctl_stop( ctmp1 )
1673         !
1674      END SELECT
1675      !
1676      rau0_rcp    = rau0 * rcp 
1677      r1_rau0     = 1._wp / rau0
1678      r1_rcp      = 1._wp / rcp
1679      r1_rau0_rcp = 1._wp / rau0_rcp 
1680      !
1681      IF(lwp) WRITE(numout,*)
1682      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
1683      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
1684      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
1685      IF(lwp) WRITE(numout,*) '          rau0 * rcp                       rau0_rcp = ', rau0_rcp
1686      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
1687      IF(lwp .AND. lflush) CALL flush(numout)
1688      !
1689   END SUBROUTINE eos_init
1690
1691   !!======================================================================
1692END MODULE eosbn2
Note: See TracBrowser for help on using the repository browser.