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 tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/PISCES – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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