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 trunk/NEMO/OFF_SRC – NEMO

source: trunk/NEMO/OFF_SRC/eosbn2.F90 @ 1326

Last change on this file since 1326 was 1326, checked in by cetlod, 15 years ago

suppression of useless CPP key key_mpp_omp, see ticket:350

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