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/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES – NEMO

source: branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 2957

Last change on this file since 2957 was 2957, checked in by cetlod, 13 years ago

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