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

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

Code new adsorption (martin)

File size: 15.9 KB
RevLine 
[4]1      SUBROUTINE ice_bio_ini(kideb,kiut,nlay_i)
2!
3!------------------------------------------------------------------------------!
4!
5! This routine initializes biogeochemical tracers
6! (c) Martin Vancoppenolle, Oct 2014
7!
8! version: 3.04
9!
10!------------------------------------------------------------------------------!
11!
12      INCLUDE 'type.com'
13      INCLUDE 'para.com'
14      INCLUDE 'const.com'
15      INCLUDE 'ice.com'
16      INCLUDE 'thermo.com'
17      INCLUDE 'bio.com'
18
19      INTEGER :: 
20     &  ji          , ! : index for space
21     &  jk          , ! : index for ice layers
22     &  jn          , ! : index for tracers
23     &  numtra = 700, ! : reference number for bio.param
24     &  n_raw         ! : number of values in the raw initial profile
25
26      REAL(4)        zini(1)  ! forcing field dummy array
27
28      REAL(8), DIMENSION(maxnlay) ::
29     &   z_raw,
30     &   zq_raw,
31     &   zc_raw,
32     &   zq1,
33     &   zdh_bio
34
35      REAL(8), DIMENSION(maxnlay+2) ::
36     &   zdh_raw 
37
38      REAL(8), DIMENSION(0:maxnlay) ::
39     &   zb_raw
40
41      CHARACTER(len=61) :: 
42     &   zblank
43
44      CHARACTER(len=10) :: 
45     &   filenc='init.nc'
46
47      CHARACTER(len=3) :: 
48     &   zchar
49
50      LOGICAL ::
51     &   ln_initfile = .FALSE.
52
53      zblank = 
54     & '                                                             '
55      ji = 1
56
57!------------------------------------------------------------------------------!
58
59      WRITE(numout,*) ' ** ice_bio_ini : '
60      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
61      WRITE(numout,*) 
62
63      IF ( c_bio_model .EQ. 'KRILL' ) THEN
[20]64         ntra_bio = 23
[4]65      ENDIF
66!
67!-----------------------------------------------------------------------------!
68! 1) Create grid and interpolate physical variables
69!-----------------------------------------------------------------------------!
70!
71
72      CALL ice_bio_grid(kideb,kiut,nlay_i,.TRUE.)             ! Biological grid
73
74      CALL ice_bio_interp_phy2bio(kideb,kiut,nlay_i, .FALSE.) ! Physical fields
75                                                              ! on the bio grid
76
77!
78!-----------------------------------------------------------------------------!
79! 2) Initialize tracer numbers and names
80!-----------------------------------------------------------------------------!
81!
82      WRITE(numout,*) ' Initialize tracers... '
83      WRITE(numout,*)
84
85      ! Set default tracer values
86      layer_00 = 1
87      cbu_i_bio(:,:)  = 0.
88      cbub_i_bio(:,:) = 0. ;  c_gtot_i(:,:)   = 0.
89      c_w_bio(:)      = 0. ;  mixr_gas(:)     = 0.
[35]90      dmol_gas(:)     = 0. ;  f_ads(:)        = 0.
[4]91
92      flag_active(:)    = .false.    ! Activated ?
93      flag_diff(:)      = .false.    ! Diffused ?
94      nn_remp(:)        = 1          ! Squeezing remaping
95      biotr_i_nam(:)    = 'xxx'//zblank      ! Name of the tracer
96      biotr_i_typ(:)    = 'xxx'//zblank      ! Type = algal
97      biotr_i_uni(:)    = 'xxx'//zblank ! Units
98      nn_init(:) = 1
99
100      ! Compute layer limits
101      IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) 
102     &   layer_00 = nlay_bio
103
104      ! Assign Tracer Numbers
105      jn_dsi = 1 
106      jn_din = 2
107      jn_dip = 3 
108      jn_aoc = 4 
109      jn_eoc = 5 
110      jn_co2 = 6  ! 20
111      jn_dic = 7  ! 12
112      jn_alk = 8  ! 13
113      jn_ika = 9  ! 14
114      jn_oxy = 10 ! 15
115      jn_nit = 11 ! 16
116      jn_arg = 12 ! 17
117      jn_cal = 13 ! 21
[20]118      jn_aon = 14
119      jn_eon = 15
120      jn_aop = 16
121      jn_eop = 17
[4]122
123      ! Assign Tracer names
124      biotr_i_nam(jn_dsi) = 'dSi'//zblank
125      biotr_i_nam(jn_din) = 'dIN'//zblank
126      biotr_i_nam(jn_dip) = 'dIP'//zblank
127      biotr_i_nam(jn_aoc) = 'AoC'//zblank
128      biotr_i_nam(jn_eoc) = 'eoC'//zblank
129      biotr_i_nam(jn_co2) = 'CO2'//zblank
130      biotr_i_nam(jn_dic) = 'DIC'//zblank 
131      biotr_i_nam(jn_alk) = 'Alk'//zblank 
132      biotr_i_nam(jn_ika) = 'Ika'//zblank
133      biotr_i_nam(jn_cal) = 'Cal'//zblank
134      biotr_i_nam(jn_arg) = 'Arg'//zblank
135      biotr_i_nam(jn_oxy) = 'Oxy'//zblank
136      biotr_i_nam(jn_nit) = 'Nit'//zblank
[20]137      biotr_i_nam(jn_aon) = 'AoN'//zblank
138      biotr_i_nam(jn_eon) = 'eoN'//zblank
139      biotr_i_nam(jn_aop) = 'AoP'//zblank
140      biotr_i_nam(jn_eop) = 'eoP'//zblank
[4]141
142!
143!-----------------------------------------------------------------------------!
144! 3) Read tracer parameters
145!-----------------------------------------------------------------------------!
146!
147
148      WRITE(numout,*) ' Tracer parameters ... '
149      WRITE(numout,*)
150
151      OPEN( unit = numtra , file='tracer.param', status='old' )
152     
153      ! initial values
154      READ(numtra,*)
155      READ(numtra,*)
156      READ(numtra,*)
157      READ(numtra,*)
158      READ(numtra,*)
159      READ(numtra,*)
160      READ(numtra,*)
161      READ(numtra,*)
162      READ(numtra,*)
163      READ(numtra,*)
164      READ(numtra,*)
165      READ(numtra,*)
166      READ(numtra,*)
167      READ(numtra,*)
168      READ(numtra,*)
169      READ(numtra,*)
170      READ(numtra,*)
171      READ(numtra,*)
172      READ(numtra,*)
173      READ(numtra,*)
174      READ(numtra,*)
[35]175      READ(numtra,*)
[4]176
[20]177      DO jn_read = 1, 17
[4]178
179         READ(numtra,*)  c_read_name(jn_read), zchar, i_dummy1, zdummy2,
[35]180     &                   zdummy3, zdummy4, zdummy5, zdummy6
[4]181
182         DO jn = 1, ntra_bio
183           IF ( c_read_name(jn_read)//zblank .EQ. biotr_i_nam(jn) ) THEN
184               c_nc_name(jn) = zchar
185               nn_init(jn)   = i_dummy1
186               cbu_i_nml(jn) = zdummy2
187               c_w_bio(jn)   = zdummy3
188               mixr_gas(jn)  = zdummy4
189               dmol_gas(jn)  = zdummy5
[35]190               f_ads(jn)     = zdummy6
[4]191           ENDIF
192         END DO
193
194      ENDDO
195
196      CLOSE( numtra )
197
198      WRITE(numout,*)
199      DO jn = 1, ntra_bio
200         WRITE(numout,*) biotr_i_nam(jn), nn_init(jn) , cbu_i_nml(jn), 
[35]201     &                   c_w_bio(jn)    , mixr_gas(jn), dmol_gas(jn),
202     &                   f_ads(jn)
[4]203      END DO
204
205!
206!-----------------------------------------------------------------------------!
207! 4) Parameterize tracer behaviour
208!-----------------------------------------------------------------------------!
209!
210      ! flag_active determines a tracer is activated in the code
211      ! flag_diff tells whether a tracer should be diffused with brine or not
212      ! nn_remp gives the type of remapping - 1 = linear remapping (very diffusive); 2 = "squeezing remapping"
213      ! biotr_i_typ can be either - only gas is useful because flag_diff makes the difference between solutes and particles
214      !    - 'sol' -> solute (like salt)
215      !    - 'prc' -> particulate (e.g. organic matter)
216      !    - 'gas' -> gas
217
218      !------------------------
219      ! DSi : Dissolved silica
220      !------------------------
221      jn = jn_dsi
222      flag_active(jn) = .true.
223      flag_diff(jn)   = .true.
224      nn_remp(jn)     = 1
225      biotr_i_typ(jn) = 'sol'//zblank
226      biotr_i_uni(jn) = 'mmol_m3'//zblank
227
228      !------------------------------------
229      ! DIN : dissolved inorganic nitrogen
230      !------------------------------------
231      jn = jn_din
232      flag_active(jn) = .true.
233      flag_diff(jn)   = .true.
234      nn_remp(jn)     = 1
235      biotr_i_typ(jn) = 'sol'//zblank
236      biotr_i_uni(jn) = 'mmol_m3'//zblank
237
238      !----------------------------
239      ! DIP : Dissolved phosphorus
240      !----------------------------
241      jn = jn_dip
242      ! Units = ( mmole PO4 m-3 ) or ( micrmol PO4 l-1 )
243      flag_active(jn) = .true.
244      flag_diff(jn)   = .true.
245      nn_remp(jn)     = 1
246      biotr_i_typ(jn) = 'sol'//zblank
247      biotr_i_uni(jn) = 'mmol_m3'//zblank
248
249      !--------------------------------------
250      ! AoC : Ice algae organic carbon
251      !--------------------------------------
252      jn = jn_aoc
253      flag_active(jn) = .true.
254      flag_diff(jn)   = .false.
[20]255      nn_remp(jn)     = 2
[4]256      biotr_i_typ(jn) = 'prc'//zblank
257      biotr_i_uni(jn) = 'mmol_m3'//zblank
258
259      !--------------------------------------
260      ! eoC : Detritus organic Carbon
261      !--------------------------------------
262      jn = jn_eoc
263      flag_active(jn) = .true.
264      flag_diff(jn)   = .false.
[20]265      nn_remp(jn)     = 2
[4]266      biotr_i_typ(jn) = 'prc'//zblank
267      biotr_i_uni(jn) = 'mmol_m3'//zblank
268
269      !--------------------------------------
[20]270      ! AoN : Ice algae organic N
271      !--------------------------------------
272      jn = jn_aon
273      flag_active(jn) = .true.
274      flag_diff(jn)   = .false.
275      nn_remp(jn)     = 2
276      biotr_i_typ(jn) = 'prc'//zblank
277      biotr_i_uni(jn) = 'mmol_m3'//zblank
278
279      !--------------------------------------
280      ! eoN : Detritus organic N
281      !--------------------------------------
282      jn = jn_eon
283      flag_active(jn) = .true.
284      flag_diff(jn)   = .false.
285      nn_remp(jn)     = 2
286      biotr_i_typ(jn) = 'prc'//zblank
287      biotr_i_uni(jn) = 'mmol_m3'//zblank
288
289      !--------------------------------------
290      ! AoP : Ice algae organic P
291      !--------------------------------------
292      jn = jn_aop
293      flag_active(jn) = .true.
294      flag_diff(jn)   = .false.
295      nn_remp(jn)     = 2
296      biotr_i_typ(jn) = 'prc'//zblank
297      biotr_i_uni(jn) = 'mmol_m3'//zblank
298
299      !--------------------------------------
300      ! eoP : Detritus organic P
301      !--------------------------------------
302      jn = jn_eop
303      flag_active(jn) = .true.
304      flag_diff(jn)   = .false.
305      nn_remp(jn)     = 2
306      biotr_i_typ(jn) = 'prc'//zblank
307      biotr_i_uni(jn) = 'mmol_m3'//zblank
308
309      !--------------------------------------
[4]310      ! DIC : Dissolved Inorganic Carbon
311      !--------------------------------------
312      jn = jn_dic
313      flag_active(jn) = .true.
314      flag_diff(jn)   = .true.
315      nn_remp(jn)     = 1
316      biotr_i_typ(jn) = 'sol'//zblank
317      biotr_i_uni(jn) = 'mmol_m3'//zblank
318
319      !--------------------------------------
320      ! Alk : Total Alkalinity
321      !--------------------------------------
322      jn = jn_alk
323      flag_active(jn) = .true.
324      flag_diff(jn)   = .true.
325      nn_remp(jn)     = 1
326      biotr_i_typ(jn) = 'sol'//zblank
327      biotr_i_uni(jn) = 'mmol_m3'//zblank
328
329      !--------------------------------------
330      ! CO2 : aqueous CO2
331      !--------------------------------------
332      jn = jn_co2
333      flag_active(jn) = .true.
334      flag_diff(jn)   = .true.
335      nn_remp(jn)     = 1
336      biotr_i_typ(jn) = 'gas'//zblank
337      biotr_i_uni(jn) = 'mmol_m3'//zblank
338
339      !--------------------------------------
340      ! Ika : Ikaite
341      !--------------------------------------
342      jn = jn_ika
343      flag_active(jn) = .true.
344      flag_diff(jn)   = .false.
345      nn_remp(jn)     = 1
346      biotr_i_typ(jn) = 'prc'//zblank
347      biotr_i_uni(jn) = 'mmol_m3'//zblank
348
349      !--------------------------------------
350      ! Cal : Ca2+
351      !--------------------------------------
352      jn = jn_cal
353      flag_active(jn) = .true.
354      flag_diff(jn)   = .true.
355      nn_remp(jn)     = 1
356      biotr_i_typ(jn) = 'prc'//zblank
357      biotr_i_uni(jn) = 'mmol_m3'//zblank
358
359      !--------------------------------------
360      ! Arg : ARgon
361      !--------------------------------------
362      jn = jn_arg
363      flag_active(jn) = .true.
364      flag_diff(jn)   = .true.
365      nn_remp(jn)     = 1
366      biotr_i_typ(jn) = 'gas'//zblank
367      biotr_i_uni(jn) = 'mmol_m3'//zblank
368
369      !--------------------------------------
370      ! Oxy : Oxygen
371      !--------------------------------------
372      jn = jn_oxy
373      flag_active(jn) = .true.
374      flag_diff(jn)   = .true.
375      nn_remp(jn)     = 1
376      biotr_i_typ(jn) = 'gas'//zblank
377      biotr_i_uni(jn) = 'mmol_m3'//zblank
378
379      !--------------------------------------
380      ! Nit : N2
381      !--------------------------------------
382      jn = jn_nit
383      flag_active(jn) = .false.
384      flag_diff(jn)   = .false.
385      nn_remp(jn)     = 1
386      biotr_i_typ(jn) = 'gas'//zblank
387      biotr_i_uni(jn) = 'mmol_m3'//zblank
388
389!
390!-----------------------------------------------------------------------------!
391! 5) Prepare Netcdf Interpolation if required
392!-----------------------------------------------------------------------------!
393!
394      DO jn = 1, ntra_bio
395
396         IF ( flag_active (jn) .AND. ( nn_init(jn) == 3 ) ) 
397     &      ln_initfile = .TRUE.
398
399      END DO
400           
401      IF ( ln_initfile ) THEN
402
403         !--- Open file
404         !--------------
405         CALL CF_OPEN  (filenc,id) ! open forcing file
406
407         CALL CF_READ1D ( filenc, 'nlay', 1, 1, zini )       ! number of layers
408         n_raw = zini(1)
409
410         DO layer = 1, n_raw                                 ! raw grid
411            CALL CF_READ1D ( filenc, 'z_i', layer, 1, zini)
412            z_raw(layer) = zini(1)
413         END DO
414
415         zb_raw(:) = 0.0                                     ! layer interfaces
416         DO layer = 1, n_raw - 1
417            zb_raw(layer)  = ( z_raw(layer) + z_raw(layer+1) ) / 2.
418         END DO
419         zb_raw(n_raw) = ht_i_b(ji)
420
421         zdh_raw(:) = 0.                                     ! layer thicknesses
422         DO layer = 1, n_raw 
423            zdh_raw(layer) = zb_raw(layer) - zb_raw(layer-1)
424         END DO
425
426         zdh_bio(:) = 0.
427         zdh_bio(1:nlay_bio) = deltaz_i_bio(:)
428
429      ENDIF
430
431      !--------------------------------------
432      ! Assign ice concentrations to cbu_i
433      !--------------------------------------
434      DO jn = 1, ntra_bio 
435
436      IF ( flag_active(jn) ) THEN
437
438         IF ( nn_init(jn) .EQ. 1 )        !--- Initialization conserves namelist concentration
439     &      cbu_i_bio(jn,layer_00:nlay_bio) = cbu_i_nml(jn)
440
441         IF ( nn_init(jn) .EQ. 2 ) THEN   !--- Initialization conserves namelist concentration times thickness
442            zstock = cbu_i_nml(jn) * ht_i_b(ji)
443            cbu_i_bio(jn,layer_00:nlay_bio) = zstock / h_bio
444         ENDIF
445
446         IF ( nn_init(jn) .EQ. 3 ) THEN   !--- Profile read from nc file
447
448            DO layer = 1, n_raw
449               CALL CF_READ1D ( filenc, c_nc_name(jn), layer, 1, zini )
450               zc_raw(layer) = zini(1)
451            END DO
452
453            zq_raw(1:n_raw) = zc_raw(1:n_raw) * zdh_raw(1:n_raw)
454
455            CALL ice_phy_relay( n_raw , nlay_bio , 1 , 1 ,
456     &                          zdh_raw, zdh_bio, zq_raw , zq1 )
457
458            cbu_i_bio(jn,1:nlay_bio) = zq1(1:nlay_bio) /
459     &                                 deltaz_i_bio(1:nlay_bio)
460
461         ENDIF
462
463      ENDIF
464
465      END DO
466
467      IF ( ln_initfile )
468     &   CALL CF_CLOSE  (filenc) ! close forcing file
469
470!
471!-----------------------------------------------------------------------------!
472! 7) Chl-a (mg/m3)
473!-----------------------------------------------------------------------------!
474!
475      c_molar = 12. ! molar mass of carbon
476      DO layer = 1, nlay_bio
477         chla_i_bio(layer) = cbu_i_bio(4,layer) * chla_c * c_molar
478      END DO
479
480!
481!-----------------------------------------------------------------------------!
482! 8) Zero un used tracers
483!-----------------------------------------------------------------------------!
484!
485      IF ( nn_bio_opt .EQ. 0 ) THEN
486         cbu_i_bio(jn_eoc,:) = 0.0
487         flag_active(jn_eoc) = .FALSE.
488         flag_diff(jn_eoc)   = .FALSE.
489         c_w_bio(jn_eoc)     = 0.
490         mixr_gas(jn_eoc)    = 0.
491         dmol_gas(jn_eoc)    = 0.
492      ENDIF
493
494!
495!-----------------------------------------------------------------------------!
496! 9) Control prints
497!-----------------------------------------------------------------------------!
498!
499      WRITE(numout,*) '    *** Initial values *** '
500      WRITE(numout,*)
501      DO jn = 1, ntra_bio
502         IF ( flag_active(jn) ) THEN
503            WRITE(numout,*) ' Tracer     : ', biotr_i_nam(jn)
504            WRITE(numout,*) ' ~~~~~~~      '
505            WRITE(numout,*) ' Type       : ', biotr_i_typ(jn)
506            WRITE(numout,*) ' cbu_i_bio  : ', ( cbu_i_bio(jn,jk), 
507     &                   jk = 1, nlay_bio )
508            WRITE(numout,*) ' c_w_bio    : ', c_w_bio(jn)
509            WRITE(numout,*)
510         ENDIF
511      END DO ! jn
512
513      WRITE(numout,*) ' Chla '
514      WRITE(numout,*) ' ~~~~ '
515      WRITE(numout,*) ' chla_i_bio : ', 
516     &                ( chla_i_bio(layer), layer = 1, nlay_bio )
517
518      WRITE(numout,*)
519      WRITE(numout,*) ' End of ice_bio_ini '
520      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
521
522!
523!-----------------------------------------------------------------------------!
524!-- End of ice_bio_ini --
525      RETURN
526 
527      END
Note: See TracBrowser for help on using the repository browser.