source: trunk/SRC/ReadWrite/readoldoparestart.pro

Last change on this file was 495, checked in by pinsard, 10 years ago

fix thanks to coding rules; typo; dupe empty lines; trailing blanks

  • Property svn:keywords set to Id
File size: 11.8 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   fileparameters = 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 fileparameters.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   close, numrst
236   free_lun, numrst
237
238   return
239end
240
241;CDIR$ LIST
242;      SUBROUTINE dtrlec
243;CCC---------------------------------------------------------------------
244;CCC
245;CCC                       ROUTINE dtrlec
246;CCC                     ******************
247;CCC
248;CCC  Purpose :
249;CCC  --------
250;CCC     Read files for restart
251;CCC
252;CC   Method :
253;CC   -------
254;CC      Read the previous fields on the file numrst
255;CC      the first record indicates previous characteristics
256;CC      after control with the present run, we read :
257;CC      - prognostic variables on the second record
258;CC      - elliptic solver arrays
259;CC     - barotropic stream function arrays (default option)
260;CC       or free surface arrays ("key_freesurf_cstvol" defined)
261;CC      - tke arrays ("key_zdftke" defined)
262;CC      for this last three records,  the previous characteristics
263;CC      could be different with those used in the present run.
264;CC
265;CC   Input :
266;CC   ------
267;CC      common
268;CC            /comrst/          : restart parameter
269;CC            /comctl/          : parameters for the control
270;CC
271;CC   Output :
272;CC   -------
273;CC      common
274;CC            /combef/          : previous fields (before)
275;CC            /comnow/          : present fields (now)
276;CC            /combsf/          : barotropic stream function
277;CC            /comspg/          : surface pressure
278;CC            /comsol/          : diagonal preconditioned conjugate
279;CC
280;CC   Modifications :
281;CC   --------------
282;CC      original  : 91-03 ()
283;CC      additions : 92-01 (M. Imbard)
284;CC                : 92-06 correction restart file (M. Imbard)
285;CC                : 98-02 (M. Guyon) FETI method
286;CC      addition  : 98-05 (G. Roullet) free surface
287;CC----------------------------------------------------------------------
288;CC parameters and commons
289;CC ======================
290;CDIR$ NOLIST
291;#include "parameter.h"
292;#include "common.h"
293;CDIR$ LIST
294;CC----------------------------------------------------------------------
295;CC local declarations
296;CC ==================
297;      INTEGER ji, jj, jk, jl
298;      INTEGER ino0, it0, ipcg0, isor0, itke0
299;      INTEGER ino1, it1, isor1, ipcg1, itke1, idast1
300;CC----------------------------------------------------------------------
301;CC statement functions
302;CC ===================
303;CDIR$ NOLIST
304;#include "stafun.h"
305;CDIR$ LIST
306;CCC---------------------------------------------------------------------
307;CCC  OPA8, LODYC (1997)
308;CCC---------------------------------------------------------------------
309;C
310;C
311;C 0. Initialisations
312;C ------------------
313;C
314;      ino0  = no
315;      it0   = nit000
316;      ipcg0 = 0
317;      isor0 = 0
318;      itke0 = 0
319;      isor0 = nsolv-1
320;      ipcg0 = 2-nsolv
321;#ifdef key_zdftke
322;      itke0 = 1
323;#endif
324;C FETI method
325;      IF (nsolv .EQ. 3) THEN
326;          isor0=2
327;          ipcg0=2
328;      ENDIF
329;C
330;      IF(lwp) THEN
331;          WRITE(numout,*) ' '
332;          WRITE(numout,*) ' *** dtrlec:  beginning of restart'
333;          WRITE(numout,*) ' '
334;          WRITE(numout,*) ' the present run :'
335;          WRITE(numout,*) '   job number : ', no
336;          WRITE(numout,*) '   with nit000 : ', nit000
337;          WRITE(numout,*) '   with pcg option ipcg0 : ', ipcg0
338;          WRITE(numout,*) '   with sor option isor0 : ', isor0
339;          WRITE(numout,*) '   with FETI solver option ipcg0 & isor0 : ',
340;     $        ipcg0,' & ',isor0
341;          WRITE(numout,*) '   with tke option itke0 : ', itke0
342;      ENDIF
343;C
344;C 1. Read numrst
345;C --------------
346;C
347;C ... First record
348;C
349;      READ(numrst,REC=1) ino1, it1, isor1, ipcg1, itke1, idast1
350;C
351;      IF(lwp) THEN
352;          WRITE(numout,*) ' '
353;          WRITE(numout,*) ' READ numrst with '
354;          WRITE(numout,*) '   job number : ', ino1
355;          WRITE(numout,*) '   with time step it : ', it1
356;          WRITE(numout,*) '   with pcg option ipcg1 : ', ipcg1
357;          WRITE(numout,*) '   with sor option isor1 : ', isor1
358;          WRITE(numout,*) '   with tke option itke1 : ', itke1
359;          WRITE(numout,*) '   with FETI solver option ipcg1 + isor1 : ',
360;     $        ipcg1 + isor1
361;          WRITE(numout,*) ' '
362;      ENDIF
363;C
364;C ... Control of date
365;C
366;      IF ( (it0-it1).NE.1 .AND. abs(nrstdt).EQ.1 ) THEN
367;          IF(lwp) THEN
368;              WRITE(numout,*) ' ===>>>> : problem with nit000 for the',
369;     $            ' restart'
370;              WRITE(numout,*) ' =======                              ',
371;     $            ' ======='
372;              WRITE(numout,*) ' we stop. verify the file'
373;              WRITE(numout,*) ' or rerun with the value  0 for the'
374;              WRITE(numout,*) ' control of time parameter  nrstdt'
375;              WRITE(numout,*) ' '
376;          ENDIF
377;          STOP 'dtrlec'
378;      ENDIF
379;      IF ( nrstdt.EQ.1 ) ndate0 = idast1
380;C
381;C ... Read prognostic variables
382;C
383;      CALL read3(numrst,ub   ,2 )
384;      CALL read3(numrst,vb   ,3 )
385;      CALL read3(numrst,tb   ,5 )
386;      CALL read3(numrst,sb   ,6 )
387;      CALL read3(numrst,rotb ,7 )
388;      CALL read3(numrst,hdivb,8 )
389;      CALL read3(numrst,un   ,9 )
390;      CALL read3(numrst,vn   ,10)
391;      CALL read3(numrst,tn   ,12)
392;      CALL read3(numrst,sn   ,13)
393;      CALL read3(numrst,rotn ,14)
394;      CALL read3(numrst,hdivn,15)
395;C
396;C ... Read elliptic solver arrays
397;C
398;      CALL read2(numrst,gcx ,jpk,17)
399;      CALL read2(numrst,gcxb,jpk,18)
400;C
401;#ifdef key_freesurf_cstvol
402;C
403;C ... free surface formulation (eta)
404;C
405;      CALL read2(numrst,etab ,jpk,4 )
406;      CALL read2(numrst,etan ,jpk,4 )
407;#  else
408;C
409;C ... Rigid-lid formulation (bsf)
410;C
411;      CALL read2(numrst,bsfb ,jpk,4 )
412;      CALL read2(numrst,bsfn ,jpk,11)
413;      CALL read2(numrst,bsfd ,jpk,16)
414;#endif
415;C
416;#ifdef key_zdftke
417;C
418;C ... Read tke arrays
419;C
420;      IF(itke1.eq.1) THEN
421;          CALL read3(numrst,en,19)
422;      ELSE
423;          IF(lwp) THEN
424;              WRITE(numout,*) ' ===>>>> : the previous restart file',
425;     $            ' did''nt used  tke scheme'
426;              WRITE(numout,*) ' =======                ======='
427;          ENDIF
428;          nrstdt=2
429;      ENDIF
430;#endif
431;C
432;C
433;      RETURN
434;      END
Note: See TracBrowser for help on using the repository browser.