source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 8462

Last change on this file since 8462 was 8462, checked in by dford, 3 years ago

Mask land if initialising pCO2/fCO2.

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