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

Last change on this file since 1715 was 1678, checked in by cetlod, 15 years ago

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