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

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

Last change on this file since 9114 was 9114, checked in by frrh, 5 years ago

Apply changes developed under Met Office GMED ticket number 351 in development
branch branches/NERC/dev_r5518_GO6_ScalingCoupledChl.

The command issued to perform the merge is:

svn merge -r 8590:9053 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/NERC/dev_r5518_GO6_ScalingCoupledChl

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