source: branches/2017/dev_v3.20_2017_transport_max/SOURCES/source_3.20/ice_gas.f

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

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