New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcrst.F90 in branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 9128

Last change on this file since 9128 was 9128, checked in by jpalmier, 6 years ago

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