1 | SUBROUTINE ice_sal_diff(nlay_i,kideb,kiut) |
---|
2 | |
---|
3 | !!------------------------------------------------------------------ |
---|
4 | !! *** ROUTINE ice_sal_diff *** |
---|
5 | !! |
---|
6 | !! ** Purpose : |
---|
7 | !! This routine computes new salinities in the ice |
---|
8 | !! |
---|
9 | !! ** Method : Vertical salinity profile computation |
---|
10 | !! Resolves brine transport equation |
---|
11 | !! |
---|
12 | !! ** Steps |
---|
13 | !! |
---|
14 | !! ** Arguments |
---|
15 | !! |
---|
16 | !! ** Inputs / Outputs |
---|
17 | !! |
---|
18 | !! ** External |
---|
19 | !! |
---|
20 | !! ** References : Vancop. et al., 2008 |
---|
21 | !! |
---|
22 | !! ** History : |
---|
23 | !! (06-2003) Martin Vancop. LIM1D |
---|
24 | !! (06-2008) Martin Vancop. BIO-LIM |
---|
25 | !! (09-2008) Martin Vancop. Explicit gravity drainage |
---|
26 | !! |
---|
27 | !!------------------------------------------------------------------ |
---|
28 | |
---|
29 | USE lib_fortran |
---|
30 | |
---|
31 | INCLUDE 'type.com' |
---|
32 | INCLUDE 'para.com' |
---|
33 | INCLUDE 'const.com' |
---|
34 | INCLUDE 'ice.com' |
---|
35 | INCLUDE 'thermo.com' |
---|
36 | |
---|
37 | REAL(8), DIMENSION(nlay_i) :: |
---|
38 | & z_ms_i , !: mass of salt times thickness |
---|
39 | & z_sbr_i !: brine salinity |
---|
40 | |
---|
41 | REAL(8), DIMENSION(nlay_i) :: !: dummy factors for tracer equation |
---|
42 | & za , !: all |
---|
43 | & zb , !: gravity drainage |
---|
44 | & zc , !: upward advective flow |
---|
45 | & ze , !: downward advective flow |
---|
46 | & zind , !: independent term in the tridiag system |
---|
47 | & zindtbis , !: |
---|
48 | & zdiagbis !: |
---|
49 | |
---|
50 | REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation |
---|
51 | & ztrid !: tridiagonal matrix |
---|
52 | |
---|
53 | REAL(8) :: |
---|
54 | & zdummy1 , !: dummy factors |
---|
55 | & zdummy2 , !: |
---|
56 | & zdummy3 , !: |
---|
57 | & zswitchs , !: switch for summer drainage |
---|
58 | & zeps = 1.0e-20 !: numerical limit |
---|
59 | |
---|
60 | ! Rayleigh number computation |
---|
61 | REAL(8) :: |
---|
62 | & ze_i_min , !: minimum brine volume |
---|
63 | & zcp , !: temporary scalar for sea ice specific heat |
---|
64 | & zk , !: temporary scalar for sea ice thermal conductivity |
---|
65 | & zalphara !: multiplicator for diffusivity |
---|
66 | |
---|
67 | REAL(8), DIMENSION(nlay_i) :: |
---|
68 | & zsigma , !: brine salinity at layer interfaces |
---|
69 | & zperm , !: permeability |
---|
70 | & zpermin , !: minimum permeability |
---|
71 | & zrhodiff , !: density difference |
---|
72 | & zlevel , !: height of the water column |
---|
73 | & zthdiff !: thermal diffusivity |
---|
74 | |
---|
75 | INTEGER :: |
---|
76 | & layer2 , !: layer loop index |
---|
77 | & indtr !: index of tridiagonal system |
---|
78 | |
---|
79 | CHARACTER(len=4) :: |
---|
80 | & bc = 'conc' !: Boundary condition 'conc' or 'flux' |
---|
81 | |
---|
82 | REAL(8) :: |
---|
83 | & z_ms_i_ini , !: initial mass of salt |
---|
84 | & z_ms_i_fin , !: final mass of salt |
---|
85 | & z_fs_b , !: basal flux of salt |
---|
86 | & z_fs_su , !: surface flux of salt |
---|
87 | & z_dms_i !: mass variation |
---|
88 | |
---|
89 | LOGICAL :: |
---|
90 | & ln_write , |
---|
91 | & ln_con , |
---|
92 | & ln_sal |
---|
93 | |
---|
94 | ln_write = .TRUE. ! write outputs |
---|
95 | ln_con = .TRUE. ! conservation check |
---|
96 | ln_sal = .TRUE. ! compute salinity variations or not |
---|
97 | |
---|
98 | IF ( ln_write ) THEN |
---|
99 | WRITE(numout,*) |
---|
100 | WRITE(numout,*) ' ** ice_sal_diff : ' |
---|
101 | WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' |
---|
102 | WRITE(numout,*) ' ln_sal = ', ln_sal |
---|
103 | WRITE(numout,*) ' ln_grd = ', ln_grd |
---|
104 | WRITE(numout,*) ' ln_flu = ', ln_flu |
---|
105 | WRITE(numout,*) ' ln_flo = ', ln_flo |
---|
106 | ENDIF |
---|
107 | WRITE(numout,*) " nlay_i : ", nlay_i |
---|
108 | |
---|
109 | IF ( ln_sal ) THEN |
---|
110 | ! |
---|
111 | !------------------------------------------------------------------------------| |
---|
112 | ! 1) Initialization |
---|
113 | !------------------------------------------------------------------------------| |
---|
114 | ! |
---|
115 | IF ( ln_write ) THEN |
---|
116 | WRITE(numout,*) ' - Initialization ... ' |
---|
117 | ENDIF |
---|
118 | |
---|
119 | DO 10 ji = kideb, kiut |
---|
120 | |
---|
121 | ! brine diffusivity |
---|
122 | diff_br(:) = 0.0 |
---|
123 | |
---|
124 | !--------------------------- |
---|
125 | ! Brine volume and salinity |
---|
126 | !--------------------------- |
---|
127 | DO layer = 1, nlay_i |
---|
128 | e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) |
---|
129 | & - tpw ) |
---|
130 | z_sbr_i(layer) = s_i_b(ji,layer) / e_i_b(layer) |
---|
131 | END DO |
---|
132 | |
---|
133 | !-------------------- |
---|
134 | ! Conservation check |
---|
135 | !-------------------- |
---|
136 | IF ( ln_con ) THEN |
---|
137 | CALL ice_sal_column( kideb , kiut , z_ms_i_ini , |
---|
138 | & s_i_b(1,1:nlay_i), |
---|
139 | & deltaz_i_phy, nlay_i, .FALSE. ) |
---|
140 | ENDIF ! ln_con |
---|
141 | |
---|
142 | IF ( ln_write ) THEN |
---|
143 | WRITE(numout,*) ' nlay_i : ', nlay_i |
---|
144 | WRITE(numout,*) ' kideb : ', kideb |
---|
145 | WRITE(numout,*) ' kiut : ', kiut |
---|
146 | WRITE(numout,*) ' ' |
---|
147 | WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer), |
---|
148 | & layer = 1, nlay_i ) |
---|
149 | WRITE(numout,*) ' z_i_phy : ', ( z_i_phy(layer), |
---|
150 | & layer = 1, nlay_i ) |
---|
151 | WRITE(numout,*) ' s_i_b : ', ( s_i_b (ji,layer), |
---|
152 | & layer = 1, nlay_i ) |
---|
153 | WRITE(numout,*) ' t_i_b : ', ( t_i_b (ji,layer), |
---|
154 | & layer = 1, nlay_i ) |
---|
155 | WRITE(numout,*) ' e_i_b : ', ( e_i_b (layer), |
---|
156 | & layer = 1, nlay_i ) |
---|
157 | WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer), |
---|
158 | & layer = 1, nlay_i ) |
---|
159 | WRITE(numout,*) |
---|
160 | ENDIF ! ln_write |
---|
161 | |
---|
162 | 10 CONTINUE |
---|
163 | |
---|
164 | ! |
---|
165 | !------------------------------------------------------------------------------| |
---|
166 | ! 2) Rayleigh-number-based diffusivity |
---|
167 | !------------------------------------------------------------------------------| |
---|
168 | ! |
---|
169 | ! Diffusivity is a function of the local Rayleigh number |
---|
170 | ! see Notz and Worster, JGR 2008 |
---|
171 | ! Diffusivity, layer represents the interface' |
---|
172 | ! between layer and layer-1 ' |
---|
173 | ! |
---|
174 | IF ( ln_write ) THEN |
---|
175 | WRITE(numout,*) ' - Rayleigh-number based diffusivity ... ' |
---|
176 | WRITE(numout,*) ' ' |
---|
177 | ENDIF |
---|
178 | |
---|
179 | DO 20 ji = kideb, kiut |
---|
180 | |
---|
181 | !----------------------------------------- |
---|
182 | ! Brine salinity between layer interfaces |
---|
183 | !----------------------------------------- |
---|
184 | DO layer = 1, nlay_i - 1 |
---|
185 | zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
---|
186 | & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / |
---|
187 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw |
---|
188 | zsigma(layer) = - zdummy1 / tmut |
---|
189 | END DO |
---|
190 | zsigma(nlay_i) = - ( t_i_b(ji,nlay_i) - tpw ) / tmut |
---|
191 | |
---|
192 | !-------------------- |
---|
193 | ! Density difference |
---|
194 | !-------------------- |
---|
195 | DO layer = 1, nlay_i |
---|
196 | zrhodiff(layer) = - beta_ocs * ( oce_sal - zsigma(layer) ) |
---|
197 | END DO |
---|
198 | |
---|
199 | !------------------------------------------ |
---|
200 | ! Minimum permeability under current level |
---|
201 | !------------------------------------------ |
---|
202 | DO layer = 1, nlay_i |
---|
203 | ze_i_min = 99999.0 |
---|
204 | DO layer2 = layer, nlay_i |
---|
205 | ze_i_min = MIN( ze_i_min , e_i_b(layer2) ) |
---|
206 | zpermin(layer) = 1.0e-17 * ( ( 1000. * ze_i_min )**3.1 ) |
---|
207 | END DO |
---|
208 | END DO ! layer |
---|
209 | |
---|
210 | !------------------------------------------------ |
---|
211 | ! length of the water column under current level |
---|
212 | !------------------------------------------------ |
---|
213 | DO layer = nlay_i, 1, -1 |
---|
214 | zlevel(layer) = 0.0 |
---|
215 | DO layer2 = layer, nlay_i |
---|
216 | zlevel(layer) = zlevel(layer) + deltaz_i_phy(layer2) |
---|
217 | END DO |
---|
218 | END DO |
---|
219 | zlevel(nlay_i) = deltaz_i_phy(nlay_i) / 2.0 |
---|
220 | |
---|
221 | !--------------------- |
---|
222 | ! Thermal diffusivity |
---|
223 | !--------------------- |
---|
224 | zkimin = 0.1 |
---|
225 | DO layer = 1, nlay_i - 1 |
---|
226 | zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
---|
227 | & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / |
---|
228 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw |
---|
229 | zdummy2 = s_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
---|
230 | & ( s_i_b(ji,layer+1) - s_i_b(ji,layer) ) / |
---|
231 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) |
---|
232 | zcp = cpg + lfus * tmut * zdummy2 / |
---|
233 | & MAX( zdummy1 * zdummy1 , zeps ) |
---|
234 | zk = xkg + betak1 * zdummy2 / |
---|
235 | & MIN( -zeps , zdummy1 ) - betak2 * zdummy1 |
---|
236 | zk = MAX( zk, zkimin ) |
---|
237 | zthdiff(layer) = zk / ( rhog * zcp ) |
---|
238 | END DO |
---|
239 | |
---|
240 | zcp = cpg + lfus * tmut * s_i_b(ji,nlay_i) / |
---|
241 | & MAX( ( t_i_b(ji,nlay_i) - tpw ) * ( t_i_b(ji,nlay_i) - tpw ), |
---|
242 | & zeps ) |
---|
243 | zk = xkg + betak1 * s_i_b(ji,nlay_i) / |
---|
244 | & MIN( -zeps , t_i_b(ji,nlay_i) - tpw ) |
---|
245 | & - betak2 * ( t_i_b(ji,nlay_i) - tpw ) |
---|
246 | zk = MAX( zk, zkimin ) |
---|
247 | zthdiff(nlay_i) = zk / ( rhog * zcp ) |
---|
248 | |
---|
249 | !----------------- |
---|
250 | ! Rayleigh number |
---|
251 | !----------------- |
---|
252 | DO layer = 1, nlay_i |
---|
253 | rayleigh(layer) = gpes * MAX(zrhodiff(layer),0.0) * |
---|
254 | & zpermin(layer) * zlevel(layer) / |
---|
255 | & ( zthdiff(layer) * visc_br ) |
---|
256 | END DO |
---|
257 | |
---|
258 | !------------------- |
---|
259 | ! Brine Diffusivity |
---|
260 | !------------------- |
---|
261 | DO layer = 1, nlay_i |
---|
262 | zalphara = ( TANH( ra_smooth * ( rayleigh(layer) - ra_c ) ) |
---|
263 | & + 1 ) / 2.0 |
---|
264 | diff_br(layer) = ( 1.0 - zalphara ) * d_br_mol + |
---|
265 | & zalphara * ( d_br_tur ) |
---|
266 | IF ( .NOT. ln_grd ) diff_br(layer) = 0. |
---|
267 | END DO |
---|
268 | |
---|
269 | |
---|
270 | IF ( ln_write ) THEN |
---|
271 | WRITE(numout,*) ' zsigma : ', ( zsigma(layer), |
---|
272 | & layer = 1, nlay_i) |
---|
273 | WRITE(numout,*) ' zrhodiff : ', ( zrhodiff(layer), |
---|
274 | & layer = 1, nlay_i ) |
---|
275 | WRITE(numout,*) ' zpermin : ', ( zpermin(layer), |
---|
276 | & layer = 1, nlay_i ) |
---|
277 | WRITE(numout,*) ' zthdiff : ', ( zthdiff(layer), |
---|
278 | & layer = 1, nlay_i ) |
---|
279 | WRITE(numout,*) ' zlevel : ', ( zlevel(layer), |
---|
280 | & layer = 1, nlay_i ) |
---|
281 | WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer), |
---|
282 | & layer = 1, nlay_i ) |
---|
283 | WRITE(numout,*) ' diff_br : ', ( diff_br(layer), |
---|
284 | & layer = 1, nlay_i ) |
---|
285 | WRITE(numout,*) |
---|
286 | ENDIF |
---|
287 | |
---|
288 | 20 CONTINUE |
---|
289 | |
---|
290 | ! |
---|
291 | !------------------------------------------------------------------------------| |
---|
292 | ! 3) Flooding and flushing velocities |
---|
293 | !------------------------------------------------------------------------------| |
---|
294 | ! |
---|
295 | IF ( ln_write ) THEN |
---|
296 | WRITE(numout,*) ' - Flooding and flushing velocities ' |
---|
297 | WRITE(numout,*) ' ' |
---|
298 | ENDIF |
---|
299 | |
---|
300 | DO 30 ji = kideb, kiut |
---|
301 | !----------------------- |
---|
302 | ! Permeability switches |
---|
303 | !----------------------- |
---|
304 | ! Permeability switch = 1 if brine volume fraction > e_thr_flu |
---|
305 | zswitch_per = 1.0 |
---|
306 | zbvmin = 1.0 |
---|
307 | DO layer = 1, nlay_i |
---|
308 | zbvmin = MIN( e_i_b(layer) , zbvmin ) ! minimum brine volume |
---|
309 | END DO |
---|
310 | IF ( zbvmin .LT. e_thr_flu ) zswitch_per = 0.0 |
---|
311 | |
---|
312 | ! Summer switch = 1 if Tsu ge tpw and min brine volume superior than e_thr_flu |
---|
313 | zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 if winter, 1 if summer |
---|
314 | zswitchs = zswitchs * zswitch_per |
---|
315 | |
---|
316 | !------------------- |
---|
317 | ! Flooding velocity |
---|
318 | !------------------- |
---|
319 | zdhitot = dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) |
---|
320 | zdhstot = dh_s_melt(ji) |
---|
321 | w_flood = ( ( rho0 - rhog ) / rho0 * zdhitot - |
---|
322 | & rhon / rho0 * zdhstot ) / ddtb * zswitch_per |
---|
323 | IF ( w_flood .GT. 0 ) THEN |
---|
324 | z_flood = 0.0 |
---|
325 | ELSE |
---|
326 | z_flood = 1.0 |
---|
327 | ENDIF |
---|
328 | IF ( .NOT. ln_flo ) w_flood = 0. |
---|
329 | |
---|
330 | !------------------ |
---|
331 | ! Percolating flux |
---|
332 | !------------------ |
---|
333 | ! Percolating flow ( rho dh * beta * switch / rho0 ) |
---|
334 | qsummer = ( - rhog * MIN ( dh_i_surf(ji) , 0.0 ) |
---|
335 | & - rhon * MIN ( dh_s_tot(ji) , 0.0 ) ) |
---|
336 | qsummer = qsummer * flu_beta * zswitchs / 1000.0 ! 1000 is a ref density for brine |
---|
337 | |
---|
338 | w_flush = qsummer / ddtb / e_i_b(1) |
---|
339 | |
---|
340 | IF ( .NOT. ln_flu ) THEN |
---|
341 | w_flush = 0. |
---|
342 | qsummer = 0. |
---|
343 | ENDIF |
---|
344 | |
---|
345 | IF ( ln_write ) THEN |
---|
346 | WRITE(numout,*) ' zswitchs : ', zswitchs |
---|
347 | WRITE(numout,*) ' w_flush : ', w_flush |
---|
348 | WRITE(numout,*) ' w_flood : ', w_flood |
---|
349 | ENDIF |
---|
350 | |
---|
351 | 30 CONTINUE |
---|
352 | ! |
---|
353 | !------------------------------------------------------------------------------| |
---|
354 | ! 4) Compute dummy factors for tracer diffusion equation |
---|
355 | !------------------------------------------------------------------------------| |
---|
356 | ! |
---|
357 | IF ( ln_write ) THEN |
---|
358 | WRITE(numout,*) ' - Compute dummy factors for tracer diffusion' |
---|
359 | WRITE(numout,*) ' ' |
---|
360 | ENDIF |
---|
361 | |
---|
362 | DO 40 ji = kideb, kiut |
---|
363 | !---------------- |
---|
364 | ! za factors |
---|
365 | !---------------- |
---|
366 | DO layer = 1, nlay_i |
---|
367 | za(layer) = ddtb / ( deltaz_i_phy(layer) * e_i_b(layer) ) |
---|
368 | END DO |
---|
369 | |
---|
370 | !-------------------- |
---|
371 | ! zb, zc, ze factors |
---|
372 | !-------------------- |
---|
373 | DO layer = 1, nlay_i - 1 |
---|
374 | ! interpolate brine volume at the interface between layers |
---|
375 | zdummy1 = ( e_i_b(layer + 1 ) - e_i_b(layer) ) / |
---|
376 | & ( z_i_phy(layer + 1) - z_i_phy(layer) ) |
---|
377 | zdummy2 = deltaz_i_phy(layer) / 2.0 |
---|
378 | zdummy3 = e_i_b(layer) + zdummy1 * zdummy2 |
---|
379 | zb(layer) = zdummy3 * diff_br(layer) / |
---|
380 | & ( z_i_phy(layer + 1) - z_i_phy(layer) ) |
---|
381 | zc(layer) = w_flood * zdummy3 * z_flood |
---|
382 | ze(layer) = ( w_flood * ( 1. - z_flood ) + w_flush ) * zdummy3 |
---|
383 | ! old qsummer scheme |
---|
384 | ! ze(layer) = ( w_flood * ( 1. - z_flood ) ) * zdummy3 + |
---|
385 | ! & zswitchs * qsummer / ddtb |
---|
386 | END DO |
---|
387 | |
---|
388 | ! Fixed boundary condition (imposed cc.) |
---|
389 | zb(nlay_i) = 2. * e_i_b(nlay_i) / |
---|
390 | & deltaz_i_phy(nlay_i) * diff_br(nlay_i) |
---|
391 | zc(nlay_i) = w_flood * e_i_b(nlay_i) * z_flood |
---|
392 | ze(nlay_i) = ( w_flood * ( 1. - z_flood ) + w_flush ) * |
---|
393 | & e_i_b(nlay_i) |
---|
394 | ! ze(nlay_i) = ( w_flood * ( 1. - z_flood ) ) * e_i_b(nlay_i) + |
---|
395 | ! & zswitchs * qsummer / ddtb |
---|
396 | |
---|
397 | IF ( ln_write ) THEN |
---|
398 | WRITE(numout,*) |
---|
399 | WRITE(numout,*) ' -Dummy factors ' |
---|
400 | WRITE(numout,*) ' za : ', ( za (layer), |
---|
401 | & layer = 1, nlay_i) |
---|
402 | WRITE(numout,*) ' zb : ', ( zb (layer), |
---|
403 | & layer = 1, nlay_i) |
---|
404 | WRITE(numout,*) ' zc : ', ( zc (layer), |
---|
405 | & layer = 1, nlay_i) |
---|
406 | WRITE(numout,*) ' ze : ', ( ze (layer), |
---|
407 | & layer = 1, nlay_i) |
---|
408 | ENDIF |
---|
409 | |
---|
410 | 40 CONTINUE |
---|
411 | ! |
---|
412 | !----------------------------------------------------------------------- |
---|
413 | ! 5) Tridiagonal system terms for tracer diffusion equation |
---|
414 | !----------------------------------------------------------------------- |
---|
415 | ! |
---|
416 | DO 50 ji = kideb, kiut |
---|
417 | |
---|
418 | !---------------- |
---|
419 | ! first equation |
---|
420 | !---------------- |
---|
421 | ztrid(1,1) = 0.0 |
---|
422 | ztrid(1,2) = 1.0 + za(1) * ( zb(1) + ze(1) ) |
---|
423 | ztrid(1,3) = za(1) * ( -zb(1) + zc(1) ) |
---|
424 | zind(1) = z_sbr_i(1) |
---|
425 | |
---|
426 | !----------------- |
---|
427 | ! inner equations |
---|
428 | !----------------- |
---|
429 | DO layer = 2, nlay_i - 1 |
---|
430 | ztrid(layer,1) = - za(layer) * ( zb(layer-1) + ze(layer-1) ) |
---|
431 | ztrid(layer,2) = 1.0 + za(layer) * ( zb(layer) + |
---|
432 | & ze(layer) + zb(layer-1) - zc(layer-1) ) |
---|
433 | ztrid(layer,3) = za(layer) * ( -zb(layer) + zc(layer) ) |
---|
434 | zind(layer) = z_sbr_i(layer) |
---|
435 | END DO |
---|
436 | |
---|
437 | !---------------- |
---|
438 | ! last equation |
---|
439 | !---------------- |
---|
440 | WRITE(numout,*) " nlay_i : ", nlay_i |
---|
441 | ztrid(nlay_i,1) = - za(nlay_i) * ( zb(nlay_i-1) + ze(nlay_i-1) ) |
---|
442 | ztrid(nlay_i,2) = 1.0 + za(nlay_i) * ( zb(nlay_i) + |
---|
443 | & ze(nlay_i) + zb(nlay_i-1) - zc(nlay_i-1) ) |
---|
444 | ztrid(nlay_i,3) = 0. |
---|
445 | zind(nlay_i) = z_sbr_i(nlay_i) + za(nlay_i) * ( zb(nlay_i) |
---|
446 | & - zc(nlay_i) ) * oce_sal |
---|
447 | |
---|
448 | IF ( ln_write ) THEN |
---|
449 | WRITE(numout,*) |
---|
450 | WRITE(numout,*) ' -Tridiag terms, winter ... ' |
---|
451 | WRITE(numout,*) |
---|
452 | DO layer = 1, nlay_i |
---|
453 | WRITE(numout,*) ' layer : ', layer |
---|
454 | WRITE(numout,*) ' ztrid : ', ztrid(layer,1), |
---|
455 | & ztrid(layer,2), ztrid(layer,3) |
---|
456 | WRITE(numout,*) ' zind : ', zind(layer) |
---|
457 | END DO |
---|
458 | ENDIF |
---|
459 | |
---|
460 | 50 CONTINUE |
---|
461 | |
---|
462 | ! |
---|
463 | !----------------------------------------------------------------------- |
---|
464 | ! 6) Solving the tridiagonal system |
---|
465 | !----------------------------------------------------------------------- |
---|
466 | ! |
---|
467 | DO 60 ji = kideb, kiut |
---|
468 | |
---|
469 | ! The tridiagonal system is solved with Gauss elimination |
---|
470 | ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, |
---|
471 | ! McGraw-Hill 1984. |
---|
472 | |
---|
473 | zindtbis(1) = zind(1) |
---|
474 | zdiagbis(1) = ztrid(1,2) |
---|
475 | |
---|
476 | DO layer = 2, nlay_i |
---|
477 | zdiagbis(layer) = ztrid(layer,2) - ztrid(layer,1) * |
---|
478 | & ztrid(layer-1,3) / zdiagbis(layer-1) |
---|
479 | zindtbis(layer) = zind(layer) - ztrid(layer,1) * |
---|
480 | & zindtbis(layer-1) / zdiagbis(layer-1) |
---|
481 | END DO |
---|
482 | |
---|
483 | ! Recover brine salinity |
---|
484 | z_sbr_i(nlay_i) = zindtbis(nlay_i) / zdiagbis(nlay_i) |
---|
485 | DO layer = nlay_i - 1 , 1 , -1 |
---|
486 | z_sbr_i(layer) = ( zindtbis(layer) - ztrid(layer,3)* |
---|
487 | & z_sbr_i(layer+1)) / zdiagbis(layer) |
---|
488 | END DO |
---|
489 | ! Recover ice salinity |
---|
490 | DO layer = 1, nlay_i |
---|
491 | sn_i_b(layer) = z_sbr_i(layer) * e_i_b(layer) |
---|
492 | ! s_i_b(ji,layer) = z_sbr_i(layer) * e_i_b(layer) |
---|
493 | END DO |
---|
494 | |
---|
495 | IF ( ln_write ) THEN |
---|
496 | WRITE(numout,*) |
---|
497 | WRITE(numout,*) ' -Solving the tridiagonal system ... ' |
---|
498 | WRITE(numout,*) |
---|
499 | WRITE(numout,*) ' zdiagbis: ', ( zdiagbis(layer) , |
---|
500 | & layer = 1, nlay_i ) |
---|
501 | WRITE(numout,*) ' zindtbis: ', ( zdiagbis(layer) , |
---|
502 | & layer = 1, nlay_i ) |
---|
503 | WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i(layer) , |
---|
504 | & layer = 1, nlay_i ) |
---|
505 | WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) , |
---|
506 | & layer = 1, nlay_i ) |
---|
507 | ENDIF |
---|
508 | |
---|
509 | 60 CONTINUE |
---|
510 | ! |
---|
511 | !----------------------------------------------------------------------- |
---|
512 | ! 7) Mass of salt conserved ? |
---|
513 | !----------------------------------------------------------------------- |
---|
514 | ! |
---|
515 | DO 70 ji = kideb, kiut |
---|
516 | |
---|
517 | ! Final mass of salt |
---|
518 | CALL ice_sal_column( kideb , kiut , z_ms_i_fin , |
---|
519 | & sn_i_b(1:nlay_i), |
---|
520 | & deltaz_i_phy, nlay_i, .FALSE. ) |
---|
521 | |
---|
522 | ! Bottom flux ( positive upwards for conservation routine ) |
---|
523 | zfb = - e_i_b(nlay_i) * |
---|
524 | & ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) * |
---|
525 | & ( z_sbr_i(nlay_i) - oce_sal ) + w_flood * ( z_flood * |
---|
526 | & oce_sal + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) ) |
---|
527 | & - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i) |
---|
528 | ! & - qsummer * z_sbr_i(nlay_i) / ddtb |
---|
529 | |
---|
530 | fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards |
---|
531 | IF ( ln_write ) THEN |
---|
532 | WRITE(numout,*) ' fsb : ', fsb |
---|
533 | WRITE(numout,*) |
---|
534 | ENDIF |
---|
535 | |
---|
536 | ! Surface flux of salt |
---|
537 | zfsu = 0.0 |
---|
538 | |
---|
539 | ! conservation check |
---|
540 | zerror = 1.0e-15 |
---|
541 | CALL ice_sal_conserv(kideb,kiut,'ice_sal_diff : ',zerror, |
---|
542 | & z_ms_i_ini,z_ms_i_fin, |
---|
543 | & zfb , zfsu , ddtb) |
---|
544 | |
---|
545 | 70 CONTINUE |
---|
546 | |
---|
547 | ENDIF ! ln_sal |
---|
548 | ! |
---|
549 | !------------------------------------------------------------------------------| |
---|
550 | ! End of la sous-routine |
---|
551 | WRITE(numout,*) |
---|
552 | END SUBROUTINE |
---|