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/trunk/src/TOP/PISCES/P4Z – NEMO

source: NEMO/trunk/src/TOP/PISCES/P4Z/p4zsms.F90 @ 14086

Last change on this file since 14086 was 14086, checked in by cetlod, 3 years ago

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