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 NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/TRA/eosbn2.F90 @ 11954

Last change on this file since 11954 was 11671, checked in by acc, 5 years ago

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