source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 @ 5707

Last change on this file since 5707 was 5707, checked in by acc, 5 years ago

JPALM —25-08-2015 — add MEDUSA in the branch. MEDUSA version already up-to-date with this trunk revision

  • Property svn:keywords set to Id
File size: 19.8 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 daymod
30   !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
31   USE sms_medusa
32   USE trcsms_medusa
33   !!
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   trc_rst_opn       ! called by ???
38   PUBLIC   trc_rst_read      ! called by ???
39   PUBLIC   trc_rst_wri       ! called by ???
40   PUBLIC   trc_rst_cal
41
42   !! * Substitutions
43#  include "top_substitute.h90"
44   
45CONTAINS
46   
47   SUBROUTINE trc_rst_opn( kt )
48      !!----------------------------------------------------------------------
49      !!                    ***  trc_rst_opn  ***
50      !!
51      !! ** purpose  :   output of sea-trc variable in a netcdf file
52      !!----------------------------------------------------------------------
53      INTEGER, INTENT(in) ::   kt       ! number of iteration
54      !
55      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
56      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
57      !!----------------------------------------------------------------------
58      !
59      IF( lk_offline ) THEN
60         IF( kt == nittrc000 ) THEN
61            lrst_trc = .FALSE.
62            nitrst = nitend
63         ENDIF
64
65         IF( MOD( kt - 1, nstock ) == 0 ) THEN
66            ! we use kt - 1 and not kt - nittrc000 to keep the same periodicity from the beginning of the experiment
67            nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
68            IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
69         ENDIF
70      ELSE
71         IF( kt == nittrc000 ) lrst_trc = .FALSE.
72      ENDIF
73
74      ! to get better performances with NetCDF format:
75      ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1)
76      ! 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
77      IF( kt == nitrst - 2*nn_dttrc .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
78         ! beware of the format used to write kt (default is i8.8, that should be large enough)
79         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
80         ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
81         ENDIF
82         ! create the file
83         IF(lwp) WRITE(numout,*)
84         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trcrst_out)
85         IF(lwp) WRITE(numout,*) '             open trc restart.output NetCDF file: '//clname
86         CALL iom_open( clname, numrtw, ldwrt = .TRUE., kiolib = jprstlib )
87         lrst_trc = .TRUE.
88      ENDIF
89      !
90   END SUBROUTINE trc_rst_opn
91
92   SUBROUTINE trc_rst_read
93      !!----------------------------------------------------------------------
94      !!                    ***  trc_read_opn  ***
95      !!
96      !! ** purpose  :   read passive tracer fields in restart files
97      !!----------------------------------------------------------------------
98      INTEGER  ::  jn     
99      !! AXY (05/11/13): temporary variables
100      REAL(wp) ::    fq0,fq1,fq2
101
102      !!----------------------------------------------------------------------
103      !
104      IF(lwp) WRITE(numout,*)
105      IF(lwp) WRITE(numout,*) 'trc_rst_read : read data in the TOP restart file'
106      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
107
108      ! READ prognostic variables and computes diagnostic variable
109      DO jn = 1, jptra
110         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
111      END DO
112
113      DO jn = 1, jptra
114         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
115      END DO
116
117      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
118      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
119      !!                 version of NEMO date significantly earlier than the current
120      !!                 version
121
122#if defined key_medusa
123      !! AXY (13/01/12): check if the restart contains sediment fields;
124      !!                 this is only relevant for simulations that include
125      !!                 biogeochemistry and are restarted from earlier runs
126      !!                 in which there was no sediment component
127      !!
128      IF( iom_varid( numrtr, 'B_SED_N', ldstop = .FALSE. ) > 0 ) THEN
129         !! YES; in which case read them
130         !!
131         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields present - reading in ...'
132         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_N',  zb_sed_n(:,:)  )
133         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_N',  zn_sed_n(:,:)  )
134         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_FE', zb_sed_fe(:,:) )
135         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_FE', zn_sed_fe(:,:) )
136         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_SI', zb_sed_si(:,:) )
137         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_SI', zn_sed_si(:,:) )
138         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_C',  zb_sed_c(:,:)  )
139         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_C',  zn_sed_c(:,:)  )
140         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_CA', zb_sed_ca(:,:) )
141         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_CA', zn_sed_ca(:,:) )
142      ELSE
143         !! NO; in which case set them to zero
144         !!
145         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields absent - setting to zero ...'
146         zb_sed_n(:,:)  = 0.0   !! organic N
147         zn_sed_n(:,:)  = 0.0
148         zb_sed_fe(:,:) = 0.0   !! organic Fe
149         zn_sed_fe(:,:) = 0.0
150         zb_sed_si(:,:) = 0.0   !! inorganic Si
151         zn_sed_si(:,:) = 0.0
152         zb_sed_c(:,:)  = 0.0   !! organic C
153         zn_sed_c(:,:)  = 0.0
154         zb_sed_ca(:,:) = 0.0   !! inorganic C
155         zn_sed_ca(:,:) = 0.0
156      ENDIF
157      !!
158      !! calculate stats on these fields
159      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
160      fq0 = MINVAL(zn_sed_n(:,:))
161      fq1 = MAXVAL(zn_sed_n(:,:))
162      fq2 = SUM(zn_sed_n(:,:))
163      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
164         &        fq0, fq1, fq2
165      fq0 = MINVAL(zn_sed_fe(:,:))
166      fq1 = MAXVAL(zn_sed_fe(:,:))
167      fq2 = SUM(zn_sed_fe(:,:))
168      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
169         &        fq0, fq1, fq2
170      fq0 = MINVAL(zn_sed_si(:,:))
171      fq1 = MAXVAL(zn_sed_si(:,:))
172      fq2 = SUM(zn_sed_si(:,:))
173      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
174         &        fq0, fq1, fq2
175      fq0 = MINVAL(zn_sed_c(:,:))
176      fq1 = MAXVAL(zn_sed_c(:,:))
177      fq2 = SUM(zn_sed_c(:,:))
178      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
179         &        fq0, fq1, fq2
180      fq0 = MINVAL(zn_sed_ca(:,:))
181      fq1 = MAXVAL(zn_sed_ca(:,:))
182      fq2 = SUM(zn_sed_ca(:,:))
183      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
184         &        fq0, fq1, fq2
185#endif
186 
187      !
188   END SUBROUTINE trc_rst_read
189
190   SUBROUTINE trc_rst_wri( kt )
191      !!----------------------------------------------------------------------
192      !!                    ***  trc_rst_wri  ***
193      !!
194      !! ** purpose  :   write passive tracer fields in restart files
195      !!----------------------------------------------------------------------
196      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
197      !!
198      INTEGER  :: jn
199      REAL(wp) :: zarak0
200      !! AXY (05/11/13): temporary variables
201      REAL(wp) ::    fq0,fq1,fq2
202      !!----------------------------------------------------------------------
203      !
204      CALL iom_rstput( kt, nitrst, numrtw, 'rdttrc1', rdttrc(1) )   ! surface passive tracer time step
205      ! prognostic variables
206      ! --------------------
207      DO jn = 1, jptra
208         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
209      END DO
210
211      DO jn = 1, jptra
212         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
213      END DO
214
215      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
216      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
217      !!                 version of NEMO date significantly earlier than the current
218      !!                 version
219
220#if defined key_medusa
221      !! AXY (13/01/12): write out "before" and "now" state of seafloor
222      !!                 sediment pools into restart; this happens
223      !!                 whether or not the pools are to be used by
224      !!                 MEDUSA (which is controlled by a switch in the
225      !!                 namelist_top file)
226      !!
227      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
228      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_N',  zb_sed_n(:,:)  )
229      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_N',  zn_sed_n(:,:)  )
230      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_FE', zb_sed_fe(:,:) )
231      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_FE', zn_sed_fe(:,:) )
232      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_SI', zb_sed_si(:,:) )
233      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_SI', zn_sed_si(:,:) )
234      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_C',  zb_sed_c(:,:)  )
235      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_C',  zn_sed_c(:,:)  )
236      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_CA', zb_sed_ca(:,:) )
237      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_CA', zn_sed_ca(:,:) )
238      !!
239      !! calculate stats on these fields
240      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
241      fq0 = MINVAL(zn_sed_n(:,:))
242      fq1 = MAXVAL(zn_sed_n(:,:))
243      fq2 = SUM(zn_sed_n(:,:))
244      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
245         &        fq0, fq1, fq2
246      fq0 = MINVAL(zn_sed_fe(:,:))
247      fq1 = MAXVAL(zn_sed_fe(:,:))
248      fq2 = SUM(zn_sed_fe(:,:))
249      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
250         &        fq0, fq1, fq2
251      fq0 = MINVAL(zn_sed_si(:,:))
252      fq1 = MAXVAL(zn_sed_si(:,:))
253      fq2 = SUM(zn_sed_si(:,:))
254      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
255         &        fq0, fq1, fq2
256      fq0 = MINVAL(zn_sed_c(:,:))
257      fq1 = MAXVAL(zn_sed_c(:,:))
258      fq2 = SUM(zn_sed_c(:,:))
259      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
260         &        fq0, fq1, fq2
261      fq0 = MINVAL(zn_sed_ca(:,:))
262      fq1 = MAXVAL(zn_sed_ca(:,:))
263      fq2 = SUM(zn_sed_ca(:,:))
264      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
265         &        fq0, fq1, fq2
266#endif
267
268      !
269      IF( kt == nitrst ) THEN
270          CALL trc_rst_stat            ! statistics
271          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
272#if ! defined key_trdmxl_trc
273          lrst_trc = .FALSE.
274#endif
275      ENDIF
276      !
277   END SUBROUTINE trc_rst_wri 
278
279
280   SUBROUTINE trc_rst_cal( kt, cdrw )
281      !!---------------------------------------------------------------------
282      !!                   ***  ROUTINE trc_rst_cal  ***
283      !!
284      !!  ** Purpose : Read or write calendar in restart file:
285      !!
286      !!  WRITE(READ) mode:
287      !!       kt        : number of time step since the begining of the experiment at the
288      !!                   end of the current(previous) run
289      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
290      !!                   end of the current(previous) run (REAL -> keep fractions of day)
291      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
292      !!
293      !!   According to namelist parameter nrstdt,
294      !!       nn_rsttr = 0  no control on the date (nittrc000 is  arbitrary).
295      !!       nn_rsttr = 1  we verify that nittrc000 is equal to the last
296      !!                   time step of previous run + 1.
297      !!       In both those options, the  exact duration of the experiment
298      !!       since the beginning (cumulated duration of all previous restart runs)
299      !!       is not stored in the restart and is assumed to be (nittrc000-1)*rdt.
300      !!       This is valid is the time step has remained constant.
301      !!
302      !!       nn_rsttr = 2  the duration of the experiment in days (adatrj)
303      !!                    has been stored in the restart file.
304      !!----------------------------------------------------------------------
305      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
306      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
307      !
308      INTEGER  ::  jlibalt = jprstlib
309      LOGICAL  ::  llok
310      REAL(wp) ::  zkt, zrdttrc1
311      REAL(wp) ::  zndastp
312
313      ! Time domain : restart
314      ! ---------------------
315
316      IF( TRIM(cdrw) == 'READ' ) THEN
317
318         IF(lwp) WRITE(numout,*)
319         IF(lwp) WRITE(numout,*) 'trc_rst_cal : read the TOP restart file for calendar'
320         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
321
322         IF ( jprstlib == jprstdimg ) THEN
323           ! eventually read netcdf file (monobloc)  for restarting on different number of processors
324           ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90
325           INQUIRE( FILE = TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
326           IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
327         ENDIF
328
329         CALL iom_open( cn_trcrst_in, numrtr, kiolib = jlibalt )
330
331         CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
332
333         IF(lwp) THEN
334            WRITE(numout,*) ' *** Info read in restart : '
335            WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
336            WRITE(numout,*) ' *** restart option'
337            SELECT CASE ( nn_rsttr )
338            CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
339            CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
340            CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
341            END SELECT
342            WRITE(numout,*)
343         ENDIF
344         ! Control of date
345         IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
346            &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
347            &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
348         IF( lk_offline ) THEN      ! set the date in offline mode
349            ! Check dynamics and tracer time-step consistency and force Euler restart if changed
350            IF( iom_varid( numrtr, 'rdttrc1', ldstop = .FALSE. ) > 0 )   THEN
351               CALL iom_get( numrtr, 'rdttrc1', zrdttrc1 )
352               IF( zrdttrc1 /= rdt * nn_dttrc )   neuler = 0
353            ENDIF
354            !                          ! define ndastp and adatrj
355            IF( nn_rsttr == 2 ) THEN
356               CALL iom_get( numrtr, 'ndastp', zndastp ) 
357               ndastp = NINT( zndastp )
358               CALL iom_get( numrtr, 'adatrj', adatrj  )
359            ELSE
360               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
361               adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
362               ! note this is wrong if time step has changed during run
363            ENDIF
364            !
365            IF(lwp) THEN
366              WRITE(numout,*) ' *** Info used values : '
367              WRITE(numout,*) '   date ndastp                                      : ', ndastp
368              WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
369              WRITE(numout,*)
370            ENDIF
371            !
372            CALL day_init          ! compute calendar
373            !
374         ENDIF
375         !
376      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
377         !
378         IF(  kt == nitrst ) THEN
379            IF(lwp) WRITE(numout,*)
380            IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
381            IF(lwp) WRITE(numout,*) '~~~~~~~'
382         ENDIF
383         CALL iom_rstput( kt, nitrst, numrtw, 'kt'     , REAL( kt    , wp) )   ! time-step
384         CALL iom_rstput( kt, nitrst, numrtw, 'ndastp' , REAL( ndastp, wp) )   ! date
385         CALL iom_rstput( kt, nitrst, numrtw, 'adatrj' , adatrj            )   ! number of elapsed days since
386         !                                                                     ! the begining of the run [s]
387      ENDIF
388
389   END SUBROUTINE trc_rst_cal
390
391
392   SUBROUTINE trc_rst_stat
393      !!----------------------------------------------------------------------
394      !!                    ***  trc_rst_stat  ***
395      !!
396      !! ** purpose  :   Compute tracers statistics
397      !!----------------------------------------------------------------------
398      INTEGER  :: jk, jn
399      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
400      !!----------------------------------------------------------------------
401
402      IF( lwp ) THEN
403         WRITE(numout,*) 
404         WRITE(numout,*) '           ----TRACER STAT----             '
405         WRITE(numout,*) 
406      ENDIF
407      !
408      DO jn = 1, jptra
409         ztraf = glob_sum( trn(:,:,:,jn) * cvol(:,:,:) )
410         zmin  = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
411         zmax  = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
412         IF( lk_mpp ) THEN
413            CALL mpp_min( zmin )      ! min over the global domain
414            CALL mpp_max( zmax )      ! max over the global domain
415         END IF
416         zmean  = ztraf / areatot
417         zdrift = ( ( ztraf - trai(jn) ) / ( trai(jn) + 1.e-12 )  ) * 100._wp
418         IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
419      END DO
420      WRITE(numout,*) 
4219000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
422      &      '    max :',e18.10,'    drift :',e18.10, ' %')
423      !
424   END SUBROUTINE trc_rst_stat
425
426#else
427   !!----------------------------------------------------------------------
428   !!  Dummy module :                                     No passive tracer
429   !!----------------------------------------------------------------------
430CONTAINS
431   SUBROUTINE trc_rst_read                      ! Empty routines
432   END SUBROUTINE trc_rst_read
433   SUBROUTINE trc_rst_wri( kt )
434      INTEGER, INTENT ( in ) :: kt
435      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
436   END SUBROUTINE trc_rst_wri   
437#endif
438
439   !!----------------------------------------------------------------------
440   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
441   !! $Id$
442   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
443   !!======================================================================
444END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.