source: branches/CMIP5_IPSL/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 1830

Last change on this file since 1830 was 1830, checked in by cetlod, 11 years ago

Computation of additional diagnostics for PISCES model ( under CPP key key_diaar5 )

  • needed for AR5 outputs (vertical inventories, passive tracers at surface,… )
  • new output file with suffix dbio_T
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 24.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
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#if defined key_diaar5
340      ! nitrogen fixation rate in ocean ( vertically integrated )
341      zw2d(:,:) = 0.
342      DO jk = 1, jpkm1
343         zw2d(:,:) = zw2d(:,:) + znitrpot(:,:,jk) * 1.e-7 * zrfact2  * fse3t(:,:,jk) * tmask(:,:,jk)
344      ENDDO
345      IF( jnt == nrdttrc ) CALL iom_put( "INTNFIX" , zw2d )
346# endif
347# endif
348# endif
349      !
350       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
351         WRITE(charout, FMT="('sed ')")
352         CALL prt_ctl_trc_info(charout)
353         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
354       ENDIF
355
356   END SUBROUTINE p4z_sed
357
358   SUBROUTINE p4z_sbc(kt)
359
360      !!----------------------------------------------------------------------
361      !!                  ***  ROUTINE p4z_sbc  ***
362      !!
363      !! ** Purpose :   Read and interpolate the external sources of
364      !!                nutrients
365      !!
366      !! ** Method  :   Read the files and interpolate the appropriate variables
367      !!
368      !! ** input   :   external netcdf files
369      !!
370      !!----------------------------------------------------------------------
371      !! * arguments
372      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
373
374      !! * Local declarations
375      INTEGER ::   &
376         imois, imois2,       &  ! temporary integers
377         i15  , iman             !    "          "
378      REAL(wp) ::   &
379         zxy                     !    "         "
380
381
382      !!---------------------------------------------------------------------
383
384      ! Initialization
385      ! --------------
386
387      i15 = nday / 16
388      iman  = INT( raamo )
389      imois = nmonth + i15 - 1
390      IF( imois == 0 ) imois = iman
391      imois2 = nmonth
392
393      ! 1. first call kt=nit000
394      ! -----------------------
395
396      IF( kt == nit000 ) THEN
397         ! initializations
398         nflx1  = 0
399         nflx11 = 0
400         ! open the file
401         IF(lwp) THEN
402            WRITE(numout,*) ' '
403            WRITE(numout,*) ' **** Routine p4z_sbc'
404         ENDIF
405         CALL iom_open ( 'dust.orca.nc', numdust )
406      ENDIF
407
408
409     ! Read monthly file
410      ! ----------------
411
412      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
413
414         ! Calendar computation
415
416         ! nflx1 number of the first file record used in the simulation
417         ! nflx2 number of the last  file record
418
419         nflx1 = imois
420         nflx2 = nflx1+1
421         nflx1 = MOD( nflx1, iman )
422         nflx2 = MOD( nflx2, iman )
423         IF( nflx1 == 0 )   nflx1 = iman
424         IF( nflx2 == 0 )   nflx2 = iman
425         IF(lwp) WRITE(numout,*) 'first record file used nflx1 ',nflx1
426         IF(lwp) WRITE(numout,*) 'last  record file used nflx2 ',nflx2
427
428         ! Read monthly fluxes data
429
430         ! humidity
431         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
432         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
433
434         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
435            WRITE(numout,*)
436            WRITE(numout,*) ' read clio flx ok'
437            WRITE(numout,*)
438               WRITE(numout,*)
439               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
440               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
441         ENDIF
442
443      ENDIF
444
445     ! 3. at every time step interpolation of fluxes
446      ! ---------------------------------------------
447
448      zxy = FLOAT( nday + 15 - 30 * i15 ) / 30
449      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
450
451      IF( kt == nitend ) CALL iom_close (numdust)
452
453   END SUBROUTINE p4z_sbc
454
455
456   SUBROUTINE p4z_sed_init
457
458      !!----------------------------------------------------------------------
459      !!                  ***  ROUTINE p4z_sed_init  ***
460      !!
461      !! ** Purpose :   Initialization of the external sources of nutrients
462      !!
463      !! ** Method  :   Read the files and compute the budget
464      !!      called at the first timestep (nittrc000)
465      !!
466      !! ** input   :   external netcdf files
467      !!
468      !!----------------------------------------------------------------------
469
470      INTEGER ::   ji, jj, jk, jm
471      INTEGER , PARAMETER ::   jpmois = 12, jpan = 1
472      INTEGER :: numriv, numbath, numdep
473
474
475      REAL(wp) ::   zcoef
476      REAL(wp) ::   expide, denitide,zmaskt
477      REAL(wp) , DIMENSION (jpi,jpj)     ::   riverdoc, river, ndepo
478      REAL(wp) , DIMENSION (jpi,jpj,jpk) ::   cmask
479      REAL(wp) , DIMENSION(jpi,jpj,12)    ::   zdustmo
480
481      NAMELIST/nampissed/ ln_dustfer, ln_river, ln_ndepo, ln_sedinput, sedfeinput, dustsolub
482
483
484      REWIND( numnat )                     ! read numnat
485      READ  ( numnat, nampissed )
486
487      IF(lwp) THEN
488         WRITE(numout,*) ' '
489         WRITE(numout,*) ' Namelist : nampissed '
490         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
491         WRITE(numout,*) '    Dust input from the atmosphere           ln_dustfer  = ', ln_dustfer
492         WRITE(numout,*) '    River input of nutrients                 ln_river    = ', ln_river
493         WRITE(numout,*) '    Atmospheric deposition of N              ln_ndepo    = ', ln_ndepo
494         WRITE(numout,*) '    Fe input from sediments                  ln_sedinput = ', ln_sedinput
495         WRITE(numout,*) '    Coastal release of Iron                  sedfeinput  =', sedfeinput
496         WRITE(numout,*) '    Solubility of the dust                   dustsolub   =', dustsolub
497      ENDIF
498
499      ! Dust input from the atmosphere
500      ! ------------------------------
501      IF( ln_dustfer ) THEN
502         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
503         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
504         CALL iom_open ( 'dust.orca.nc', numdust )
505         DO jm = 1, jpmois
506            CALL iom_get( numdust, jpdom_data, 'dust', zdustmo(:,:,jm), jm )
507         END DO
508         CALL iom_close( numdust )
509      ELSE
510         zdustmo(:,:,:) = 0.e0
511         dust(:,:) = 0.0
512      ENDIF
513
514      ! Nutrient input from rivers
515      ! --------------------------
516      IF( ln_river ) THEN
517         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
518         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
519         CALL iom_open ( 'river.orca.nc', numriv )
520         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
521         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
522         CALL iom_close( numriv )
523      ELSE
524         river   (:,:) = 0.e0
525         riverdoc(:,:) = 0.e0
526      endif
527
528      ! Nutrient input from dust
529      ! ------------------------
530      IF( ln_ndepo ) THEN
531         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
532         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
533         CALL iom_open ( 'ndeposition.orca.nc', numdep )
534         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
535         CALL iom_close( numdep )
536      ELSE
537         ndepo(:,:) = 0.e0
538      ENDIF
539
540      ! Coastal and island masks
541      ! ------------------------
542      IF( ln_sedinput ) THEN     
543         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron'
544         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
545         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
546         CALL iom_open ( 'bathy.orca.nc', numbath )
547         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
548         CALL iom_close( numbath )
549         !
550         DO jk = 1, 5
551            DO jj = 2, jpjm1
552               DO ji = fs_2, fs_jpim1
553                  IF( tmask(ji,jj,jk) /= 0. ) THEN
554                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
555                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
556                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
557                  ENDIF
558               END DO
559            END DO
560         END DO
561         DO jk = 1, jpk
562            DO jj = 1, jpj
563               DO ji = 1, jpi
564                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
565                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
566                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
567               END DO
568            END DO
569         END DO
570      ELSE
571         cmask(:,:,:) = 0.e0
572      ENDIF
573
574      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
575
576
577      ! Number of seconds per year and per month
578      ryyss = nyear_len(1) * rday
579      rmtss = ryyss / raamo
580
581      ! total atmospheric supply of Si
582      ! ------------------------------
583      sumdepsi = 0.e0
584      DO jm = 1, jpmois
585         DO jj = 2, jpjm1
586            DO ji = fs_2, fs_jpim1
587               sumdepsi = sumdepsi + zdustmo(ji,jj,jm) / (12.*rmtss) * 8.8        &
588                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
589            END DO
590         END DO
591      END DO
592      IF( lk_mpp )  CALL mpp_sum( sumdepsi )  ! sum over the global domain
593
594      ! N/P and Si releases due to coastal rivers
595      ! -----------------------------------------
596      DO jj = 1, jpj
597         DO ji = 1, jpi
598            zcoef = ryyss * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1) * tmask_i(ji,jj)
599            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn )
600            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn )
601            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*ryyss*fse3t(ji,jj,1) + rtrn )
602         END DO
603      END DO
604      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
605      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
606
607      rivpo4input = 0.e0
608      rivalkinput = 0.e0
609      nitdepinput = 0.e0
610      DO jj = 2 , jpjm1
611         DO ji = fs_2, fs_jpim1
612            zcoef = cvol(ji,jj,1) * ryyss
613            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
614            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
615            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
616         END DO
617     END DO
618      IF( lk_mpp ) THEN
619         CALL mpp_sum( rivpo4input )  ! sum over the global domain
620         CALL mpp_sum( rivalkinput )  ! sum over the global domain
621         CALL mpp_sum( nitdepinput )  ! sum over the global domain
622      ENDIF
623
624
625      ! Coastal supply of iron
626      ! -------------------------
627      DO jk = 1, jpkm1
628         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
629      END DO
630      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
631
632
633   END SUBROUTINE p4z_sed_init
634
635#else
636   !!======================================================================
637   !!  Dummy module :                                   No PISCES bio-model
638   !!======================================================================
639CONTAINS
640   SUBROUTINE p4z_sed                         ! Empty routine
641   END SUBROUTINE p4z_sed
642#endif 
643
644   !!======================================================================
645END MODULE  p4zsed
Note: See TracBrowser for help on using the repository browser.