source: tags/LIM1D_v3.20/SOURCES/source_3.20/ice_bio_sms.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: 16.0 KB
Line 
1      SUBROUTINE ice_bio_sms(nlay_i,kideb,kiut)
2!
3!-----------------------------------------------------------------------------!
4!                           *** ice_bio_sms ***
5!     Biological sources minus sinks
6!     (c) Martin Vancoppenolle, LOCEAN, October 2012
7!     
8!     version : source_3.03
9!     status  : in permanent revision :-)
10!
11!=============================================================================!
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, DIMENSION(3) :: 
21     &   zindex
22
23      REAL, DIMENSION(nlay_bio) ::
24     &   zek,
25     &   zlim_nut
26
27      LOGICAL ::
28     &   ln_write_bio
29
30      ln_write_bio = .TRUE.
31      zeps = 1.0e-14
32      ji = 1
33
34!
35!=============================================================================!
36!
37!-----------------------------------------------------------------------------!
38! 1) Control prints
39!-----------------------------------------------------------------------------!
40!
41      IF ( ln_write_bio ) THEN
42
43         WRITE(numout,*) 
44         WRITE(numout,*) ' ** ice_bio_sms : '
45         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
46         WRITE(numout,*)
47
48         CALL ice_brine
49
50         WRITE(numout,*) ' Initial tracer values   '
51         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~  '
52
53         DO jn = 1, ntra_bio
54            IF ( flag_active(jn) ) THEN
55               WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
56               WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), 
57     &                                         layer = 1, nlay_bio )
58            ENDIF ! flag_active
59         END DO ! jn
60
61         WRITE(numout,*)
62         WRITE(numout,*) ' Parameters '
63         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~  '
64         WRITE(numout,*) ' si_c      : ', si_c
65         WRITE(numout,*) ' no3_c     : ', no3_c
66         WRITE(numout,*) ' po4_c     : ', po4_c
67         WRITE(numout,*) ' mumax_bio : ', mumax_bio   
68         WRITE(numout,*) ' klys_bio  : ', klys_bio
69         WRITE(numout,*) ' khs_si_bio: ', khs_si_bio
70         WRITE(numout,*) ' ek_bio    : ', ek_bio
71         WRITE(numout,*) ' astar_alg : ', astar_alg
72         WRITE(numout,*) ' rg_bio    : ', rg_bio
73         WRITE(numout,*) ' nn_lim_sal: ', nn_lim_sal 
74         WRITE(numout,*) ' lim_sal_wid :', lim_sal_wid
75         WRITE(numout,*) ' lim_sal_smax:', lim_sal_smax
76         WRITE(numout,*) 
77
78      ENDIF ! ln_write_bio
79!
80!-----------------------------------------------------------------------------!
81! 2) Source Minus Sink terms
82!-----------------------------------------------------------------------------!
83!
84
85      ! Reference of the upper layer
86      IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) THEN
87
88         lim_dsi(:) = 0.; lim_no3(:) = 0.; lim_po4(:) = 0.
89         lim_lig(:) = 0.; lim_tem(:) = 0.; lim_sal(:) = 0.
90         lys_bio(:) = 0.; rem_bio(:) = 0.; syn_bio(:) = 0.
91         rsp_bio(:) = 0.; exu_bio(:) = 0.
92
93         lim_dsi(nlay_bio) = 1.; lim_no3(nlay_bio) = 1.
94         lim_po4(nlay_bio) = 1.
95         lim_lig(nlay_bio) = 1.; lim_tem(nlay_bio) = 1.
96         lim_sal(nlay_bio) = 1.
97
98      ELSE
99
100         lim_dsi(:) = 1.; lim_no3(:) = 1.; lim_po4(:) = 1.
101         lim_lig(:) = 1.; lim_tem(:) = 1.; lim_sal(:) = 1.
102         lys_bio(:) = 0.; rem_bio(:) = 0.; syn_bio(:) = 0.
103         rsp_bio(:) = 0.; exu_bio(:) = 0.
104
105      ENDIF
106
107      !------------------------------------------------------------------------
108      ! 2.1 Limitation factors
109      !------------------------------------------------------------------------
110
111      ! 2.1.1 Nutrients & temperature
112      !--------------------------------
113      DO layer = layer_00, nlay_bio
114
115         IF ( ln_lim_dsi ) ! DSi limitation
116     &      lim_dsi(layer) = c_i_bio(1,layer) 
117     &                     / ( khs_si_bio + c_i_bio(1,layer) )
118
119         IF ( ln_lim_no3 ) ! NO3 limitation
120     &      lim_no3(layer) = c_i_bio(2,layer) 
121     &                     / ( khs_n_bio + c_i_bio(2,layer) )
122
123         IF ( ln_lim_po4 ) ! PO4 limitation
124     &      lim_po4(layer) = c_i_bio(3,layer) 
125     &                     / ( khs_p_bio + c_i_bio(3,layer) )
126
127         IF ( ln_lim_dsi .OR. ln_lim_no3 .OR. ln_lim_po4 ) ! Nut limitation
128!    &      zlim_nut(layer) = lim_dsi(layer) * lim_no3(layer)
129!    &                      * lim_po4(layer)
130     &      zlim_nut(layer) = MIN ( lim_dsi(layer) , lim_no3(layer)
131     &                      , lim_po4(layer) )
132
133         IF ( ln_lim_tem ) ! Temperature factor
134     &      lim_tem(layer) = EXP( rg_bio * tc_bio(layer) )
135
136      END DO ! layer
137
138      ! 2.1.2 Brine salinity
139      !----------------------
140      IF ( ln_lim_sal ) THEN
141
142         SELECT CASE ( nn_lim_sal )
143
144            CASE (1) ! Arrigo and Sullivan (1992)
145
146            DO layer = layer_00, nlay_bio
147               zarg = 2.16 - 8.30e-5 * sbr_bio(layer)**2.11
148     &              - 0.55 * LOG( sbr_bio(layer) )
149               lim_sal(layer) = EXP ( - zarg * zarg )
150            END DO
151
152            CASE (2) ! Simple Gaussian ( Martouf and Grozny, 2010 )
153           
154            DO layer = layer_00, nlay_bio
155                zarg = ( sbr_bio(layer) - lim_sal_smax ) /  lim_sal_wid
156                lim_sal(layer) = EXP( - zarg * zarg )
157            END DO
158
159         END SELECT
160
161      ENDIF ! ln_lim_sal
162
163      ! 2.1.3 Light
164      !-------------
165      IF ( ln_lim_lig ) THEN
166
167         SELECT CASE ( nn_phs )
168
169            CASE (1) ! prescribed Ek
170
171               zek(:) = ek_bio
172               chlC_bio(:) = chla_c
173
174            CASE (2) ! Ek depends on chla/C
175
176               zek(:) = mumax_bio / ( alpha_bio * chla_c )
177               chlC_bio(:) = chla_c
178
179            CASE (3) ! full account for limitation factors
180
181               zek(:) = mumax_bio * lim_tem(:) * lim_sal(:) 
182     &                * zlim_nut(:) / ( alpha_bio * chla_c ) 
183               chlC_bio(:) = chla_c
184
185            CASE (4) ! varying chl-a carbon ratio
186
187               zthmin = 0.01  ! min
188               zthmax = 0.05  ! max
189
190               DO layer = layer_00, nlay_bio
191                  chlC_bio (layer) = ( zthmax - ( zthmax - zthmin ) * 
192     &            MIN( par_bio(layer) / Estar, 1.0 ) ) * zlim_nut(layer)
193               END DO
194
195               zek(:) = mumax_bio * lim_tem(:) * lim_sal(:)
196     &                * zlim_nut(:) / ( alpha_bio * chlC_bio(:) )
197
198         END SELECT
199         
200         DO layer = layer_00, nlay_bio ! Jassby & Platt 76 formulation
201            lim_lig(layer) = TANH ( par_bio(layer) / zek(layer) )
202         END DO
203
204         WRITE(numout,*)
205         WRITE(numout,*) ' lim_lig : ', ( lim_lig(layer), 
206     &                   layer = layer_00, nlay_bio )
207         WRITE(numout,*) ' zlim_nut: ', ( zlim_nut(layer), 
208     &                   layer = layer_00, nlay_bio )
209         WRITE(numout,*) ' par_bio : ', ( par_bio(layer), 
210     &                   layer = layer_00, nlay_bio )
211         WRITE(numout,*) ' zek     : ', ( zek(layer)    , 
212     &                   layer = layer_00, nlay_bio )
213         WRITE(numout,*) ' chlC_bio: ', ( chlC_bio(layer), 
214     &                   layer = layer_00, nlay_bio )
215         WRITE(numout,*)
216
217      ENDIF ! ln_lim_lig
218
219      IF ( .NOT. ln_lim_lig ) chlC_bio(:) = chla_c
220
221      !------------------------------------------------------------------------
222      ! 2.2 Loss & remineralization
223      !------------------------------------------------------------------------
224
225      ! Exudation is zero
226      exu_bio(:) = 0.0
227
228      DO layer = layer_00, nlay_bio
229
230         ! Respiration (NP, NPDr)
231         ! IF ( ln_rsp ) ...
232         rsp_bio(layer) = krsp_bio * cbu_i_bio(jn_aoc,layer) * 
233     &                    lim_tem(layer)
234
235         ! Lysis (NP, NPDr)
236         IF ( ln_lys ) THEN
237
238            IF ( nn_lys .GE. 1 )   ! no T-dependence
239     &         lys_bio(layer) = klys_bio * cbu_i_bio(jn_aoc,layer)
240
241            IF ( nn_lys .EQ. 2 )   ! T-dependence
242     &         lys_bio(layer) = lys_bio(layer) * lim_tem(layer)
243
244         ENDIF
245
246      END DO
247
248      ! Remineralization
249      IF ( ln_rem ) THEN
250
251         SELECT CASE ( nn_bio_opt ) 
252
253         CASE(0) ! NP
254
255            rem_bio(layer_00:nlay_bio) = frem_bio * 
256     &                                   lys_bio(layer_00:nlay_bio)
257
258         CASE(1) ! NPDr
259
260            IF ( nn_rem .GE. 1 )   ! no T-dependence
261     &         rem_bio(layer_00:nlay_bio) = krem_bio * 
262     &                          cbu_i_bio(jn_eoc,layer_00:nlay_bio)
263
264         
265            IF ( nn_rem .EQ. 2 )   ! T-dependence
266     &         rem_bio(layer_00:nlay_bio) = rem_bio(layer_00:nlay_bio) *
267     &                       EXP( rg_bac * tc_bio(layer_00:nlay_bio) )
268
269         END SELECT
270
271      ENDIF
272
273
274      !------------------------------------------------------------------------
275      ! 2.3 Diatom Synthesis
276      !------------------------------------------------------------------------
277
278      DO layer = layer_00, nlay_bio
279   
280         IF ( ln_syn ) THEN
281            ! raw synthesis (in s-1)
282            z_syn = mumax_bio
283     &            * lim_tem(layer) * lim_sal(layer) * lim_lig(layer)
284     &            * zlim_nut(layer)
285
286            ! preclude excessive growth if insufficient nutrient stock
287            zsyn1 = z_syn * cbu_i_bio(jn_aoc,layer) ! normal carbon production rate (mmolC/m3/s)
288            zsyn2 = zsyn1;    zsyn3 = zsyn1;    zsyn4 = zsyn1
289   
290            IF ( ln_lim_dsi )                  ! prevent exhaustion of Si
291     &         zsyn2 = cbu_i_bio(jn_dsi,layer) / ddtb / si_c - zeps 
292 
293            IF ( ln_lim_no3 )                  ! prevent exhaustion of N
294     &         zsyn3 = cbu_i_bio(jn_din,layer) / ddtb / no3_c - zeps
295 
296            IF ( ln_lim_po4 )                  ! prevent exhaustion of P
297     &         zsyn4 = cbu_i_bio(jn_dip,layer) / ddtb / po4_c - zeps
298   
299            ! nutrient stock-limited PP rate
300            syn_bio(layer) = MIN( zsyn1, zsyn2, zsyn3, zsyn4 ) 
301
302         ELSE  ! ln_syn is .FALSE.
303               ! PP is prescribed to zero everywhere
304               ! except in the bottom layer
305               ! where it equals pp_presc
306               ! if sun is shining
307
308            IF ( ( layer .EQ. nlay_bio ) .AND. ( fsolgb(ji) .GT. 5.0 ) )
309     &         syn_bio(layer) = pp_presc 
310
311         ENDIF ! ln_syn
312   
313
314      END DO ! layer
315
316      !------------------------------------------------------------------------
317      ! 2.4 Update reservoirs
318      !------------------------------------------------------------------------
319
320      DO layer = layer_00, nlay_bio
321
322         ! DSi uptake
323         zd_dsi1 = si_c * ( rem_bio(layer) + 
324     &             rsp_bio(layer) - syn_bio(layer) ) * ddtb
325         zd_dsi2 = - cbu_i_bio(jn_dsi,layer)        ! Maximum DSi uptake
326         zd_dsi  = MAX ( zd_dsi1, zd_dsi2 )
327
328         ! NO3 uptake
329         zd_no31 = no3_c * ( rem_bio(layer) +
330     &             rsp_bio(layer) - syn_bio(layer) ) * ddtb
331         zd_no32 = - cbu_i_bio(jn_din,layer)        ! Max NO3 uptake
332         zd_no3  = MAX ( zd_no31, zd_no32 )
333
334         ! update PO4
335         zd_po41 = po4_c * ( rem_bio(layer) +
336     &             rsp_bio(layer) - syn_bio(layer) ) * ddtb
337         zd_po42 = - cbu_i_bio(jn_dip,layer)        ! Max PO4 uptake
338         zd_po4  = MAX ( zd_po41, zd_po42 )
339
340         ! update algae
341         zd_daf1 = ( syn_bio(layer) - lys_bio(layer) -
342     &               rsp_bio(layer) - exu_bio(layer) ) * ddtb
343         zd_daf2 = - cbu_i_bio(jn_aoc,layer)        ! Max algal carbon loss
344         zd_daf = MAX ( zd_daf1, zd_daf2 )
345
346         ! update detritus
347         IF ( nn_bio_opt .GE. 1 ) THEN
348            zd_eoc1 = ( lys_bio(layer) + exu_bio(layer) - 
349     &                  rem_bio(layer) ) * ddtb
350            zd_eoc2 = - cbu_i_bio(jn_eoc,layer)     ! Max detrital carbon loss
351            zd_eoc  = MAX( zd_eoc1, zd_eoc2 )
352         ENDIF
353
354         cbu_i_bio(jn_dsi,layer) = cbu_i_bio(jn_dsi,layer) + zd_dsi
355         cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_no3
356         cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_po4
357         cbu_i_bio(jn_aoc,layer) = cbu_i_bio(jn_aoc,layer) + zd_daf
358
359         IF ( nn_bio_opt .GE. 1 ) 
360     &      cbu_i_bio(jn_eoc,layer) = cbu_i_bio(jn_eoc,layer) + zd_eoc
361
362         IF ( ln_write_bio ) THEN
363
364            WRITE(numout,*) ' Terms for layer: ', layer
365            WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
366            WRITE(numout,*) ' syn_bio   : ', syn_bio(layer)
367            WRITE(numout,*) ' lys_bio   : ', lys_bio(layer)
368            WRITE(numout,*) ' rem_bio   : ', rem_bio(layer)
369            WRITE(numout,*) ' exu_bio   : ', exu_bio(layer)
370            WRITE(numout,*) ' rsp_bio   : ', rsp_bio(layer)
371            WRITE(numout,*) ' lim_dsi   : ', lim_dsi(layer)
372            WRITE(numout,*) ' lim_no3   : ', lim_no3(layer)
373            WRITE(numout,*) ' lim_po4   : ', lim_po4(layer)
374            WRITE(numout,*) ' lim_lig   : ', lim_lig(layer)
375            WRITE(numout,*) ' lim_tem   : ', lim_tem(layer)
376            WRITE(numout,*) ' lim_sal   : ', lim_sal(layer)
377            WRITE(numout,*)
378            WRITE(numout,*) ' tc_bio    : ', tc_bio(layer) 
379            WRITE(numout,*) ' sbr_bio   : ', sbr_bio(layer) 
380            WRITE(numout,*) ' par_bio   : ', par_bio(layer)
381            WRITE(numout,*)
382            WRITE(numout,*) ' zd_dsi    : ', zd_dsi
383            WRITE(numout,*) ' zd_no3    : ', zd_no3
384            WRITE(numout,*) ' zd_po4    : ', zd_po4
385            WRITE(numout,*) ' zd_daf    : ', zd_daf
386            WRITE(numout,*) ' zd_eoc    : ', zd_eoc
387            WRITE(numout,*)
388            WRITE(numout,*) ' dSi       : ', cbu_i_bio(jn_dsi,layer)
389            WRITE(numout,*) ' dIN       : ', cbu_i_bio(jn_din,layer)
390            WRITE(numout,*) ' dIP       : ', cbu_i_bio(jn_dip,layer)
391            WRITE(numout,*) ' AoC       : ', cbu_i_bio(jn_aoc,layer)
392            WRITE(numout,*) ' eoC       : ', cbu_i_bio(jn_eoc,layer)
393            WRITE(numout,*)
394
395         ENDIF ! ln_write_bio
396
397      END DO ! layer
398
399      !------------------------------------------------------------------------
400      ! 2.5 Update DIC, Alk and Oxy
401      !------------------------------------------------------------------------
402
403      DO layer = layer_00, nlay_bio
404
405         zncp  =  - syn_bio(layer) + rsp_bio(layer) + rem_bio(layer)
406
407         IF ( ln_carbon .AND. flag_active(jn_dic) 
408     &                  .AND. flag_active(jn_alk) ) THEN
409
410            zd_dic1   = - zncp * ddtb
411            zd_dic2   = - cbu_i_bio(jn_dic,layer)
412            zd_dic    = MAX( zd_dic1 , zd_dic2 )
413
414            zd_alk1   =   no3_c * zncp * ddtb
415            zd_alk2   = - cbu_i_bio(jn_alk,layer)
416            zd_alk    = MAX( zd_alk1 , zd_alk2 )
417
418            cbu_i_bio(jn_dic,layer) = cbu_i_bio(jn_dic,layer) + zd_dic 
419            cbu_i_bio(jn_alk,layer) = cbu_i_bio(jn_alk,layer) + zd_alk 
420
421            IF ( ln_write_bio ) THEN
422               WRITE(numout,*) ' DIC       : ', cbu_i_bio(jn_dic,layer)
423               WRITE(numout,*) ' Alk       : ', cbu_i_bio(jn_alk,layer)
424            ENDIF
425
426         ENDIF
427
428         IF ( flag_active(jn_oxy) ) THEN
429
430            zd_oxy1   =   oxy_c * zncp * ddtb
431            zd_oxy2   = - cbu_i_bio(jn_oxy,layer)
432            zd_oxy    = MAX( zd_oxy1, zd_oxy2 )
433
434            cbu_i_bio(jn_oxy,layer) = cbu_i_bio(jn_oxy,layer) + zd_oxy
435
436            IF ( ln_write_bio ) THEN
437               WRITE(numout,*) ' Oxy       : ', cbu_i_bio(jn_oxy,layer)
438            ENDIF
439
440         ENDIF
441
442      END DO
443
444!
445!------------------------------------------------------------------------------!
446! X) End of the routine
447!------------------------------------------------------------------------------!
448!
449      WRITE(numout,*)
450      WRITE(numout,*) '    *** After biological sources and sinks *** '
451      WRITE(numout,*) '    model output '
452      DO jn = 1, ntra_bio
453         IF ( flag_active(jn) ) THEN
454           WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn)
455           WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn, jk), jk = 1,
456     &                                      nlay_bio )
457         ENDIF
458      END DO
459      WRITE(numout,*)
460
461      WRITE(numout,*)
462      WRITE(numout,*) ' End of ice_bio_sms '
463      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
464     
465
466!==============================================================================|
467! Fin de la routine ice_bio_sms
468      WRITE(numout,*)
469
470      END
Note: See TracBrowser for help on using the repository browser.