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_dimg.h90 in trunk/NEMO/OPA_SRC – NEMO

source: trunk/NEMO/OPA_SRC/restart_dimg.h90 @ 359

Last change on this file since 359 was 359, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1   !!---------------------------------------------------------------------
2   !!                     ***  restart_dimg.h90  ***
3   !!---------------------------------------------------------------------
4   !!  OPA 9.0 , LOCEAN-IPSL (2005)
5   !! $Header$
6   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
7   !!----------------------------------------------------------------------
8
9   SUBROUTINE rst_write(kt)
10     !!---------------------------------------------------------------------
11     !!                  ***  ROUTINE rst_write  ***
12     !!
13     !! ** Purpose :   Write restart fields in direct access format in mpp.
14     !!      one per process
15     !!
16     !! ** Method  :   each nstock time step , save fields which are necessary
17     !!      for restart
18     !!      Record #1 hold general information on the state of the run
19     !!      Data fields (either 3D or 2D ) starts ar record #2
20     !!
21     !! History :
22     !!        !  91-03  ()  original code
23     !!        !  91-11  (G. Madec)
24     !!        !  92-06  (M. Imbard)  correction restart file
25     !!        !  92-07  (M. Imbard)  split into diawri and rstwri
26     !!        !  98-02  (M. Guyon)  FETI method
27     !!        !  98-05  (G. Roullet)  free surface
28     !!        !  99-11  (M. Imbard)  NetCDF FORMAT with ioipsl
29     !!   8.5  !  03-06  (J.M. Molines)  F90: Free form, mpp support
30     !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
31     !!----------------------------------------------------------------------
32     !! * Arguments
33     INTEGER, INTENT( in ) ::   kt         ! ocean time-step
34
35     !! * Local declarations
36     INTEGER :: ino0, it0, ipcg0, isor0, itke0
37     INTEGER :: irecl8, irec
38     INTEGER :: jk               ! dummy loop indices
39     INTEGER :: inum = 11        ! temporary logical unit
40     INTEGER :: ios1 , ios2      ! flag for ice and bulk in the current run
41     INTEGER :: ios3             ! flag for free surface.  0 = none 1 = yes.  0 = none 1 = yes
42     INTEGER :: ios4             ! flag for coupled (1) or not (0)
43
44     CHARACTER(LEN=80)  :: clres
45
46     REAL(wp), DIMENSION( 1) ::   zfice, zfblk   ! used only in case of ice & bulk
47     !!----------------------------------------------------------------------
48
49    IF( MOD( kt, nstock ) == 0 .OR. kt == nitend ) THEN
50
51       ! 0. Initializations
52       ! ------------------
53
54       IF(lwp) THEN
55          WRITE(numout,*) ' '
56          WRITE(numout,*) ' rst_write: output done in inum = ',   &
57               inum,' at it= ',kt,' date= ',ndastp
58          WRITE(numout,*) ' -------'
59       ENDIF
60
61       ! Open direct access file, with reclength for 2D wp fields
62       irecl8= jpi * jpj * wp
63       WRITE(clres,'(a,i3.3)') 'restart.output.',narea
64       OPEN(inum,FILE=clres,FORM='UNFORMATTED', ACCESS='DIRECT', RECL=irecl8 )
65
66
67       ino0  = no
68       it0   = kt
69       ipcg0 = 0
70       isor0 = 0
71       itke0 = 0
72       isor0 = nsolv - 1
73       ipcg0 = 2 - nsolv
74       IF ( lk_zdftke ) itke0=1
75       ! FETI method
76       IF (nsolv == 3) THEN
77          isor0 = 2
78          ipcg0 = 2
79       ENDIF
80
81       ! 1. Write in inum
82       ! ------------------
83
84       ! first record
85       ios1 = 0
86       ios2 = 0
87       ios3 = 0
88       ios4 = 0
89       IF ( lk_ice_lim                           )   ios1 = 1
90       IF ( l_bulk                               )   ios2 = 1
91       IF ( lk_dynspg_flt                        )   ios3 = 1
92       IF ( lk_cpl                               )   ios4 = 1
93
94       WRITE(inum,REC=1) irecl8, ino0, it0, isor0, ipcg0, itke0, &
95          &              nfice, nfbulk , ios1, ios2, ios3, ios4, &
96          &              ndastp, adatrj, jpi, jpj, jpk,  &
97          &              jpni, jpnj, jpnij, narea, jpiglo, jpjglo, &
98          &              nlcit, nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt
99
100       ! prognostic variables
101
102       irec=2
103
104       ! 'before' fields
105       DO jk = 1, jpk
106          WRITE(inum,REC=irec) ub(:,:,jk)   ;    irec = irec +1
107       END DO
108       DO jk = 1, jpk
109          WRITE(inum,REC=irec) vb(:,:,jk)   ;    irec = irec +1
110       END DO
111       DO jk = 1, jpk
112          WRITE(inum,REC=irec) tb(:,:,jk)   ;    irec = irec +1
113       END DO
114       DO jk = 1, jpk
115          WRITE(inum,REC=irec) sb(:,:,jk)   ;    irec = irec +1
116       END DO
117       DO jk = 1, jpk
118          WRITE(inum,REC=irec) rotb(:,:,jk)   ;    irec = irec +1
119       END DO
120       DO jk = 1, jpk
121          WRITE(inum,REC=irec) hdivb(:,:,jk)   ;    irec = irec +1
122       END DO
123
124       ! 'now' fields
125       DO jk = 1, jpk
126          WRITE(inum,REC=irec) un(:,:,jk)   ;   irec = irec +1
127       END DO
128       DO jk = 1, jpk
129          WRITE(inum,REC=irec) vn(:,:,jk)   ;   irec = irec +1
130       END DO
131       DO jk = 1, jpk
132          WRITE(inum,REC=irec) tn(:,:,jk)   ;   irec = irec +1
133       END DO
134       DO jk = 1, jpk
135          WRITE(inum,REC=irec) sn(:,:,jk)   ;   irec = irec +1
136       END DO
137       DO jk = 1, jpk
138          WRITE(inum,REC=irec) rotn(:,:,jk)   ;   irec = irec +1
139       END DO
140       DO jk = 1, jpk
141          WRITE(inum,REC=irec) hdivn(:,:,jk)   ;   irec = irec +1
142       END DO
143
144       ! elliptic solver arrays
145       WRITE(inum,REC=irec ) gcx(1:jpi,1:jpj)   ;   irec = irec +1
146       WRITE(inum,REC=irec ) gcxb(1:jpi,1:jpj)   ;   irec = irec +1
147#if defined key_dynspg_rl
148       ! Rigid-lid formulation (bsf)
149       WRITE(inum,REC=irec ) bsfb(:,:)   ;   irec = irec +1
150       WRITE(inum,REC=irec ) bsfn(:,:)   ;   irec = irec +1
151       WRITE(inum,REC=irec ) bsfd(:,:)   ;   irec = irec +1
152# else
153       ! free surface formulation (ssh)
154       WRITE(inum,REC=irec ) sshb(:,:)   ;   irec = irec +1
155       WRITE(inum,REC=irec ) sshn(:,:)   ;   irec = irec +1
156# if defined key_dynspg_ts
157       ! free surface formulation issued from barotropic loop
158       WRITE(inum,REC=irec ) sshb_b(:,:)   ;   irec = irec +1
159       WRITE(inum,REC=irec ) sshn_b(:,:)   ;   irec = irec +1
160
161       ! horizontal transports issued from barotropic loop
162       WRITE(inum,REC=irec) un_b(:,:)   ;   irec = irec +1
163       WRITE(inum,REC=irec) vn_b(:,:)   ;   irec = irec +1
164# endif
165#endif
166
167       ! TKE arrays
168#if defined key_zdftke
169         DO jk = 1, jpk
170            WRITE(inum,REC=irec) en(:,:,jk)   ;   irec = irec + 1
171         END DO
172#endif
173
174#if defined key_ice_lim
175          zfice(1) = FLOAT( nfice )                                      ! Louvain La Neuve Sea Ice Model
176          WRITE(inum,REC=irec) zfice(:)      ;   irec = irec + 1
177          WRITE(inum,REC=irec) sst_io(:,:)   ;   irec = irec + 1
178          WRITE(inum,REC=irec) sss_io(:,:)   ;   irec = irec + 1
179          WRITE(inum,REC=irec) u_io  (:,:)   ;   irec = irec + 1
180          WRITE(inum,REC=irec) v_io  (:,:)   ;   irec = irec + 1
181#    if defined key_coupled
182          WRITE(inum,REC=irec) alb_ice(:,:)  ;   irec = irec + 1
183#    endif
184#endif
185# if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
186          zfblk(1) = FLOAT( nfbulk )                                 ! Bulk
187          WRITE(inum,REC=irec) zfblk(:)   ;   irec = irec + 1
188          WRITE(inum,REC=irec) gsst(:,:)  ;   irec = irec + 1
189# endif
190
191       CLOSE(inum)
192    ENDIF
193
194  END SUBROUTINE rst_write
195
196
197  SUBROUTINE rst_read
198    !!---------------------------------------------------------------------
199    !!                  ***  ROUTINE rst_read  ***
200    !! ** Purpose :
201    !!        Read restart fields in direct access format, one per process
202    !!
203    !! ** Method :   Just does the opposit than rst_wri
204    !!
205    !! History :
206    !!        !  91-03  ()  original code
207    !!        !  91-11  (G. Madec)
208    !!        !  92-06  (M. Imbard)  correction restart file
209    !!        !  92-07  (M. Imbard)  split into diawri and rstwri
210    !!        !  98-02  (M. Guyon)  FETI method
211    !!        !  98-05  (G. Roullet)  free surface
212    !!        !  99-11  (M. Imbard)  NetCDF FORMAT with ioipsl
213    !!   8.5  !  03-06  (J.M. Molines)  F90: Free form, mpp support
214    !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
215    !!----------------------------------------------------------------------
216    USE lib_mpp
217    !! * Local declarations
218    INTEGER :: ino0, it0, ipcg0, isor0, itke0
219    INTEGER :: ino1, it1, isor1, ipcg1, itke1, idast1
220    INTEGER :: iice1, ibulk1
221    INTEGER :: ipi,ipj,ipk, ipni,ipnj,ipnij,iarea
222    INTEGER :: irecl8, irec
223    INTEGER :: ji,jj,jk
224    INTEGER :: ick, inum
225    INTEGER :: ios1, ios2, ios3, ios4
226
227    CHARACTER(LEN=80)  :: clres
228
229    LOGICAL   :: lstop
230
231    REAL(wp), DIMENSION( 1) ::   zfice, zfblk   ! used only in case of ice & bulk
232    !!----------------------------------------------------------------------
233
234
235    ! 0. Initialisations
236    ! ------------------
237
238    inum = 11
239    ino0  = no
240    it0   = nit000
241    ipcg0 = 0
242    isor0 = 0
243    itke0 = 0
244    isor0 = nsolv-1
245    ipcg0 = 2-nsolv
246    IF (lk_zdftke ) itke0 = 1
247    ! FETI method
248    IF( nsolv == 3 ) THEN
249       isor0=2
250       ipcg0=2
251    ENDIF
252
253    IF(lwp) THEN
254       WRITE(numout,*)
255       WRITE(numout,*) ' *** rst_read:  beginning of restart'
256       WRITE(numout,*) ' '
257       WRITE(numout,*) ' the present run :'
258       WRITE(numout,*) '   job number : ', no
259       WRITE(numout,*) '   with nit000 : ', nit000
260       WRITE(numout,*) '   with pcg option ipcg0 : ', ipcg0
261       WRITE(numout,*) '   with sor option isor0 : ', isor0
262       WRITE(numout,*) '   with FETI solver option ipcg0 & isor0 : ', ipcg0,' & ',isor0
263       WRITE(numout,*) '   with tke option itke0 : ', itke0
264       WRITE(numout,*) '   with nfice            : ', nfice
265       WRITE(numout,*) '   with nfbulk           : ', nfbulk
266    ENDIF
267
268    ! Open direct access file, with reclength for 2D wp fields
269    WRITE(clres,'(a,i3.3)') 'restart.',narea
270
271    OPEN(inum,FILE=clres,FORM='UNFORMATTED', ACCESS='DIRECT', RECL=8 )
272    READ(inum,REC=1)irecl8
273    CLOSE(inum)
274
275    OPEN(inum,FILE=clres,FORM='UNFORMATTED', ACCESS='DIRECT', RECL=irecl8 )
276
277    ! 1. Read inum
278    ! --------------
279
280    READ(inum,REC=1) irecl8, ino1, it1, isor1, ipcg1, itke1,   &
281       &             iice1, ibulk1, ios1, ios2, ios3, ios4,    &
282       &             idast1, adatrj0,  ipi,ipj,ipk,ipni,ipnj,ipnij,iarea
283
284    ! Performs checks on the file
285
286    ! processor layout changed ? check only on lwp
287    lstop =.FALSE.
288    IF ( ipni /= jpni ) THEN
289       lstop=.TRUE.
290       IF (lwp)  WRITE(numout,*) ' E R R O R : Processor splitting change along I '
291       IF (lwp)  WRITE(numout,*) ' =========='
292    END IF
293
294    IF ( ipnj /= jpnj ) THEN
295       lstop=.TRUE.
296       IF (lwp)  WRITE(numout,*) ' E R R O R : Processor splitting change along J '
297       IF (lwp)  WRITE(numout,*) ' =========='
298    END IF
299
300    IF ( ipnij /= jpnij ) THEN
301       lstop=.TRUE.
302       IF (lwp)  WRITE(numout,*) ' E R R O R : Total number of processors changed '
303       IF (lwp)  WRITE(numout,*) ' =========='
304    END IF
305
306    ick = narea -iarea
307    CALL mpp_sum( ick )
308
309    IF (ick /= 0 ) THEN
310       lstop=.TRUE.
311       IF (lwp)  WRITE(numout,*) ' E R R O R : mismatch in area numbering ...'
312       IF (lwp)  WRITE(numout,*) ' =========='
313    END IF
314
315    IF(lwp) THEN
316       WRITE(numout,*)
317       WRITE(numout,*) ' READ inum with '
318       WRITE(numout,*) '   job number : ', ino1
319       WRITE(numout,*) '   with time step it : ', it1
320       WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
321       WRITE(numout,*) '   with sor option isor1 : ', isor1
322       WRITE(numout,*) '   with tke option itke1 : ', itke1
323       WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ', ipcg1 + isor1
324       WRITE(numout,*)
325    ENDIF
326
327    ! Control of date
328
329    IF( (it0-it1) /= 1 .AND. nrstdt /= 0 ) THEN
330       lstop=.TRUE.
331       IF(lwp) THEN
332          WRITE(numout,*) ' ===>>>> : problem with nit000 for the restart'
333          WRITE(numout,*) ' =======                               ======='
334          WRITE(numout,*) ' we stop. verify the file'
335          WRITE(numout,*) ' or rerun with the value  0 for the'
336          WRITE(numout,*) ' control of time parameter  nrstdt'
337          WRITE(numout,*)
338       ENDIF
339    ENDIF
340
341    IF ( nrstdt /= 2 ) THEN  !  Compatibility with OPA8
342                             !  the beginning of the new run is ndate0 read in the namelist
343                             !  adatrj0 is recalculated assuming constant time step.
344
345       adatrj0 =  ( FLOAT( nit000-1 ) * rdttra(1) ) / rday
346    ELSE                 !  restart option nrstdt = 2
347                         !  both adatrj0 and ndastp are read in the restart file.
348       ndastp = idast1
349    ENDIF
350
351    IF (lstop ) STOP 'rst_read'
352
353    irec=2
354
355    ! 'before' fields
356    DO jk = 1, jpk
357       READ(inum,REC=irec) ub(:,:,jk)   ;   irec = irec +1
358    END DO
359    DO jk = 1, jpk
360       READ(inum,REC=irec) vb(:,:,jk)   ;   irec = irec +1
361    END DO
362    DO jk = 1, jpk
363       READ(inum,REC=irec) tb(:,:,jk)   ;   irec = irec +1
364    END DO
365    DO jk = 1, jpk
366       READ(inum,REC=irec) sb(:,:,jk)   ;   irec = irec +1
367    END DO
368    DO jk = 1, jpk
369       READ(inum,REC=irec) rotb(:,:,jk)   ;   irec = irec +1
370    END DO
371    DO jk = 1, jpk
372       READ(inum,REC=irec) hdivb(:,:,jk)   ;   irec = irec +1
373    END DO
374
375    ! 'now' fields
376    DO jk = 1, jpk
377       READ(inum,REC=irec) un(:,:,jk)   ;   irec = irec +1
378    END DO
379    DO jk = 1, jpk
380       READ(inum,REC=irec) vn(:,:,jk)   ;   irec = irec +1
381    END DO
382    DO jk = 1, jpk
383       READ(inum,REC=irec) tn(:,:,jk)   ;   irec = irec +1
384    END DO
385    DO jk = 1, jpk
386       READ(inum,REC=irec) sn(:,:,jk)   ;   irec = irec +1
387    END DO
388    DO jk = 1, jpk
389       READ(inum,REC=irec) rotn(:,:,jk)   ;   irec = irec +1
390    END DO
391    DO jk = 1, jpk
392       READ(inum,REC=irec) hdivn(:,:,jk)   ;   irec = irec +1
393    END DO
394
395    ! elliptic solver arrays
396    READ(inum,REC=irec ) gcx(1:jpi,1:jpj)   ;   irec = irec +1
397    READ(inum,REC=irec ) gcxb(1:jpi,1:jpj)   ;   irec = irec +1
398#if defined key_dynspg_rl
399    ! Rigid-lid formulation (bsf)
400    READ(inum,REC=irec ) bsfb(:,:)   ;   irec = irec +1
401    READ(inum,REC=irec ) bsfn(:,:)   ;   irec = irec +1
402    READ(inum,REC=irec ) bsfd(:,:)   ;   irec = irec +1
403#else
404    ! free surface formulation (eta)
405    READ(inum,REC=irec ) sshb(:,:)   ;   irec = irec +1
406    READ(inum,REC=irec ) sshn(:,:)   ;   irec = irec +1
407# if defined key_dynspg_ts
408    ! free surface formulation issued from barotropic loop
409    READ(inum,REC=irec ) sshb_b(:,:)   ;   irec = irec +1
410    READ(inum,REC=irec ) sshn_b(:,:)   ;   irec = irec +1
411    ! horizontal transports issued from barotropic loop
412    READ(inum,REC=irec) un_b(:,:)   ;   irec = irec +1
413    READ(inum,REC=irec) vn_b(:,:)   ;   irec = irec +1
414# endif
415#endif
416
417    ! TKE arrays
418#if defined key_zdftke
419    IF ( itke1 == 1 ) THEN
420       DO jk = 1, jpk
421          READ(inum,REC=irec) en(:,:,jk)   ;   irec = irec +1
422       END DO
423    ELSE
424       IF(lwp) THEN
425          WRITE(numout,*) ' ===>>>> : the previous restart file didnot used  tke scheme'
426          WRITE(numout,*) ' =======                ======='
427       ENDIF
428       nrstdt = 2
429    ENDIF
430#endif
431#if defined key_ice_lim
432    ! Louvain La Neuve Sea Ice Model
433    ! check if it was in the previous run
434    IF ( ios1 == 1 ) THEN
435       READ(inum,REC=irec) zfice(:)      ;   irec = irec + 1
436       READ(inum,REC=irec) sst_io(:,:)   ;   irec = irec + 1
437       READ(inum,REC=irec) sss_io(:,:)   ;   irec = irec + 1
438       READ(inum,REC=irec) u_io  (:,:)   ;   irec = irec + 1
439       READ(inum,REC=irec) v_io  (:,:)   ;   irec = irec + 1
440# if defined key_coupled
441       READ(inum,REC=irec) alb_ice(:,:)   ;   irec = irec + 1
442# endif
443    ENDIF
444    IF ( zfice(1) /= FLOAT(nfice) .OR. ios1 == 0 ) THEN
445         IF(lwp) WRITE(numout,*)
446         IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization'
447         IF(lwp) WRITE(numout,*)
448         sst_io(:,:) = sst_io(:,:) + ( nfice-1 )*( tn(:,:,1) + rt0 )
449         sss_io(:,:) = sss_io(:,:) + ( nfice-1 )*  sn(:,:,1)
450         DO jj = 2, jpj
451            DO ji = 2, jpi
452               u_io(ji,jj) = u_io(ji,jj) + (nfice-1)*0.5*( un(ji-1,jj,1)+un(ji-1,jj-1,1) )
453               v_io(ji,jj) = v_io(ji,jj) + (nfice-1)*0.5*( vn(ji,jj-1,1)+vn(ji-1,jj-1,1) )
454            END DO
455         END DO
456# if defined key_coupled
457         alb_ice(:,:) = 0.8 * tmask(:,:,1)
458# endif
459    ENDIF
460#endif
461#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
462      ! bulk forcing
463      IF( ios2 == 1 ) THEN
464         READ(inum,REC=irec) zfblk(:)   ; irec = irec + 1
465         READ(inum,REC=irec) gsst (:,:) ; irec = irec + 1
466      ENDIF
467      IF( zfblk(1) /= FLOAT(nfbulk)  .OR. ios2 == 0 ) THEN
468         IF(lwp) WRITE(numout,*)
469         IF(lwp) WRITE(numout,*) 'rst_read :  Bulk forcing ==> Initialization '
470         IF(lwp) WRITE(numout,*)
471         gsst(:,:) = 0.e0
472         gsst(:,:) = gsst(:,:) + ( nfbulk-1 )*( tn(:,:,1) + rt0 )
473      ENDIF
474#endif
475    CLOSE(inum)
476  ! In case of restart with neuler = 0 then put all before fields = to now fields
477    IF ( neuler == 0 ) THEN
478      tb(:,:,:)=tn(:,:,:)
479      sb(:,:,:)=sn(:,:,:)
480      ub(:,:,:)=un(:,:,:)
481      vb(:,:,:)=vn(:,:,:)
482      rotb(:,:,:)=rotn(:,:,:)
483      hdivb(:,:,:)=hdivn(:,:,:)
484#if defined key_dynspg_rl
485      bsfb(:,:)=bsfn(:,:)      ! rigid lid
486#else
487      sshb(:,:)=sshn(:,:)      ! free surface formulation (eta)
488#endif
489    ENDIF
490
491  END SUBROUTINE rst_read
Note: See TracBrowser for help on using the repository browser.