source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 13354

Last change on this file since 13354 was 13354, checked in by jwhile, 7 weeks ago

Merged in trunk

File size: 44.1 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         IF( iom_varid( numrtr, 'pgrow_avg_3d', ldstop = .FALSE. ) > 0 ) THEN
371            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...'
372            CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg_3d',  pgrow_avg_3d(:,:,:)  )
373            CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg_3d',  ploss_avg_3d(:,:,:)  )
374            CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg_3d',   phyt_avg_3d(:,:,:)   )
375         ELSE
376            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg_3d absent - setting to zero ...'
377            pgrow_avg_3d(:,:,:) = 0.0
378            ploss_avg_3d(:,:,:) = 0.0
379            phyt_avg_3d(:,:,:)  = 0.0
380         ENDIF
381      ENDIF
382
383#endif
384      !
385#if defined key_idtra
386      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
387      !!                        writting here undre their key.
388      !!                        problems in CFC restart, maybe because of this...
389      !!                        and pb in idtra diag or diad-restart writing.
390      !!----------------------------------------------------------------------
391      IF( iom_varid( numrtr, 'qint_IDTRA', ldstop = .FALSE. ) > 0 ) THEN
392         !! YES; in which case read them
393         !!
394         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties present - reading in ...'
395         CALL iom_get( numrtr, jpdom_autoglo, 'qint_IDTRA',  qint_idtra(:,:,1)  )
396      ELSE
397         !! NO; in which case set them to zero
398         !!
399         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties absent - setting to zero ...'
400         qint_idtra(:,:,1)  = 0.0   !! CHN
401      ENDIF
402      !!
403      !! calculate stats on these fields
404      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
405      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
406#endif
407      !
408#if defined key_cfc
409      DO jl = 1, jp_cfc
410         jn = jp_cfc0 + jl - 1
411         IF( iom_varid( numrtr, 'qint_'//ctrcnm(jn), ldstop = .FALSE. ) > 0 ) THEN
412            !! YES; in which case read them
413            !!
414            IF(lwp) WRITE(numout,*) ' CFC averaged properties present - reading in ...'
415            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
416         ELSE
417            !! NO; in which case set them to zero
418            !!
419            IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...'
420            qint_cfc(:,:,jl)  = 0.0   !! CHN
421         ENDIF
422         !!
423         !! calculate stats on these fields
424         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
425         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
426      END DO
427#endif
428      !
429   END SUBROUTINE trc_rst_read
430
431   SUBROUTINE trc_rst_wri( kt )
432      !!----------------------------------------------------------------------
433      !!                    ***  trc_rst_wri  ***
434      !!
435      !! ** purpose  :   write passive tracer fields in restart files
436      !!----------------------------------------------------------------------
437      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
438      !!
439      INTEGER  :: jn, jl
440      REAL(wp) :: zarak0
441      !! AXY (05/11/13): temporary variables
442      REAL(wp) ::    fq0,fq1,fq2
443      !!----------------------------------------------------------------------
444      !
445      CALL iom_rstput( kt, nitrst, numrtw, 'rdttrc1', rdttrc(1) )   ! surface passive tracer time step
446      ! prognostic variables
447      ! --------------------
448      DO jn = 1, jptra
449         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
450      END DO
451
452      DO jn = 1, jptra
453         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
454      END DO
455
456      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
457      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
458      !!                 version of NEMO date significantly earlier than the current
459      !!                 version
460
461#if defined key_medusa
462      !! AXY (13/01/12): write out "before" and "now" state of seafloor
463      !!                 sediment pools into restart; this happens
464      !!                 whether or not the pools are to be used by
465      !!                 MEDUSA (which is controlled by a switch in the
466      !!                 namelist_top file)
467      !!
468      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
469      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_N',  zb_sed_n(:,:)  )
470      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_N',  zn_sed_n(:,:)  )
471      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_FE', zb_sed_fe(:,:) )
472      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_FE', zn_sed_fe(:,:) )
473      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_SI', zb_sed_si(:,:) )
474      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_SI', zn_sed_si(:,:) )
475      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_C',  zb_sed_c(:,:)  )
476      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_C',  zn_sed_c(:,:)  )
477      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_CA', zb_sed_ca(:,:) )
478      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_CA', zn_sed_ca(:,:) )
479      !!
480      !! calculate stats on these fields
481      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
482      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
483      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
484      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
485      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
486      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
487      !!
488      !! AXY (07/07/15): write out temporally averaged fields for DMS
489      !!                 calculations
490      !!
491      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS - writing out ...'
492      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
493      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
494      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
495      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
496      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
497      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
498      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
499      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
500      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_DIN',  zb_dms_din(:,:)  )
501      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_DIN',  zn_dms_din(:,:)  )
502         !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
503         !!                  -- needed for the coupling with atm
504      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_srf',  zb_dms_srf(:,:)  )
505      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_srf',  zn_dms_srf(:,:)  )
506      CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
507      CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
508      !! JPALM 07-12-2017 -- To make things cleaner, we want to store an 
509      !!                     unscaled Chl field in the restart and only
510      !!                     scale it when reading it in.
511      CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf',  zn_chl_srf(:,:)  )
512      !!
513      !! calculate stats on these fields
514      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
515      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
516      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
517      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
518      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
519      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
520      call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
521      call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
522      call trc_rst_dia_stat(zn_chl_srf(:,:), 'unscaled CHL cpl')
523      !!
524      IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
525      call trc_rst_dia_stat(dust(:,:), 'Dust dep')
526      call trc_rst_dia_stat(zirondep(:,:), 'Iron dep')
527      !!
528      !! 
529      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
530      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
531# if defined key_roam
532      IF(lwp) WRITE(numout,*) 'Carbonate chem variable - writing out ...'
533      CALL iom_rstput( kt, nitrst, numrtw, 'pH_3D'   ,  f3_pH(:,:,:)     )
534      CALL iom_rstput( kt, nitrst, numrtw, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
535      CALL iom_rstput( kt, nitrst, numrtw, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
536      CALL iom_rstput( kt, nitrst, numrtw, 'CO3_3D'  ,  f3_co3(:,:,:)    )
537      CALL iom_rstput( kt, nitrst, numrtw, 'omcal_3D',  f3_omcal(:,:,:)  )
538      CALL iom_rstput( kt, nitrst, numrtw, 'omarg_3D',  f3_omarg(:,:,:)  )
539      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
540      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
541      !!
542      IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
543      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
544      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
545      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
546      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
547      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
548      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
549      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
550      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
551      !!
552      IF ( ln_foam_medusa ) THEN
553         !! Fields for observation operator on first timestep
554         IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...'
555         CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:)  )
556         CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:)  )
557      ENDIF
558# endif
559      IF ( ln_foam_medusa ) THEN
560         !! Fields for assimilation on first timestep
561         IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...'
562         CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) )
563         CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) )
564         CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  )
565         CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   )
566         CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg_3d', pgrow_avg_3d(:,:,:) )
567         CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg_3d', ploss_avg_3d(:,:,:) )
568         CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg_3d',  phyt_avg_3d(:,:,:)  )
569      ENDIF
570!!
571#endif
572      !
573#if defined key_idtra
574      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
575      !!                        writting here undre their key.
576      !!                        problems in CFC restart, maybe because of this...
577      !!                        and pb in idtra diag or diad-restart writing.
578      !!----------------------------------------------------------------------
579      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties - writing out ...'
580      CALL iom_rstput( kt, nitrst, numrtw, 'qint_IDTRA',  qint_idtra(:,:,1) )
581      !!
582      !! calculate stats on these fields
583      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
584      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
585#endif
586      !
587#if defined key_cfc
588      DO jl = 1, jp_cfc
589         jn = jp_cfc0 + jl - 1
590         IF(lwp) WRITE(numout,*) ' CFC averaged properties - writing out ...'
591         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
592         !!
593         !! calculate stats on these fields
594         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
595         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
596      END DO
597#endif
598      !
599
600      IF( kt == nitrst ) THEN
601          CALL trc_rst_stat            ! statistics
602#if defined key_medusa && defined key_roam
603          CALL trc_rst_conserve        ! conservation check
604#endif
605          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
606#if ! defined key_trdmxl_trc
607          lrst_trc = .FALSE.
608#endif
609          IF( lk_offline .AND. ln_rst_list ) THEN
610             nrst_lst = nrst_lst + 1
611             nitrst = nstocklist( nrst_lst )
612          ENDIF
613      ENDIF
614      !
615   END SUBROUTINE trc_rst_wri 
616
617
618   SUBROUTINE trc_rst_cal( kt, cdrw )
619      !!---------------------------------------------------------------------
620      !!                   ***  ROUTINE trc_rst_cal  ***
621      !!
622      !!  ** Purpose : Read or write calendar in restart file:
623      !!
624      !!  WRITE(READ) mode:
625      !!       kt        : number of time step since the begining of the experiment at the
626      !!                   end of the current(previous) run
627      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
628      !!                   end of the current(previous) run (REAL -> keep fractions of day)
629      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
630      !!
631      !!   According to namelist parameter nrstdt,
632      !!       nn_rsttr = 0  no control on the date (nittrc000 is  arbitrary).
633      !!       nn_rsttr = 1  we verify that nittrc000 is equal to the last
634      !!                   time step of previous run + 1.
635      !!       In both those options, the  exact duration of the experiment
636      !!       since the beginning (cumulated duration of all previous restart runs)
637      !!       is not stored in the restart and is assumed to be (nittrc000-1)*rdt.
638      !!       This is valid is the time step has remained constant.
639      !!
640      !!       nn_rsttr = 2  the duration of the experiment in days (adatrj)
641      !!                    has been stored in the restart file.
642      !!----------------------------------------------------------------------
643      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
644      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
645      !
646      INTEGER  ::  jlibalt = jprstlib
647      LOGICAL  ::  llok
648      REAL(wp) ::  zkt, zrdttrc1
649      REAL(wp) ::  zndastp
650
651      ! Time domain : restart
652      ! ---------------------
653
654      IF( TRIM(cdrw) == 'READ' ) THEN
655
656         IF(lwp) WRITE(numout,*)
657         IF(lwp) WRITE(numout,*) 'trc_rst_cal : read the TOP restart file for calendar'
658         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
659
660         IF ( jprstlib == jprstdimg ) THEN
661           ! eventually read netcdf file (monobloc)  for restarting on different number of processors
662           ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90
663           INQUIRE( FILE = TRIM(cn_trcrst_indir)//'/'//TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
664           IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
665         ENDIF
666
667         IF( ln_rsttr ) THEN
668            CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
669            CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
670
671            IF(lwp) THEN
672               WRITE(numout,*) ' *** Info read in restart : '
673               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
674               WRITE(numout,*) ' *** restart option'
675               SELECT CASE ( nn_rsttr )
676               CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
677               CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
678               CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
679               END SELECT
680               WRITE(numout,*)
681            ENDIF
682            ! Control of date
683            IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
684               &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
685               &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
686         ENDIF
687         !
688         IF( lk_offline ) THEN   
689            !                                          ! set the date in offline mode
690            IF( ln_rsttr .AND. nn_rsttr == 2 ) THEN
691               CALL iom_get( numrtr, 'ndastp', zndastp ) 
692               ndastp = NINT( zndastp )
693               CALL iom_get( numrtr, 'adatrj', adatrj  )
694             ELSE
695               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
696               adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
697               ! note this is wrong if time step has changed during run
698            ENDIF
699            !
700            IF(lwp) THEN
701              WRITE(numout,*) ' *** Info used values : '
702              WRITE(numout,*) '   date ndastp                                      : ', ndastp
703              WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
704              WRITE(numout,*)
705            ENDIF
706            !
707            IF( ln_rsttr )  THEN   ;    neuler = 1
708            ELSE                   ;    neuler = 0
709            ENDIF
710            !
711            CALL day_init          ! compute calendar
712            !
713         ENDIF
714         !
715      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
716         !
717         IF(  kt == nitrst ) THEN
718            IF(lwp) WRITE(numout,*)
719            IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
720            IF(lwp) WRITE(numout,*) '~~~~~~~'
721         ENDIF
722         CALL iom_rstput( kt, nitrst, numrtw, 'kt'     , REAL( kt    , wp) )   ! time-step
723         CALL iom_rstput( kt, nitrst, numrtw, 'ndastp' , REAL( ndastp, wp) )   ! date
724         CALL iom_rstput( kt, nitrst, numrtw, 'adatrj' , adatrj            )   ! number of elapsed days since
725         !                                                                     ! the begining of the run [s]
726      ENDIF
727
728   END SUBROUTINE trc_rst_cal
729
730
731   SUBROUTINE trc_rst_stat
732      !!----------------------------------------------------------------------
733      !!                    ***  trc_rst_stat  ***
734      !!
735      !! ** purpose  :   Compute tracers statistics
736      !!----------------------------------------------------------------------
737      INTEGER  :: jk, jn
738      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
739      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
740      !!----------------------------------------------------------------------
741
742      IF( lwp ) THEN
743         WRITE(numout,*) 
744         WRITE(numout,*) '           ----TRACER STAT----             '
745         WRITE(numout,*) 
746      ENDIF
747      !
748      DO jk = 1, jpk
749         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
750      END DO
751      !
752      DO jn = 1, jptra
753         ztraf = glob_sum( trn(:,:,:,jn) * zvol(:,:,:) )
754         zmin  = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
755         zmax  = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
756         IF( lk_mpp ) THEN
757            CALL mpp_min( zmin )      ! min over the global domain
758            CALL mpp_max( zmax )      ! max over the global domain
759         END IF
760         zmean  = ztraf / areatot
761         zdrift = ( ( ztraf - trai(jn) ) / ( trai(jn) + 1.e-12 )  ) * 100._wp
762         IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
763      END DO
764      IF(lwp) WRITE(numout,*) 
7659000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
766      &      '    max :',e18.10,'    drift :',e18.10, ' %')
767      !
768   END SUBROUTINE trc_rst_stat
769
770
771# if defined key_medusa && defined key_roam
772   SUBROUTINE trc_rst_conserve
773      !!----------------------------------------------------------------------
774      !!                    ***  trc_rst_conserve  ***
775      !!
776      !! ** purpose  :   Compute tracers conservation statistics
777      !!
778      !! AXY (17/11/2017)
779      !! This routine calculates the "now" inventories of the elemental
780      !! cycles of MEDUSA and compares them to those calculate when the
781      !! model was initialised / restarted; the cycles calculated are:
782      !!    nitrogen, silicon, iron, carbon, alkalinity and oxygen
783      !!----------------------------------------------------------------------
784      INTEGER  :: ji, jj, jk, jn
785      REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio
786      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol
787      REAL(wp), DIMENSION(jpi,jpj)     :: z2d, zarea
788      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
789      !!----------------------------------------------------------------------
790      !
791      IF( lwp ) THEN
792         WRITE(numout,*) 
793         WRITE(numout,*) '           ----TRACER CONSERVATION----             '
794         WRITE(numout,*) 
795      ENDIF
796      !
797      ! ocean volume
798      DO jk = 1, jpk
799         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
800      END DO
801      !
802      ! ocean area (for sediments)
803      zarea(:,:)      = e1e2t(:,:) * tmask(:,:,1)
804      !
805      !----------------------------------------------------------------------
806      ! nitrogen
807      z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
808                   trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
809      z2d(:,:)   = zn_sed_n(:,:)
810      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
811      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
812      ! total tracer, and delta
813      zinvt      = zsum3d + zsum2d
814      zdelta     = zinvt - cycletot(1)
815      zratio     = 1.0e2 * zdelta / cycletot(1)
816      !
817      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt,   &
818         cycletot(1), zdelta, zratio
819      IF ( lwp ) WRITE(numout,*) 
820      !
821      !----------------------------------------------------------------------
822      ! silicon
823      z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
824      z2d(:,:)   = zn_sed_si(:,:)
825      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
826      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
827      ! total tracer, and delta
828      zinvt      = zsum3d + zsum2d
829      zdelta     = zinvt - cycletot(2)
830      zratio     = 1.0e2 * zdelta / cycletot(2)
831      !
832      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt,    &
833         cycletot(2), zdelta, zratio
834      IF ( lwp ) WRITE(numout,*) 
835      !
836      !----------------------------------------------------------------------
837      ! iron
838      z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
839            trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
840      z2d(:,:)   = zn_sed_fe(:,:)
841      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
842      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
843      ! total tracer, and delta
844      zinvt      = zsum3d + zsum2d
845      zdelta     = zinvt - cycletot(3)
846      zratio     = 1.0e2 * zdelta / cycletot(3)
847      !
848      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt,       &
849         cycletot(3), zdelta, zratio
850      IF ( lwp ) WRITE(numout,*) 
851      !
852      !----------------------------------------------------------------------
853      ! carbon
854      z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  + &
855                   (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + &
856                   trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
857      z2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
858      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
859      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
860      ! total tracer, and delta
861      zinvt      = zsum3d + zsum2d
862      zdelta     = zinvt - cycletot(4)
863      zratio     = 1.0e2 * zdelta / cycletot(4)
864      !
865      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt,     &
866         cycletot(4), zdelta, zratio
867      IF ( lwp ) WRITE(numout,*) 
868      !
869      !----------------------------------------------------------------------
870      ! alkalinity
871      z3d(:,:,:) = trn(:,:,:,jpalk)
872      z2d(:,:)   = zn_sed_ca(:,:) * 2.0
873      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
874      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
875      ! total tracer, and delta
876      zinvt      = zsum3d + zsum2d
877      zdelta     = zinvt - cycletot(5)
878      zratio     = 1.0e2 * zdelta / cycletot(5)
879      !
880      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, &
881         cycletot(5), zdelta, zratio
882      IF ( lwp ) WRITE(numout,*) 
883      !
884      !----------------------------------------------------------------------
885      ! oxygen
886      z3d(:,:,:) = trn(:,:,:,jpoxy)
887      z2d(:,:)   = 0.0
888      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
889      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
890      ! total tracer, and delta
891      zinvt      = zsum3d + zsum2d
892      zdelta     = zinvt - cycletot(6)
893      zratio     = 1.0e2 * zdelta / cycletot(6)
894      !
895      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt,     &
896         cycletot(6), zdelta, zratio
897      !
898      !----------------------------------------------------------------------
899      ! Check
900      zsum3d        = glob_sum( zvol(:,:,:) )
901      zsum2d        = glob_sum( zarea(:,:) )
902      IF ( lwp ) THEN
903         WRITE(numout,*)
904         WRITE(numout,*) ' check : cvol    : ', zsum3d
905         WRITE(numout,*) ' check : carea   : ', zsum2d
906         WRITE(numout,*)
907      ENDIF
908      !
9099010  FORMAT(' element:',a10,                     &
910             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
911             ' total:',e18.10,' initial:',e18.10, &
912             ' delta:',e18.10,' %:',e18.10)
913      !
914   END SUBROUTINE trc_rst_conserve 
915# endif
916
917
918#else
919   !!----------------------------------------------------------------------
920   !!  Dummy module :                                     No passive tracer
921   !!----------------------------------------------------------------------
922CONTAINS
923   SUBROUTINE trc_rst_read                      ! Empty routines
924   END SUBROUTINE trc_rst_read
925   SUBROUTINE trc_rst_wri( kt )
926      INTEGER, INTENT ( in ) :: kt
927      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
928   END SUBROUTINE trc_rst_wri   
929#endif
930
931   !!----------------------------------------------------------------------
932   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
933   !! $Id$
934   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
935   !!======================================================================
936END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.