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

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

Ludivine source files
File size: 11.9 KB
Line 
1      SUBROUTINE ice_gas(nlay_i,kideb,kiut)
2
3!
4!-----------------------------------------------------------------------------!
5!                           *** ice_gas ***
6!     
7!     Nov 2015
8!
9!     Version : source_3.10
10!
11!     to do
12!     - put the computation of solubilities in a subroutine
13!        - check conservation
14!        - brine characteristics to subroutinize   
15!        - missing carbonate chemistry in the end
16!-----------------------------------------------------------------------------!
17
18      USE lib_fortran
19 
20      INCLUDE 'type.com'
21      INCLUDE 'para.com'
22      INCLUDE 'const.com'
23      INCLUDE 'ice.com'
24      INCLUDE 'thermo.com'
25      INCLUDE 'bio.com'
26
27       REAL(8) :: 
28     &   zsat,
29     &   zsat_diff
30
31       LOGICAL ::
32     &   ln_write_gas = .TRUE.
33
34!==============================================================================!
35
36      IF ( ln_write_gas ) THEN
37
38         WRITE(numout,*) 
39         WRITE(numout,*) ' ** ice_gas : '
40         WRITE(numout,*) ' ~~~~~~~~~~~~ '
41   
42      ENDIF
43     
44!
45!------------------------------------------------------------------------------!
46! 1) Initialize conservation test and gas concentrations
47!------------------------------------------------------------------------------!
48!
49
50
51      !--- Gas concentrations
52      DO jn = 1, ntra_bio
53
54         ! Brine concentrations
55         IF ( ( flag_active(jn)            ) .AND.
56     &        ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
57            DO layer = 1, nlay_bio
58              c_i_bio(jn,layer)   = cbu_i_bio(jn,layer) / e_i_bio(layer)
59            END DO
60         ENDIF
61
62      END DO
63
64      ! Carbonate chemistry
65      IF ( ln_carbon ) CALL ice_carb_chem
66
67      DO jn = 1, ntra_bio
68
69         IF ( ( flag_active(jn)            ) .AND.
70     &        ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
71            DO layer = 1, nlay_bio
72              c_gtot_i(jn,layer)  = cbu_i_bio(jn,layer) + 
73     &                              cbub_i_bio(jn,layer)
74            END DO
75
76         ENDIF
77
78      END DO
79
80      !--- Conservation test
81      CALL ice_gas_column( ntra_bio , c_gtot_i, deltaz_i_bio, .TRUE., 
82     &                     mt_i_gas_init )
83
84      DO jn = 1, ntra_bio
85
86         IF ( ( flag_active(jn)            ) .AND.
87     &        ( biotr_i_typ(jn) .EQ. 'gas' ) .AND. 
88     &        ln_write_gas                   ) THEN
89
90           WRITE(numout,*)
91           WRITE(numout,*) ' *** Initial values  *** '
92           WRITE(numout,*) ' --- time step --- : ', numit
93
94           WRITE(numout,*) ' --- Tracer    --- :   ', biotr_i_nam(jn)
95           WRITE(numout,*) 'c_i_bio  : ',c_i_bio(jn,layer_00:nlay_bio)
96           WRITE(numout,*) 'cbu_i_bio: ',cbu_i_bio(jn,layer_00:nlay_bio)
97           WRITE(numout,*) 'cbub_i_bio: ',
98     &                     cbub_i_bio(jn,layer_00:nlay_bio)
99           WRITE(numout,*) 'c_gtot_i : ',c_gtot_i(jn,layer_00:nlay_bio)
100           WRITE(numout,*)
101           WRITE(numout,*) ' mt_i_gas_init : ', mt_i_gas_init(jn)
102
103         ENDIF
104
105      END DO
106     
107!
108!------------------------------------------------------------------------------!
109! 2) Bubble formation and dissolution
110!------------------------------------------------------------------------------!
111!
112
113                          !---------------------------
114      CALL ice_gas_solu   ! Saturation concentrations
115                          !---------------------------
116                 
117      IF ( ln_bubform ) THEN
118     
119         IF ( ln_write_gas ) THEN
120            WRITE(numout,*)
121            WRITE(numout,*) ' *** Bubble formation and dissolution *** '
122            WRITE(numout,*)
123         ENDIF
124         
125          DO jn = 1, ntra_bio
126         
127             IF ( flag_active(jn) .AND. 
128     &         ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
129     
130                IF ( ln_write_gas ) THEN
131                   WRITE(numout,*) ' Tracer number : ', jn
132                   WRITE(numout,*) ' --- ', biotr_i_nam(jn)
133                ENDIF
134     
135                DO layer = 1, nlay_bio
136               
137                   zsat      = csat_gas(jn,layer) * sursat_gas
138                   zsat_diff = zsat - c_i_bio(jn,layer)    ! undersaturation
139                   
140                   ! bubble formation or dissolution rate
141                   IF ( zsat_diff .LE. 0 ) THEN  ! sursaturation
142                      zbub_rate = bub_form_rate
143                   ELSE
144                      zbub_rate = bub_diss_rate
145                   ENDIF
146                   
147                   ! change in bulk bubble concentration
148                   zd_cbri = zsat_diff * zbub_rate * ddtb
149                   zd_cbri = MAX( -c_i_bio(jn,layer), zd_cbri )         ! cannot exhaust brines
150                   zd_cbri = MIN( cbub_i_bio(jn,layer)/e_i_bio(layer),  ! cannot exhaust gas bubbles
151     &                            zd_cbri )
152                   c_i_bio(jn,layer) = c_i_bio(jn,layer) + zd_cbri
153               
154                   ! change in bulk bubble concentration
155                   zd_cbub = - zd_cbri*e_i_bio(layer)
156                   cbub_i_bio(jn,layer) = cbub_i_bio(jn,layer) + zd_cbub
157 
158                   ! Update bulk concentration
159                   cbu_i_bio(jn,layer) =e_i_bio(layer)*c_i_bio(jn,layer)
160 
161                   ! Diagnostic % saturation ratio
162                   pc_sat(jn,layer) = c_i_bio(jn,layer) / 
163     &                                csat_gas(jn,layer) * 100.
164 
165                   ! Update DIC for CO2
166                   IF ( jn .EQ. jn_co2 ) 
167     &               cbu_i_bio(jn_dic,layer) = 
168     &                  cbu_i_bio(jn_dic,layer) - zd_cbub
169             
170                END DO ! layer
171         
172             ENDIF ! flags, biotr_i_typ
173             
174          END DO ! ntra_bio
175!
176       ENDIF ! ln_bubform
177
178!
179!------------------------------------------------------------------------------!
180! 3) Bubble rise
181!------------------------------------------------------------------------------!
182!
183      IF ( ln_bubrise ) THEN
184     
185         IF ( ln_write_gas ) THEN
186            WRITE(numout,*)
187            WRITE(numout,*) ' *** Bubble rise *** '
188         ENDIF
189         
190         f_bub(:) = 0.
191         
192         IF ( ln_write_gas ) WRITE(numout,*) ' *** immediate rise'
193         
194             DO jn = 1, ntra_bio
195         
196                IF ( flag_active(jn) .AND. 
197     &             ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
198     
199                   IF ( ln_write_gas ) THEN
200                      WRITE(numout,*) ' Tracer number : ', jn
201                      WRITE(numout,*) ' --- ', biotr_i_nam(jn)
202                   ENDIF
203                   
204                   DO layer = nlay_bio, 1, -1
205                   
206                   zswitch = 0. ! by default, bubbles rise
207                   
208                   ! i_imper_layer: identifies the first impermeable layer above "layer"
209                   ! if no such layer, then i_imper_layer = 0 and gas get out of the ice
210                   ! via f_bub
211                   
212                   i_imper_layer = 0
213                   
214                   IF ( e_i_bio(layer) .LE. e_thr_bubrise ) THEN
215                      i_imper_layer = layer
216                   ELSE
217                      DO layer_2 = layer - 1, 1, -1 ! find a layer above
218                       IF ( ( zswitch .LE. 0 ) .AND. 
219     &                    ( e_i_bio(layer_2) .LT. e_thr_bubrise ) ) THEN
220                            zswitch = 1.
221                            i_imper_layer = MAX(i_imper_layer, layer_2)
222                       ENDIF ! e_i_bio
223                      END DO ! layer2
224                   ENDIF
225
226                   ! Transfer gas among layers
227                   IF ( ( i_imper_layer .GT. 0 ) .AND. 
228     &                  ( i_imper_layer .NE. layer ) ) THEN ! in-ice bubble transfer
229                      cbub_i_bio(jn,i_imper_layer) = ! add gas in the impermeable layer
230     &                cbub_i_bio(jn,i_imper_layer)+cbub_i_bio(jn,layer)
231                      cbub_i_bio(jn,layer) = 0.      ! remove gas from the permeable layer
232                   ENDIF
233                   
234                   IF ( i_imper_layer .EQ. 0 ) THEN ! ice-atmosphere bubble transfer
235                      f_bub(jn) = f_bub(jn) + cbub_i_bio(jn,layer) *
236     &                deltaz_i_bio(layer) / ddtb
237                      cbub_i_bio(jn,layer) = 0.      ! remove gas from the permeable layer
238                   ENDIF ! i_imper_layer
239                 
240                 END DO
241               
242                ENDIF ! flag_active, biotr_i_typ
243               
244             END DO ! jn
245         
246      ENDIF ! ln_bubrise
247
248!
249!------------------------------------------------------------------------------!
250! 4) Update total gas concentration
251!------------------------------------------------------------------------------!
252!
253
254      DO jn = 1, ntra_bio
255         
256         IF ( ( flag_active(jn)            ) .AND. 
257     &        ( biotr_i_typ(jn) .EQ. 'gas' ) .AND.
258     &        ln_write_gas ) THEN
259
260         DO layer = 1, nlay_bio  ! dégazage des couches 2 à nlay_bio
261            c_gtot_i(jn,layer) = cbu_i_bio(jn,layer) + 
262     &                           cbub_i_bio(jn,layer) 
263         END DO
264
265         WRITE(numout,*) '    *** After ice_gas *** '
266
267         WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
268         WRITE(numout,*) ' c_i_bio  : ', ( c_i_bio(jn,layer),    ! concentration of tracers in brines (mmol m-3)
269     &                                   layer = 1, nlay_bio )
270         WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer),  ! concentration of tracers in bulk ice
271     &                                   layer = 1, nlay_bio )
272         WRITE(numout,*) ' cbub_i_bio : ', ( cbub_i_bio(jn,layer), ! concentration of tracers in bubbles in bulk ice
273     &                                   layer = 1, nlay_bio )
274         WRITE(numout,*) ' c_gtot_i : ', ( c_gtot_i(jn,layer),   ! concentration of tracers in bubbles + bulk ice
275     &                                   layer = 1, nlay_bio )
276         WRITE(numout,*)
277         
278         ENDIF
279         
280      ENDDO
281
282!
283!------------------------------------------------------------------------------!
284! 4) Conservation test and output
285!------------------------------------------------------------------------------!
286!
287
288      CALL ice_gas_column( ntra_bio , c_gtot_i, deltaz_i_bio, .TRUE., 
289     &                     mt_i_gas_final )
290
291      IF ( ln_write_gas ) THEN
292         DO jn = 1, ntra_bio
293            WRITE(numout,*) ' tracer         : ', biotr_i_nam(jn)
294            WRITE(numout,*) ' mt_i_gas_init  : ', mt_i_gas_init(jn)
295            WRITE(numout,*) ' mt_i_gas_final : ', mt_i_gas_final(jn)
296         END DO
297      ENDIF
298
299      DO jn = 1, ntra_bio
300         zfbub = f_bub(jn) * ddtb
301         WRITE(numout,*) ' tracer : ', biotr_i_nam(jn)
302         WRITE(numout,*) ' zfbub  : ', zfbub
303
304      IF ( ABS( mt_i_gas_init(jn) - mt_i_gas_final(jn) - zfbub ) 
305     &                                            .GT. 1.0e-10 )
306     &   THEN
307         WRITE(numout,*)
308         WRITE(numout,*) ' ALERT!!! TOTAL GAS NOT CONSERVED '
309         WRITE(numout,*) '      at time step : ', numit
310         WRITE(numout,*) '         tracer    : ', jn
311         WRITE(numout,*) ' mt_i_gas_init     : ', mt_i_gas_init(jn)
312         WRITE(numout,*) ' mt_i_gas_final    : ', mt_i_gas_final(jn)
313         WRITE(numout,*) ' diff    : ', mt_i_gas_init(jn) - 
314     &                                  mt_i_gas_final(jn) - zfbub
315         WRITE(numout,*) ' zfbub   : ', zfbub
316         WRITE(numout,*)
317      ENDIF
318      ENDDO
319     
320!
321!------------------------------------------------------------------------------!
322! x) Update DIC and carbonate chemistry
323!------------------------------------------------------------------------------!
324!
325      IF ( ln_carbon ) CALL ice_carb_chem
326
327!------------------------------------------------------------------------------!
328! X) End of the routine
329!------------------------------------------------------------------------------!
330!
331      WRITE(numout,*)
332      WRITE(numout,*) ' End of ice_gas '
333      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
334     
335
336!==============================================================================|
337! Fin de la routine ice_gas
338      WRITE(numout,*)
339
340      END
Note: See TracBrowser for help on using the repository browser.