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.F90 in branches/dev_001_GM/NEMO/TOP_SRC/PISCES – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90 @ 771

Last change on this file since 771 was 771, checked in by gm, 16 years ago

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