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

source: trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 1180

Last change on this file since 1180 was 1180, checked in by cetlod, 16 years ago

update PISCES modules to couple with the sediment model, see ticket:249
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 23.3 KB
Line 
1MODULE p4zsed
2   !!======================================================================
3   !!                         ***  MODULE p4sed  ***
4   !! TOP :   PISCES Compute loss of organic matter in the sediments
5   !!======================================================================
6   !! History :   1.0  !  2004-03 (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!----------------------------------------------------------------------
9#if defined key_pisces
10   !!----------------------------------------------------------------------
11   !!   'key_pisces'                                       PISCES bio-model
12   !!----------------------------------------------------------------------
13   !!   p4z_sed        :  Compute loss of organic matter in the sediments
14   !!   p4z_sbc        :  Read and interpolate time-varying nutrients fluxes
15   !!   p4z_sed_init   :  Initialization of p4z_sed
16   !!----------------------------------------------------------------------
17   USE trc
18   USE oce_trc         !
19   USE sms_pisces
20   USE lib_mpp
21   USE prtctl_trc
22   USE p4zbio
23   USE p4zint
24   USE p4zopt
25   USE p4zsink
26   USE p4zrem
27   USE p4zlim
28   USE lbclnk
29   USE iom
30
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   p4z_sed   
36
37   !! * Shared module variables
38   LOGICAL, PUBLIC ::    &
39     bdustfer  = .FALSE.      ,  &  !:
40     briver    = .FALSE.      ,  &  !:
41     bndepo    = .FALSE.      ,  &  !:
42     bsedinput = .FALSE.            !:
43
44   REAL(wp), PUBLIC ::   &
45     sedfeinput = 1.E-9_wp   ,  &  !:
46     dustsolub  = 0.014_wp         !:
47
48   !! * Module variables
49   INTEGER ::          &
50      numdust,          &  ! logical unit for surface fluxes data
51      nflx1 , nflx2,   &  !  first and second record used
52      nflx11, nflx12      ! ???
53   REAL(wp), DIMENSION(jpi,jpj,2) ::    &  !:
54     dustmo                                !: 2 consecutive set of dust fields
55   REAL(wp), DIMENSION(jpi,jpj)   ::    &
56     rivinp, cotdep, nitdep, dust
57   REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   &
58     ironsed
59   REAL(wp) :: sumdepsi, rivalkinput, rivpo4input, nitdepinput
60
61   !!* Substitution
62#  include "domzgr_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
65   !! $Header:$
66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68
69CONTAINS
70
71   SUBROUTINE p4z_sed(kt, jnt)
72      !!---------------------------------------------------------------------
73      !!                     ***  ROUTINE p4z_sed  ***
74      !!
75      !! ** Purpose :   Compute loss of organic matter in the sediments. This
76      !!              is by no way a sediment model. The loss is simply
77      !!              computed to balance the inout from rivers and dust
78      !!
79      !! ** Method  : - ???
80      !!---------------------------------------------------------------------
81      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
82      INTEGER  ::   ji, jj, jk
83      INTEGER  ::   ikt
84#if ! defined key_sed
85      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
86#endif
87      REAL(wp) ::   zconctmp , zdenitot  , znitrpottot
88      REAL(wp) ::   zlim, zconctmp2, zstep, zfact
89      REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
90      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
91      CHARACTER (len=25) :: charout
92      !!---------------------------------------------------------------------
93
94
95      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_sed_init      ! Initialization (first time-step only)
96
97      IF( (jnt == 1) .and. (bdustfer) )  CALL p4z_sbc( kt )
98
99      zstep = rfact2 / rjjss      ! Time step duration for the biology
100
101      zirondep(:,:,:) = 0.e0          ! Initialisation of variables used to compute deposition
102      zsidep  (:,:)   = 0.e0
103
104      ! Iron and Si deposition at the surface
105      ! -------------------------------------
106
107      DO jj = 1, jpj
108         DO ji = 1, jpi
109            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmoss ) + 3.e-10 / raass )   &
110               &             * rfact2 / fse3t(ji,jj,1)
111            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / ( fse3t(ji,jj,1) * 28.1 * rmoss )
112         END DO
113      END DO
114
115      ! Iron solubilization of particles in the water column
116      ! ----------------------------------------------------
117
118      DO jk = 2, jpkm1
119         zirondep(:,:,jk) = dust(:,:) / ( 10. * 55.85 * rmoss ) * rfact2 * 1.e-4
120      END DO
121
122      ! Add the external input of nutrients, carbon and alkalinity
123      ! ----------------------------------------------------------
124
125      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2 
126      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2
127      trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivinp(:,:) * 3.e-5 * rfact2
128      trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep (:,:) + cotdep(:,:)   * rfact2 / 6.
129      trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivinp(:,:) * 2.631 * rfact2
130      trn(:,:,1,jptal) = trn(:,:,1,jptal) + (cotdep(:,:) - rno3*(rivinp(:,:) +  nitdep(:,:) ) ) * rfact2
131
132
133      ! Add the external input of iron which is 3D distributed
134      ! (dust, river and sediment mobilization)
135      ! ------------------------------------------------------
136
137      DO jk = 1, jpkm1
138         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer)   &
139         &       + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
140      END DO
141
142
143#if ! defined key_sed
144
145      ! Initialisation of variables used to compute Sinking Speed
146      ! ---------------------------------------------------------
147
148      zsumsedsi  = 0.e0
149      zsumsedpo4 = 0.e0
150      zsumsedcal = 0.e0
151
152      ! Loss of biogenic silicon, Caco3 organic carbon in the sediments.
153      ! First, the total loss is computed.
154      ! The factor for calcite comes from the alkalinity effect
155      ! -------------------------------------------------------------
156
157      DO jj = 1, jpj
158         DO ji = 1, jpi
159            ikt = MAX( mbathy(ji,jj)-1, 1 )
160            zfact = e1t(ji,jj) * e2t(ji,jj) / rjjss * tmask_i(ji,jj)
161
162# if defined key_kriest
163            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
164            zsumsedpo4 = zsumsedpo4 + zfact *  trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
165# else
166            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt)
167            zsumsedpo4 = zsumsedpo4 + zfact *( trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt)   &
168               &       + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )
169# endif
170
171            zsumsedcal = zsumsedcal + zfact *  trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0
172
173         END DO
174      END DO
175
176      IF( lk_mpp ) THEN
177         CALL mpp_sum( zsumsedsi  )   ! sums over the global domain
178         CALL mpp_sum( zsumsedcal )   ! sums over the global domain
179         CALL mpp_sum( zsumsedpo4 )   ! sums over the global domain
180      ENDIF
181
182#endif
183
184      ! Then this loss is scaled at each bottom grid cell for
185      ! equilibrating the total budget of silica in the ocean.
186      ! Thus, the amount of silica lost in the sediments equal
187      ! the supply at the surface (dust+rivers)
188      ! ------------------------------------------------------
189
190      DO jj = 1, jpj
191         DO ji = 1, jpi
192            ikt = MAX( mbathy(ji,jj) - 1, 1 )
193            zconctmp = trn(ji,jj,ikt,jpdsi) * zstep / fse3t(ji,jj,ikt)   &
194# if ! defined key_kriest
195     &             * wscal (ji,jj,ikt)
196# else
197     &             * wsbio4(ji,jj,ikt)
198# endif
199
200            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zconctmp
201
202#if ! defined key_sed
203            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zconctmp   &
204            &      * ( 1.- ( sumdepsi + rivalkinput / raass / 6. ) / zsumsedsi )
205#endif
206         END DO
207      END DO
208
209      DO jj = 1, jpj
210         DO ji = 1, jpi
211            ikt = MAX( mbathy(ji,jj) - 1, 1 )
212            zconctmp = trn(ji,jj,ikt,jpcal) * wscal(ji,jj,ikt) * zstep / fse3t(ji,jj,ikt)
213            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zconctmp
214
215#if ! defined key_sed
216            trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + zconctmp   &
217               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal ) * 2.e0
218            trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zconctmp   &
219               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal )
220#endif
221         END DO
222      END DO
223
224      DO jj = 1, jpj
225         DO ji = 1, jpi
226            ikt = MAX( mbathy(ji,jj) - 1, 1 )
227            zfact = zstep / fse3t(ji,jj,ikt)
228
229# if ! defined key_kriest
230            zconctmp  = trn(ji,jj,ikt,jpgoc)
231            zconctmp2 = trn(ji,jj,ikt,jppoc)
232            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - zconctmp  * wsbio4(ji,jj,ikt) * zfact
233            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
234#if ! defined key_sed
235            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
236            &      + ( zconctmp  * wsbio4(ji,jj,ikt) + zconctmp2 * wsbio3(ji,jj,ikt) ) * zfact   &
237            &      * ( 1.- rivpo4input / (raass * zsumsedpo4 ) )
238#endif
239            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * wsbio4(ji,jj,ikt) * zfact
240            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
241
242# else
243
244            zconctmp  = trn(ji,jj,ikt,jpnum)
245            zconctmp2 = trn(ji,jj,ikt,jppoc)
246            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum)   &
247            &      - zconctmp  * wsbio4(ji,jj,ikt) * zfact
248            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc)   &
249            &      - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
250#if ! defined key_sed
251            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
252            &      + ( zconctmp2 * wsbio3(ji,jj,ikt) )   &
253            &      * zfact * ( 1.- rivpo4input / ( raass * zsumsedpo4 ) )
254#endif
255            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe)   &
256            &      - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
257
258# endif
259         END DO
260      END DO
261
262      ! Nitrogen fixation (simple parameterization). The total gain
263      ! from nitrogen fixation is scaled to balance the loss by
264      ! denitrification
265      ! -------------------------------------------------------------
266
267      zdenitot = 0.e0
268      DO jk = 1, jpkm1
269         DO jj = 1,jpj
270            DO ji = 1,jpi
271               zdenitot = zdenitot + denitr(ji,jj,jk) * rdenit * e1t(ji,jj) * e2t(ji,jj)   &
272               &        *fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) * xnegtr(ji,jj,jk)
273            END DO
274         END DO
275      END DO
276
277      IF( lk_mpp )   CALL mpp_sum( zdenitot )      ! sum over the global domain
278
279      ! Potential nitrogen fication dependant on temperature and iron
280      ! -------------------------------------------------------------
281
282!CDIR NOVERRCHK
283      DO jk = 1, jpk
284!CDIR NOVERRCHK
285         DO jj = 1, jpj
286!CDIR NOVERRCHK
287            DO ji = 1, jpi
288               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
289               IF( zlim <= 0.2 )   zlim = 0.01
290               znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) / rjjss )   &
291# if defined key_off_degrad
292               &                  * facvol(ji,jj,jk)   &
293# endif
294               &                  * zlim * rfact2 * trn(ji,jj,jk,jpfer)   &
295               &                  / ( conc3 + trn(ji,jj,jk,jpfer) ) * ( 1.- EXP( -etot(ji,jj,jk) / 50.) )
296            END DO
297         END DO
298      END DO
299
300      znitrpottot = 0.e0
301      DO jk = 1, jpkm1
302         DO jj = 1, jpj
303            DO ji = 1, jpi
304               znitrpottot = znitrpottot + znitrpot(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)   &
305                  &                                           * tmask(ji,jj,jk) * tmask_i(ji,jj) 
306            END DO
307         END DO
308      END DO
309
310      IF( lk_mpp )   CALL mpp_sum( znitrpottot )  ! sum over the global domain
311
312      ! Nitrogen change due to nitrogen fixation
313      ! ----------------------------------------
314
315      DO jk = 1, jpk
316         DO jj = 1, jpj
317            DO ji = 1, jpi
318# if ! defined key_c1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 )
319!!             zfact = znitrpot(ji,jj,jk) * zdenitot / znitrpottot
320               zfact = znitrpot(ji,jj,jk) * 1.e-7
321# else
322               zfact = znitrpot(ji,jj,jk) * 1.e-7
323# endif
324               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
325               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
326               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30./ 46.* zfact
327            END DO
328         END DO
329      END DO
330
331# if defined key_trc_diaadd
332      DO jj = 1,jpj
333         DO ji = 1,jpi
334            trc2d(ji,jj,jp_pcs0_2d + 11) = zirondep(ji,jj,1) * 1.e+3 * rfact2r * fse3t(ji,jj,1)
335            trc2d(ji,jj,jp_pcs0_2d + 12) = znitrpot(ji,jj,1) * 1.e-7 * fse3t(ji,jj,1) * 1.e+3 / rfact2
336         END DO
337      END DO
338# endif
339      !
340       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
341         WRITE(charout, FMT="('sed ')")
342         CALL prt_ctl_trc_info(charout)
343         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
344       ENDIF
345
346   END SUBROUTINE p4z_sed
347
348   SUBROUTINE p4z_sbc(kt)
349
350      !!----------------------------------------------------------------------
351      !!                  ***  ROUTINE p4z_sbc  ***
352      !!
353      !! ** Purpose :   Read and interpolate the external sources of
354      !!                nutrients
355      !!
356      !! ** Method  :   Read the files and interpolate the appropriate variables
357      !!
358      !! ** input   :   external netcdf files
359      !!
360      !!----------------------------------------------------------------------
361      !! * arguments
362      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
363
364      !! * Local declarations
365      INTEGER ::   &
366         imois, imois2,       &  ! temporary integers
367         i15  , iman             !    "          "
368      REAL(wp) ::   &
369         zxy                     !    "         "
370
371
372      !!---------------------------------------------------------------------
373
374      ! Initialization
375      ! --------------
376
377      i15 = nday / 16
378      iman  = INT( raamo )
379      imois = nmonth + i15 - 1
380      IF( imois == 0 ) imois = iman
381      imois2 = nmonth
382
383      ! 1. first call kt=nit000
384      ! -----------------------
385
386      IF( kt == nit000 ) THEN
387         ! initializations
388         nflx1  = 0
389         nflx11 = 0
390         ! open the file
391         IF(lwp) THEN
392            WRITE(numout,*) ' '
393            WRITE(numout,*) ' **** Routine p4z_sbc'
394         ENDIF
395         CALL iom_open ( 'dust.orca.nc', numdust )
396      ENDIF
397
398
399     ! Read monthly file
400      ! ----------------
401
402      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
403
404         ! Calendar computation
405
406         ! nflx1 number of the first file record used in the simulation
407         ! nflx2 number of the last  file record
408
409         nflx1 = imois
410         nflx2 = nflx1+1
411         nflx1 = MOD( nflx1, iman )
412         nflx2 = MOD( nflx2, iman )
413         IF( nflx1 == 0 )   nflx1 = iman
414         IF( nflx2 == 0 )   nflx2 = iman
415         IF(lwp) WRITE(numout,*) 'first record file used nflx1 ',nflx1
416         IF(lwp) WRITE(numout,*) 'last  record file used nflx2 ',nflx2
417
418         ! Read monthly fluxes data
419
420         ! humidity
421         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
422         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
423
424         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
425            WRITE(numout,*)
426            WRITE(numout,*) ' read clio flx ok'
427            WRITE(numout,*)
428               WRITE(numout,*)
429               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
430               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
431         ENDIF
432
433      ENDIF
434
435     ! 3. at every time step interpolation of fluxes
436      ! ---------------------------------------------
437
438      zxy = FLOAT( nday + 15 - 30 * i15 ) / 30
439      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
440
441      IF( kt == nitend ) CALL iom_close (numdust)
442
443   END SUBROUTINE p4z_sbc
444
445
446   SUBROUTINE p4z_sed_init
447
448      !!----------------------------------------------------------------------
449      !!                  ***  ROUTINE p4z_sed_init  ***
450      !!
451      !! ** Purpose :   Initialization of the external sources of nutrients
452      !!
453      !! ** Method  :   Read the files and compute the budget
454      !!      called at the first timestep (nittrc000)
455      !!
456      !! ** input   :   external netcdf files
457      !!
458      !!----------------------------------------------------------------------
459
460      INTEGER ::   ji, jj, jk, jm
461      INTEGER , PARAMETER ::   jpmois = 12, jpan = 1
462      INTEGER :: numriv, numbath, numdep
463
464
465      REAL(wp) ::   zcoef
466      REAL(wp) ::   expide, denitide,zmaskt
467      REAL(wp) , DIMENSION (jpi,jpj)     ::   riverdoc, river, ndepo
468      REAL(wp) , DIMENSION (jpi,jpj,jpk) ::   cmask
469      REAL(wp), DIMENSION(jpi,jpj,12)    ::   zdustmo
470
471      NAMELIST/nampissed/ bdustfer, briver, bndepo, bsedinput, sedfeinput, dustsolub
472
473
474      REWIND( numnat )                     ! read numnat
475      READ  ( numnat, nampissed )
476
477      IF(lwp) THEN
478         WRITE(numout,*) ' '
479         WRITE(numout,*) ' Namelist : nampissed '
480         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
481         WRITE(numout,*) '    Dust input from the atmosphere           bdustfer  = ', bdustfer
482         WRITE(numout,*) '    River input of nutrients                 briver    = ', briver
483         WRITE(numout,*) '    Atmospheric deposition of N              bndepo    = ', bndepo
484         WRITE(numout,*) '    Fe input from sediments                  bsedinput = ', bsedinput
485         WRITE(numout,*) '    Coastal release of Iron                  sedfeinput =', sedfeinput
486         WRITE(numout,*) '    Solubility of the dust                   dustsolub  =', dustsolub
487      ENDIF
488
489      ! Dust input from the atmosphere
490      ! ------------------------------
491      IF( bdustfer ) THEN
492         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
493         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
494         CALL iom_open ( 'dust.orca.nc', numdust )
495         DO jm = 1, jpmois
496            CALL iom_get( numdust, jpdom_data, 'dust', zdustmo(:,:,jm), jm )
497         END DO
498         CALL iom_close( numdust )
499      ELSE
500         zdustmo(:,:,:) = 0.e0
501         dust(:,:) = 0.0
502      ENDIF
503
504      ! Nutrient input from rivers
505      ! --------------------------
506      IF( briver ) THEN
507         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
508         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
509         CALL iom_open ( 'river.orca.nc', numriv )
510         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
511         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
512         CALL iom_close( numriv )
513      ELSE
514         river   (:,:) = 0.e0
515         riverdoc(:,:) = 0.e0
516      endif
517
518      ! Nutrient input from dust
519      ! ------------------------
520      IF( bndepo ) THEN
521         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
522         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
523         CALL iom_open ( 'ndeposition.orca.nc', numdep )
524         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
525         CALL iom_close( numdep )
526      ELSE
527         ndepo(:,:) = 0.e0
528      ENDIF
529
530      ! Coastal and island masks
531      ! ------------------------
532      IF( bsedinput ) THEN     
533         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron'
534         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
535         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
536         CALL iom_open ( 'bathy.orca.nc', numbath )
537         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
538         CALL iom_close( numbath )
539         !
540         DO jk = 1, 5
541            DO jj = 2, jpjm1
542               DO ji = 2, jpim1
543                  IF( tmask(ji,jj,jk) /= 0. ) THEN
544                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
545                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
546                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
547                  ENDIF
548               END DO
549            END DO
550         END DO
551         DO jk = 1, jpk
552            DO jj = 1, jpj
553               DO ji = 1, jpi
554                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
555                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
556                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
557               END DO
558            END DO
559         END DO
560      ELSE
561         cmask(:,:,:) = 0.e0
562      ENDIF
563
564      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
565
566
567      ! total atmospheric supply of Si
568      ! ------------------------------
569      sumdepsi = 0.e0
570      DO jm = 1, jpmois
571         DO jj = 2, jpjm1
572            DO ji = 2, jpim1
573               sumdepsi = sumdepsi + zdustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
574                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
575            END DO
576         END DO
577      END DO
578      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
579
580      ! N/P and Si releases due to coastal rivers
581      ! -----------------------------------------
582      DO jj = 1, jpj
583         DO ji = 1, jpi
584            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
585            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
586            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
587            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
588         END DO
589      END DO
590      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
591      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
592
593      rivpo4input=0.e0
594      rivalkinput=0.e0
595      nitdepinput=0.e0
596      DO jj = 2 , jpjm1
597         DO ji = 2, jpim1
598            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
599            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
600            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
601            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
602         END DO
603     END DO
604      IF( lk_mpp ) THEN
605         CALL mpp_sum( rivpo4input )  ! sum over the global domain
606         CALL mpp_sum( rivalkinput )  ! sum over the global domain
607         CALL mpp_sum( nitdepinput )  ! sum over the global domain
608      ENDIF
609
610
611      ! Coastal supply of iron
612      ! -------------------------
613      DO jk = 1, jpkm1
614         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
615      END DO
616      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
617
618
619   END SUBROUTINE p4z_sed_init
620
621
622
623#else
624   !!======================================================================
625   !!  Dummy module :                                   No PISCES bio-model
626   !!======================================================================
627CONTAINS
628   SUBROUTINE p4z_sed                         ! Empty routine
629   END SUBROUTINE p4z_sed
630#endif 
631
632   !!======================================================================
633END MODULE  p4zsed
Note: See TracBrowser for help on using the repository browser.