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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsed.F90 @ 2082

Last change on this file since 2082 was 2082, checked in by cetlod, 14 years ago

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