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.
restart.F90 in NEMO/trunk/src/OCE/IOM – NEMO

source: NEMO/trunk/src/OCE/IOM/restart.F90 @ 14179

Last change on this file since 14179 was 14179, checked in by techene, 3 years ago

#2385 restart with qco bug fix for euler 1st and RK3

  • Property svn:keywords set to Id
File size: 20.5 KB
Line 
1MODULE restart
2   !!======================================================================
3   !!                     ***  MODULE  restart  ***
4   !! Ocean restart :  write the ocean restart file
5   !!======================================================================
6   !! History :  OPA  !  1999-11  (M. Imbard)  Original code
7   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form
8   !!            2.0  !  2006-07  (S. Masson)  use IOM for restart
9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
10   !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
11   !!            3.7  !  2014-01  (G. Madec) suppression of curl and hdiv from the restart
12   !!             -   !  2014-12  (G. Madec) remove KPP scheme
13   !!            4.1  !  2020-11  (S. Techene, G. Madec)  move ssh initiatlisation in rst_read_ssh
14   !!             -   !                                   add restart in Shallow Water Eq. case
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   rst_opn       : open the ocean restart file for writting
19   !!   rst_write     : write the ocean restart file
20   !!   rst_read_open : open the restart file for reading
21   !!   rst_read      : read the ocean restart file
22   !!   rst_read_ssh  : ssh set from restart or domcfg.nc file or usr_def_istat_ssh
23   !!----------------------------------------------------------------------
24   USE oce            ! ocean dynamics and tracers
25   USE dom_oce        ! ocean space and time domain
26   USE sbc_ice        ! only lk_si3
27   USE phycst         ! physical constants
28   USE eosbn2         ! equation of state
29   USE wet_dry        ! Wetting/Drying flux limiting
30   USE usrdef_istate, ONLY : usr_def_istate_ssh   ! user defined ssh initial state
31   USE trdmxl_oce     ! ocean active mixed layer tracers trends variables
32   USE diu_bulk       ! ???
33   !
34   USE in_out_manager ! I/O manager
35   USE iom            ! I/O module
36   USE lib_mpp        ! distribued memory computing library
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   rst_opn         ! called by step.F90
42   PUBLIC   rst_write       ! called by step.F90
43   PUBLIC   rst_read_open   ! called in rst_read_ssh
44   PUBLIC   rst_read        ! called by istate.F90
45   PUBLIC   rst_read_ssh    ! called by domain.F90
46   
47   !! * Substitutions
48#  include "do_loop_substitute.h90"
49#  include "domzgr_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE rst_opn( kt )
58      !!---------------------------------------------------------------------
59      !!                   ***  ROUTINE rst_opn  ***
60      !!
61      !! ** Purpose : + initialization (should be read in the namelist) of nitrst
62      !!              + open the restart when we are one time step before nitrst
63      !!                   - restart header is defined when kt = nitrst-1
64      !!                   - restart data  are written when kt = nitrst
65      !!              + define lrst_oce to .TRUE. when we need to define or write the restart
66      !!----------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   kt     ! ocean time-step
68      !!
69      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
70      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
71      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file
72      CHARACTER(LEN=52)   ::   clpname   ! ocean output restart file name including prefix for AGRIF
73      CHARACTER(LEN=256)  ::   clinfo    ! info character
74      !!----------------------------------------------------------------------
75      !
76      IF( kt == nit000 ) THEN   ! default definitions
77         lrst_oce = .FALSE.
78         IF( ln_rst_list ) THEN
79            nrst_lst = 1
80            nitrst = nn_stocklist( nrst_lst )
81         ELSE
82            nitrst = nitend
83         ENDIF
84      ENDIF
85
86      IF( .NOT. ln_rst_list .AND. nn_stock == -1 )   RETURN   ! we will never do any restart
87
88      ! frequency-based restart dumping (nn_stock)
89      IF( .NOT. ln_rst_list .AND. MOD( kt - 1, nn_stock ) == 0 ) THEN
90         ! we use kt - 1 and not kt - nit000 to keep the same periodicity from the beginning of the experiment
91         nitrst = kt + nn_stock - 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      ! to get better performances with NetCDF format:
95      ! we open and define the ocean restart file one time step before writing the data (-> at nitrst - 1)
96      ! except if we write ocean restart files every time step or if an ocean restart file was writen at nitend - 1
97      IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
98         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN
99            ! beware of the format used to write kt (default is i8.8, that should be large enough...)
100            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
101            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
102            ENDIF
103            ! create the file
104            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
105            clpath = TRIM(cn_ocerst_outdir)
106            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
107            IF(lwp) THEN
108               WRITE(numout,*)
109               IF(.NOT.lwxios) THEN
110                  WRITE(numout,*) '             open ocean restart NetCDF file: ',TRIM(clpath)//TRIM(clname)
111                  IF ( snc4set%luse )      WRITE(numout,*) '             opened for NetCDF4 chunking and compression'
112                  IF( kt == nitrst - 1 ) THEN   ;   WRITE(numout,*) '             kt = nitrst - 1 = ', kt
113                  ELSE                          ;   WRITE(numout,*) '             kt = '             , kt
114                  ENDIF
115               ENDIF
116            ENDIF
117            !
118            IF(.NOT.lwxios) THEN
119               CALL iom_open( TRIM(clpath)//TRIM(clname), numrow, ldwrt = .TRUE. )
120            ELSE
121#if defined key_iomput
122               cw_ocerst_cxt = "rstw_"//TRIM(ADJUSTL(clkt))
123               IF( TRIM(Agrif_CFixed()) == '0' ) THEN
124                  clpname = clname
125               ELSE
126                  clpname = TRIM(Agrif_CFixed())//"_"//clname
127               ENDIF
128               numrow = iom_xios_setid(TRIM(clpath)//TRIM(clpname))
129               CALL iom_init( cw_ocerst_cxt, kdid = numrow, ld_closedef = .false. )
130               CALL iom_swap(      cxios_context          )
131#else
132               clinfo = 'Can not use XIOS in rst_opn'
133               CALL ctl_stop(TRIM(clinfo))
134#endif
135            ENDIF
136            lrst_oce = .TRUE.
137         ENDIF
138      ENDIF
139      !
140   END SUBROUTINE rst_opn
141
142
143   SUBROUTINE rst_write( kt, Kbb, Kmm )
144      !!---------------------------------------------------------------------
145      !!                   ***  ROUTINE rstwrite  ***
146      !!
147      !! ** Purpose :   Write restart fields in NetCDF format
148      !!
149      !! ** Method  :   Write in numrow when kt == nitrst in NetCDF
150      !!              file, save fields which are necessary for restart
151      !!
152      !!                NB: ssh is written here (rst_write)
153      !!                    but is read or set in rst_read_ssh
154      !!----------------------------------------------------------------------
155      INTEGER, INTENT(in) ::   kt         ! ocean time-step
156      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
157      !!----------------------------------------------------------------------
158      !
159         CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rn_Dt       )   ! dynamics time step
160      !
161      IF( .NOT.lwxios )   CALL iom_delay_rst( 'WRITE', 'OCE', numrow )   ! save only ocean delayed global communication variables
162      !
163      IF( .NOT.ln_diurnal_only ) THEN
164         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , ssh(:,:        ,Kbb) )     ! before fields
165         CALL iom_rstput( kt, nitrst, numrow, 'ub'     , uu(:,:,:       ,Kbb) )
166         CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vv(:,:,:       ,Kbb) )
167         IF( .NOT.lk_SWE ) THEN
168            CALL iom_rstput( kt, nitrst, numrow, 'tb'  , ts(:,:,:,jp_tem,Kbb) )
169            CALL iom_rstput( kt, nitrst, numrow, 'sb'  , ts(:,:,:,jp_sal,Kbb) )
170         ENDIF
171         !
172#if ! defined key_RK3
173         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , ssh(:,:        ,Kmm) )     ! now fields
174         CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:       ,Kmm) )
175         CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vv(:,:,:       ,Kmm) )
176         IF( .NOT.lk_SWE ) THEN
177            CALL iom_rstput( kt, nitrst, numrow, 'tn'  , ts(:,:,:,jp_tem,Kmm) )
178            CALL iom_rstput( kt, nitrst, numrow, 'sn'  , ts(:,:,:,jp_sal,Kmm) )
179            CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop                 )
180         ENDIF
181#endif
182      ENDIF
183
184      IF( ln_diurnal )   CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst )
185      IF( kt == nitrst ) THEN
186         IF( .NOT.lwxios ) THEN
187            CALL iom_close( numrow )     ! close the restart file (only at last time step)
188         ELSE
189            CALL iom_context_finalize(      cw_ocerst_cxt          )
190            iom_file(numrow)%nfid       = 0
191            numrow = 0
192         ENDIF
193!!gm         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
194!!gm  not sure what to do here   ===>>>  ask to Sebastian
195         lrst_oce = .FALSE.
196         IF( ln_rst_list ) THEN
197            nrst_lst = MIN(nrst_lst + 1, SIZE(nn_stocklist,1))
198            nitrst   = nn_stocklist( nrst_lst )
199         ENDIF
200      ENDIF
201      !
202   END SUBROUTINE rst_write
203
204
205   SUBROUTINE rst_read_open
206      !!----------------------------------------------------------------------
207      !!                   ***  ROUTINE rst_read_open  ***
208      !!
209      !! ** Purpose :   Open read files for NetCDF restart
210      !!
211      !! ** Method  :   Use a non-zero, positive value of numror to assess whether or not
212      !!                the file has already been opened
213      !!----------------------------------------------------------------------
214      LOGICAL             ::   llok
215      CHARACTER(len=lc)   ::   clpath   ! full path to ocean output restart file
216      CHARACTER(len=lc+2) ::   clpname  ! file name including agrif prefix
217      !!----------------------------------------------------------------------
218      !
219      IF( numror <= 0 ) THEN
220         IF(lwp) THEN                                             ! Contol prints
221            WRITE(numout,*)
222            WRITE(numout,*) 'rst_read : read oce NetCDF restart file'
223            IF ( snc4set%luse )      WRITE(numout,*) 'rst_read : configured with NetCDF4 support'
224            WRITE(numout,*) '~~~~~~~~'
225         ENDIF
226         lxios_sini = .FALSE.
227         clpath = TRIM(cn_ocerst_indir)
228         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
229         CALL iom_open( TRIM(clpath)//cn_ocerst_in, numror )
230! are we using XIOS to read the data? Part above will have to modified once XIOS
231! can handle checking if variable is in the restart file (there will be no need to open
232! restart)
233         lrxios = lrxios.AND.lxios_sini
234
235         IF( lrxios) THEN
236             cr_ocerst_cxt = 'oce_rst'
237             IF(lwp) WRITE(numout,*) 'Enable restart reading by XIOS'
238!            IF( TRIM(Agrif_CFixed()) == '0' ) THEN
239!               clpname = cn_ocerst_in
240!            ELSE
241!               clpname = TRIM(Agrif_CFixed())//"_"//cn_ocerst_in
242!            ENDIF
243             CALL iom_init( cr_ocerst_cxt, kdid = numror, ld_closedef = .TRUE. )
244             CALL iom_swap(      cxios_context          )
245         ENDIF
246
247      ENDIF
248
249   END SUBROUTINE rst_read_open
250
251
252   SUBROUTINE rst_read( Kbb, Kmm )
253      !!----------------------------------------------------------------------
254      !!                   ***  ROUTINE rst_read  ***
255      !!
256      !! ** Purpose :   Read velocity and T-S fields in the restart file
257      !!
258      !! ** Method  :   Read in restart.nc fields which are necessary for restart
259      !!
260      !!                NB: restart file openned           in DOM/domain.F90:dom_init
261      !!                    before field in restart tested in DOM/domain.F90:dom_init
262      !!                    (sshb)
263      !!
264      !!                NB: ssh is read or set in rst_read_ssh
265      !!----------------------------------------------------------------------
266      INTEGER          , INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
267      INTEGER  ::   jk
268      REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
269      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zgdept       ! 3D workspace for QCO
270      !!----------------------------------------------------------------------
271      !
272      IF(.NOT.lrxios )   CALL iom_delay_rst( 'READ', 'OCE', numror )   ! read only ocean delayed global communication variables
273      !
274      !                             !*  Diurnal DSST
275      IF( ln_diurnal )   CALL iom_get( numror, jpdom_auto, 'Dsst' , x_dsst )
276      IF( ln_diurnal_only ) THEN
277         IF(lwp) WRITE( numout, * ) &
278         &   "rst_read:- ln_diurnal_only set, setting rhop=rho0"
279         rhop = rho0
280         CALL iom_get( numror, jpdom_auto, 'tn'     , w3d )
281         ts(:,:,1,jp_tem,Kmm) = w3d(:,:,1)
282         RETURN
283      ENDIF
284      !
285#if defined key_RK3
286      !                             !*  Read Kbb fields   (NB: in RK3 Kmm = Kbb = Nbb)
287      IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
288      CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
289      CALL iom_get( numror, jpdom_auto, 'vb'   , vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
290      IF( .NOT.lk_SWE ) THEN
291         CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
292         CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
293      ENDIF
294#else
295      !                             !*  Read Kmm fields   (MLF only)
296      IF(lwp) WRITE(numout,*)    '           Kmm u, v and T-S fields read in the restart file'
297      CALL iom_get( numror, jpdom_auto, 'un'   , uu(:,:,:       ,Kmm), cd_type = 'U', psgn = -1._wp )
298      CALL iom_get( numror, jpdom_auto, 'vn'   , vv(:,:,:       ,Kmm), cd_type = 'V', psgn = -1._wp )
299      IF( .NOT.lk_SWE ) THEN
300         CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
301         CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
302      ENDIF
303      !
304      IF( l_1st_euler ) THEN        !*  Euler restart   (MLF only)
305         IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields set to Kmm values'
306         uu(:,:,:  ,Kbb) = uu(:,:,:  ,Kmm)         ! all before fields set to now values
307         vv(:,:,:  ,Kbb) = vv(:,:,:  ,Kmm)
308         IF( .NOT.lk_SWE ) ts(:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
309         !
310      ELSE                          !* Leap frog restart   (MLF only)
311         IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
312         CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
313         CALL iom_get( numror, jpdom_auto, 'vb'   , vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
314         IF( .NOT.lk_SWE ) THEN
315            CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
316            CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
317         ENDIF
318      ENDIF
319#endif
320      !
321      IF( .NOT.lk_SWE ) THEN
322         IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
323            CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density
324         ELSE
325#if defined key_qco
326            ALLOCATE( zgdept(jpi,jpj,jpk) )
327            DO jk = 1, jpk
328               zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
329            END DO
330            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, zgdept )
331            DEALLOCATE( zgdept )
332#else
333            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept(:,:,:,Kmm) )
334#endif
335         ENDIF
336      ENDIF
337      !
338   END SUBROUTINE rst_read
339
340
341   SUBROUTINE rst_read_ssh( Kbb, Kmm, Kaa )
342      !!---------------------------------------------------------------------
343      !!                   ***  ROUTINE rst_read_ssh  ***
344      !!
345      !! ** Purpose :   ssh initialization of the sea surface height (ssh)
346      !!
347      !! ** Method  :   set ssh from restart or read configuration, or user_def
348      !!              * ln_rstart = T
349      !!                   USE of IOM library to read ssh in the restart file
350      !!                   Leap-Frog: Kbb and Kmm are read except for l_1st_euler=T
351      !!
352      !!              * otherwise
353      !!                   call user defined ssh or
354      !!                   set to -ssh_ref in wet and drying case with domcfg.nc
355      !!
356      !!              NB: ssh_b/n are written by restart.F90
357      !!----------------------------------------------------------------------
358      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time level indices
359      !
360      INTEGER ::   ji, jj, jk
361      !!----------------------------------------------------------------------
362      !
363      IF(lwp) THEN
364         WRITE(numout,*)
365         WRITE(numout,*) 'rst_read_ssh : ssh initialization'
366         WRITE(numout,*) '~~~~~~~~~~~~ '
367      ENDIF
368      !
369      !                            !=============================!
370      IF( ln_rstart ) THEN         !==  Read the restart file  ==!
371         !                         !=============================!
372         !
373#if defined key_RK3
374         !                                     !*  RK3: Read ssh at Kbb
375         IF(lwp) WRITE(numout,*)
376         IF(lwp) WRITE(numout,*)    '      Kbb sea surface height read in the restart file'
377         CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb) )
378         !
379         !                                     !*  RK3: Set ssh at Kmm for AGRIF
380         ssh(:,:,Kmm) = 0._wp
381#else
382         !                                     !*  MLF: Read ssh at Kmm
383         IF(lwp) WRITE(numout,*)
384         IF(lwp) WRITE(numout,*)    '      Kmm sea surface height read in the restart file'
385         CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm) )
386         !
387         IF( l_1st_euler ) THEN                !*  MLF: Euler at first time-step
388            IF(lwp) WRITE(numout,*)
389            IF(lwp) WRITE(numout,*) '      Euler first time step : ssh(Kbb) = ssh(Kmm)'
390            ssh(:,:,Kbb) = ssh(:,:,Kmm)
391            !
392         ELSE                                  !*  MLF: read ssh at Kbb
393            IF(lwp) WRITE(numout,*)
394            IF(lwp) WRITE(numout,*) '      Kbb sea surface height read in the restart file'
395            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) )
396         ENDIF
397#endif
398         !                         !============================!
399      ELSE                         !==  Initialize at "rest"  ==!
400         !                         !============================!
401         !
402         IF(lwp) WRITE(numout,*)
403         IF(lwp) WRITE(numout,*)    '      initialization at rest'
404         !
405         IF( ll_wd ) THEN                      !* wet and dry
406            !
407            IF( ln_read_cfg  ) THEN                 ! read configuration : ssh_ref is read in domain_cfg file
408!!st  why ssh is not masked : i.e. ssh(:,:,Kmm) = -ssh_ref*ssmask(:,:),
409!!st  since at the 1st time step lbclnk will be applied on ssh at Kaa but not initially at Kbb and Kmm
410               ssh(:,:,Kbb) = -ssh_ref
411               !
412               DO_2D( 1, 1, 1, 1 )
413                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN   ! if total depth is less than min depth
414                     ssh(ji,jj,Kbb) = rn_wdmin1 - ht_0(ji,jj)
415                  ENDIF
416               END_2D
417            ELSE                                    ! user define configuration case 
418               CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
419            ENDIF
420            !
421         ELSE                                  !* user defined configuration
422            CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
423            !
424         ENDIF
425         !
426#if defined key_RK3
427         ssh(:,:,Kmm) = 0._wp                  !* RK3: set Kmm to 0 for AGRIF
428#else
429         ssh(:,:,Kmm) = ssh(:,:,Kbb)           !* MLF: set now values from to before ones
430#endif
431      ENDIF
432      !
433      !                            !==========================!
434      ssh(:,:,Kaa) = 0._wp         !==  Set to 0 for AGRIF  ==!
435      !                            !==========================!
436      !
437   END SUBROUTINE rst_read_ssh
438
439   !!=====================================================================
440END MODULE restart
Note: See TracBrowser for help on using the repository browser.