New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcrst.F90 in branches/NERC/dev_r5518_GO6_ScalingCoupledChl/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/NERC/dev_r5518_GO6_ScalingCoupledChl/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 8955

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

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