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

Last change on this file since 269 was 262, checked in by pinsard, 17 years ago

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