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.
p4zsms.F90 in NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90 @ 12759

Last change on this file since 12759 was 12759, checked in by aumont, 4 years ago

make parameterizations in PISCES-operationnal more similar to thos of PISCES-QUOTA (prey switching, optimal allocation, size, ...)
  • Property svn:keywords set to Id
File size: 27.5 KB
Line 
1MODULE p4zsms
2   !!======================================================================
3   !!                         ***  MODULE p4zsms  ***
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   !!----------------------------------------------------------------------
9   !!   p4z_sms        : Time loop of passive tracers sms
10   !!----------------------------------------------------------------------
11   USE oce_trc         ! shared variables between ocean and passive tracers
12   USE trc             ! passive tracers common variables
13   USE trcdta          !
14   USE sms_pisces      ! PISCES Source Minus Sink variables
15   USE p4zbio          ! Biological model
16   USE p4zche          ! Chemical model
17   USE p4zlys          ! Calcite saturation
18   USE p4zflx          ! Gas exchange
19   USE p4zsbc          ! External source of nutrients
20   USE p4zsed          ! Sedimentation
21   USE p4zint          ! time interpolation
22   USE p4zrem          ! remineralisation
23   USE iom             ! I/O manager
24   USE trd_oce         ! Ocean trends variables
25   USE trdtrc          ! TOP trends variables
26   USE sedmodel        ! Sediment model
27   USE prtctl_trc      ! print control for debugging
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   p4z_sms_init   ! called in trcini_pisces.F90
33   PUBLIC   p4z_sms        ! called in trcsms_pisces.F90
34
35   INTEGER ::    numco2, numnut, numnit      ! logical unit for co2 budget
36   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget ! total budget of the different conservative elements
37   REAL(wp) ::   xfact1, xfact2, xfact3
38
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     ! Array used to indicate negative tracer values
40
41   !!----------------------------------------------------------------------
42   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE p4z_sms( kt )
49      !!---------------------------------------------------------------------
50      !!                     ***  ROUTINE p4z_sms  ***
51      !!
52      !! ** Purpose :   Managment of the call to Biological sources and sinks
53      !!              routines of PISCES bio-model
54      !!
55      !! ** Method  : - calls the various SMS subroutines
56      !!              - calls the sediment module (if ln_sediment) 
57      !!              - several calls of bio and sed (possible time-splitting)
58      !!              - handles the potential negative concentrations (xnegtr)
59      !!---------------------------------------------------------------------
60      !
61      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
62      !!
63      INTEGER ::   ji, jj, jk, jnt, jn, jl
64      REAL(wp) ::  ztra
65      CHARACTER (len=25) :: charout
66      !!---------------------------------------------------------------------
67      !
68      IF( ln_timing )   CALL timing_start('p4z_sms')
69      !
70      IF( kt == nittrc000 ) THEN
71        !
72        ALLOCATE( xnegtr(jpi,jpj,jpk) )
73        !
74        IF( .NOT. ln_rsttr ) THEN
75            CALL p4z_che            ! initialize the chemical constants
76            CALL ahini_for_at(hi)   !  set PH at kt=nit000
77            t_oce_co2_flx_cum = 0._wp
78        ELSE
79            CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
80        ENDIF
81        !
82      ENDIF
83      !
84      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
85      !
86      rfact = r2dttrc
87      !
88      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
89         rfactr  = 1. / rfact
90         rfact2  = rfact / REAL( nrdttrc, wp )
91         rfact2r = 1. / rfact2
92         xstep = rfact2 / rday         ! Time step duration for biology
93         IF(lwp) WRITE(numout,*) 
94         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
95         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
96         IF(lwp) WRITE(numout,*)
97      ENDIF
98
99      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
100         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
101            trb(:,:,:,jn) = trn(:,:,:,jn)
102         END DO
103      ENDIF
104      !
105      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients
106      !
107#if ! defined key_sed_off
108      CALL p4z_che              ! computation of chemical constants
109      CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
110      !
111      DO jnt = 1, nrdttrc          ! Potential time splitting if requested
112         !
113         CALL p4z_bio( kt, jnt )   ! Biological SMS
114         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation and dissolution
115         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
116         CALL p4z_flx( kt, jnt )   ! Compute surface air-sea fluxes
117         !
118         ! Handling of the negative concentrations
119         ! The biological SMS may generate negative concentrations
120         ! Trends are tested at each grid cell. If a negative concentrations
121         ! is created at a grid cell, all the sources and sinks at that grid
122         ! cell are scale to avoid that negative concentration. This approach
123         ! is quite simplistic but it conserves mass.
124         ! ------------------------------------------------------------------
125         xnegtr(:,:,:) = 1.e0
126         DO jn = jp_pcs0, jp_pcs1
127            DO jk = 1, jpk
128               DO jj = 1, jpj
129                  DO ji = 1, jpi
130                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
131                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
132                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
133                     ENDIF
134                 END DO
135               END DO
136            END DO
137         END DO
138         !                                ! where at least 1 tracer concentration becomes negative
139         !                                !
140         ! Concentrations are updated
141         DO jn = jp_pcs0, jp_pcs1
142           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
143         END DO
144        !
145        ! Trends are are reset to 0
146         DO jn = jp_pcs0, jp_pcs1
147            tra(:,:,:,jn) = 0._wp
148         END DO
149         !
150         IF( ln_top_euler ) THEN
151            DO jn = jp_pcs0, jp_pcs1
152               trn(:,:,:,jn) = trb(:,:,:,jn)
153            END DO
154         ENDIF
155      END DO
156
157      !
158      IF( l_trdtrc ) THEN
159         DO jn = jp_pcs0, jp_pcs1
160           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
161         END DO
162      END IF
163#endif
164      !
165      ! If ln_sediment is set to .true. then the sediment module is called
166      IF( ln_sediment ) THEN 
167         !
168         CALL sed_model( kt )     !  Main program of Sediment model
169         !
170         IF( ln_top_euler ) THEN
171            DO jn = jp_pcs0, jp_pcs1
172               trn(:,:,:,jn) = trb(:,:,:,jn)
173            END DO
174         ENDIF
175         !
176      ENDIF
177      !
178      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file
179      !
180
181      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt )    ! Mass conservation checking
182
183      IF( lwm .AND. kt == nittrc000    )  CALL FLUSH( numonp )       ! flush output namelist PISCES
184      !
185      IF( ln_timing )  CALL timing_stop('p4z_sms')
186      !
187   END SUBROUTINE p4z_sms
188
189
190   SUBROUTINE p4z_sms_init
191      !!----------------------------------------------------------------------
192      !!                     ***  p4z_sms_init  *** 
193      !!
194      !! ** Purpose :   read the general PISCES namelist
195      !!
196      !! ** input   :   file 'namelist_pisces' containing the following
197      !!                namelist: nampisbio, nampisdmp, nampismass
198      !!----------------------------------------------------------------------
199      INTEGER :: ios                 ! Local integer output status for namelist read
200      !!
201      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
202         &                   ldocp, ldocz, lthet, no3rat3, po4rat3
203         !
204      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
205      NAMELIST/nampismass/ ln_check_mass
206      !!----------------------------------------------------------------------
207      !
208      IF(lwp) THEN
209         WRITE(numout,*)
210         WRITE(numout,*) 'p4z_sms_init : PISCES initialization'
211         WRITE(numout,*) '~~~~~~~~~~~~'
212      ENDIF
213
214      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
215      READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
216901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbio in reference namelist' )
217      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
218      READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
219902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbio in configuration namelist' )
220      IF(lwm) WRITE( numonp, nampisbio )
221      !
222      IF(lwp) THEN                         ! control print
223         WRITE(numout,*) '   Namelist : nampisbio'
224         WRITE(numout,*) '      frequency for the biology                 nrdttrc     =', nrdttrc
225         WRITE(numout,*) '      POC sinking speed                         wsbio       =', wsbio
226         WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort 
227         IF( ln_p5z ) THEN
228            WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3
229            WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3
230         ENDIF
231         WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3
232         WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2
233         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
234         WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
235         IF( ln_ligand ) THEN
236            IF( ln_p4z ) THEN
237               WRITE(numout,*) '      Phyto ligand production per unit doc           ldocp  =', ldocp
238               WRITE(numout,*) '      Zoo ligand production per unit doc             ldocz  =', ldocz
239               WRITE(numout,*) '      Proportional loss of ligands due to Fe uptake  lthet  =', lthet
240            ENDIF
241         ENDIF
242      ENDIF
243
244
245      REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
246      READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
247905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisdmp in reference namelist' )
248      REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
249      READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
250906   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisdmp in configuration namelist' )
251      IF(lwm) WRITE( numonp, nampisdmp )
252      !
253      IF(lwp) THEN                         ! control print
254         WRITE(numout,*)
255         WRITE(numout,*) '   Namelist : nampisdmp --- relaxation to GLODAP'
256         WRITE(numout,*) '      Relaxation of tracer to glodap mean value   ln_pisdmp =', ln_pisdmp
257         WRITE(numout,*) '      Frequency of Relaxation                     nn_pisdmp =', nn_pisdmp
258      ENDIF
259
260      REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
261      READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
262907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampismass in reference namelist' )
263      REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check
264      READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
265908   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampismass in configuration namelist' )
266      IF(lwm) WRITE( numonp, nampismass )
267
268      IF(lwp) THEN                         ! control print
269         WRITE(numout,*)
270         WRITE(numout,*) '   Namelist : nampismass  --- mass conservation checking'
271         WRITE(numout,*) '      Flag to check mass conservation of NO3/Si/TALK   ln_check_mass = ', ln_check_mass
272      ENDIF
273      !
274   END SUBROUTINE p4z_sms_init
275
276
277   SUBROUTINE p4z_rst( kt, cdrw )
278      !!---------------------------------------------------------------------
279      !!                   ***  ROUTINE p4z_rst  ***
280      !!
281      !!  ** Purpose : Read or write specific PISCES variables in restart file:
282      !!
283      !!  WRITE(READ) mode:
284      !!       kt        : number of time step since the begining of the experiment at the
285      !!                   end of the current(previous) run
286      !!---------------------------------------------------------------------
287      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
288      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
289      !!---------------------------------------------------------------------
290      !
291      IF( TRIM(cdrw) == 'READ' ) THEN
292         !
293         ! Read the specific variable of PISCES
294         IF(lwp) WRITE(numout,*)
295         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
296         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
297         !
298         ! Read the pH. If not in the restart file, then it is initialized from
299         ! the initial conditions
300         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
301            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
302         ELSE
303            CALL p4z_che                              ! initialize the chemical constants
304            CALL ahini_for_at(hi)
305         ENDIF
306
307         ! Read the Si half saturation constant and the maximum Silica concentration
308         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
309         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
310            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
311         ELSE
312            xksimax(:,:) = xksi(:,:)
313         ENDIF
314
315         ! Read the cumulative total flux. If not in the restart file, it is set to 0         
316         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
317            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
318         ELSE
319            t_oce_co2_flx_cum = 0._wp
320         ENDIF
321         !
322         ! PISCES size proxy
323         IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
324            CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
325         ELSE
326            sized(:,:,:) = 1.
327         ENDIF
328         IF( iom_varid( numrtr, 'sizen', ldstop = .FALSE. ) > 0 ) THEN
329            CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
330         ELSE
331            sizen(:,:,:) = 1.
332         ENDIF
333
334         ! PISCES-QUOTA specific part
335         IF( ln_p5z ) THEN
336            ! Read the size of the different phytoplankton groups
337            ! If not in the restart file, they are set to 1
338            IF( iom_varid( numrtr, 'sizep', ldstop = .FALSE. ) > 0 ) THEN
339               CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
340            ELSE
341               sizep(:,:,:) = 1.
342            ENDIF
343        ENDIF
344        !
345      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
346         ! write the specific variables of PISCES
347         IF( kt == nitrst ) THEN
348            IF(lwp) WRITE(numout,*)
349            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
350            IF(lwp) WRITE(numout,*) '~~~~~~~'
351         ENDIF
352         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )           ! pH
353         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )    ! Si 1/2 saturation constant
354         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
355         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux
356         CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
357         CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms
358         IF( ln_p5z ) THEN
359            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )  ! Size of picophytoplankton
360         ENDIF
361      ENDIF
362      !
363   END SUBROUTINE p4z_rst
364
365
366   SUBROUTINE p4z_dmp( kt )
367      !!----------------------------------------------------------------------
368      !!                    ***  p4z_dmp  ***
369      !!
370      !! ** purpose  : Relaxation of the total budget of some elements
371      !!               This routine avoids the model to drift far from the
372      !!               observed content in various elements
373      !!               Elements that may be relaxed : Alk, P, N, Si
374      !!----------------------------------------------------------------------
375      !
376      INTEGER, INTENT( in )  ::     kt ! time step
377      !
378      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
379      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphate
380      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
381      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
382      !
383      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
384      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
385      !!---------------------------------------------------------------------
386
387      IF(lwp)  WRITE(numout,*)
388      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
389      IF(lwp)  WRITE(numout,*)
390
391      IF( cn_cfg == "ORCA" .OR. cn_cfg == "orca") THEN
392         IF( .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
393            !                                                ! --------------------------- !
394            ! set total alkalinity, phosphate, nitrate & silicate
395            zarea          = 1._wp / glob_sum( 'p4zsms', cvol(:,:,:) ) * 1e6             
396
397            zalksumn = glob_sum( 'p4zsms', trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
398            zpo4sumn = glob_sum( 'p4zsms', trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
399            zno3sumn = glob_sum( 'p4zsms', trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
400            zsilsumn = glob_sum( 'p4zsms', trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
401 
402            ! Correct the trn mean content of alkalinity
403            IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
404            trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
405
406            ! Correct the trn mean content of PO4
407            IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
408            trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
409
410            ! Correct the trn mean content of NO3
411            IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
412            trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
413
414            ! Correct the trn mean content of SiO3
415            IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
416            trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
417            !
418            !
419            IF( .NOT. ln_top_euler ) THEN
420               zalksumb = glob_sum( 'p4zsms', trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
421               zpo4sumb = glob_sum( 'p4zsms', trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
422               zno3sumb = glob_sum( 'p4zsms', trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
423               zsilsumb = glob_sum( 'p4zsms', trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
424 
425               IF(lwp) WRITE(numout,*) ' '
426               ! Correct the trb mean content of alkalinity
427               IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
428               trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
429
430               ! Correct the trb mean content of PO4
431               IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
432               trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
433
434               ! Correct the trb mean content of NO3
435               IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
436               trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
437
438               ! Correct the trb mean content of SiO3
439               IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
440               trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
441           ENDIF
442        ENDIF
443        !
444      ENDIF
445        !
446   END SUBROUTINE p4z_dmp
447
448
449   SUBROUTINE p4z_chk_mass( kt )
450      !!----------------------------------------------------------------------
451      !!                  ***  ROUTINE p4z_chk_mass  ***
452      !!
453      !! ** Purpose :  Mass conservation check
454      !!
455      !!---------------------------------------------------------------------
456      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
457      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
458      CHARACTER(LEN=100)   ::   cltxt
459      INTEGER :: jk
460      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork
461      !!----------------------------------------------------------------------
462      !
463      IF( kt == nittrc000 ) THEN
464         xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
465         xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
466         xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
467         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
468            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
469            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
470            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
471            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
472            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
473            IF( lwp ) WRITE(numnut,*) 
474         ENDIF
475      ENDIF
476
477      ! Compute the budget of NO3
478      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
479         IF( ln_p4z ) THEN
480            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
481               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
482               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
483               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
484        ELSE
485            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4) + trn(:,:,:,jpnph)   &
486               &          +   trn(:,:,:,jpndi) + trn(:,:,:,jpnpi)                      & 
487               &          +   trn(:,:,:,jppon) + trn(:,:,:,jpgon) + trn(:,:,:,jpdon)   &
488               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * no3rat3 
489        ENDIF
490        !
491        no3budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
492        no3budget = no3budget / areatot
493        CALL iom_put( "pno3tot", no3budget )
494      ENDIF
495      !
496      ! Compute the budget of PO4
497      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
498         IF( ln_p4z ) THEN
499            zwork(:,:,:) =    trn(:,:,:,jppo4)                                         &
500               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
501               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
502               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
503        ELSE
504            zwork(:,:,:) =    trn(:,:,:,jppo4) + trn(:,:,:,jppph)                      &
505               &          +   trn(:,:,:,jppdi) + trn(:,:,:,jpppi)                      & 
506               &          +   trn(:,:,:,jppop) + trn(:,:,:,jpgop) + trn(:,:,:,jpdop)   &
507               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * po4rat3 
508        ENDIF
509        !
510        po4budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
511        po4budget = po4budget / areatot
512        CALL iom_put( "ppo4tot", po4budget )
513      ENDIF
514      !
515      ! Compute the budget of SiO3
516      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
517         zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
518         !
519         silbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
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         zwork(:,:,:) =  trn(:,:,:,jpno3) * rno3 + trn(:,:,:,jptal) + trn(:,:,:,jpcal) * 2.             
526         !
527         alkbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )         !
528         alkbudget = alkbudget / areatot
529         CALL iom_put( "palktot", alkbudget )
530      ENDIF
531      !
532      ! Compute the budget of Iron
533      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
534         zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
535            &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
536            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3   
537         !
538         ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
539         ferbudget = ferbudget / areatot
540         CALL iom_put( "pfertot", ferbudget )
541      ENDIF
542      !
543      ! Global budget of N SMS : denitrification in the water column and in the sediment
544      !                          nitrogen fixation by the diazotrophs
545      ! --------------------------------------------------------------------------------
546      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
547         znitrpottot  = glob_sum ( 'p4zsms', nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
548         CALL iom_put( "tnfix"  , znitrpottot * xfact3 )  ! Global  nitrogen fixation molC/l  to molN/m3
549      ENDIF
550      !
551      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
552         zrdenittot = glob_sum ( 'p4zsms', denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
553         zsdenittot = glob_sum ( 'p4zsms', sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )
554         CALL iom_put( "tdenit" , ( zrdenittot + zsdenittot ) * xfact3 )  ! Total denitrification molC/l to molN/m3
555      ENDIF
556      !
557      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
558         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( 'p4zsms', e1e2t(:,:) )
559         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
560         tpp            = tpp           * 1000. * xfact1
561         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
562         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
563         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
564             &                                no3budget * rno3 * 1.e+06, &
565             &                                po4budget * po4r * 1.e+06, &
566             &                                silbudget        * 1.e+06, &
567             &                                ferbudget        * 1.e+09
568         !
569         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
570            &                             zrdenittot  * xfact2  , &
571            &                             zsdenittot  * xfact2
572      ENDIF
573      !
574 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
575 9100  FORMAT(i8,5e18.10)
576 9200  FORMAT(i8,3f10.5)
577       !
578   END SUBROUTINE p4z_chk_mass
579
580   !!======================================================================
581END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.