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

source: trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 1073

Last change on this file since 1073 was 1073, checked in by cetlod, 16 years ago

update PISCES model, see ticket:190
  • Property svn:executable set to *
File size: 23.2 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   !!----------------------------------------------------------------------
9#if defined key_pisces
10   !!----------------------------------------------------------------------
11   !!   'key_pisces'                                       PISCES bio-model
12   !!----------------------------------------------------------------------
13   !!   p4z_sed        :  Compute loss of organic matter in the sediments
14   !!   p4z_sbc        :  Read and interpolate time-varying nutrients fluxes
15   !!   p4z_sed_init   :  Initialization of p4z_sed
16   !!----------------------------------------------------------------------
17   USE trc
18   USE oce_trc         !
19   USE trp_trc
20   USE sms_pisces
21   USE lib_mpp
22   USE prtctl_trc
23   USE p4zbio
24   USE p4zint
25   USE p4zopt
26   USE p4zsink
27   USE p4zrem
28   USE p4zlim
29   USE lbclnk
30   USE iom
31
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   p4z_sed   
37
38   !! * Shared module variables
39   LOGICAL, PUBLIC ::    &
40     bdustfer  = .FALSE.      ,  &  !:
41     briver    = .FALSE.      ,  &  !:
42     bndepo    = .FALSE.      ,  &  !:
43     bsedinput = .FALSE.            !:
44
45   REAL(wp), PUBLIC ::   &
46     sedfeinput = 1.E-9_wp   ,  &  !:
47     dustsolub  = 0.014_wp         !:
48
49   !! * Module variables
50   INTEGER ::          &
51      numdust,          &  ! logical unit for surface fluxes data
52      nflx1 , nflx2,   &  !  first and second record used
53      nflx11, nflx12      ! ???
54   REAL(wp), DIMENSION(jpi,jpj,2) ::    &  !:
55     dustmo                                !: 2 consecutive set of dust fields
56   REAL(wp), DIMENSION(jpi,jpj)   ::    &
57     rivinp, cotdep, nitdep, dust
58   REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   &
59     ironsed
60   REAL(wp) :: sumdepsi, rivalkinput, rivpo4input, nitdepinput
61
62   !!* Substitution
63#  include "domzgr_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
66   !! $Header:$
67   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69
70CONTAINS
71
72   SUBROUTINE p4z_sed(kt, jnt)
73      !!---------------------------------------------------------------------
74      !!                     ***  ROUTINE p4z_sed  ***
75      !!
76      !! ** Purpose :   Compute loss of organic matter in the sediments. This
77      !!              is by no way a sediment model. The loss is simply
78      !!              computed to balance the inout from rivers and dust
79      !!
80      !! ** Method  : - ???
81      !!---------------------------------------------------------------------
82      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
83      INTEGER  ::   ji, jj, jk
84      INTEGER  ::   ikt
85      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
86      REAL(wp) ::   zconctmp , zdenitot  , znitrpottot
87      REAL(wp) ::   zlim, zconctmp2, zstep, zfact
88      REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
90      CHARACTER (len=25) :: charout
91      !!---------------------------------------------------------------------
92
93
94      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_sed_init      ! Initialization (first time-step only)
95
96      IF( (jnt == 1) .and. (bdustfer) )  CALL p4z_sbc( kt )
97
98      zstep = rfact2 / rjjss      ! Time step duration for the biology
99
100      zirondep(:,:,:) = 0.e0          ! Initialisation of variables used to compute deposition
101      zsidep  (:,:)   = 0.e0
102
103      ! Iron and Si deposition at the surface
104      ! -------------------------------------
105
106      DO jj = 1, jpj
107         DO ji = 1, jpi
108            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmoss ) + 3.e-10 / raass )   &
109               &             * rfact2 / fse3t(ji,jj,1)
110            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / ( fse3t(ji,jj,1) * 28.1 * rmoss )
111         END DO
112      END DO
113
114      ! Iron solubilization of particles in the water column
115      ! ----------------------------------------------------
116
117      DO jk = 2, jpkm1
118         zirondep(:,:,jk) = dust(:,:) / ( 10. * 55.85 * rmoss ) * rfact2 * 1.e-4
119      END DO
120
121      ! Add the external input of nutrients, carbon and alkalinity
122      ! ----------------------------------------------------------
123
124      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2 
125      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2
126      trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivinp(:,:) * 3.e-5 * rfact2
127      trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep (:,:) + cotdep(:,:)   * rfact2 / 6.
128      trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivinp(:,:) * 2.631 * rfact2
129      trn(:,:,1,jptal) = trn(:,:,1,jptal) + (cotdep(:,:) - rno3*(rivinp(:,:) +  nitdep(:,:) ) ) * rfact2
130
131
132      ! Add the external input of iron which is 3D distributed
133      ! (dust, river and sediment mobilization)
134      ! ------------------------------------------------------
135
136      DO jk = 1, jpkm1
137         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer)   &
138         &       + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
139      END DO
140
141      ! Initialisation of variables used to compute Sinking Speed
142      ! ---------------------------------------------------------
143
144      zsumsedsi  = 0.e0
145      zsumsedpo4 = 0.e0
146      zsumsedcal = 0.e0
147
148      ! Loss of biogenic silicon, Caco3 organic carbon in the sediments.
149      ! First, the total loss is computed.
150      ! The factor for calcite comes from the alkalinity effect
151      ! -------------------------------------------------------------
152
153      DO jj = 1, jpj
154         DO ji = 1, jpi
155            ikt = MAX( mbathy(ji,jj)-1, 1 )
156            zfact = e1t(ji,jj) * e2t(ji,jj) / rjjss * tmask_i(ji,jj)
157
158# if defined key_kriest
159            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
160            zsumsedpo4 = zsumsedpo4 + zfact *  trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
161# else
162            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt)
163            zsumsedpo4 = zsumsedpo4 + zfact *( trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt)   &
164               &       + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )
165# endif
166
167            zsumsedcal = zsumsedcal + zfact *  trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0
168
169         END DO
170      END DO
171
172      IF( lk_mpp ) THEN
173         CALL mpp_sum( zsumsedsi  )   ! sums over the global domain
174         CALL mpp_sum( zsumsedcal )   ! sums over the global domain
175         CALL mpp_sum( zsumsedpo4 )   ! sums over the global domain
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
184      DO jj = 1, jpj
185         DO ji = 1, jpi
186            ikt = MAX( mbathy(ji,jj) - 1, 1 )
187            zconctmp = trn(ji,jj,ikt,jpdsi) * zstep / fse3t(ji,jj,ikt)   &
188# if ! defined key_kriest
189     &             * wscal (ji,jj,ikt)
190# else
191     &             * wsbio4(ji,jj,ikt)
192# endif
193
194            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zconctmp
195            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zconctmp   &
196            &      * ( 1.- ( sumdepsi + rivalkinput / raass / 6. ) / zsumsedsi )
197         END DO
198      END DO
199
200      DO jj = 1, jpj
201         DO ji = 1, jpi
202            ikt = MAX( mbathy(ji,jj) - 1, 1 )
203            zconctmp = trn(ji,jj,ikt,jpcal) * wscal(ji,jj,ikt) * zstep / fse3t(ji,jj,ikt)
204            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zconctmp
205            trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + zconctmp   &
206               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal ) * 2.e0
207            trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zconctmp   &
208               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal )
209         END DO
210      END DO
211
212      DO jj = 1, jpj
213         DO ji = 1, jpi
214            ikt = MAX( mbathy(ji,jj) - 1, 1 )
215            zfact = zstep / fse3t(ji,jj,ikt)
216
217# if ! defined key_kriest
218            zconctmp  = trn(ji,jj,ikt,jpgoc)
219            zconctmp2 = trn(ji,jj,ikt,jppoc)
220            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - zconctmp  * wsbio4(ji,jj,ikt) * zfact
221            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
222            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
223            &      + ( zconctmp  * wsbio4(ji,jj,ikt) + zconctmp2 * wsbio3(ji,jj,ikt) ) * zfact   &
224            &      * ( 1.- rivpo4input / (raass * zsumsedpo4 ) )
225            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * wsbio4(ji,jj,ikt) * zfact
226            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
227
228# else
229
230            zconctmp  = trn(ji,jj,ikt,jpnum)
231            zconctmp2 = trn(ji,jj,ikt,jppoc)
232            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum)   &
233            &      - zconctmp  * wsbio4(ji,jj,ikt) * zfact
234            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc)   &
235            &      - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
236            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
237            &      + ( zconctmp2 * wsbio3(ji,jj,ikt) )   &
238            &      * zfact * ( 1.- rivpo4input / ( raass * zsumsedpo4 ) )
239            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe)   &
240            &      - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
241
242# endif
243         END DO
244      END DO
245
246      ! Nitrogen fixation (simple parameterization). The total gain
247      ! from nitrogen fixation is scaled to balance the loss by
248      ! denitrification
249      ! -------------------------------------------------------------
250
251      zdenitot = 0.e0
252      DO jk = 1, jpkm1
253         DO jj = 1,jpj
254            DO ji = 1,jpi
255               zdenitot = zdenitot + denitr(ji,jj,jk) * rdenit * e1t(ji,jj) * e2t(ji,jj)   &
256               &        *fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) * xnegtr(ji,jj,jk)
257            END DO
258         END DO
259      END DO
260
261      IF( lk_mpp )   CALL mpp_sum( zdenitot )      ! sum over the global domain
262
263      ! Potential nitrogen fication dependant on temperature and iron
264      ! -------------------------------------------------------------
265
266!CDIR NOVERRCHK
267      DO jk = 1, jpk
268!CDIR NOVERRCHK
269         DO jj = 1, jpj
270!CDIR NOVERRCHK
271            DO ji = 1, jpi
272               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
273               IF( zlim <= 0.2 )   zlim = 0.01
274               znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) / rjjss )   &
275# if defined key_off_degrad
276               &                  * facvol(ji,jj,jk)   &
277# endif
278               &                  * zlim * rfact2 * trn(ji,jj,jk,jpfer)   &
279               &                  / ( conc3 + trn(ji,jj,jk,jpfer) ) * ( 1.- EXP( -etot(ji,jj,jk) / 50.) )
280            END DO
281         END DO
282      END DO
283
284      znitrpottot = 0.e0
285      DO jk = 1, jpkm1
286         DO jj = 1, jpj
287            DO ji = 1, jpi
288               znitrpottot = znitrpottot + znitrpot(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)   &
289                  &                                           * tmask(ji,jj,jk) * tmask_i(ji,jj) 
290            END DO
291         END DO
292      END DO
293
294      IF( lk_mpp )   CALL mpp_sum( znitrpottot )  ! sum over the global domain
295
296      ! Nitrogen change due to nitrogen fixation
297      ! ----------------------------------------
298
299      DO jk = 1, jpk
300         DO jj = 1, jpj
301            DO ji = 1, jpi
302# if ! defined key_c1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 )
303!!             zfact = znitrpot(ji,jj,jk) * zdenitot / znitrpottot
304               zfact = znitrpot(ji,jj,jk) * 1.e-7
305# else
306               zfact = znitrpot(ji,jj,jk) * 1.e-7
307# endif
308               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
309               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
310               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30./ 46.* zfact
311            END DO
312         END DO
313      END DO
314
315# if defined key_trc_diaadd
316      DO jj = 1,jpj
317         DO ji = 1,jpi
318            trc2d(ji,jj,13) = znitrpot(ji,jj,1) * 1.e-7 * fse3t(ji,jj,1) * 1.e+3 / rfact2
319            trc2d(ji,jj,12) = zirondep(ji,jj,1) * 1.e+3 * rfact2r * fse3t(ji,jj,1)
320         END DO
321      END DO
322# endif
323      !
324       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
325         WRITE(charout, FMT="('sed ')")
326         CALL prt_ctl_trc_info(charout)
327         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
328       ENDIF
329
330   END SUBROUTINE p4z_sed
331
332   SUBROUTINE p4z_sbc(kt)
333
334      !!----------------------------------------------------------------------
335      !!                  ***  ROUTINE p4z_sbc  ***
336      !!
337      !! ** Purpose :   Read and interpolate the external sources of
338      !!                nutrients
339      !!
340      !! ** Method  :   Read the files and interpolate the appropriate variables
341      !!
342      !! ** input   :   external netcdf files
343      !!
344      !!----------------------------------------------------------------------
345      !! * arguments
346      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
347
348      !! * Local declarations
349      INTEGER ::   &
350         imois, imois2,       &  ! temporary integers
351         i15  , iman             !    "          "
352      REAL(wp) ::   &
353         zxy                     !    "         "
354
355
356      !!---------------------------------------------------------------------
357
358      ! Initialization
359      ! --------------
360
361      i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
362      iman  = INT( raamo )
363      imois = nmonth + i15 - 1
364      IF( imois == 0 ) imois = iman
365      imois2 = nmonth
366
367      ! 1. first call kt=nit000
368      ! -----------------------
369
370      IF( kt == nit000 ) THEN
371         ! initializations
372         nflx1  = 0
373         nflx11 = 0
374         ! open the file
375         IF(lwp) THEN
376            WRITE(numout,*) ' '
377            WRITE(numout,*) ' **** Routine p4z_sbc'
378         ENDIF
379         CALL iom_open ( 'dust.orca.nc', numdust )
380      ENDIF
381
382
383     ! Read monthly file
384      ! ----------------
385
386      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
387
388         ! Calendar computation
389
390         ! nflx1 number of the first file record used in the simulation
391         ! nflx2 number of the last  file record
392
393         nflx1 = imois
394         nflx2 = nflx1+1
395         nflx1 = MOD( nflx1, iman )
396         nflx2 = MOD( nflx2, iman )
397         IF( nflx1 == 0 )   nflx1 = iman
398         IF( nflx2 == 0 )   nflx2 = iman
399         IF(lwp) WRITE(numout,*) 'first record file used nflx1 ',nflx1
400         IF(lwp) WRITE(numout,*) 'last  record file used nflx2 ',nflx2
401
402         ! Read monthly fluxes data
403
404         ! humidity
405         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
406         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
407
408         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
409            WRITE(numout,*)
410            WRITE(numout,*) ' read clio flx ok'
411            WRITE(numout,*)
412               WRITE(numout,*)
413               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
414               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
415         ENDIF
416
417      ENDIF
418
419     ! 3. at every time step interpolation of fluxes
420      ! ---------------------------------------------
421
422      zxy = FLOAT( nday ) / FLOAT( nobis(nflx1) ) + 0.5 - i15
423
424      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
425
426      IF( kt == nitend ) CALL iom_close (numdust)
427
428   END SUBROUTINE p4z_sbc
429
430
431   SUBROUTINE p4z_sed_init
432
433      !!----------------------------------------------------------------------
434      !!                  ***  ROUTINE p4z_sed_init  ***
435      !!
436      !! ** Purpose :   Initialization of the external sources of nutrients
437      !!
438      !! ** Method  :   Read the files and compute the budget
439      !!      called at the first timestep (nittrc000)
440      !!
441      !! ** input   :   external netcdf files
442      !!
443      !!----------------------------------------------------------------------
444
445      INTEGER ::   ji, jj, jk, jm
446      INTEGER , PARAMETER ::   jpmois = 12, jpan = 1
447      INTEGER :: numriv, numbath, numdep
448
449
450      REAL(wp) ::   zcoef
451      REAL(wp) ::   expide, denitide,zmaskt
452      REAL(wp) , DIMENSION (jpi,jpj)     ::   riverdoc, river, ndepo
453      REAL(wp) , DIMENSION (jpi,jpj,jpk) ::   cmask
454      REAL(wp), DIMENSION(jpi,jpj,12)    ::   zdustmo
455
456      NAMELIST/natsms/ bdustfer, briver, bndepo, bsedinput, sedfeinput, dustsolub
457
458
459      REWIND( numnat )                     ! read numnat
460      READ  ( numnat, natsms )
461
462      IF(lwp) THEN
463         WRITE(numout,*) ' '
464         WRITE(numout,*) ' Namelist : natsms '
465         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
466         WRITE(numout,*) '    Dust input from the atmosphere           bdustfer  = ', bdustfer
467         WRITE(numout,*) '    River input of nutrients                 briver    = ', briver
468         WRITE(numout,*) '    Atmospheric deposition of N              bndepo    = ', bndepo
469         WRITE(numout,*) '    Fe input from sediments                  bsedinput = ', bsedinput
470         WRITE(numout,*) '    Coastal release of Iron                  sedfeinput =', sedfeinput
471         WRITE(numout,*) '    Solubility of the dust                   dustsolub  =', dustsolub
472      ENDIF
473
474      ! Dust input from the atmosphere
475      ! ------------------------------
476      IF( bdustfer ) THEN
477         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
478         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
479         CALL iom_open ( 'dust.orca.nc', numdust )
480         DO jm = 1, jpmois
481            CALL iom_get( numdust, jpdom_data, 'dust', zdustmo(:,:,jm), jm )
482         END DO
483         CALL iom_close( numdust )
484      ELSE
485         zdustmo(:,:,:) = 0.e0
486         dust(:,:) = 0.0
487      ENDIF
488
489      ! Nutrient input from rivers
490      ! --------------------------
491      IF( briver ) THEN
492         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
493         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
494         CALL iom_open ( 'river.orca.nc', numriv )
495         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
496         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
497         CALL iom_close( numriv )
498      ELSE
499         river   (:,:) = 0.e0
500         riverdoc(:,:) = 0.e0
501      endif
502
503      ! Nutrient input from dust
504      ! ------------------------
505      IF( bndepo ) THEN
506         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
507         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
508         CALL iom_open ( 'ndeposition.orca.nc', numdep )
509         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
510         CALL iom_close( numdep )
511      ELSE
512         ndepo(:,:) = 0.e0
513      ENDIF
514
515      ! Coastal and island masks
516      ! ------------------------
517      IF( bsedinput ) THEN     
518         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron'
519         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
520         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
521         CALL iom_open ( 'bathy.orca.nc', numbath )
522         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
523         CALL iom_close( numbath )
524         !
525         DO jk = 1, 5
526            DO jj = 2, jpjm1
527               DO ji = 2, jpim1
528                  IF( tmask(ji,jj,jk) /= 0. ) THEN
529                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
530                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
531                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
532                  ENDIF
533               END DO
534            END DO
535         END DO
536         DO jk = 1, jpk
537            DO jj = 1, jpj
538               DO ji = 1, jpi
539                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
540                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
541                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
542               END DO
543            END DO
544         END DO
545      ELSE
546         cmask(:,:,:) = 0.e0
547      ENDIF
548
549      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
550
551
552      ! total atmospheric supply of Si
553      ! ------------------------------
554      sumdepsi = 0.e0
555      DO jm = 1, jpmois
556         DO jj = 2, jpjm1
557            DO ji = 2, jpim1
558               sumdepsi = sumdepsi + zdustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
559                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
560            END DO
561         END DO
562      END DO
563      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
564
565      ! N/P and Si releases due to coastal rivers
566      ! -----------------------------------------
567      DO jj = 1, jpj
568         DO ji = 1, jpi
569            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
570            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
571            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
572            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
573         END DO
574      END DO
575      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
576      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
577
578      rivpo4input=0.e0
579      rivalkinput=0.e0
580      nitdepinput=0.e0
581      DO jj = 2 , jpjm1
582         DO ji = 2, jpim1
583            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
584            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
585            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
586            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
587         END DO
588     END DO
589      IF( lk_mpp ) THEN
590         CALL mpp_sum( rivpo4input )  ! sum over the global domain
591         CALL mpp_sum( rivalkinput )  ! sum over the global domain
592         CALL mpp_sum( nitdepinput )  ! sum over the global domain
593      ENDIF
594
595
596      ! Coastal supply of iron
597      ! -------------------------
598      DO jk = 1, jpkm1
599         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
600      END DO
601      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
602
603
604   END SUBROUTINE p4z_sed_init
605
606
607
608#else
609   !!======================================================================
610   !!  Dummy module :                                   No PISCES bio-model
611   !!======================================================================
612CONTAINS
613   SUBROUTINE p4z_sed                         ! Empty routine
614   END SUBROUTINE p4z_sed
615#endif 
616
617   !!======================================================================
618END MODULE  p4zsed
Note: See TracBrowser for help on using the repository browser.