source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 9309

Last change on this file since 9309 was 9262, checked in by frrh, 3 years ago

Add code from Julien's Met Office GMED ticket 364 to output
details relating to tracer conservation checks.

Command:
svn merge -r 9204:9260
svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/NERC/dev_r5518_GO6_conserv_check_up

File size: 41.0 KB
Line 
1MODULE trcrst
2   !!======================================================================
3   !!                         ***  MODULE trcrst  ***
4   !! TOP :   Manage the passive tracer restart
5   !!======================================================================
6   !! History :    -   !  1991-03  ()  original code
7   !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90
8   !!              -   !  2005-10 (C. Ethe) print control
9   !!             2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture
10   !!----------------------------------------------------------------------
11#if defined key_top
12   !!----------------------------------------------------------------------
13   !!   'key_top'                                                TOP models
14   !!----------------------------------------------------------------------
15   !!----------------------------------------------------------------------
16   !!   trc_rst :   Restart for passive tracer
17   !!----------------------------------------------------------------------
18   !!----------------------------------------------------------------------
19   !!   'key_top'                                                TOP models
20   !!----------------------------------------------------------------------
21   !!   trc_rst_opn    : open  restart file
22   !!   trc_rst_read   : read  restart file
23   !!   trc_rst_wri    : write restart file
24   !!----------------------------------------------------------------------
25   USE oce_trc
26   USE trc
27   USE trcnam_trp
28   USE iom
29   USE ioipsl, ONLY : ju2ymds    ! for calendar
30   USE daymod
31   !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
32   USE par_medusa
33   USE sms_medusa
34   USE trcsms_medusa
35   !!
36#if defined key_idtra
37   USE trcsms_idtra
38#endif
39   !!
40#if defined key_cfc
41   USE trcsms_cfc
42#endif
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE sbc_oce, ONLY: lk_oasis 
45   USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
46   USE trcstat
47
48   IMPLICIT NONE
49   PRIVATE
50
51   PUBLIC   trc_rst_opn       ! called by ???
52   PUBLIC   trc_rst_read      ! called by ???
53   PUBLIC   trc_rst_wri       ! called by ???
54   PUBLIC   trc_rst_cal
55   PUBLIC   trc_rst_stat
56#if defined key_medusa && defined key_roam
57   PUBLIC   trc_rst_conserve
58#endif
59
60   !! * Substitutions
61#  include "top_substitute.h90"
62   
63CONTAINS
64   
65   SUBROUTINE trc_rst_opn( kt )
66      !!----------------------------------------------------------------------
67      !!                    ***  trc_rst_opn  ***
68      !!
69      !! ** purpose  :   output of sea-trc variable in a netcdf file
70      !!----------------------------------------------------------------------
71      INTEGER, INTENT(in) ::   kt       ! number of iteration
72      INTEGER             ::   iyear, imonth, iday
73      REAL (wp)           ::   zsec
74      REAL (wp)           ::   zfjulday
75      !
76      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
77      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
78      CHARACTER(LEN=256)  ::   clpath   ! full path to ocean output restart file
79      !!----------------------------------------------------------------------
80      !
81      IF( lk_offline ) THEN
82         IF( kt == nittrc000 ) THEN
83            lrst_trc = .FALSE.
84            IF( ln_rst_list ) THEN
85               nrst_lst = 1
86               nitrst = nstocklist( nrst_lst )
87            ELSE
88               nitrst = nitend
89            ENDIF
90         ENDIF
91
92         IF( .NOT. ln_rst_list .AND. MOD( kt - 1, nstock ) == 0 ) THEN
93            ! we use kt - 1 and not kt - nittrc000 to keep the same periodicity from the beginning of the experiment
94            nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
95            IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
96         ENDIF
97      ELSE
98         IF( kt == nittrc000 ) lrst_trc = .FALSE.
99      ENDIF
100
101      ! to get better performances with NetCDF format:
102      ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1)
103      ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1
104      IF( kt == nitrst - 2*nn_dttrc .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
105         IF ( ln_rstdate ) THEN
106            !! JPALM -- 22-12-2015 -- modif to get the good date on restart trc file name
107            !!                     -- the condition to open the rst file is not the same than for the dynamic rst.
108            !!                     -- here it - for an obscure reason - is open 2 time-step before the restart writing process
109            !!                     instead of 1.
110            !!                     -- i am not sure if someone forgot +1 in the if loop condition as
111            !!                     it is writen in all comments nitrst -2*nn_dttrc + 1 and the condition is
112            !!                     nitrst - 2*nn_dttrc
113            !!                     -- nevertheless we didn't wanted to broke something already working
114            !!                     and just adapted the part we added.
115            !!                     -- So instead of calling ju2ymds( fjulday + (rdttra(1))
116            !!                     we call ju2ymds( fjulday + (2*rdttra(1))
117            !!--------------------------------------------------------------------     
118            zfjulday = fjulday + (2*rdttra(1)) / rday
119            IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
120            CALL ju2ymds( zfjulday + (2*rdttra(1)) / rday, iyear, imonth, iday, zsec )
121            WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
122         ELSE
123            ! beware of the format used to write kt (default is i8.8, that should be large enough)
124            IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
125            ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
126            ENDIF
127         ENDIF
128         ! create the file
129         IF(lwp) WRITE(numout,*)
130         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trcrst_out)
131         clpath = TRIM(cn_trcrst_outdir)
132         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
133         IF(lwp) WRITE(numout,*) &
134             '             open trc restart.output NetCDF file: ',TRIM(clpath)//clname
135         CALL iom_open( TRIM(clpath)//TRIM(clname), numrtw, ldwrt = .TRUE., kiolib = jprstlib )
136         lrst_trc = .TRUE.
137      ENDIF
138      !
139   END SUBROUTINE trc_rst_opn
140
141   SUBROUTINE trc_rst_read
142      !!----------------------------------------------------------------------
143      !!                    ***  trc_rst_opn  ***
144      !!
145      !! ** purpose  :   read passive tracer fields in restart files
146      !!----------------------------------------------------------------------
147      INTEGER  ::  jn, jl     
148      !! AXY (05/11/13): temporary variables
149      REAL(wp) ::    fq0,fq1,fq2
150
151      !!----------------------------------------------------------------------
152      !
153      IF(lwp) WRITE(numout,*)
154      IF(lwp) WRITE(numout,*) 'trc_rst_read : read data in the TOP restart file'
155      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
156
157      ! READ prognostic variables and computes diagnostic variable
158      DO jn = 1, jptra
159         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
160         trn(:,:,:,jn) = trn(:,:,:,jn) * tmask(:,:,:)
161      END DO
162
163      DO jn = 1, jptra
164         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
165         trb(:,:,:,jn) = trb(:,:,:,jn) * tmask(:,:,:)
166      END DO
167      !
168      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
169      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
170      !!                 version of NEMO date significantly earlier than the current
171      !!                 version
172
173#if defined key_medusa
174      !! AXY (13/01/12): check if the restart contains sediment fields;
175      !!                 this is only relevant for simulations that include
176      !!                 biogeochemistry and are restarted from earlier runs
177      !!                 in which there was no sediment component
178      !!
179      IF( iom_varid( numrtr, 'B_SED_N', ldstop = .FALSE. ) > 0 ) THEN
180         !! YES; in which case read them
181         !!
182         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields present - reading in ...'
183         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_N',  zb_sed_n(:,:)  )
184         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_N',  zn_sed_n(:,:)  )
185         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_FE', zb_sed_fe(:,:) )
186         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_FE', zn_sed_fe(:,:) )
187         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_SI', zb_sed_si(:,:) )
188         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_SI', zn_sed_si(:,:) )
189         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_C',  zb_sed_c(:,:)  )
190         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_C',  zn_sed_c(:,:)  )
191         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_CA', zb_sed_ca(:,:) )
192         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_CA', zn_sed_ca(:,:) )
193      ELSE
194         !! NO; in which case set them to zero
195         !!
196         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields absent - setting to zero ...'
197         zb_sed_n(:,:)  = 0.0   !! organic N
198         zn_sed_n(:,:)  = 0.0
199         zb_sed_fe(:,:) = 0.0   !! organic Fe
200         zn_sed_fe(:,:) = 0.0
201         zb_sed_si(:,:) = 0.0   !! inorganic Si
202         zn_sed_si(:,:) = 0.0
203         zb_sed_c(:,:)  = 0.0   !! organic C
204         zn_sed_c(:,:)  = 0.0
205         zb_sed_ca(:,:) = 0.0   !! inorganic C
206         zn_sed_ca(:,:) = 0.0
207      ENDIF
208      !!
209      !! calculate stats on these fields
210      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
211      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
212      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
213      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
214      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
215      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
216      !!
217      !! AXY (07/07/15): read in temporally averaged fields for DMS
218      !!                 calculations
219      !!
220      IF( iom_varid( numrtr, 'B_DMS_CHN', ldstop = .FALSE. ) > 0 ) THEN
221         !! YES; in which case read them
222         !!
223         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS present - reading in ...'
224         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
225         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
226         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
227         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
228         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
229         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
230         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
231         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
232         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_DIN',  zb_dms_din(:,:)  )
233         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_DIN',  zn_dms_din(:,:)  )
234      ELSE
235         !! NO; in which case set them to zero
236         !!
237         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS absent - setting to zero ...'
238         zb_dms_chn(:,:)  = 0.0   !! CHN
239         zn_dms_chn(:,:)  = 0.0
240         zb_dms_chd(:,:)  = 0.0   !! CHD
241         zn_dms_chd(:,:)  = 0.0
242         zb_dms_mld(:,:)  = 0.0   !! MLD
243         zn_dms_mld(:,:)  = 0.0
244         zb_dms_qsr(:,:)  = 0.0   !! QSR
245         zn_dms_qsr(:,:)  = 0.0
246         zb_dms_din(:,:)  = 0.0   !! DIN
247         zn_dms_din(:,:)  = 0.0
248      ENDIF
249      !! 
250      !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
251      !!                  -- needed for the coupling with atm
252      IF( iom_varid( numrtr, 'N_DMS_srf', ldstop = .FALSE. ) > 0 ) THEN
253         IF(lwp) WRITE(numout,*) 'DMS surf concentration - reading in ...'
254         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_srf',  zb_dms_srf(:,:)  )
255         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_srf',  zn_dms_srf(:,:)  )
256      ELSE
257         IF(lwp) WRITE(numout,*) 'DMS surf concentration - setting to zero ...'
258         zb_dms_srf(:,:)  = 0.0   !! DMS
259         zn_dms_srf(:,:)  = 0.0
260      ENDIF
261      IF (lk_oasis) THEN
262         DMS_out_cpl(:,:) = zn_dms_srf(:,:)        !! Coupling variable
263      END IF
264      !!
265      IF( iom_varid( numrtr, 'B_CO2_flx', ldstop = .FALSE. ) > 0 ) THEN
266         IF(lwp) WRITE(numout,*) 'CO2 air-sea flux - reading in ...'
267         CALL iom_get( numrtr, jpdom_autoglo, 'B_CO2_flx',  zb_co2_flx(:,:)  )
268         CALL iom_get( numrtr, jpdom_autoglo, 'N_CO2_flx',  zn_co2_flx(:,:)  )
269      ELSE
270         IF(lwp) WRITE(numout,*) 'CO2 air-sea flux - setting to zero ...'
271         zb_co2_flx(:,:)  = 0.0   !! CO2 flx
272         zn_co2_flx(:,:)  = 0.0
273      ENDIF
274      IF (lk_oasis) THEN
275         CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable
276      END IF
277      !!
278      !! JPALM 02-06-2017 -- in complement to DMS surf
279      !!                  -- the atm model needs surf Chl
280      !!                     as proxy of org matter from the ocean
281      !!                  -- needed for the coupling with atm
282      !!       07-12-2017 -- To make things cleaner, we want to store an 
283      !!                     unscaled Chl field in the restart and only
284      !!                     scale it when reading it in.
285
286      IF( iom_varid( numrtr, 'N_CHL_srf', ldstop = .FALSE. ) > 0 ) THEN
287         IF(lwp) WRITE(numout,*) 'Chl cpl concentration - reading in ... - scale by ', scl_chl
288         CALL iom_get( numrtr, jpdom_autoglo, 'N_CHL_srf',  zn_chl_srf(:,:)  )
289      ELSE
290         IF(lwp) WRITE(numout,*) 'set Chl coupled concentration - scaled by ', scl_chl
291         zn_chl_srf(:,:)  = MAX( 0.0, (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6 )
292      ENDIF
293      IF (lk_oasis) THEN
294         chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling variable
295      END IF
296      !!
297      !! calculate stats on these fields
298      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
299      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
300      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
301      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
302      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
303      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
304      call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
305      call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
306      IF (lk_oasis) THEN
307         call trc_rst_dia_stat(chloro_out_cpl(:,:), 'CHL  cpl')
308      END IF
309      !! 
310      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
311      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
312# if defined key_roam
313      IF( iom_varid( numrtr, 'pH_3D', ldstop = .FALSE. ) > 0 ) THEN
314         IF(lwp) WRITE(numout,*) 'Carbonate chem variable - reading in ...'
315         CALL iom_get( numrtr, jpdom_autoglo, 'pH_3D'   ,  f3_pH(:,:,:)     )
316         CALL iom_get( numrtr, jpdom_autoglo, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
317         CALL iom_get( numrtr, jpdom_autoglo, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
318         CALL iom_get( numrtr, jpdom_autoglo, 'CO3_3D'  ,  f3_co3(:,:,:)    )
319         CALL iom_get( numrtr, jpdom_autoglo, 'omcal_3D',  f3_omcal(:,:,:)  )
320         CALL iom_get( numrtr, jpdom_autoglo, 'omarg_3D',  f3_omarg(:,:,:)  )
321         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
322         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
323         !!
324         IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
325      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
326      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
327      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
328      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
329      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
330      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
331      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
332      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
333
334      ELSE
335         IF(lwp) WRITE(numout,*) 'WARNING : No Carbonate-chem variable in the restart.... '
336         IF(lwp) WRITE(numout,*) 'Is not a problem if start a month, but may be very problematic if not '
337         IF(lwp) WRITE(numout,*) 'Check if   mod(kt*rdt,2592000) == rdt' 
338         IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
339      ENDIF
340# endif
341
342
343#endif
344      !
345#if defined key_idtra
346      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
347      !!                        writting here undre their key.
348      !!                        problems in CFC restart, maybe because of this...
349      !!                        and pb in idtra diag or diad-restart writing.
350      !!----------------------------------------------------------------------
351      IF( iom_varid( numrtr, 'qint_IDTRA', ldstop = .FALSE. ) > 0 ) THEN
352         !! YES; in which case read them
353         !!
354         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties present - reading in ...'
355         CALL iom_get( numrtr, jpdom_autoglo, 'qint_IDTRA',  qint_idtra(:,:,1)  )
356      ELSE
357         !! NO; in which case set them to zero
358         !!
359         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties absent - setting to zero ...'
360         qint_idtra(:,:,1)  = 0.0   !! CHN
361      ENDIF
362      !!
363      !! calculate stats on these fields
364      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
365      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
366#endif
367      !
368#if defined key_cfc
369      DO jl = 1, jp_cfc
370         jn = jp_cfc0 + jl - 1
371         IF( iom_varid( numrtr, 'qint_'//ctrcnm(jn), ldstop = .FALSE. ) > 0 ) THEN
372            !! YES; in which case read them
373            !!
374            IF(lwp) WRITE(numout,*) ' CFC averaged properties present - reading in ...'
375            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
376         ELSE
377            !! NO; in which case set them to zero
378            !!
379            IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...'
380            qint_cfc(:,:,jn)  = 0.0   !! CHN
381         ENDIF
382         !!
383         !! calculate stats on these fields
384         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
385         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
386      END DO
387#endif
388      !
389   END SUBROUTINE trc_rst_read
390
391   SUBROUTINE trc_rst_wri( kt )
392      !!----------------------------------------------------------------------
393      !!                    ***  trc_rst_wri  ***
394      !!
395      !! ** purpose  :   write passive tracer fields in restart files
396      !!----------------------------------------------------------------------
397      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
398      !!
399      INTEGER  :: jn, jl
400      REAL(wp) :: zarak0
401      !! AXY (05/11/13): temporary variables
402      REAL(wp) ::    fq0,fq1,fq2
403      !!----------------------------------------------------------------------
404      !
405      CALL iom_rstput( kt, nitrst, numrtw, 'rdttrc1', rdttrc(1) )   ! surface passive tracer time step
406      ! prognostic variables
407      ! --------------------
408      DO jn = 1, jptra
409         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
410      END DO
411
412      DO jn = 1, jptra
413         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
414      END DO
415
416      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
417      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
418      !!                 version of NEMO date significantly earlier than the current
419      !!                 version
420
421#if defined key_medusa
422      !! AXY (13/01/12): write out "before" and "now" state of seafloor
423      !!                 sediment pools into restart; this happens
424      !!                 whether or not the pools are to be used by
425      !!                 MEDUSA (which is controlled by a switch in the
426      !!                 namelist_top file)
427      !!
428      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
429      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_N',  zb_sed_n(:,:)  )
430      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_N',  zn_sed_n(:,:)  )
431      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_FE', zb_sed_fe(:,:) )
432      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_FE', zn_sed_fe(:,:) )
433      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_SI', zb_sed_si(:,:) )
434      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_SI', zn_sed_si(:,:) )
435      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_C',  zb_sed_c(:,:)  )
436      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_C',  zn_sed_c(:,:)  )
437      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_CA', zb_sed_ca(:,:) )
438      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_CA', zn_sed_ca(:,:) )
439      !!
440      !! calculate stats on these fields
441      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
442      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
443      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
444      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
445      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
446      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
447      !!
448      !! AXY (07/07/15): write out temporally averaged fields for DMS
449      !!                 calculations
450      !!
451      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS - writing out ...'
452      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
453      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
454      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
455      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
456      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
457      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
458      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
459      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
460      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_DIN',  zb_dms_din(:,:)  )
461      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_DIN',  zn_dms_din(:,:)  )
462         !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
463         !!                  -- needed for the coupling with atm
464      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_srf',  zb_dms_srf(:,:)  )
465      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_srf',  zn_dms_srf(:,:)  )
466      CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
467      CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
468      !! JPALM 07-12-2017 -- To make things cleaner, we want to store an 
469      !!                     unscaled Chl field in the restart and only
470      !!                     scale it when reading it in.
471      CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf',  zn_chl_srf(:,:)  )
472      !!
473      !! calculate stats on these fields
474      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
475      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
476      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
477      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
478      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
479      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
480      call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
481      call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
482      call trc_rst_dia_stat(zn_chl_srf(:,:), 'unscaled CHL cpl')
483      !!
484      IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
485      call trc_rst_dia_stat(dust(:,:), 'Dust dep')
486      call trc_rst_dia_stat(zirondep(:,:), 'Iron dep')
487      !!
488      !! 
489      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
490      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
491# if defined key_roam
492      IF(lwp) WRITE(numout,*) 'Carbonate chem variable - writing out ...'
493      CALL iom_rstput( kt, nitrst, numrtw, 'pH_3D'   ,  f3_pH(:,:,:)     )
494      CALL iom_rstput( kt, nitrst, numrtw, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
495      CALL iom_rstput( kt, nitrst, numrtw, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
496      CALL iom_rstput( kt, nitrst, numrtw, 'CO3_3D'  ,  f3_co3(:,:,:)    )
497      CALL iom_rstput( kt, nitrst, numrtw, 'omcal_3D',  f3_omcal(:,:,:)  )
498      CALL iom_rstput( kt, nitrst, numrtw, 'omarg_3D',  f3_omarg(:,:,:)  )
499      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
500      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
501      !!
502      IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
503      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
504      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
505      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
506      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
507      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
508      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
509      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
510      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
511      !!
512# endif
513!!
514#endif
515      !
516#if defined key_idtra
517      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
518      !!                        writting here undre their key.
519      !!                        problems in CFC restart, maybe because of this...
520      !!                        and pb in idtra diag or diad-restart writing.
521      !!----------------------------------------------------------------------
522      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties - writing out ...'
523      CALL iom_rstput( kt, nitrst, numrtw, 'qint_IDTRA',  qint_idtra(:,:,1) )
524      !!
525      !! calculate stats on these fields
526      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
527      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
528#endif
529      !
530#if defined key_cfc
531      DO jl = 1, jp_cfc
532         jn = jp_cfc0 + jl - 1
533         IF(lwp) WRITE(numout,*) ' CFC averaged properties - writing out ...'
534         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
535         !!
536         !! calculate stats on these fields
537         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
538         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
539      END DO
540#endif
541      !
542
543      IF( kt == nitrst ) THEN
544          CALL trc_rst_stat            ! statistics
545#if defined key_medusa && defined key_roam
546          CALL trc_rst_conserve        ! conservation check
547#endif
548          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
549#if ! defined key_trdmxl_trc
550          lrst_trc = .FALSE.
551#endif
552          IF( lk_offline .AND. ln_rst_list ) THEN
553             nrst_lst = nrst_lst + 1
554             nitrst = nstocklist( nrst_lst )
555          ENDIF
556      ENDIF
557      !
558   END SUBROUTINE trc_rst_wri 
559
560
561   SUBROUTINE trc_rst_cal( kt, cdrw )
562      !!---------------------------------------------------------------------
563      !!                   ***  ROUTINE trc_rst_cal  ***
564      !!
565      !!  ** Purpose : Read or write calendar in restart file:
566      !!
567      !!  WRITE(READ) mode:
568      !!       kt        : number of time step since the begining of the experiment at the
569      !!                   end of the current(previous) run
570      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
571      !!                   end of the current(previous) run (REAL -> keep fractions of day)
572      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
573      !!
574      !!   According to namelist parameter nrstdt,
575      !!       nn_rsttr = 0  no control on the date (nittrc000 is  arbitrary).
576      !!       nn_rsttr = 1  we verify that nittrc000 is equal to the last
577      !!                   time step of previous run + 1.
578      !!       In both those options, the  exact duration of the experiment
579      !!       since the beginning (cumulated duration of all previous restart runs)
580      !!       is not stored in the restart and is assumed to be (nittrc000-1)*rdt.
581      !!       This is valid is the time step has remained constant.
582      !!
583      !!       nn_rsttr = 2  the duration of the experiment in days (adatrj)
584      !!                    has been stored in the restart file.
585      !!----------------------------------------------------------------------
586      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
587      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
588      !
589      INTEGER  ::  jlibalt = jprstlib
590      LOGICAL  ::  llok
591      REAL(wp) ::  zkt, zrdttrc1
592      REAL(wp) ::  zndastp
593
594      ! Time domain : restart
595      ! ---------------------
596
597      IF( TRIM(cdrw) == 'READ' ) THEN
598
599         IF(lwp) WRITE(numout,*)
600         IF(lwp) WRITE(numout,*) 'trc_rst_cal : read the TOP restart file for calendar'
601         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
602
603         IF ( jprstlib == jprstdimg ) THEN
604           ! eventually read netcdf file (monobloc)  for restarting on different number of processors
605           ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90
606           INQUIRE( FILE = TRIM(cn_trcrst_indir)//'/'//TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
607           IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
608         ENDIF
609
610         IF( ln_rsttr ) THEN
611            CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
612            CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
613
614            IF(lwp) THEN
615               WRITE(numout,*) ' *** Info read in restart : '
616               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
617               WRITE(numout,*) ' *** restart option'
618               SELECT CASE ( nn_rsttr )
619               CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
620               CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
621               CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
622               END SELECT
623               WRITE(numout,*)
624            ENDIF
625            ! Control of date
626            IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
627               &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
628               &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
629         ENDIF
630         !
631         IF( lk_offline ) THEN   
632            !                                          ! set the date in offline mode
633            IF( ln_rsttr .AND. nn_rsttr == 2 ) THEN
634               CALL iom_get( numrtr, 'ndastp', zndastp ) 
635               ndastp = NINT( zndastp )
636               CALL iom_get( numrtr, 'adatrj', adatrj  )
637             ELSE
638               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
639               adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
640               ! note this is wrong if time step has changed during run
641            ENDIF
642            !
643            IF(lwp) THEN
644              WRITE(numout,*) ' *** Info used values : '
645              WRITE(numout,*) '   date ndastp                                      : ', ndastp
646              WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
647              WRITE(numout,*)
648            ENDIF
649            !
650            IF( ln_rsttr )  THEN   ;    neuler = 1
651            ELSE                   ;    neuler = 0
652            ENDIF
653            !
654            CALL day_init          ! compute calendar
655            !
656         ENDIF
657         !
658      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
659         !
660         IF(  kt == nitrst ) THEN
661            IF(lwp) WRITE(numout,*)
662            IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
663            IF(lwp) WRITE(numout,*) '~~~~~~~'
664         ENDIF
665         CALL iom_rstput( kt, nitrst, numrtw, 'kt'     , REAL( kt    , wp) )   ! time-step
666         CALL iom_rstput( kt, nitrst, numrtw, 'ndastp' , REAL( ndastp, wp) )   ! date
667         CALL iom_rstput( kt, nitrst, numrtw, 'adatrj' , adatrj            )   ! number of elapsed days since
668         !                                                                     ! the begining of the run [s]
669      ENDIF
670
671   END SUBROUTINE trc_rst_cal
672
673
674   SUBROUTINE trc_rst_stat
675      !!----------------------------------------------------------------------
676      !!                    ***  trc_rst_stat  ***
677      !!
678      !! ** purpose  :   Compute tracers statistics
679      !!----------------------------------------------------------------------
680      INTEGER  :: jk, jn
681      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
682      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
683      !!----------------------------------------------------------------------
684
685      IF( lwp ) THEN
686         WRITE(numout,*) 
687         WRITE(numout,*) '           ----TRACER STAT----             '
688         WRITE(numout,*) 
689      ENDIF
690      !
691      DO jk = 1, jpk
692         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
693      END DO
694      !
695      DO jn = 1, jptra
696         ztraf = glob_sum( trn(:,:,:,jn) * zvol(:,:,:) )
697         zmin  = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
698         zmax  = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
699         IF( lk_mpp ) THEN
700            CALL mpp_min( zmin )      ! min over the global domain
701            CALL mpp_max( zmax )      ! max over the global domain
702         END IF
703         zmean  = ztraf / areatot
704         zdrift = ( ( ztraf - trai(jn) ) / ( trai(jn) + 1.e-12 )  ) * 100._wp
705         IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
706      END DO
707      IF(lwp) WRITE(numout,*) 
7089000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
709      &      '    max :',e18.10,'    drift :',e18.10, ' %')
710      !
711   END SUBROUTINE trc_rst_stat
712
713
714# if defined key_medusa && defined key_roam
715   SUBROUTINE trc_rst_conserve
716      !!----------------------------------------------------------------------
717      !!                    ***  trc_rst_conserve  ***
718      !!
719      !! ** purpose  :   Compute tracers conservation statistics
720      !!
721      !! AXY (17/11/2017)
722      !! This routine calculates the "now" inventories of the elemental
723      !! cycles of MEDUSA and compares them to those calculate when the
724      !! model was initialised / restarted; the cycles calculated are:
725      !!    nitrogen, silicon, iron, carbon, alkalinity and oxygen
726      !!----------------------------------------------------------------------
727      INTEGER  :: ji, jj, jk, jn
728      REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio, loc_vol, loc_are
729      REAL(wp) :: zq1, zq2, loc_vol, loc_area
730      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol
731      REAL(wp), DIMENSION(jpi,jpj)     :: z2d, zarea
732      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
733      !!----------------------------------------------------------------------
734      !
735      IF( lwp ) THEN
736         WRITE(numout,*) 
737         WRITE(numout,*) '           ----TRACER CONSERVATION----             '
738         WRITE(numout,*) 
739      ENDIF
740      !
741      ! ocean volume
742      DO jk = 1, jpk
743         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
744      END DO
745      !
746      ! ocean area (for sediments)
747      zarea(:,:)      = e1e2t(:,:) * tmask(:,:,1)
748      !
749      !----------------------------------------------------------------------
750      ! nitrogen
751      z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
752                   trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
753      z2d(:,:)   = zn_sed_n(:,:)
754      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
755      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
756      ! total tracer, and delta
757      zinvt      = zsum3d + zsum2d
758      zdelta     = zinvt - cycletot(1)
759      zratio     = 1.0e2 * zdelta / cycletot(1)
760      !
761      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt,   &
762         cycletot(1), zdelta, zratio
763      IF ( lwp ) WRITE(numout,*) 
764      !
765      !----------------------------------------------------------------------
766      ! silicon
767      z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
768      z2d(:,:)   = zn_sed_si(:,:)
769      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
770      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
771      ! total tracer, and delta
772      zinvt      = zsum3d + zsum2d
773      zdelta     = zinvt - cycletot(2)
774      zratio     = 1.0e2 * zdelta / cycletot(2)
775      !
776      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt,    &
777         cycletot(2), zdelta, zratio
778      IF ( lwp ) WRITE(numout,*) 
779      !
780      !----------------------------------------------------------------------
781      ! iron
782      z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
783            trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
784      z2d(:,:)   = zn_sed_fe(:,:)
785      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
786      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
787      ! total tracer, and delta
788      zinvt      = zsum3d + zsum2d
789      zdelta     = zinvt - cycletot(3)
790      zratio     = 1.0e2 * zdelta / cycletot(3)
791      !
792      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt,       &
793         cycletot(3), zdelta, zratio
794      IF ( lwp ) WRITE(numout,*) 
795      !
796      !----------------------------------------------------------------------
797      ! carbon
798      z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  + &
799                   (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + &
800                   trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
801      z2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
802      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
803      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
804      ! total tracer, and delta
805      zinvt      = zsum3d + zsum2d
806      zdelta     = zinvt - cycletot(4)
807      zratio     = 1.0e2 * zdelta / cycletot(4)
808      !
809      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt,     &
810         cycletot(4), zdelta, zratio
811      IF ( lwp ) WRITE(numout,*) 
812      !
813      !----------------------------------------------------------------------
814      ! alkalinity
815      z3d(:,:,:) = trn(:,:,:,jpalk)
816      z2d(:,:)   = zn_sed_ca(:,:) * 2.0
817      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
818      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
819      ! total tracer, and delta
820      zinvt      = zsum3d + zsum2d
821      zdelta     = zinvt - cycletot(5)
822      zratio     = 1.0e2 * zdelta / cycletot(5)
823      !
824      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, &
825         cycletot(5), zdelta, zratio
826      IF ( lwp ) WRITE(numout,*) 
827      !
828      !----------------------------------------------------------------------
829      ! oxygen
830      z3d(:,:,:) = trn(:,:,:,jpoxy)
831      z2d(:,:)   = 0.0
832      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
833      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
834      ! total tracer, and delta
835      zinvt      = zsum3d + zsum2d
836      zdelta     = zinvt - cycletot(6)
837      zratio     = 1.0e2 * zdelta / cycletot(6)
838      !
839      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt,     &
840         cycletot(6), zdelta, zratio
841      !
842      !----------------------------------------------------------------------
843      ! Check
844      zsum3d        = glob_sum( zvol(:,:,:) )
845      zsum2d        = glob_sum( zarea(:,:) )
846      IF ( lwp ) THEN
847         WRITE(numout,*)
848         WRITE(numout,*) ' check : cvol    : ', zsum3d
849         WRITE(numout,*) ' check : carea   : ', zsum2d
850         WRITE(numout,*)
851      ENDIF
852      !
8539010  FORMAT(' element:',a10,                     &
854             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
855             ' total:',e18.10,' initial:',e18.10, &
856             ' delta:',e18.10,' %:',e18.10)
857      !
858   END SUBROUTINE trc_rst_conserve 
859# endif
860
861
862#else
863   !!----------------------------------------------------------------------
864   !!  Dummy module :                                     No passive tracer
865   !!----------------------------------------------------------------------
866CONTAINS
867   SUBROUTINE trc_rst_read                      ! Empty routines
868   END SUBROUTINE trc_rst_read
869   SUBROUTINE trc_rst_wri( kt )
870      INTEGER, INTENT ( in ) :: kt
871      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
872   END SUBROUTINE trc_rst_wri   
873#endif
874
875   !!----------------------------------------------------------------------
876   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
877   !! $Id$
878   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
879   !!======================================================================
880END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.