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

Last change on this file since 8427 was 8427, checked in by timgraham, 3 years ago

GMED ticket 341. Fixes found through use of preset Na Ns? to make the code more robust.

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