source: trunk/ReadWrite/readoldoparestart.pro @ 44

Last change on this file since 44 was 44, checked in by pinsard, 18 years ago

upgrade of LECTURE/ReadWrite according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files
  • 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   19L*jpk
81   numrec = if size NE filepamameters.size then begin
82      print, 'The size of the file is not the expected one!'
83      print, 'Check your file or the values of ibloc, jpiglo,'
84      print, 'jpjglo, jpk, jpbyt, numrec in this program'
85      return
86   endif
87; first record: six 64-bit integer to read.
88; default definition
89   ino1 = long64(9999)
90   it1 = long64(9999)
91   isor1 = long64(9999)
92   ipcg1 = long64(9999)
93   itke1 = long64(9999)
94   idast1 = long64(9999)
95; read
96   readu, numrst, ino1, it1, isor1, ipcg1, itke1, idast1
97   print, ino1, it1, isor1, ipcg1, itke1, idast1
98; other records
99   params = {jpiglo:jpiglo, jpjglo:jpjglo, jpk:jpk, reclen:reclen}
100;      CALL read3(numrst,ub   ,2 )
101   IF arg_present(ub) THEN ub = read3fromopa(numrst, params, 2)
102;      CALL read3(numrst,vb   ,3 )
103   IF arg_present(vb) THEN vb = read3fromopa(numrst, params, 3)
104;      CALL read3(numrst,tb   ,5 )
105   IF arg_present(tb) THEN tb = read3fromopa(numrst, params, 5)
106;      CALL read3(numrst,sb   ,6 )
107   IF arg_present(sb) THEN sb = read3fromopa(numrst, params, 6)
108;      CALL read3(numrst,rotb ,7 )
109   IF arg_present(rotb) THEN rotb = read3fromopa(numrst, params, 7)
110;      CALL read3(numrst,hdivb,8 )
111   IF arg_present(hdivb) THEN hdivb = read3fromopa(numrst, params, 8)
112;      CALL read3(numrst,un   ,9 )
113   IF arg_present(un) THEN un = read3fromopa(numrst, params, 9)
114;      CALL read3(numrst,vn   ,10)
115   IF arg_present(vn) THEN vn = read3fromopa(numrst, params, 10)
116;      CALL read3(numrst,tn   ,12)
117   IF arg_present(tn) THEN tn = read3fromopa(numrst, params, 12)
118;      CALL read3(numrst,sn   ,13)
119   IF arg_present(sn) THEN sn = read3fromopa(numrst, params, 13)
120;      CALL read3(numrst,rotn ,14)
121   IF arg_present(rotn) THEN rotn = read3fromopa(numrst, params, 14)
122;      CALL read3(numrst,hdivn,15)
123   IF arg_present(hdivn) THEN hdivn = read3fromopa(numrst, params, 15)
124;C
125;C ... Read elliptic solver arrays
126;C
127;      CALL read2(numrst,gcx ,jpk,17)
128   IF arg_present(gcx) THEN gcx = read2fromopa(numrst, params, 17)
129;      CALL read2(numrst,gcxb,jpk,18)
130   IF arg_present(gcxb) THEN gcxb = read2fromopa(numrst, params, 18)
131;C
132;#ifdef key_freesurf_cstvol
133;C
134;C ... free surface formulation (eta)
135;C
136;      CALL read2(numrst,etab ,jpk,4 )
137   IF arg_present(etab) THEN etab = read2fromopa(numrst, params, 4)
138;      CALL read2(numrst,etan ,jpk,4 )
139   IF arg_present(etan) THEN etan = read2fromopa(numrst, params, 4)
140;#  else
141;C
142;C ... Rigid-lid formulation (bsf)
143;C
144;      CALL read2(numrst,bsfb ,jpk,4 )
145   IF arg_present(bsfb) THEN bsfb = read2fromopa(numrst, params, 4)
146;      CALL read2(numrst,bsfn ,jpk,11)
147   IF arg_present(bsfn) THEN bsfn = read2fromopa(numrst, params, 11)
148;      CALL read2(numrst,bsfd ,jpk,16)
149   IF arg_present(bsfd) THENbsfd  = read2fromopa(numrst, params, 16)
150;#endif
151;#ifdef key_zdftke
152;          CALL read3(numrst,en,19)
153   IF arg_present(en) THEN en = read3fromopa(numrst, params, 19)
154
155
156
157   close, numrst
158   free_lun, numrst
159
160   return
161end
162
163
164
165;CDIR$ LIST
166;      SUBROUTINE dtrlec
167;CCC---------------------------------------------------------------------
168;CCC
169;CCC                       ROUTINE dtrlec
170;CCC                     ******************
171;CCC
172;CCC  Purpose :
173;CCC  --------
174;CCC     Read files for restart
175;CCC
176;CC   Method :
177;CC   -------
178;CC      Read the previous fields on the file numrst
179;CC      the first record indicates previous characterics
180;CC      after control with the present run, we read :
181;CC      - prognostic variables on the second record
182;CC      - elliptic solver arrays
183;CC     - barotropic stream function arrays (default option)
184;CC       or free surface arrays ("key_freesurf_cstvol" defined)
185;CC      - tke arrays ("key_zdftke" defined)
186;CC      for this last three records,  the previous characteristics
187;CC      could be different with those used in the present run.
188;CC
189;CC   Input :
190;CC   ------
191;CC      common
192;CC            /comrst/          : restart parameter
193;CC            /comctl/          : parameters for the control
194;CC
195;CC   Output :
196;CC   -------
197;CC      common
198;CC            /combef/          : previous fields (before)
199;CC            /comnow/          : present fields (now)
200;CC            /combsf/          : barotropic stream function
201;CC            /comspg/          : surface pressure
202;CC            /comsol/          : diagonal preconditioned conjugate
203;CC
204;CC   Modifications :
205;CC   --------------
206;CC      original  : 91-03 ()
207;CC      additions : 92-01 (M. Imbard)
208;CC                : 92-06 correction restart file (M. Imbard)
209;CC                : 98-02 (M. Guyon) FETI method
210;CC      addition  : 98-05 (G. Roullet) free surface
211;CC----------------------------------------------------------------------
212;CC parameters and commons
213;CC ======================
214;CDIR$ NOLIST
215;#include "parameter.h"
216;#include "common.h"
217;CDIR$ LIST
218;CC----------------------------------------------------------------------
219;CC local declarations
220;CC ==================
221;      INTEGER ji, jj, jk, jl
222;      INTEGER ino0, it0, ipcg0, isor0, itke0
223;      INTEGER ino1, it1, isor1, ipcg1, itke1, idast1
224;CC----------------------------------------------------------------------
225;CC statement functions
226;CC ===================
227;CDIR$ NOLIST
228;#include "stafun.h"
229;CDIR$ LIST
230;CCC---------------------------------------------------------------------
231;CCC  OPA8, LODYC (1997)
232;CCC---------------------------------------------------------------------
233;C
234;C
235;C 0. Initialisations
236;C ------------------
237;C
238;      ino0  = no
239;      it0   = nit000
240;      ipcg0 = 0
241;      isor0 = 0
242;      itke0 = 0
243;      isor0 = nsolv-1
244;      ipcg0 = 2-nsolv
245;#ifdef key_zdftke
246;      itke0 = 1
247;#endif
248;C FETI method
249;      IF (nsolv .EQ. 3) THEN
250;          isor0=2
251;          ipcg0=2
252;      ENDIF
253;C
254;      IF(lwp) THEN
255;          WRITE(numout,*) ' '
256;          WRITE(numout,*) ' *** dtrlec:  beginning of restart'
257;          WRITE(numout,*) ' '
258;          WRITE(numout,*) ' the present run :'
259;          WRITE(numout,*) '   job number : ', no
260;          WRITE(numout,*) '   with nit000 : ', nit000
261;          WRITE(numout,*) '   with pcg option ipcg0 : ', ipcg0
262;          WRITE(numout,*) '   with sor option isor0 : ', isor0
263;          WRITE(numout,*) '   with FETI solver option ipcg0 & isor0 : ',
264;     $        ipcg0,' & ',isor0
265;          WRITE(numout,*) '   with tke option itke0 : ', itke0
266;      ENDIF
267;C
268;C 1. Read numrst
269;C --------------
270;C
271;C ... First record
272;C
273;      READ(numrst,REC=1) ino1, it1, isor1, ipcg1, itke1, idast1
274;C
275;      IF(lwp) THEN
276;          WRITE(numout,*) ' '
277;          WRITE(numout,*) ' READ numrst with '
278;          WRITE(numout,*) '   job number : ', ino1
279;          WRITE(numout,*) '   with time step it : ', it1
280;          WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
281;          WRITE(numout,*) '   with sor option isor1 : ', isor1
282;          WRITE(numout,*) '   with tke option itke1 : ', itke1
283;          WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ',
284;     $        ipcg1 + isor1
285;          WRITE(numout,*) ' '
286;      ENDIF
287;C
288;C ... Control of date
289;C
290;      IF ( (it0-it1).NE.1 .AND. abs(nrstdt).EQ.1 ) THEN
291;          IF(lwp) THEN
292;              WRITE(numout,*) ' ===>>>> : problem with nit000 for the',
293;     $            ' restart'
294;              WRITE(numout,*) ' =======                              ',
295;     $            ' ======='
296;              WRITE(numout,*) ' we stop. verify the file'
297;              WRITE(numout,*) ' or rerun with the value  0 for the'
298;              WRITE(numout,*) ' control of time parameter  nrstdt'
299;              WRITE(numout,*) ' '
300;          ENDIF
301;          STOP 'dtrlec'
302;      ENDIF
303;      IF ( nrstdt.EQ.1 ) ndate0 = idast1
304;C
305;C ... Read prognostic variables
306;C
307;      CALL read3(numrst,ub   ,2 )
308;      CALL read3(numrst,vb   ,3 )
309;      CALL read3(numrst,tb   ,5 )
310;      CALL read3(numrst,sb   ,6 )
311;      CALL read3(numrst,rotb ,7 )
312;      CALL read3(numrst,hdivb,8 )
313;      CALL read3(numrst,un   ,9 )
314;      CALL read3(numrst,vn   ,10)
315;      CALL read3(numrst,tn   ,12)
316;      CALL read3(numrst,sn   ,13)
317;      CALL read3(numrst,rotn ,14)
318;      CALL read3(numrst,hdivn,15)
319;C
320;C ... Read elliptic solver arrays
321;C
322;      CALL read2(numrst,gcx ,jpk,17)
323;      CALL read2(numrst,gcxb,jpk,18)
324;C
325;#ifdef key_freesurf_cstvol
326;C
327;C ... free surface formulation (eta)
328;C
329;      CALL read2(numrst,etab ,jpk,4 )
330;      CALL read2(numrst,etan ,jpk,4 )
331;#  else
332;C
333;C ... Rigid-lid formulation (bsf)
334;C
335;      CALL read2(numrst,bsfb ,jpk,4 )
336;      CALL read2(numrst,bsfn ,jpk,11)
337;      CALL read2(numrst,bsfd ,jpk,16)
338;#endif
339;C
340;#ifdef key_zdftke
341;C
342;C ... Read tke arrays
343;C
344;      IF(itke1.eq.1) THEN
345;          CALL read3(numrst,en,19)
346;      ELSE
347;          IF(lwp) THEN
348;              WRITE(numout,*) ' ===>>>> : the previous restart file',
349;     $            ' didnt used  tke scheme'
350;              WRITE(numout,*) ' =======                ======='
351;          ENDIF
352;          nrstdt=2
353;      ENDIF
354;#endif
355;C
356;C
357;      RETURN
358;      END
Note: See TracBrowser for help on using the repository browser.