Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ice_bio_diff.f @ 67

Last change on this file since 67 was 4, checked in by vancop, 8 years ago
initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 25.0 KB

Line
1	SUBROUTINE ice_bio_diff(kideb,kiut,nlay_i)
2
3	!------------------------------------------------------------------------------!
4	!
5	! --- ice_bio_diff ---
6	!
7	! Transport and diffusion of tracers
8	! (c) Martin Vancoppenolle, May 2007
9	! 1.1 Rayleigh number based diffusivity, Oct 2008
10	! 1.2 Includes open upper boundary condition for gas exchange Jul 2010 - Mar 2011
11	! 1.3 Simplification and implementation of flooding velocity, Feb 2012
12	!
13	! Rewriting of surface gas boundary condition (M. Vancoppenolle & S. Moreau, 2015)
14	! - routine for gas transfer velocity (ice_gas_trvel)
15	! - add solubility and partial pressure for each gas
16	!
17	! ! Use the real surface pCO2 and not the 1st layer value
18	! ! remove ji
19	!
20	!------------------------------------------------------------------------------!
21
22	USE lib_fortran
23
24	INCLUDE 'type.com'
25	INCLUDE 'para.com'
26	INCLUDE 'const.com'
27	INCLUDE 'ice.com'
28	INCLUDE 'thermo.com'
29	INCLUDE 'bio.com'
30
31	INTEGER ::
32	& ji , ! : horizontal space index
33	& jn ! : horizontal space index jn
34
35	REAL(8), DIMENSION( maxnlay ) :: !: dummy factors for tracer equation
36	& za , !: winter
37	& zb ,
38	& zc ,
39	& ze , !: summer
40	& zind , !: independent term in the tridiag system
41	& zindtbis , !:
42	& zdiagbis
43
44	REAL(8), DIMENSION(maxnlay,3) ::!: dummy factors for tracer equation
45	& ztrid !: tridiagonal matrix
46
47	REAL(8), DIMENSION(nlay_bio) ::
48	& z_sbr_i !: brine salinity
49
50	REAL(8), DIMENSION(ntra_bio_max) ::
51	& zpp_gas !: partial pressure
52
53	REAL(8) ::
54	& zdummy1 , !: dummy factors
55	& zdummy2 , !:
56	& zdummy3 , !:
57	& zswitchs , !: switch for summer drainage
58	& f_sn_rat ,
59	& wspd_trs ,
60	& zsat_arg ,
61	& zsat_oxy ,
62	& zsat_CO2 ,
63	& zsat_nit ,
64	& mol_diff
65
66	INTEGER ::
67	& indtr , !: index of tridiagonal system
68	& iter !: time step index
69
70	LOGICAL ::
71	& ln_write_bio ,
72	& ln_con_bio ,
73	& ln_flood
74
75	ln_write_bio = .TRUE.
76	ln_con_bio = .TRUE.
77	ln_flood = .TRUE.
78
79	zsol = 0.
80	zb0 = 0.
81	zpatm_gas = 0.
82	zpp_gas(:) = 0.0
83
84	!=======================================================================
85
86	WRITE(numout,*)
87	WRITE(numout,) ' * ice_bio_diff : '
88	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
89	WRITE(numout,*)
90
91	DO ji = kideb, kiut
92
93	!
94	!-----------------------------------------------------------------------
95	! 1) Initialization
96	!-----------------------------------------------------------------------
97	!
98	IF ( ln_write_bio ) THEN
99	WRITE(numout,*) ' Initialization '
100	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
101	WRITE(numout,*)
102	ENDIF
103
104	!---------------
105	! Interpolation
106	!---------------
107	CALL ice_bio_interp_phy2bio(kideb,kiut,nlay_i,.FALSE.)
108	! interpolation of physical variables
109	! on the biological grid
110	! mass of salt, heat content, brine volume, Rb, PAR
111
112	CALL ice_bio_interp_diffus(kideb,kiut,nlay_i,.TRUE.)
113	IF ( ln_write_bio ) THEN
114	WRITE(numout,*) ' diff_br_bio : ', ( diff_br_bio(layer),
115	& layer = 1, nlay_bio )
116	ENDIF
117
118	!--------------------------------
119	! Brine concentration of tracers
120	!--------------------------------
121	DO jn = 1, ntra_bio
122	IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
123	DO jk = 1, nlay_bio
124	c_i_bio(jn,jk) = cbu_i_bio(jn,jk) / e_i_bio(jk)
125	END DO
126	ENDIF
127	IF ( flag_adsorb(jn) .AND. flag_active(jn) ) THEN
128	DO jk = 1, nlay_bio
129	c_i_bio(jn,jk) = 0.
130	END DO
131	ENDIF
132	WRITE(numout,) ' 01 ** jn = ', jn
133	WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:)
134	END DO
135
136
137	!---------------------------------------------------------------
138	! Equilibrate carbonate system to get aqueous CO2 concentration
139	!---------------------------------------------------------------
140	IF ( ln_carbon ) CALL ice_carb_chem
141
142	DO jn = 1, ntra_bio
143	WRITE(numout,) ' 02 ** jn = ', jn
144	WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:)
145	END DO
146
147	!--------------------
148	! Conservation check
149	!--------------------
150
151	CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_init,cbu_i_bio,
152	& deltaz_i_bio, .FALSE.)
153
154	IF ( ln_write_bio ) THEN
155	DO jn = 1, ntra_bio
156	WRITE(numout,*) ' mt_i_bio_init : ', mt_i_bio_init(jn)
157	END DO
158	ENDIF
159
160	! layer by layer
161	DO jn = 1, ntra_bio
162	DO jk = 1, nlay_bio
163	m_i_bio_init(jn,jk) = cbu_i_bio(jn,jk)*deltaz_i_bio(jk)
164	END DO
165	END DO
166
167	!
168	!-----------------------------------------------------------------------
169	! 3) Surface boundary condition for gases
170	!-----------------------------------------------------------------------
171	!
172
173	zaeff = e_i_bio(1) * brines_ar ! proportionality factor for surface flux
174
175	zpp_gas(:) = mixr_gas(:) * psbqb(ji) / ! atmospheric partial pressure
176	& 101325.
177	!----------------------------
178	CALL ice_gas_solu ! Solubilities (sol_gas(jn))
179	!----------------------------
180
181	!-------------------------
182	gas_trvel(:) = 0. ! Gas transfer velocities
183	!-------------------------
184
185	DO jn = 1, ntra_bio
186
187	IF ( flag_active(jn) .AND. ( biotr_i_typ(jn) .EQ. 'gas' ) .AND. ! select gases only
188	& ( flag_diff(jn) .OR. flag_adsorb(jn) ) .AND. ln_gasflux ) THEN
189
190	SELECT CASE ( i_gasflux )
191
192	CASE (1) ! No flux
193
194	CASE (2) ! Gas diffusion between ice and atm
195
196	IF ( e_i_bio(1) .GT. e_thr_gasflux )
197	& gas_trvel(jn) = dmol_gas(jn) / h_bl_gas
198
199	END SELECT
200
201	IF ( ln_write_bio ) THEN
202
203	WRITE(numout,) ' ** Gas transfer velocity --- '
204	WRITE(numout,*) ' --- Tracer --- : ', biotr_i_nam(jn)
205	WRITE(numout,*) ' gas_trvel : ', gas_trvel(jn)
206
207	ENDIF ! ln_write_bio
208
209	ENDIF ! flags
210
211	!
212	!-----------------------------------------------------------------------
213	! 3) Switches, flushing and flooding velocities
214	!-----------------------------------------------------------------------
215	!
216	!----------
217	! Switches
218	!----------
219	! summer switch
220	zbvmin = 1.0
221	DO layer = 1, nlay_bio
222	zbvmin = MIN( e_i_bio(layer) , zbvmin ) ! minimum brine volume
223	END DO
224	zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 si hiver 1 si ete
225	IF ( zbvmin .LT. e_thr_flu ) zswitchs = 0.0
226
227	IF ( ln_write_bio ) THEN
228	WRITE(numout,*) ' zswitchs : ', zswitchs
229	WRITE(numout,*)
230	ENDIF
231
232	! flood switch
233	IF ( w_flood .GT. 0 ) THEN
234	z_flood = 0.0
235	ELSE
236	z_flood = 1.0
237	ENDIF
238
239	! recompute flushing velocity
240	w_flush = qsummer / ddtb / e_i_bio(1)
241
242	!
243	!-----------------------------------------------------------------------
244	! 4) Compute dummy factors for tracer diffusion equation
245	!-----------------------------------------------------------------------
246	!
247	IF ( ( flag_diff(jn) .OR. flag_adsorb(jn) )
248	& .AND. flag_active(jn) ) THEN
249
250	IF ( ln_write_bio ) THEN
251	WRITE(numout,*) ' --------------------------------- '
252	WRITE(numout,*) ' Diffusion for ', biotr_i_nam(jn)
253	WRITE(numout,*) ' --------------------------------- '
254
255	WRITE(numout,*)
256	WRITE(numout,*) ' Dummy factors '
257	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
258	WRITE(numout,*)
259	ENDIF
260
261	!--------------------
262	! za factors
263	!--------------------
264	DO layer = 1, nlay_bio
265	za(layer) = ddtb / ( deltaz_i_bio(layer) * e_i_bio(layer) )
266	END DO
267
268	!--------------------
269	! zb, zc, ze factors
270	!--------------------
271	DO layer = 1, nlay_bio - 1
272	! interpolate brine volume at the interface between layers
273	zdummy1 = ( e_i_bio(layer + 1) - e_i_bio(layer) ) /
274	& ( z_i_bio(layer + 1) - z_i_bio(layer) )
275	zdummy2 = deltaz_i_bio(layer) / 2.0
276	zdummy3 = e_i_bio(layer) + zdummy1 * zdummy2
277	zb(layer) = zdummy3 * diff_br_bio(layer) /
278	& ( z_i_bio(layer+1) - z_i_bio(layer) )
279	zc(layer) = w_flood * zdummy3 * z_flood
280	ze(layer) = ( w_flood * ( 1. - z_flood ) + w_flush ) * zdummy3
281	END DO
282
283	! Basal values
284	zb(nlay_bio) = 2. * e_i_bio(nlay_bio) /
285	& deltaz_i_bio(nlay_bio) * diff_br_bio(nlay_bio)
286	zc(nlay_bio) = w_flood * e_i_bio(nlay_bio) * z_flood
287	ze(nlay_bio) = ( w_flood * ( 1. - z_flood ) + w_flush ) *
288	& e_i_bio(nlay_bio)
289
290	! Open upper boundary condition for gas
291	zb0 = 0.
292	IF ( biotr_i_typ(jn) .EQ. 'gas' ) zb0 = gas_trvel(jn) * zaeff
293
294	! Block fluxes above the biologically active layer (SL & BAL cases)
295	IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) THEN
296	zb(1:nlay_bio-1) = 0.
297	zc(1:nlay_bio-1) = 0.
298	ze(1:nlay_bio-1) = 0.
299	ENDIF
300
301	IF ( ln_write_bio ) THEN
302	WRITE(numout,*) ' Winter factors '
303	WRITE(numout,*) ' za : ', ( za (layer),
304	& layer = 1, nlay_bio)
305	WRITE(numout,*) ' zb : ', ( zb (layer),
306	& layer = 1, nlay_bio)
307	WRITE(numout,*) ' zc : ', ( zc (layer),
308	& layer = 1, nlay_bio)
309	WRITE(numout,*) ' ze : ', ( ze (layer),
310	& layer = 1, nlay_bio)
311	WRITE(numout,*)
312	ENDIF
313
314	!
315	!-----------------------------------------------------------------------
316	! 5) Tridiagonal system terms for tracer diffusion equation, winter
317	!-----------------------------------------------------------------------
318	!
319	!----------------
320	! first equation
321	!----------------
322	ztrid(1,1) = 0.0
323	ztrid(1,2) = 1.0 + za(1) * ( zb(1) + ze(1) + zb0 )
324	ztrid(1,3) = za(1) * ( -zb(1) + zc(1) )
325	zind(1) = c_i_bio(jn,1) + za(1)zb0sol_gas(jn,1)*zpp_gas(jn)
326
327	! IF ( jn .NE. jn_dic ) THEN
328	! zind(1) = c_i_bio(jn,1) + za(1)zb0sol_gas(jn,1)*zpp_gas(jn)
329	! ELSE
330	! CALL ice_carb_chem
331	! zind(1) = c_i_bio(jn,1) + za(1)zb0( sol_gas(jn,1)*zpp_gas(jn) - co2aq(1) )
332	! ENDIF
333
334	! for dic would read c_i_bio(jn,1) + za(1)zb0(sol_gas(jn,1)*zpp_gas(jn) - zco2 )
335	! where zco2 = co2aq(1)
336
337	!-----------------
338	! inner equations
339	!-----------------
340	DO layer = 2, nlay_bio - 1
341	ztrid(layer,1) = - za(layer) * ( zb(layer-1) + ze(layer-1) )
342	ztrid(layer,2) = 1.0 + za(layer) * ( zb(layer) +
343	& ze(layer) + zb(layer-1) - zc(layer-1) )
344	ztrid(layer,3) = za(layer) * ( -zb(layer) + zc(layer) )
345	zind(layer) = c_i_bio(jn,layer)
346	END DO
347
348	!----------------
349	! last equation
350	!----------------
351	ztrid(nlay_bio,1) = -za(nlay_bio) * ( zb(nlay_bio-1) +
352	& ze(nlay_bio-1) )
353	ztrid(nlay_bio,2) = 1.0 + za(nlay_bio) * ( zb(nlay_bio) +
354	& ze(nlay_bio) + zb(nlay_bio-1) -
355	& zc(nlay_bio-1) )
356	ztrid(nlay_bio,3) = 0.
357	zind(nlay_bio) = c_i_bio(jn,nlay_bio) + za(nlay_bio) *
358	& ( zb(nlay_bio) - zc(nlay_bio) )*c_w_bio(jn)
359
360	IF ( ln_write_bio ) THEN
361	WRITE(numout,*)
362	WRITE(numout,*) ' Tridiag terms : '
363	WRITE(numout,*) ' ~~~~~~~~~~~~~~~ '
364	WRITE(numout,*)
365	DO layer = 1, nlay_bio
366	WRITE(numout,*) ' layer : ', layer
367	WRITE(numout,*) ' ztrid : ', ztrid(layer,1),
368	& ztrid(layer,2),
369	& ztrid(layer,3)
370	WRITE(numout,*) ' zind : ',zind(layer)
371	END DO
372	ENDIF
373	!
374	!-----------------------------------------------------------------------
375	! 6) Solving the tridiagonal system
376	!-----------------------------------------------------------------------
377	!
378	! The tridiagonal system is solved with Gauss elimination
379	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
380	! McGraw-Hill 1984.
381	zindtbis(1) = zind(1)
382	zdiagbis(1) = ztrid(1,2)
383	DO layer = 2, nlay_bio
384	zdiagbis(layer) = ztrid(layer,2) - ztrid(layer,1) *
385	& ztrid(layer-1,3) / zdiagbis(layer-1)
386	zindtbis(layer) = zind(layer) - ztrid(layer,1) *
387	& zindtbis(layer-1) / zdiagbis(layer-1)
388	END DO
389
390	!-----------------------------
391	! Tracer brine concentrations
392	!-----------------------------
393	c_i_bio(jn,nlay_bio) = zindtbis(nlay_bio) / zdiagbis(nlay_bio)
394	DO layer = nlay_bio - 1 , 1 , -1
395	c_i_bio(jn,layer) = (zindtbis(layer) - ztrid(layer,3)*
396	& c_i_bio(jn,layer+1)) / zdiagbis(layer)
397	END DO
398
399	IF ( ln_write_bio ) THEN
400	WRITE(numout,*)
401	WRITE(numout,*) ' Resolution '
402	WRITE(numout,*) ' ~~~~~~~~~~~ '
403	WRITE(numout,*)
404	DO layer = 1, nlay_bio
405	WRITE(numout,*) ' layer : ', layer
406	WRITE(numout,*) ' zdiagbis : ', zdiagbis(layer)
407	WRITE(numout,*) ' zindtbis : ', zindtbis(layer)
408	WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,layer)
409	END DO
410	ENDIF
411
412	ENDIF ! flag
413
414	! Update gas flux
415	fgas(jn) = 0.
416
417	IF ( ( biotr_i_typ(jn) .EQ. 'gas' ) .AND.
418	& flag_active(jn) ) THEN
419
420	fgas(jn) = - zb0 * ( c_i_bio(jn,1) -
421	& sol_gas(jn,1)*zpp_gas(jn) ) ! mmol/m2/s, positive to the ice
422
423	ENDIF
424
425	WRITE(numout,*)
426	WRITE(numout,*) ' tracer : ', biotr_i_nam(jn)
427	WRITE(numout,*) ' fgas : ', fgas(jn)
428	WRITE(numout,*)
429
430
431	END DO ! jn
432
433	!
434	!-----------------------------------------------------------------------
435	! 7) Recover bulk tracer concentrations
436	!-----------------------------------------------------------------------
437	!
438	!--------------------------------
439	! Tracer bulk ice concentrations
440	!--------------------------------
441
442	IF ( ln_write_bio ) THEN
443	WRITE(numout,*)
444	WRITE(numout,*) ' fgas Argon : ', fgas(jn_arg)
445	WRITE(numout,*) ' fgas Oxygen : ', fgas(jn_oxy)
446	WRITE(numout,*) ' fgas CO2 : ', fgas(jn_co2)
447	WRITE(numout,*)
448	ENDIF
449
450	DO jn = 1, ntra_bio
451
452	IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
453
454	DO layer = 1, nlay_bio
455	cbu_i_bio(jn,layer) = c_i_bio(jn,layer) * e_i_bio(layer)
456	END DO
457
458	WRITE(numout,) ' 03 ** jn = ', jn
459	WRITE(numout,*) ' cbu_i_bio : ', cbu_i_bio(jn,:)
460	WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:)
461
462	IF ( jn .EQ. jn_dic ) THEN
463
464	WRITE(numout,*)
465	WRITE(numout,*) ' Change in DIC :'
466	WRITE(numout,*) ' c_i_bio(11,1) before :',c_i_bio(jn,1)
467	WRITE(numout,*) ' cbu_i_bio(11,1) befe :',cbu_i_bio(jn,1)
468
469	cbu_i_bio(jn,1) = cbu_i_bio(jn,1) + fgas(jn_co2) * ! change in bulk DIC due to fgas
470	& ddtb / deltaz_i_bio(1)
471	c_i_bio(jn,1) = cbu_i_bio(jn,1) / e_i_bio(1) ! update surface brine DIC after CO2 flux
472
473	WRITE(numout,*) ' c_i_bio(11,1) after :',c_i_bio(jn,1)
474	WRITE(numout,*) ' cbu_i_bio(11,1) after:',cbu_i_bio(jn,1)
475	WRITE(numout,*) ' fgas(20) :', fgas(jn_co2)
476	WRITE(numout,*)
477
478	ENDIF
479
480	ENDIF ! flag_diff
481
482	IF ( flag_adsorb(jn) .AND. flag_active(jn) ) THEN
483	DO layer = 1, nlay_bio
484	cbu_i_bio(jn,layer) = cbu_i_bio(jn,layer)
485	& + c_i_bio(jn,layer) * e_i_bio(layer)
486	END DO
487	ENDIF ! flag_adsorb
488
489	END DO ! jn
490
491	IF ( ln_write_bio ) THEN
492	DO jn = 1, ntra_bio
493	IF ( ( flag_diff(jn) .OR. flag_adsorb(jn) )
494	& .AND. flag_active(jn) ) THEN
495
496	WRITE(numout,*)
497	WRITE(numout,*) ' Tracer concentrations '
498	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ '
499	WRITE(numout,*)
500	WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn)
501	WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer),
502	& layer = 1, nlay_bio )
503	WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn,layer),
504	& layer = 1, nlay_bio )
505	WRITE(numout,*)
506	WRITE(numout,*) ' The ', biotr_i_nam(jn), 'fluxdwn is :',
507	& fgas(jn), ' mol / m2 / s '
508	ENDIF ! flag_diff
509	END DO ! jn
510	ENDIF ! ln_write_bio
511	!
512	!-----------------------------------------------------------------------
513	! 8) Conservation check
514	!-----------------------------------------------------------------------
515	!
516	IF ( ln_con_bio ) THEN
517
518	IF ( ln_write_bio ) THEN
519	WRITE(numout,*)
520	WRITE(numout,*) ' Conservation check : '
521	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~ '
522	WRITE(numout,*)
523	ENDIF ! ln_write_bio
524
525	CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_final,cbu_i_bio,
526	& deltaz_i_bio, .FALSE.)
527
528	zerror = 1.0d-15
529
530	DO jn = 1, ntra_bio
531
532	IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
533
534	IF ( ln_write_bio ) THEN
535	WRITE(numout,*) ' mt_i_bio_final : ', mt_i_bio_final(jn)
536	ENDIF
537
538	zfb = - e_i_bio(nlay_bio) * ( diff_br_bio(nlay_bio) * 2.0 /
539	& deltaz_i_bio(nlay_bio) * ( c_i_bio(jn,nlay_bio) -
540	& c_w_bio(jn) ) + w_flood * ( z_flood * c_w_bio(jn) +
541	& ( 1. - z_flood ) * c_i_bio(jn,nlay_bio) ) + w_flush *
542	& c_i_bio(jn,nlay_bio) )
543	f_bo_tra(jn) = zfb
544	fcb(jn) = - zfb ! ice-ocean tracer flux
545
546	f_su_tra(jn) = fgas(jn)
547	& + zswitchs * ( qsummer * c_s_bio(jn) ) / ddtb
548
549	! fcb_max = idealized flux if c_i_bio(nlay_bio) .EQ. 0.0 all the time
550	fcb_max(jn)= e_i_bio(nlay_bio) * ( diff_br_bio(nlay_bio) * 2.0 /
551	& deltaz_i_bio(nlay_bio) * ( 0.0 -
552	& c_w_bio(jn) ) + w_flood * ( z_flood * c_w_bio(jn) +
553	& ( 1. - z_flood ) * 0. ) + w_flush * 0. )
554
555
556	IF ( ln_write_bio ) THEN
557	WRITE(numout,*) ' f_bo_tra : ', f_bo_tra(jn)
558	WRITE(numout,*) ' f_su_tra : ', f_su_tra(jn)
559	WRITE(numout,*) ' zswitchs : ', zswitchs
560	WRITE(numout,*) ' qsummer : ', qsummer
561	WRITE(numout,*) ' c_i_bio(nlay_bio) : ',c_i_bio(jn,nlay_bio)
562	WRITE(numout,*) ' mt_i_bio_init : ', mt_i_bio_init(jn)
563	WRITE(numout,*) ' mt_i_bio_final : ', mt_i_bio_final(jn)
564	WRITE(numout,*)
565	ENDIF ! ln_write_bio
566
567	ENDIF ! flag_diff
568
569	END DO ! jn
570
571	CALL ice_bio_conserv(kideb,kiut,ntra_bio,'ice_bio_diff ',
572	& biotr_i_nam, zerror,
573	& mt_i_bio_init,mt_i_bio_final,
574	& f_bo_tra, f_su_tra, ddtb)
575
576	ENDIF ! ln_con_bio
577
578	!
579	!-----------------------------------------------------------------------
580	! X) Control Print
581	!-----------------------------------------------------------------------
582	!
583	WRITE(numout,*)
584	WRITE(numout,) ' * After CO2 fluxes with atmosphere : '
585	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
586	WRITE(numout,*)
587
588	DO jn = 1, ntra_bio
589
590	IF ( biotr_i_nam(jn) .EQ. 'CO2' ) THEN
591	WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn)
592	WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3)
593	& layer = 1, nlay_bio )
594	WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice
595	& layer = 1, nlay_bio )
596	ENDIF
597
598	IF ( biotr_i_nam(jn) .EQ. 'DIC' ) THEN
599	WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn)
600	WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3)
601	& layer = 1, nlay_bio )
602	WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice
603	& layer = 1, nlay_bio )
604	ENDIF
605
606	IF ( biotr_i_nam(jn) .EQ. 'Alk' ) THEN
607	WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn)
608	WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3)
609	& layer = 1, nlay_bio )
610	WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice
611	& layer = 1, nlay_bio )
612	ENDIF
613
614	IF ( biotr_i_nam(jn) .EQ. 'Ika' ) THEN
615	WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn)
616	WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3)
617	& layer = 1, nlay_bio )
618	WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice
619	& layer = 1, nlay_bio )
620	ENDIF
621
622	END DO
623
624	!
625	!-----------------------------------------------------------------------
626	! 8) Flux divergence diagnostics
627	!-----------------------------------------------------------------------
628	!
629	! not sure if those are good anymore
630	WRITE(numout,*)
631	WRITE(numout,*) ' Diagnostics '
632	WRITE(numout,*)
633	WRITE(numout,*)
634
635	! Right now, does not include summer
636	DO jn = 1, ntra_bio
637
638	IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
639
640	fdiff(jn,0) = 0.0 ! top flux
641
642	DO layer = 1, nlay_bio - 1 ! inner fluxes
643	! interpolate brine volume at the interface between layers
644	zdummy1 = ( e_i_bio(layer + 1 ) - e_i_bio(layer) ) /
645	& ( z_i_bio(layer + 1) - z_i_bio(layer) )
646	zdummy2 = deltaz_i_bio(layer) / 2.0
647	zdummy3 = e_i_bio(layer) + zdummy1 * zdummy2
648	fdiff(jn,layer) = - zdummy3 *
649	& diff_br_bio(layer) *
650	& ( c_i_bio(jn,layer+1) - c_i_bio(jn,layer) ) /
651	& ( z_i_bio(layer + 1) - z_i_bio(layer) )
652	END DO
653
654	! lower flux
655	fdiff(jn,nlay_bio) = - 2.0 * e_i_bio(nlay_bio) *
656	& diff_br_bio(layer) *
657	& ( c_w_bio(jn) - c_i_bio(jn,nlay_bio) ) /
658	& deltaz_i_bio(nlay_bio)
659	WRITE(numout,*) ' c_w_bio : ', c_w_bio(jn)
660	WRITE(numout,*) ' c_i_bio(N) : ', c_i_bio(jn,nlay_bio)
661	WRITE(numout,*) ' deltaz_i_bio : ', deltaz_i_bio(nlay_bio)
662
663	! divergence of the fluxes
664	DO layer = 1, nlay_bio
665	diag_divf_bio(jn,layer) =
666	& - ( fdiff(jn,layer) - fdiff(jn,layer-1) ) /
667	& deltaz_i_bio(layer)
668	END DO
669
670	DO layer = 1, nlay_bio - 1
671	WRITE(numout,*) ' fdiff(layer) : ', fdiff(jn,layer)
672	END DO
673	DO layer = 1, nlay_bio - 1
674	WRITE(numout,*) ' div_F(layer) : ', diag_divf_bio(jn,layer)
675	END DO
676
677	ENDIF ! flag_diff
678
679	END DO ! jn
680	!
681	!-----------------------------------------------------------------------
682	! 9) End of the routine
683	!-----------------------------------------------------------------------
684	!
685	END DO ! ji
686
687	IF ( ln_write_bio ) THEN
688
689	WRITE(numout,*)
690	WRITE(numout,) ' After diffusion of tracers * '
691	WRITE(numout,*) ' model output '
692
693	DO jn = 1, ntra_bio
694	IF ( flag_active(jn) ) THEN
695	WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn)
696	WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn, jk),
697	& jk = 1, nlay_bio )
698	ENDIF ! flag_active
699	END DO ! jn
700	WRITE(numout,*)
701
702	ENDIF ! ln_write_bio
703
704	IF ( ln_carbon ) CALL ice_carb_chem
705
706	WRITE(numout,*)
707	WRITE(numout,*) ' End of ice_bio_diff '
708	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
709	!
710	!=============================================================================!
711	!-- End of ice_bio_diff --
712
713	END

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

Context Navigation

source: trunk/SOURCES/source_3.20/ice_bio_diff.f @ 67

Download in other formats: