Ticket #1390: p4zsed.F90

File p4zsed.F90, 30.4 KB (added by lambertn, 6 years ago)
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   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) USE of fldread
9   !!----------------------------------------------------------------------
10#if defined key_pisces
11   !!----------------------------------------------------------------------
12   !!   'key_pisces'                                       PISCES bio-model
13   !!----------------------------------------------------------------------
14   !!   p4z_sed        :  Compute loss of organic matter in the sediments
15   !!   p4z_sbc        :  Read and interpolate time-varying nutrients fluxes
16   !!   p4z_sed_init   :  Initialization of p4z_sed
17   !!----------------------------------------------------------------------
18   USE oce_trc         !  shared variables between ocean and passive tracers
19   USE trc             !  passive tracers common variables
20   USE sms_pisces      !  PISCES Source Minus Sink variables
21   USE p4zsink         !  vertical flux of particulate matter due to sinking
22   USE p4zopt          !  optical model
23   USE p4zlim          !  Co-limitations of differents nutrients
24   USE p4zrem          !  Remineralisation of organic matter
25   USE p4zint          !  interpolation and computation of various fields
26   USE iom             !  I/O manager
27   USE fldread         !  time interpolation
28   USE prtctl_trc      !  print control for debugging
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   p4z_sed   
34   PUBLIC   p4z_sed_init   
35   PUBLIC   p4z_sed_alloc
36
37   !! * Shared module variables
38   LOGICAL  :: ln_dust     = .FALSE.    !: boolean for dust input from the atmosphere
39   LOGICAL  :: ln_river    = .FALSE.    !: boolean for river input of nutrients
40   LOGICAL  :: ln_ndepo    = .FALSE.    !: boolean for atmospheric deposition of N
41   LOGICAL  :: ln_ironsed  = .FALSE.    !: boolean for Fe input from sediments
42
43   REAL(wp) :: sedfeinput  = 1.E-9_wp   !: Coastal release of Iron
44   REAL(wp) :: dustsolub   = 0.014_wp   !: Solubility of the dust
45   REAL(wp) :: wdust       = 2.0_wp     !: Sinking speed of the dust
46   REAL(wp) :: nitrfix     = 1E-7_wp    !: Nitrogen fixation rate   
47   REAL(wp) :: diazolight  = 50._wp     !: Nitrogen fixation sensitivty to light
48   REAL(wp) :: concfediaz  = 1.E-10_wp  !: Fe half-saturation Cste for diazotrophs
49
50
51   !! * Module variables
52   REAL(wp) :: ryyss                  !: number of seconds per year
53   REAL(wp) :: r1_ryyss                 !: inverse of ryyss
54   REAL(wp) :: rmtss                  !: number of seconds per month
55   REAL(wp) :: r1_rday                  !: inverse of rday
56   LOGICAL  :: ll_sbc
57
58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdic  ! structure of input riverdic
60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_riverdoc  ! structure of input riverdoc
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ndepo     ! structure of input nitrogen deposition
62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
63
64   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file
65   INTEGER  :: ntimes_dust, ntimes_riv, ntimes_ndep       ! number of time steps in a file
66
67   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: dust      !: dust fields
68   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivinp, cotdep    !: river input fields
69   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: nitdep    !: atmospheric N deposition
70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ironsed   !: Coastal supply of iron
71
72   REAL(wp) :: sumdepsi, rivalkinput, rivpo4input, nitdepinput
73
74   !!* Substitution
75#  include "top_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
78   !! $Header:$
79   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81
82CONTAINS
83
84
85   SUBROUTINE p4z_sed( kt, jnt )
86      !!---------------------------------------------------------------------
87      !!                     ***  ROUTINE p4z_sed  ***
88      !!
89      !! ** Purpose :   Compute loss of organic matter in the sediments. This
90      !!              is by no way a sediment model. The loss is simply
91      !!              computed to balance the inout from rivers and dust
92      !!
93      !! ** Method  : - ???
94      !!---------------------------------------------------------------------
95      !
96      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
97      INTEGER  ::   ji, jj, jk, ikt
98#if ! defined key_sed
99      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
100      REAL(wp) ::   zrivalk, zrivsil, zrivpo4
101#endif
102      REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact, zfactcal
103      REAL(wp) ::   zsiloss, zcaloss, zwsbio3, zwsbio4, zwscal, zdep
104      CHARACTER (len=25) :: charout
105      REAL(wp), POINTER, DIMENSION(:,:  ) :: zsidep, zwork1, zwork2, zwork3
106      REAL(wp), POINTER, DIMENSION(:,:,:) :: znitrpot, zirondep
107      !!---------------------------------------------------------------------
108      !
109      IF( nn_timing == 1 )  CALL timing_start('p4z_sed')
110      !
111      ! Allocate temporary workspace
112      CALL wrk_alloc( jpi, jpj,      zsidep, zwork1, zwork2, zwork3 )
113      CALL wrk_alloc( jpi, jpj, jpk, znitrpot, zirondep             )
114
115      IF( jnt == 1 .AND. ll_sbc ) CALL p4z_sbc( kt )
116
117      zirondep(:,:,:) = 0.e0          ! Initialisation of variables USEd to compute deposition
118      zsidep  (:,:)   = 0.e0
119
120      ! Iron and Si deposition at the surface
121      ! -------------------------------------
122      DO jj = 1, jpj
123         DO ji = 1, jpi
124            zdep  = rfact2 / fse3t(ji,jj,1)
125            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 * r1_ryyss ) * zdep
126            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * zdep / ( 28.1 * rmtss )
127         END DO
128      END DO
129
130      ! Iron solubilization of particles in the water column
131      ! ----------------------------------------------------
132      DO jk = 2, jpkm1
133         zirondep(:,:,jk) = dust(:,:) / ( wdust * 55.85 * rmtss ) * rfact2 * 1.e-4 * EXP( -fsdept(:,:,jk) / 1000. )
134      END DO
135
136      ! Add the external input of nutrients, carbon and alkalinity
137      ! ----------------------------------------------------------
138      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2 
139      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2
140      trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivinp(:,:) * 3.e-5 * rfact2
141      trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep (:,:) + cotdep(:,:)   * rfact2 / 6.
142      trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivinp(:,:) * 2.631 * rfact2
143      trn(:,:,1,jptal) = trn(:,:,1,jptal) + (cotdep(:,:) - rno3*(rivinp(:,:) +  nitdep(:,:) ) ) * rfact2
144
145
146      ! Add the external input of iron which is 3D distributed
147      ! (dust, river and sediment mobilization)
148      ! ------------------------------------------------------
149      DO jk = 1, jpkm1
150         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
151      END DO
152
153#if ! defined key_sed
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 = mbkt(ji,jj) 
161# if defined key_kriest
162            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * wscal (ji,jj,ikt)
163            zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
164# else
165            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * wsbio4(ji,jj,ikt)
166            zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) 
167# endif
168            ! For calcite, burial efficiency is made a function of saturation
169            zfactcal      = MIN( excess(ji,jj,ikt), 0.2 )
170            zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
171            zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0 * zfactcal
172         END DO
173      END DO
174      zsumsedsi  = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * r1_rday
175      zsumsedpo4 = glob_sum( zwork2(:,:) * e1e2t(:,:) ) * r1_rday
176      zsumsedcal = glob_sum( zwork3(:,:) * e1e2t(:,:) ) * r1_rday
177#endif
178
179      ! THEN this loss is scaled at each bottom grid cell for
180      ! equilibrating the total budget of silica in the ocean.
181      ! Thus, the amount of silica lost in the sediments equal
182      ! the supply at the surface (dust+rivers)
183      ! ------------------------------------------------------
184#if ! defined key_sed
185      zrivsil =  1._wp - ( sumdepsi + rivalkinput * r1_ryyss / 6. ) / zsumsedsi 
186      zrivpo4 =  1._wp - ( rivpo4input * r1_ryyss ) / zsumsedpo4 
187#endif
188
189      DO jj = 1, jpj
190         DO ji = 1, jpi
191            ikt  = mbkt(ji,jj)
192            zdep = xstep / fse3t(ji,jj,ikt)
193            zwsbio4 = wsbio4(ji,jj,ikt) * zdep
194            zwscal  = wscal (ji,jj,ikt) * zdep
195# if defined key_kriest
196            zsiloss = trn(ji,jj,ikt,jpgsi) * zwsbio4
197# else
198            zsiloss = trn(ji,jj,ikt,jpgsi) * zwscal
199# endif
200            zcaloss = trn(ji,jj,ikt,jpcal) * zwscal
201            !
202            trn(ji,jj,ikt,jpgsi) = trn(ji,jj,ikt,jpgsi) - zsiloss
203            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zcaloss
204#if ! defined key_sed
205            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
206            zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
207            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
208            zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / zsumsedcal 
209            trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
210            trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + zcaloss * zrivalk
211#endif
212         END DO
213      END DO
214
215      DO jj = 1, jpj
216         DO ji = 1, jpi
217            ikt  = mbkt(ji,jj)
218            zdep = xstep / fse3t(ji,jj,ikt)
219            zwsbio4 = wsbio4(ji,jj,ikt) * zdep
220            zwsbio3 = wsbio3(ji,jj,ikt) * zdep
221# if ! defined key_kriest
222            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zwsbio4
223            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3
224            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zwsbio4
225            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3
226#if ! defined key_sed
227            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) &
228               &               + ( trn(ji,jj,ikt,jpgoc) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4
229#endif
230
231# else
232            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zwsbio4
233            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3
234            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3
235#if ! defined key_sed
236            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) &
237               &               + ( trn(ji,jj,ikt,jpnum) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 ) * zrivpo4
238#endif
239
240# endif
241         END DO
242      END DO
243
244
245      ! Nitrogen fixation (simple parameterization). The total gain
246      ! from nitrogen fixation is scaled to balance the loss by
247      ! denitrification
248      ! -------------------------------------------------------------
249
250      zdenitot = glob_sum(  ( denitr(:,:,:) * rdenit + denitnh4(:,:,:) * rdenita ) * cvol(:,:,:) ) 
251
252      ! Potential nitrogen fixation dependant on temperature and iron
253      ! -------------------------------------------------------------
254
255!CDIR NOVERRCHK
256      DO jk = 1, jpk
257!CDIR NOVERRCHK
258         DO jj = 1, jpj
259!CDIR NOVERRCHK
260            DO ji = 1, jpi
261               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
262               IF( zlim <= 0.2 )   zlim = 0.01
263#if defined key_degrad
264               zfact = zlim * rfact2 * facvol(ji,jj,jk)
265#else
266               zfact = zlim * rfact2 
267#endif
268               znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
269                 &                 *  zfact * trn(ji,jj,jk,jpfer) / ( concfediaz + trn(ji,jj,jk,jpfer) ) &
270                 &                 * ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) )
271            END DO
272         END DO
273      END DO
274
275      znitrpottot = glob_sum( znitrpot(:,:,:) * cvol(:,:,:) )
276
277      ! Nitrogen change due to nitrogen fixation
278      ! ----------------------------------------
279      DO jk = 1, jpk
280         DO jj = 1, jpj
281            DO ji = 1, jpi
282               zfact = znitrpot(ji,jj,jk) * nitrfix
283               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
284               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3 * zfact
285               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
286               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30. / 46. * zfact
287           !    trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + zfact
288           END DO
289         END DO
290      END DO
291      !
292      IF( ln_diatrc ) THEN
293         zfact = 1.e+3 * rfact2r
294         IF( lk_iomput ) THEN
295            zwork1(:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1) 
296            zwork2(:,:)  =    znitrpot(:,:,1) * nitrfix                   * zfact * fse3t(:,:,1) * tmask(:,:,1)
297            IF( jnt == nrdttrc ) THEN
298               CALL iom_put( "Irondep", zwork1  )  ! surface downward net flux of iron
299               CALL iom_put( "Nfix"   , zwork2 )  ! nitrogen fixation at surface
300            ENDIF
301         ELSE
302            trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)           * zfact * fse3t(:,:,1) * tmask(:,:,1)
303            trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1)
304         ENDIF
305      ENDIF
306      !
307      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging)
308         WRITE(charout, fmt="('sed ')")
309         CALL prt_ctl_trc_info(charout)
310         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
311      ENDIF
312      !
313      CALL wrk_dealloc( jpi, jpj,      zsidep, zwork1, zwork2, zwork3 )
314      CALL wrk_dealloc( jpi, jpj, jpk, znitrpot, zirondep             )
315      !
316      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed')
317      !
318   END SUBROUTINE p4z_sed
319
320   SUBROUTINE p4z_sbc( kt )
321      !!----------------------------------------------------------------------
322      !!                  ***  routine p4z_sbc  ***
323      !!
324      !! ** purpose :   read and interpolate the external sources of
325      !!                nutrients
326      !!
327      !! ** method  :   read the files and interpolate the appropriate variables
328      !!
329      !! ** input   :   external netcdf files
330      !!
331      !!----------------------------------------------------------------------
332      !! * arguments
333      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
334
335      !! * local declarations
336      INTEGER  :: ji,jj 
337      REAL(wp) :: zcoef
338      !!---------------------------------------------------------------------
339      !
340      IF( nn_timing == 1 )  CALL timing_start('p4z_sbc')
341      !
342      ! Compute dust at nit000 or only if there is more than 1 time record in dust file
343      IF( ln_dust ) THEN
344         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
345            CALL fld_read( kt, 1, sf_dust )
346            dust(:,:) = sf_dust(1)%fnow(:,:,1)
347         ENDIF
348      ENDIF
349
350      ! N/P and Si releases due to coastal rivers
351      ! Compute river at nit000 or only if there is more than 1 time record in river file
352      ! -----------------------------------------
353      IF( ln_river ) THEN
354         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN
355            CALL fld_read( kt, 1, sf_riverdic )
356            CALL fld_read( kt, 1, sf_riverdoc )
357            DO jj = 1, jpj
358               DO ji = 1, jpi
359                  zcoef = ryyss * cvol(ji,jj,1) 
360                  cotdep(ji,jj) =   sf_riverdic(1)%fnow(ji,jj,1)                                  * 1E9 / ( 12. * zcoef + rtrn )
361                  rivinp(ji,jj) = ( sf_riverdic(1)%fnow(ji,jj,1) + sf_riverdoc(1)%fnow(ji,jj,1) ) * 1E9 / ( 31.6* zcoef + rtrn )
362               END DO
363            END DO
364         ENDIF
365      ENDIF
366
367      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file
368      IF( ln_ndepo ) THEN
369         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
370            CALL fld_read( kt, 1, sf_ndepo )
371            DO jj = 1, jpj
372               DO ji = 1, jpi
373                  nitdep(ji,jj) = 7.6 * sf_ndepo(1)%fnow(ji,jj,1) / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn )
374               END DO
375            END DO
376         ENDIF
377      ENDIF
378      !
379      IF( nn_timing == 1 )  CALL timing_stop('p4z_sbc')
380      !
381   END SUBROUTINE p4z_sbc
382
383   SUBROUTINE p4z_sed_init
384
385      !!----------------------------------------------------------------------
386      !!                  ***  routine p4z_sed_init  ***
387      !!
388      !! ** purpose :   initialization of the external sources of nutrients
389      !!
390      !! ** method  :   read the files and compute the budget
391      !!                called at the first timestep (nittrc000)
392      !!
393      !! ** input   :   external netcdf files
394      !!
395      !!----------------------------------------------------------------------
396      !
397      INTEGER  :: ji, jj, jk, jm
398      INTEGER  :: numdust, numriv, numiron, numdepo
399      INTEGER  :: ierr, ierr1, ierr2, ierr3
400      REAL(wp) :: zexpide, zdenitide, zmaskt
401      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
402      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriverdic, zriverdoc, zcmask
403      !
404      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
405      TYPE(FLD_N) ::   sn_dust, sn_riverdoc, sn_riverdic, sn_ndepo, sn_ironsed        ! informations about the fields to be read
406      NAMELIST/nampissed/cn_dir, sn_dust, sn_riverdic, sn_riverdoc, sn_ndepo, sn_ironsed, &
407        &                ln_dust, ln_river, ln_ndepo, ln_ironsed,         &
408        &                sedfeinput, dustsolub, wdust, nitrfix, diazolight, concfediaz 
409      !!----------------------------------------------------------------------
410      !
411      IF( nn_timing == 1 )  CALL timing_start('p4z_sed_init')
412      !
413      !                                    ! number of seconds per year and per month
414      ryyss    = nyear_len(1) * rday
415      rmtss    = ryyss / raamo
416      r1_rday  = 1. / rday
417      r1_ryyss = 1. / ryyss
418      !                            !* set file information
419      cn_dir  = './'            ! directory in which the model is executed
420      ! ... default values (NB: frequency positive => hours, negative => months)
421      !                  !   file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
422      !                  !   name       !  (hours)  !   name      !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
423      sn_dust     = FLD_N( 'dust'       ,    -1     ,  'dust'     ,  .true.    , .true.  ,   'yearly'  , ''       , ''         )
424      sn_riverdic = FLD_N( 'river'      ,   -12     ,  'riverdic' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
425      sn_riverdoc = FLD_N( 'river'      ,   -12     ,  'riverdoc' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
426      sn_ndepo    = FLD_N( 'ndeposition',   -12     ,  'ndep'     ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
427      sn_ironsed  = FLD_N( 'ironsed'    ,   -12     ,  'bathy'    ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
428
429      REWIND( numnatp )                     ! read numnatp
430      READ  ( numnatp, nampissed )
431
432      IF(lwp) THEN
433         WRITE(numout,*) ' '
434         WRITE(numout,*) ' namelist : nampissed '
435         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
436         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust
437         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river
438         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
439         WRITE(numout,*) '    fe input from sediments                  ln_sedinput = ', ln_ironsed
440         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput
441         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub
442         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust
443         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix
444         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight
445         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
446       END IF
447
448      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN
449          ll_sbc = .TRUE.
450      ELSE
451          ll_sbc = .FALSE.
452      ENDIF
453
454      ! dust input from the atmosphere
455      ! ------------------------------
456      IF( ln_dust ) THEN
457         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
458         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
459         !
460         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
461         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
462         !
463         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
464                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
465         IF( sf_dust(1)%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
466         !
467         ! Get total input dust ; need to compute total atmospheric supply of Si in a year
468         CALL iom_open (  TRIM( sf_dust(1)%clrootname ) , numdust )
469         CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust)  ! get number of record in file
470         ALLOCATE( zdust(jpi,jpj,ntimes_dust) )
471         DO jm = 1, ntimes_dust
472            CALL iom_get( numdust, jpdom_data, TRIM( sf_dust(1)%clvar ), zdust(:,:,jm), jm )
473         END DO
474         CALL iom_close( numdust )
475         sumdepsi = 0.e0
476         DO jm = 1, ntimes_dust
477            sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
478         ENDDO
479         sumdepsi = sumdepsi * r1_ryyss * 8.8 * 0.075 / 28.1 
480         DEALLOCATE( zdust)
481      ELSE
482         dust(:,:) = 0._wp
483         sumdepsi  = 0._wp
484      END IF
485
486      ! nutrient input from rivers
487      ! --------------------------
488      IF( ln_river ) THEN
489         ALLOCATE( sf_riverdic(1), STAT=ierr1 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
490         ALLOCATE( sf_riverdoc(1), STAT=ierr2 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
491         IF( ierr1 + ierr2 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
492         !
493         CALL fld_fill( sf_riverdic, (/ sn_riverdic /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
494         CALL fld_fill( sf_riverdoc, (/ sn_riverdoc /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
495                                   ALLOCATE( sf_riverdic(1)%fnow(jpi,jpj,1)   )
496                                   ALLOCATE( sf_riverdoc(1)%fnow(jpi,jpj,1)   )
497         IF( sf_riverdic(1)%ln_tint ) ALLOCATE( sf_riverdic(1)%fdta(jpi,jpj,1,2) )
498         IF( sf_riverdoc(1)%ln_tint ) ALLOCATE( sf_riverdoc(1)%fdta(jpi,jpj,1,2) )
499         ! Get total input rivers ; need to compute total river supply in a year
500         CALL iom_open ( TRIM( sf_riverdic(1)%clrootname ), numriv )
501         CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv)
502         ALLOCATE( zriverdic(jpi,jpj,ntimes_riv) )   ;     ALLOCATE( zriverdoc(jpi,jpj,ntimes_riv) )
503         DO jm = 1, ntimes_riv
504            CALL iom_get( numriv, jpdom_data, TRIM( sf_riverdic(1)%clvar ), zriverdic(:,:,jm), jm )
505            CALL iom_get( numriv, jpdom_data, TRIM( sf_riverdoc(1)%clvar ), zriverdoc(:,:,jm), jm )
506         END DO
507         CALL iom_close( numriv )
508         ! N/P and Si releases due to coastal rivers
509         ! -----------------------------------------
510         rivpo4input = 0._wp 
511         rivalkinput = 0._wp 
512         DO jm = 1, ntimes_riv
513            rivpo4input = rivpo4input + glob_sum( ( zriverdic(:,:,jm) + zriverdoc(:,:,jm) ) * tmask(:,:,1) ) 
514            rivalkinput = rivalkinput + glob_sum(   zriverdic(:,:,jm)                       * tmask(:,:,1) ) 
515         END DO
516         rivpo4input = rivpo4input * 1E9 / 31.6_wp
517         rivalkinput = rivalkinput * 1E9 / 12._wp 
518         DEALLOCATE( zriverdic)   ;    DEALLOCATE( zriverdoc) 
519      ELSE
520         rivinp(:,:) = 0._wp
521         cotdep(:,:) = 0._wp
522         rivpo4input = 0._wp
523         rivalkinput = 0._wp
524      END IF 
525
526      ! nutrient input from dust
527      ! ------------------------
528      IF( ln_ndepo ) THEN
529         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc'
530         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
531         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
532         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
533         !
534         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
535                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   )
536         IF( sf_ndepo(1)%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) )
537         ! Get total input dust ; need to compute total atmospheric supply of N in a year
538         CALL iom_open ( TRIM( sf_ndepo(1)%clrootname ), numdepo )
539         CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep)
540         ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) )
541         DO jm = 1, ntimes_ndep
542            CALL iom_get( numdepo, jpdom_data, TRIM( sf_ndepo(1)%clvar ), zndepo(:,:,jm), jm )
543         END DO
544
545         CALL iom_close( numdepo )
546         nitdepinput = 0._wp
547         DO jm = 1, ntimes_ndep
548           nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
549         ENDDO
550         nitdepinput = nitdepinput * 7.6 / 14E6 
551         DEALLOCATE( zndepo)
552      ELSE
553         nitdep(:,:) = 0._wp
554         nitdepinput = 0._wp
555      ENDIF
556      ! coastal and island masks
557      ! ------------------------
558      IF( ln_ironsed ) THEN     
559         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron'
560         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
561         ALLOCATE( sf_ironsed(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
562         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
563         !
564         CALL fld_fill( sf_ironsed, (/ sn_ironsed /), cn_dir, 'p4z_sed_init', 'Iron from coastal supply ', 'nampissed' )
565
566         CALL iom_open ( TRIM( sf_ironsed(1)%clrootname ), numiron )
567         ALLOCATE( zcmask(jpi,jpj,jpk) )
568         CALL iom_get  ( numiron, jpdom_data, TRIM( sf_ironsed(1)%clvar ), zcmask(:,:,:), 1 )
569         CALL iom_close( numiron )
570         !
571         DO jk = 1, 5
572            DO jj = 2, jpjm1
573               DO ji = fs_2, fs_jpim1
574                  IF( tmask(ji,jj,jk) /= 0. ) THEN
575                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
576                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
577                     IF( zmaskt == 0. )   zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) ) 
578                  END IF
579               END DO
580            END DO
581         END DO
582         CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
583         DO jk = 1, jpk
584            DO jj = 1, jpj
585               DO ji = 1, jpi
586                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
587                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
588                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
589               END DO
590            END DO
591         END DO
592         ! Coastal supply of iron
593         ! -------------------------
594         ironsed(:,:,jpk) = 0._wp
595         DO jk = 1, jpkm1
596            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
597         END DO
598         DEALLOCATE( zcmask)
599      ELSE
600         ironsed(:,:,:) = 0._wp
601      ENDIF
602      !
603      IF( ll_sbc ) CALL p4z_sbc( nit000 ) 
604      !
605      IF(lwp) THEN
606         WRITE(numout,*)
607         WRITE(numout,*) '    Total input of elements from river supply'
608         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
609         WRITE(numout,*) '    N Supply   : ', rivpo4input/7.6*1E3/1E12*14.,' TgN/yr'
610         WRITE(numout,*) '    Si Supply  : ', rivalkinput/6.*1E3/1E12*32.,' TgSi/yr'
611         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr'
612         WRITE(numout,*) '    DIC Supply : ', rivpo4input*2.631*1E3*12./1E12,'TgC/yr'
613         WRITE(numout,*) 
614         WRITE(numout,*) '    Total input of elements from atmospheric supply'
615         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
616         WRITE(numout,*) '    N Supply   : ', nitdepinput/7.6*1E3/1E12*14.,' TgN/yr'
617         WRITE(numout,*) 
618      ENDIF
619      !
620      IF( nn_timing == 1 )  CALL timing_stop('p4z_sed_init')
621      !
622   END SUBROUTINE p4z_sed_init
623
624   INTEGER FUNCTION p4z_sed_alloc()
625      !!----------------------------------------------------------------------
626      !!                     ***  ROUTINE p4z_sed_alloc  ***
627      !!----------------------------------------------------------------------
628
629      ALLOCATE( dust  (jpi,jpj), rivinp(jpi,jpj)     , cotdep(jpi,jpj),      &
630        &       nitdep(jpi,jpj), ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc ) 
631
632      IF( p4z_sed_alloc /= 0 ) CALL ctl_warn('p4z_sed_alloc : failed to allocate arrays.')
633
634   END FUNCTION p4z_sed_alloc
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