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

source: branches/2011/dev_LOCEAN_2011/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 2977

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

Add in branch 2011/dev_LOCEAN_2011 changes from 2011/dev_r2787_PISCES_improvment, 2011/dev_r2787_LOCEAN_offline_fldread and 2011/dev_r2787_LOCEAN3_TRA_TRP branches, see ticket #877
  • Property svn:keywords set to Id
File size: 29.7 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   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      USE wrk_nemo, ONLY: wrk_in_USE, wrk_not_released
96      USE wrk_nemo, ONLY: zsidep   => wrk_2d_11
97      USE wrk_nemo, ONLY: zwork1   => wrk_2d_12, zwork2 => wrk_2d_13, zwork3 => wrk_2d_14
98      USE wrk_nemo, ONLY: znitrpot => wrk_3d_2, zirondep => wrk_3d_3
99      !
100      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
101      INTEGER  ::   ji, jj, jk, ikt
102#if ! defined key_sed
103      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
104      REAL(wp) ::   zrivalk, zrivsil, zrivpo4
105#endif
106      REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact, zfactcal
107      REAL(wp) ::   zsiloss, zcaloss, zwsbio3, zwsbio4, zwscal, zdep
108      CHARACTER (len=25) :: charout
109      !!---------------------------------------------------------------------
110
111      IF( ( wrk_in_USE(2, 11,12,13,14) ) .OR. ( wrk_in_USE(3, 2,3) ) ) THEN
112         CALL ctl_stop('p4z_sed: requested workspace arrays unavailable')  ;  RETURN
113      END IF
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,jpdsi) * 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,jpdsi) * 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,jpdsi) * zwsbio4
197# else
198            zsiloss = trn(ji,jj,ikt,jpdsi) * zwscal
199# endif
200            zcaloss = trn(ji,jj,ikt,jpcal) * zwscal
201            !
202            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - 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      IF( ( wrk_not_released(2, 11,12,13,14) ) .OR. ( wrk_not_released(3, 2,3) ) )   &
314        &         CALL ctl_stop('p4z_sed: failed to release workspace arrays')
315
316   END SUBROUTINE p4z_sed
317
318   SUBROUTINE p4z_sbc( kt )
319      !!----------------------------------------------------------------------
320      !!                  ***  routine p4z_sbc  ***
321      !!
322      !! ** purpose :   read and interpolate the external sources of
323      !!                nutrients
324      !!
325      !! ** method  :   read the files and interpolate the appropriate variables
326      !!
327      !! ** input   :   external netcdf files
328      !!
329      !!----------------------------------------------------------------------
330      !! * arguments
331      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
332
333      !! * local declarations
334      INTEGER  :: ji,jj 
335      REAL(wp) :: zcoef
336
337      !!---------------------------------------------------------------------
338
339      ! Compute dust at nit000 or only if there is more than 1 time record in dust file
340      IF( ln_dust ) THEN
341         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
342            CALL fld_read( kt, 1, sf_dust )
343            dust(:,:) = sf_dust(1)%fnow(:,:,1)
344         ENDIF
345      ENDIF
346
347      ! N/P and Si releases due to coastal rivers
348      ! Compute river at nit000 or only if there is more than 1 time record in river file
349      ! -----------------------------------------
350      IF( ln_river ) THEN
351         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN
352            CALL fld_read( kt, 1, sf_riverdic )
353            CALL fld_read( kt, 1, sf_riverdoc )
354            DO jj = 1, jpj
355               DO ji = 1, jpi
356                  zcoef = ryyss * cvol(ji,jj,1) 
357                  cotdep(ji,jj) =   sf_riverdic(1)%fnow(ji,jj,1)                                  * 1E9 / ( 12. * zcoef + rtrn )
358                  rivinp(ji,jj) = ( sf_riverdic(1)%fnow(ji,jj,1) + sf_riverdoc(1)%fnow(ji,jj,1) ) * 1E9 / ( 31.6* zcoef + rtrn )
359               END DO
360            END DO
361         ENDIF
362      ENDIF
363
364      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file
365      IF( ln_ndepo ) THEN
366         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
367            CALL fld_read( kt, 1, sf_ndepo )
368            DO jj = 1, jpj
369               DO ji = 1, jpi
370                  nitdep(ji,jj) = 7.6 * sf_ndepo(1)%fnow(ji,jj,1) / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn )
371               END DO
372            END DO
373         ENDIF
374      ENDIF
375      !
376   END SUBROUTINE p4z_sbc
377
378   SUBROUTINE p4z_sed_init
379
380      !!----------------------------------------------------------------------
381      !!                  ***  routine p4z_sed_init  ***
382      !!
383      !! ** purpose :   initialization of the external sources of nutrients
384      !!
385      !! ** method  :   read the files and compute the budget
386      !!                called at the first timestep (nit000)
387      !!
388      !! ** input   :   external netcdf files
389      !!
390      !!----------------------------------------------------------------------
391      !
392      INTEGER  :: ji, jj, jk, jm
393      INTEGER  :: numdust, numriv, numiron, numdepo
394      INTEGER  :: ierr, ierr1, ierr2, ierr3
395      REAL(wp) :: zexpide, zdenitide, zmaskt
396      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
397      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriverdic, zriverdoc, zcmask
398      !
399      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
400      TYPE(FLD_N) ::   sn_dust, sn_riverdoc, sn_riverdic, sn_ndepo, sn_ironsed        ! informations about the fields to be read
401      NAMELIST/nampissed/cn_dir, sn_dust, sn_riverdic, sn_riverdoc, sn_ndepo, sn_ironsed, &
402        &                ln_dust, ln_river, ln_ndepo, ln_ironsed,         &
403        &                sedfeinput, dustsolub, wdust, nitrfix, diazolight, concfediaz 
404      !!----------------------------------------------------------------------
405      !                                    ! number of seconds per year and per month
406      ryyss    = nyear_len(1) * rday
407      rmtss    = ryyss / raamo
408      r1_rday  = 1. / rday
409      r1_ryyss = 1. / ryyss
410      !                            !* set file information
411      cn_dir  = './'            ! directory in which the model is executed
412      ! ... default values (NB: frequency positive => hours, negative => months)
413      !                  !   file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
414      !                  !   name       !  (hours)  !   name      !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
415      sn_dust     = FLD_N( 'dust'       ,    -1     ,  'dust'     ,  .true.    , .true.  ,   'yearly'  , ''       , ''         )
416      sn_riverdic = FLD_N( 'river'      ,   -12     ,  'riverdic' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
417      sn_riverdoc = FLD_N( 'river'      ,   -12     ,  'riverdoc' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
418      sn_ndepo    = FLD_N( 'ndeposition',   -12     ,  'ndep'     ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
419      sn_ironsed  = FLD_N( 'ironsed'    ,   -12     ,  'bathy'    ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
420
421      REWIND( numnatp )                     ! read numnatp
422      READ  ( numnatp, nampissed )
423
424      IF(lwp) THEN
425         WRITE(numout,*) ' '
426         WRITE(numout,*) ' namelist : nampissed '
427         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
428         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust
429         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river
430         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
431         WRITE(numout,*) '    fe input from sediments                  ln_sedinput = ', ln_ironsed
432         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput
433         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub
434         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust
435         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix
436         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight
437         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
438       END IF
439
440      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN
441          ll_sbc = .TRUE.
442      ELSE
443          ll_sbc = .FALSE.
444      ENDIF
445
446      ! dust input from the atmosphere
447      ! ------------------------------
448      IF( ln_dust ) THEN
449         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
450         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
451         !
452         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
453         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
454         !
455         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
456                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
457         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
458         !
459         ! Get total input dust ; need to compute total atmospheric supply of Si in a year
460         CALL iom_open (  TRIM( sn_dust%clname ) , numdust )
461         CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust)  ! get number of record in file
462         ALLOCATE( zdust(jpi,jpj,ntimes_dust) )
463         DO jm = 1, ntimes_dust
464            CALL iom_get( numdust, jpdom_data, TRIM( sn_dust%clvar ), zdust(:,:,jm), jm )
465         END DO
466         CALL iom_close( numdust )
467         sumdepsi = 0.e0
468         DO jm = 1, ntimes_dust
469            sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
470         ENDDO
471         sumdepsi = sumdepsi * r1_ryyss * 8.8 * 0.075 / 28.1 
472         DEALLOCATE( zdust)
473      ELSE
474         dust(:,:) = 0._wp
475         sumdepsi  = 0._wp
476      END IF
477
478      ! nutrient input from rivers
479      ! --------------------------
480      IF( ln_river ) THEN
481         ALLOCATE( sf_riverdic(1), STAT=ierr1 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
482         ALLOCATE( sf_riverdoc(1), STAT=ierr2 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
483         IF( ierr1 + ierr2 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
484         !
485         CALL fld_fill( sf_riverdic, (/ sn_riverdic /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
486         CALL fld_fill( sf_riverdoc, (/ sn_riverdoc /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
487                                   ALLOCATE( sf_riverdic(1)%fnow(jpi,jpj,1)   )
488                                   ALLOCATE( sf_riverdoc(1)%fnow(jpi,jpj,1)   )
489         IF( sn_riverdic%ln_tint ) ALLOCATE( sf_riverdic(1)%fdta(jpi,jpj,1,2) )
490         IF( sn_riverdoc%ln_tint ) ALLOCATE( sf_riverdoc(1)%fdta(jpi,jpj,1,2) )
491         ! Get total input rivers ; need to compute total river supply in a year
492         CALL iom_open ( TRIM( sn_riverdic%clname ), numriv )
493         CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv)
494         ALLOCATE( zriverdic(jpi,jpj,ntimes_riv) )   ;     ALLOCATE( zriverdoc(jpi,jpj,ntimes_riv) )
495         DO jm = 1, ntimes_riv
496            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdic%clvar ), zriverdic(:,:,jm), jm )
497            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdoc%clvar ), zriverdoc(:,:,jm), jm )
498         END DO
499         CALL iom_close( numriv )
500         ! N/P and Si releases due to coastal rivers
501         ! -----------------------------------------
502         rivpo4input = 0._wp 
503         rivalkinput = 0._wp 
504         DO jm = 1, ntimes_riv
505            rivpo4input = rivpo4input + glob_sum( ( zriverdic(:,:,jm) + zriverdoc(:,:,jm) ) * tmask(:,:,1) ) 
506            rivalkinput = rivalkinput + glob_sum(   zriverdic(:,:,jm)                       * tmask(:,:,1) ) 
507         END DO
508         rivpo4input = rivpo4input * 1E9 / 31.6_wp
509         rivalkinput = rivalkinput * 1E9 / 12._wp 
510         DEALLOCATE( zriverdic)   ;    DEALLOCATE( zriverdoc) 
511      ELSE
512         rivinp(:,:) = 0._wp
513         cotdep(:,:) = 0._wp
514         rivpo4input = 0._wp
515         rivalkinput = 0._wp
516      END IF 
517
518      ! nutrient input from dust
519      ! ------------------------
520      IF( ln_ndepo ) THEN
521         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc'
522         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
523         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
524         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
525         !
526         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
527                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   )
528         IF( sn_ndepo%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) )
529         !
530         ! Get total input dust ; need to compute total atmospheric supply of N in a year
531         CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo )
532         CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep)
533         ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) )
534         DO jm = 1, ntimes_ndep
535            CALL iom_get( numdepo, jpdom_data, TRIM( sn_ndepo%clvar ), zndepo(:,:,jm), jm )
536         END DO
537         CALL iom_close( numdepo )
538         nitdepinput = 0._wp
539         DO jm = 1, ntimes_ndep
540           nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
541         ENDDO
542         nitdepinput = nitdepinput * 7.6 / 14E6 
543         DEALLOCATE( zndepo)
544      ELSE
545         nitdep(:,:) = 0._wp
546         nitdepinput = 0._wp
547      ENDIF
548
549      ! coastal and island masks
550      ! ------------------------
551      IF( ln_ironsed ) THEN     
552         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron'
553         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
554         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
555         ALLOCATE( zcmask(jpi,jpj,jpk) )
556         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
557         CALL iom_close( numiron )
558         !
559         DO jk = 1, 5
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1
562                  IF( tmask(ji,jj,jk) /= 0. ) THEN
563                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
564                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
565                     IF( zmaskt == 0. )   zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) ) 
566                  END IF
567               END DO
568            END DO
569         END DO
570         CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
571         DO jk = 1, jpk
572            DO jj = 1, jpj
573               DO ji = 1, jpi
574                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
575                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
576                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
577               END DO
578            END DO
579         END DO
580         ! Coastal supply of iron
581         ! -------------------------
582         ironsed(:,:,jpk) = 0._wp
583         DO jk = 1, jpkm1
584            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
585         END DO
586         DEALLOCATE( zcmask)
587      ELSE
588         ironsed(:,:,:) = 0._wp
589      ENDIF
590      !
591      IF(lwp) THEN
592         WRITE(numout,*)
593         WRITE(numout,*) '    Total input of elements from river supply'
594         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
595         WRITE(numout,*) '    N Supply   : ', rivpo4input/7.6*1E3/1E12*14.,' TgN/yr'
596         WRITE(numout,*) '    Si Supply  : ', rivalkinput/6.*1E3/1E12*32.,' TgSi/yr'
597         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr'
598         WRITE(numout,*) '    DIC Supply : ', rivpo4input*2.631*1E3*12./1E12,'TgC/yr'
599         WRITE(numout,*) 
600         WRITE(numout,*) '    Total input of elements from atmospheric supply'
601         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
602         WRITE(numout,*) '    N Supply   : ', nitdepinput/7.6*1E3/1E12*14.,' TgN/yr'
603         WRITE(numout,*) 
604      ENDIF
605       !
606   END SUBROUTINE p4z_sed_init
607
608   INTEGER FUNCTION p4z_sed_alloc()
609      !!----------------------------------------------------------------------
610      !!                     ***  ROUTINE p4z_sed_alloc  ***
611      !!----------------------------------------------------------------------
612
613      ALLOCATE( dust  (jpi,jpj), rivinp(jpi,jpj)     , cotdep(jpi,jpj),      &
614        &       nitdep(jpi,jpj), ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc ) 
615
616      IF( p4z_sed_alloc /= 0 ) CALL ctl_warn('p4z_sed_alloc : failed to allocate arrays.')
617
618   END FUNCTION p4z_sed_alloc
619#else
620   !!======================================================================
621   !!  Dummy module :                                   No PISCES bio-model
622   !!======================================================================
623CONTAINS
624   SUBROUTINE p4z_sed                         ! Empty routine
625   END SUBROUTINE p4z_sed
626#endif 
627
628   !!======================================================================
629END MODULE  p4zsed
Note: See TracBrowser for help on using the repository browser.