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 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • 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,jpk,jpts)   ! 1 : potential temperature  [Celcius]
300      !                                                    ! 2 : salinity               [psu]
301      REAL(wp), INTENT(  out) ::   prd(jpi,jpj,jpk)        ! in situ density            [-]
302      REAL(wp), INTENT(  out) ::   prhop(jpi,jpj,jpk)      ! 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,jpk,jpts)   ! 1 : potential temperature  [Celcius]
602      !                                                    ! 2 : salinity               [psu]
603      REAL(wp), INTENT(  out) ::   pn2(jpi,jpj,jpk)        ! 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,jpk,jpts)            ! pot. temperature & salinity
767      REAL(wp), INTENT(  out) ::   palph(jpi,jpj,jpk)               ! thermal expansion coeff.
768      REAL(wp), INTENT(  out) ::   pbeta(jpi,jpj,jpk)               ! 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.