Context Navigation

source: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_sal_adv.f @ 10

Last change on this file since 10 was 10, checked in by vancop, 8 years ago

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

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

Download in other formats: