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

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

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