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

source: branches/DEV_r1784_mid_year_merge_2010/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 1970

Last change on this file since 1970 was 1970, checked in by acc, 14 years ago

ticket #684 step 5: Add in changes from the trunk between revisions 1821 and 1879.

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