New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcini.pisces.h90 in trunk/NEMO/TOP_SRC/SMS – NEMO

source: trunk/NEMO/TOP_SRC/SMS/trcini.pisces.h90 @ 728

Last change on this file since 728 was 728, checked in by cetlod, 17 years ago

add tmask_i in global sum, see ticket:17

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