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

Last change on this file since 433 was 341, checked in by opalod, 19 years ago

nemo_v1_update_028 : CT : add missing headers
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 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$
40   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
41      !!----------------------------------------------------------------------
42      !! local declarations
43      !! ==================
44      INTEGER :: ji,jj,jk
45      INTEGER :: ichl,iband,mo
46      INTEGER , PARAMETER :: jpmois = 12,      &
47         jpan   = 1
48
49      REAL(wp) :: xtoto,expide,denitide,ztra,zmaskt
50      REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo
51      CHARACTER (len=34) :: clname
52
53      INTEGER :: ipi,ipj,ipk,itime
54      INTEGER , DIMENSION (jpmois) :: istep
55      INTEGER , DIMENSION (jpan) :: istep0
56      REAL(wp) :: zsecond, zdate0
57      REAL(wp) , DIMENSION (jpi,jpj) :: zlon,zlat
58      REAL(wp), DIMENSION (jpk) :: zlev
59      INTEGER :: numriv,numdust,numbath,numdep
60
61      !! 1. initialization
62      !! -----------------
63
64      !! computation of the record length for direct access FILE
65      !! this length depend of 512 for the t3d machine
66      !!
67      rfact = rdttra(1) * float(ndttrc)
68      rfactr = 1./rfact
69      IF(lwp) WRITE(numout,*) ' Tracer time step=',rfact,' rdt=',rdt
70      rfact2= rfact / float(nrdttrc)
71      rfact2r = 1./rfact2
72      IF(lwp) write(numout,*) ' Biology time step=',rfact2
73
74      !!    INITIALISE DUST INPUT FROM ATMOSPHERE
75      !!    -------------------------------------
76
77      IF (bdustfer) THEN
78         clname='dust.orca.nc'
79         CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0        &
80            &      ,zlon,zlat,zlev,itime,istep,zdate0,zsecond,numdust)
81         CALL flinget(numdust,'dust',jpidta,jpjdta,0,jpmois,1,                 &
82            &        12,mig(1),nlci,mjg(1),nlcj,dustmo(1:nlci,1:nlcj,:) )
83         CALL flinclo(numdust)
84
85         ! Extra-halo initialization in MPP
86         IF( lk_mpp ) THEN
87            DO ji = nlci+1, jpi
88               dustmo(ji,:,:) = dustmo(1,:,:)
89            ENDDO
90            DO jj = nlcj+1, jpj
91               dustmo(:,jj,:)=dustmo(:,1,:)
92            ENDDO
93         ENDIF
94      ELSE
95         dustmo(:,:,:)=0.
96      ENDIF
97
98      !!    INITIALISE THE NUTRIENT INPUT BY RIVERS
99      !!    ---------------------------------------
100
101      IF (briver) THEN
102         clname='river.orca.nc'
103         CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0        &
104            &      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numriv)
105         CALL flinget(numriv,'riverdic',jpidta,jpjdta,0,jpan,1,                &
106            &        1,mig(1),nlci,mjg(1),nlcj,river(1:nlci,1:nlcj) )
107         CALL flinget(numriv,'riverdoc',jpidta,jpjdta,0,jpan,1,                &
108            &        1,mig(1),nlci,mjg(1),nlcj,riverdoc(1:nlci,1:nlcj) )
109         CALL flinclo(numriv)
110
111         ! Extra-halo initialization in MPP
112         IF( lk_mpp ) THEN
113            DO ji = nlci+1, jpi
114               river(ji,:) = river(1,:)
115               riverdoc(ji,:) = riverdoc(1,:)
116            ENDDO
117            DO jj = nlcj+1, jpj
118               river(:,jj)=river(:,1)
119               riverdoc(:,jj) = riverdoc(:,1)
120            ENDDO
121         ENDIF
122
123      ELSE
124         river(:,:)=0.
125         riverdoc(:,:)=0.
126      endif
127
128      !!    INITIALISE THE N INPUT BY DUST
129      !!  ---------------------------------------
130
131      IF (bndepo) THEN
132         clname='ndeposition.orca.nc'
133         CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,0        &
134            &      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numdep)
135         CALL flinget(numdep,'ndep',jpidta,jpjdta,0,jpan,1,                   &
136            &        1,mig(1),nlci,mjg(1),nlcj,ndepo(1:nlci,1:nlcj) )
137         CALL flinclo(numdep)
138
139         ! Extra-halo initialization in MPP
140         IF( lk_mpp ) THEN
141            DO ji = nlci+1, jpi
142               ndepo(ji,:) = ndepo(1,:)
143            ENDDO
144            DO jj = nlcj+1, jpj
145               ndepo(:,jj)=ndepo(:,1)
146            ENDDO
147         ENDIF
148
149      ELSE
150         ndepo(:,:)=0.
151      ENDIF
152
153      !!    Computation of the coastal mask.
154      !!    Computation of an island mask to enhance coastal supply
155      !!    of iron
156      !!    -------------------------------------------------------
157
158      IF (bsedinput) THEN
159         clname='bathy.orca.nc'
160         CALL flinopen(clname,mig(1),nlci,mjg(1),nlcj,.false.,ipi,ipj,ipk      &
161            &      ,zlon,zlat,zlev,itime,istep0,zdate0,zsecond,numbath)
162         CALL flinget(numbath,'bathy',jpidta,jpjdta,jpk,jpan,1,                &
163            &        1,mig(1),nlci,mjg(1),nlcj,cmask(1:nlci,1:nlcj,1:jpk) )
164         CALL flinclo(numbath)
165
166         do jk=1,5
167            do jj=2,jpj-1
168               do ji=2,jpi-1
169                  if (tmask(ji,jj,jk).ne.0) then
170                     zmaskt=tmask(ji+1,jj,jk)*tmask(ji-1,jj,jk)*tmask(ji,jj+1,jk)    &
171                        &          *tmask(ji,jj-1,jk)*tmask(ji,jj,jk+1)
172                     if (zmaskt.eq.0) then
173                        cmask(ji,jj,jk)=0.1
174                     endif
175                  endif
176               end do
177            end do
178         end do
179
180
181         ! Extra-halo initialization in MPP
182         IF( lk_mpp ) THEN
183            DO ji = nlci+1, jpi
184               cmask(ji,:,:) = cmask(1,:,:)
185            ENDDO
186            DO jj = nlcj+1, jpj
187               cmask(:,jj,:)=cmask(:,1,:)
188            ENDDO
189         ENDIF
190
191         DO jk = 1, jpk
192            DO jj = 1, jpj
193               DO ji = 1, jpi
194                  expide=min(8.,(fsdept(ji,jj,jk)/500.)**(-1.5))
195                  denitide=-0.9543+0.7662*log(expide)-0.235*log(expide)**2
196                  cmask(ji,jj,jk)=cmask(ji,jj,jk)*min(1.,exp(denitide)/0.5)
197               END DO
198            END DO
199         END DO
200
201      ELSE
202         cmask(:,:,:)=0.
203      ENDIF
204
205      ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
206      CALL lbc_lnk( cmask , 'T', 1. )
207
208      !!     Computation of the total atmospheric supply of Si
209      !!     -------------------------------------------------
210
211      sumdepsi=0.
212      DO mo=1,12
213         DO jj=2,jpjm1
214            DO ji=2,jpim1
215               sumdepsi=sumdepsi+dustmo(ji,jj,mo)/(12.*rmoss)*8.8        &
216                  *0.075/28.1*e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,1)
217            END DO
218         END DO
219      END DO
220
221      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
222
223      !!    COMPUTATION OF THE N/P RELEASE DUE TO COASTAL RIVERS
224      !!    COMPUTATION OF THE Si RELEASE DUE TO COASTAL RIVERS
225      !!    ---------------------------------------------------
226
227      DO jj=1,jpj
228         DO ji=1,jpi
229            cotdep(ji,jj)=river(ji,jj)*1E9/(12.*raass                          &
230               *e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1)+rtrn)*tmask(ji,jj,1)
231            rivinp(ji,jj)=(river(ji,jj)+riverdoc(ji,jj))*1E9                   &
232               /(31.6*raass*e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1)+rtrn)   &
233               *tmask(ji,jj,1)
234            nitdep(ji,jj)=7.6*ndepo(ji,jj)*tmask(ji,jj,1)/(14E6*raass          &
235               *fse3t(ji,jj,1)+rtrn)
236         END DO
237      END DO
238      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
239      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
240
241      rivpo4input=0.
242      rivalkinput=0.
243      nitdepinput=0.
244      DO jj=2,jpjm1
245         DO ji=2,jpim1
246            rivpo4input=rivpo4input+rivinp(ji,jj)*(e1t(ji,jj)*e2t(ji,jj)    &
247               *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
248            rivalkinput=rivalkinput+cotdep(ji,jj)*(e1t(ji,jj)*e2t(ji,jj)    &
249               *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
250            nitdepinput=nitdepinput+nitdep(ji,jj)*(e1t(ji,jj)*e2t(ji,jj)    &
251               *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass
252         END DO
253      END DO
254
255      IF( lk_mpp ) THEN
256         CALL mpp_sum( rivpo4input )  ! sum over the global domain
257         CALL mpp_sum( rivalkinput )  ! sum over the global domain
258         CALL mpp_sum( nitdepinput )  ! sum over the global domain
259      ENDIF
260
261
262      !!    Coastal supply of iron
263      !!    ----------------------
264
265      DO jk=1,jpkm1
266         ironsed(:,:,jk)=sedfeinput*cmask(:,:,jk)         &
267            /(fse3t(:,:,jk)*rjjss)
268      END DO
269
270      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
271      CALL lbc_lnk( ironsed , 'T', 1. )
272      !!----------------------------------------------------------------------
273      !!
274      !! Initialize biological variables
275      !!
276      !!----------------------------------------------------------------------
277      !! Set biological ratios
278      !! ---------------------
279
280      rno3   = (16.+2.)/122.
281      po4r   = 1./122.
282      o2nit  = 32./122.
283      rdenit = 97.6/16.
284      o2ut   = 140./122.
285
286      !!----------------------------------------------------------------------
287      !!
288      !! Initialize chemical variables
289      !!
290      !!----------------------------------------------------------------------
291
292      !! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
293      !! ----------------------------------------------------------
294
295      atcox = 0.20946
296
297      !! Set lower/upper limits for temperature and salinity
298      !! ---------------------------------------------------
299
300      salchl = 1./1.80655
301      calcon = 1.03E-2
302
303      !! Set coefficients for apparent solubility equilibrium
304      !!   of calcite (Ingle, 1800, eq. 6)
305      !! ----------------------------------------------------
306
307      akcc1 = -34.452
308      akcc2 = -39.866
309      akcc3 = 110.21
310      akcc4 = -7.5752E-6
311
312
313      !! Set coefficients for seawater pressure correction
314      !! -------------------------------------------------
315
316      devk1(1) = -25.5
317      devk2(1) = 0.1271
318      devk3(1) = 0.
319      devk4(1) = -3.08E-3
320      devk5(1) = 0.0877E-3
321
322      devk1(2) = -15.82
323      devk2(2) = -0.0219
324      devk3(2) = 0.
325      devk4(2) = 1.13E-3
326      devk5(2) = -0.1475E-3
327
328      devk1(3) = -29.48
329      devk2(3) = 0.1622
330      devk3(3) = 2.608E-3
331      devk4(3) = -2.84E-3
332      devk5(3) = 0.
333
334      devk1(4) = -25.60
335      devk2(4) = 0.2324
336      devk3(4) = -3.6246E-3
337      devk4(4) = -5.13E-3
338      devk5(4) = 0.0794E-3
339
340      devkst = 0.23
341      devks  = 35.4
342
343      !! Set universal gas constants
344      !! ---------------------------
345
346      rgas = 83.143
347      oxyco = 1./22.4144
348
349      !! Set boron constants
350      !! -------------------
351
352      bor1 = 0.00023
353      bor2 = 1./10.82
354
355      !! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
356      !! -----------------------------------------------------------------
357
358      c00 = -60.2409
359      c01 = 93.4517
360      c02 = 23.3585
361      c03 = 0.023517
362      c04 = -0.023656
363      c05 = 0.0047036
364
365      ca0 = -162.8301
366      ca1 = 218.2968
367      ca2 = 90.9241
368      ca3 = -1.47696
369      ca4 = 0.025695
370      ca5 = -0.025225
371      ca6 = 0.0049867
372
373      !! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
374      !! ---------------------------------------------------------------------
375
376      c10 = -3670.7
377      c11 =  62.008
378      c12 = -9.7944
379      c13 = 0.0118
380      c14 = -0.000116
381
382      !! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
383      !! ---------------------------------------------------------------------
384
385      c20 = -1394.7
386      c21 = -4.777
387      c22 = 0.0184
388      c23 = -0.000118
389
390      !! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
391      !! ------------------------------------------------------------------
392
393      cb0  = -8966.90
394      cb1  = -2890.53
395      cb2  = -77.942
396      cb3  = 1.728
397      cb4  = -0.0996
398      cb5  = 148.0248
399      cb6  = 137.1942
400      cb7  = 1.62142
401      cb8  = -24.4344
402      cb9  = -25.085
403      cb10 = -0.2474
404      cb11 = 0.053105
405
406      !! Set coeff. for dissoc. of water (Dickson and Riley, 1979,
407      !!   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415
408      !!   after pers. commun. to B. Bacastow, 1988)
409      !! ---------------------------------------------------------
410
411      cw0 = -13847.26
412      cw1 = 148.9652
413      cw2 = -23.6521
414      cw3 = 118.67
415      cw4 = -5.977
416      cw5 = 1.0495
417      cw6 = -0.01615
418
419      !
420      ! Set coeff. for dissoc. of phosphate (Millero (1974)
421      ! ---------------------------------------------------
422      !
423      cp10 = 115.525
424      cp11 = -4576.752
425      cp12 = -18.453
426      cp13 = -106.736
427      cp14 = 0.69171
428      cp15 = -0.65643
429      cp16 = -0.01844
430
431      cp20 = 172.0883
432      cp21 = -8814.715
433      cp22 = -27.927
434      cp23 = -160.340
435      cp24 = 1.3566
436      cp25 = 0.37335
437      cp26 = -0.05778
438
439
440      cp30 = -18.141
441      cp31 = -3070.75
442      cp32 = 17.27039
443      cp33 = 2.81197
444      cp34 = -44.99486
445      cp35 = -0.09984
446      !
447      ! Set coeff. for dissoc. of phosphate (Millero (1974)
448      ! ---------------------------------------------------
449      !
450      cs10 = 117.385
451      cs11 = -8904.2
452      cs12 = -19.334
453      cs13 = -458.79
454      cs14 =  3.5913
455      cs15 = 188.74
456      cs16 = -1.5998
457      cs17 = -12.1652
458      cs18 = 0.07871
459      cs19 = -0.001005
460
461      !! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
462      !! ----------------------------------------------------------------
463
464      ox0 = -58.3877
465      ox1 = 85.8079
466      ox2 = 23.8439
467      ox3 = -0.034892
468      ox4 = 0.015568
469      ox5 = -0.0019387
470
471      !!  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
472      !!  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
473
474      open(49,file='kRGB61.txt',form='formatted')
475      do ichl=1,61
476         READ(49,*) xtoto,(xkrgb(iband,ichl),iband = 1,3)
477      end do
478      close(49)
479
480#if defined key_off_degrad
481
482      !! Read volume for degraded regions (DEGINIT)
483      !! ------------------------------------------
484
485#    if defined key_vpp
486      CALL READ3S(902,facvol,jpi,jpj,jpk)
487#    else
488      READ (902) facvol
489#    endif
490#endif
491
492
493      !!  Call p4zche to initialize the chemical constants
494      !!  ------------------------------------------------
495
496      CALL p4zche
497      !!
498      !!  Initialize a counter for the computation of chemistry
499      !!
500      ndayflxtr=0
501
502      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
503
504   END SUBROUTINE trc_ini
Note: See TracBrowser for help on using the repository browser.