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

Last change on this file since 12537 was 12537, checked in by aumont, 5 years ago

Comments in routines have been revised and significantly augmented
  • Property svn:keywords set to Id
File size: 27.3 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-QUOTA specific part
323         IF( ln_p5z ) THEN
324            ! Read the size of the different phytoplankton groups
325            ! If not in the restart file, they are set to 1
326            IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
327               CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
328               CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
329               CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
330            ELSE
331               sizep(:,:,:) = 1.
332               sizen(:,:,:) = 1.
333               sized(:,:,:) = 1.
334            ENDIF
335        ENDIF
336        !
337      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
338         ! write the specific variables of PISCES
339         IF( kt == nitrst ) THEN
340            IF(lwp) WRITE(numout,*)
341            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
342            IF(lwp) WRITE(numout,*) '~~~~~~~'
343         ENDIF
344         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )           ! pH
345         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )    ! Si 1/2 saturation constant
346         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
347         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux
348         IF( ln_p5z ) THEN
349            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )  ! Size of picophytoplankton
350            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
351            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms
352         ENDIF
353      ENDIF
354      !
355   END SUBROUTINE p4z_rst
356
357
358   SUBROUTINE p4z_dmp( kt )
359      !!----------------------------------------------------------------------
360      !!                    ***  p4z_dmp  ***
361      !!
362      !! ** purpose  : Relaxation of the total budget of some elements
363      !!               This routine avoids the model to drift far from the
364      !!               observed content in various elements
365      !!               Elements that may be relaxed : Alk, P, N, Si
366      !!----------------------------------------------------------------------
367      !
368      INTEGER, INTENT( in )  ::     kt ! time step
369      !
370      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
371      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphate
372      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
373      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
374      !
375      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
376      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
377      !!---------------------------------------------------------------------
378
379      IF(lwp)  WRITE(numout,*)
380      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
381      IF(lwp)  WRITE(numout,*)
382
383      IF( cn_cfg == "ORCA" .OR. cn_cfg == "orca") THEN
384         IF( .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
385            !                                                ! --------------------------- !
386            ! set total alkalinity, phosphate, nitrate & silicate
387            zarea          = 1._wp / glob_sum( 'p4zsms', cvol(:,:,:) ) * 1e6             
388
389            zalksumn = glob_sum( 'p4zsms', trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
390            zpo4sumn = glob_sum( 'p4zsms', trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
391            zno3sumn = glob_sum( 'p4zsms', trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
392            zsilsumn = glob_sum( 'p4zsms', trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
393 
394            ! Correct the trn mean content of alkalinity
395            IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
396            trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
397
398            ! Correct the trn mean content of PO4
399            IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
400            trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
401
402            ! Correct the trn mean content of NO3
403            IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
404            trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
405
406            ! Correct the trn mean content of SiO3
407            IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
408            trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
409            !
410            !
411            IF( .NOT. ln_top_euler ) THEN
412               zalksumb = glob_sum( 'p4zsms', trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
413               zpo4sumb = glob_sum( 'p4zsms', trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
414               zno3sumb = glob_sum( 'p4zsms', trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
415               zsilsumb = glob_sum( 'p4zsms', trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
416 
417               IF(lwp) WRITE(numout,*) ' '
418               ! Correct the trb mean content of alkalinity
419               IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
420               trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
421
422               ! Correct the trb mean content of PO4
423               IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
424               trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
425
426               ! Correct the trb mean content of NO3
427               IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
428               trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
429
430               ! Correct the trb mean content of SiO3
431               IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
432               trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
433           ENDIF
434        ENDIF
435        !
436      ENDIF
437        !
438   END SUBROUTINE p4z_dmp
439
440
441   SUBROUTINE p4z_chk_mass( kt )
442      !!----------------------------------------------------------------------
443      !!                  ***  ROUTINE p4z_chk_mass  ***
444      !!
445      !! ** Purpose :  Mass conservation check
446      !!
447      !!---------------------------------------------------------------------
448      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
449      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
450      CHARACTER(LEN=100)   ::   cltxt
451      INTEGER :: jk
452      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork
453      !!----------------------------------------------------------------------
454      !
455      IF( kt == nittrc000 ) THEN
456         xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
457         xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
458         xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
459         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
460            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
461            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
462            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
463            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
464            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
465            IF( lwp ) WRITE(numnut,*) 
466         ENDIF
467      ENDIF
468
469      ! Compute the budget of NO3
470      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
471         IF( ln_p4z ) THEN
472            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
473               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
474               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
475               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
476        ELSE
477            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4) + trn(:,:,:,jpnph)   &
478               &          +   trn(:,:,:,jpndi) + trn(:,:,:,jpnpi)                      & 
479               &          +   trn(:,:,:,jppon) + trn(:,:,:,jpgon) + trn(:,:,:,jpdon)   &
480               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * no3rat3 
481        ENDIF
482        !
483        no3budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
484        no3budget = no3budget / areatot
485        CALL iom_put( "pno3tot", no3budget )
486      ENDIF
487      !
488      ! Compute the budget of PO4
489      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
490         IF( ln_p4z ) THEN
491            zwork(:,:,:) =    trn(:,:,:,jppo4)                                         &
492               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
493               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
494               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
495        ELSE
496            zwork(:,:,:) =    trn(:,:,:,jppo4) + trn(:,:,:,jppph)                      &
497               &          +   trn(:,:,:,jppdi) + trn(:,:,:,jpppi)                      & 
498               &          +   trn(:,:,:,jppop) + trn(:,:,:,jpgop) + trn(:,:,:,jpdop)   &
499               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * po4rat3 
500        ENDIF
501        !
502        po4budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
503        po4budget = po4budget / areatot
504        CALL iom_put( "ppo4tot", po4budget )
505      ENDIF
506      !
507      ! Compute the budget of SiO3
508      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
509         zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
510         !
511         silbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
512         silbudget = silbudget / areatot
513         CALL iom_put( "psiltot", silbudget )
514      ENDIF
515      !
516      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
517         zwork(:,:,:) =  trn(:,:,:,jpno3) * rno3 + trn(:,:,:,jptal) + trn(:,:,:,jpcal) * 2.             
518         !
519         alkbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )         !
520         alkbudget = alkbudget / areatot
521         CALL iom_put( "palktot", alkbudget )
522      ENDIF
523      !
524      ! Compute the budget of Iron
525      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
526         zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
527            &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
528            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3   
529         !
530         ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  ) 
531         ferbudget = ferbudget / areatot
532         CALL iom_put( "pfertot", ferbudget )
533      ENDIF
534      !
535      ! Global budget of N SMS : denitrification in the water column and in the sediment
536      !                          nitrogen fixation by the diazotrophs
537      ! --------------------------------------------------------------------------------
538      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
539         znitrpottot  = glob_sum ( 'p4zsms', nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
540         CALL iom_put( "tnfix"  , znitrpottot * xfact3 )  ! Global  nitrogen fixation molC/l  to molN/m3
541      ENDIF
542      !
543      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
544         zrdenittot = glob_sum ( 'p4zsms', denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
545         zsdenittot = glob_sum ( 'p4zsms', sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )
546         CALL iom_put( "tdenit" , ( zrdenittot + zsdenittot ) * xfact3 )  ! Total denitrification molC/l to molN/m3
547      ENDIF
548      !
549      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
550         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( 'p4zsms', e1e2t(:,:) )
551         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
552         tpp            = tpp           * 1000. * xfact1
553         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
554         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
555         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
556             &                                no3budget * rno3 * 1.e+06, &
557             &                                po4budget * po4r * 1.e+06, &
558             &                                silbudget        * 1.e+06, &
559             &                                ferbudget        * 1.e+09
560         !
561         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
562            &                             zrdenittot  * xfact2  , &
563            &                             zsdenittot  * xfact2
564      ENDIF
565      !
566 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
567 9100  FORMAT(i8,5e18.10)
568 9200  FORMAT(i8,3f10.5)
569       !
570   END SUBROUTINE p4z_chk_mass
571
572   !!======================================================================
573END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.