1 |
C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_init.F,v 1.1 2001/09/13 17:43:55 adcroft Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "FLT_CPPOPTIONS.h" |
5 |
|
6 |
subroutine flt_init ( myCurrentIter,myCurrentTime, myThid ) |
7 |
|
8 |
c ================================================================== |
9 |
c SUBROUTINE flt_init |
10 |
c ================================================================== |
11 |
c |
12 |
c o This routine initializes the start/restart positions. |
13 |
c It does the following: |
14 |
c o First it checks for local files. These are supposed to be restarts |
15 |
c from a previous integration. The floats can therefore be read in |
16 |
c without any further check (because they exist on the specific tile). |
17 |
c o If no local files are available the routine assumes that this |
18 |
c is an initialization. In that case it reads a global file |
19 |
c (that has the same format as local files) and sorts those floats |
20 |
c that exist on the specific tile into the local array. |
21 |
c o At the end the float positions are written to the trajectory file |
22 |
c |
23 |
c ================================================================== |
24 |
c SUBROUTINE flt_init |
25 |
c ================================================================== |
26 |
|
27 |
#include "EEPARAMS.h" |
28 |
#include "SIZE.h" |
29 |
#include "FLT.h" |
30 |
#include "GRID.h" |
31 |
#include "PARAMS.h" |
32 |
|
33 |
c == routine arguments == |
34 |
|
35 |
c mythid - thread number for this instance of the routine. |
36 |
|
37 |
INTEGER myCurrentIter, myThid |
38 |
INTEGER ip, iG, jG |
39 |
_RL myCurrentTime |
40 |
|
41 |
c == local variables == |
42 |
|
43 |
INTEGER imax |
44 |
parameter(imax=9) |
45 |
_RL tmp(imax) |
46 |
integer jtlo,jthi,itlo,ithi |
47 |
INTEGER bi, bj, xx, yy |
48 |
_RL xlo, xhi, ylo, yhi |
49 |
|
50 |
logical globalFile |
51 |
|
52 |
c number of active record in the file (might be lower than the |
53 |
c total number of records because the tile could have contained |
54 |
c more floats at an earlier restart |
55 |
_RL npart_read, npart_dist |
56 |
|
57 |
character*(max_len_mbuf) msgbuf |
58 |
INTEGER K, I, J, IL, iUnit |
59 |
INTEGER errIO |
60 |
INTEGER IFNBLNK |
61 |
EXTERNAL IFNBLNK |
62 |
INTEGER ILNBLNK |
63 |
EXTERNAL ILNBLNK |
64 |
CHARACTER*(MAX_LEN_PREC) record |
65 |
|
66 |
namelist /flt_nml/ flt_int_traj, flt_int_prof, flt_noise |
67 |
& ,flt_file |
68 |
|
69 |
c == end of interface == |
70 |
|
71 |
_BEGIN_MASTER(mythid) |
72 |
|
73 |
c Set default values. |
74 |
flt_int_traj = 3600. |
75 |
flt_int_prof = 43200. |
76 |
flt_noise = 0.0 |
77 |
flt_file = 'float_pos' |
78 |
|
79 |
c call nml_filter( 'data.flt', scrunit1, myThid ) |
80 |
c if (scrunit1 .eq. 0) then |
81 |
c stop 'flt_init: reading namelist failed' |
82 |
c end if |
83 |
c read( scrunit1, nml = flt_nml ) |
84 |
c close( scrunit1 ) |
85 |
|
86 |
C-- Open the parameter file |
87 |
#ifdef TARGET_BGL |
88 |
OPEN(UNIT=scrUnit1,FILE='scratch1',STATUS='UNKNOWN') |
89 |
OPEN(UNIT=scrUnit2,FILE='scratch2',STATUS='UNKNOWN') |
90 |
#else |
91 |
OPEN(UNIT=scrUnit1,STATUS='SCRATCH') |
92 |
OPEN(UNIT=scrUnit2,STATUS='SCRATCH') |
93 |
#endif |
94 |
OPEN(UNIT=modelDataUnit,FILE='data.flt',STATUS='OLD', |
95 |
& IOSTAT=errIO) |
96 |
IF ( errIO .LT. 0 ) THEN |
97 |
WRITE(msgBuf,'(A)') |
98 |
& 'S/R FLT_INIT' |
99 |
CALL PRINT_ERROR( msgBuf , 1) |
100 |
WRITE(msgBuf,'(A)') |
101 |
& 'Unable to open flt parameter' |
102 |
CALL PRINT_ERROR( msgBuf , 1) |
103 |
WRITE(msgBuf,'(A)') |
104 |
& 'file "data.flt"' |
105 |
CALL PRINT_ERROR( msgBuf , 1) |
106 |
CALL MODELDATA_EXAMPLE( myThid ) |
107 |
STOP 'ABNORMAL END: S/R FLT_INIT' |
108 |
ENDIF |
109 |
|
110 |
DO WHILE ( .TRUE. ) |
111 |
READ(modelDataUnit,FMT='(A)',END=1001) RECORD |
112 |
IL = MAX(ILNBLNK(RECORD),1) |
113 |
IF ( RECORD(1:1) .NE. commentCharacter ) |
114 |
& WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL) |
115 |
WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL) |
116 |
ENDDO |
117 |
1001 CONTINUE |
118 |
CLOSE(modelDataUnit) |
119 |
|
120 |
C-- Report contents of model parameter file |
121 |
WRITE(msgBuf,'(A)') |
122 |
&'// =======================================================' |
123 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
124 |
& SQUEEZE_RIGHT , 1) |
125 |
WRITE(msgBuf,'(A)') '// Float parameter file "data.flt"' |
126 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
127 |
& SQUEEZE_RIGHT , 1) |
128 |
WRITE(msgBuf,'(A)') |
129 |
&'// =======================================================' |
130 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
131 |
& SQUEEZE_RIGHT , 1) |
132 |
iUnit = scrUnit2 |
133 |
REWIND(iUnit) |
134 |
DO WHILE ( .TRUE. ) |
135 |
READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD |
136 |
IL = MAX(ILNBLNK(RECORD),1) |
137 |
WRITE(msgBuf,'(A,A)') '>',RECORD(:IL) |
138 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
139 |
& SQUEEZE_RIGHT , 1) |
140 |
ENDDO |
141 |
2001 CONTINUE |
142 |
CLOSE(iUnit) |
143 |
WRITE(msgBuf,'(A)') ' ' |
144 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
145 |
& SQUEEZE_RIGHT , 1) |
146 |
|
147 |
_END_MASTER( mythid ) |
148 |
|
149 |
C-- Read settings from model parameter file "data.flt". |
150 |
iUnit = scrUnit1 |
151 |
REWIND(iUnit) |
152 |
|
153 |
READ(UNIT=iUnit,NML=FLT_NML) !,IOSTAT=errIO) |
154 |
IF ( errIO .LT. 0 ) THEN |
155 |
WRITE(msgBuf,'(A)') |
156 |
& 'S/R FLT_INIT' |
157 |
CALL PRINT_ERROR( msgBuf , 1) |
158 |
WRITE(msgBuf,'(A)') |
159 |
& 'Error reading float parameter file ' |
160 |
CALL PRINT_ERROR( msgBuf , 1) |
161 |
WRITE(msgBuf,'(A)') |
162 |
& 'parameter file "data.flt"' |
163 |
CALL PRINT_ERROR( msgBuf , 1) |
164 |
WRITE(msgBuf,'(A)') |
165 |
& 'Problem in namelist FLT_NML' |
166 |
CALL PRINT_ERROR( msgBuf , 1) |
167 |
CALL MODELDATA_EXAMPLE( myThid ) |
168 |
STOP 'ABNORMAL END: S/R FLT_INIT' |
169 |
ENDIF |
170 |
|
171 |
c do some checks |
172 |
|
173 |
IF ( useFLT .AND. useOBCS ) THEN |
174 |
WRITE(msgBuf,'(A,A)') |
175 |
& 'S/R FLT_INIT: Integrating floats is currently not possible', |
176 |
& 'in combination with open boundaries.' |
177 |
CALL PRINT_ERROR( msgBuf , myThid) |
178 |
STOP 'ABNORMAL END: S/R FLT_INIT' |
179 |
ENDIF |
180 |
|
181 |
_BARRIER |
182 |
|
183 |
c This might be faster, since the assignment is only done once. |
184 |
jtlo = mybylo(mythid) |
185 |
jthi = mybyhi(mythid) |
186 |
itlo = mybxlo(mythid) |
187 |
ithi = mybxhi(mythid) |
188 |
|
189 |
|
190 |
do bj = jtlo,jthi |
191 |
do bi = itlo,ithi |
192 |
c |
193 |
c (1) read actual number floats from file |
194 |
call mdsreadvector_flt(flt_file,globalFile,64,'RL', |
195 |
& imax,tmp,bi,bj,1,mythid) |
196 |
npart_read = tmp(1) |
197 |
max_npart = tmp(6) |
198 |
|
199 |
if (globalFile) then |
200 |
npart_tile(bi,bj) = 0 |
201 |
else |
202 |
npart_tile(bi,bj) = INT(npart_read) |
203 |
endif |
204 |
|
205 |
do ip=1,INT(npart_read) |
206 |
|
207 |
call mdsreadvector_flt(flt_file,globalFile,64,'RL', |
208 |
& imax,tmp,bi,bj,ip+1,mythid) |
209 |
|
210 |
if (globalFile) then |
211 |
c |
212 |
c check if floats are existing on tile. If not, set to zero |
213 |
c use southern/western side for axis information |
214 |
c |
215 |
|
216 |
c note: The possible area for a float has to extended to the |
217 |
c space "between" two T points, i.e. xc(sNx) of one tile |
218 |
c and xc(1) of the neighboring tile. This cannot be solved |
219 |
c by simply using xc(sNx+1) or xc(0) because periodicity |
220 |
c could imply wrong values |
221 |
c |
222 |
iG = myXGlobalLo + (bi-1)*sNx |
223 |
jG = myYGlobalLo + (bj-1)*sNy |
224 |
|
225 |
xlo = xc(1, 1, bi,bj) - delX(iG)*0.5 |
226 |
xhi = xc(sNx,1,bi,bj) + delX(iG+sNx-1)*0.5 |
227 |
ylo = yc(1, 1, bi,bj) - delY(jG)*0.5 |
228 |
yhi = yc(1,sNy,bi,bj) + delY(jG+sNy-1)*0.5 |
229 |
|
230 |
if (tmp(3) .ge. xlo .and. tmp(3) .le. xhi .and. |
231 |
& tmp(4) .ge. ylo .and. tmp(4) .le. yhi) then |
232 |
|
233 |
npart_tile(bi,bj) = npart_tile(bi,bj) + 1 |
234 |
if (npart_tile(bi,bj) .gt. max_npart_tile) |
235 |
& stop ' max_npart_tile too low. stop in flt_init' |
236 |
|
237 |
npart(npart_tile(bi,bj),bi,bj) = tmp(1) |
238 |
tstart(npart_tile(bi,bj),bi,bj) = tmp(2) |
239 |
xpart(npart_tile(bi,bj),bi,bj) = tmp(3) |
240 |
ypart(npart_tile(bi,bj),bi,bj) = tmp(4) |
241 |
kpart(npart_tile(bi,bj),bi,bj) = tmp(5) |
242 |
kfloat(npart_tile(bi,bj),bi,bj) = tmp(6) |
243 |
iup(npart_tile(bi,bj),bi,bj) = tmp(7) |
244 |
itop(npart_tile(bi,bj),bi,bj) = tmp(8) |
245 |
tend(npart_tile(bi,bj),bi,bj) = tmp(9) |
246 |
endif |
247 |
|
248 |
c else |
249 |
|
250 |
c npart(ip,bi,bj) = tmp(1) |
251 |
c tstart(ip,bi,bj) = tmp(2) |
252 |
c xpart(ip,bi,bj) = tmp(3) |
253 |
c ypart(ip,bi,bj) = tmp(4) |
254 |
c kpart(ip,bi,bj) = tmp(5) |
255 |
c kfloat(ip,bi,bj) = tmp(6) |
256 |
c iup(ip,bi,bj) = tmp(7) |
257 |
c itop(ip,bi,bj) = tmp(8) |
258 |
c tend(ip,bi,bj) = tmp(9) |
259 |
|
260 |
endif |
261 |
|
262 |
|
263 |
enddo |
264 |
|
265 |
_BARRIER |
266 |
_BEGIN_MASTER( myThid ) |
267 |
npart_dist = DBLE(npart_tile(bi,bj)) |
268 |
_GLOBAL_SUM_R8( npart_dist, myThid ) |
269 |
if (.not. globalFile) _GLOBAL_SUM_R8(npart_read,myThid) |
270 |
if (myProcId .eq. 0) then |
271 |
write(errormessageunit,*) ' max _npart: ',max_npart |
272 |
write(errormessageunit,*) 'sum npart_read: ',npart_read |
273 |
write(errormessageunit,*) 'sum npart_tile: ',npart_dist |
274 |
endif |
275 |
_END_MASTER( myThid ) |
276 |
_BARRIER |
277 |
|
278 |
enddo |
279 |
enddo |
280 |
|
281 |
return |
282 |
end |
283 |
|