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

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

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