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

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

improvements/corrections of some *.pro headers
  • 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      print, '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      print, 'The size of the file is not the expected one!'
168      print, 'Check your file or the values of ibloc, jpiglo,'
169      print, '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.