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

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

Code new adsorption (martin)
File size: 10.2 KB
Line 
1      SUBROUTINE ice_bio_adv(kideb,kiut,nlay_i)
2
3!------------------------------------------------------------------------------!
4!
5!                              --- ice_bio_adv ---
6!
7!    Gravity Drainage for tracers (Advection)
8!
9!------------------------------------------------------------------------------!
10
11      USE lib_fortran
12 
13      INCLUDE 'type.com'
14      INCLUDE 'para.com'
15      INCLUDE 'const.com'
16      INCLUDE 'ice.com'
17      INCLUDE 'thermo.com'
18      INCLUDE 'bio.com'
19
20      INTEGER :: 
21     &   ji                 ,    ! : horizontal space index
22     &   jn                      ! : horizontal space index jn
23
24      REAL(8), DIMENSION( maxnlay ) ::  !: dummy factors for tracer equation
25     &   za                 ,    !: winter
26     &   zb                 ,   
27     &   zc                 ,   
28     &   ze                 ,    !: summer
29     &   zind               ,    !: independent term in the tridiag system
30     &   zindtbis           ,    !:
31     &   zdiagbis
32
33      REAL(8), DIMENSION(maxnlay,3) ::!: dummy factors for tracer equation
34     &   ztrid                   !: tridiagonal matrix
35
36      REAL(8), DIMENSION(nlay_bio) :: 
37     &   z_sbr_i                 !: brine salinity
38
39      REAL(8), DIMENSION(ntra_bio_max) :: 
40     &   zpp_gas                 !: partial pressure
41
42      REAL(8) :: 
43     &   zdummy1            ,    !: dummy factors
44     &   zdummy2            ,    !:
45     &   zdummy3            ,    !:
46     &   zswitchs           ,    !: switch for summer drainage
47     &   f_sn_rat           ,
48     &   wspd_trs           ,           
49     &   zsat_arg           ,
50     &   zsat_oxy           ,
51     &   zsat_CO2           ,
52     &   zsat_nit           ,
53     &   mol_diff
54
55      INTEGER :: 
56     &   indtr              ,    !: index of tridiagonal system
57     &   iter                    !: time step  index
58
59      REAL(8) ::
60     &   int_DSI_n            ,
61     &   int_DIN_n            ,
62     &   int_DIP_n            ,
63     &   int_DSI_n1           ,
64     &   int_DIN_n1           ,
65     &   int_DIP_n1
66     
67      LOGICAL ::
68     &   ln_write_bio       ,
69     &   ln_con_bio         ,
70     &   ln_flood           
71
72      ln_write_bio = .TRUE.
73      ln_con_bio   = .TRUE.
74      ln_flood     = .TRUE.
75
76      zsol      = 0.
77      zb0       = 0.
78      zpatm_gas = 0. 
79      zpp_gas(:) = 0.0
80
81!=======================================================================
82
83      WRITE(numout,*) 
84      WRITE(numout,*) ' ** ice_bio_adv : '
85      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
86      WRITE(numout,*)
87
88      DO ji = kideb, kiut
89
90!
91!-----------------------------------------------------------------------
92! 1) Initialization
93!-----------------------------------------------------------------------
94!
95      IF ( ln_write_bio ) THEN
96         WRITE(numout,*) ' Initialization '
97         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
98         WRITE(numout,*)
99      ENDIF
100     
101      !---------------
102      ! Interpolation
103      !---------------
104      CALL ice_bio_interp_phy2bio(kideb,kiut,nlay_i,.FALSE.) 
105                             ! interpolation of physical variables
106                             ! on the biological grid
107                             ! mass of salt, heat content, brine volume, Rb, PAR
108
109      CALL ice_bio_interp_diffus(kideb,kiut,nlay_i,.TRUE.) 
110
111      !--------------------------------
112      ! Brine concentration of tracers
113      !--------------------------------
114      DO jn = 1, ntra_bio
115         ! must implement adsorption here
116         IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
117
118            zf = 1 - f_ads(jn)
119            c_i_bio(jn,:) = cbu_i_bio(jn,:) * zf / e_i_bio(:)
120
121            WRITE(numout,*) ' 01 *** jn = ', jn
122            WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:)
123         ENDIF
124      END DO
125
126
127      !--------------------------------
128      ! Integrated old bulk concentration
129      !--------------------------------
130
131      ! jn = 1 pour => DSI
132      ! jn = 2 pour => DIN
133      ! jn = 3 pour => DIP
134
135      int_DSI_n = SUM (cbu_i_bio(1,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
136      int_DIN_n = SUM (cbu_i_bio(2,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
137      int_DIP_n = SUM (cbu_i_bio(3,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
138
139      WRITE(numout,*) ' int_DSI_n : ', int_DSI_n
140      WRITE(numout,*) ' int_DIN_n : ', int_DIN_n
141      WRITE(numout,*) ' int_DIP_n : ', int_DIP_n
142     
143      !---------------------------------------------------------------
144      ! Equilibrate carbonate system to get aqueous CO2 concentration
145      !---------------------------------------------------------------
146!      IF ( ln_carbon ) CALL ice_carb_chem
147
148!      DO jn = 1, ntra_bio
149!         WRITE(numout,*) ' 02 *** jn = ', jn
150!         WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:)
151!      END DO
152
153      !--------------------
154      ! Conservation check
155      !--------------------
156
157!      CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_init,cbu_i_bio,
158!     &                    deltaz_i_bio, .FALSE.)
159
160!      IF ( ln_write_bio ) THEN
161!         DO jn = 1, ntra_bio
162!            WRITE(numout,*) ' mt_i_bio_init : ', mt_i_bio_init(jn)
163!         END DO
164!      ENDIF
165
166      ! layer by layer
167!      DO jn = 1, ntra_bio
168!         DO jk = 1, nlay_bio
169!            m_i_bio_init(jn,jk) = cbu_i_bio(jn,jk)*deltaz_i_bio(jk)
170!         END DO
171!      END DO
172     
173
174!------------------------------------------------------------------------------|
175! Update ice concentrations
176!------------------------------------------------------------------------------|
177!
178 
179      DO jn = 1, ntra_bio
180         IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN
181
182            DO layer = 1, nlay_bio
183               za(layer) = w_adv_br(layer) * ddtb / deltaz_i_bio(layer)
184            END DO
185
186      ! first layer
187            cbu_i_bio(jn,1) = c_i_bio(jn,1) * ( e_i_b(1) + za(1) ) +
188     &                        c_i_bio(jn,2) * ( - za(1) )
189
190      ! inner layers
191            DO layer = 2, nlay_bio - 1
192               cbu_i_bio(jn,layer) = c_i_bio(jn,layer-1) *
193     &                               (za(layer)/2.) +
194     &                               c_i_bio(jn,layer) * e_i_b(layer)  +
195     &                               c_i_bio(jn,layer+1)*(-za(layer)/2.)
196            END DO
197
198         ! lowermost layer
199            cbu_i_bio(jn,nlay_bio)  = c_i_bio(jn,nlay_bio-1) *
200     &                                ( za(nlay_bio)/2. ) +
201     &                                c_i_bio(jn,nlay_bio)   * 
202     &                                (e_i_b(nlay_bio) +
203     &                                za(nlay_bio)/2. ) - 
204     &                                za(nlay_bio) * c_w_bio(jn) 
205
206      IF ( ln_write_bio ) THEN
207         WRITE(numout,*)
208         WRITE(numout,*) ' cbu_i_bio   : ', ( cbu_i_bio(jn,layer) ,
209     &                   layer = 1, nlay_bio )
210      ENDIF
211
212         ENDIF ! flag_diif flag_active
213   
214      END DO  ! ntrabio
215
216
217!--------------------------------
218! Integrated NEW bulk concentration
219!--------------------------------
220
221      ! jn = 1 pour => DSI
222      ! jn = 2 pour => DIN
223      ! jn = 3 pour => DIP
224     
225      int_DSI_n1 = SUM(cbu_i_bio(1,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
226      int_DIN_n1 = SUM(cbu_i_bio(2,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
227      int_DIP_n1 = SUM(cbu_i_bio(3,1:nlay_bio)*deltaz_i_bio(1:nlay_bio))
228
229      WRITE(numout,*) ' int_DSI_n1 : ', int_DSI_n1
230      WRITE(numout,*) ' int_DIN_n1 : ', int_DIN_n1
231      WRITE(numout,*) ' int_DIP_n1 : ', int_DIP_n1
232
233!-------------------------------------------------
234! Flux nutrient with Sea Water (positive donward)
235!-------------------------------------------------
236
237      FDSI_AD = - (int_DSI_n1 - int_DSI_n)/ddtb
238      FDIN_AD = - (int_DIN_n1 - int_DIN_n)/ddtb
239      FDIP_AD = - (int_DIP_n1 - int_DIP_n)/ddtb
240
241      WRITE(numout,*) ' FDIN_AD : ', FDIN_AD
242      WRITE(numout,*) ' FDIP_AD : ', FDIP_AD
243      WRITE(numout,*) ' FDSI_AD : ', FDSI_AD
244     
245      WRITE(numout,*) 'TESTBIO  w_adv_br : ', ( w_adv_br(layer),
246     &                   layer = 1, nlay_bio )
247
248!
249!-----------------------------------------------------------------------
250! 8) Mass of salt conserved ?
251!-----------------------------------------------------------------------
252!
253      ! Final mass of salt
254!      CALL ice_sal_column( kideb , kiut , z_ms_i_fin ,
255!     &                    sn_i_b(1:nlay_i),
256!     &                    deltaz_i_phy, nlay_i, .FALSE. )
257
258      ! Bottom flux ( positive upwards for conservation routine )
259!      zfb      =  - e_i_b(nlay_bio) *
260!     &           ( diff_br(nlay_bio) * 2.0 / deltaz_i_bio(nlay_bio) *
261!     &           ( c_i_bio(nlay_bio) - oce_sal ) + w_flood * ( z_flood *
262!     &           oce_sal + ( 1. - z_flood ) * c_i_bio(nlay_bio) ) )
263!     &           - e_i_b(nlay_bio) * w_flush * c_i_bio(nlay_bio)
264!     &            - qsummer * z_sbr_i(nlay_i) / ddtb
265
266
267!      fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards
268!      IF ( ln_write ) THEN
269!         WRITE(numout,*) ' fsb : ', fsb
270!         WRITE(numout,*)
271!      ENDIF
272
273      ! Surface flux of salt
274!      zfsu     = 0.0
275
276      ! conservation check
277!      zerror = 1.0e-15
278!      CALL ice_sal_conserv(kideb,kiut,'ice_sal_adv  : ',zerror,
279!     &                           z_ms_i_ini,z_ms_i_fin,
280!     &                           zfb   , zfsu  , ddtb)
281
282!      ENDIF ! ln_sal
283!
284!------------------------------------------------------------------------------|
285
286
287!
288!-----------------------------------------------------------------------
289! 9) End of the routine
290!-----------------------------------------------------------------------
291!
292      END DO ! ji
293
294      IF ( ln_write_bio ) THEN
295
296         WRITE(numout,*)
297         WRITE(numout,*) ' *** After advection of tracers *** '
298         WRITE(numout,*) '     model output '
299
300         DO jn = 1, ntra_bio
301            IF ( flag_active(jn) ) THEN
302               WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn)
303               WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn, jk), 
304     &                         jk = 1, nlay_bio )
305            ENDIF ! flag_active
306         END DO ! jn
307         WRITE(numout,*)
308
309      ENDIF ! ln_write_bio
310
311      IF ( ln_carbon ) CALL ice_carb_chem
312
313      WRITE(numout,*)
314      WRITE(numout,*) ' End of ice_bio_adv '
315      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
316!
317!=============================================================================!
318!-- End of ice_bio_adv --
319 
320      END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.