source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 10149

Last change on this file since 10149 was 10149, checked in by frrh, 23 months ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

svn merge -r 10054:10141 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3

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