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 trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
File size: 28.8 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   !!   p4zsms         :  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 p4zsms.F90
33   PUBLIC   p4z_sms        ! called in p4zsms.F90
34
35   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget
36   REAL(wp) :: xfact1, xfact2, xfact3
37   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget
38
39   !!* Array used to indicate negative tracer values
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
41
42
43   !!----------------------------------------------------------------------
44   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
45   !! $Id: p4zsms.F90 3320 2012-03-05 16:37:52Z cetlod $
46   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48
49CONTAINS
50
51   SUBROUTINE p4z_sms( kt )
52      !!---------------------------------------------------------------------
53      !!                     ***  ROUTINE p4z_sms  ***
54      !!
55      !! ** Purpose :   Managment of the call to Biological sources and sinks
56      !!              routines of PISCES bio-model
57      !!
58      !! ** Method  : - at each new day ...
59      !!              - several calls of bio and sed ???
60      !!              - ...
61      !!---------------------------------------------------------------------
62      !
63      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
64      !!
65      INTEGER ::   ji, jj, jk, jnt, jn, jl
66      REAL(wp) ::  ztra
67      CHARACTER (len=25) :: charout
68      !!---------------------------------------------------------------------
69      !
70      IF( nn_timing == 1 )  CALL timing_start('p4z_sms')
71      !
72      IF( kt == nittrc000 ) THEN
73        !
74        ALLOCATE( xnegtr(jpi,jpj,jpk) )
75        !
76        CALL p4z_che                              ! initialize the chemical constants
77        !
78        IF( .NOT. ln_rsttr ) THEN  ;   CALL ahini_for_at(hi)   !  set PH at kt=nit000
79        ELSE                       ;   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!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
102            DO jk = 1, jpk
103               DO jj = 1, jpj
104                  DO ji = 1, jpi
105                     trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
106                  END DO
107               END DO
108            END DO
109         END DO
110      ENDIF
111      !
112      IF( ndayflxtr /= nday_year ) THEN      ! New days
113         !
114         ndayflxtr = nday_year
115
116         IF(lwp) write(numout,*)
117         IF(lwp) write(numout,*) ' New chemical constants and various rates for biogeochemistry at new day : ', nday_year
118         IF(lwp) write(numout,*) '~~~~~~'
119
120         CALL p4z_che              ! computation of chemical constants
121         CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
122         !
123      ENDIF
124
125      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients
126
127      DO jnt = 1, nrdttrc          ! Potential time splitting if requested
128         !
129         CALL p4z_bio( kt, jnt )   ! Biology
130         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
131         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
132         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
133         !
134!$OMP PARALLEL
135!$OMP DO schedule(static) private(jk, jj, ji)
136         DO jk = 1, jpk
137            DO jj = 1, jpj
138               DO ji = 1, jpi
139                  xnegtr(ji,jj,jk) = 1.e0
140               END DO
141            END DO
142         END DO
143         DO jn = jp_pcs0, jp_pcs1
144!$OMP DO schedule(static) private(jk, jj, ji, ztra)
145            DO jk = 1, jpk
146               DO jj = 1, jpj
147                  DO ji = 1, jpi
148                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
149                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
150                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
151                     ENDIF
152                 END DO
153               END DO
154            END DO
155         END DO
156         !                                ! where at least 1 tracer concentration becomes negative
157         !                                !
158         DO jn = jp_pcs0, jp_pcs1
159!$OMP DO schedule(static) private(jk, jj, ji)
160            DO jk = 1, jpk
161               DO jj = 1, jpj
162                  DO ji = 1, jpi
163                     trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) + xnegtr(ji,jj,jk) * tra(ji,jj,jk,jn)
164                  END DO
165               END DO
166            END DO
167         END DO
168        !
169         DO jn = jp_pcs0, jp_pcs1
170!$OMP DO schedule(static) private(jk, jj, ji)
171            DO jk = 1, jpk
172               DO jj = 1, jpj
173                  DO ji = 1, jpi
174                     tra(ji,jj,jk,jn) = 0._wp
175                  END DO
176               END DO
177            END DO
178         END DO
179!$OMP END PARALLEL
180         !
181         IF( ln_top_euler ) THEN
182            DO jn = jp_pcs0, jp_pcs1
183!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
184               DO jk = 1, jpk
185                  DO jj = 1, jpj
186                     DO ji = 1, jpi
187                        trn(ji,jj,jk,jn) = trb(ji,jj,jk,jn)
188                     END DO
189                  END DO
190               END DO
191            END DO
192         ENDIF
193      END DO
194
195      !
196      IF( l_trdtrc ) THEN
197         DO jn = jp_pcs0, jp_pcs1
198           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
199         END DO
200      END IF
201      !
202      IF( lk_sed ) THEN 
203         !
204         CALL sed_model( kt )     !  Main program of Sediment model
205         !
206         DO jn = jp_pcs0, jp_pcs1
207           CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
208         END DO
209         !
210      ENDIF
211      !
212      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file
213      !
214
215      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt ) ! Mass conservation checking
216
217      IF ( lwm .AND. kt == nittrc000 ) CALL FLUSH    ( numonp )     ! flush output namelist PISCES
218      IF( nn_timing == 1 )  CALL timing_stop('p4z_sms')
219      !
220      !
221   END SUBROUTINE p4z_sms
222
223   SUBROUTINE p4z_sms_init
224      !!----------------------------------------------------------------------
225      !!                     ***  p4z_sms_init  *** 
226      !!
227      !! ** Purpose :   read PISCES namelist
228      !!
229      !! ** input   :   file 'namelist.trc.s' containing the following
230      !!             namelist: natext, natbio, natsms
231      !!----------------------------------------------------------------------
232      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
233         &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet,  &
234         &                   no3rat3, po4rat3
235
236      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
237      NAMELIST/nampismass/ ln_check_mass
238      INTEGER :: ios                 ! Local integer output status for namelist read
239      !!----------------------------------------------------------------------
240
241      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
242      READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
243901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
244
245      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
246      READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
247902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
248      IF(lwm) WRITE ( numonp, nampisbio )
249
250      IF(lwp) THEN                         ! control print
251         WRITE(numout,*) ' Namelist : nampisbio'
252         WRITE(numout,*) '    frequence pour la biologie                nrdttrc    =', nrdttrc
253         WRITE(numout,*) '    POC sinking speed                         wsbio      =', wsbio
254         WRITE(numout,*) '    half saturation constant for mortality    xkmort     =', xkmort 
255         IF( ln_p5z ) THEN
256            WRITE(numout,*) '    N/C in zooplankton                        no3rat3    =', no3rat3
257            WRITE(numout,*) '    P/C in zooplankton                        po4rat3    =', po4rat3
258         ENDIF
259         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3     =', ferat3
260         WRITE(numout,*) '    Big particles sinking speed               wsbio2     =', wsbio2
261         WRITE(numout,*) '    Big particles maximum sinking speed       wsbio2max  =', wsbio2max
262         WRITE(numout,*) '    Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
263         WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
264         WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
265         IF( ln_ligand ) THEN
266            WRITE(numout,*) '    FeP sinking speed                             wfep   =', wfep
267            IF( ln_p4z ) THEN
268              WRITE(numout,*) '    Phyto ligand production per unit doc          ldocp  =', ldocp
269              WRITE(numout,*) '    Zoo ligand production per unit doc            ldocz  =', ldocz
270              WRITE(numout,*) '    Proportional loss of ligands due to Fe uptake lthet  =', lthet
271            ENDIF
272         ENDIF
273      ENDIF
274
275
276      REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
277      READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
278905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
279
280      REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
281      READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
282906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
283      IF(lwm) WRITE ( numonp, nampisdmp )
284
285      IF(lwp) THEN                         ! control print
286         WRITE(numout,*)
287         WRITE(numout,*) ' Namelist : nampisdmp'
288         WRITE(numout,*) '    Relaxation of tracer to glodap mean value             ln_pisdmp      =', ln_pisdmp
289         WRITE(numout,*) '    Frequency of Relaxation                               nn_pisdmp      =', nn_pisdmp
290         WRITE(numout,*) ' '
291      ENDIF
292
293      REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
294      READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
295907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
296
297      REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check
298      READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
299908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
300      IF(lwm) WRITE ( numonp, nampismass )
301
302      IF(lwp) THEN                         ! control print
303         WRITE(numout,*) ' '
304         WRITE(numout,*) ' Namelist parameter for mass conservation checking'
305         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
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   SUBROUTINE p4z_rst( kt, cdrw )
312      !!---------------------------------------------------------------------
313      !!                   ***  ROUTINE p4z_rst  ***
314      !!
315      !!  ** Purpose : Read or write variables in restart file:
316      !!
317      !!  WRITE(READ) mode:
318      !!       kt        : number of time step since the begining of the experiment at the
319      !!                   end of the current(previous) run
320      !!---------------------------------------------------------------------
321      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
322      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
323      !
324      !!---------------------------------------------------------------------
325
326      IF( TRIM(cdrw) == 'READ' ) THEN
327         !
328         IF(lwp) WRITE(numout,*)
329         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
330         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
331         !
332         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
333            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
334         ELSE
335            CALL ahini_for_at(hi)
336         ENDIF
337         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
338         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
339            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
340         ELSE
341            xksimax(:,:) = xksi(:,:)
342         ENDIF
343         !
344         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
345            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
346         ELSE
347            t_oce_co2_flx_cum = 0._wp
348         ENDIF
349         !
350         IF( ln_p5z ) THEN
351            IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
352               CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sized(:,:,:)  )
353               CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sized(:,:,:)  )
354               CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
355            ELSE
356               sizep(:,:,:) = 1.
357               sizen(:,:,:) = 1.
358               sized(:,:,:) = 1.
359            ENDIF
360        ENDIF
361        !
362      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
363         IF( kt == nitrst ) THEN
364            IF(lwp) WRITE(numout,*)
365            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
366            IF(lwp) WRITE(numout,*) '~~~~~~~'
367         ENDIF
368         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
369         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
370         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
371         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum )
372         IF( ln_p5z ) THEN
373            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sized(:,:,:) )
374            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sized(:,:,:) )
375            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )
376         ENDIF
377      ENDIF
378      !
379   END SUBROUTINE p4z_rst
380
381   SUBROUTINE p4z_dmp( kt )
382      !!----------------------------------------------------------------------
383      !!                    ***  p4z_dmp  ***
384      !!
385      !! ** purpose  : Relaxation of some tracers
386      !!----------------------------------------------------------------------
387      !
388      INTEGER, INTENT( in )  ::     kt ! time step
389      INTEGER ::   ji, jj, jk
390      !
391      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
392      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
393      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
394      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
395      !
396      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
397      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
398      REAL(wp), POINTER, DIMENSION(:,:,:) :: zctrn_jptal, zctrn_jppo4, zctrn_jppo3, zctrn_jpsil !workspace arrays
399      REAL(wp), POINTER, DIMENSION(:,:,:) :: zctrb_jptal, zctrb_jppo4, zctrb_jppo3, zctrb_jpsil !workspace arrays
400      !!---------------------------------------------------------------------
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" .AND. .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
408         !                                                ! --------------------------- !
409         CALL wrk_alloc( jpi, jpj, jpk, zctrn_jptal, zctrn_jppo4, zctrn_jppo3, zctrn_jpsil )
410         CALL wrk_alloc( jpi, jpj, jpk, zctrb_jptal, zctrb_jppo4, zctrb_jppo3, zctrb_jpsil )
411
412         ! set total alkalinity, phosphate, nitrate & silicate
413         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6             
414
415!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
416         DO jk = 1, jpk
417            DO jj = 1, jpj
418               DO ji = 1, jpi
419                  zctrn_jptal(ji,jj,jk) = trn(ji,jj,jk,jptal) * cvol(ji,jj,jk)
420                  zctrn_jppo4(ji,jj,jk) = trn(ji,jj,jk,jppo4) * cvol(ji,jj,jk)
421                  zctrn_jppo3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * cvol(ji,jj,jk)
422                  zctrn_jpsil(ji,jj,jk) = trn(ji,jj,jk,jpsil) * cvol(ji,jj,jk)
423               END DO
424            END DO
425         END DO
426
427         zalksumn = glob_sum( zctrn_jptal(:,:,:)  ) * zarea
428         zpo4sumn = glob_sum( zctrn_jppo4(:,:,:)  ) * zarea * po4r
429         zno3sumn = glob_sum( zctrn_jppo3(:,:,:)  ) * zarea * rno3
430         zsilsumn = glob_sum( zctrn_jpsil(:,:,:)  ) * zarea
431
432!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
433         DO jk = 1, jpk
434            DO jj = 1, jpj
435               DO ji = 1, jpi
436                  trn(ji,jj,jk,jpsil) = MIN( 400.e-6,trn(ji,jj,jk,jpsil) * silmean / zsilsumn )
437                  trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) * alkmean / zalksumn
438                  trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) * po4mean / zpo4sumn
439                  trn(ji,jj,jk,jpno3) = trn(ji,jj,jk,jpno3) * no3mean / zno3sumn
440               END DO
441            END DO
442         END DO
443
444         IF(lwp) THEN
445                WRITE(numout,*) '       TALKN mean : ', zalksumn
446                WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
447                WRITE(numout,*) '       NO3N  mean : ', zno3sumn
448                WRITE(numout,*) '       SiO3N mean : ', zsilsumn
449         END IF
450         !
451         !
452         IF( .NOT. ln_top_euler ) THEN
453!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
454            DO jk = 1, jpk
455               DO jj = 1, jpj
456                  DO ji = 1, jpi
457                     zctrb_jptal(ji,jj,jk) = trb(ji,jj,jk,jptal) * cvol(ji,jj,jk)
458                     zctrb_jppo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) * cvol(ji,jj,jk)
459                     zctrb_jppo3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * cvol(ji,jj,jk)
460                     zctrb_jpsil(ji,jj,jk) = trb(ji,jj,jk,jpsil) * cvol(ji,jj,jk)
461                  END DO
462               END DO
463            END DO
464
465            zalksumb = glob_sum( zctrb_jptal(:,:,:)  ) * zarea
466            zpo4sumb = glob_sum( zctrb_jppo4(:,:,:)  ) * zarea * po4r
467            zno3sumb = glob_sum( zctrb_jppo3(:,:,:)  ) * zarea * rno3
468            zsilsumb = glob_sum( zctrb_jpsil(:,:,:)  ) * zarea
469
470!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
471            DO jk = 1, jpk
472               DO jj = 1, jpj
473                  DO ji = 1, jpi
474                     trb(ji,jj,jk,jpsil) = MIN( 400.e-6,trb(ji,jj,jk,jpsil) * silmean / zsilsumb )
475                     trb(ji,jj,jk,jptal) = trb(ji,jj,jk,jptal) * alkmean / zalksumb
476                     trb(ji,jj,jk,jppo4) = trb(ji,jj,jk,jppo4) * po4mean / zpo4sumb
477                     trb(ji,jj,jk,jpno3) = trb(ji,jj,jk,jpno3) * no3mean / zno3sumb
478                  END DO
479               END DO
480            END DO
481
482            IF(lwp) THEN
483                WRITE(numout,*) ' '
484                WRITE(numout,*) '       TALKB mean : ', zalksumb
485                WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
486                WRITE(numout,*) '       NO3B  mean : ', zno3sumb
487                WRITE(numout,*) '       SiO3B mean : ', zsilsumb
488            END IF
489        ENDIF
490        !
491        CALL wrk_dealloc( jpi, jpj, jpk, zctrb_jptal, zctrb_jppo4, zctrb_jppo3, zctrb_jpsil )
492        CALL wrk_dealloc( jpi, jpj, jpk, zctrn_jptal, zctrn_jppo4, zctrn_jppo3, zctrn_jpsil )
493        !
494      ENDIF
495        !
496   END SUBROUTINE p4z_dmp
497
498
499   SUBROUTINE p4z_chk_mass( kt )
500      !!----------------------------------------------------------------------
501      !!                  ***  ROUTINE p4z_chk_mass  ***
502      !!
503      !! ** Purpose :  Mass conservation check
504      !!
505      !!---------------------------------------------------------------------
506      !
507      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
508      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
509      CHARACTER(LEN=100)   ::   cltxt
510      INTEGER :: jk
511      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwork
512      !!----------------------------------------------------------------------
513
514      !
515      !!---------------------------------------------------------------------
516
517      IF( kt == nittrc000 ) THEN
518         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
519            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
520            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
521            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
522            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
523            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
524            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
525            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
526            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
527            IF( lwp ) WRITE(numnut,*) 
528         ENDIF
529      ENDIF
530
531      CALL wrk_alloc( jpi, jpj, jpk, zwork )
532      !
533      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
534         !   Compute the budget of NO3, ALK, Si, Fer
535         IF( ln_p4z ) THEN
536            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
537               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
538               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
539               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
540        ELSE
541            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4) + trn(:,:,:,jpnph)   &
542               &          +   trn(:,:,:,jpndi) + trn(:,:,:,jpnpi)                      & 
543               &          +   trn(:,:,:,jppon) + trn(:,:,:,jpgon) + trn(:,:,:,jpdon)   &
544               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * no3rat3 
545        ENDIF
546        !
547        no3budget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
548        no3budget = no3budget / areatot
549        CALL iom_put( "pno3tot", no3budget )
550      ENDIF
551      !
552      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
553         IF( ln_p4z ) THEN
554            zwork(:,:,:) =    trn(:,:,:,jppo4)                                         &
555               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
556               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &       
557               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
558        ELSE
559            zwork(:,:,:) =    trn(:,:,:,jppo4) + trn(:,:,:,jppph)                      &
560               &          +   trn(:,:,:,jppdi) + trn(:,:,:,jpppi)                      & 
561               &          +   trn(:,:,:,jppop) + trn(:,:,:,jpgop) + trn(:,:,:,jpdop)   &
562               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * po4rat3 
563        ENDIF
564        !
565        po4budget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
566        po4budget = po4budget / areatot
567        CALL iom_put( "ppo4tot", po4budget )
568      ENDIF
569      !
570      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
571         zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
572         !
573         silbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
574         silbudget = silbudget / areatot
575         CALL iom_put( "psiltot", silbudget )
576      ENDIF
577      !
578      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
579         zwork(:,:,:) =  trn(:,:,:,jpno3) * rno3 + trn(:,:,:,jptal) + trn(:,:,:,jpcal) * 2.             
580         !
581         alkbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  )         !
582         alkbudget = alkbudget / areatot
583         CALL iom_put( "palktot", alkbudget )
584      ENDIF
585      !
586      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
587         zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
588            &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
589            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3   
590         IF( ln_ligand)  zwork(:,:,:) = zwork(:,:,:) + trn(:,:,:,jpfep)               
591         !
592         ferbudget = glob_sum( zwork(:,:,:) * cvol(:,:,:)  ) 
593         ferbudget = ferbudget / areatot
594         CALL iom_put( "pfertot", ferbudget )
595      ENDIF
596      !
597      CALL wrk_dealloc( jpi, jpj, jpk, zwork )
598      !
599      ! Global budget of N SMS : denitrification in the water column and in the sediment
600      !                          nitrogen fixation by the diazotrophs
601      ! --------------------------------------------------------------------------------
602      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
603         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
604         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3
605      ENDIF
606      !
607      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
608         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
609         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3
610      ENDIF
611      !
612      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
613         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) )
614         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments
615      ENDIF
616
617      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
618         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
619         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
620         tpp            = tpp           * 1000. * xfact1
621         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
622         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
623         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
624             &                                no3budget * rno3 * 1.e+06, &
625             &                                po4budget * po4r * 1.e+06, &
626             &                                silbudget        * 1.e+06, &
627             &                                ferbudget        * 1.e+09
628         !
629         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
630         &                             zrdenittot  * xfact2  , &
631         &                             zsdenittot  * xfact2
632
633      ENDIF
634      !
635 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
636 9100  FORMAT(i8,5e18.10)
637 9200  FORMAT(i8,3f10.5)
638
639       !
640   END SUBROUTINE p4z_chk_mass
641
642   !!======================================================================
643END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.