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

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

replace some print by some report in some .pro (continuation)
  • 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;
143  compile_opt idl2, strictarrsubs
144;
145   iname_file = findfile(filename)
146   if iname_file[0] EQ '' then begin
147      ras = report('Bad file name')
148      return
149   ENDIF ELSE iname_file = iname_file[0]
150; open the file
151   openr,numrst , iname_file, /get_lun, /swap_if_little_endian
152; check the size of the file
153   filepamameters = fstat(numrst)
154; parameter definition
155   IF keyword_set(ibloc) THEN ibloc = long(ibloc) ELSE ibloc = 4096L
156   jpiglo = long(jpiglo)
157   jpjglo = long(jpjglo)
158   jpk = long(jpk)
159   IF keyword_set(jpbyt) THEN jpbyt = long(jpbyt) ELSE jpbyt = 8L
160; record length computation
161   reclen = ibloc*((jpiglo*jpjglo*jpbyt-1 )/ibloc+1)
162   IF keyword_set(numrec) THEN numrec = long(numrec) ELSE numrec = 19L*jpk
163   toomuch = reclen-jpiglo*jpjglo*jpbyt
164; expected size computation
165   size = numrec*reclen-toomuch
166   if size NE filepamameters.size then begin
167      ras = report(['The size of the file is not the expected one!', $
168      'Check your file or the values of ibloc, jpiglo,', $
169      'jpjglo, jpk, jpbyt, numrec in this program'])
170      return
171   endif
172; first record: six 64-bit integer to read.
173; default definition
174   ino1 = long64(9999)
175   it1 = long64(9999)
176   isor1 = long64(9999)
177   ipcg1 = long64(9999)
178   itke1 = long64(9999)
179   idast1 = long64(9999)
180; read
181   readu, numrst, ino1, it1, isor1, ipcg1, itke1, idast1
182   print, ino1, it1, isor1, ipcg1, itke1, idast1
183; other records
184   params = {jpiglo:jpiglo, jpjglo:jpjglo, jpk:jpk, reclen:reclen}
185;      CALL read3(numrst,ub   ,2 )
186   IF arg_present(ub) THEN ub = read3fromopa(numrst, params, 2)
187;      CALL read3(numrst,vb   ,3 )
188   IF arg_present(vb) THEN vb = read3fromopa(numrst, params, 3)
189;      CALL read3(numrst,tb   ,5 )
190   IF arg_present(tb) THEN tb = read3fromopa(numrst, params, 5)
191;      CALL read3(numrst,sb   ,6 )
192   IF arg_present(sb) THEN sb = read3fromopa(numrst, params, 6)
193;      CALL read3(numrst,rotb ,7 )
194   IF arg_present(rotb) THEN rotb = read3fromopa(numrst, params, 7)
195;      CALL read3(numrst,hdivb,8 )
196   IF arg_present(hdivb) THEN hdivb = read3fromopa(numrst, params, 8)
197;      CALL read3(numrst,un   ,9 )
198   IF arg_present(un) THEN un = read3fromopa(numrst, params, 9)
199;      CALL read3(numrst,vn   ,10)
200   IF arg_present(vn) THEN vn = read3fromopa(numrst, params, 10)
201;      CALL read3(numrst,tn   ,12)
202   IF arg_present(tn) THEN tn = read3fromopa(numrst, params, 12)
203;      CALL read3(numrst,sn   ,13)
204   IF arg_present(sn) THEN sn = read3fromopa(numrst, params, 13)
205;      CALL read3(numrst,rotn ,14)
206   IF arg_present(rotn) THEN rotn = read3fromopa(numrst, params, 14)
207;      CALL read3(numrst,hdivn,15)
208   IF arg_present(hdivn) THEN hdivn = read3fromopa(numrst, params, 15)
209;C
210;C ... Read elliptic solver arrays
211;C
212;      CALL read2(numrst,gcx ,jpk,17)
213   IF arg_present(gcx) THEN gcx = read2fromopa(numrst, params, 17)
214;      CALL read2(numrst,gcxb,jpk,18)
215   IF arg_present(gcxb) THEN gcxb = read2fromopa(numrst, params, 18)
216;C
217;#ifdef key_freesurf_cstvol
218;C
219;C ... free surface formulation (eta)
220;C
221;      CALL read2(numrst,etab ,jpk,4 )
222   IF arg_present(etab) THEN etab = read2fromopa(numrst, params, 4)
223;      CALL read2(numrst,etan ,jpk,4 )
224   IF arg_present(etan) THEN etan = read2fromopa(numrst, params, 4)
225;#  else
226;C
227;C ... Rigid-lid formulation (bsf)
228;C
229;      CALL read2(numrst,bsfb ,jpk,4 )
230   IF arg_present(bsfb) THEN bsfb = read2fromopa(numrst, params, 4)
231;      CALL read2(numrst,bsfn ,jpk,11)
232   IF arg_present(bsfn) THEN bsfn = read2fromopa(numrst, params, 11)
233;      CALL read2(numrst,bsfd ,jpk,16)
234   IF arg_present(bsfd) THEN bsfd  = read2fromopa(numrst, params, 16)
235;#endif
236;#ifdef key_zdftke
237;          CALL read3(numrst,en,19)
238   IF arg_present(en) THEN en = read3fromopa(numrst, params, 19)
239
240
241
242   close, numrst
243   free_lun, numrst
244
245   return
246end
247
248
249
250;CDIR$ LIST
251;      SUBROUTINE dtrlec
252;CCC---------------------------------------------------------------------
253;CCC
254;CCC                       ROUTINE dtrlec
255;CCC                     ******************
256;CCC
257;CCC  Purpose :
258;CCC  --------
259;CCC     Read files for restart
260;CCC
261;CC   Method :
262;CC   -------
263;CC      Read the previous fields on the file numrst
264;CC      the first record indicates previous characteristics
265;CC      after control with the present run, we read :
266;CC      - prognostic variables on the second record
267;CC      - elliptic solver arrays
268;CC     - barotropic stream function arrays (default option)
269;CC       or free surface arrays ("key_freesurf_cstvol" defined)
270;CC      - tke arrays ("key_zdftke" defined)
271;CC      for this last three records,  the previous characteristics
272;CC      could be different with those used in the present run.
273;CC
274;CC   Input :
275;CC   ------
276;CC      common
277;CC            /comrst/          : restart parameter
278;CC            /comctl/          : parameters for the control
279;CC
280;CC   Output :
281;CC   -------
282;CC      common
283;CC            /combef/          : previous fields (before)
284;CC            /comnow/          : present fields (now)
285;CC            /combsf/          : barotropic stream function
286;CC            /comspg/          : surface pressure
287;CC            /comsol/          : diagonal preconditioned conjugate
288;CC
289;CC   Modifications :
290;CC   --------------
291;CC      original  : 91-03 ()
292;CC      additions : 92-01 (M. Imbard)
293;CC                : 92-06 correction restart file (M. Imbard)
294;CC                : 98-02 (M. Guyon) FETI method
295;CC      addition  : 98-05 (G. Roullet) free surface
296;CC----------------------------------------------------------------------
297;CC parameters and commons
298;CC ======================
299;CDIR$ NOLIST
300;#include "parameter.h"
301;#include "common.h"
302;CDIR$ LIST
303;CC----------------------------------------------------------------------
304;CC local declarations
305;CC ==================
306;      INTEGER ji, jj, jk, jl
307;      INTEGER ino0, it0, ipcg0, isor0, itke0
308;      INTEGER ino1, it1, isor1, ipcg1, itke1, idast1
309;CC----------------------------------------------------------------------
310;CC statement functions
311;CC ===================
312;CDIR$ NOLIST
313;#include "stafun.h"
314;CDIR$ LIST
315;CCC---------------------------------------------------------------------
316;CCC  OPA8, LODYC (1997)
317;CCC---------------------------------------------------------------------
318;C
319;C
320;C 0. Initialisations
321;C ------------------
322;C
323;      ino0  = no
324;      it0   = nit000
325;      ipcg0 = 0
326;      isor0 = 0
327;      itke0 = 0
328;      isor0 = nsolv-1
329;      ipcg0 = 2-nsolv
330;#ifdef key_zdftke
331;      itke0 = 1
332;#endif
333;C FETI method
334;      IF (nsolv .EQ. 3) THEN
335;          isor0=2
336;          ipcg0=2
337;      ENDIF
338;C
339;      IF(lwp) THEN
340;          WRITE(numout,*) ' '
341;          WRITE(numout,*) ' *** dtrlec:  beginning of restart'
342;          WRITE(numout,*) ' '
343;          WRITE(numout,*) ' the present run :'
344;          WRITE(numout,*) '   job number : ', no
345;          WRITE(numout,*) '   with nit000 : ', nit000
346;          WRITE(numout,*) '   with pcg option ipcg0 : ', ipcg0
347;          WRITE(numout,*) '   with sor option isor0 : ', isor0
348;          WRITE(numout,*) '   with FETI solver option ipcg0 & isor0 : ',
349;     $        ipcg0,' & ',isor0
350;          WRITE(numout,*) '   with tke option itke0 : ', itke0
351;      ENDIF
352;C
353;C 1. Read numrst
354;C --------------
355;C
356;C ... First record
357;C
358;      READ(numrst,REC=1) ino1, it1, isor1, ipcg1, itke1, idast1
359;C
360;      IF(lwp) THEN
361;          WRITE(numout,*) ' '
362;          WRITE(numout,*) ' READ numrst with '
363;          WRITE(numout,*) '   job number : ', ino1
364;          WRITE(numout,*) '   with time step it : ', it1
365;          WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
366;          WRITE(numout,*) '   with sor option isor1 : ', isor1
367;          WRITE(numout,*) '   with tke option itke1 : ', itke1
368;          WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ',
369;     $        ipcg1 + isor1
370;          WRITE(numout,*) ' '
371;      ENDIF
372;C
373;C ... Control of date
374;C
375;      IF ( (it0-it1).NE.1 .AND. abs(nrstdt).EQ.1 ) THEN
376;          IF(lwp) THEN
377;              WRITE(numout,*) ' ===>>>> : problem with nit000 for the',
378;     $            ' restart'
379;              WRITE(numout,*) ' =======                              ',
380;     $            ' ======='
381;              WRITE(numout,*) ' we stop. verify the file'
382;              WRITE(numout,*) ' or rerun with the value  0 for the'
383;              WRITE(numout,*) ' control of time parameter  nrstdt'
384;              WRITE(numout,*) ' '
385;          ENDIF
386;          STOP 'dtrlec'
387;      ENDIF
388;      IF ( nrstdt.EQ.1 ) ndate0 = idast1
389;C
390;C ... Read prognostic variables
391;C
392;      CALL read3(numrst,ub   ,2 )
393;      CALL read3(numrst,vb   ,3 )
394;      CALL read3(numrst,tb   ,5 )
395;      CALL read3(numrst,sb   ,6 )
396;      CALL read3(numrst,rotb ,7 )
397;      CALL read3(numrst,hdivb,8 )
398;      CALL read3(numrst,un   ,9 )
399;      CALL read3(numrst,vn   ,10)
400;      CALL read3(numrst,tn   ,12)
401;      CALL read3(numrst,sn   ,13)
402;      CALL read3(numrst,rotn ,14)
403;      CALL read3(numrst,hdivn,15)
404;C
405;C ... Read elliptic solver arrays
406;C
407;      CALL read2(numrst,gcx ,jpk,17)
408;      CALL read2(numrst,gcxb,jpk,18)
409;C
410;#ifdef key_freesurf_cstvol
411;C
412;C ... free surface formulation (eta)
413;C
414;      CALL read2(numrst,etab ,jpk,4 )
415;      CALL read2(numrst,etan ,jpk,4 )
416;#  else
417;C
418;C ... Rigid-lid formulation (bsf)
419;C
420;      CALL read2(numrst,bsfb ,jpk,4 )
421;      CALL read2(numrst,bsfn ,jpk,11)
422;      CALL read2(numrst,bsfd ,jpk,16)
423;#endif
424;C
425;#ifdef key_zdftke
426;C
427;C ... Read tke arrays
428;C
429;      IF(itke1.eq.1) THEN
430;          CALL read3(numrst,en,19)
431;      ELSE
432;          IF(lwp) THEN
433;              WRITE(numout,*) ' ===>>>> : the previous restart file',
434;     $            ' did''nt used  tke scheme'
435;              WRITE(numout,*) ' =======                ======='
436;          ENDIF
437;          nrstdt=2
438;      ENDIF
439;#endif
440;C
441;C
442;      RETURN
443;      END
Note: See TracBrowser for help on using the repository browser.