source: trunk/SRC/ReadWrite/readoldoparestart.pro @ 74

Last change on this file since 74 was 74, checked in by smasson, 18 years ago

debug xxx and cie + clean data file + rm perldoc_idl
  • Property svn:executable set to *
File size: 11.6 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:readoldoparestart
6;      based on the OPA subroutine dtrlec included at the end of the
7;      file.
8;
9; PURPOSE:read the old restart files of OPA (before NetCDF)
10;
11; CATEGORY:for OPA before NetCDF
12;
13; CALLING SEQUENCE:readoldoparestart, filename, jpiglo, jpjglo, jpk
14;
15; INPUTS:
16;     filename: with the whole path if necessary
17;     jpiglo, jpjglo, jpk: dimensions of the opa grid
18;
19; KEYWORD PARAMETERS:
20;     IBLOC: ibloc size, default: ibloc = 4096L
21;     JPBYT: jpbyt size, defalut: jpbyt = 8L
22;     NUMREC: number of records in the file. defalut: numrec = 19L*jpk
23;     UB, VB, TB, SB, ROTB, HDIVB, UN, VN, TN, SN, ROTN, HDIVN, GCX,
24;     GCXB, ETAB, TAN, BSFB, BSFN, BSFD, EN: the variable we want to
25;     read.
26;
27; OUTPUTS:according to the given keywords.
28;
29; COMMON BLOCKS:none
30;
31; SIDE EFFECTS:
32;
33; RESTRICTIONS:bug for etab and etan written on the same record???
34;
35; EXAMPLE:
36;
37; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
38;                      June 2002
39;-
40;------------------------------------------------------------
41;------------------------------------------------------------
42;------------------------------------------------------------
43
44
45FUNCTION read2fromopa, unit, params, num
46   offset=params.reclen*params.jpk*(num-1L)
47   a=assoc(unit,dblarr(params.jpiglo,params.jpjglo,/nozero),offset)
48   return, a[0]
49end
50FUNCTION read3fromopa, unit, params, num
51   offset=params.reclen*params.jpk*(num-1L)
52   a=assoc(unit,dblarr(params.jpiglo,params.jpjglo,params.jpk,/nozero),offset)
53   return, a[0]
54end
55
56
57PRO readoldoparestart, filename, jpiglo, jpjglo, jpk, IBLOC = ibloc, JPBYT = jpbyt, NUMREC = numrec, ub = ub, vb = vb, tb = tb, sb = sb, rotb = rotb, hdivb = hdivb, un = un, vn = vn, tn = tn, sn = sn, rotn = rotn, hdivn = hdivn, gcx = gcx, gcxb = gcxb, etab = etab, etan = etan, bsfb = bsfb, bsfn = bsfn, bsfd = bsfd, en = en
58;
59   iname_file = findfile(filename)
60   if iname_file[0] EQ '' then begin
61      print, 'Bad file name'
62      return
63   ENDIF ELSE iname_file = iname_file[0]
64; open the file
65   openr,numrst , iname_file, /get_lun, /swap_if_little_endian
66; check the size of the file
67   filepamameters = fstat(numrst)
68; parameter definition
69   IF keyword_set(ibloc) THEN ibloc = long(ibloc) ELSE ibloc = 4096L
70   jpiglo = long(jpiglo)
71   jpjglo = long(jpjglo)
72   jpk = long(jpk)
73   IF keyword_set(jpbyt) THEN jpbyt = long(jpbyt) ELSE jpbyt = 8L
74; record length computation
75   reclen = ibloc*((jpiglo*jpjglo*jpbyt-1 )/ibloc+1)
76   IF keyword_set(numrec) THEN numrec = long(numrec) ELSE numrec = 19L*jpk
77   toomuch = reclen-jpiglo*jpjglo*jpbyt
78; expected size computation
79   size = numrec*reclen-toomuch
80   if size NE filepamameters.size then begin
81      print, 'The size of the file is not the expected one!'
82      print, 'Check your file or the values of ibloc, jpiglo,'
83      print, 'jpjglo, jpk, jpbyt, numrec in this program'
84      return
85   endif
86; first record: six 64-bit integer to read.
87; default definition
88   ino1 = long64(9999)
89   it1 = long64(9999)
90   isor1 = long64(9999)
91   ipcg1 = long64(9999)
92   itke1 = long64(9999)
93   idast1 = long64(9999)
94; read
95   readu, numrst, ino1, it1, isor1, ipcg1, itke1, idast1
96   print, ino1, it1, isor1, ipcg1, itke1, idast1
97; other records
98   params = {jpiglo:jpiglo, jpjglo:jpjglo, jpk:jpk, reclen:reclen}
99;      CALL read3(numrst,ub   ,2 )
100   IF arg_present(ub) THEN ub = read3fromopa(numrst, params, 2)
101;      CALL read3(numrst,vb   ,3 )
102   IF arg_present(vb) THEN vb = read3fromopa(numrst, params, 3)
103;      CALL read3(numrst,tb   ,5 )
104   IF arg_present(tb) THEN tb = read3fromopa(numrst, params, 5)
105;      CALL read3(numrst,sb   ,6 )
106   IF arg_present(sb) THEN sb = read3fromopa(numrst, params, 6)
107;      CALL read3(numrst,rotb ,7 )
108   IF arg_present(rotb) THEN rotb = read3fromopa(numrst, params, 7)
109;      CALL read3(numrst,hdivb,8 )
110   IF arg_present(hdivb) THEN hdivb = read3fromopa(numrst, params, 8)
111;      CALL read3(numrst,un   ,9 )
112   IF arg_present(un) THEN un = read3fromopa(numrst, params, 9)
113;      CALL read3(numrst,vn   ,10)
114   IF arg_present(vn) THEN vn = read3fromopa(numrst, params, 10)
115;      CALL read3(numrst,tn   ,12)
116   IF arg_present(tn) THEN tn = read3fromopa(numrst, params, 12)
117;      CALL read3(numrst,sn   ,13)
118   IF arg_present(sn) THEN sn = read3fromopa(numrst, params, 13)
119;      CALL read3(numrst,rotn ,14)
120   IF arg_present(rotn) THEN rotn = read3fromopa(numrst, params, 14)
121;      CALL read3(numrst,hdivn,15)
122   IF arg_present(hdivn) THEN hdivn = read3fromopa(numrst, params, 15)
123;C
124;C ... Read elliptic solver arrays
125;C
126;      CALL read2(numrst,gcx ,jpk,17)
127   IF arg_present(gcx) THEN gcx = read2fromopa(numrst, params, 17)
128;      CALL read2(numrst,gcxb,jpk,18)
129   IF arg_present(gcxb) THEN gcxb = read2fromopa(numrst, params, 18)
130;C
131;#ifdef key_freesurf_cstvol
132;C
133;C ... free surface formulation (eta)
134;C
135;      CALL read2(numrst,etab ,jpk,4 )
136   IF arg_present(etab) THEN etab = read2fromopa(numrst, params, 4)
137;      CALL read2(numrst,etan ,jpk,4 )
138   IF arg_present(etan) THEN etan = read2fromopa(numrst, params, 4)
139;#  else
140;C
141;C ... Rigid-lid formulation (bsf)
142;C
143;      CALL read2(numrst,bsfb ,jpk,4 )
144   IF arg_present(bsfb) THEN bsfb = read2fromopa(numrst, params, 4)
145;      CALL read2(numrst,bsfn ,jpk,11)
146   IF arg_present(bsfn) THEN bsfn = read2fromopa(numrst, params, 11)
147;      CALL read2(numrst,bsfd ,jpk,16)
148   IF arg_present(bsfd) THEN bsfd  = read2fromopa(numrst, params, 16)
149;#endif
150;#ifdef key_zdftke
151;          CALL read3(numrst,en,19)
152   IF arg_present(en) THEN en = read3fromopa(numrst, params, 19)
153
154
155
156   close, numrst
157   free_lun, numrst
158
159   return
160end
161
162
163
164;CDIR$ LIST
165;      SUBROUTINE dtrlec
166;CCC---------------------------------------------------------------------
167;CCC
168;CCC                       ROUTINE dtrlec
169;CCC                     ******************
170;CCC
171;CCC  Purpose :
172;CCC  --------
173;CCC     Read files for restart
174;CCC
175;CC   Method :
176;CC   -------
177;CC      Read the previous fields on the file numrst
178;CC      the first record indicates previous characterics
179;CC      after control with the present run, we read :
180;CC      - prognostic variables on the second record
181;CC      - elliptic solver arrays
182;CC     - barotropic stream function arrays (default option)
183;CC       or free surface arrays ("key_freesurf_cstvol" defined)
184;CC      - tke arrays ("key_zdftke" defined)
185;CC      for this last three records,  the previous characteristics
186;CC      could be different with those used in the present run.
187;CC
188;CC   Input :
189;CC   ------
190;CC      common
191;CC            /comrst/          : restart parameter
192;CC            /comctl/          : parameters for the control
193;CC
194;CC   Output :
195;CC   -------
196;CC      common
197;CC            /combef/          : previous fields (before)
198;CC            /comnow/          : present fields (now)
199;CC            /combsf/          : barotropic stream function
200;CC            /comspg/          : surface pressure
201;CC            /comsol/          : diagonal preconditioned conjugate
202;CC
203;CC   Modifications :
204;CC   --------------
205;CC      original  : 91-03 ()
206;CC      additions : 92-01 (M. Imbard)
207;CC                : 92-06 correction restart file (M. Imbard)
208;CC                : 98-02 (M. Guyon) FETI method
209;CC      addition  : 98-05 (G. Roullet) free surface
210;CC----------------------------------------------------------------------
211;CC parameters and commons
212;CC ======================
213;CDIR$ NOLIST
214;#include "parameter.h"
215;#include "common.h"
216;CDIR$ LIST
217;CC----------------------------------------------------------------------
218;CC local declarations
219;CC ==================
220;      INTEGER ji, jj, jk, jl
221;      INTEGER ino0, it0, ipcg0, isor0, itke0
222;      INTEGER ino1, it1, isor1, ipcg1, itke1, idast1
223;CC----------------------------------------------------------------------
224;CC statement functions
225;CC ===================
226;CDIR$ NOLIST
227;#include "stafun.h"
228;CDIR$ LIST
229;CCC---------------------------------------------------------------------
230;CCC  OPA8, LODYC (1997)
231;CCC---------------------------------------------------------------------
232;C
233;C
234;C 0. Initialisations
235;C ------------------
236;C
237;      ino0  = no
238;      it0   = nit000
239;      ipcg0 = 0
240;      isor0 = 0
241;      itke0 = 0
242;      isor0 = nsolv-1
243;      ipcg0 = 2-nsolv
244;#ifdef key_zdftke
245;      itke0 = 1
246;#endif
247;C FETI method
248;      IF (nsolv .EQ. 3) THEN
249;          isor0=2
250;          ipcg0=2
251;      ENDIF
252;C
253;      IF(lwp) THEN
254;          WRITE(numout,*) ' '
255;          WRITE(numout,*) ' *** dtrlec:  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 : ',
263;     $        ipcg0,' & ',isor0
264;          WRITE(numout,*) '   with tke option itke0 : ', itke0
265;      ENDIF
266;C
267;C 1. Read numrst
268;C --------------
269;C
270;C ... First record
271;C
272;      READ(numrst,REC=1) ino1, it1, isor1, ipcg1, itke1, idast1
273;C
274;      IF(lwp) THEN
275;          WRITE(numout,*) ' '
276;          WRITE(numout,*) ' READ numrst with '
277;          WRITE(numout,*) '   job number : ', ino1
278;          WRITE(numout,*) '   with time step it : ', it1
279;          WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
280;          WRITE(numout,*) '   with sor option isor1 : ', isor1
281;          WRITE(numout,*) '   with tke option itke1 : ', itke1
282;          WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ',
283;     $        ipcg1 + isor1
284;          WRITE(numout,*) ' '
285;      ENDIF
286;C
287;C ... Control of date
288;C
289;      IF ( (it0-it1).NE.1 .AND. abs(nrstdt).EQ.1 ) THEN
290;          IF(lwp) THEN
291;              WRITE(numout,*) ' ===>>>> : problem with nit000 for the',
292;     $            ' restart'
293;              WRITE(numout,*) ' =======                              ',
294;     $            ' ======='
295;              WRITE(numout,*) ' we stop. verify the file'
296;              WRITE(numout,*) ' or rerun with the value  0 for the'
297;              WRITE(numout,*) ' control of time parameter  nrstdt'
298;              WRITE(numout,*) ' '
299;          ENDIF
300;          STOP 'dtrlec'
301;      ENDIF
302;      IF ( nrstdt.EQ.1 ) ndate0 = idast1
303;C
304;C ... Read prognostic variables
305;C
306;      CALL read3(numrst,ub   ,2 )
307;      CALL read3(numrst,vb   ,3 )
308;      CALL read3(numrst,tb   ,5 )
309;      CALL read3(numrst,sb   ,6 )
310;      CALL read3(numrst,rotb ,7 )
311;      CALL read3(numrst,hdivb,8 )
312;      CALL read3(numrst,un   ,9 )
313;      CALL read3(numrst,vn   ,10)
314;      CALL read3(numrst,tn   ,12)
315;      CALL read3(numrst,sn   ,13)
316;      CALL read3(numrst,rotn ,14)
317;      CALL read3(numrst,hdivn,15)
318;C
319;C ... Read elliptic solver arrays
320;C
321;      CALL read2(numrst,gcx ,jpk,17)
322;      CALL read2(numrst,gcxb,jpk,18)
323;C
324;#ifdef key_freesurf_cstvol
325;C
326;C ... free surface formulation (eta)
327;C
328;      CALL read2(numrst,etab ,jpk,4 )
329;      CALL read2(numrst,etan ,jpk,4 )
330;#  else
331;C
332;C ... Rigid-lid formulation (bsf)
333;C
334;      CALL read2(numrst,bsfb ,jpk,4 )
335;      CALL read2(numrst,bsfn ,jpk,11)
336;      CALL read2(numrst,bsfd ,jpk,16)
337;#endif
338;C
339;#ifdef key_zdftke
340;C
341;C ... Read tke arrays
342;C
343;      IF(itke1.eq.1) THEN
344;          CALL read3(numrst,en,19)
345;      ELSE
346;          IF(lwp) THEN
347;              WRITE(numout,*) ' ===>>>> : the previous restart file',
348;     $            ' didnt used  tke scheme'
349;              WRITE(numout,*) ' =======                ======='
350;          ENDIF
351;          nrstdt=2
352;      ENDIF
353;#endif
354;C
355;C
356;      RETURN
357;      END
Note: See TracBrowser for help on using the repository browser.