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.
trcrst.F90 in branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 9292

Last change on this file since 9292 was 9292, checked in by dford, 6 years ago

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