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 @ 2852

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

dev_r2787_PISCES_improvment:minor bug correction in nitrogen fixation diagnostic
  • Property svn:keywords set to Id
File size: 29.6 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 ) * cvol(:,:,:) * xnegtr(:,:,:) )
253
254      ! Potential nitrogen fixation dependant on temperature and iron
255      ! -------------------------------------------------------------
256
257!CDIR NOVERRCHK
258      DO jk = 1, jpk
259!CDIR NOVERRCHK
260         DO jj = 1, jpj
261!CDIR NOVERRCHK
262            DO ji = 1, jpi
263               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
264               IF( zlim <= 0.2 )   zlim = 0.01
265#if defined key_degrad
266               zfact = zlim * rfact2 * facvol(ji,jj,jk)
267#else
268               zfact = zlim * rfact2 
269#endif
270               znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
271                 &                 *  zfact * trn(ji,jj,jk,jpfer) / ( concfediaz + trn(ji,jj,jk,jpfer) ) &
272                 &                 * ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) )
273            END DO
274         END DO
275      END DO
276
277      znitrpottot = glob_sum( znitrpot(:,:,:) * cvol(:,:,:) )
278
279      ! Nitrogen change due to nitrogen fixation
280      ! ----------------------------------------
281      DO jk = 1, jpk
282         DO jj = 1, jpj
283            DO ji = 1, jpi
284               zfact = znitrpot(ji,jj,jk) * nitrfix
285               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
286               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3 * zfact
287               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
288               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30. / 46. * zfact
289           !    trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + zfact
290           END DO
291         END DO
292      END DO
293
294#if defined key_diatrc
295      zfact = 1.e+3 * rfact2r
296#if defined key_iomput
297      zwork1(:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1) 
298      zwork2(:,:)  =    znitrpot(:,:,1) * nitrfix                   * zfact * fse3t(:,:,1) * tmask(:,:,1)
299      IF( jnt == nrdttrc ) THEN
300          CALL iom_put( "Irondep", zwork1  )  ! surface downward net flux of iron
301          CALL iom_put( "Nfix"   , zwork2 )  ! nitrogen fixation at surface
302      ENDIF
303#else
304      trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)           * zfact * fse3t(:,:,1) * tmask(:,:,1)
305      trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1)
306#endif
307      !
308      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging)
309         WRITE(charout, fmt="('sed ')")
310         CALL prt_ctl_trc_info(charout)
311         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
312      ENDIF
313
314      IF( ( wrk_not_released(2, 11,12,13,14) ) .OR. ( wrk_not_released(3, 2,3) ) )   &
315        &         CALL ctl_stop('p4z_sed: failed to release workspace arrays')
316
317   END SUBROUTINE p4z_sed
318
319   SUBROUTINE p4z_sbc( kt )
320      !!----------------------------------------------------------------------
321      !!                  ***  routine p4z_sbc  ***
322      !!
323      !! ** purpose :   read and interpolate the external sources of
324      !!                nutrients
325      !!
326      !! ** method  :   read the files and interpolate the appropriate variables
327      !!
328      !! ** input   :   external netcdf files
329      !!
330      !!----------------------------------------------------------------------
331      !! * arguments
332      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
333
334      !! * local declarations
335      INTEGER  :: ji,jj 
336      REAL(wp) :: zcoef
337
338      !!---------------------------------------------------------------------
339
340      ! Compute dust at nit000 or only if there is more than 1 time record in dust file
341      IF( ln_dust ) THEN
342         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
343            CALL fld_read( kt, 1, sf_dust )
344            dust(:,:) = sf_dust(1)%fnow(:,:,1)
345         ENDIF
346      ENDIF
347
348      ! N/P and Si releases due to coastal rivers
349      ! Compute river at nit000 or only if there is more than 1 time record in river file
350      ! -----------------------------------------
351      IF( ln_river ) THEN
352         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN
353            CALL fld_read( kt, 1, sf_riverdic )
354            CALL fld_read( kt, 1, sf_riverdoc )
355            DO jj = 1, jpj
356               DO ji = 1, jpi
357                  zcoef = ryyss * cvol(ji,jj,1) 
358                  cotdep(ji,jj) =   sf_riverdic(1)%fnow(ji,jj,1)                                  * 1E9 / ( 12. * zcoef + rtrn )
359                  rivinp(ji,jj) = ( sf_riverdic(1)%fnow(ji,jj,1) + sf_riverdoc(1)%fnow(ji,jj,1) ) * 1E9 / ( 31.6* zcoef + rtrn )
360               END DO
361            END DO
362         ENDIF
363      ENDIF
364
365      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file
366      IF( ln_ndepo ) THEN
367         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
368            CALL fld_read( kt, 1, sf_ndepo )
369            DO jj = 1, jpj
370               DO ji = 1, jpi
371                  nitdep(ji,jj) = 7.6 * sf_ndepo(1)%fnow(ji,jj,1) / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn )
372               END DO
373            END DO
374         ENDIF
375      ENDIF
376      !
377   END SUBROUTINE p4z_sbc
378
379
380   SUBROUTINE p4z_sed_init
381
382      !!----------------------------------------------------------------------
383      !!                  ***  routine p4z_sed_init  ***
384      !!
385      !! ** purpose :   initialization of the external sources of nutrients
386      !!
387      !! ** method  :   read the files and compute the budget
388      !!                called at the first timestep (nit000)
389      !!
390      !! ** input   :   external netcdf files
391      !!
392      !!----------------------------------------------------------------------
393      !
394      INTEGER  :: ji, jj, jk, jm
395      INTEGER  :: numdust, numriv, numiron, numdepo
396      INTEGER  :: ierr, ierr1, ierr2, ierr3
397      REAL(wp) :: zexpide, zdenitide, zmaskt
398      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
399      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriverdic, zriverdoc, zcmask
400      !
401      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
402      TYPE(FLD_N) ::   sn_dust, sn_riverdoc, sn_riverdic, sn_ndepo, sn_ironsed        ! informations about the fields to be read
403      NAMELIST/nampissed/cn_dir, sn_dust, sn_riverdic, sn_riverdoc, sn_ndepo, sn_ironsed, &
404        &                ln_dust, ln_river, ln_ndepo, ln_ironsed,         &
405        &                sedfeinput, dustsolub, wdust, nitrfix, diazolight, concfediaz 
406      !!----------------------------------------------------------------------
407      !                                    ! number of seconds per year and per month
408      ryyss    = nyear_len(1) * rday
409      rmtss    = ryyss / raamo
410      r1_rday  = 1. / rday
411      r1_ryyss = 1. / ryyss
412      !                            !* set file information
413      cn_dir  = './'            ! directory in which the model is executed
414      ! ... default values (NB: frequency positive => hours, negative => months)
415      !                  !   file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
416      !                  !   name       !  (hours)  !   name      !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
417      sn_dust     = FLD_N( 'dust'       ,    -1     ,  'dust'     ,  .true.    , .true.  ,   'yearly'  , ''       , ''         )
418      sn_riverdic = FLD_N( 'river'      ,   -12     ,  'riverdic' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
419      sn_riverdoc = FLD_N( 'river'      ,   -12     ,  'riverdoc' ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
420      sn_ndepo    = FLD_N( 'ndeposition',   -12     ,  'ndep'     ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
421      sn_ironsed  = FLD_N( 'ironsed'    ,   -12     ,  'bathy'    ,  .false.   , .true.  ,   'yearly'  , ''       , ''         )
422
423
424      REWIND( numnatp )                     ! read numnat
425      READ  ( numnatp, nampissed )
426
427      IF(lwp) THEN
428         WRITE(numout,*) ' '
429         WRITE(numout,*) ' namelist : nampissed '
430         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
431         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust
432         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river
433         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
434         WRITE(numout,*) '    fe input from sediments                  ln_sedinput = ', ln_ironsed
435         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput
436         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub
437         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust
438         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix
439         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight
440         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
441       END IF
442
443      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN
444          ll_sbc = .TRUE.
445      ELSE
446          ll_sbc = .FALSE.
447      ENDIF
448
449      ! dust input from the atmosphere
450      ! ------------------------------
451      IF( ln_dust ) THEN
452         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
453         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
454         !
455         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
456         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
457         !
458         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
459                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
460         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
461         !
462         ! Get total input dust ; need to compute total atmospheric supply of Si in a year
463         CALL iom_open (  TRIM( sn_dust%clname ) , numdust )
464         CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust)  ! get number of record in file
465         ALLOCATE( zdust(jpi,jpj,ntimes_dust) )
466         DO jm = 1, ntimes_dust
467            CALL iom_get( numdust, jpdom_data, TRIM( sn_dust%clvar ), zdust(:,:,jm), jm )
468         END DO
469         CALL iom_close( numdust )
470         sumdepsi = 0.e0
471         DO jm = 1, ntimes_dust
472            sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
473         ENDDO
474         sumdepsi = sumdepsi * r1_ryyss * 8.8 * 0.075 / 28.1 
475         DEALLOCATE( zdust)
476      ELSE
477         dust(:,:) = 0._wp
478         sumdepsi  = 0._wp
479      END IF
480
481      ! nutrient input from rivers
482      ! --------------------------
483      IF( ln_river ) THEN
484         ALLOCATE( sf_riverdic(1), STAT=ierr1 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
485         ALLOCATE( sf_riverdoc(1), STAT=ierr2 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
486         IF( ierr1 + ierr2 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
487         !
488         CALL fld_fill( sf_riverdic, (/ sn_riverdic /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
489         CALL fld_fill( sf_riverdoc, (/ sn_riverdoc /), cn_dir, 'p4z_sed_init', 'Input DOC from river ', 'nampissed' )
490                                   ALLOCATE( sf_riverdic(1)%fnow(jpi,jpj,1)   )
491                                   ALLOCATE( sf_riverdoc(1)%fnow(jpi,jpj,1)   )
492         IF( sn_riverdic%ln_tint ) ALLOCATE( sf_riverdic(1)%fdta(jpi,jpj,1,2) )
493         IF( sn_riverdoc%ln_tint ) ALLOCATE( sf_riverdoc(1)%fdta(jpi,jpj,1,2) )
494         ! Get total input rivers ; need to compute total river supply in a year
495         CALL iom_open ( TRIM( sn_riverdic%clname ), numriv )
496         CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv)
497         ALLOCATE( zriverdic(jpi,jpj,ntimes_riv) )   ;     ALLOCATE( zriverdoc(jpi,jpj,ntimes_riv) )
498         DO jm = 1, ntimes_riv
499            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdic%clvar ), zriverdic(:,:,jm), jm )
500            CALL iom_get( numriv, jpdom_data, TRIM( sn_riverdoc%clvar ), zriverdoc(:,:,jm), jm )
501         END DO
502         CALL iom_close( numriv )
503         ! N/P and Si releases due to coastal rivers
504         ! -----------------------------------------
505         rivpo4input = 0._wp 
506         rivalkinput = 0._wp 
507         DO jm = 1, ntimes_riv
508            rivpo4input = rivpo4input + glob_sum( ( zriverdic(:,:,jm) + zriverdoc(:,:,jm) ) * tmask(:,:,1) ) 
509            rivalkinput = rivalkinput + glob_sum(   zriverdic(:,:,jm)                       * tmask(:,:,1) ) 
510         END DO
511         rivpo4input = rivpo4input * 1E9 / 31.6_wp
512         rivalkinput = rivalkinput * 1E9 / 12._wp 
513         DEALLOCATE( zriverdic)   ;    DEALLOCATE( zriverdoc) 
514      ELSE
515         rivinp(:,:) = 0._wp
516         cotdep(:,:) = 0._wp
517         rivpo4input = 0._wp
518         rivalkinput = 0._wp
519      END IF 
520
521      ! nutrient input from dust
522      ! ------------------------
523      IF( ln_ndepo ) THEN
524         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc'
525         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
526         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
527         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_apr structure' )
528         !
529         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Iron from sediment ', 'nampissed' )
530                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   )
531         IF( sn_ndepo%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) )
532         !
533         ! Get total input dust ; need to compute total atmospheric supply of N in a year
534         CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo )
535         CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep)
536         ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) )
537         DO jm = 1, ntimes_ndep
538            CALL iom_get( numdepo, jpdom_data, TRIM( sn_ndepo%clvar ), zndepo(:,:,jm), jm )
539         END DO
540         CALL iom_close( numdepo )
541         nitdepinput = 0._wp
542         DO jm = 1, ntimes_ndep
543           nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) ) 
544         ENDDO
545         nitdepinput = nitdepinput * 7.6 / 14E6 
546         DEALLOCATE( zndepo)
547      ELSE
548         nitdep(:,:) = 0._wp
549         nitdepinput = 0._wp
550      ENDIF
551
552      ! coastal and island masks
553      ! ------------------------
554      IF( ln_ironsed ) THEN     
555         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron'
556         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
557         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
558         ALLOCATE( zcmask(jpi,jpj,jpk) )
559         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
560         CALL iom_close( numiron )
561         !
562         DO jk = 1, 5
563            DO jj = 2, jpjm1
564               DO ji = fs_2, fs_jpim1
565                  IF( tmask(ji,jj,jk) /= 0. ) THEN
566                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
567                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
568                     IF( zmaskt == 0. )   zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) ) 
569                  END IF
570               END DO
571            END DO
572         END DO
573         CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
574         DO jk = 1, jpk
575            DO jj = 1, jpj
576               DO ji = 1, jpi
577                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
578                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
579                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
580               END DO
581            END DO
582         END DO
583         ! Coastal supply of iron
584         ! -------------------------
585         ironsed(:,:,jpk) = 0._wp
586         DO jk = 1, jpkm1
587            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
588         END DO
589         DEALLOCATE( zcmask)
590      ELSE
591         ironsed(:,:,:) = 0._wp
592      ENDIF
593      !
594      IF(lwp) THEN
595         WRITE(numout,*)
596         WRITE(numout,*) '    Total input of elements from river supply'
597         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
598         WRITE(numout,*) '    N Supply   : ', rivpo4input/7.6*1E3/1E12*14.,' TgN/yr'
599         WRITE(numout,*) '    Si Supply  : ', rivalkinput/6.*1E3/1E12*32.,' TgSi/yr'
600         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr'
601         WRITE(numout,*) '    DIC Supply : ', rivpo4input*2.631*1E3*12./1E12,'TgC/yr'
602         WRITE(numout,*) 
603         WRITE(numout,*) '    Total input of elements from atmospheric supply'
604         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
605         WRITE(numout,*) '    N Supply   : ', nitdepinput/7.6*1E3/1E12*14.,' TgN/yr'
606         WRITE(numout,*) 
607      ENDIF
608       !
609   END SUBROUTINE p4z_sed_init
610
611   INTEGER FUNCTION p4z_sed_alloc()
612      !!----------------------------------------------------------------------
613      !!                     ***  ROUTINE p4z_sed_alloc  ***
614      !!----------------------------------------------------------------------
615
616      ALLOCATE( dust  (jpi,jpj), rivinp(jpi,jpj)     , cotdep(jpi,jpj),      &
617        &       nitdep(jpi,jpj), ironsed(jpi,jpj,jpk), STAT=p4z_sed_alloc ) 
618
619      IF( p4z_sed_alloc /= 0 ) CALL ctl_warn('p4z_sed_alloc : failed to allocate arrays.')
620
621   END FUNCTION p4z_sed_alloc
622#else
623   !!======================================================================
624   !!  Dummy module :                                   No PISCES bio-model
625   !!======================================================================
626CONTAINS
627   SUBROUTINE p4z_sed                         ! Empty routine
628   END SUBROUTINE p4z_sed
629#endif 
630
631   !!======================================================================
632END MODULE  p4zsed
Note: See TracBrowser for help on using the repository browser.