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 @ 9129

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

JPALM -- merge with GO6_package rev 8638-9114

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