New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
eosbn2.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
  • Property svn:keywords set to Id
File size: 44.4 KB
Line 
1MODULE eosbn2
2   !!==============================================================================
3   !!                       ***  MODULE  eosbn2  ***
4   !! Ocean diagnostic variable : equation of state - in situ and potential density
5   !!                                               - Brunt-Vaisala frequency
6   !!==============================================================================
7   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
8   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
9   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
10   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
11   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
12   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
13   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
14   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
15   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
16   !!             -   ! 2003-08  (G. Madec)  F90, free form
17   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function
18   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   eos            : generic interface of the equation of state
24   !!   eos_insitu     : Compute the in situ density
25   !!   eos_insitu_pot : Compute the insitu and surface referenced potential
26   !!                    volumic mass
27   !!   eos_insitu_2d  : Compute the in situ density for 2d fields
28   !!   eos_bn2        : Compute the Brunt-Vaisala frequency
29   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff.
30   !!   tfreez         : Compute the surface freezing temperature
31   !!   eos_init       : set eos parameters (namelist)
32   !!----------------------------------------------------------------------
33   USE dom_oce         ! ocean space and time domain
34   USE phycst          ! physical constants
35   USE zdfddm          ! vertical physics: double diffusion
36   USE in_out_manager  ! I/O manager
37   USE lib_mpp         ! MPP library
38   USE prtctl          ! Print control
39
40   IMPLICIT NONE
41   PRIVATE
42
43   !                   !! * Interface
44   INTERFACE eos
45      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
46   END INTERFACE
47   INTERFACE bn2
48      MODULE PROCEDURE eos_bn2
49   END INTERFACE
50
51   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules
52   PUBLIC   eos_init   ! called by istate module
53   PUBLIC   bn2        ! called by step module
54   PUBLIC   eos_alpbet ! called by ldfslp module
55   PUBLIC   tfreez     ! called by sbcice_... modules
56
57   !                                          !!* Namelist (nameos) *
58   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
59   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state)
60   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state)
61
62   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio
63
64   !! * Control permutation of array indices
65#  include "dom_oce_ftrans.h90"
66#  include "zdfddm_ftrans.h90"
67   
68   !! * Substitutions
69#  include "domzgr_substitute.h90"
70#  include "vectopt_loop_substitute.h90"
71   !!----------------------------------------------------------------------
72   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
73   !! $Id$
74   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   SUBROUTINE eos_insitu( pts, prd )
79      !!----------------------------------------------------------------------
80      !!                   ***  ROUTINE eos_insitu  ***
81      !!
82      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
83      !!       potential temperature and salinity using an equation of state
84      !!       defined through the namelist parameter nn_eos.
85      !!
86      !! ** Method  :   3 cases:
87      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
88      !!         the in situ density is computed directly as a function of
89      !!         potential temperature relative to the surface (the opa t
90      !!         variable), salt and pressure (assuming no pressure variation
91      !!         along geopotential surfaces, i.e. the pressure p in decibars
92      !!         is approximated by the depth in meters.
93      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
94      !!         with pressure                      p        decibars
95      !!              potential temperature         t        deg celsius
96      !!              salinity                      s        psu
97      !!              reference volumic mass        rau0     kg/m**3
98      !!              in situ volumic mass          rho      kg/m**3
99      !!              in situ density anomalie      prd      no units
100      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
101      !!          t = 40 deg celcius, s=40 psu
102      !!      nn_eos = 1 : linear equation of state function of temperature only
103      !!              prd(t) = 0.0285 - rn_alpha * t
104      !!      nn_eos = 2 : linear equation of state function of temperature and
105      !!               salinity
106      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
107      !!      Note that no boundary condition problem occurs in this routine
108      !!      as pts are defined over the whole domain.
109      !!
110      !! ** Action  :   compute prd , the in situ density (no units)
111      !!
112      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
113      !!----------------------------------------------------------------------
114      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
115      USE wrk_nemo, ONLY:   zws => wrk_3d_1   ! 3D workspace
116      !!
117
118!FTRANS zws :I :I :z
119!FTRANS pts :I :I :z :I
120!FTRANS prd :I :I :z
121
122      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
123      !                                                      ! 2 : salinity               [psu]
124      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-]
125      !!
126      INTEGER  ::   ji, jj, jk           ! dummy loop indices
127      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
128      REAL(wp) ::   zr1, zr2, zr3, zr4   !   -      -
129      REAL(wp) ::   zrhop, ze, zbw, zb   !   -      -
130      REAL(wp) ::   zd , zc , zaw, za    !   -      -
131      REAL(wp) ::   zb1, za1, zkw, zk0   !   -      -
132      REAL(wp) ::   zrau0r               !   -      -
133      !!----------------------------------------------------------------------
134
135      IF( wrk_in_use(3, 1) ) THEN
136         CALL ctl_stop('eos_insitu: requested workspace array unavailable')   ;   RETURN
137      ENDIF
138
139      SELECT CASE( nn_eos )
140      !
141      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
142         zrau0r = 1.e0 / rau0
143!CDIR NOVERRCHK
144         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
145         ! 
146#if defined key_z_first
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               DO jk = 1, jpkm1
150#else
151         DO jk = 1, jpkm1
152            DO jj = 1, jpj
153               DO ji = 1, jpi
154#endif
155                  zt = pts   (ji,jj,jk,jp_tem)
156                  zs = pts   (ji,jj,jk,jp_sal)
157                  zh = fsdept(ji,jj,jk)        ! depth
158                  zsr= zws   (ji,jj,jk)        ! square root salinity
159                  !
160                  ! compute volumic mass pure water at atm pressure
161                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
162                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
163                  ! seawater volumic mass atm pressure
164                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
165                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
166                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
167                  zr4= 4.8314e-4_wp
168                  !
169                  ! potential volumic mass (reference to the surface)
170                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
171                  !
172                  ! add the compression terms
173                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
174                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
175                  zb = zbw + ze * zs
176                  !
177                  zd = -2.042967e-2_wp
178                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
179                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
180                  za = ( zd*zsr + zc ) *zs + zaw
181                  !
182                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
183                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
184                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
185                  zk0= ( zb1*zsr + za1 )*zs + zkw
186                  !
187                  ! masked in situ density anomaly
188                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
189                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
190               END DO
191            END DO
192         END DO
193         !
194      CASE( 1 )                !==  Linear formulation function of temperature only  ==!
195#if defined key_z_first
196         DO jj = 1, jpj
197            DO ji = 1, jpi
198               DO jk = 1, jpkm1
199                  prd(ji,jj,jk) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jk,jp_tem) ) * tmask(ji,jj,jk)
200               END DO
201            END DO
202         END DO
203#else
204         DO jk = 1, jpkm1
205            prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
206         END DO
207#endif
208         !
209      CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==!
210#if defined key_z_first
211         DO jj = 1, jpj
212            DO ji = 1, jpi
213               DO jk = 1, jpkm1
214                  prd(ji,jj,jk) = ( rn_beta  * pts(ji,jj,jk,jp_sal) - rn_alpha * pts(ji,jj,jk,jp_tem) ) * tmask(ji,jj,jk)
215               END DO
216            END DO
217         END DO
218#else
219         DO jk = 1, jpkm1
220            prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
221         END DO
222#endif
223         !
224      END SELECT
225      !
226      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk )
227      !
228      IF( wrk_not_released(3, 1) )   CALL ctl_stop('eos_insitu: failed to release workspace array')
229      !
230
231!! * Reset control of array index permutation
232!FTRANS CLEAR
233#  include "dom_oce_ftrans.h90"
234#  include "zdfddm_ftrans.h90"
235
236   END SUBROUTINE eos_insitu
237
238
239   SUBROUTINE eos_insitu_pot( pts, prd, prhop )
240      !!----------------------------------------------------------------------
241      !!                  ***  ROUTINE eos_insitu_pot  ***
242      !!           
243      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
244      !!      potential volumic mass (Kg/m3) from potential temperature and
245      !!      salinity fields using an equation of state defined through the
246      !!     namelist parameter nn_eos.
247      !!
248      !! ** Method  :
249      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
250      !!         the in situ density is computed directly as a function of
251      !!         potential temperature relative to the surface (the opa t
252      !!         variable), salt and pressure (assuming no pressure variation
253      !!         along geopotential surfaces, i.e. the pressure p in decibars
254      !!         is approximated by the depth in meters.
255      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
256      !!              rhop(t,s)  = rho(t,s,0)
257      !!         with pressure                      p        decibars
258      !!              potential temperature         t        deg celsius
259      !!              salinity                      s        psu
260      !!              reference volumic mass        rau0     kg/m**3
261      !!              in situ volumic mass          rho      kg/m**3
262      !!              in situ density anomalie      prd      no units
263      !!
264      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
265      !!          t = 40 deg celcius, s=40 psu
266      !!
267      !!      nn_eos = 1 : linear equation of state function of temperature only
268      !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
269      !!              rhop(t,s)  = rho(t,s)
270      !!
271      !!      nn_eos = 2 : linear equation of state function of temperature and
272      !!               salinity
273      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0
274      !!                       = rn_beta * s - rn_alpha * tn - 1.
275      !!              rhop(t,s)  = rho(t,s)
276      !!      Note that no boundary condition problem occurs in this routine
277      !!      as (tn,sn) or (ta,sa) are defined over the whole domain.
278      !!
279      !! ** Action  : - prd  , the in situ density (no units)
280      !!              - prhop, the potential volumic mass (Kg/m3)
281      !!
282      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
283      !!                Brown and Campana, Mon. Weather Rev., 1978
284      !!----------------------------------------------------------------------
285      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
286      USE wrk_nemo, ONLY:   zws => wrk_3d_1 ! 3D workspace
287      !!
288
289!FTRANS zws :I :I :z
290!FTRANS pts :I :I :z :I
291!FTRANS prd :I :I :z
292!FTRANS prhop :I :I :z
293
294!!DCSE NEMO: This style defeats ftrans
295!     REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
296!     !                                                                ! 2 : salinity               [psu]
297!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
298!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
299      REAL(wp), INTENT(in   ) ::   pts(jpi,jpj,jpkorig,jpts)   ! 1 : potential temperature  [Celcius]
300      !                                                    ! 2 : salinity               [psu]
301      REAL(wp), INTENT(  out) ::   prd(jpi,jpj,jpkorig)        ! in situ density            [-]
302      REAL(wp), INTENT(  out) ::   prhop(jpi,jpj,jpkorig)      ! potential density (surface referenced)
303      !
304      INTEGER  ::   ji, jj, jk   ! dummy loop indices
305      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars
306      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r       !   -      -
307      !!----------------------------------------------------------------------
308
309      IF( wrk_in_use(3, 1) ) THEN
310         CALL ctl_stop('eos_insitu_pot: requested workspace array unavailable')   ;   RETURN
311      ENDIF
312
313      SELECT CASE ( nn_eos )
314      !
315      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
316         zrau0r = 1.e0 / rau0
317!CDIR NOVERRCHK
318         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
319         ! 
320#if defined key_z_first
321         DO jj = 1, jpj
322            DO ji = 1, jpi
323               DO jk = 1, jpkm1
324#else
325         DO jk = 1, jpkm1
326            DO jj = 1, jpj
327               DO ji = 1, jpi
328#endif
329                  zt = pts   (ji,jj,jk,jp_tem)
330                  zs = pts   (ji,jj,jk,jp_sal)
331                  zh = fsdept(ji,jj,jk)        ! depth
332                  zsr= zws   (ji,jj,jk)        ! square root salinity
333                  !
334                  ! compute volumic mass pure water at atm pressure
335                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
336                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
337                  ! seawater volumic mass atm pressure
338                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
339                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp
340                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
341                  zr4= 4.8314e-4_wp
342                  !
343                  ! potential volumic mass (reference to the surface)
344                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
345                  !
346                  ! save potential volumic mass
347                  prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk)
348                  !
349                  ! add the compression terms
350                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
351                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
352                  zb = zbw + ze * zs
353                  !
354                  zd = -2.042967e-2_wp
355                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
356                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
357                  za = ( zd*zsr + zc ) *zs + zaw
358                  !
359                  zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp
360                  za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp
361                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
362                  zk0= ( zb1*zsr + za1 )*zs + zkw
363                  !
364                  ! masked in situ density anomaly
365                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
366                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
367               END DO
368            END DO
369         END DO
370         !
371      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
372#if defined key_z_first
373         DO jj = 1, jpj
374            DO ji = 1, jpi
375               DO jk = 1, jpkm1
376                  prd  (ji,jj,jk) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jk,jp_tem) )        * tmask(ji,jj,jk)
377                  prhop(ji,jj,jk) = ( 1.e0_wp   +            prd(ji,jj,jk)        ) * rau0 * tmask(ji,jj,jk)
378               END DO
379            END DO
380         END DO
381#else
382         DO jk = 1, jpkm1
383            prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
384            prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
385         END DO
386#endif
387         !
388      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
389#if defined key_z_first
390         DO jj = 1, jpj
391            DO ji = 1, jpi
392               DO jk = 1, jpkm1
393                  prd  (ji,jj,jk) = ( rn_beta  * pts(ji,jj,jk,jp_sal) - rn_alpha * pts(ji,jj,jk,jp_tem) )        * tmask(ji,jj,jk)
394                  prhop(ji,jj,jk) = ( 1.e0_wp  + prd(ji,jj,jk)                                          ) * rau0 * tmask(ji,jj,jk)
395               END DO
396            END DO
397         END DO
398#else
399         DO jk = 1, jpkm1
400            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
401            prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
402         END DO
403#endif
404         !
405      END SELECT
406      !
407      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
408      !
409      IF( wrk_not_released(3, 1) )   CALL ctl_stop('eos_insitu_pot: failed to release workspace array')
410      !
411
412!! * Reset control of array index permutation
413!FTRANS CLEAR
414#  include "dom_oce_ftrans.h90"
415#  include "zdfddm_ftrans.h90"
416
417   END SUBROUTINE eos_insitu_pot
418
419
420   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
421      !!----------------------------------------------------------------------
422      !!                  ***  ROUTINE eos_insitu_2d  ***
423      !!
424      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
425      !!      potential temperature and salinity using an equation of state
426      !!      defined through the namelist parameter nn_eos. * 2D field case
427      !!
428      !! ** Method :
429      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
430      !!         the in situ density is computed directly as a function of
431      !!         potential temperature relative to the surface (the opa t
432      !!         variable), salt and pressure (assuming no pressure variation
433      !!         along geopotential surfaces, i.e. the pressure p in decibars
434      !!         is approximated by the depth in meters.
435      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
436      !!         with pressure                      p        decibars
437      !!              potential temperature         t        deg celsius
438      !!              salinity                      s        psu
439      !!              reference volumic mass        rau0     kg/m**3
440      !!              in situ volumic mass          rho      kg/m**3
441      !!              in situ density anomalie      prd      no units
442      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
443      !!          t = 40 deg celcius, s=40 psu
444      !!      nn_eos = 1 : linear equation of state function of temperature only
445      !!              prd(t) = 0.0285 - rn_alpha * t
446      !!      nn_eos = 2 : linear equation of state function of temperature and
447      !!               salinity
448      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
449      !!      Note that no boundary condition problem occurs in this routine
450      !!      as pts are defined over the whole domain.
451      !!
452      !! ** Action  : - prd , the in situ density (no units)
453      !!
454      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
455      !!----------------------------------------------------------------------
456      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
457      USE wrk_nemo, ONLY:   zws => wrk_2d_5 ! 2D workspace
458      !!
459      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
460      !                                                           ! 2 : salinity               [psu]
461      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m]
462      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
463      !!
464      INTEGER  ::   ji, jj                    ! dummy loop indices
465      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary scalars
466      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask        !    -         -
467      !!----------------------------------------------------------------------
468
469      IF( wrk_in_use(2, 5) ) THEN
470         CALL ctl_stop('eos_insitu_2d: requested workspace array unavailable')   ;   RETURN
471      ENDIF
472
473      prd(:,:) = 0._wp
474
475      SELECT CASE( nn_eos )
476      !
477      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
478      !
479!CDIR NOVERRCHK
480         DO jj = 1, jpjm1
481!CDIR NOVERRCHK
482            DO ji = 1, fs_jpim1   ! vector opt.
483               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
484            END DO
485         END DO
486         DO jj = 1, jpjm1
487            DO ji = 1, fs_jpim1   ! vector opt.
488#if defined key_z_first
489               zmask = tmask_1(ji,jj)          ! land/sea bottom mask = surf. mask
490#else
491               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
492#endif
493               zt    = pts  (ji,jj,jp_tem)            ! interpolated T
494               zs    = pts  (ji,jj,jp_sal)            ! interpolated S
495               zsr   = zws  (ji,jj)            ! square root of interpolated S
496               zh    = pdep (ji,jj)            ! depth at the partial step level
497               !
498               ! compute volumic mass pure water at atm pressure
499               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
500                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
501               ! seawater volumic mass atm pressure
502               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   &
503                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp
504               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
505               zr4 = 4.8314e-4_wp
506               !
507               ! potential volumic mass (reference to the surface)
508               zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
509               !
510               ! add the compression terms
511               ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
512               zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
513               zb = zbw + ze * zs
514               !
515               zd =    -2.042967e-2_wp
516               zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
517               zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp
518               za = ( zd*zsr + zc ) *zs + zaw
519               !
520               zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp
521               za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp
522               zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   &
523                  &                             +2098.925_wp ) *zt+190925.6_wp
524               zk0= ( zb1*zsr + za1 )*zs + zkw
525               !
526               ! masked in situ density anomaly
527               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
528            END DO
529         END DO
530         !
531      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
532         DO jj = 1, jpjm1
533            DO ji = 1, fs_jpim1   ! vector opt.
534#if defined key_z_first
535               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask_1(ji,jj)
536#else
537               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
538#endif
539            END DO
540         END DO
541         !
542      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
543         DO jj = 1, jpjm1
544            DO ji = 1, fs_jpim1   ! vector opt.
545#if defined key_z_first
546               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask_1(ji,jj) 
547#else
548               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
549#endif
550            END DO
551         END DO
552         !
553      END SELECT
554
555      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
556      !
557      IF( wrk_not_released(2, 5) )   CALL ctl_stop('eos_insitu_2d: failed to release workspace array')
558      !
559   END SUBROUTINE eos_insitu_2d
560
561
562   SUBROUTINE eos_bn2( pts, pn2 )
563      !!----------------------------------------------------------------------
564      !!                  ***  ROUTINE eos_bn2  ***
565      !!
566      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time-
567      !!      step of the input arguments
568      !!     
569      !! ** Method :
570      !!       * nn_eos = 0  : UNESCO sea water properties
571      !!         The brunt-vaisala frequency is computed using the polynomial
572      !!      polynomial expression of McDougall (1987):
573      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
574      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
575      !!      computed and used in zdfddm module :
576      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
577      !!       * nn_eos = 1  : linear equation of state (temperature only)
578      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
579      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
580      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
581      !!      The use of potential density to compute N^2 introduces e r r o r
582      !!      in the sign of N^2 at great depths. We recommand the use of
583      !!      nn_eos = 0, except for academical studies.
584      !!        Macro-tasked on horizontal slab (jk-loop)
585      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
586      !!      and is never used at this level.
587      !!
588      !! ** Action  : - pn2 : the brunt-vaisala frequency
589      !!
590      !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987.
591      !!----------------------------------------------------------------------
592
593!FTRANS pts :I :I :z :I
594!FTRANS pn2 :I :I :z
595
596!!DCSE_NEMO: This style defeats ftrans
597!     REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
598!     !                                                               ! 2 : salinity               [psu]
599!     REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1]
600
601      REAL(wp), INTENT(in   ) ::   pts(jpi,jpj,jpkorig,jpts)   ! 1 : potential temperature  [Celcius]
602      !                                                    ! 2 : salinity               [psu]
603      REAL(wp), INTENT(  out) ::   pn2(jpi,jpj,jpkorig)        ! Brunt-Vaisala frequency    [s-1]
604      !!
605      INTEGER  ::   ji, jj, jk   ! dummy loop indices
606      REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars
607#if defined key_zdfddm
608      REAL(wp) ::   zds   ! local scalars
609#endif
610      !!----------------------------------------------------------------------
611
612      ! pn2 : interior points only (2=< jk =< jpkm1 )
613      ! --------------------------
614      !
615      SELECT CASE( nn_eos )
616      !
617      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
618#if defined key_z_first
619         DO jj = 1, jpj
620            DO ji = 1, jpi
621               DO jk = 2, jpkm1
622#else
623         DO jk = 2, jpkm1
624            DO jj = 1, jpj
625               DO ji = 1, jpi
626#endif
627                  zgde3w = grav / fse3w(ji,jj,jk)
628                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt
629                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0  ! salinity anomaly (s-35) at w-pt
630                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point
631                  !
632                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta
633                     &                                  - 0.203814e-03_wp ) * zt   &
634                     &                                  + 0.170907e-01_wp ) * zt   &
635                     &   +         0.665157e-01_wp                                 &
636                     &   +     ( - 0.678662e-05_wp * zs                            &
637                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
638                     &   +   ( ( - 0.302285e-13_wp * zh                            &
639                     &           - 0.251520e-11_wp * zs                            &
640                     &           + 0.512857e-12_wp * zt * zt              ) * zh   &
641                     &           - 0.164759e-06_wp * zs                            &
642                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
643                     &                                  + 0.380374e-04_wp ) * zh
644                     !
645                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta
646                     &                               - 0.301985e-05_wp ) * zt      &
647                     &   +       0.785567e-03_wp                                   &
648                     &   + (     0.515032e-08_wp * zs                              &
649                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     &
650                     &   + ( (   0.121551e-17_wp * zh                              &
651                     &         - 0.602281e-15_wp * zs                              &
652                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     &
653                     &                                + 0.408195e-10_wp   * zs     &
654                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     &
655                     &                                - 0.121555e-07_wp ) * zh
656                     !
657                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2
658                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
659                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
660#if defined key_zdfddm
661                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
662                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
663                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
664                  rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
665#endif
666               END DO
667            END DO
668         END DO
669         !
670      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
671#if defined key_z_first
672         DO jj = 1, jpj
673            DO ji = 1, jpi
674               DO jk = 2, jpkm1
675                  pn2(ji,jj,jk) = grav * rn_alpha * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
676                     &                 / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
677               END DO
678            END DO
679         END DO
680#else
681         DO jk = 2, jpkm1
682            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
683         END DO
684#endif
685         !
686      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
687#if defined key_z_first
688         DO jj = 1, jpj
689            DO ji = 1, jpi
690               DO jk = 2, jpkm1
691                  pn2(ji,jj,jk) = grav * (  rn_alpha * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )      &
692                     &                    - rn_beta  * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )   &
693                     &                 / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
694               END DO
695            END DO
696         END DO 
697#else
698         DO jk = 2, jpkm1
699            pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      &
700               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   &
701               &               / fse3w(:,:,jk) * tmask(:,:,jk)
702         END DO 
703#endif
704#if defined key_zdfddm
705#if defined key_z_first
706         DO jj = 1, jpj                                   ! Rrau = (alpha / beta) (dk[t] / dk[s])
707            DO ji = 1, jpi
708               DO jk = 2, jpkm1
709#else
710         DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s])
711            DO jj = 1, jpj
712               DO ji = 1, jpi
713#endif
714                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
715                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
716                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
717               END DO
718            END DO
719         END DO
720#endif
721      END SELECT
722
723      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
724#if defined key_zdfddm
725      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
726#endif
727      !
728
729!! * Reset control of array index permutation
730!FTRANS CLEAR
731#  include "dom_oce_ftrans.h90"
732#  include "zdfddm_ftrans.h90"
733
734   END SUBROUTINE eos_bn2
735
736
737   SUBROUTINE eos_alpbet( pts, palph, pbeta )
738      !!----------------------------------------------------------------------
739      !!                 ***  ROUTINE ldf_slp_grif  ***
740      !!
741      !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points.
742      !!
743      !! ** Method  :   calculates alpha and beta at T-points
744      !!       * nn_eos = 0  : UNESCO sea water properties
745      !!         The brunt-vaisala frequency is computed using the polynomial
746      !!      polynomial expression of McDougall (1987):
747      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
748      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
749      !!      computed and used in zdfddm module :
750      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
751      !!       * nn_eos = 1  : linear equation of state (temperature only)
752      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
753      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
754      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
755      !!       * nn_eos = 3  : Jackett JAOT 2003 ???
756      !!
757      !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point
758      !!----------------------------------------------------------------------
759
760!FTRANS pts :I :I :z :I
761!FTRANS palph :I :I :z
762!FTRANS pbeta :I :I :z
763!!DCSE_NEMO: This style defeats ftrans
764!     REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity
765!     REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff.
766      REAL(wp), INTENT(in   ) ::   pts(jpi,jpj,jpkorig,jpts)            ! pot. temperature & salinity
767      REAL(wp), INTENT(  out) ::   palph(jpi,jpj,jpkorig)               ! thermal expansion coeff.
768      REAL(wp), INTENT(  out) ::   pbeta(jpi,jpj,jpkorig)               ! haline  expansion coeff.
769      !
770      INTEGER  ::   ji, jj, jk   ! dummy loop indices
771      REAL(wp) ::   zt, zs, zh   ! local scalars
772      !!----------------------------------------------------------------------
773      !
774      SELECT CASE ( nn_eos )
775      !
776      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
777#if defined key_z_first
778         DO jj = 1, jpj
779            DO ji = 1, jpi
780               DO jk = 1, jpk
781#else
782         DO jk = 1, jpk
783            DO jj = 1, jpj
784               DO ji = 1, jpi
785#endif
786                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
787                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
788                  zh = fsdept(ji,jj,jk)              ! depth in meters
789                  !
790                  pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   &
791                     &                                        - 0.301985e-05_wp ) * zt   &
792                     &           + 0.785567e-03_wp                                       &
793                     &           + (     0.515032e-08_wp * zs                            &
794                     &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   &
795                     &           + ( (   0.121551e-17_wp * zh                            &
796                     &                 - 0.602281e-15_wp * zs                            &
797                     &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   &
798                     &                                        + 0.408195e-10_wp   * zs   &
799                     &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   &
800                     &                                        - 0.121555e-07_wp ) * zh
801                     !
802                  palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             &
803                      &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &
804                      &                                  - 0.203814e-03_wp ) * zt   &
805                      &                                  + 0.170907e-01_wp ) * zt   &
806                      &   + 0.665157e-01_wp                                         &
807                      &   +     ( - 0.678662e-05_wp * zs                            &
808                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
809                      &   +   ( ( - 0.302285e-13_wp * zh                            &
810                      &           - 0.251520e-11_wp * zs                            &
811                      &           + 0.512857e-12_wp * zt * zt              ) * zh   &
812                      &           - 0.164759e-06_wp * zs                            &
813                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
814                      &                                  + 0.380374e-04_wp ) * zh)
815               END DO
816            END DO
817         END DO
818         !
819      CASE ( 1 )
820         palph(:,:,:) = - rn_alpha
821         pbeta(:,:,:) =   0._wp
822         !
823      CASE ( 2 )
824         palph(:,:,:) = - rn_alpha
825         pbeta(:,:,:) =   rn_beta
826         !
827      CASE DEFAULT
828         IF(lwp) WRITE(numout,cform_err)
829         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
830         nstop = nstop + 1
831         !
832      END SELECT
833      !
834
835!! * Reset control of array index permutation
836!FTRANS CLEAR
837#  include "dom_oce_ftrans.h90"
838#  include "zdfddm_ftrans.h90"
839
840   END SUBROUTINE eos_alpbet
841
842
843   FUNCTION tfreez( psal ) RESULT( ptf )
844      !!----------------------------------------------------------------------
845      !!                 ***  ROUTINE eos_init  ***
846      !!
847      !! ** Purpose :   Compute the sea surface freezing temperature [Celcius]
848      !!
849      !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???)
850      !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
851      !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
852      !!
853      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
854      !!----------------------------------------------------------------------
855      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
856      ! Leave result array automatic rather than making explicitly allocated
857      REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius]
858      !!----------------------------------------------------------------------
859      !
860      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
861         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
862      !
863   END FUNCTION tfreez
864
865
866   SUBROUTINE eos_init
867      !!----------------------------------------------------------------------
868      !!                 ***  ROUTINE eos_init  ***
869      !!
870      !! ** Purpose :   initializations for the equation of state
871      !!
872      !! ** Method  :   Read the namelist nameos and control the parameters
873      !!----------------------------------------------------------------------
874      NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta
875      !!----------------------------------------------------------------------
876      !
877      REWIND( numnam )            ! Read Namelist nameos : equation of state
878      READ  ( numnam, nameos )
879      !
880      IF(lwp) THEN                ! Control print
881         WRITE(numout,*)
882         WRITE(numout,*) 'eos_init : equation of state'
883         WRITE(numout,*) '~~~~~~~~'
884         WRITE(numout,*) '          Namelist nameos : set eos parameters'
885         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
886         WRITE(numout,*) '             thermal exp. coef. (linear)    rn_alpha = ', rn_alpha
887         WRITE(numout,*) '             saline  exp. coef. (linear)    rn_beta  = ', rn_beta
888      ENDIF
889      !
890      SELECT CASE( nn_eos )         ! check option
891      !
892      CASE( 0 )                        !==  Jackett and McDougall (1994) formulation  ==!
893         IF(lwp) WRITE(numout,*)
894         IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1994) equation of state and'
895         IF(lwp) WRITE(numout,*) '                 McDougall (1987) Brunt-Vaisala frequency'
896         !
897      CASE( 1 )                        !==  Linear formulation = F( temperature )  ==!
898         IF(lwp) WRITE(numout,*)
899         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )'
900         IF( lk_zdfddm ) CALL ctl_stop( '          double diffusive mixing parameterization requires',   &
901              &                         ' that T and S are used as state variables' )
902         !
903      CASE( 2 )                        !==  Linear formulation = F( temperature , salinity )  ==!
904         ralpbet = rn_alpha / rn_beta
905         IF(lwp) WRITE(numout,*)
906         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )'
907         !
908      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
909         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
910         CALL ctl_stop( ctmp1 )
911         !
912      END SELECT
913      !
914   END SUBROUTINE eos_init
915
916   !!======================================================================
917END MODULE eosbn2 
Note: See TracBrowser for help on using the repository browser.