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_crs.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2_crs.F90 @ 7795

Last change on this file since 7795 was 6772, checked in by cbricaud, 8 years ago

clean in coarsening branch
  • Property svn:executable set to *
File size: 49.5 KB
Line 
1MODULE eosbn2_crs
2   !!==============================================================================
3   !!                       ***  MODULE  eosbn2  ***
4   !! Ocean diagnostic variable : equation of state - in situ and potential density
5   !!                                               - Brunt-Vaisala frequency
6   !!==============================================================================
7   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
8   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
9   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
10   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
11   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
12   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
13   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
14   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
15   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
16   !!             -   ! 2003-08  (G. Madec)  F90, free form
17   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function (now eos_fzp function)
18   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add alpha/beta used in ldfslp
20   !!            3.7  ! 2012-03  (F. Roquet, G. Madec)  add primitive of alpha and beta used in PE computation
21   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state
22   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module
23   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80
24   !!----------------------------------------------------------------------
25
26   !!----------------------------------------------------------------------
27   !!   eos            : generic interface of the equation of state
28   !!   eos_insitu     : Compute the in situ density
29   !!   eos_insitu_pot : Compute the insitu and surface referenced potential volumic mass
30   !!   eos_insitu_2d  : Compute the in situ density for 2d fields
31   !!   bn2            : Compute the Brunt-Vaisala frequency
32   !!   eos_rab        : generic interface of in situ thermal/haline expansion ratio
33   !!   eos_rab_3d     : compute in situ thermal/haline expansion ratio
34   !!   eos_rab_2d     : compute in situ thermal/haline expansion ratio for 2d fields
35   !!   eos_fzp_2d     : freezing temperature for 2d fields
36   !!   eos_fzp_0d     : freezing temperature for scalar
37   !!   eos_init       : set eos parameters (namelist)
38   !!----------------------------------------------------------------------
39   USE crs         ! ocean space and time domain
40   USE phycst          ! physical constants
41   !
42   !USE in_out_manager  ! I/O manager
43   !USE lib_mpp         ! MPP library
44   !USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
45   USE prtctl          ! Print control
46   USE wrk_nemo        ! Memory Allocation
47   USE crslbclnk         ! ocean lateral boundary conditions
48   USE timing          ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   !                   !! * Interface
54   INTERFACE eos_crs
55      MODULE PROCEDURE eos_insitu_pot , eos_insitu_2d
56   END INTERFACE
57   !
58   INTERFACE eos_rab_crs
59      MODULE PROCEDURE rab_crs_3d, rab_crs_2d, rab_crs_0d
60   END INTERFACE
61   !
62   PUBLIC   eos_crs        ! called by step, istate, tranpc and zpsgrd modules
63   PUBLIC   eos_rab_crs    ! called by ldfslp, zdfddm, trabbl
64   PUBLIC   eos_init_crs   ! called by istate module
65
66   !                                          !!* Namelist (nameos) *
67   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
68   LOGICAL , PUBLIC ::   ln_useCT  = .FALSE.  ! determine if eos_pt_from_ct is used to compute sst_m
69
70   !                                   !!!  simplified eos coefficients
71   ! default value: Vallis 2006
72   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.
73   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.
74   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2       
75   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2       
76   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T 
77   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S 
78   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt 
79   
80   ! TEOS10/EOS80 parameters
81   REAL(wp) ::   r1_S0, r1_T0, r1_Z0, rdeltaS
82   
83   ! EOS parameters
84   REAL(wp) ::   EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
85   REAL(wp) ::   EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510
86   REAL(wp) ::   EOS020 , EOS120 , EOS220 , EOS320 , EOS420
87   REAL(wp) ::   EOS030 , EOS130 , EOS230 , EOS330
88   REAL(wp) ::   EOS040 , EOS140 , EOS240
89   REAL(wp) ::   EOS050 , EOS150
90   REAL(wp) ::   EOS060
91   REAL(wp) ::   EOS001 , EOS101 , EOS201 , EOS301 , EOS401
92   REAL(wp) ::   EOS011 , EOS111 , EOS211 , EOS311
93   REAL(wp) ::   EOS021 , EOS121 , EOS221
94   REAL(wp) ::   EOS031 , EOS131
95   REAL(wp) ::   EOS041
96   REAL(wp) ::   EOS002 , EOS102 , EOS202
97   REAL(wp) ::   EOS012 , EOS112
98   REAL(wp) ::   EOS022
99   REAL(wp) ::   EOS003 , EOS103
100   REAL(wp) ::   EOS013 
101   
102   ! ALPHA parameters
103   REAL(wp) ::   ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
104   REAL(wp) ::   ALP010 , ALP110 , ALP210 , ALP310 , ALP410
105   REAL(wp) ::   ALP020 , ALP120 , ALP220 , ALP320
106   REAL(wp) ::   ALP030 , ALP130 , ALP230
107   REAL(wp) ::   ALP040 , ALP140
108   REAL(wp) ::   ALP050
109   REAL(wp) ::   ALP001 , ALP101 , ALP201 , ALP301
110   REAL(wp) ::   ALP011 , ALP111 , ALP211
111   REAL(wp) ::   ALP021 , ALP121
112   REAL(wp) ::   ALP031
113   REAL(wp) ::   ALP002 , ALP102
114   REAL(wp) ::   ALP012
115   REAL(wp) ::   ALP003
116   
117   ! BETA parameters
118   REAL(wp) ::   BET000 , BET100 , BET200 , BET300 , BET400 , BET500
119   REAL(wp) ::   BET010 , BET110 , BET210 , BET310 , BET410
120   REAL(wp) ::   BET020 , BET120 , BET220 , BET320
121   REAL(wp) ::   BET030 , BET130 , BET230
122   REAL(wp) ::   BET040 , BET140
123   REAL(wp) ::   BET050
124   REAL(wp) ::   BET001 , BET101 , BET201 , BET301
125   REAL(wp) ::   BET011 , BET111 , BET211
126   REAL(wp) ::   BET021 , BET121
127   REAL(wp) ::   BET031
128   REAL(wp) ::   BET002 , BET102
129   REAL(wp) ::   BET012
130   REAL(wp) ::   BET003
131
132   ! PEN parameters
133   REAL(wp) ::   PEN000 , PEN100 , PEN200 , PEN300 , PEN400
134   REAL(wp) ::   PEN010 , PEN110 , PEN210 , PEN310
135   REAL(wp) ::   PEN020 , PEN120 , PEN220
136   REAL(wp) ::   PEN030 , PEN130
137   REAL(wp) ::   PEN040
138   REAL(wp) ::   PEN001 , PEN101 , PEN201
139   REAL(wp) ::   PEN011 , PEN111
140   REAL(wp) ::   PEN021
141   REAL(wp) ::   PEN002 , PEN102
142   REAL(wp) ::   PEN012
143   
144   ! ALPHA_PEN parameters
145   REAL(wp) ::   APE000 , APE100 , APE200 , APE300
146   REAL(wp) ::   APE010 , APE110 , APE210
147   REAL(wp) ::   APE020 , APE120
148   REAL(wp) ::   APE030
149   REAL(wp) ::   APE001 , APE101
150   REAL(wp) ::   APE011
151   REAL(wp) ::   APE002
152
153   ! BETA_PEN parameters
154   REAL(wp) ::   BPE000 , BPE100 , BPE200 , BPE300
155   REAL(wp) ::   BPE010 , BPE110 , BPE210
156   REAL(wp) ::   BPE020 , BPE120
157   REAL(wp) ::   BPE030
158   REAL(wp) ::   BPE001 , BPE101
159   REAL(wp) ::   BPE011
160   REAL(wp) ::   BPE002
161
162   !! * Substitutions
163#  include "domzgr_substitute.h90"
164#  include "vectopt_loop_substitute.h90"
165   !!----------------------------------------------------------------------
166   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
167   !! $Id: eosbn2.F90 4990 2014-12-15 16:42:49Z timgraham $
168   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
169   !!----------------------------------------------------------------------
170CONTAINS
171
172   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
173      !!----------------------------------------------------------------------
174      !!                  ***  ROUTINE eos_insitu_pot  ***
175      !!
176      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
177      !!      potential volumic mass (Kg/m3) from potential temperature and
178      !!      salinity fields using an equation of state defined through the
179      !!     namelist parameter nn_eos.
180      !!
181      !! ** Action  : - prd  , the in situ density (no units)
182      !!              - prhop, the potential volumic mass (Kg/m3)
183      !!
184      !!----------------------------------------------------------------------
185      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
186      !                                                                ! 2 : salinity               [psu]
187      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
188      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
189      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
190      !
191      INTEGER  ::   ji, jj, jk                ! dummy loop indices
192      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
193      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
194      !!----------------------------------------------------------------------
195      !
196      IF( nn_timing == 1 )   CALL timing_start('eos-pot_crs')
197      !
198      SELECT CASE ( nn_eos )
199      !
200      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
201         !
202         DO jk = 1, jpkm1
203            DO jj = 1, jpj_crs
204               DO ji = 1, jpi_crs
205                  !
206                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
207                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
208                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
209                  ztm = tmask_crs(ji,jj,jk)                                     ! tmask
210                  !
211                  zn3 = EOS013*zt   &
212                     &   + EOS103*zs+EOS003
213                     !
214                  zn2 = (EOS022*zt   &
215                     &   + EOS112*zs+EOS012)*zt   &
216                     &   + (EOS202*zs+EOS102)*zs+EOS002
217                     !
218                  zn1 = (((EOS041*zt   &
219                     &   + EOS131*zs+EOS031)*zt   &
220                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
221                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
222                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
223                     !
224                  zn0 = (((((EOS060*zt   &
225                     &   + EOS150*zs+EOS050)*zt   &
226                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
227                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
228                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
229                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
230                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
231                     !
232                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
233                  !
234                  prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
235                  !
236                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
237               END DO
238            END DO
239         END DO
240         !
241      CASE( 1 )                !==  simplified EOS  ==!
242         !
243         DO jk = 1, jpkm1
244            DO jj = 1, jpj_crs
245               DO ji = 1, jpi_crs
246                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
247                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
248                  zh  = pdep (ji,jj,jk)
249                  ztm = tmask_crs(ji,jj,jk)
250                  !                                                     ! potential density referenced at the surface
251                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
252                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
253                     &  - rn_nu * zt * zs
254                  prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
255                  !                                                     ! density anomaly (masked)
256                  zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
257                  prd(ji,jj,jk) = zn * r1_rau0 * ztm
258                  !
259               END DO
260            END DO
261         END DO
262         !
263      END SELECT
264      !
265      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
266      !
267      IF( nn_timing == 1 )   CALL timing_stop('eos-pot_crs')
268      !
269   END SUBROUTINE eos_insitu_pot
270
271   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
272      !!----------------------------------------------------------------------
273      !!                  ***  ROUTINE eos_insitu_2d  ***
274      !!
275      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
276      !!      potential temperature and salinity using an equation of state
277      !!      defined through the namelist parameter nn_eos. * 2D field case
278      !!
279      !! ** Action  : - prd , the in situ density (no units) (unmasked)
280      !!
281      !!----------------------------------------------------------------------
282      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
283      !                                                           ! 2 : salinity               [psu]
284      REAL(wp), DIMENSION(jpi_crs,jpj_crs)     , INTENT(in   ) ::   pdep  ! depth                      [m]
285      REAL(wp), DIMENSION(jpi_crs,jpj_crs)     , INTENT(  out) ::   prd   ! in situ density
286      !
287      INTEGER  ::   ji, jj, jk                ! dummy loop indices
288      REAL(wp) ::   zt , zh , zs              ! local scalars
289      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
290      !!----------------------------------------------------------------------
291      !
292      IF( nn_timing == 1 )   CALL timing_start('eos2d')
293      !
294      prd(:,:) = 0._wp
295      !
296      SELECT CASE( nn_eos )
297      !
298      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
299         !
300         DO jj = 1, jpj_crsm1
301            DO ji = 1, jpi_crsm1   ! vector opt.
302               !
303               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
304               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
305               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
306               !
307               zn3 = EOS013*zt   &
308                  &   + EOS103*zs+EOS003
309                  !
310               zn2 = (EOS022*zt   &
311                  &   + EOS112*zs+EOS012)*zt   &
312                  &   + (EOS202*zs+EOS102)*zs+EOS002
313                  !
314               zn1 = (((EOS041*zt   &
315                  &   + EOS131*zs+EOS031)*zt   &
316                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
317                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
318                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
319                  !
320               zn0 = (((((EOS060*zt   &
321                  &   + EOS150*zs+EOS050)*zt   &
322                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
323                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
324                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
325                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
326                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
327                  !
328               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
329               !
330               prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
331               !
332            END DO
333         END DO
334         !
335         CALL crs_lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions
336         !
337      CASE( 1 )                !==  simplified EOS  ==!
338         !
339         DO jj = 1, jpj_crsm1
340            DO ji = 1, jpi_crsm1   ! vector opt.
341               !
342               zt    = pts  (ji,jj,jp_tem)  - 10._wp
343               zs    = pts  (ji,jj,jp_sal)  - 35._wp
344               zh    = pdep (ji,jj)                         ! depth at the partial step level
345               !
346               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
347                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
348                  &  - rn_nu * zt * zs
349                  !
350               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
351               !
352            END DO
353         END DO
354         !
355         CALL crs_lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions
356         !
357      END SELECT
358      !
359      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
360      !
361      IF( nn_timing == 1 )   CALL timing_stop('eos2d')
362      !
363   END SUBROUTINE eos_insitu_2d
364
365   SUBROUTINE rab_crs_3d( pts, pab )
366      !!----------------------------------------------------------------------
367      !!                 ***  ROUTINE rab_3d  ***
368      !!
369      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points
370      !!
371      !! ** Method  :   calculates alpha / beta at T-points
372      !!
373      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
374      !!----------------------------------------------------------------------
375      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
376      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
377      !
378      INTEGER  ::   ji, jj, jk                ! dummy loop indices
379      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
380      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
381      !!----------------------------------------------------------------------
382      !
383      IF( nn_timing == 1 )   CALL timing_start('rab_3d')
384      !
385      SELECT CASE ( nn_eos )
386      !
387      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
388         !
389         DO jk = 1, jpkm1
390            DO jj = 1, jpj_crs
391               DO ji = 1, jpi_crs
392                  !
393                  zh  = fsdept_crs(ji,jj,jk) * r1_Z0                                ! depth
394                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
395                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
396                  ztm = tmask_crs(ji,jj,jk)                                         ! tmask
397                  !
398                  ! alpha
399                  zn3 = ALP003
400                  !
401                  zn2 = ALP012*zt + ALP102*zs+ALP002
402                  !
403                  zn1 = ((ALP031*zt   &
404                     &   + ALP121*zs+ALP021)*zt   &
405                     &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
406                     &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
407                     !
408                  zn0 = ((((ALP050*zt   &
409                     &   + ALP140*zs+ALP040)*zt   &
410                     &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
411                     &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
412                     &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
413                     &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
414                     !
415                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
416                  !
417                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
418                  !
419                  ! beta
420                  zn3 = BET003
421                  !
422                  zn2 = BET012*zt + BET102*zs+BET002
423                  !
424                  zn1 = ((BET031*zt   &
425                     &   + BET121*zs+BET021)*zt   &
426                     &   + (BET211*zs+BET111)*zs+BET011)*zt   &
427                     &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
428                     !
429                  zn0 = ((((BET050*zt   &
430                     &   + BET140*zs+BET040)*zt   &
431                     &   + (BET230*zs+BET130)*zs+BET030)*zt   &
432                     &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
433                     &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
434                     &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
435                     !
436                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
437                  !
438                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
439                  !
440               END DO
441            END DO
442         END DO
443         !
444      CASE( 1 )                  !==  simplified EOS  ==!
445         !
446         DO jk = 1, jpkm1
447            DO jj = 1, jpj_crs
448               DO ji = 1, jpi_crs
449                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
450                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
451                  zh  = fsdept_crs(ji,jj,jk)                 ! depth in meters at t-point
452                  ztm = tmask_crs(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
453                  !
454                  zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
455                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
456                  !
457                  zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
458                  pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
459                  !
460               END DO
461            END DO
462         END DO
463         !
464      CASE DEFAULT
465         IF(lwp) WRITE(numout,cform_err)
466         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
467         nstop = nstop + 1
468         !
469      END SELECT
470      !
471      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
472         &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk )
473      !
474      IF( nn_timing == 1 )   CALL timing_stop('rab_3d')
475      !
476   END SUBROUTINE rab_crs_3d
477
478   SUBROUTINE rab_crs_2d( pts, pdep, pab )
479      !!----------------------------------------------------------------------
480      !!                 ***  ROUTINE rab_2d  ***
481      !!
482      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
483      !!
484      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
485      !!----------------------------------------------------------------------
486      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
487      REAL(wp), DIMENSION(jpi_crs,jpj_crs)         , INTENT(in   ) ::   pdep   ! depth                  [m]
488      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
489      !
490      INTEGER  ::   ji, jj, jk                ! dummy loop indices
491      REAL(wp) ::   zt , zh , zs              ! local scalars
492      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
493      !!----------------------------------------------------------------------
494      !
495      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
496      !
497      pab(:,:,:) = 0._wp
498      !
499      SELECT CASE ( nn_eos )
500      !
501      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
502         !
503         DO jj = 1, jpj_crsm1
504            DO ji = 1, jpi_crsm1   ! vector opt.
505               !
506               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
507               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
508               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
509               !
510               ! alpha
511               zn3 = ALP003
512               !
513               zn2 = ALP012*zt + ALP102*zs+ALP002
514               !
515               zn1 = ((ALP031*zt   &
516                  &   + ALP121*zs+ALP021)*zt   &
517                  &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
518                  &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
519                  !
520               zn0 = ((((ALP050*zt   &
521                  &   + ALP140*zs+ALP040)*zt   &
522                  &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
523                  &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
524                  &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
525                  &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
526                  !
527               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
528               !
529               pab(ji,jj,jp_tem) = zn * r1_rau0
530               !
531               ! beta
532               zn3 = BET003
533               !
534               zn2 = BET012*zt + BET102*zs+BET002
535               !
536               zn1 = ((BET031*zt   &
537                  &   + BET121*zs+BET021)*zt   &
538                  &   + (BET211*zs+BET111)*zs+BET011)*zt   &
539                  &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
540                  !
541               zn0 = ((((BET050*zt   &
542                  &   + BET140*zs+BET040)*zt   &
543                  &   + (BET230*zs+BET130)*zs+BET030)*zt   &
544                  &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
545                  &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
546                  &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
547                  !
548               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
549               !
550               pab(ji,jj,jp_sal) = zn / zs * r1_rau0
551               !
552               !
553            END DO
554         END DO
555         !
556         CALL crs_lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions
557         CALL crs_lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                   
558         !
559      CASE( 1 )                  !==  simplified EOS  ==!
560         !
561         DO jj = 1, jpj_crsm1
562            DO ji = 1, jpi_crsm1   ! vector opt.
563               !
564               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
565               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
566               zh    = pdep (ji,jj)                   ! depth at the partial step level
567               !
568               zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
569               pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
570               !
571               zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
572               pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
573               !
574            END DO
575         END DO
576         !
577         CALL crs_lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions
578         CALL crs_lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                   
579         !
580      CASE DEFAULT
581         IF(lwp) WRITE(numout,cform_err)
582         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
583         nstop = nstop + 1
584         !
585      END SELECT
586      !
587      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
588         &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
589      !
590      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
591      !
592   END SUBROUTINE rab_crs_2d
593
594
595   SUBROUTINE rab_crs_0d( pts, pdep, pab )
596      !!----------------------------------------------------------------------
597      !!                 ***  ROUTINE rab_0d  ***
598      !!
599      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
600      !!
601      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
602      !!----------------------------------------------------------------------
603      REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
604      REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m]
605      REAL(wp), DIMENSION(jpts)    , INTENT(  out) ::   pab    ! thermal/haline expansion ratio
606      !
607      REAL(wp) ::   zt , zh , zs              ! local scalars
608      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
609      !!----------------------------------------------------------------------
610      !
611      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
612      !
613      pab(:) = 0._wp
614      !
615      SELECT CASE ( nn_eos )
616      !
617      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
618         !
619         !
620         zh  = pdep * r1_Z0                                  ! depth
621         zt  = pts (jp_tem) * r1_T0                           ! temperature
622         zs  = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
623         !
624         ! alpha
625         zn3 = ALP003
626         !
627         zn2 = ALP012*zt + ALP102*zs+ALP002
628         !
629         zn1 = ((ALP031*zt   &
630            &   + ALP121*zs+ALP021)*zt   &
631            &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
632            &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
633            !
634         zn0 = ((((ALP050*zt   &
635            &   + ALP140*zs+ALP040)*zt   &
636            &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
637            &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
638            &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
639            &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
640            !
641         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
642         !
643         pab(jp_tem) = zn * r1_rau0
644         !
645         ! beta
646         zn3 = BET003
647         !
648         zn2 = BET012*zt + BET102*zs+BET002
649         !
650         zn1 = ((BET031*zt   &
651            &   + BET121*zs+BET021)*zt   &
652            &   + (BET211*zs+BET111)*zs+BET011)*zt   &
653            &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
654            !
655         zn0 = ((((BET050*zt   &
656            &   + BET140*zs+BET040)*zt   &
657            &   + (BET230*zs+BET130)*zs+BET030)*zt   &
658            &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
659            &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
660            &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
661            !
662         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
663         !
664         pab(jp_sal) = zn / zs * r1_rau0
665         !
666         !
667         !
668      CASE( 1 )                  !==  simplified EOS  ==!
669         !
670         zt    = pts(jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
671         zs    = pts(jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
672         zh    = pdep                    ! depth at the partial step level
673         !
674         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
675         pab(jp_tem) = zn * r1_rau0   ! alpha
676         !
677         zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
678         pab(jp_sal) = zn * r1_rau0   ! beta
679         !
680      CASE DEFAULT
681         IF(lwp) WRITE(numout,cform_err)
682         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
683         nstop = nstop + 1
684         !
685      END SELECT
686      !
687      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
688      !
689   END SUBROUTINE rab_crs_0d
690
691   SUBROUTINE eos_init_crs
692      !!----------------------------------------------------------------------
693      !!                 ***  ROUTINE eos_init  ***
694      !!
695      !! ** Purpose :   initializations for the equation of state
696      !!
697      !! ** Method  :   Read the namelist nameos and control the parameters
698      !!----------------------------------------------------------------------
699      INTEGER  ::   ios   ! local integer
700      !!
701      NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
702         &                                             rn_lambda2, rn_mu2, rn_nu
703      !!----------------------------------------------------------------------
704      !
705      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
706      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
707901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
708      !
709      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
710      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
711902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
712      IF(lwm) WRITE( numond, nameos )
713      !
714      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
715      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
716      !
717      IF(lwp) THEN                ! Control print
718         WRITE(numout,*)
719         WRITE(numout,*) 'eos_init : equation of state'
720         WRITE(numout,*) '~~~~~~~~'
721         WRITE(numout,*) '          Namelist nameos : set eos parameters'
722         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
723         IF( ln_useCT )   THEN
724            WRITE(numout,*) '             model uses Conservative Temperature'
725            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
726         ENDIF
727      ENDIF
728      !
729      SELECT CASE( nn_eos )         ! check option
730      !
731      CASE( -1 )                       !==  polynomial TEOS-10  ==!
732         IF(lwp) WRITE(numout,*)
733         IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
734         !
735         rdeltaS = 32._wp
736         r1_S0  = 0.875_wp/35.16504_wp
737         r1_T0  = 1._wp/40._wp
738         r1_Z0  = 1.e-4_wp
739         !
740         EOS000 = 8.0189615746e+02_wp
741         EOS100 = 8.6672408165e+02_wp
742         EOS200 = -1.7864682637e+03_wp
743         EOS300 = 2.0375295546e+03_wp
744         EOS400 = -1.2849161071e+03_wp
745         EOS500 = 4.3227585684e+02_wp
746         EOS600 = -6.0579916612e+01_wp
747         EOS010 = 2.6010145068e+01_wp
748         EOS110 = -6.5281885265e+01_wp
749         EOS210 = 8.1770425108e+01_wp
750         EOS310 = -5.6888046321e+01_wp
751         EOS410 = 1.7681814114e+01_wp
752         EOS510 = -1.9193502195_wp
753         EOS020 = -3.7074170417e+01_wp
754         EOS120 = 6.1548258127e+01_wp
755         EOS220 = -6.0362551501e+01_wp
756         EOS320 = 2.9130021253e+01_wp
757         EOS420 = -5.4723692739_wp
758         EOS030 = 2.1661789529e+01_wp
759         EOS130 = -3.3449108469e+01_wp
760         EOS230 = 1.9717078466e+01_wp
761         EOS330 = -3.1742946532_wp
762         EOS040 = -8.3627885467_wp
763         EOS140 = 1.1311538584e+01_wp
764         EOS240 = -5.3563304045_wp
765         EOS050 = 5.4048723791e-01_wp
766         EOS150 = 4.8169980163e-01_wp
767         EOS060 = -1.9083568888e-01_wp
768         EOS001 = 1.9681925209e+01_wp
769         EOS101 = -4.2549998214e+01_wp
770         EOS201 = 5.0774768218e+01_wp
771         EOS301 = -3.0938076334e+01_wp
772         EOS401 = 6.6051753097_wp
773         EOS011 = -1.3336301113e+01_wp
774         EOS111 = -4.4870114575_wp
775         EOS211 = 5.0042598061_wp
776         EOS311 = -6.5399043664e-01_wp
777         EOS021 = 6.7080479603_wp
778         EOS121 = 3.5063081279_wp
779         EOS221 = -1.8795372996_wp
780         EOS031 = -2.4649669534_wp
781         EOS131 = -5.5077101279e-01_wp
782         EOS041 = 5.5927935970e-01_wp
783         EOS002 = 2.0660924175_wp
784         EOS102 = -4.9527603989_wp
785         EOS202 = 2.5019633244_wp
786         EOS012 = 2.0564311499_wp
787         EOS112 = -2.1311365518e-01_wp
788         EOS022 = -1.2419983026_wp
789         EOS003 = -2.3342758797e-02_wp
790         EOS103 = -1.8507636718e-02_wp
791         EOS013 = 3.7969820455e-01_wp
792         !
793         ALP000 = -6.5025362670e-01_wp
794         ALP100 = 1.6320471316_wp
795         ALP200 = -2.0442606277_wp
796         ALP300 = 1.4222011580_wp
797         ALP400 = -4.4204535284e-01_wp
798         ALP500 = 4.7983755487e-02_wp
799         ALP010 = 1.8537085209_wp
800         ALP110 = -3.0774129064_wp
801         ALP210 = 3.0181275751_wp
802         ALP310 = -1.4565010626_wp
803         ALP410 = 2.7361846370e-01_wp
804         ALP020 = -1.6246342147_wp
805         ALP120 = 2.5086831352_wp
806         ALP220 = -1.4787808849_wp
807         ALP320 = 2.3807209899e-01_wp
808         ALP030 = 8.3627885467e-01_wp
809         ALP130 = -1.1311538584_wp
810         ALP230 = 5.3563304045e-01_wp
811         ALP040 = -6.7560904739e-02_wp
812         ALP140 = -6.0212475204e-02_wp
813         ALP050 = 2.8625353333e-02_wp
814         ALP001 = 3.3340752782e-01_wp
815         ALP101 = 1.1217528644e-01_wp
816         ALP201 = -1.2510649515e-01_wp
817         ALP301 = 1.6349760916e-02_wp
818         ALP011 = -3.3540239802e-01_wp
819         ALP111 = -1.7531540640e-01_wp
820         ALP211 = 9.3976864981e-02_wp
821         ALP021 = 1.8487252150e-01_wp
822         ALP121 = 4.1307825959e-02_wp
823         ALP031 = -5.5927935970e-02_wp
824         ALP002 = -5.1410778748e-02_wp
825         ALP102 = 5.3278413794e-03_wp
826         ALP012 = 6.2099915132e-02_wp
827         ALP003 = -9.4924551138e-03_wp
828         !
829         BET000 = 1.0783203594e+01_wp
830         BET100 = -4.4452095908e+01_wp
831         BET200 = 7.6048755820e+01_wp
832         BET300 = -6.3944280668e+01_wp
833         BET400 = 2.6890441098e+01_wp
834         BET500 = -4.5221697773_wp
835         BET010 = -8.1219372432e-01_wp
836         BET110 = 2.0346663041_wp
837         BET210 = -2.1232895170_wp
838         BET310 = 8.7994140485e-01_wp
839         BET410 = -1.1939638360e-01_wp
840         BET020 = 7.6574242289e-01_wp
841         BET120 = -1.5019813020_wp
842         BET220 = 1.0872489522_wp
843         BET320 = -2.7233429080e-01_wp
844         BET030 = -4.1615152308e-01_wp
845         BET130 = 4.9061350869e-01_wp
846         BET230 = -1.1847737788e-01_wp
847         BET040 = 1.4073062708e-01_wp
848         BET140 = -1.3327978879e-01_wp
849         BET050 = 5.9929880134e-03_wp
850         BET001 = -5.2937873009e-01_wp
851         BET101 = 1.2634116779_wp
852         BET201 = -1.1547328025_wp
853         BET301 = 3.2870876279e-01_wp
854         BET011 = -5.5824407214e-02_wp
855         BET111 = 1.2451933313e-01_wp
856         BET211 = -2.4409539932e-02_wp
857         BET021 = 4.3623149752e-02_wp
858         BET121 = -4.6767901790e-02_wp
859         BET031 = -6.8523260060e-03_wp
860         BET002 = -6.1618945251e-02_wp
861         BET102 = 6.2255521644e-02_wp
862         BET012 = -2.6514181169e-03_wp
863         BET003 = -2.3025968587e-04_wp
864         !
865         PEN000 = -9.8409626043_wp
866         PEN100 = 2.1274999107e+01_wp
867         PEN200 = -2.5387384109e+01_wp
868         PEN300 = 1.5469038167e+01_wp
869         PEN400 = -3.3025876549_wp
870         PEN010 = 6.6681505563_wp
871         PEN110 = 2.2435057288_wp
872         PEN210 = -2.5021299030_wp
873         PEN310 = 3.2699521832e-01_wp
874         PEN020 = -3.3540239802_wp
875         PEN120 = -1.7531540640_wp
876         PEN220 = 9.3976864981e-01_wp
877         PEN030 = 1.2324834767_wp
878         PEN130 = 2.7538550639e-01_wp
879         PEN040 = -2.7963967985e-01_wp
880         PEN001 = -1.3773949450_wp
881         PEN101 = 3.3018402659_wp
882         PEN201 = -1.6679755496_wp
883         PEN011 = -1.3709540999_wp
884         PEN111 = 1.4207577012e-01_wp
885         PEN021 = 8.2799886843e-01_wp
886         PEN002 = 1.7507069098e-02_wp
887         PEN102 = 1.3880727538e-02_wp
888         PEN012 = -2.8477365341e-01_wp
889         !
890         APE000 = -1.6670376391e-01_wp
891         APE100 = -5.6087643219e-02_wp
892         APE200 = 6.2553247576e-02_wp
893         APE300 = -8.1748804580e-03_wp
894         APE010 = 1.6770119901e-01_wp
895         APE110 = 8.7657703198e-02_wp
896         APE210 = -4.6988432490e-02_wp
897         APE020 = -9.2436260751e-02_wp
898         APE120 = -2.0653912979e-02_wp
899         APE030 = 2.7963967985e-02_wp
900         APE001 = 3.4273852498e-02_wp
901         APE101 = -3.5518942529e-03_wp
902         APE011 = -4.1399943421e-02_wp
903         APE002 = 7.1193413354e-03_wp
904         !
905         BPE000 = 2.6468936504e-01_wp
906         BPE100 = -6.3170583896e-01_wp
907         BPE200 = 5.7736640125e-01_wp
908         BPE300 = -1.6435438140e-01_wp
909         BPE010 = 2.7912203607e-02_wp
910         BPE110 = -6.2259666565e-02_wp
911         BPE210 = 1.2204769966e-02_wp
912         BPE020 = -2.1811574876e-02_wp
913         BPE120 = 2.3383950895e-02_wp
914         BPE030 = 3.4261630030e-03_wp
915         BPE001 = 4.1079296834e-02_wp
916         BPE101 = -4.1503681096e-02_wp
917         BPE011 = 1.7676120780e-03_wp
918         BPE002 = 1.7269476440e-04_wp
919         !
920      CASE( 0 )                        !==  polynomial EOS-80 formulation  ==!
921         !
922         IF(lwp) WRITE(numout,*)
923         IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)'
924         !
925         rdeltaS = 20._wp
926         r1_S0  = 1._wp/40._wp
927         r1_T0  = 1._wp/40._wp
928         r1_Z0  = 1.e-4_wp
929         !
930         EOS000 = 9.5356891948e+02_wp
931         EOS100 = 1.7136499189e+02_wp
932         EOS200 = -3.7501039454e+02_wp
933         EOS300 = 5.1856810420e+02_wp
934         EOS400 = -3.7264470465e+02_wp
935         EOS500 = 1.4302533998e+02_wp
936         EOS600 = -2.2856621162e+01_wp
937         EOS010 = 1.0087518651e+01_wp
938         EOS110 = -1.3647741861e+01_wp
939         EOS210 = 8.8478359933_wp
940         EOS310 = -7.2329388377_wp
941         EOS410 = 1.4774410611_wp
942         EOS510 = 2.0036720553e-01_wp
943         EOS020 = -2.5579830599e+01_wp
944         EOS120 = 2.4043512327e+01_wp
945         EOS220 = -1.6807503990e+01_wp
946         EOS320 = 8.3811577084_wp
947         EOS420 = -1.9771060192_wp
948         EOS030 = 1.6846451198e+01_wp
949         EOS130 = -2.1482926901e+01_wp
950         EOS230 = 1.0108954054e+01_wp
951         EOS330 = -6.2675951440e-01_wp
952         EOS040 = -8.0812310102_wp
953         EOS140 = 1.0102374985e+01_wp
954         EOS240 = -4.8340368631_wp
955         EOS050 = 1.2079167803_wp
956         EOS150 = 1.1515380987e-01_wp
957         EOS060 = -2.4520288837e-01_wp
958         EOS001 = 1.0748601068e+01_wp
959         EOS101 = -1.7817043500e+01_wp
960         EOS201 = 2.2181366768e+01_wp
961         EOS301 = -1.6750916338e+01_wp
962         EOS401 = 4.1202230403_wp
963         EOS011 = -1.5852644587e+01_wp
964         EOS111 = -7.6639383522e-01_wp
965         EOS211 = 4.1144627302_wp
966         EOS311 = -6.6955877448e-01_wp
967         EOS021 = 9.9994861860_wp
968         EOS121 = -1.9467067787e-01_wp
969         EOS221 = -1.2177554330_wp
970         EOS031 = -3.4866102017_wp
971         EOS131 = 2.2229155620e-01_wp
972         EOS041 = 5.9503008642e-01_wp
973         EOS002 = 1.0375676547_wp
974         EOS102 = -3.4249470629_wp
975         EOS202 = 2.0542026429_wp
976         EOS012 = 2.1836324814_wp
977         EOS112 = -3.4453674320e-01_wp
978         EOS022 = -1.2548163097_wp
979         EOS003 = 1.8729078427e-02_wp
980         EOS103 = -5.7238495240e-02_wp
981         EOS013 = 3.8306136687e-01_wp
982         !
983         ALP000 = -2.5218796628e-01_wp
984         ALP100 = 3.4119354654e-01_wp
985         ALP200 = -2.2119589983e-01_wp
986         ALP300 = 1.8082347094e-01_wp
987         ALP400 = -3.6936026529e-02_wp
988         ALP500 = -5.0091801383e-03_wp
989         ALP010 = 1.2789915300_wp
990         ALP110 = -1.2021756164_wp
991         ALP210 = 8.4037519952e-01_wp
992         ALP310 = -4.1905788542e-01_wp
993         ALP410 = 9.8855300959e-02_wp
994         ALP020 = -1.2634838399_wp
995         ALP120 = 1.6112195176_wp
996         ALP220 = -7.5817155402e-01_wp
997         ALP320 = 4.7006963580e-02_wp
998         ALP030 = 8.0812310102e-01_wp
999         ALP130 = -1.0102374985_wp
1000         ALP230 = 4.8340368631e-01_wp
1001         ALP040 = -1.5098959754e-01_wp
1002         ALP140 = -1.4394226233e-02_wp
1003         ALP050 = 3.6780433255e-02_wp
1004         ALP001 = 3.9631611467e-01_wp
1005         ALP101 = 1.9159845880e-02_wp
1006         ALP201 = -1.0286156825e-01_wp
1007         ALP301 = 1.6738969362e-02_wp
1008         ALP011 = -4.9997430930e-01_wp
1009         ALP111 = 9.7335338937e-03_wp
1010         ALP211 = 6.0887771651e-02_wp
1011         ALP021 = 2.6149576513e-01_wp
1012         ALP121 = -1.6671866715e-02_wp
1013         ALP031 = -5.9503008642e-02_wp
1014         ALP002 = -5.4590812035e-02_wp
1015         ALP102 = 8.6134185799e-03_wp
1016         ALP012 = 6.2740815484e-02_wp
1017         ALP003 = -9.5765341718e-03_wp
1018         !
1019         BET000 = 2.1420623987_wp
1020         BET100 = -9.3752598635_wp
1021         BET200 = 1.9446303907e+01_wp
1022         BET300 = -1.8632235232e+01_wp
1023         BET400 = 8.9390837485_wp
1024         BET500 = -1.7142465871_wp
1025         BET010 = -1.7059677327e-01_wp
1026         BET110 = 2.2119589983e-01_wp
1027         BET210 = -2.7123520642e-01_wp
1028         BET310 = 7.3872053057e-02_wp
1029         BET410 = 1.2522950346e-02_wp
1030         BET020 = 3.0054390409e-01_wp
1031         BET120 = -4.2018759976e-01_wp
1032         BET220 = 3.1429341406e-01_wp
1033         BET320 = -9.8855300959e-02_wp
1034         BET030 = -2.6853658626e-01_wp
1035         BET130 = 2.5272385134e-01_wp
1036         BET230 = -2.3503481790e-02_wp
1037         BET040 = 1.2627968731e-01_wp
1038         BET140 = -1.2085092158e-01_wp
1039         BET050 = 1.4394226233e-03_wp
1040         BET001 = -2.2271304375e-01_wp
1041         BET101 = 5.5453416919e-01_wp
1042         BET201 = -6.2815936268e-01_wp
1043         BET301 = 2.0601115202e-01_wp
1044         BET011 = -9.5799229402e-03_wp
1045         BET111 = 1.0286156825e-01_wp
1046         BET211 = -2.5108454043e-02_wp
1047         BET021 = -2.4333834734e-03_wp
1048         BET121 = -3.0443885826e-02_wp
1049         BET031 = 2.7786444526e-03_wp
1050         BET002 = -4.2811838287e-02_wp
1051         BET102 = 5.1355066072e-02_wp
1052         BET012 = -4.3067092900e-03_wp
1053         BET003 = -7.1548119050e-04_wp
1054         !
1055         PEN000 = -5.3743005340_wp
1056         PEN100 = 8.9085217499_wp
1057         PEN200 = -1.1090683384e+01_wp
1058         PEN300 = 8.3754581690_wp
1059         PEN400 = -2.0601115202_wp
1060         PEN010 = 7.9263222935_wp
1061         PEN110 = 3.8319691761e-01_wp
1062         PEN210 = -2.0572313651_wp
1063         PEN310 = 3.3477938724e-01_wp
1064         PEN020 = -4.9997430930_wp
1065         PEN120 = 9.7335338937e-02_wp
1066         PEN220 = 6.0887771651e-01_wp
1067         PEN030 = 1.7433051009_wp
1068         PEN130 = -1.1114577810e-01_wp
1069         PEN040 = -2.9751504321e-01_wp
1070         PEN001 = -6.9171176978e-01_wp
1071         PEN101 = 2.2832980419_wp
1072         PEN201 = -1.3694684286_wp
1073         PEN011 = -1.4557549876_wp
1074         PEN111 = 2.2969116213e-01_wp
1075         PEN021 = 8.3654420645e-01_wp
1076         PEN002 = -1.4046808820e-02_wp
1077         PEN102 = 4.2928871430e-02_wp
1078         PEN012 = -2.8729602515e-01_wp
1079         !
1080         APE000 = -1.9815805734e-01_wp
1081         APE100 = -9.5799229402e-03_wp
1082         APE200 = 5.1430784127e-02_wp
1083         APE300 = -8.3694846809e-03_wp
1084         APE010 = 2.4998715465e-01_wp
1085         APE110 = -4.8667669469e-03_wp
1086         APE210 = -3.0443885826e-02_wp
1087         APE020 = -1.3074788257e-01_wp
1088         APE120 = 8.3359333577e-03_wp
1089         APE030 = 2.9751504321e-02_wp
1090         APE001 = 3.6393874690e-02_wp
1091         APE101 = -5.7422790533e-03_wp
1092         APE011 = -4.1827210323e-02_wp
1093         APE002 = 7.1824006288e-03_wp
1094         !
1095         BPE000 = 1.1135652187e-01_wp
1096         BPE100 = -2.7726708459e-01_wp
1097         BPE200 = 3.1407968134e-01_wp
1098         BPE300 = -1.0300557601e-01_wp
1099         BPE010 = 4.7899614701e-03_wp
1100         BPE110 = -5.1430784127e-02_wp
1101         BPE210 = 1.2554227021e-02_wp
1102         BPE020 = 1.2166917367e-03_wp
1103         BPE120 = 1.5221942913e-02_wp
1104         BPE030 = -1.3893222263e-03_wp
1105         BPE001 = 2.8541225524e-02_wp
1106         BPE101 = -3.4236710714e-02_wp
1107         BPE011 = 2.8711395266e-03_wp
1108         BPE002 = 5.3661089288e-04_wp
1109         !
1110      CASE( 1 )                        !==  Simplified EOS     ==!
1111         IF(lwp) THEN
1112            WRITE(numout,*)
1113            WRITE(numout,*) '          use of simplified eos:    rhd(dT=T-10,dS=S-35,Z) = '
1114            WRITE(numout,*) '             [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0'
1115            WRITE(numout,*)
1116            WRITE(numout,*) '             thermal exp. coef.    rn_a0      = ', rn_a0
1117            WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0
1118            WRITE(numout,*) '             cabbeling coef.       rn_lambda1 = ', rn_lambda1
1119            WRITE(numout,*) '             cabbeling coef.       rn_lambda2 = ', rn_lambda2
1120            WRITE(numout,*) '             thermobar. coef.      rn_mu1     = ', rn_mu1
1121            WRITE(numout,*) '             thermobar. coef.      rn_mu2     = ', rn_mu2
1122            WRITE(numout,*) '             2nd cabbel. coef.     rn_nu      = ', rn_nu
1123            WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization '
1124         ENDIF
1125         !
1126      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
1127         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
1128         CALL ctl_stop( ctmp1 )
1129         !
1130      END SELECT
1131      !
1132      r1_rau0     = 1._wp / rau0
1133      r1_rcp      = 1._wp / rcp
1134      r1_rau0_rcp = 1._wp / ( rau0 * rcp )
1135      !
1136      IF(lwp) WRITE(numout,*)
1137      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
1138      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
1139      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
1140      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
1141      !
1142   END SUBROUTINE eos_init_crs
1143
1144   !!======================================================================
1145END MODULE eosbn2_crs
Note: See TracBrowser for help on using the repository browser.