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

Last change on this file since 163 was 163, checked in by navarro, 18 years ago

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