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_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 7766

Last change on this file since 7766 was 7766, checked in by jpalmier, 7 years ago

JPALM -- adapted MEDUSA_SBC to coupled mode, and merged this branch with RICHARD s fixes_part2 branch

File size: 36.5 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  !! 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, 'B_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      !! calculate stats on these fields
276      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
277      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
278      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
279      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
280      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
281      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
282      call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
283      call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
284      !! 
285      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
286      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
287# if defined key_roam
288      IF( iom_varid( numrtr, 'pH_3D', ldstop = .FALSE. ) > 0 ) THEN
289         IF(lwp) WRITE(numout,*) 'Carbonate chem variable - reading in ...'
290         CALL iom_get( numrtr, jpdom_autoglo, 'pH_3D'   ,  f3_pH(:,:,:)     )
291         CALL iom_get( numrtr, jpdom_autoglo, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
292         CALL iom_get( numrtr, jpdom_autoglo, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
293         CALL iom_get( numrtr, jpdom_autoglo, 'CO3_3D'  ,  f3_co3(:,:,:)    )
294         CALL iom_get( numrtr, jpdom_autoglo, 'omcal_3D',  f3_omcal(:,:,:)  )
295         CALL iom_get( numrtr, jpdom_autoglo, 'omarg_3D',  f3_omarg(:,:,:)  )
296         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
297         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
298         !!
299         IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
300      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
301      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
302      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
303      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
304      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
305      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
306      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
307      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
308
309      ELSE
310         IF(lwp) WRITE(numout,*) 'WARNING : No Carbonate-chem variable in the restart.... '
311         IF(lwp) WRITE(numout,*) 'Is not a problem if start a month, but may be very problematic if not '
312         IF(lwp) WRITE(numout,*) 'Check if   mod(kt*rdt,2592000) == rdt' 
313         IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
314      ENDIF
315# endif
316
317
318#endif
319      !
320#if defined key_idtra
321      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
322      !!                        writting here undre their key.
323      !!                        problems in CFC restart, maybe because of this...
324      !!                        and pb in idtra diag or diad-restart writing.
325      !!----------------------------------------------------------------------
326      IF( iom_varid( numrtr, 'qint_IDTRA', ldstop = .FALSE. ) > 0 ) THEN
327         !! YES; in which case read them
328         !!
329         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties present - reading in ...'
330         CALL iom_get( numrtr, jpdom_autoglo, 'qint_IDTRA',  qint_idtra(:,:,1)  )
331      ELSE
332         !! NO; in which case set them to zero
333         !!
334         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties absent - setting to zero ...'
335         qint_idtra(:,:,1)  = 0.0   !! CHN
336      ENDIF
337      !!
338      !! calculate stats on these fields
339      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
340      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
341#endif
342      !
343#if defined key_cfc
344      DO jl = 1, jp_cfc
345         jn = jp_cfc0 + jl - 1
346         IF( iom_varid( numrtr, 'qint_'//ctrcnm(jn), ldstop = .FALSE. ) > 0 ) THEN
347            !! YES; in which case read them
348            !!
349            IF(lwp) WRITE(numout,*) ' CFC averaged properties present - reading in ...'
350            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
351         ELSE
352            !! NO; in which case set them to zero
353            !!
354            IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...'
355            qint_cfc(:,:,jn)  = 0.0   !! CHN
356         ENDIF
357         !!
358         !! calculate stats on these fields
359         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
360         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
361      END DO
362#endif
363      !
364   END SUBROUTINE trc_rst_read
365
366   SUBROUTINE trc_rst_wri( kt )
367      !!----------------------------------------------------------------------
368      !!                    ***  trc_rst_wri  ***
369      !!
370      !! ** purpose  :   write passive tracer fields in restart files
371      !!----------------------------------------------------------------------
372      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
373      !!
374      INTEGER  :: jn, jl
375      REAL(wp) :: zarak0
376      !! AXY (05/11/13): temporary variables
377      REAL(wp) ::    fq0,fq1,fq2
378      !!----------------------------------------------------------------------
379      !
380      CALL iom_rstput( kt, nitrst, numrtw, 'rdttrc1', rdttrc(1) )   ! surface passive tracer time step
381      ! prognostic variables
382      ! --------------------
383      DO jn = 1, jptra
384         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
385      END DO
386
387      DO jn = 1, jptra
388         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
389      END DO
390
391      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
392      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
393      !!                 version of NEMO date significantly earlier than the current
394      !!                 version
395
396#if defined key_medusa
397      !! AXY (13/01/12): write out "before" and "now" state of seafloor
398      !!                 sediment pools into restart; this happens
399      !!                 whether or not the pools are to be used by
400      !!                 MEDUSA (which is controlled by a switch in the
401      !!                 namelist_top file)
402      !!
403      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
404      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_N',  zb_sed_n(:,:)  )
405      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_N',  zn_sed_n(:,:)  )
406      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_FE', zb_sed_fe(:,:) )
407      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_FE', zn_sed_fe(:,:) )
408      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_SI', zb_sed_si(:,:) )
409      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_SI', zn_sed_si(:,:) )
410      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_C',  zb_sed_c(:,:)  )
411      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_C',  zn_sed_c(:,:)  )
412      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_CA', zb_sed_ca(:,:) )
413      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_CA', zn_sed_ca(:,:) )
414      !!
415      !! calculate stats on these fields
416      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
417      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
418      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
419      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
420      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
421      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
422      !!
423      !! AXY (07/07/15): write out temporally averaged fields for DMS
424      !!                 calculations
425      !!
426      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS - writing out ...'
427      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
428      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
429      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
430      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
431      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
432      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
433      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
434      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
435      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_DIN',  zb_dms_din(:,:)  )
436      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_DIN',  zn_dms_din(:,:)  )
437         !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
438         !!                  -- needed for the coupling with atm
439      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_srf',  zb_dms_srf(:,:)  )
440      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_srf',  zn_dms_srf(:,:)  )
441      CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
442      CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
443      !!
444      !! calculate stats on these fields
445      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
446      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
447      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
448      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
449      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
450      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
451      call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
452      call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
453      !!
454      IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
455      call trc_rst_dia_stat(dust(:,:), 'Dust dep')
456      call trc_rst_dia_stat(zirondep(:,:), 'Iron dep')
457      !!
458      !! 
459      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
460      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
461# if defined key_roam
462      IF(lwp) WRITE(numout,*) 'Carbonate chem variable - writing out ...'
463      CALL iom_rstput( kt, nitrst, numrtw, 'pH_3D'   ,  f3_pH(:,:,:)     )
464      CALL iom_rstput( kt, nitrst, numrtw, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
465      CALL iom_rstput( kt, nitrst, numrtw, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
466      CALL iom_rstput( kt, nitrst, numrtw, 'CO3_3D'  ,  f3_co3(:,:,:)    )
467      CALL iom_rstput( kt, nitrst, numrtw, 'omcal_3D',  f3_omcal(:,:,:)  )
468      CALL iom_rstput( kt, nitrst, numrtw, 'omarg_3D',  f3_omarg(:,:,:)  )
469      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
470      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
471      !!
472      IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
473      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
474      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
475      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
476      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
477      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
478      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
479      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
480      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
481      !!
482# endif
483!!
484#endif
485      !
486#if defined key_idtra
487      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and
488      !!                        writting here undre their key.
489      !!                        problems in CFC restart, maybe because of this...
490      !!                        and pb in idtra diag or diad-restart writing.
491      !!----------------------------------------------------------------------
492      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties - writing out ...'
493      CALL iom_rstput( kt, nitrst, numrtw, 'qint_IDTRA',  qint_idtra(:,:,1) )
494      !!
495      !! calculate stats on these fields
496      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
497      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
498#endif
499      !
500#if defined key_cfc
501      DO jl = 1, jp_cfc
502         jn = jp_cfc0 + jl - 1
503         IF(lwp) WRITE(numout,*) ' CFC averaged properties - writing out ...'
504         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
505         !!
506         !! calculate stats on these fields
507         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
508         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
509      END DO
510#endif
511      !
512
513      IF( kt == nitrst ) THEN
514          CALL trc_rst_stat            ! statistics
515          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
516#if ! defined key_trdmxl_trc
517          lrst_trc = .FALSE.
518#endif
519          IF( lk_offline .AND. ln_rst_list ) THEN
520             nrst_lst = nrst_lst + 1
521             nitrst = nstocklist( nrst_lst )
522          ENDIF
523      ENDIF
524      !
525   END SUBROUTINE trc_rst_wri 
526
527
528   SUBROUTINE trc_rst_cal( kt, cdrw )
529      !!---------------------------------------------------------------------
530      !!                   ***  ROUTINE trc_rst_cal  ***
531      !!
532      !!  ** Purpose : Read or write calendar in restart file:
533      !!
534      !!  WRITE(READ) mode:
535      !!       kt        : number of time step since the begining of the experiment at the
536      !!                   end of the current(previous) run
537      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
538      !!                   end of the current(previous) run (REAL -> keep fractions of day)
539      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
540      !!
541      !!   According to namelist parameter nrstdt,
542      !!       nn_rsttr = 0  no control on the date (nittrc000 is  arbitrary).
543      !!       nn_rsttr = 1  we verify that nittrc000 is equal to the last
544      !!                   time step of previous run + 1.
545      !!       In both those options, the  exact duration of the experiment
546      !!       since the beginning (cumulated duration of all previous restart runs)
547      !!       is not stored in the restart and is assumed to be (nittrc000-1)*rdt.
548      !!       This is valid is the time step has remained constant.
549      !!
550      !!       nn_rsttr = 2  the duration of the experiment in days (adatrj)
551      !!                    has been stored in the restart file.
552      !!----------------------------------------------------------------------
553      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
554      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
555      !
556      INTEGER  ::  jlibalt = jprstlib
557      LOGICAL  ::  llok
558      REAL(wp) ::  zkt, zrdttrc1
559      REAL(wp) ::  zndastp
560
561      ! Time domain : restart
562      ! ---------------------
563
564      IF( TRIM(cdrw) == 'READ' ) THEN
565
566         IF(lwp) WRITE(numout,*)
567         IF(lwp) WRITE(numout,*) 'trc_rst_cal : read the TOP restart file for calendar'
568         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
569
570         IF ( jprstlib == jprstdimg ) THEN
571           ! eventually read netcdf file (monobloc)  for restarting on different number of processors
572           ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90
573           INQUIRE( FILE = TRIM(cn_trcrst_indir)//'/'//TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
574           IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
575         ENDIF
576
577         IF( ln_rsttr ) THEN
578            CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
579            CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
580
581            IF(lwp) THEN
582               WRITE(numout,*) ' *** Info read in restart : '
583               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
584               WRITE(numout,*) ' *** restart option'
585               SELECT CASE ( nn_rsttr )
586               CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
587               CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
588               CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
589               END SELECT
590               WRITE(numout,*)
591            ENDIF
592            ! Control of date
593            IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
594               &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
595               &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
596         ENDIF
597         !
598         IF( lk_offline ) THEN   
599            !                                          ! set the date in offline mode
600            IF( ln_rsttr .AND. nn_rsttr == 2 ) THEN
601               CALL iom_get( numrtr, 'ndastp', zndastp ) 
602               ndastp = NINT( zndastp )
603               CALL iom_get( numrtr, 'adatrj', adatrj  )
604             ELSE
605               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
606               adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
607               ! note this is wrong if time step has changed during run
608            ENDIF
609            !
610            IF(lwp) THEN
611              WRITE(numout,*) ' *** Info used values : '
612              WRITE(numout,*) '   date ndastp                                      : ', ndastp
613              WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
614              WRITE(numout,*)
615            ENDIF
616            !
617            IF( ln_rsttr )  THEN   ;    neuler = 1
618            ELSE                   ;    neuler = 0
619            ENDIF
620            !
621            CALL day_init          ! compute calendar
622            !
623         ENDIF
624         !
625      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
626         !
627         IF(  kt == nitrst ) THEN
628            IF(lwp) WRITE(numout,*)
629            IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
630            IF(lwp) WRITE(numout,*) '~~~~~~~'
631         ENDIF
632         CALL iom_rstput( kt, nitrst, numrtw, 'kt'     , REAL( kt    , wp) )   ! time-step
633         CALL iom_rstput( kt, nitrst, numrtw, 'ndastp' , REAL( ndastp, wp) )   ! date
634         CALL iom_rstput( kt, nitrst, numrtw, 'adatrj' , adatrj            )   ! number of elapsed days since
635         !                                                                     ! the begining of the run [s]
636      ENDIF
637
638   END SUBROUTINE trc_rst_cal
639
640
641   SUBROUTINE trc_rst_stat
642      !!----------------------------------------------------------------------
643      !!                    ***  trc_rst_stat  ***
644      !!
645      !! ** purpose  :   Compute tracers statistics
646      !!----------------------------------------------------------------------
647      INTEGER  :: jk, jn
648      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
649      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
650      !!----------------------------------------------------------------------
651
652      IF( lwp ) THEN
653         WRITE(numout,*) 
654         WRITE(numout,*) '           ----TRACER STAT----             '
655         WRITE(numout,*) 
656      ENDIF
657      !
658      DO jk = 1, jpk
659         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
660      END DO
661      !
662      DO jn = 1, jptra
663         ztraf = glob_sum( trn(:,:,:,jn) * zvol(:,:,:) )
664         zmin  = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
665         zmax  = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
666         IF( lk_mpp ) THEN
667            CALL mpp_min( zmin )      ! min over the global domain
668            CALL mpp_max( zmax )      ! max over the global domain
669         END IF
670         zmean  = ztraf / areatot
671         zdrift = ( ( ztraf - trai(jn) ) / ( trai(jn) + 1.e-12 )  ) * 100._wp
672         IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
673      END DO
674      IF(lwp) WRITE(numout,*) 
6759000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
676      &      '    max :',e18.10,'    drift :',e18.10, ' %')
677      !
678   END SUBROUTINE trc_rst_stat
679
680
681   SUBROUTINE trc_rst_tra_stat
682      !!----------------------------------------------------------------------
683      !!                    ***  trc_rst_tra_stat  ***
684      !!
685      !! ** purpose  :   Compute tracers statistics - check where crazy values appears
686      !!----------------------------------------------------------------------
687      INTEGER  :: jk, jn
688      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift, areasf
689      REAL(wp), DIMENSION(jpi,jpj) :: zvol
690      !!----------------------------------------------------------------------
691
692      IF( lwp ) THEN
693         WRITE(numout,*)
694         WRITE(numout,*) '           ----SURFACE TRA STAT----             '
695         WRITE(numout,*)
696      ENDIF
697      !
698      zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
699      areasf = glob_sum(zvol(:,:))
700      DO jn = 1, jptra
701         ztraf = glob_sum( tra(:,:,1,jn) * zvol(:,:) )
702         zmin  = MINVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) )
703         zmax  = MAXVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) )
704         IF( lk_mpp ) THEN
705            CALL mpp_min( zmin )      ! min over the global domain
706            CALL mpp_max( zmax )      ! max over the global domain
707         END IF
708         zmean  = ztraf / areasf
709         IF(lwp) WRITE(numout,9001) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax
710      END DO
711      IF(lwp) WRITE(numout,*)
7129001  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
713      &      '    max :',e18.10)
714      !
715   END SUBROUTINE trc_rst_tra_stat
716
717
718
719   SUBROUTINE trc_rst_dia_stat( dgtr, names)
720      !!----------------------------------------------------------------------
721      !!                    ***  trc_rst_dia_stat  ***
722      !!
723      !! ** purpose  :   Compute tracers statistics
724      !!----------------------------------------------------------------------
725      REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) ::   dgtr      ! 2D diag var
726      CHARACTER(len=*)             , INTENT(in) ::   names     ! 2D diag name
727      !!---------------------------------------------------------------------
728      INTEGER  :: jk, jn
729      REAL(wp) :: ztraf, zmin, zmax, zmean, areasf
730      REAL(wp), DIMENSION(jpi,jpj) :: zvol
731      !!----------------------------------------------------------------------
732
733      IF( lwp )  WRITE(numout,*) 'STAT- ', names
734      !
735      zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
736      ztraf = glob_sum( dgtr(:,:) * zvol(:,:) )
737      !! areasf = glob_sum(e1e2t(:,:) * tmask(:,:,1) )
738      areasf = glob_sum(zvol(:,:))
739      zmin  = MINVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
740      zmax  = MAXVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
741      IF( lk_mpp ) THEN
742         CALL mpp_min( zmin )      ! min over the global domain
743         CALL mpp_max( zmax )      ! max over the global domain
744      END IF
745      zmean  = ztraf / areasf
746      IF(lwp) WRITE(numout,9002) TRIM( names ), zmean, zmin, zmax
747      !
748      IF(lwp) WRITE(numout,*)
7499002  FORMAT(' tracer name :',a10,'    mean :',e18.10,'    min :',e18.10, &
750      &      '    max :',e18.10 )
751      !
752   END SUBROUTINE trc_rst_dia_stat
753
754
755#else
756   !!----------------------------------------------------------------------
757   !!  Dummy module :                                     No passive tracer
758   !!----------------------------------------------------------------------
759CONTAINS
760   SUBROUTINE trc_rst_read                      ! Empty routines
761   END SUBROUTINE trc_rst_read
762   SUBROUTINE trc_rst_wri( kt )
763      INTEGER, INTENT ( in ) :: kt
764      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
765   END SUBROUTINE trc_rst_wri   
766#endif
767
768   !!----------------------------------------------------------------------
769   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
770   !! $Id$
771   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
772   !!======================================================================
773END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.