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.
p5zsms.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z – NEMO

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90 @ 7190

Last change on this file since 7190 was 7190, checked in by aumont, 7 years ago

various important bug fixes
  • Property svn:executable set to *
File size: 27.6 KB
Line 
1MODULE p5zsms
2   !!======================================================================
3   !!                         ***  MODULE p5zsms  ***
4   !! TOP :   PISCES Source Minus Sink manager
5   !!======================================================================
6   !! History :   1.0  !  2004-03 (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
9   !!----------------------------------------------------------------------
10#if defined key_pisces_quota
11   !!----------------------------------------------------------------------
12   !!   'key_pisces'                                       PISCES bio-model
13   !!----------------------------------------------------------------------
14   !!   p5zsms        :  Time loop of passive tracers sms
15   !!----------------------------------------------------------------------
16   USE oce_trc         !  shared variables between ocean and passive tracers
17   USE trc             !  passive tracers common variables
18   USE trcdta
19   USE sms_pisces      !  PISCES Source Minus Sink variables
20   USE p5zbio          !  Biological model
21   USE p4zche          !  Chemical model
22   USE p4zlys          !  Calcite saturation
23   USE p4zflx          !  Gas exchange
24   USE p4zsbc          !  External source of nutrients
25   USE p5zsed          !  Sedimentation
26   USE p4zint          !  time interpolation
27   USE p5zrem          !  remineralisation
28   USE iom             !  I/O manager
29   USE trd_oce         !  Ocean trends variables
30   USE trdtrc          !  TOP trends variables
31   USE sedmodel        !  Sediment model
32   USE prtctl_trc      !  print control for debugging
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   p5z_sms_init    ! called in p5zsms.F90
38   PUBLIC   p5z_sms    ! called in p5zsms.F90
39
40   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget
41   REAL(wp) :: xfact1, xfact2, xfact3
42   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget
43
44   !!* Array used to indicate negative tracer values
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
46
47   !! * Substitutions
48#  include "top_substitute.h90"
49
50   !!----------------------------------------------------------------------
51   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
52   !! $Id: p4zsms.F90 3320 2012-03-05 16:37:52Z cetlod $
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE p5z_sms( kt )
59      !!---------------------------------------------------------------------
60      !!                     ***  ROUTINE p5z_sms  ***
61      !!
62      !! ** Purpose :   Managment of the call to Biological sources and sinks
63      !!              routines of PISCES bio-model
64      !!
65      !! ** Method  : - at each new day ...
66      !!              - several calls of bio and sed ???
67      !!              - ...
68      !!---------------------------------------------------------------------
69      !
70      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
71      !!
72      INTEGER ::   ji,jj,jk,knt, jn, jl
73      REAL(wp) ::  ztra
74#if defined key_kriest
75      REAL(wp) ::  zcoef1, zcoef2
76#endif
77      CHARACTER (len=25) :: charout
78      !!---------------------------------------------------------------------
79      !
80      IF( nn_timing == 1 )  CALL timing_start('p5z_sms')
81      !
82      IF( kt == nittrc000 ) THEN
83        !
84        ALLOCATE( xnegtr(jpi,jpj,jpk) )
85        !
86        CALL p4z_che                              ! initialize the chemical constants
87        !
88        IF( .NOT. ln_rsttr ) THEN  ;   CALL ahini_for_at(hi)   !  set PH at kt=nit000
89        ELSE                       ;   CALL p5z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
90        ENDIF
91        !
92      ENDIF
93      !
94      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
95      !
96      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN   ;    rfact = rdttrc(1)     !  at nittrc000
97      ELSEIF( kt <= nittrc000 + nn_dttrc )                          THEN   ;    rfact = 2. * rdttrc(1)   ! at nittrc000 or nittrc000+nn_dttrc (Leapfrog)
98      ENDIF
99      !
100      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
101         rfactr  = 1. / rfact
102         rfact2  = rfact / FLOAT( nrdttrc )
103         rfact2r = 1. / rfact2
104         xstep = rfact2 / rday         ! Time step duration for biology
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1)
107         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
108         IF(lwp) WRITE(numout,*)
109      ENDIF
110
111      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
112         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
113            trb(:,:,:,jn) = trn(:,:,:,jn)
114         END DO
115      ENDIF
116      !
117      IF( ndayflxtr /= nday_year ) THEN      ! New days
118         !
119         ndayflxtr = nday_year
120
121         IF(lwp) write(numout,*)
122         IF(lwp) write(numout,*) ' New chemical constants and various rates for biogeochemistry at new day : ', nday_year
123         IF(lwp) write(numout,*) '~~~~~~'
124
125         CALL p4z_che              ! computation of chemical constants
126         CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
127         !
128      ENDIF
129
130      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients
131
132      DO knt = 1, nrdttrc          ! Potential time splitting if requested
133         !
134         CALL p5z_bio (kt, knt)    ! Biology
135         CALL p4z_lys( kt, knt )   ! Compute CaCO3 saturation
136         CALL p4z_flx( kt, knt )   ! Compute surface fluxes
137         CALL p5z_sed (kt, knt)    ! Sedimentation
138         !
139         xnegtr(:,:,:) = 1.e0
140         DO jn = jp_pcs0, jp_pcs1
141            DO jk = 1, jpk
142               DO jj = 1, jpj
143                  DO ji = 1, jpi
144                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
145                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
146                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
147                     ENDIF
148                 END DO
149               END DO
150            END DO
151         END DO
152         !                                ! where at least 1 tracer concentration becomes negative
153         !                                !
154         DO jn = jp_pcs0, jp_pcs1
155           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
156         END DO
157        !
158         DO jn = jp_pcs0, jp_pcs1
159            tra(:,:,:,jn) = 0._wp
160         END DO
161         !
162         IF( ln_top_euler ) THEN
163            DO jn = jp_pcs0, jp_pcs1
164               trn(:,:,:,jn) = trb(:,:,:,jn)
165            END DO
166         ENDIF
167      END DO
168
169#if defined key_kriest
170      !
171      zcoef1 = 1.e0 / xkr_massp
172      zcoef2 = 1.e0 / xkr_massp / 1.1
173      DO jk = 1,jpkm1
174         trb(:,:,jk,jpnum) = MAX(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  )
175         trb(:,:,jk,jpnum) = MIN(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef2              )
176      END DO
177      !
178#endif
179      !
180      !
181      IF( l_trdtrc )  THEN
182         DO jn = jp_pcs0, jp_pcs1
183           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
184         END DO
185      ENDIF
186      !
187      IF( lk_sed ) THEN 
188         !
189         CALL sed_model( kt )     !  Main program of Sediment model
190         !
191         DO jn = jp_pcs0, jp_pcs1
192           CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
193         END DO
194         !
195      ENDIF
196      !
197      IF( lrst_trc )  CALL p5z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file
198      !
199      IF( lk_iomput .OR. ln_check_mass )  CALL p5z_chk_mass( kt ) ! Mass conservation checking
200
201      IF ( lwm .AND. kt == nittrc000 ) CALL FLUSH    ( numonp )     ! flush output namelist PISCES
202      IF( nn_timing == 1 )  CALL timing_stop('p5z_sms')
203      !
204      !
205   END SUBROUTINE p5z_sms
206
207   SUBROUTINE p5z_sms_init
208      !!----------------------------------------------------------------------
209      !!                     ***  p5z_sms_init  *** 
210      !!
211      !! ** Purpose :   read PISCES namelist
212      !!
213      !! ** input   :   file 'namelist.trc.s' containing the following
214      !!             namelist: natext, natbio, natsms
215      !!                       natkriest ("key_kriest")
216      !!----------------------------------------------------------------------
217#if defined key_ligand
218      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, no3rat3, po4rat3, wsbio2,    &
219      &                   wsbio2max, wsbio2scale, niter1max, niter2max, wfep, ldocp, ldocz, lthet
220#else
221      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, no3rat3, po4rat3, wsbio2,    &
222      &                   wsbio2max, wsbio2scale, niter1max, niter2max
223#endif
224#if defined key_kriest
225      NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_ncontent, xkr_mass_min, xkr_mass_max
226#endif
227      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
228      NAMELIST/nampismass/ ln_check_mass
229      INTEGER :: ios                 ! Local integer output status for namelist read
230      !!----------------------------------------------------------------------
231
232      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
233      READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
234901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
235
236      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
237      READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
238902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
239      IF(lwm) WRITE ( numonp, nampisbio )
240
241      IF(lwp) THEN                         ! control print
242         WRITE(numout,*) ' Namelist : nampisbio'
243         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
244         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
245         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
246         WRITE(numout,*) '    N/C in zooplankton                        no3rat3   =', no3rat3
247         WRITE(numout,*) '    P/C in zooplankton                        po4rat3   =', po4rat3
248         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
249         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
250         WRITE(numout,*) '    Big particles maximum sinking speed       wsbio2max =', wsbio2max
251         WRITE(numout,*) '    Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
252         WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
253         WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
254#if defined key_ligand
255         WRITE(numout,*) '    FeP sinking speed                         wfep       =', wfep
256#endif
257      ENDIF
258
259#if defined key_kriest
260
261      !                               ! nampiskrp : kriest parameters
262      !                               ! -----------------------------
263      REWIND( numnatp_ref )              ! Namelist nampiskrp in reference namelist : Pisces Kriest
264      READ  ( numnatp_ref, nampiskrp, IOSTAT = ios, ERR = 903)
265903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrp in reference namelist', lwp )
266
267      REWIND( numnatp_cfg )              ! Namelist nampiskrp in configuration namelist : Pisces Kriest
268      READ  ( numnatp_cfg, nampiskrp, IOSTAT = ios, ERR = 904 )
269904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrp in configuration namelist', lwp )
270      IF(lwm) WRITE ( numonp, nampiskrp )
271
272      IF(lwp) THEN
273         WRITE(numout,*)
274         WRITE(numout,*) ' Namelist : nampiskrp'
275         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
276         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
277         WRITE(numout,*) '    N content factor                         xkr_ncontent = ', xkr_ncontent
278         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
279         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
280         WRITE(numout,*)
281     ENDIF
282
283
284     ! Computation of some variables
285     xkr_massp = xkr_ncontent * 7.625 * xkr_mass_min**xkr_zeta
286
287#endif
288
289      REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
290      READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
291905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
292
293      REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
294      READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
295906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
296      IF(lwm) WRITE ( numonp, nampisdmp )
297
298      IF(lwp) THEN                         ! control print
299         WRITE(numout,*)
300         WRITE(numout,*) ' Namelist : nampisdmp'
301         WRITE(numout,*) '    Relaxation of tracer to glodap mean value             ln_pisdmp      =', ln_pisdmp
302         WRITE(numout,*) '    Frequency of Relaxation                               nn_pisdmp      =', nn_pisdmp
303         WRITE(numout,*) ' '
304      ENDIF
305
306      REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
307      READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
308907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
309
310      REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check
311      READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
312908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
313      IF(lwm) WRITE ( numonp, nampismass )
314
315      IF(lwp) THEN                         ! control print
316         WRITE(numout,*) ' '
317         WRITE(numout,*) ' Namelist parameter for mass conservation checking'
318         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
319         WRITE(numout,*) '    Flag to check mass conservation of NO3/Si/TALK ln_check_mass = ', ln_check_mass
320      ENDIF
321
322   END SUBROUTINE p5z_sms_init
323
324   SUBROUTINE p5z_rst( kt, cdrw )
325      !!---------------------------------------------------------------------
326      !!                   ***  ROUTINE p5z_rst  ***
327      !!
328      !!  ** Purpose : Read or write variables in restart file:
329      !!
330      !!  WRITE(READ) mode:
331      !!       kt        : number of time step since the begining of the experiment at the
332      !!                   end of the current(previous) run
333      !!---------------------------------------------------------------------
334      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
335      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
336      !!---------------------------------------------------------------------
337
338      IF( TRIM(cdrw) == 'READ' ) THEN
339         !
340         IF(lwp) WRITE(numout,*)
341         IF(lwp) WRITE(numout,*) ' p5z_rst : Read specific variables from pisces model '
342         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
343         !
344         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
345            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
346         ELSE
347            CALL ahini_for_at(hi)
348         ENDIF
349         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
350         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
351            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
352         ELSE
353            xksimax(:,:) = xksi(:,:)
354         ENDIF
355         !
356         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
357            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
358         ELSE
359            t_oce_co2_flx_cum = 0._wp
360         ENDIF
361         !
362         IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
363            CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sized(:,:,:)  )
364            CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sized(:,:,:)  )
365            CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
366         ELSE
367            sizep(:,:,:) = 1.
368            sizen(:,:,:) = 1.
369            sized(:,:,:) = 1.
370         ENDIF
371         !
372      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
373         IF( kt == nitrst ) THEN
374            IF(lwp) WRITE(numout,*)
375            IF(lwp) WRITE(numout,*) 'p5z_rst : write pisces restart file  kt =', kt
376            IF(lwp) WRITE(numout,*) '~~~~~~~'
377         ENDIF
378         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
379         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
380         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
381         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum )
382         CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sized(:,:,:) )
383         CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sized(:,:,:) )
384         CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )
385      ENDIF
386      !
387   END SUBROUTINE p5z_rst
388
389   SUBROUTINE p4z_dmp( kt )
390      !!----------------------------------------------------------------------
391      !!                    ***  p4z_dmp  ***
392      !!
393      !! ** purpose  : Relaxation of some tracers
394      !!----------------------------------------------------------------------
395      !
396      INTEGER, INTENT( in )  ::     kt ! time step
397      !
398      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
399      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
400      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
401      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
402      !
403      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
404      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
405      !!---------------------------------------------------------------------
406
407
408      IF(lwp)  WRITE(numout,*)
409      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
410      IF(lwp)  WRITE(numout,*)
411
412      IF( cp_cfg == "orca" .AND. .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
413         !                                                    ! --------------------------- !
414         ! set total alkalinity, phosphate, nitrate & silicate
415         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6             
416
417         zalksumn = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
418         zpo4sumn = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
419         zno3sumn = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
420         zsilsumn = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
421 
422         IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
423         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
424
425         IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
426         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
427
428         IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
429         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
430
431         IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
432         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
433         !
434         IF( .NOT. ln_top_euler ) THEN
435            zalksumb = glob_sum( trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
436            zpo4sumb = glob_sum( trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
437            zno3sumb = glob_sum( trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
438            zsilsumb = glob_sum( trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
439
440            IF(lwp) WRITE(numout,*) ' '
441            IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
442            trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
443
444            IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
445            trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
446
447            IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
448            trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
449
450            IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
451            trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
452        ENDIF
453        !
454      ENDIF
455
456   END SUBROUTINE p4z_dmp
457
458
459   SUBROUTINE p5z_chk_mass( kt )
460      !!----------------------------------------------------------------------
461      !!                  ***  ROUTINE p5z_chk_mass  ***
462      !!
463      !! ** Purpose :  Mass conservation check
464      !!
465      !!---------------------------------------------------------------------
466      !
467      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
468      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
469      CHARACTER(LEN=100)   ::   cltxt
470      INTEGER :: jk
471      !!
472      !!---------------------------------------------------------------------
473
474      IF( kt == nittrc000 ) THEN
475         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
476            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
477            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
478            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
479            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
480            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
481            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
482            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
483            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
484            IF( lwp ) WRITE(numnut,*)
485         ENDIF
486      ENDIF
487
488      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
489         !   Compute the budget of NO3, ALK, Si, Fer
490         no3budget = glob_sum( (   trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)  &
491            &                    + trn(:,:,:,jpnph) + trn(:,:,:,jpndi)  &
492            &                    + (trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)) * no3rat3  &
493            &                    + trn(:,:,:,jppon) + trn(:,:,:,jpnpi) &
494#if ! defined key_kriest
495            &                    + trn(:,:,:,jpgon)                     &
496#endif
497            &                    + trn(:,:,:,jpdon)                     ) * cvol(:,:,:)  ) 
498         !
499         no3budget = no3budget / areatot
500         CALL iom_put( "pno3tot", no3budget )
501      ENDIF
502      !
503      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
504         po4budget = glob_sum( (   trn(:,:,:,jppo4)                     &
505            &                    + trn(:,:,:,jppph) + trn(:,:,:,jppdi)  &
506            &                    + (trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)) * po4rat3  &
507            &                    + trn(:,:,:,jppop) + trn(:,:,:,jpppi)  &
508#if ! defined key_kriest
509            &                    + trn(:,:,:,jpgop)                     &
510#endif
511            &                    + trn(:,:,:,jpdop)                     ) * cvol(:,:,:)  )
512         po4budget = po4budget / areatot
513         CALL iom_put( "ppo4tot", po4budget )
514      ENDIF
515      !
516      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
517         silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
518            &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
519         !
520         silbudget = silbudget / areatot
521         CALL iom_put( "psiltot", silbudget )
522      ENDIF
523      !
524      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
525         !
526         alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
527            &                    + trn(:,:,:,jptal)                     &
528            &                    + trn(:,:,:,jpcal) * 2.                ) * cvol(:,:,:)  )
529         !
530         alkbudget = alkbudget / areatot
531         CALL iom_put( "palktot", alkbudget )
532      ENDIF
533      !
534      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
535         ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
536            &                    + trn(:,:,:,jpdfe) + trn(:,:,:,jppfe)  &
537#if ! defined key_kriest
538            &                    + trn(:,:,:,jpbfe)                     &
539#endif
540#if defined key_ligand
541            &                    + trn(:,:,:,jpfep)                     &
542#endif
543            &                    + trn(:,:,:,jpsfe)                     &
544            &                    + trn(:,:,:,jpzoo) * ferat3            &
545            &                    + trn(:,:,:,jpmes) * ferat3            ) * cvol(:,:,:)  )
546         !
547         ferbudget = ferbudget / areatot
548         CALL iom_put( "pfertot", ferbudget )
549      ENDIF
550      !
551
552      ! Global budget of N SMS : denitrification in the water column and in the sediment
553      !                          nitrogen fixation by the diazotrophs
554      ! --------------------------------------------------------------------------------
555      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
556         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
557         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3
558      ENDIF
559      !
560      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
561         zrdenittot   = glob_sum ( denitrc(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
562         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3
563      ENDIF
564      !
565      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
566         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) )
567         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments
568      ENDIF
569
570      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
571         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
572         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
573         tpp            = tpp           * 1000. * xfact1
574         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
575         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
576         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
577             &                                no3budget * rno3 * 1.e+06, &
578             &                                po4budget * po4r * 1.e+06, &
579             &                                silbudget        * 1.e+06, &
580             &                                ferbudget        * 1.e+09
581         !
582         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
583         &                             zrdenittot  * xfact2  , &
584         &                             zsdenittot  * xfact2
585
586      ENDIF
587       !
588 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
589 9100  FORMAT(i8,5e18.10)
590 9200  FORMAT(i8,3f10.5)
591       !
592   END SUBROUTINE p5z_chk_mass
593
594#else
595   !!======================================================================
596   !!  Dummy module :                                   No PISCES bio-model
597   !!======================================================================
598CONTAINS
599   SUBROUTINE p5z_sms( kt )                   ! Empty routine
600      INTEGER, INTENT( in ) ::   kt
601      WRITE(*,*) 'p5z_sms: You should not have seen this print! error?', kt
602   END SUBROUTINE p5z_sms
603#endif 
604
605   !!======================================================================
606END MODULE p5zsms 
Note: See TracBrowser for help on using the repository browser.