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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 @ 9559

Last change on this file since 9559 was 9559, checked in by cetlod, 6 years ago

Minor corrections to compute PISCES chemical constant at every time-step rather than every day

File size: 24.9 KB
RevLine 
[3443]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   !!----------------------------------------------------------------------
[9169]9   !!   p4z_sms        : Time loop of passive tracers sms
[3443]10   !!----------------------------------------------------------------------
[9169]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
[3443]28
29   IMPLICIT NONE
30   PRIVATE
31
[4990]32   PUBLIC   p4z_sms_init   ! called in p4zsms.F90
33   PUBLIC   p4z_sms        ! called in p4zsms.F90
[3443]34
[9169]35   INTEGER ::    numco2, numnut, numnit      ! logical unit for co2 budget
36   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget
37   REAL(wp) ::   xfact1, xfact2, xfact3
[3443]38
[9169]39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     ! Array used to indicate negative tracer values
[5385]40
[3443]41   !!----------------------------------------------------------------------
42   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
43   !! $Id: p4zsms.F90 3320 2012-03-05 16:37:52Z cetlod $
44   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
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  : - at each new day ...
56      !!              - several calls of bio and sed ???
57      !!              - ...
58      !!---------------------------------------------------------------------
59      !
60      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
61      !!
[5385]62      INTEGER ::   ji, jj, jk, jnt, jn, jl
63      REAL(wp) ::  ztra
[3443]64      CHARACTER (len=25) :: charout
65      !!---------------------------------------------------------------------
66      !
[9124]67      IF( ln_timing )   CALL timing_start('p4z_sms')
[3443]68      !
[4152]69      IF( kt == nittrc000 ) THEN
70        !
[5385]71        ALLOCATE( xnegtr(jpi,jpj,jpk) )
72        !
[9559]73        IF( .NOT. ln_rsttr ) THEN
74            CALL p4z_che                              ! initialize the chemical constants
75            CALL ahini_for_at(hi)   !  set PH at kt=nit000
76        ELSE
77            CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
[4152]78        ENDIF
79        !
80      ENDIF
81      !
[3443]82      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
[3496]83      !
[7646]84      rfact = r2dttrc
[5385]85      !
86      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
87         rfactr  = 1. / rfact
[7646]88         rfact2  = rfact / REAL( nrdttrc, wp )
[5385]89         rfact2r = 1. / rfact2
90         xstep = rfact2 / rday         ! Time step duration for biology
91         IF(lwp) WRITE(numout,*) 
[6140]92         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
[5385]93         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
94         IF(lwp) WRITE(numout,*)
95      ENDIF
96
97      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
98         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
[7753]99            trb(:,:,:,jn) = trn(:,:,:,jn)
[5385]100         END DO
101      ENDIF
102      !
[9559]103      CALL p4z_che              ! computation of chemical constants
104      CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
105      !
[3443]106      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients
107
108      DO jnt = 1, nrdttrc          ! Potential time splitting if requested
109         !
[5385]110         CALL p4z_bio( kt, jnt )   ! Biology
111         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
[6421]112         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
[5385]113         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
[3443]114         !
[7753]115         xnegtr(:,:,:) = 1.e0
[3443]116         DO jn = jp_pcs0, jp_pcs1
[5385]117            DO jk = 1, jpk
118               DO jj = 1, jpj
119                  DO ji = 1, jpi
120                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
121                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
122                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
123                     ENDIF
124                 END DO
125               END DO
126            END DO
127         END DO
128         !                                ! where at least 1 tracer concentration becomes negative
129         !                                !
130         DO jn = jp_pcs0, jp_pcs1
[7753]131           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
[5385]132         END DO
133        !
134         DO jn = jp_pcs0, jp_pcs1
[7753]135            tra(:,:,:,jn) = 0._wp
[5385]136         END DO
[3443]137         !
[5385]138         IF( ln_top_euler ) THEN
139            DO jn = jp_pcs0, jp_pcs1
[7753]140               trn(:,:,:,jn) = trb(:,:,:,jn)
[5385]141            END DO
142         ENDIF
[3443]143      END DO
144
145      !
[5385]146      IF( l_trdtrc ) THEN
147         DO jn = jp_pcs0, jp_pcs1
148           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
149         END DO
150      END IF
151      !
[3443]152      IF( lk_sed ) THEN 
153         !
154         CALL sed_model( kt )     !  Main program of Sediment model
155         !
156      ENDIF
157      !
158      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file
159      !
[5385]160
[9124]161      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt )    ! Mass conservation checking
[3481]162
[9124]163      IF( lwm .AND. kt == nittrc000    )  CALL FLUSH( numonp )       ! flush output namelist PISCES
[3443]164      !
[9124]165      IF( ln_timing )  CALL timing_stop('p4z_sms')
[3481]166      !
[3443]167   END SUBROUTINE p4z_sms
168
[9124]169
[3443]170   SUBROUTINE p4z_sms_init
171      !!----------------------------------------------------------------------
172      !!                     ***  p4z_sms_init  *** 
173      !!
174      !! ** Purpose :   read PISCES namelist
175      !!
176      !! ** input   :   file 'namelist.trc.s' containing the following
177      !!             namelist: natext, natbio, natsms
178      !!----------------------------------------------------------------------
[9124]179      INTEGER :: ios                 ! Local integer output status for namelist read
180      !!
[7646]181      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
182         &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet,  &
183         &                   no3rat3, po4rat3
[9124]184         !
[4148]185      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
[3531]186      NAMELIST/nampismass/ ln_check_mass
187      !!----------------------------------------------------------------------
[9169]188      !
189      IF(lwp) THEN
190         WRITE(numout,*)
191         WRITE(numout,*) 'p4z_sms_init : PISCES initialization'
192         WRITE(numout,*) '~~~~~~~~~~~~'
193      ENDIF
[3443]194
[4147]195      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
196      READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
[9169]197901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
[4147]198      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
199      READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
[9169]200902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
201      IF(lwm) WRITE( numonp, nampisbio )
202      !
[3443]203      IF(lwp) THEN                         ! control print
[9169]204         WRITE(numout,*) '   Namelist : nampisbio'
205         WRITE(numout,*) '      frequency for the biology                 nrdttrc     =', nrdttrc
206         WRITE(numout,*) '      POC sinking speed                         wsbio       =', wsbio
207         WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort 
[7646]208         IF( ln_p5z ) THEN
[9169]209            WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3
210            WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3
[7646]211         ENDIF
[9169]212         WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3
213         WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2
214         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
215         WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
216         WRITE(numout,*) '      Maximum number of iterations for POC      niter1max   =', niter1max
217         WRITE(numout,*) '      Maximum number of iterations for GOC      niter2max   =', niter2max
[7646]218         IF( ln_ligand ) THEN
[9169]219            WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep
[7646]220            IF( ln_p4z ) THEN
[9169]221               WRITE(numout,*) '      Phyto ligand production per unit doc           ldocp  =', ldocp
222               WRITE(numout,*) '      Zoo ligand production per unit doc             ldocz  =', ldocz
223               WRITE(numout,*) '      Proportional loss of ligands due to Fe uptake  lthet  =', lthet
[7646]224            ENDIF
225         ENDIF
[3443]226      ENDIF
227
228
[4147]229      REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
230      READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
[9169]231905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
[4147]232      REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
233      READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
[9169]234906   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
235      IF(lwm) WRITE( numonp, nampisdmp )
236      !
[3443]237      IF(lwp) THEN                         ! control print
238         WRITE(numout,*)
[9169]239         WRITE(numout,*) '   Namelist : nampisdmp --- relaxation to GLODAP'
240         WRITE(numout,*) '      Relaxation of tracer to glodap mean value   ln_pisdmp =', ln_pisdmp
241         WRITE(numout,*) '      Frequency of Relaxation                     nn_pisdmp =', nn_pisdmp
[3443]242      ENDIF
243
[4147]244      REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
245      READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
[9169]246907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
[4147]247      REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check
248      READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
[9169]249908   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
250      IF(lwm) WRITE( numonp, nampismass )
[4147]251
[3531]252      IF(lwp) THEN                         ! control print
[9169]253         WRITE(numout,*)
254         WRITE(numout,*) '   Namelist : nampismass  --- mass conservation checking'
255         WRITE(numout,*) '      Flag to check mass conservation of NO3/Si/TALK   ln_check_mass = ', ln_check_mass
[3531]256      ENDIF
[9124]257      !
[3443]258   END SUBROUTINE p4z_sms_init
259
[9124]260
[3443]261   SUBROUTINE p4z_rst( kt, cdrw )
262      !!---------------------------------------------------------------------
263      !!                   ***  ROUTINE p4z_rst  ***
264      !!
265      !!  ** Purpose : Read or write variables in restart file:
266      !!
267      !!  WRITE(READ) mode:
268      !!       kt        : number of time step since the begining of the experiment at the
269      !!                   end of the current(previous) run
270      !!---------------------------------------------------------------------
271      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
272      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
[9124]273      !!---------------------------------------------------------------------
[3443]274      !
275      IF( TRIM(cdrw) == 'READ' ) THEN
276         !
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
279         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
280         !
281         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
282            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
283         ELSE
[7646]284            CALL ahini_for_at(hi)
[3443]285         ENDIF
286         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
287         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
288            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
289         ELSE
290            xksimax(:,:) = xksi(:,:)
291         ENDIF
292         !
[4996]293         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
294            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
295         ELSE
296            t_oce_co2_flx_cum = 0._wp
297         ENDIF
298         !
[7646]299         IF( ln_p5z ) THEN
300            IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
301               CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sized(:,:,:)  )
302               CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sized(:,:,:)  )
303               CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
304            ELSE
305               sizep(:,:,:) = 1.
306               sizen(:,:,:) = 1.
307               sized(:,:,:) = 1.
308            ENDIF
309        ENDIF
310        !
[3443]311      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
312         IF( kt == nitrst ) THEN
313            IF(lwp) WRITE(numout,*)
314            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
315            IF(lwp) WRITE(numout,*) '~~~~~~~'
316         ENDIF
317         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
318         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
319         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
[4996]320         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum )
[7646]321         IF( ln_p5z ) THEN
322            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sized(:,:,:) )
323            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sized(:,:,:) )
324            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )
325         ENDIF
[3443]326      ENDIF
327      !
328   END SUBROUTINE p4z_rst
329
[9124]330
[3443]331   SUBROUTINE p4z_dmp( kt )
332      !!----------------------------------------------------------------------
333      !!                    ***  p4z_dmp  ***
334      !!
335      !! ** purpose  : Relaxation of some tracers
336      !!----------------------------------------------------------------------
337      !
338      INTEGER, INTENT( in )  ::     kt ! time step
339      !
340      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
341      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
342      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
343      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
344      !
[5385]345      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
346      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
[3443]347      !!---------------------------------------------------------------------
348
349      IF(lwp)  WRITE(numout,*)
[3557]350      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
[3443]351      IF(lwp)  WRITE(numout,*)
352
[7646]353      IF( cn_cfg == "orca" .AND. .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
354         !                                                ! --------------------------- !
[3443]355         ! set total alkalinity, phosphate, nitrate & silicate
356         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6             
357
[7753]358         zalksumn = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
359         zpo4sumn = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
360         zno3sumn = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
361         zsilsumn = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
362 
363         IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
364         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
[3443]365
[7753]366         IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
367         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
[3443]368
[7753]369         IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
370         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
[3443]371
[7753]372         IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
373         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
[3443]374         !
[5385]375         !
376         IF( .NOT. ln_top_euler ) THEN
[7753]377            zalksumb = glob_sum( trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
378            zpo4sumb = glob_sum( trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
379            zno3sumb = glob_sum( trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
380            zsilsumb = glob_sum( trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
381 
382            IF(lwp) WRITE(numout,*) ' '
383            IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
384            trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
[5385]385
[7753]386            IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
387            trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
[5385]388
[7753]389            IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
390            trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
[5385]391
[7753]392            IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
393            trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
[5385]394        ENDIF
395        !
[3443]396      ENDIF
[5385]397        !
[3443]398   END SUBROUTINE p4z_dmp
399
400
401   SUBROUTINE p4z_chk_mass( kt )
402      !!----------------------------------------------------------------------
403      !!                  ***  ROUTINE p4z_chk_mass  ***
404      !!
405      !! ** Purpose :  Mass conservation check
406      !!
407      !!---------------------------------------------------------------------
[5547]408      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
409      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
[5385]410      CHARACTER(LEN=100)   ::   cltxt
411      INTEGER :: jk
[9125]412      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork
[5385]413      !!----------------------------------------------------------------------
414      !
[3443]415      IF( kt == nittrc000 ) THEN
[8533]416         xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
417         xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
418         xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
[3451]419         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
[3496]420            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
421            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
[5385]422            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
423            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
424            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
425            IF( lwp ) WRITE(numnut,*) 
[3443]426         ENDIF
427      ENDIF
428
[4996]429      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
430         !   Compute the budget of NO3, ALK, Si, Fer
[7646]431         IF( ln_p4z ) THEN
432            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
433               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
434               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
435               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
436        ELSE
437            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4) + trn(:,:,:,jpnph)   &
438               &          +   trn(:,:,:,jpndi) + trn(:,:,:,jpnpi)                      & 
439               &          +   trn(:,:,:,jppon) + trn(:,:,:,jpgon) + trn(:,:,:,jpdon)   &
440               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * no3rat3 
441        ENDIF
442        !
443        no3budget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
444        no3budget = no3budget / areatot
445        CALL iom_put( "pno3tot", no3budget )
[4996]446      ENDIF
447      !
[5385]448      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
[7646]449         IF( ln_p4z ) THEN
450            zwork(:,:,:) =    trn(:,:,:,jppo4)                                         &
451               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
452               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
453               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
454        ELSE
455            zwork(:,:,:) =    trn(:,:,:,jppo4) + trn(:,:,:,jppph)                      &
456               &          +   trn(:,:,:,jppdi) + trn(:,:,:,jpppi)                      & 
457               &          +   trn(:,:,:,jppop) + trn(:,:,:,jpgop) + trn(:,:,:,jpdop)   &
458               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * po4rat3 
459        ENDIF
460        !
461        po4budget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
462        po4budget = po4budget / areatot
463        CALL iom_put( "ppo4tot", po4budget )
[5385]464      ENDIF
465      !
466      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
[7646]467         zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
[4996]468         !
[7646]469         silbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
[4996]470         silbudget = silbudget / areatot
471         CALL iom_put( "psiltot", silbudget )
472      ENDIF
473      !
[5385]474      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
[7646]475         zwork(:,:,:) =  trn(:,:,:,jpno3) * rno3 + trn(:,:,:,jptal) + trn(:,:,:,jpcal) * 2.             
[4996]476         !
[7646]477         alkbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  )         !
[4996]478         alkbudget = alkbudget / areatot
479         CALL iom_put( "palktot", alkbudget )
480      ENDIF
481      !
[5385]482      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
[7646]483         zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
484            &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
485            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3   
486         IF( ln_ligand)  zwork(:,:,:) = zwork(:,:,:) + trn(:,:,:,jpfep)               
[3496]487         !
[7646]488         ferbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
[3496]489         ferbudget = ferbudget / areatot
[4996]490         CALL iom_put( "pfertot", ferbudget )
[3496]491      ENDIF
[4996]492      !
[5385]493      ! Global budget of N SMS : denitrification in the water column and in the sediment
494      !                          nitrogen fixation by the diazotrophs
495      ! --------------------------------------------------------------------------------
496      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
497         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
[8533]498         CALL iom_put( "tnfix"  , znitrpottot * xfact3 )  ! Global  nitrogen fixation molC/l  to molN/m3
[5385]499      ENDIF
500      !
501      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
[8533]502         zrdenittot = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
503         zsdenittot = glob_sum ( sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )
504         CALL iom_put( "tdenit" , ( zrdenittot + zsdenittot ) * xfact3 )  ! Total denitrification molC/l to molN/m3
[5385]505      ENDIF
506      !
[4996]507      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
508         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
[5385]509         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
510         tpp            = tpp           * 1000. * xfact1
511         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
[4996]512         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
[5385]513         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
[4996]514             &                                no3budget * rno3 * 1.e+06, &
[5385]515             &                                po4budget * po4r * 1.e+06, &
[4996]516             &                                silbudget        * 1.e+06, &
517             &                                ferbudget        * 1.e+09
[5385]518         !
519         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
[9124]520            &                             zrdenittot  * xfact2  , &
521            &                             zsdenittot  * xfact2
[4996]522      ENDIF
523      !
[3496]524 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
[5385]525 9100  FORMAT(i8,5e18.10)
526 9200  FORMAT(i8,3f10.5)
[3443]527       !
528   END SUBROUTINE p4z_chk_mass
529
530   !!======================================================================
531END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.