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/OPA_SRC – NEMO

source: trunk/NEMO/OPA_SRC/eosbn2.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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