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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OFF_SRC/eosbn2.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

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