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 @ 776

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

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