Context Navigation

source: branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_sal_adv.f @ 29

Last change on this file since 29 was 29, checked in by vancop, 8 years ago
Add ice_sal_adv routine for GN13
File size: 14.8 KB

Line
1	SUBROUTINE ice_sal_adv(nlay_i,kideb,kiut)
2
3	!!------------------------------------------------------------------
4	!! * ROUTINE ice_sal_adv *
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	& zRae , !: effective Ra
46	& ze , !: downward advective flow
47	& zind , !: independent term in the tridiag system
48	& zindtbis , !:
49	& zdiagbis !:
50
51
52	REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation
53	& ztrid !: tridiagonal matrix
54
55	REAL(8) ::
56	& zdummy1 , !: dummy factors
57	& zdummy2 , !:
58	& zdummy3 , !:
59	& zswitchs , !: switch for summer drainage
60	& zeps = 1.0e-20 !: numerical limit
61
62	! Rayleigh number computation
63	REAL(8) ::
64	& ze_i_min , !: minimum brine volume
65	& zcp , !: temporary scalar for sea ice specific heat
66	& zk , !: temporary scalar for sea ice thermal conductivity
67	& zalphara !: multiplicator for diffusivity
68
69	REAL(8), DIMENSION(nlay_i) ::
70	& zsigma , !: brine salinity at layer interfaces
71	& zperm , !: permeability
72	& zpermin , !: minimum permeability
73	& zrhodiff , !: density difference
74	& zlevel , !: height of the water column
75	& zthdiff !: thermal diffusivity
76
77	REAL(8), DIMENSION(nlay_i+1) ::
78	& z_sbr_int !: brine salinity at layer interfaces
79
80	INTEGER ::
81	& layer2 , !: layer loop index
82	& indtr !: index of tridiagonal system
83
84	CHARACTER(len=4) ::
85	& bc = 'conc' !: Boundary condition 'conc' or 'flux'
86
87	REAL(8) ::
88	& z_ms_i_ini , !: initial mass of salt
89	& z_ms_i_fin , !: final mass of salt
90	& z_fs_b , !: basal flux of salt
91	& z_fs_su , !: surface flux of salt
92	& z_dms_i !: mass variation
93
94	LOGICAL ::
95	& ln_write ,
96	& ln_con ,
97	& ln_sal
98
99	ln_write = .TRUE. ! write outputs
100	ln_con = .TRUE. ! conservation check
101	ln_sal = .TRUE. ! compute salinity variations or not
102
103	IF ( ln_write ) THEN
104	WRITE(numout,*)
105	WRITE(numout,) ' * ice_sal_adv : '
106	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ '
107	WRITE(numout,*) ' ln_sal = ', ln_sal
108	WRITE(numout,*) ' ln_grd = ', ln_grd
109	WRITE(numout,*) ' ln_flu = ', ln_flu
110	WRITE(numout,*) ' ln_flo = ', ln_flo
111	WRITE(numout,*) ' c_gravdr = ', c_gravdr
112	WRITE(numout,*) ' c_sbr = ', c_sbr
113	WRITE(numout,*) ' c_perm = ', c_perm
114	WRITE(numout,*) ' c_permeff= ', c_permeff
115	ENDIF
116
117	ji = 1
118
119	IF ( ln_sal ) THEN
120	!
121	!------------------------------------------------------------------------------\|
122	! 1) Initialization
123	!------------------------------------------------------------------------------\|
124	!
125	IF ( ln_write ) THEN
126	WRITE(numout,*) ' - Initialization ... '
127	ENDIF
128
129	! ! brine diffusivity
130	! diff_br(:) = 0.0
131
132	! Darcy velocity
133	w_adv_br(:) = 0.0
134
135	!--------------------
136	! Conservation check
137	!--------------------
138	IF ( ln_con ) THEN
139	CALL ice_sal_column( kideb , kiut , z_ms_i_ini ,
140	& s_i_b(1,1:nlay_i),
141	& deltaz_i_phy, nlay_i, .FALSE. )
142	ENDIF ! ln_con
143
144	IF ( ln_write ) THEN
145	WRITE(numout,*) ' '
146	WRITE(numout,*) ' nlay_i : ', nlay_i
147	WRITE(numout,*) ' kideb : ', kideb
148	WRITE(numout,*) ' kiut : ', kiut
149	WRITE(numout,*) ' '
150	WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer),
151	& layer = 1, nlay_i )
152	WRITE(numout,*) ' z_i_phy : ', ( z_i_phy(layer),
153	& layer = 1, nlay_i )
154	WRITE(numout,*) ' s_i_b : ', ( s_i_b (ji,layer),
155	& layer = 1, nlay_i )
156	WRITE(numout,*) ' t_i_b : ', ( t_i_b (ji,layer),
157	& layer = 1, nlay_i )
158	WRITE(numout,*) ' t_i_int : ', ( t_i_int (ji,layer),
159	& layer = 1, nlay_i+1 )
160	WRITE(numout,*)
161	ENDIF ! ln_write
162
163	!
164	!------------------------------------------------------------------------------\|
165	! 2) Brine salinity at layer mid points and interfaces, brine fraction
166	!------------------------------------------------------------------------------\|
167	!
168	DO layer = 1, nlay_i
169	e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer)
170	& - tpw )
171	END DO
172
173	IF ( c_sbr .EQ. 'LIN' ) THEN ! Linear liquidus
174
175	DO layer = 1, nlay_i + 1
176	zTc = t_i_int(ji,layer) - tpw
177	z_sbr_int(layer) = - zTc / tmut !--- interfacial value
178	END DO
179
180	DO layer = 1, nlay_i
181	zTc = t_i_b(ji,layer) - tpw
182	z_sbr_i(layer) = - zTc / tmut !--- mid-point value
183	END DO
184
185	ENDIF
186
187	IF ( c_sbr .EQ. 'WEA' ) THEN ! Weast (1971) 3rd order liquidus
188
189	DO layer = 1, nlay_i + 1
190	zTc = t_i_int(ji,layer) - tpw
191	z_sbr_int(layer) = -17.6zTc -0.389zTc*2. -0.00362zTc**3.
192	END DO
193
194	DO layer = 1, nlay_i
195	zTc = t_i_b(ji,layer) - tpw
196	z_sbr_i(layer) = -17.6zTc -0.389zTc*2. -0.00362zTc**3.
197	END DO
198
199	ENDIF
200
201	IF ( c_sbr .EQ. 'NTZ' ) THEN ! Notz (2005) 3rd order liquidus
202
203	DO layer = 1, nlay_i + 1
204	zTc = t_i_int(ji,layer) - tpw
205	z_sbr_int(layer) = -21.4zTc - 0.886zTc*2. - 0.017zTc**3.
206	END DO
207
208	DO layer = 1, nlay_i
209	zTc = t_i_b(ji,layer) - tpw
210	z_sbr_i(layer) = -21.4zTc - 0.886zTc*2. - 0.017zTc**3.
211	END DO
212
213	ENDIF
214
215	IF ( ln_write ) THEN
216	WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer),
217	& layer = 1, nlay_i )
218	WRITE(numout,*) ' z_sbr_int : ', ( z_sbr_int (layer),
219	& layer = 1, nlay_i + 1 )
220	WRITE(numout,*) ' e_i_b : ', ( e_i_b (layer),
221	& layer = 1, nlay_i )
222	ENDIF ! ln_write
223	!
224	!------------------------------------------------------------------------------\|
225	! 3) Effective permeability at layer mid-points
226	!------------------------------------------------------------------------------\|
227	!
228	IF ( c_permeff .EQ. 'HAR' ) THEN ! Harmonic mean
229	WRITE(numout,*)
230	ENDIF
231
232	IF ( c_permeff .EQ. 'MIN' ) THEN ! Minimum permeability
233
234	DO layer = 1, nlay_i
235	ze_i_min = 99999.0
236	DO layer2 = layer, nlay_i
237	ze_i_min = MIN( ze_i_min , e_i_b(layer2) )
238	IF ( c_perm .EQ. 'FRE' ) ! Freitag (1999)
239	& perm_eff(layer) = 1.995e-8 * ze_i_min**3.1
240	IF ( c_perm .EQ. 'RJW' ) ! Rees-Jones and Worster (2014)
241	& perm_eff(layer) = 1.0e-8 * ze_i_min**3
242	END DO
243	END DO ! layer
244
245	ENDIF
246
247	!
248	!------------------------------------------------------------------------------\|
249	! 4) Rayleigh number at mid-points (see Vancoppenolle et al TCD2013)
250	!------------------------------------------------------------------------------\|
251	!
252	z1 = gpes / ( thdiff_br * visc_br );
253
254	DO layer = 1, nlay_i
255
256	z2 = beta_ocs * ( z_sbr_i(layer) - s_w ) ! Delta rho
257	z3 = ht_i_b(ji) - z_i_phy(layer)
258	z4 = perm_eff(layer)
259
260	rayleigh(layer) = MAX( z1 * z2 * z3 * z4, 0.0)
261
262	END DO
263
264	! IF ( ln_write ) THEN
265	! WRITE(numout,*)
266	! WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer),
267	! & layer = 1, nlay_i )
268	! WRITE(numout,*) ' diff_br : ', ( diff_br(layer),
269	! & layer = 1, nlay_i )
270	! WRITE(numout,*)
271
272
273	!
274	!------------------------------------------------------------------------------\|
275	! 6) Brine velocity at layer mid-points
276	!------------------------------------------------------------------------------\|
277	!
278	! c_wadv = 'GN' or 'RW'
279
280	alpha_GN = 1.56e-3
281	rho_br_GN = 1020.
282	Rc_GN = 1.01
283	!alpha_GN = 1.0e-3
284	!Rc_GN = 1.0
285
286	w_adv_br(:) = 0.0
287	zRae(:) = 0.
288	DO layer = 1, nlay_i
289	zRae(layer) = MAX( rayleigh(layer) - Rc_GN, 0.) ! correction from the orig scheme
290	END DO
291
292	DO layer = 1, nlay_i
293	w_adv_br(layer) = - alpha_GN / rho_br_GN *
294	& SUM ( zRae(1:layer)*deltaz_i_phy(1:layer) )
295	END DO
296
297	IF ( ln_write ) THEN
298	WRITE(numout,*)
299	WRITE(numout,*) ' w_adv_br : ', ( w_adv_br(layer),
300	& layer = 1, nlay_i )
301	ENDIF
302	!
303	!------------------------------------------------------------------------------\|
304	! 7) New salinities
305	!------------------------------------------------------------------------------\|
306	!
307
308	DO layer = 1, nlay_i
309	za(layer) = w_adv_br(layer) * ddtb / deltaz_i_phy(layer)
310	END DO
311
312	! first layer
313	sn_i_b(1) = z_sbr_i(1) * ( e_i_b(1) + za(1) ) +
314	& z_sbr_i(2) * ( - za(1) )
315
316
317	! inner layers
318	DO layer = 2, nlay_i - 1
319	sn_i_b(layer) = z_sbr_i(layer-1) * ( za(layer)/2. ) +
320	& z_sbr_i(layer) * e_i_b(layer) +
321	& z_sbr_i(layer+1) * ( - za(layer)/2. )
322	END DO
323
324	! lowermost layer
325	sn_i_b(nlay_i) = z_sbr_i(nlay_i-1) * ( za(nlay_i)/2. ) +
326	& z_sbr_i(nlay_i) * ( e_i_b(nlay_i) +
327	& za(nlay_i)/2. ) - za(nlay_i) * s_w
328
329	IF ( ln_write ) THEN
330	WRITE(numout,*)
331	WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) ,
332	& layer = 1, nlay_i )
333	ENDIF
334
335	!
336	!-----------------------------------------------------------------------
337	! 8) Salt flux diagnostics
338	!-----------------------------------------------------------------------
339	!
340	! diagnostics to compute a salt budget
341	! fsupw2 upwelling salt flux into layer k
342	! fsupw1 upwelling salt flux from layer k into layer k-1
343	! fsdwn downwelling salt flux from layer k (brine channels)
344	! dzdsdt = rate of change in salt content
345
346	!--- First layer ---
347	fsupw2(1) = -1./4 * ( w_adv_br(1) + w_adv_br(2) ) *
348	& ( z_sbr_i(1) + z_sbr_i(2) )
349
350	fsupw1(1) = 0.
351
352	dzdsdt(1) = ( sn_i_b(1) - s_i_b(ji,1) ) * deltaz_i_phy(1) / ddtb
353
354	fsdwn(1) = 0.25( z_sbr_i(1)( -5.*w_adv_br(1) - w_adv_br(2) ) +
355	& z_sbr_i(2)( 3.w_adv_br(1) - w_adv_br(2) ) )
356
357	!--- Layers 2->N-1 ---
358	DO layer = 2, nlay_i-1
359
360	fsupw2(layer) = -0.25*
361	& ( w_adv_br(layer) + w_adv_br(layer+1) )*
362	& ( z_sbr_i(layer) + z_sbr_i(layer+1) )
363
364	fsupw1(layer) = -0.25*
365	& ( w_adv_br(layer-1) + w_adv_br(layer) )*
366	& ( z_sbr_i(layer-1) + z_sbr_i(layer) )
367
368	dzdsdt(layer) = ( sn_i_b(layer) - s_i_b(ji,layer) )*
369	& deltaz_i_phy(layer) / ddtb
370
371	fsdwn(layer) = 0.25*(
372	& z_sbr_i(layer-1) * ( w_adv_br(layer-1) - w_adv_br(layer) ) +
373	& z_sbr_i(layer) * ( w_adv_br(layer-1) - w_adv_br(layer+1) ) +
374	& z_sbr_i(layer+1) * ( w_adv_br(layer) - w_adv_br(layer+1) ) )
375
376	END DO
377
378	!--- Last layer ---
379	fsupw2(nlay_i) = 0.5s_w( -3.*w_adv_br(nlay_i) +
380	& w_adv_br(nlay_i-1) )
381
382	fsupw1(nlay_i) = -0.25( w_adv_br(nlay_i) + w_adv_br(nlay_i-1) )
383	& ( z_sbr_i(nlay_i) + z_sbr_i(nlay_i-1) )
384
385	dzdsdt(nlay_i) = ( sn_i_b(nlay_i) - s_i_b(ji,nlay_i) ) *
386	& deltaz_i_phy(nlay_i) / ddtb
387
388	fsdwn(nlay_i) = ( w_adv_br(nlay_i-1) - w_adv_br(nlay_i) ) *
389	& ( s_w + 0.5*( z_sbr_i(nlay_i) + z_sbr_i(nlay_i-1) ) )
390
391	!
392	!-----------------------------------------------------------------------
393	! 9) Mass of salt conserved ?
394	!-----------------------------------------------------------------------
395	!
396	! Final mass of salt
397	CALL ice_sal_column( kideb, kiut, z_ms_i_fin ,
398	& sn_i_b(1:nlay_i),
399	& deltaz_i_phy, nlay_i, .FALSE. )
400
401	! Bottom flux (positive upwards for conservation routine)
402	zfb = - SUM(fsdwn(1:nlay_i)) + fsupw2(nlay_i)
403
404	!! ! Bottom flux ( positive upwards for conservation routine )
405	!! zfb = - e_i_b(nlay_i) *
406	!! & ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) *
407	!! & ( z_sbr_i(nlay_i) - s_w ) + w_flood * ( z_flood *
408	!! & s_w + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) )
409	!! & - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i)
410	!! & - qsummer * z_sbr_i(nlay_i) / ddtb
411	!
412	fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards
413	IF ( ln_write ) THEN
414	WRITE(numout,*) ' fsb : ', fsb
415	WRITE(numout,*)
416	ENDIF
417
418	! ! Surface flux of salt
419	zfsu = 0.0
420
421	! ! conservation check
422	zerror = 1.0e-15
423	CALL ice_sal_conserv(kideb,kiut,'ice_sal_adv : ',zerror,
424	& z_ms_i_ini,z_ms_i_fin,
425	& zfb , zfsu , ddtb)
426
427	ENDIF ! ln_sal
428
429
430	!
431	!------------------------------------------------------------------------------\|
432	! End of la sous-routine
433	WRITE(numout,*)
434	END SUBROUTINE

Note: See TracBrowser for help on using the repository browser.

Download in other formats: