Line data Source code
1 : !=========================================================================
2 : !
3 : ! Update ice and snow internal temperatures
4 : ! using zero-layer thermodynamics
5 : !
6 : ! authors: Alison McLaren, UK MetOffice
7 : ! Elizabeth C. Hunke, LANL
8 : !
9 : ! 2012: Split from icepack_therm_vertical.F90
10 :
11 : module icepack_therm_0layer
12 :
13 : use icepack_kinds
14 : use icepack_parameters, only: c0, c1, p5, puny
15 : use icepack_parameters, only: kseaice, ksno
16 : use icepack_therm_bl99, only: surface_fluxes
17 : use icepack_warnings, only: warnstr, icepack_warnings_add
18 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
19 :
20 : implicit none
21 :
22 : private
23 : public :: zerolayer_temperature
24 :
25 : !=======================================================================
26 :
27 : contains
28 :
29 : !=======================================================================
30 : !
31 : ! Compute new surface temperature using zero layer model of Semtner
32 : ! (1976).
33 : !
34 : ! New temperatures are computed iteratively by solving a
35 : ! surface flux balance equation (i.e. net surface flux from atmos
36 : ! equals conductive flux from the top to the bottom surface).
37 : !
38 : ! author: Alison McLaren, Met Office
39 : ! (but largely taken from temperature_changes)
40 :
41 108633465 : subroutine zerolayer_temperature(nilyr, nslyr, &
42 : rhoa, flw, &
43 : potT, Qa, &
44 : shcoef, lhcoef, &
45 : fswsfc, &
46 : hilyr, hslyr, &
47 : Tsf, Tbot, &
48 : fsensn, flatn, &
49 : flwoutn, fsurfn, &
50 : fcondtopn,fcondbot )
51 :
52 : integer (kind=int_kind), intent(in) :: &
53 : nilyr , & ! number of ice layers
54 : nslyr ! number of snow layers
55 :
56 : real (kind=dbl_kind), intent(in) :: &
57 : rhoa , & ! air density (kg/m^3)
58 : flw , & ! incoming longwave radiation (W/m^2)
59 : potT , & ! air potential temperature (K)
60 : Qa , & ! specific humidity (kg/kg)
61 : shcoef , & ! transfer coefficient for sensible heat
62 : lhcoef , & ! transfer coefficient for latent heat
63 : Tbot , & ! ice bottom surface temperature (deg C)
64 : fswsfc ! SW absorbed at ice/snow surface (W m-2)
65 :
66 : real (kind=dbl_kind), intent(in) :: &
67 : hilyr , & ! ice layer thickness (m)
68 : hslyr ! snow layer thickness (m)
69 :
70 : real (kind=dbl_kind), intent(inout):: &
71 : fsensn , & ! surface downward sensible heat (W m-2)
72 : flatn , & ! surface downward latent heat (W m-2)
73 : flwoutn , & ! upward LW at surface (W m-2)
74 : fsurfn , & ! net flux to top surface, excluding fcondtopn
75 : fcondtopn ! downward cond flux at top surface (W m-2)
76 :
77 : real (kind=dbl_kind), intent(out):: &
78 : fcondbot ! downward cond flux at bottom surface (W m-2)
79 :
80 : real (kind=dbl_kind), &
81 : intent(inout) :: &
82 : Tsf ! ice/snow surface temperature, Tsfcn
83 :
84 : ! local variables
85 :
86 : logical (kind=log_kind), parameter :: &
87 : l_zerolayerchecks = .true.
88 :
89 : integer (kind=int_kind), parameter :: &
90 : nitermax = 50 ! max number of iterations in temperature solver
91 :
92 : real (kind=dbl_kind), parameter :: &
93 : Tsf_errmax = 5.e-4_dbl_kind ! max allowed error in Tsf
94 : ! recommend Tsf_errmax < 0.01 K
95 :
96 : integer (kind=int_kind) :: &
97 : niter ! iteration counter in temperature solver
98 :
99 : real (kind=dbl_kind) :: &
100 17906208 : Tsf_start , & ! Tsf at start of iteration
101 17906208 : dTsf , & ! Tsf - Tsf_start
102 17906208 : dfsurf_dT ! derivative of fsurfn wrt Tsf
103 :
104 : real (kind=dbl_kind) :: &
105 17906208 : dTsf_prev , & ! dTsf from previous iteration
106 17906208 : dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1)
107 17906208 : dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1)
108 17906208 : dflwout_dT ! deriv of flwout wrt Tsf (W m-2 deg-1)
109 :
110 : real (kind=dbl_kind) :: &
111 17906208 : kh , & ! effective conductivity
112 17906208 : diag , & ! diagonal matrix elements
113 17906208 : rhs ! rhs of tri-diagonal matrix equation
114 :
115 : real (kind=dbl_kind) :: &
116 17906208 : heff , & ! effective ice thickness (m)
117 : ! ( hice + hsno*kseaice/ksnow)
118 17906208 : kratio , & ! ratio of ice and snow conductivies
119 17906208 : avg_Tsf ! = 1. if Tsf averaged w/Tsf_start, else = 0.
120 :
121 : logical (kind=log_kind) :: &
122 : converged ! = true when local solution has converged
123 :
124 : character(len=*),parameter :: subname='(zerolayer_temperature)'
125 :
126 : !-----------------------------------------------------------------
127 : ! Initialize
128 : !-----------------------------------------------------------------
129 :
130 108633465 : fcondbot = c0
131 :
132 108633465 : converged = .false.
133 :
134 108633465 : dTsf_prev = c0
135 :
136 : !-----------------------------------------------------------------
137 : ! Solve for new temperatures.
138 : ! Iterate until temperatures converge with minimal temperature
139 : ! change.
140 : !-----------------------------------------------------------------
141 :
142 5540306715 : do niter = 1, nitermax
143 :
144 5540306715 : if (.not. converged) then
145 :
146 : !-----------------------------------------------------------------
147 : ! Update radiative and turbulent fluxes and their derivatives
148 : ! with respect to Tsf.
149 : !-----------------------------------------------------------------
150 :
151 : call surface_fluxes (Tsf, fswsfc, &
152 : rhoa, flw, &
153 : potT, Qa, &
154 : shcoef, lhcoef, &
155 : flwoutn, fsensn, &
156 : flatn, fsurfn, &
157 : dflwout_dT, dfsens_dT, &
158 205765486 : dflat_dT, dfsurf_dT)
159 205765486 : if (icepack_warnings_aborted(subname)) return
160 :
161 : !-----------------------------------------------------------------
162 : ! Compute effective ice thickness (includes snow) and thermal
163 : ! conductivity
164 : !-----------------------------------------------------------------
165 :
166 205765486 : kratio = kseaice/ksno
167 :
168 205765486 : heff = hilyr + kratio * hslyr
169 205765486 : kh = kseaice / heff
170 :
171 : !-----------------------------------------------------------------
172 : ! Compute conductive flux at top surface, fcondtopn.
173 : ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less
174 : ! than zero (but not less than -puny).
175 : !-----------------------------------------------------------------
176 :
177 205765486 : fcondtopn = kh * (Tsf - Tbot)
178 :
179 205765486 : if (fsurfn < fcondtopn) &
180 194092600 : Tsf = min (Tsf, -puny)
181 :
182 : !-----------------------------------------------------------------
183 : ! Save surface temperature at start of iteration
184 : !-----------------------------------------------------------------
185 :
186 205765486 : Tsf_start = Tsf
187 :
188 : !-----------------------------------------------------------------
189 : ! Solve surface balance equation to obtain the new temperatures.
190 : !-----------------------------------------------------------------
191 :
192 205765486 : diag = dfsurf_dT - kh
193 : rhs = dfsurf_dT*Tsf - fsurfn &
194 205765486 : - kh*Tbot
195 205765486 : Tsf = rhs / diag
196 :
197 : !-----------------------------------------------------------------
198 : ! Determine whether the computation has converged to an acceptable
199 : ! solution. Four conditions must be satisfied:
200 : !
201 : ! (1) Tsf <= 0 C.
202 : ! (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
203 : ! dTsf(niter-1) have magnitudes greater than puny, then
204 : ! dTsf(niter)/dTsf(niter-1) cannot be a negative number
205 : ! with magnitude greater than 0.5.
206 : ! (3) abs(dTsf) < Tsf_errmax
207 : ! (4) If Tsf = 0 C, then the downward turbulent/radiative
208 : ! flux, fsurfn, must be greater than or equal to the downward
209 : ! conductive flux, fcondtopn.
210 : !-----------------------------------------------------------------
211 :
212 : !-----------------------------------------------------------------
213 : ! Initialize convergence flag (true until proven false), dTsf,
214 : ! and temperature-averaging coefficients.
215 : ! Average only if test 1 or 2 fails.
216 : ! Initialize energy.
217 : !-----------------------------------------------------------------
218 :
219 205765486 : converged = .true.
220 205765486 : dTsf = Tsf - Tsf_start
221 205765486 : avg_Tsf = c0
222 :
223 : !-----------------------------------------------------------------
224 : ! Condition 1: check for Tsf > 0
225 : ! If Tsf > 0, set Tsf = 0 and leave converged=.true.
226 : !-----------------------------------------------------------------
227 :
228 205765486 : if (Tsf > puny) then
229 0 : Tsf = c0
230 0 : dTsf = -Tsf_start
231 :
232 : !-----------------------------------------------------------------
233 : ! Condition 2: check for oscillating Tsf
234 : ! If oscillating, average all temps to increase rate of convergence.
235 : ! It is possible that this may never occur.
236 : !-----------------------------------------------------------------
237 :
238 : elseif (niter > 1 & ! condition (2)
239 : .and. Tsf_start <= -puny &
240 : .and. abs(dTsf) > puny &
241 : .and. abs(dTsf_prev) > puny &
242 205765486 : .and. -dTsf/(dTsf_prev+puny*puny) > p5) then
243 :
244 0 : avg_Tsf = c1 ! average with starting temp
245 0 : dTsf = p5 * dTsf
246 0 : converged = .false.
247 : endif
248 :
249 : !-----------------------------------------------------------------
250 : ! If condition 2 failed, average new surface temperature with
251 : ! starting value.
252 : !-----------------------------------------------------------------
253 : Tsf = Tsf &
254 205765486 : + avg_Tsf * p5 * (Tsf_start - Tsf)
255 :
256 : !-----------------------------------------------------------------
257 : ! Condition 3: check for large change in Tsf
258 : !-----------------------------------------------------------------
259 :
260 205765486 : if (abs(dTsf) > Tsf_errmax) then
261 97132021 : converged = .false.
262 : endif
263 :
264 : !-----------------------------------------------------------------
265 : ! Condition 4: check for fsurfn < fcondtopn with Tsf > 0
266 : !-----------------------------------------------------------------
267 :
268 205765486 : fsurfn = fsurfn + dTsf*dfsurf_dT
269 205765486 : fcondtopn = kh * (Tsf-Tbot)
270 :
271 205765486 : if (Tsf > -puny .and. fsurfn < fcondtopn) then
272 0 : converged = .false.
273 : endif
274 :
275 205765486 : fcondbot = fcondtopn
276 :
277 205765486 : dTsf_prev = dTsf
278 :
279 : endif ! converged
280 :
281 : enddo ! temperature iteration niter
282 :
283 : !-----------------------------------------------------------------
284 : ! Check for convergence failures.
285 : !-----------------------------------------------------------------
286 108633465 : if (.not.converged) then
287 0 : write(warnstr,*) subname, 'Thermo iteration does not converge,'
288 0 : call icepack_warnings_add(warnstr)
289 0 : write(warnstr,*) subname, 'Ice thickness:', hilyr*nilyr
290 0 : call icepack_warnings_add(warnstr)
291 0 : write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr
292 0 : call icepack_warnings_add(warnstr)
293 0 : write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, &
294 0 : Tsf_errmax
295 0 : call icepack_warnings_add(warnstr)
296 0 : write(warnstr,*) subname, 'Tsf:', Tsf
297 0 : call icepack_warnings_add(warnstr)
298 0 : write(warnstr,*) subname, 'fsurfn:', fsurfn
299 0 : call icepack_warnings_add(warnstr)
300 0 : write(warnstr,*) subname, 'fcondtopn, fcondbot', &
301 0 : fcondtopn, fcondbot
302 0 : call icepack_warnings_add(warnstr)
303 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
304 0 : call icepack_warnings_add(subname//" zerolayer_temperature: Thermo iteration does not converge" )
305 0 : return
306 : endif
307 :
308 : !-----------------------------------------------------------------
309 : ! Check that if Tsfc < 0, then fcondtopn = fsurfn
310 : !-----------------------------------------------------------------
311 :
312 : if (l_zerolayerchecks) then
313 108633465 : if (Tsf < c0 .and. &
314 : abs(fcondtopn-fsurfn) > puny) then
315 :
316 0 : write(warnstr,*) subname, 'fcondtopn does not equal fsurfn,'
317 0 : call icepack_warnings_add(warnstr)
318 0 : write(warnstr,*) subname, 'Tsf=',Tsf
319 0 : call icepack_warnings_add(warnstr)
320 0 : write(warnstr,*) subname, 'fcondtopn=',fcondtopn
321 0 : call icepack_warnings_add(warnstr)
322 0 : write(warnstr,*) subname, 'fsurfn=',fsurfn
323 0 : call icepack_warnings_add(warnstr)
324 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
325 0 : call icepack_warnings_add(subname//" zerolayer_temperature: fcondtopn /= fsurfn" )
326 0 : return
327 : endif
328 : endif ! l_zerolayerchecks
329 :
330 : ! update fluxes that depend on Tsf
331 108633465 : flwoutn = flwoutn + dTsf_prev * dflwout_dT
332 108633465 : fsensn = fsensn + dTsf_prev * dfsens_dT
333 108633465 : flatn = flatn + dTsf_prev * dflat_dT
334 :
335 : end subroutine zerolayer_temperature
336 :
337 : !=======================================================================
338 :
339 : end module icepack_therm_0layer
340 :
341 : !=======================================================================
|
