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

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

dev_001_GM - create 1 trcini_ module by trc model (CFC, LOBSTER, PISCES..) - never compiled

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