source: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini.pisces.h90 @ 764

Last change on this file since 764 was 764, checked in by gm, 13 years ago

dev_001_GM - create new directory and move files only

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.4 KB
Line 
1   !!======================================================================
2   !!                         ***  trcini.pisces.h90  ***
3   !! TOP :   Initialisation of PISCES biological model
4   !!======================================================================
5   !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
6   !!              -   !  1999-10  (O. Aumont, C. Le Quere)
7   !!              -   !  2002     (O. Aumont)  PISCES
8   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
9   !!----------------------------------------------------------------------
10
11#  include "domzgr_substitute.h90"
12#  include "passivetrc_substitute.h90"
13   !!----------------------------------------------------------------------
14   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005)
15   !! $Id$
16   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
17   !!----------------------------------------------------------------------
18
19CONTAINS
20
21   SUBROUTINE trc_ini
22      !!----------------------------------------------------------------------
23      !!                   ***  ROUTINE trc_ini ***
24      !!
25      !! ** Purpose :   Initialisation of PISCES biological and chemical variables
26      !!----------------------------------------------------------------------
27      USE iom
28      !!
29      INTEGER :: ji,jj,jk
30      INTEGER :: ichl,iband,jm
31      INTEGER , PARAMETER :: jpmois = 12, jpan   = 1
32
33      REAL(wp) :: zcoef
34      REAL(wp) :: ztoto,expide,denitide,ztra,zmaskt
35      REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo
36      REAL(wp) , DIMENSION (jpi,jpj,jpk) :: cmask
37
38      INTEGER :: numriv, numdust, numbath, numdep
39      INTEGER :: numlight
40
41#if defined key_trc_kriest
42      INTEGER  ::   jn, kiter
43      REAL(wp) ::   znum, zdiv
44      REAL(wp) ::   zws,zwr, zwl,wmax, xnummax, &
45      REAL(wp) ::   zmin, zmax, zl, zr, xacc
46#endif
47      !!----------------------------------------------------------------------
48
49
50      IF(lwp) WRITE(numout,*)
51      IF(lwp) WRITE(numout,*) ' trc_ini :   PISCES biological and chemical initialisation'
52      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
53
54
55      !                                            ! Time-step
56      rfact   = rdttra(1) * float(ndttrc)          ! ---------
57      rfactr  = 1. / rfact
58      rfact2  = rfact / float(nrdttrc)
59      rfact2r = 1. / rfact2
60
61      IF(lwp) WRITE(numout,*) '    Tracer  time step=', rfact, ' rdt = ', rdt
62      IF(lwp) write(numout,*) '    Biology time step=', rfact2
63
64
65      !                                            ! Dust input from the atmosphere
66      IF( bdustfer ) THEN                          ! ------------------------------
67         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
68         CALL iom_open ( 'dust.orca.nc', numdust )
69         DO jm = 1, jpmois
70            CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm )
71         END DO
72         CALL iom_close( numdust )
73      ELSE
74         dustmo(:,:,:) = 0.e0
75      ENDIF
76
77
78      !                                            ! Nutrient input from rivers
79      IF( briver ) THEN                            ! --------------------------
80         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
81         CALL iom_open ( 'river.orca.nc', numriv )
82         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
83         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
84         CALL iom_close( numriv )
85      ELSE
86         river   (:,:) = 0.e0
87         riverdoc(:,:) = 0.e0
88      endif
89
90      !                                            ! Nutrient input from dust
91      IF( bndepo ) THEN                            ! ------------------------
92         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
93         CALL iom_open ( 'ndeposition.orca.nc', numdep )
94         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
95         CALL iom_close( numdep )
96      ELSE
97         ndepo(:,:) = 0.e0
98      ENDIF
99
100      !                                            ! Coastal and island masks
101      IF( bsedinput ) THEN                         ! ------------------------
102         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron '
103         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
104         CALL iom_open ( 'bathy.orca.nc', numbath )
105         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
106         CALL iom_close( numbath )
107         !
108         DO jk = 1, 5
109            DO jj = 2, jpjm1
110               DO ji = 2, jpim1
111                  IF( tmask(ji,jj,jk) /= 0. ) THEN
112                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
113                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
114                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
115                  ENDIF
116               END DO
117            END DO
118         END DO
119         DO jk = 1, jpk
120            DO jj = 1, jpj
121               DO ji = 1, jpi
122                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
123                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
124                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
125               END DO
126            END DO
127         END DO   
128      ELSE
129         cmask(:,:,:) = 0.e0
130      ENDIF
131
132      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
133
134
135      !                                            ! total atmospheric supply of Si
136      !                                            ! ------------------------------
137      sumdepsi = 0.e0
138      DO jm = 1, jpmois
139         DO jj = 2, jpjm1
140            DO ji = 2, jpim1
141               sumdepsi = sumdepsi + dustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
142                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
143            END DO
144         END DO
145      END DO
146      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
147
148      !                                            ! N/P and Si releases due to coastal rivers
149      !                                            ! -----------------------------------------
150      DO jj = 1, jpj
151         DO ji = 1, jpi
152            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
153            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
154            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
155            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
156         END DO
157      END DO
158      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
159      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
160
161      rivpo4input=0.e0
162      rivalkinput=0.e0
163      nitdepinput=0.e0
164      DO jj = 2 , jpjm1
165         DO ji = 2, jpim1
166            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
167            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
168            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
169            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
170         END DO
171      END DO
172      IF( lk_mpp ) THEN
173         CALL mpp_sum( rivpo4input )  ! sum over the global domain
174         CALL mpp_sum( rivalkinput )  ! sum over the global domain
175         CALL mpp_sum( nitdepinput )  ! sum over the global domain
176      ENDIF
177
178
179      !                                            !    Coastal supply of iron
180      !                                            ! -------------------------
181      DO jk = 1, jpkm1
182         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
183      END DO
184      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
185
186
187
188#if defined key_trc_kriest
189      !----------------------------------------------------------------------
190      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
191      !  Search of the maximum number of particles in aggregates for each k-level. 
192      !  Bissection Method
193      !--------------------------------------------------------------------
194      WRITE(numout,*)
195      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
196           
197      xacc     =  0.001
198      kiter    = 50
199      zmin     =  1.10
200      zmax     = xkr_mass_max / xkr_mass_min
201      xkr_frac = zmax
202
203      DO jk =1,jpk
204         zl = zmin
205         zr = zmax
206         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
207         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
208         znum = zl - 1.
209         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
210            & - ( xkr_wsbio_max * xkr_eta * znum * &
211            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
212            & - wmax
213 
214         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
215         znum = zr - 1.
216         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
217            & - ( xkr_wsbio_max * xkr_eta * znum * &
218            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
219            & - wmax
220
221iflag:  DO jn = 1, kiter               
222           IF( zwl == 0.e0 ) THEN
223              xnummax = zl
224           ELSE IF ( zwr == 0.e0 ) THEN
225              xnummax = zr
226           ELSE
227              xnummax = ( zr + zl ) / 2.
228              zdiv = xkr_zeta + xkr_eta - xkr_eta * xnummax
229              znum = xnummax - 1.
230              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
231                 & - ( xkr_wsbio_max * xkr_eta * znum * &
232                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
233                 & - wmax
234              IF( zws * zwl < 0. ) THEN
235                 zr = xnummax
236              ELSE
237                 zl = xnummax
238              ENDIF
239              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
240              znum = zl - 1.
241              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
242                 & - ( xkr_wsbio_max * xkr_eta * znum * &
243                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
244                 & - wmax
245             
246              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
247              znum = zr - 1.
248              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
249                 & - ( xkr_wsbio_max * xkr_eta * znum * &
250                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
251                 & - wmax
252
253              IF ( ABS ( zws )  <= xacc ) EXIT iflag
254
255           ENDIF
256                   
257        END DO iflag
258               
259        xnumm(jk) = xnummax
260        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
261       
262     END DO
263
264#endif
265
266      !----------------------------------------------------------------------
267      ! Initialize biological variables
268      !----------------------------------------------------------------------
269      ! Set biological ratios
270      ! ---------------------
271      rno3   = (16.+2.) / 122.
272      po4r   =   1.e0   / 122.
273      o2nit  =  32.     / 122.
274      rdenit =  97.6    /  16.
275      o2ut   = 140.     / 122.
276
277      !----------------------------------------------------------------------
278      ! Initialize chemical variables
279      !----------------------------------------------------------------------
280
281      ! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
282      ! ----------------------------------------------------------
283      atcox = 0.20946
284
285      ! Set lower/upper limits for temperature and salinity
286      ! ---------------------------------------------------
287      salchl = 1.e0 / 1.80655
288      calcon = 1.03e-2
289
290      ! Set coefficients for apparent solubility equilibrium of calcite
291      ! Millero et al. 1995 from Mucci 1983
292      ! --------------------------------------------------------------
293      akcc1 = -171.9065
294      akcc2 =   -0.077993
295      akcc3 = 2839.319
296      akcc4 =   71.595
297      akcc5 =   -0.77712
298      akcc6 =    0.0028426
299      akcc7 =  178.34
300      akcc8 =   -0.07711
301      akcc9 =    0.0041249
302
303      ! Set coefficients for seawater pressure correction
304      ! -------------------------------------------------
305      devk1(1) = -25.5
306      devk2(1) =   0.1271
307      devk3(1) =   0.e0
308      devk4(1) =  -3.08E-3
309      devk5(1) =   0.0877E-3
310      !
311      devk1(2) = -15.82
312      devk2(2) =  -0.0219
313      devk3(2) =   0.e0
314      devk4(2) =   1.13E-3
315      devk5(2) =  -0.1475E-3
316      !
317      devk1(3) = -29.48
318      devk2(3) =   0.1622
319      devk3(3) =   2.608E-3
320      devk4(3) =  -2.84E-3
321      devk5(3) =   0.e0
322      !
323      devk1(4) = -14.51
324      devk2(4) =   0.1211
325      devk3(4) =  -0.321E-3
326      devk4(4) =  -2.67E-3
327      devk5(4) =   0.0427E-3
328      !
329      devk1(5) = -23.12
330      devk2(5) =   0.1758
331      devk3(5) =  -2.647E-3
332      devk4(5) =  -5.15E-3
333      devk5(5) =   0.09E-3
334      !
335      devk1(6) = -26.57
336      devk2(6) =   0.2020
337      devk3(6) =  -3.042E-3
338      devk4(6) =  -4.08E-3
339      devk5(6) =   0.0714E-3
340      !
341      devk1(7) = -25.60
342      devk2(7) =   0.2324
343      devk3(7) =  -3.6246E-3
344      devk4(7) =  -5.13E-3
345      devk5(7) =   0.0794E-3
346      !
347      ! For calcite with Edmond and Gieske 1970
348      !     devkst = 0.23
349      !     devks  = 35.4
350      ! Millero 95 takes this depth dependance for calcite
351      devk1(8) = -48.76
352      devk2(8) =   0.5304
353      devk3(8) =   0.e0
354      devk4(8) = -11.76E-3
355      devk5(8) =   0.3692E-3
356      !
357      ! Coefficients for sulfate and fluoride
358      devk1(9) = -18.03
359      devk2(9) =   0.0466
360      devk3(9) =   0.316e-3
361      devk4(9) =  -4.53e-3
362      devk5(9) =   0.09e-3
363
364      devk1(10) = -9.78
365      devk2(10) = -0.0090
366      devk3(10) = -0.942e-3
367      devk4(10) = -3.91e-3
368      devk5(10) =  0.054e-3
369
370
371      ! Set universal gas constants
372      ! ---------------------------
373      rgas  = 83.143
374      oxyco =  1.e0 / 22.4144
375
376      ! Set boron constants
377      ! -------------------
378      bor1 = 0.00023
379      bor2 = 1.e0 / 10.82
380
381      ! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
382      ! -----------------------------------------------------------------
383      c00 = -60.2409
384      c01 =  93.4517
385      c02 =  23.3585
386      c03 =   0.023517
387      c04 =  -0.023656
388      c05 =   0.0047036
389      !
390      ca0 = -162.8301
391      ca1 =  218.2968
392      ca2 =   90.9241
393      ca3 =   -1.47696
394      ca4 =    0.025695
395      ca5 =   -0.025225
396      ca6 =    0.0049867
397
398      ! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
399      ! ---------------------------------------------------------------------
400      c10 = -3670.7
401      c11 =    62.008
402      c12 =    -9.7944
403      c13 =     0.0118
404      c14 =    -0.000116
405
406      ! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
407      ! ---------------------------------------------------------------------
408      c20 = -1394.7
409      c21 =    -4.777
410      c22 =     0.0184
411      c23 =    -0.000118
412
413      ! Set constants for calculate concentrations for sulfate and fluoride
414      ! sulfates (Morris & Riley 1966)
415      !----------------------------------------------------------------------
416      st1 = 0.14
417      st2 = 1.e0 / 96.062
418
419      ! fluoride
420      ! --------
421      ft1 = 0.000067
422      ft2 = 1.e0 / 18.9984
423
424      ! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP)
425      !----------------------------------------------------------
426      ks0  =    141.328
427      ks1  =  -4276.1
428      ks2  =    -23.093
429      ks3  = -13856.
430      ks4  =    324.57
431      ks5  =    -47.986
432      ks6  =  35474.
433      ks7  =   -771.54
434      ks8  =    114.723
435      ks9  =  -2698.
436      ks10 =   1776.
437      ks11 =      1.
438      ks12 =     -0.001005
439
440      ! fluorides (Dickson & Riley 1979 change to mol/kg soln)
441      !-------------------------------------------------------
442      kf0 =  -12.641
443      kf1 = 1590.2
444      kf2 =    1.525
445      kf3 =    1.0
446      kf4 =   -0.001005
447
448      ! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
449      ! ------------------------------------------------------------------
450      cb0  = -8966.90
451      cb1  = -2890.53
452      cb2  =   -77.942
453      cb3  =     1.728
454      cb4  =    -0.0996
455      cb5  =   148.0248
456      cb6  =   137.1942
457      cb7  =     1.62142
458      cb8  =   -24.4344
459      cb9  =   -25.085
460      cb10 =    -0.2474
461      cb11 =     0.053105
462
463      ! Set coeff. for dissoc. of water (Dickson and Riley, 1979,
464      !   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415
465      !   after pers. commun. to B. Bacastow, 1988)
466      ! ---------------------------------------------------------
467      cw0 = -13847.26
468      cw1 =    148.9652
469      cw2 =    -23.6521
470      cw3 =    118.67
471      cw4 =     -5.977
472      cw5 =      1.0495
473      cw6 =     -0.01615
474
475      ! Set coeff. for dissoc. of phosphate (Millero (1974)
476      ! ---------------------------------------------------
477      cp10 =   115.54
478      cp11 = -4576.752
479      cp12 =   -18.453
480      cp13 =  -106.736
481      cp14 =     0.69171
482      cp15 =    -0.65643
483      cp16 =    -0.01844
484      !
485      cp20 =   172.1033
486      cp21 = -8814.715
487      cp22 =   -27.927
488      cp23 =  -160.340
489      cp24 =     1.3566
490      cp25 =     0.37335
491      cp26 =    -0.05778
492      !
493      cp30 =   -18.126
494      cp31 = -3070.75
495      cp32 =    17.27039
496      cp33 =     2.81197
497      cp34 =   -44.99486
498      cp35 =    -0.09984
499
500      ! Set coeff. for dissoc. of phosphate (Millero (1974)
501      ! ---------------------------------------------------
502      cs10 =   117.385
503      cs11 = -8904.2
504      cs12 =   -19.334
505      cs13 =  -458.79
506      cs14 =     3.5913
507      cs15 =   188.74
508      cs16 =    -1.5998
509      cs17 =   -12.1652
510      cs18 =     0.07871
511      cs19 =     0.e0
512      cs20 =     1.e0
513      cs21 =    -0.001005
514
515
516      ! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
517      ! ----------------------------------------------------------------
518      ox0 = -58.3877
519      ox1 =  85.8079
520      ox2 =  23.8439
521      ox3 =  -0.034892
522      ox4 =   0.015568
523      ox5 =  -0.0019387
524
525      !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
526      !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
527
528      CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
529         &           1, numout, .TRUE., 1 )
530      DO ichl = 1,61
531         READ(numlight,*) ztoto, ( xkrgb(iband,ichl), iband = 1,3 )
532      END DO
533      CLOSE(numlight)
534
535
536      CALL p4zche        ! initialize the chemical constants
537
538
539      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry
540
541
542      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
543      !
544   END SUBROUTINE trc_ini
Note: See TracBrowser for help on using the repository browser.