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

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

new compilation options (compile_opt idl2, strictarrsubs) in each routine

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