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

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

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