FV3 Bundle
gsw_frazil_properties_potential.f90
Go to the documentation of this file.
1 !==========================================================================
2 elemental subroutine gsw_frazil_properties_potential (sa_bulk, h_pot_bulk,&
3  p, sa_final, ct_final, w_ih_final)
4 !==========================================================================
5 !
6 ! Calculates the mass fraction of ice (mass of ice divided by mass of ice
7 ! plus seawater), w_Ih_final, which results from given values of the bulk
8 ! Absolute Salinity, SA_bulk, bulk potential enthalpy, h_pot_bulk,
9 ! occuring at pressure p. The final equilibrium values of Absolute
10 ! Salinity, SA_final, and Conservative Temperature, CT_final, of the
11 ! interstitial seawater phase are also returned. This code assumes that
12 ! there is no dissolved air in the seawater (that is, saturation_fraction
13 ! is assumed to be zero thoughout the code).
14 !
15 ! When the mass fraction w_Ih_final is calculated as being a positive
16 ! value, the seawater-ice mixture is at thermodynamic equlibrium.
17 !
18 ! This code returns w_Ih_final = 0 when the input bulk enthalpy, h_bulk,
19 ! is sufficiently large (i.e. sufficiently "warm") so that there is no ice
20 ! present in the final state. In this case the final state consists of
21 ! only seawater rather than being an equlibrium mixture of seawater and
22 ! ice which occurs when w_Ih_final is positive. Note that when
23 ! w_Ih_final = 0, the final seawater is not at the freezing temperature.
24 !
25 ! Note that this code uses the exact forms of CT_freezing and
26 ! pot_enthalpy_ice_freezing.
27 !
28 ! SA_bulk = bulk Absolute Salinity of the seawater and ice mixture
29 ! [ g/kg ]
30 ! h_pot_bulk = bulk potential enthalpy of the seawater and ice mixture
31 ! [ J/kg ]
32 ! p = sea pressure [ dbar ]
33 ! ( i.e. absolute pressure - 10.1325 dbar )
34 !
35 ! SA_final = Absolute Salinity of the seawater in the final state,
36 ! whether or not any ice is present. [ g/kg ]
37 ! CT_final = Conservative Temperature of the seawater in the the final
38 ! state, whether or not any ice is present. [ deg C ]
39 ! w_Ih_final = mass fraction of ice in the final seawater-ice mixture.
40 ! If this ice mass fraction is positive, the system is at
41 ! thermodynamic equilibrium. If this ice mass fraction is
42 ! zero there is no ice in the final state which consists
43 ! only of seawater which is warmer than the freezing
44 ! temperature. [unitless]
45 !--------------------------------------------------------------------------
46 
48 
53 
55 
56 use gsw_mod_kinds
57 
58 implicit none
59 
60 real (r8), intent(in) :: sa_bulk, h_pot_bulk, p
61 real (r8), intent(out) :: sa_final, ct_final, w_ih_final
62 
63 integer :: iterations, max_iterations
64 
65 real (r8) :: ctf_sa, ctf, dfunc_dw_ih, dfunc_dw_ih_mean_poly, dpot_h_ihf_dsa
66 real (r8) :: func, func0, h_pot_ihf, sa, w_ih_old, w_ih, x, xa, y, z
67 
68 real (r8), parameter :: f01 = -9.041191886754806e-1_r8
69 real (r8), parameter :: f02 = 4.169608567309818e-2_r8
70 real (r8), parameter :: f03 = -9.325971761333677e-3_r8
71 real (r8), parameter :: f04 = 4.699055851002199e-2_r8
72 real (r8), parameter :: f05 = -3.086923404061666e-2_r8
73 real (r8), parameter :: f06 = 1.057761186019000e-2_r8
74 real (r8), parameter :: f07 = -7.349302346007727e-2_r8
75 real (r8), parameter :: f08 = 1.444842576424337e-1_r8
76 real (r8), parameter :: f09 = -1.408425967872030e-1_r8
77 real (r8), parameter :: f10 = 1.070981398567760e-1_r8
78 real (r8), parameter :: f11 = -1.768451760854797e-2_r8
79 real (r8), parameter :: f12 = -4.013688314067293e-1_r8
80 real (r8), parameter :: f13 = 7.209753205388577e-1_r8
81 real (r8), parameter :: f14 = -1.807444462285120e-1_r8
82 real (r8), parameter :: f15 = 1.362305015808993e-1_r8
83 real (r8), parameter :: f16 = -9.500974920072897e-1_r8
84 real (r8), parameter :: f17 = 1.192134856624248_r8
85 real (r8), parameter :: f18 = -9.191161283559850e-2_r8
86 real (r8), parameter :: f19 = -1.008594411490973_r8
87 real (r8), parameter :: f20 = 8.020279271484482e-1_r8
88 real (r8), parameter :: f21 = -3.930534388853466e-1_r8
89 real (r8), parameter :: f22 = -2.026853316399942e-2_r8
90 real (r8), parameter :: f23 = -2.722731069001690e-2_r8
91 real (r8), parameter :: f24 = 5.032098120548072e-2_r8
92 real (r8), parameter :: f25 = -2.354888890484222e-2_r8
93 real (r8), parameter :: f26 = -2.454090179215001e-2_r8
94 real (r8), parameter :: f27 = 4.125987229048937e-2_r8
95 real (r8), parameter :: f28 = -3.533404753585094e-2_r8
96 real (r8), parameter :: f29 = 3.766063025852511e-2_r8
97 real (r8), parameter :: f30 = -3.358409746243470e-2_r8
98 real (r8), parameter :: f31 = -2.242158862056258e-2_r8
99 real (r8), parameter :: f32 = 2.102254738058931e-2_r8
100 real (r8), parameter :: f33 = -3.048635435546108e-2_r8
101 real (r8), parameter :: f34 = -1.996293091714222e-2_r8
102 real (r8), parameter :: f35 = 2.577703068234217e-2_r8
103 real (r8), parameter :: f36 = -1.292053030649309e-2_r8
104 
105 real (r8), parameter :: g01 = 3.332286683867741e5_r8
106 real (r8), parameter :: g02 = 1.416532517833479e4_r8
107 real (r8), parameter :: g03 = -1.021129089258645e4_r8
108 real (r8), parameter :: g04 = 2.356370992641009e4_r8
109 real (r8), parameter :: g05 = -8.483432350173174e3_r8
110 real (r8), parameter :: g06 = 2.279927781684362e4_r8
111 real (r8), parameter :: g07 = 1.506238790315354e4_r8
112 real (r8), parameter :: g08 = 4.194030718568807e3_r8
113 real (r8), parameter :: g09 = -3.146939594885272e5_r8
114 real (r8), parameter :: g10 = -7.549939721380912e4_r8
115 real (r8), parameter :: g11 = 2.790535212869292e6_r8
116 real (r8), parameter :: g12 = 1.078851928118102e5_r8
117 real (r8), parameter :: g13 = -1.062493860205067e7_r8
118 real (r8), parameter :: g14 = 2.082909703458225e7_r8
119 real (r8), parameter :: g15 = -2.046810820868635e7_r8
120 real (r8), parameter :: g16 = 8.039606992745191e6_r8
121 real (r8), parameter :: g17 = -2.023984705844567e4_r8
122 real (r8), parameter :: g18 = 2.871769638352535e4_r8
123 real (r8), parameter :: g19 = -1.444841553038544e4_r8
124 real (r8), parameter :: g20 = 2.261532522236573e4_r8
125 real (r8), parameter :: g21 = -2.090579366221046e4_r8
126 real (r8), parameter :: g22 = -1.128417003723530e4_r8
127 real (r8), parameter :: g23 = 3.222965226084112e3_r8
128 real (r8), parameter :: g24 = -1.226388046175992e4_r8
129 real (r8), parameter :: g25 = 1.506847628109789e4_r8
130 real (r8), parameter :: g26 = -4.584670946447444e4_r8
131 real (r8), parameter :: g27 = 1.596119496322347e4_r8
132 real (r8), parameter :: g28 = -6.338852410446789e4_r8
133 real (r8), parameter :: g29 = 8.951570926106525e4_r8
134 
135 real (r8), parameter :: saturation_fraction = 0.0_r8
136 
137 character (*), parameter :: func_name = "gsw_frazil_properties_potential"
138 
139 !--------------------------------------------------------------------------
140 ! Finding func0. This is the value of the function, func, that would
141 ! result in the output w_Ih_final being exactly zero.
142 !--------------------------------------------------------------------------
143 func0 = h_pot_bulk - &
144  gsw_cp0*gsw_ct_freezing_exact(sa_bulk,p,saturation_fraction)
145 
146 !--------------------------------------------------------------------------
147 ! Setting the three outputs for data points that have func0 non-negative
148 !--------------------------------------------------------------------------
149 if (func0 >= 0.0_r8) then
150  ! When func0 is zero or positive then the final answer will contain
151  ! no frazil ice; that is, it will be pure seawater that is warmer
152  ! han the freezing temperature. If func0 >= 0 we do not need to go
153  ! through the modified Newton-Raphson procedure and we can simply
154  ! write down the answer, as in the following 4 lines of code.
155  sa_final = sa_bulk
156  ct_final = h_pot_bulk/gsw_cp0
157  w_ih_final = 0.0_r8
158  return
159 end if
160 
161 !--------------------------------------------------------------------------
162 ! Begin finding the solution for data points that have func0 < 0, so that
163 ! the output will have a positive ice mass fraction w_Ih_final.
164 !--------------------------------------------------------------------------
165 
166 ! Evalaute a polynomial for w_Ih in terms of SA_bulk, func0 and p
167 x = sa_bulk*1e-2_r8
168 y = func0/3e5_r8
169 z = p*1e-4_r8
170 
171 w_ih = y*(f01 + x*(f02 + x*(f03 + x*(f04 + x*(f05 + f06*x)))) &
172  + y*(f07 + x*(f08 + x*(f09 + x*(f10 + f11*x))) + y*(f12 + x*(f13 &
173  + x*(f14 + f15*x)) + y*(f16 + x*(f17 + f18*x) + y*(f19 + f20*x &
174  + f21*y)))) + z*(f22 + x*(f23 + x*(f24 + f25*x)) + y*(x*(f26 + f27*x) &
175  + y*(f28 + f29*x + f30*y)) + z*(f31 + x*(f32 + f33*x) + y*(f34 &
176  + f35*x + f36*y))))
177 
178 if (w_ih .gt. 0.9_r8) then
179  ! The ice mass fraction out of this code is restricted to be less than 0.9.
180  sa_final = gsw_error_code(1,func_name)
181  ct_final = sa_final
182  w_ih_final = sa_final
183  return
184 end if
185 
186 ! The initial guess at the absolute salinity of the interstitial seawater
187 sa = sa_bulk/(1.0_r8 - w_ih)
188 
189 !--------------------------------------------------------------------------
190 ! Doing a Newton step with a separate polynomial estimate of the mean
191 ! derivative dfunc_dw_Ih_mean_poly.
192 !--------------------------------------------------------------------------
193 ctf = gsw_ct_freezing_exact(sa,p,saturation_fraction)
194 h_pot_ihf = gsw_pot_enthalpy_ice_freezing(sa,p)
195 func = h_pot_bulk - (1.0_r8 - w_ih)*gsw_cp0*ctf - w_ih*h_pot_ihf
196 
197 xa = sa*1e-2_r8
198 
199 dfunc_dw_ih_mean_poly = g01 + xa*(g02 + xa*(g03 + xa*(g04 + g05*xa))) &
200  + w_ih*(xa*(g06 + xa*(g07 + g08*xa)) + w_ih*(xa*(g09 + g10*xa) &
201  + w_ih*xa*(g11 + g12*xa + w_ih*(g13 + w_ih*(g14 + w_ih*(g15 &
202  + g16*w_ih)))))) + z*(g17 + xa*(g18 + g19*xa) + w_ih*(g20 &
203  + w_ih*(g21 + g22*w_ih) + xa*(g23 + g24*xa*w_ih)) &
204  + z*(g25 + xa*(g26 + g27*xa) + w_ih*(g28 + g29*w_ih)))
205 
206 w_ih_old = w_ih
207 w_ih = w_ih_old - func/dfunc_dw_ih_mean_poly
208 
209 sa = sa_bulk/(1.0_r8 - w_ih)
210 
211 !--------------------------------------------------------------------------
212 ! Calculating the estimate of the derivative of func, dfunc_dw_Ih, to be
213 ! fed into Newton's Method.
214 !--------------------------------------------------------------------------
215 ctf = gsw_ct_freezing_exact(sa,p,saturation_fraction)
216 
217 h_pot_ihf = gsw_pot_enthalpy_ice_freezing(sa,p)
218 
219 call gsw_ct_freezing_first_derivatives(sa,p,saturation_fraction,ctf_sa)
220 call gsw_pot_enthalpy_ice_freezing_first_derivatives(sa,p,dpot_h_ihf_dsa)
221 
222 dfunc_dw_ih = gsw_cp0*ctf - h_pot_ihf - &
223  sa*(gsw_cp0*ctf_sa + w_ih*dpot_h_ihf_dsa/(1.0_r8 - w_ih))
224 
225 if (w_ih .ge. 0.0_r8 .and. w_ih .le. 0.20_r8 .and. sa .gt. 15.0_r8 &
226  .and. sa .lt. 60.0_r8 .and. p .le. 3000.0_r8) then
227  max_iterations = 1
228 else if (w_ih .ge. 0.0_r8 .and. w_ih .le. 0.85_r8 .and. sa .gt. 0.0_r8 &
229  .and. sa .lt. 120.0_r8 .and. p .le. 3500.0_r8) then
230  max_iterations = 2
231 else
232  max_iterations = 3
233 end if
234 
235 do iterations = 1, max_iterations
236 
237  if (iterations .gt. 1) then
238  ! On the first iteration ctf and h_pot_ihf are both known
239  ctf = gsw_ct_freezing_exact(sa,p,saturation_fraction)
240  h_pot_ihf = gsw_pot_enthalpy_ice_freezing(sa,p)
241  end if
242 
243  ! This is the function, func, whose zero we seek ...
244  func = h_pot_bulk - (1.0_r8 - w_ih)*gsw_cp0*ctf - w_ih*h_pot_ihf
245 
246  w_ih_old = w_ih
247  w_ih = w_ih_old - func/dfunc_dw_ih
248 
249  if (w_ih .gt. 0.9_r8) then
250  sa_final = gsw_error_code(1+iterations,func_name)
251  ct_final = sa_final
252  w_ih_final = sa_final
253  return
254  end if
255 
256  sa = sa_bulk/(1.0_r8 - w_ih)
257 
258 end do
259 
260 if (w_ih .lt. 0.0_r8) then
261  sa_final = sa_bulk
262  ct_final = h_pot_bulk/gsw_cp0
263  w_ih_final = 0.0_r8
264 else
265  sa_final = sa
266  ct_final = gsw_ct_freezing_exact(sa,p,saturation_fraction)
267  w_ih_final = w_ih
268 end if
269 
270 return
271 end subroutine
272 
273 !--------------------------------------------------------------------------
real(r8), parameter, public gsw_error_limit
elemental subroutine gsw_frazil_properties_potential(sa_bulk, h_pot_bulk, p, sa_final, ct_final, w_ih_final)
elemental real(r8) function, public gsw_error_code(err_num, func_name, error_code)