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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/PISCES – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 13 years ago

update licence of all NEMO files...

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