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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 3147

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

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