source: trunk/SOURCES/source_3.20/ice_bio_ini.f @ 4

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

initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 14.9 KB
Line 
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
64         ntra_bio = 21
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.
90      dmol_gas(:)     = 0.
91
92      flag_active(:)    = .false.    ! Activated ?
93      flag_adsorb(:)    = .false.    ! Adsorbed ?
94      flag_diff(:)      = .false.    ! Diffused ?
95      nn_remp(:)        = 1          ! Squeezing remaping
96      biotr_i_nam(:)    = 'xxx'//zblank      ! Name of the tracer
97      biotr_i_typ(:)    = 'xxx'//zblank      ! Type = algal
98      biotr_i_uni(:)    = 'xxx'//zblank ! Units
99      nn_init(:) = 1
100
101      ! Compute layer limits
102      IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) 
103     &   layer_00 = nlay_bio
104
105      ! Assign Tracer Numbers
106      jn_dsi = 1 
107      jn_din = 2
108      jn_dip = 3 
109      jn_aoc = 4 
110      jn_eoc = 5 
111      jn_co2 = 6  ! 20
112      jn_dic = 7  ! 12
113      jn_alk = 8  ! 13
114      jn_ika = 9  ! 14
115      jn_oxy = 10 ! 15
116      jn_nit = 11 ! 16
117      jn_arg = 12 ! 17
118      jn_cal = 13 ! 21
119
120      ! Assign Tracer names
121      biotr_i_nam(jn_dsi) = 'dSi'//zblank
122      biotr_i_nam(jn_din) = 'dIN'//zblank
123      biotr_i_nam(jn_dip) = 'dIP'//zblank
124      biotr_i_nam(jn_aoc) = 'AoC'//zblank
125      biotr_i_nam(jn_eoc) = 'eoC'//zblank
126      biotr_i_nam(jn_co2) = 'CO2'//zblank
127      biotr_i_nam(jn_dic) = 'DIC'//zblank 
128      biotr_i_nam(jn_alk) = 'Alk'//zblank 
129      biotr_i_nam(jn_ika) = 'Ika'//zblank
130      biotr_i_nam(jn_cal) = 'Cal'//zblank
131      biotr_i_nam(jn_arg) = 'Arg'//zblank
132      biotr_i_nam(jn_oxy) = 'Oxy'//zblank
133      biotr_i_nam(jn_nit) = 'Nit'//zblank
134
135!
136!-----------------------------------------------------------------------------!
137! 3) Read tracer parameters
138!-----------------------------------------------------------------------------!
139!
140
141      WRITE(numout,*) ' Tracer parameters ... '
142      WRITE(numout,*)
143
144      OPEN( unit = numtra , file='tracer.param', status='old' )
145     
146      ! initial values
147      READ(numtra,*)
148      READ(numtra,*)
149      READ(numtra,*)
150      READ(numtra,*)
151      READ(numtra,*)
152      READ(numtra,*)
153      READ(numtra,*)
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
169      DO jn_read = 1, 13
170
171         READ(numtra,*)  c_read_name(jn_read), zchar, i_dummy1, zdummy2,
172     &                   zdummy3, zdummy4, zdummy5
173
174         DO jn = 1, ntra_bio
175           IF ( c_read_name(jn_read)//zblank .EQ. biotr_i_nam(jn) ) THEN
176               c_nc_name(jn) = zchar
177               nn_init(jn)   = i_dummy1
178               cbu_i_nml(jn) = zdummy2
179               c_w_bio(jn)   = zdummy3
180               mixr_gas(jn)  = zdummy4
181               dmol_gas(jn)  = zdummy5
182           ENDIF
183         END DO
184
185      ENDDO
186
187      CLOSE( numtra )
188
189      WRITE(numout,*)
190      DO jn = 1, ntra_bio
191         WRITE(numout,*) biotr_i_nam(jn), nn_init(jn) , cbu_i_nml(jn), 
192     &                   c_w_bio(jn)    , mixr_gas(jn), dmol_gas(jn)
193      END DO
194
195!
196!-----------------------------------------------------------------------------!
197! 4) Parameterize tracer behaviour
198!-----------------------------------------------------------------------------!
199!
200      ! flag_active determines a tracer is activated in the code
201      ! flag_adsorb is just there for one specific test and should not be activated
202      ! flag_diff tells whether a tracer should be diffused with brine or not
203      ! nn_remp gives the type of remapping - 1 = linear remapping (very diffusive); 2 = "squeezing remapping"
204      ! biotr_i_typ can be either - only gas is useful because flag_diff makes the difference between solutes and particles
205      !    - 'sol' -> solute (like salt)
206      !    - 'prc' -> particulate (e.g. organic matter)
207      !    - 'gas' -> gas
208
209      !------------------------
210      ! DSi : Dissolved silica
211      !------------------------
212      jn = jn_dsi
213      flag_active(jn) = .true.
214      flag_adsorb(jn) = .false.
215      flag_diff(jn)   = .true.
216      nn_remp(jn)     = 1
217      biotr_i_typ(jn) = 'sol'//zblank
218      biotr_i_uni(jn) = 'mmol_m3'//zblank
219
220      !------------------------------------
221      ! DIN : dissolved inorganic nitrogen
222      !------------------------------------
223      jn = jn_din
224      flag_active(jn) = .true.
225      flag_adsorb(jn) = .false.
226      flag_diff(jn)   = .true.
227      nn_remp(jn)     = 1
228      biotr_i_typ(jn) = 'sol'//zblank
229      biotr_i_uni(jn) = 'mmol_m3'//zblank
230
231      !----------------------------
232      ! DIP : Dissolved phosphorus
233      !----------------------------
234      jn = jn_dip
235      ! Units = ( mmole PO4 m-3 ) or ( micrmol PO4 l-1 )
236      flag_active(jn) = .true.
237      flag_adsorb(jn) = .false.
238      flag_diff(jn)   = .true.
239      nn_remp(jn)     = 1
240      biotr_i_typ(jn) = 'sol'//zblank
241      biotr_i_uni(jn) = 'mmol_m3'//zblank
242
243      !--------------------------------------
244      ! AoC : Ice algae organic carbon
245      !--------------------------------------
246      jn = jn_aoc
247      flag_active(jn) = .true.
248      flag_adsorb(jn) = .false.
249      flag_diff(jn)   = .false.
250      nn_remp(jn)     = 1
251      biotr_i_typ(jn) = 'prc'//zblank
252      biotr_i_uni(jn) = 'mmol_m3'//zblank
253
254      !--------------------------------------
255      ! eoC : Detritus organic Carbon
256      !--------------------------------------
257      jn = jn_eoc
258      flag_active(jn) = .true.
259      flag_adsorb(jn) = .false.
260      flag_diff(jn)   = .false.
261      nn_remp(jn)     = 1
262      biotr_i_typ(jn) = 'prc'//zblank
263      biotr_i_uni(jn) = 'mmol_m3'//zblank
264
265      !--------------------------------------
266      ! DIC : Dissolved Inorganic Carbon
267      !--------------------------------------
268      jn = jn_dic
269      flag_active(jn) = .true.
270      flag_adsorb(jn) = .false.
271      flag_diff(jn)   = .true.
272      nn_remp(jn)     = 1
273      biotr_i_typ(jn) = 'sol'//zblank
274      biotr_i_uni(jn) = 'mmol_m3'//zblank
275
276      !--------------------------------------
277      ! Alk : Total Alkalinity
278      !--------------------------------------
279      jn = jn_alk
280      flag_active(jn) = .true.
281      flag_adsorb(jn) = .false.
282      flag_diff(jn)   = .true.
283      nn_remp(jn)     = 1
284      biotr_i_typ(jn) = 'sol'//zblank
285      biotr_i_uni(jn) = 'mmol_m3'//zblank
286
287      !--------------------------------------
288      ! CO2 : aqueous CO2
289      !--------------------------------------
290      jn = jn_co2
291      flag_active(jn) = .true.
292      flag_adsorb(jn) = .false.
293      flag_diff(jn)   = .true.
294      nn_remp(jn)     = 1
295      biotr_i_typ(jn) = 'gas'//zblank
296      biotr_i_uni(jn) = 'mmol_m3'//zblank
297
298      !--------------------------------------
299      ! Ika : Ikaite
300      !--------------------------------------
301      jn = jn_ika
302      flag_active(jn) = .true.
303      flag_adsorb(jn) = .false.
304      flag_diff(jn)   = .false.
305      nn_remp(jn)     = 1
306      biotr_i_typ(jn) = 'prc'//zblank
307      biotr_i_uni(jn) = 'mmol_m3'//zblank
308
309      !--------------------------------------
310      ! Cal : Ca2+
311      !--------------------------------------
312      jn = jn_cal
313      flag_active(jn) = .true.
314      flag_adsorb(jn) = .false.
315      flag_diff(jn)   = .true.
316      nn_remp(jn)     = 1
317      biotr_i_typ(jn) = 'prc'//zblank
318      biotr_i_uni(jn) = 'mmol_m3'//zblank
319
320      !--------------------------------------
321      ! Arg : ARgon
322      !--------------------------------------
323      jn = jn_arg
324      flag_active(jn) = .true.
325      flag_adsorb(jn) = .false.
326      flag_diff(jn)   = .true.
327      nn_remp(jn)     = 1
328      biotr_i_typ(jn) = 'gas'//zblank
329      biotr_i_uni(jn) = 'mmol_m3'//zblank
330
331      !--------------------------------------
332      ! Oxy : Oxygen
333      !--------------------------------------
334      jn = jn_oxy
335      flag_active(jn) = .true.
336      flag_adsorb(jn) = .false.
337      flag_diff(jn)   = .true.
338      nn_remp(jn)     = 1
339      biotr_i_typ(jn) = 'gas'//zblank
340      biotr_i_uni(jn) = 'mmol_m3'//zblank
341
342      !--------------------------------------
343      ! Nit : N2
344      !--------------------------------------
345      jn = jn_nit
346      flag_active(jn) = .false.
347      flag_adsorb(jn) = .false.
348      flag_diff(jn)   = .false.
349      nn_remp(jn)     = 1
350      biotr_i_typ(jn) = 'gas'//zblank
351      biotr_i_uni(jn) = 'mmol_m3'//zblank
352
353!
354!-----------------------------------------------------------------------------!
355! 5) Prepare Netcdf Interpolation if required
356!-----------------------------------------------------------------------------!
357!
358      DO jn = 1, ntra_bio
359
360         IF ( flag_active (jn) .AND. ( nn_init(jn) == 3 ) ) 
361     &      ln_initfile = .TRUE.
362
363      END DO
364           
365      IF ( ln_initfile ) THEN
366
367         !--- Open file
368         !--------------
369         CALL CF_OPEN  (filenc,id) ! open forcing file
370
371         CALL CF_READ1D ( filenc, 'nlay', 1, 1, zini )       ! number of layers
372         n_raw = zini(1)
373
374         DO layer = 1, n_raw                                 ! raw grid
375            CALL CF_READ1D ( filenc, 'z_i', layer, 1, zini)
376            z_raw(layer) = zini(1)
377         END DO
378
379         zb_raw(:) = 0.0                                     ! layer interfaces
380         DO layer = 1, n_raw - 1
381            zb_raw(layer)  = ( z_raw(layer) + z_raw(layer+1) ) / 2.
382         END DO
383         zb_raw(n_raw) = ht_i_b(ji)
384
385         zdh_raw(:) = 0.                                     ! layer thicknesses
386         DO layer = 1, n_raw 
387            zdh_raw(layer) = zb_raw(layer) - zb_raw(layer-1)
388         END DO
389
390         zdh_bio(:) = 0.
391         zdh_bio(1:nlay_bio) = deltaz_i_bio(:)
392
393      ENDIF
394
395      !--------------------------------------
396      ! Assign ice concentrations to cbu_i
397      !--------------------------------------
398      DO jn = 1, ntra_bio 
399
400      IF ( flag_active(jn) ) THEN
401
402         IF ( nn_init(jn) .EQ. 1 )        !--- Initialization conserves namelist concentration
403     &      cbu_i_bio(jn,layer_00:nlay_bio) = cbu_i_nml(jn)
404
405         IF ( nn_init(jn) .EQ. 2 ) THEN   !--- Initialization conserves namelist concentration times thickness
406            zstock = cbu_i_nml(jn) * ht_i_b(ji)
407            cbu_i_bio(jn,layer_00:nlay_bio) = zstock / h_bio
408         ENDIF
409
410         IF ( nn_init(jn) .EQ. 3 ) THEN   !--- Profile read from nc file
411
412            DO layer = 1, n_raw
413               CALL CF_READ1D ( filenc, c_nc_name(jn), layer, 1, zini )
414               zc_raw(layer) = zini(1)
415            END DO
416
417            zq_raw(1:n_raw) = zc_raw(1:n_raw) * zdh_raw(1:n_raw)
418
419            CALL ice_phy_relay( n_raw , nlay_bio , 1 , 1 ,
420     &                          zdh_raw, zdh_bio, zq_raw , zq1 )
421
422            cbu_i_bio(jn,1:nlay_bio) = zq1(1:nlay_bio) /
423     &                                 deltaz_i_bio(1:nlay_bio)
424
425         ENDIF
426
427      ENDIF
428
429      END DO
430
431      IF ( ln_initfile )
432     &   CALL CF_CLOSE  (filenc) ! close forcing file
433
434!
435!-----------------------------------------------------------------------------!
436! 7) Chl-a (mg/m3)
437!-----------------------------------------------------------------------------!
438!
439      c_molar = 12. ! molar mass of carbon
440      DO layer = 1, nlay_bio
441         chla_i_bio(layer) = cbu_i_bio(4,layer) * chla_c * c_molar
442      END DO
443
444!
445!-----------------------------------------------------------------------------!
446! 8) Zero un used tracers
447!-----------------------------------------------------------------------------!
448!
449      IF ( nn_bio_opt .EQ. 0 ) THEN
450         cbu_i_bio(jn_eoc,:) = 0.0
451         flag_active(jn_eoc) = .FALSE.
452         flag_adsorb(jn_eoc) = .FALSE.
453         flag_diff(jn_eoc)   = .FALSE.
454         c_w_bio(jn_eoc)     = 0.
455         mixr_gas(jn_eoc)    = 0.
456         dmol_gas(jn_eoc)    = 0.
457      ENDIF
458
459!
460!-----------------------------------------------------------------------------!
461! 9) Control prints
462!-----------------------------------------------------------------------------!
463!
464      WRITE(numout,*) '    *** Initial values *** '
465      WRITE(numout,*)
466      DO jn = 1, ntra_bio
467         IF ( flag_active(jn) ) THEN
468            WRITE(numout,*) ' Tracer     : ', biotr_i_nam(jn)
469            WRITE(numout,*) ' ~~~~~~~      '
470            WRITE(numout,*) ' Type       : ', biotr_i_typ(jn)
471            WRITE(numout,*) ' cbu_i_bio  : ', ( cbu_i_bio(jn,jk), 
472     &                   jk = 1, nlay_bio )
473            WRITE(numout,*) ' c_w_bio    : ', c_w_bio(jn)
474            WRITE(numout,*)
475         ENDIF
476      END DO ! jn
477
478      WRITE(numout,*) ' Chla '
479      WRITE(numout,*) ' ~~~~ '
480      WRITE(numout,*) ' chla_i_bio : ', 
481     &                ( chla_i_bio(layer), layer = 1, nlay_bio )
482
483      WRITE(numout,*)
484      WRITE(numout,*) ' End of ice_bio_ini '
485      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
486
487!
488!-----------------------------------------------------------------------------!
489!-- End of ice_bio_ini --
490      RETURN
491 
492      END
Note: See TracBrowser for help on using the repository browser.