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/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90 @ 15066

Last change on this file since 15066 was 15066, checked in by sparonuz, 12 months ago

Added interface in eosbn2.F90 to ease the weight of casts necessary in mixed precision + Removed unnecessary casts

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