17 use gsw_mod_toolbox,
only : gsw_rho_first_derivatives
19 implicit none ;
private
23 public calculate_density_scalar_nemo, calculate_density_array_nemo
29 module procedure calculate_density_scalar_nemo, calculate_density_array_nemo
35 module procedure calculate_density_derivs_scalar_nemo, calculate_density_derivs_array_nemo
38 real,
parameter :: pa2db = 1.e-4
40 real,
parameter :: rdeltas = 32.
41 real,
parameter :: r1_s0 = 0.875/35.16504
42 real,
parameter :: r1_t0 = 1./40.
43 real,
parameter :: r1_p0 = 1.e-4
44 real,
parameter :: r00 = 4.6494977072e+01
45 real,
parameter :: r01 = -5.2099962525
46 real,
parameter :: r02 = 2.2601900708e-01
47 real,
parameter :: r03 = 6.4326772569e-02
48 real,
parameter :: r04 = 1.5616995503e-02
49 real,
parameter :: r05 = -1.7243708991e-03
50 real,
parameter :: eos000 = 8.0189615746e+02
51 real,
parameter :: eos100 = 8.6672408165e+02
52 real,
parameter :: eos200 = -1.7864682637e+03
53 real,
parameter :: eos300 = 2.0375295546e+03
54 real,
parameter :: eos400 = -1.2849161071e+03
55 real,
parameter :: eos500 = 4.3227585684e+02
56 real,
parameter :: eos600 = -6.0579916612e+01
57 real,
parameter :: eos010 = 2.6010145068e+01
58 real,
parameter :: eos110 = -6.5281885265e+01
59 real,
parameter :: eos210 = 8.1770425108e+01
60 real,
parameter :: eos310 = -5.6888046321e+01
61 real,
parameter :: eos410 = 1.7681814114e+01
62 real,
parameter :: eos510 = -1.9193502195
63 real,
parameter :: eos020 = -3.7074170417e+01
64 real,
parameter :: eos120 = 6.1548258127e+01
65 real,
parameter :: eos220 = -6.0362551501e+01
66 real,
parameter :: eos320 = 2.9130021253e+01
67 real,
parameter :: eos420 = -5.4723692739
68 real,
parameter :: eos030 = 2.1661789529e+01
69 real,
parameter :: eos130 = -3.3449108469e+01
70 real,
parameter :: eos230 = 1.9717078466e+01
71 real,
parameter :: eos330 = -3.1742946532
72 real,
parameter :: eos040 = -8.3627885467
73 real,
parameter :: eos140 = 1.1311538584e+01
74 real,
parameter :: eos240 = -5.3563304045
75 real,
parameter :: eos050 = 5.4048723791e-01
76 real,
parameter :: eos150 = 4.8169980163e-01
77 real,
parameter :: eos060 = -1.9083568888e-01
78 real,
parameter :: eos001 = 1.9681925209e+01
79 real,
parameter :: eos101 = -4.2549998214e+01
80 real,
parameter :: eos201 = 5.0774768218e+01
81 real,
parameter :: eos301 = -3.0938076334e+01
82 real,
parameter :: eos401 = 6.6051753097
83 real,
parameter :: eos011 = -1.3336301113e+01
84 real,
parameter :: eos111 = -4.4870114575
85 real,
parameter :: eos211 = 5.0042598061
86 real,
parameter :: eos311 = -6.5399043664e-01
87 real,
parameter :: eos021 = 6.7080479603
88 real,
parameter :: eos121 = 3.5063081279
89 real,
parameter :: eos221 = -1.8795372996
90 real,
parameter :: eos031 = -2.4649669534
91 real,
parameter :: eos131 = -5.5077101279e-01
92 real,
parameter :: eos041 = 5.5927935970e-01
93 real,
parameter :: eos002 = 2.0660924175
94 real,
parameter :: eos102 = -4.9527603989
95 real,
parameter :: eos202 = 2.5019633244
96 real,
parameter :: eos012 = 2.0564311499
97 real,
parameter :: eos112 = -2.1311365518e-01
98 real,
parameter :: eos022 = -1.2419983026
99 real,
parameter :: eos003 = -2.3342758797e-02
100 real,
parameter :: eos103 = -1.8507636718e-02
101 real,
parameter :: eos013 = 3.7969820455e-01
102 real,
parameter :: alp000 = -6.5025362670e-01
103 real,
parameter :: alp100 = 1.6320471316
104 real,
parameter :: alp200 = -2.0442606277
105 real,
parameter :: alp300 = 1.4222011580
106 real,
parameter :: alp400 = -4.4204535284e-01
107 real,
parameter :: alp500 = 4.7983755487e-02
108 real,
parameter :: alp010 = 1.8537085209
109 real,
parameter :: alp110 = -3.0774129064
110 real,
parameter :: alp210 = 3.0181275751
111 real,
parameter :: alp310 = -1.4565010626
112 real,
parameter :: alp410 = 2.7361846370e-01
113 real,
parameter :: alp020 = -1.6246342147
114 real,
parameter :: alp120 = 2.5086831352
115 real,
parameter :: alp220 = -1.4787808849
116 real,
parameter :: alp320 = 2.3807209899e-01
117 real,
parameter :: alp030 = 8.3627885467e-01
118 real,
parameter :: alp130 = -1.1311538584
119 real,
parameter :: alp230 = 5.3563304045e-01
120 real,
parameter :: alp040 = -6.7560904739e-02
121 real,
parameter :: alp140 = -6.0212475204e-02
122 real,
parameter :: alp050 = 2.8625353333e-02
123 real,
parameter :: alp001 = 3.3340752782e-01
124 real,
parameter :: alp101 = 1.1217528644e-01
125 real,
parameter :: alp201 = -1.2510649515e-01
126 real,
parameter :: alp301 = 1.6349760916e-02
127 real,
parameter :: alp011 = -3.3540239802e-01
128 real,
parameter :: alp111 = -1.7531540640e-01
129 real,
parameter :: alp211 = 9.3976864981e-02
130 real,
parameter :: alp021 = 1.8487252150e-01
131 real,
parameter :: alp121 = 4.1307825959e-02
132 real,
parameter :: alp031 = -5.5927935970e-02
133 real,
parameter :: alp002 = -5.1410778748e-02
134 real,
parameter :: alp102 = 5.3278413794e-03
135 real,
parameter :: alp012 = 6.2099915132e-02
136 real,
parameter :: alp003 = -9.4924551138e-03
137 real,
parameter :: bet000 = 1.0783203594e+01
138 real,
parameter :: bet100 = -4.4452095908e+01
139 real,
parameter :: bet200 = 7.6048755820e+01
140 real,
parameter :: bet300 = -6.3944280668e+01
141 real,
parameter :: bet400 = 2.6890441098e+01
142 real,
parameter :: bet500 = -4.5221697773
143 real,
parameter :: bet010 = -8.1219372432e-01
144 real,
parameter :: bet110 = 2.0346663041
145 real,
parameter :: bet210 = -2.1232895170
146 real,
parameter :: bet310 = 8.7994140485e-01
147 real,
parameter :: bet410 = -1.1939638360e-01
148 real,
parameter :: bet020 = 7.6574242289e-01
149 real,
parameter :: bet120 = -1.5019813020
150 real,
parameter :: bet220 = 1.0872489522
151 real,
parameter :: bet320 = -2.7233429080e-01
152 real,
parameter :: bet030 = -4.1615152308e-01
153 real,
parameter :: bet130 = 4.9061350869e-01
154 real,
parameter :: bet230 = -1.1847737788e-01
155 real,
parameter :: bet040 = 1.4073062708e-01
156 real,
parameter :: bet140 = -1.3327978879e-01
157 real,
parameter :: bet050 = 5.9929880134e-03
158 real,
parameter :: bet001 = -5.2937873009e-01
159 real,
parameter :: bet101 = 1.2634116779
160 real,
parameter :: bet201 = -1.1547328025
161 real,
parameter :: bet301 = 3.2870876279e-01
162 real,
parameter :: bet011 = -5.5824407214e-02
163 real,
parameter :: bet111 = 1.2451933313e-01
164 real,
parameter :: bet211 = -2.4409539932e-02
165 real,
parameter :: bet021 = 4.3623149752e-02
166 real,
parameter :: bet121 = -4.6767901790e-02
167 real,
parameter :: bet031 = -6.8523260060e-03
168 real,
parameter :: bet002 = -6.1618945251e-02
169 real,
parameter :: bet102 = 6.2255521644e-02
170 real,
parameter :: bet012 = -2.6514181169e-03
171 real,
parameter :: bet003 = -2.3025968587e-04
180 subroutine calculate_density_scalar_nemo(T, S, pressure, rho, rho_ref)
181 real,
intent(in) :: t
182 real,
intent(in) :: s
183 real,
intent(in) :: pressure
184 real,
intent(out) :: rho
185 real,
optional,
intent(in) :: rho_ref
187 real :: al0, p0, lambda
189 real,
dimension(1) :: t0, s0, pressure0
190 real,
dimension(1) :: rho0
194 pressure0(1) = pressure
196 call calculate_density_array_nemo(t0, s0, pressure0, rho0, 1, 1, rho_ref)
199 end subroutine calculate_density_scalar_nemo
205 subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts, rho_ref)
206 real,
dimension(:),
intent(in) :: t
207 real,
dimension(:),
intent(in) :: s
208 real,
dimension(:),
intent(in) :: pressure
209 real,
dimension(:),
intent(out) :: rho
210 integer,
intent(in) :: start
211 integer,
intent(in) :: npts
212 real,
optional,
intent(in) :: rho_ref
215 real :: zp, zt, zh, zs, zr0, zn, zn0, zn1, zn2, zn3, zs0
218 do j=start,start+npts-1
222 zp = pressure(j)* pa2db
228 zs = sqrt( abs( zs + rdeltas ) * r1_s0 )
234 & + eos112*zs+eos012)*zt &
235 & + (eos202*zs+eos102)*zs+eos002
238 & + eos131*zs+eos031)*zt &
239 & + (eos221*zs+eos121)*zs+eos021)*zt &
240 & + ((eos311*zs+eos211)*zs+eos111)*zs+eos011)*zt &
241 & + (((eos401*zs+eos301)*zs+eos201)*zs+eos101)*zs+eos001
243 zn0 = (((((eos060*zt &
244 & + eos150*zs+eos050)*zt &
245 & + (eos240*zs+eos140)*zs+eos040)*zt &
246 & + ((eos330*zs+eos230)*zs+eos130)*zs+eos030)*zt &
247 & + (((eos420*zs+eos320)*zs+eos220)*zs+eos120)*zs+eos020)*zt &
248 & + ((((eos510*zs+eos410)*zs+eos310)*zs+eos210)*zs+eos110)*zs+eos010)*zt
250 zs0 = (((((eos600*zs+eos500)*zs+eos400)*zs+eos300)*zs+eos200)*zs+eos100)*zs + eos000
252 zr0 = (((((r05 * zp+r04) * zp+r03 ) * zp+r02 ) * zp+r01) * zp+r00) * zp
254 if (
present(rho_ref))
then
255 zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + (zn0 + (zs0 - rho_ref))
256 rho(j) = ( zn + zr0 )
258 zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + (zn0 + zs0)
259 rho(j) = ( zn + zr0 )
263 end subroutine calculate_density_array_nemo
267 subroutine calculate_density_derivs_array_nemo(T, S, pressure, drho_dT, drho_dS, start, npts)
268 real,
intent(in),
dimension(:) :: T
269 real,
intent(in),
dimension(:) :: S
270 real,
intent(in),
dimension(:) :: pressure
271 real,
intent(out),
dimension(:) :: drho_dT
273 real,
intent(out),
dimension(:) :: drho_dS
275 integer,
intent(in) :: start
276 integer,
intent(in) :: npts
279 real :: zp,zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3
282 do j=start,start+npts-1
286 zp = pressure(j)* pa2db
292 zs = sqrt( abs( zs + rdeltas ) * r1_s0 )
297 zn2 = alp012*zt + alp102*zs+alp002
300 & + alp121*zs+alp021)*zt &
301 & + (alp211*zs+alp111)*zs+alp011)*zt &
302 & + ((alp301*zs+alp201)*zs+alp101)*zs+alp001
304 zn0 = ((((alp050*zt &
305 & + alp140*zs+alp040)*zt &
306 & + (alp230*zs+alp130)*zs+alp030)*zt &
307 & + ((alp320*zs+alp220)*zs+alp120)*zs+alp020)*zt &
308 & + (((alp410*zs+alp310)*zs+alp210)*zs+alp110)*zs+alp010)*zt &
309 & + ((((alp500*zs+alp400)*zs+alp300)*zs+alp200)*zs+alp100)*zs+alp000
311 zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0
319 zn2 = bet012*zt + bet102*zs+bet002
322 & + bet121*zs+bet021)*zt &
323 & + (bet211*zs+bet111)*zs+bet011)*zt &
324 & + ((bet301*zs+bet201)*zs+bet101)*zs+bet001
326 zn0 = ((((bet050*zt &
327 & + bet140*zs+bet040)*zt &
328 & + (bet230*zs+bet130)*zs+bet030)*zt &
329 & + ((bet320*zs+bet220)*zs+bet120)*zs+bet020)*zt &
330 & + (((bet410*zs+bet310)*zs+bet210)*zs+bet110)*zs+bet010)*zt &
331 & + ((((bet500*zs+bet400)*zs+bet300)*zs+bet200)*zs+bet100)*zs+bet000
333 zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0
338 end subroutine calculate_density_derivs_array_nemo
341 subroutine calculate_density_derivs_scalar_nemo(T, S, pressure, drho_dt, drho_ds)
342 real,
intent(in) :: T
343 real,
intent(in) :: S
344 real,
intent(in) :: pressure
345 real,
intent(out) :: drho_dT
347 real,
intent(out) :: drho_dS
350 real :: al0, p0, lambda
352 real,
dimension(1) :: T0, S0, pressure0
353 real,
dimension(1) :: drdt0, drds0
357 pressure0(1) = pressure
359 call calculate_density_derivs_array_nemo(t0, s0, pressure0, drdt0, drds0, 1, 1)
362 end subroutine calculate_density_derivs_scalar_nemo
368 subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts)
369 real,
intent(in),
dimension(:) :: t
370 real,
intent(in),
dimension(:) :: s
371 real,
intent(in),
dimension(:) :: pressure
372 real,
intent(out),
dimension(:) :: rho
373 real,
intent(out),
dimension(:) :: drho_dp
376 integer,
intent(in) :: start
377 integer,
intent(in) :: npts
383 call calculate_density_array_nemo(t, s, pressure, rho, start, npts)
388 do j=start,start+npts-1
392 zp = pressure(j)* pa2db
393 call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j))
395 end subroutine calculate_compress_nemo