| 3 |
michaesp |
1 |
PROGRAM reformat
|
|
|
2 |
|
|
|
3 |
c ***********************************************************************
|
|
|
4 |
c * Change format of a trajectory file *
|
|
|
5 |
c * Michael Sprenger / Spring, summer 2010 *
|
|
|
6 |
c ***********************************************************************
|
|
|
7 |
|
|
|
8 |
implicit none
|
|
|
9 |
|
|
|
10 |
c ----------------------------------------------------------------------
|
|
|
11 |
c Declaration of variables
|
|
|
12 |
c ----------------------------------------------------------------------
|
|
|
13 |
|
|
|
14 |
c Mode & Special parameters depending on mode
|
|
|
15 |
character*80 mode
|
|
|
16 |
real clon,clat,radius
|
|
|
17 |
|
|
|
18 |
c Input and output files
|
|
|
19 |
character*80 inpfile ! Input filename
|
|
|
20 |
character*80 outfile ! Output filename
|
|
|
21 |
|
|
|
22 |
c Trajectories
|
|
|
23 |
integer ntra ! Number of trajectories
|
|
|
24 |
integer ntim ! Number of times
|
|
|
25 |
integer ncol ! Number of columns
|
|
|
26 |
real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)
|
|
|
27 |
real,allocatable, dimension (:,:) :: proxy ! Footprint
|
|
|
28 |
character*80 vars(100) ! Variable names
|
|
|
29 |
integer refdate(6) ! Reference date
|
|
|
30 |
real,allocatable, dimension (:) :: num ! Output number
|
|
|
31 |
|
|
|
32 |
c Auxiliary variables
|
|
|
33 |
integer inpmode
|
|
|
34 |
integer stat
|
|
|
35 |
integer fid
|
|
|
36 |
integer i,j,k
|
|
|
37 |
real dist
|
|
|
38 |
real lon0,lat0,lon1,lat1
|
|
|
39 |
|
|
|
40 |
c Externals
|
|
|
41 |
real sdis
|
|
|
42 |
external sdis
|
|
|
43 |
|
|
|
44 |
c ----------------------------------------------------------------------
|
|
|
45 |
c Preparations
|
|
|
46 |
c ----------------------------------------------------------------------
|
|
|
47 |
|
|
|
48 |
c Read parameters
|
|
|
49 |
open(10,file='traj2num.param')
|
|
|
50 |
read(10,*) inpfile
|
|
|
51 |
read(10,*) outfile
|
|
|
52 |
read(10,*) ntra,ntim,ncol
|
|
|
53 |
read(10,*) mode
|
|
|
54 |
if ( mode.eq.'proxy' ) then
|
|
|
55 |
read(10,*) clon
|
|
|
56 |
read(10,*) clat
|
|
|
57 |
read(10,*) radius
|
|
|
58 |
endif
|
|
|
59 |
close(10)
|
|
|
60 |
|
|
|
61 |
c Check that a valid mode is selected
|
|
|
62 |
if ( mode.eq.'boost' ) goto 10
|
|
|
63 |
if ( mode.eq.'proxy' ) goto 10
|
|
|
64 |
print*,' Unknown mode ',trim(mode)
|
|
|
65 |
stop
|
|
|
66 |
10 continue
|
|
|
67 |
|
|
|
68 |
c Determine the formats
|
|
|
69 |
call mode_tra(inpmode,inpfile)
|
|
|
70 |
if (inpmode.eq.-1) inpmode=1
|
|
|
71 |
|
|
|
72 |
c Allocate memory
|
|
|
73 |
allocate(tra(ntra,ntim,ncol),stat=stat)
|
|
|
74 |
if (stat.ne.0) print*,'*** error allocating array tra ***'
|
|
|
75 |
allocate(num(ntra),stat=stat)
|
|
|
76 |
if (stat.ne.0) print*,'*** error allocating array num ***'
|
|
|
77 |
if ( mode.eq. 'proxy' ) then
|
|
|
78 |
allocate(proxy(ntra,ncol+1),stat=stat)
|
|
|
79 |
if (stat.ne.0) print*,'*** error allocating array proxy ***'
|
|
|
80 |
endif
|
|
|
81 |
|
|
|
82 |
c Read inpufile
|
|
|
83 |
call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
|
|
|
84 |
call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
|
|
|
85 |
call close_tra(fid,inpmode)
|
|
|
86 |
|
|
|
87 |
c ----------------------------------------------------------------------
|
|
|
88 |
c Mode = 'boost': get the maximum distance traveled in one time step
|
|
|
89 |
c ----------------------------------------------------------------------
|
|
|
90 |
|
|
|
91 |
if ( mode.ne.'boost') goto 100
|
|
|
92 |
|
|
|
93 |
do i=1,ntra
|
|
|
94 |
|
|
|
95 |
num(i) = 0.
|
|
|
96 |
|
|
|
97 |
do j=2,ntim
|
|
|
98 |
|
|
|
99 |
c Get spherical distance between data points
|
|
|
100 |
lon0 = tra(i,j-1,2)
|
|
|
101 |
lat0 = tra(i,j-1,3)
|
|
|
102 |
lon1 = tra(i,j ,2)
|
|
|
103 |
lat1 = tra(i,j ,3)
|
|
|
104 |
dist = sdis( lon1,lat1,lon0,lat0 )
|
|
|
105 |
|
|
|
106 |
if ( dist.gt.num(i) ) num(i) = dist
|
|
|
107 |
|
|
|
108 |
enddo
|
|
|
109 |
|
|
|
110 |
enddo
|
|
|
111 |
|
|
|
112 |
100 continue
|
|
|
113 |
|
|
|
114 |
c ----------------------------------------------------------------------
|
|
|
115 |
c Mode = 'proxy': get nearest distance to a lat/lon point
|
|
|
116 |
c ----------------------------------------------------------------------
|
|
|
117 |
|
|
|
118 |
if ( mode.ne.'proxy') goto 200
|
|
|
119 |
|
|
|
120 |
call get_proxy (proxy,clon,clat,tra,ntra,ntim,ncol,radius)
|
|
|
121 |
do i=1,ntra
|
|
|
122 |
num(i) = proxy(i,ncol+1)
|
|
|
123 |
enddo
|
|
|
124 |
|
|
|
125 |
200 continue
|
|
|
126 |
|
|
|
127 |
c ----------------------------------------------------------------------
|
|
|
128 |
c Write output file
|
|
|
129 |
c ----------------------------------------------------------------------
|
|
|
130 |
|
|
|
131 |
open(10,file=outfile)
|
|
|
132 |
do i=1,ntra
|
|
|
133 |
write(10,*) num(i)
|
|
|
134 |
enddo
|
|
|
135 |
close(10)
|
|
|
136 |
|
|
|
137 |
end
|
|
|
138 |
|
|
|
139 |
c ***********************************************************************
|
|
|
140 |
c * SUBROUTINES *
|
|
|
141 |
c ***********************************************************************
|
|
|
142 |
|
|
|
143 |
c -----------------------------------------------------------------------
|
|
|
144 |
c Spherical distance between lat/lon points
|
|
|
145 |
c -----------------------------------------------------------------------
|
|
|
146 |
|
|
|
147 |
real function sdis(xp,yp,xq,yq)
|
|
|
148 |
c
|
|
|
149 |
c calculates spherical distance (in km) between two points given
|
|
|
150 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
151 |
c
|
|
|
152 |
real re
|
|
|
153 |
parameter (re=6370.)
|
|
|
154 |
real pi180
|
|
|
155 |
parameter (pi180=3.14159/180.)
|
|
|
156 |
real xp,yp,xq,yq,arg
|
|
|
157 |
|
|
|
158 |
arg=sin(pi180*yp)*sin(pi180*yq)+
|
|
|
159 |
> cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
|
|
|
160 |
if (arg.lt.-1.) arg=-1.
|
|
|
161 |
if (arg.gt.1.) arg=1.
|
|
|
162 |
|
|
|
163 |
sdis=re*acos(arg)
|
|
|
164 |
|
|
|
165 |
end
|
|
|
166 |
|
|
|
167 |
c -----------------------------------------------------------------------
|
|
|
168 |
c Nearest distance to a lat/lon point
|
|
|
169 |
c -----------------------------------------------------------------------
|
|
|
170 |
|
|
|
171 |
subroutine get_proxy (proxy,clon,clat,tra,ntra,ntim,ncol,radius)
|
|
|
172 |
|
|
|
173 |
c Given a set of trajectories tra(ntra,ntim,ncol), find the closest point
|
|
|
174 |
c to clon/clat and return the footprint of this trajectory at this
|
|
|
175 |
c point: proxy(ntra,ncol+1); the parameter <ncol+1> of the trajectory will
|
|
|
176 |
c become the minimum distance.
|
|
|
177 |
|
|
|
178 |
implicit none
|
|
|
179 |
|
|
|
180 |
c Declaration of parameters
|
|
|
181 |
integer ntra,ntim,ncol
|
|
|
182 |
real tra(ntra,ntim,ncol)
|
|
|
183 |
real proxy(ntra,ncol+1)
|
|
|
184 |
real radius
|
|
|
185 |
real clon,clat
|
|
|
186 |
|
|
|
187 |
c Flag for tests
|
|
|
188 |
integer test
|
|
|
189 |
parameter (test=0)
|
|
|
190 |
character*80 testfile
|
|
|
191 |
parameter (testfile='TEST.nc')
|
|
|
192 |
|
|
|
193 |
c Number of grid points for the radius mask
|
|
|
194 |
integer nmax
|
|
|
195 |
parameter (nmax = 400)
|
|
|
196 |
real degkm
|
|
|
197 |
parameter (degkm = 111.)
|
|
|
198 |
|
|
|
199 |
c Number of binary search steps
|
|
|
200 |
integer nbin
|
|
|
201 |
parameter (nbin=10)
|
|
|
202 |
|
|
|
203 |
c Auxiliry variables
|
|
|
204 |
integer i,j,k
|
|
|
205 |
real lon,lat,rlon,rlat
|
|
|
206 |
real dist
|
|
|
207 |
character*80 varname
|
|
|
208 |
integer rnx,rny
|
|
|
209 |
real rmask (nmax,nmax)
|
|
|
210 |
real rcount(nmax,nmax)
|
|
|
211 |
real rxmin,rymin,rdx,rdy
|
|
|
212 |
integer indx,indy
|
|
|
213 |
integer isnear
|
|
|
214 |
real near,near0,near1,nearm
|
|
|
215 |
integer j0,j1
|
|
|
216 |
real r0,r1,rm,frac
|
|
|
217 |
real rlon0,rlon1,rlonm,rlat0,rlat1,rlatm
|
|
|
218 |
|
|
|
219 |
c Externals
|
|
|
220 |
real sdis
|
|
|
221 |
external sdis
|
|
|
222 |
|
|
|
223 |
c Transform into clon/clat centered system
|
|
|
224 |
do i=1,ntra
|
|
|
225 |
do j=1,ntim
|
|
|
226 |
lon = tra(i,j,2)
|
|
|
227 |
lat = tra(i,j,3)
|
|
|
228 |
call getenvir_f (clon,clat,0.,rlon,rlat,lon,lat,1)
|
|
|
229 |
tra(i,j,2) = rlon
|
|
|
230 |
tra(i,j,3) = rlat
|
|
|
231 |
enddo
|
|
|
232 |
enddo
|
|
|
233 |
|
|
|
234 |
c Init the radius mask - define the grid resolution
|
|
|
235 |
rdx = 2.5 * radius / degkm / real(nmax)
|
|
|
236 |
rdy = 2.5 * radius / degkm / real(nmax)
|
|
|
237 |
rnx = nmax
|
|
|
238 |
rny = nmax
|
|
|
239 |
rxmin = -real(rnx/2)*rdx
|
|
|
240 |
rymin = -real(rny/2)*rdy
|
|
|
241 |
do i=1,rnx
|
|
|
242 |
do j=1,rny
|
|
|
243 |
rlon = rxmin + real(i-1) * rdx
|
|
|
244 |
rlat = rymin + real(j-1) * rdy
|
|
|
245 |
dist = sdis(rlon,rlat,0.,0.)
|
|
|
246 |
if ( dist.lt.radius ) then
|
|
|
247 |
rmask(i,j) = dist
|
|
|
248 |
else
|
|
|
249 |
rmask(i,j) = radius
|
|
|
250 |
endif
|
|
|
251 |
enddo
|
|
|
252 |
enddo
|
|
|
253 |
|
|
|
254 |
c Init counter for test
|
|
|
255 |
if ( test.eq.1 ) then
|
|
|
256 |
do i=1,rnx
|
|
|
257 |
do j=1,rny
|
|
|
258 |
rcount(i,j) = 0.
|
|
|
259 |
enddo
|
|
|
260 |
enddo
|
|
|
261 |
endif
|
|
|
262 |
|
|
|
263 |
|
|
|
264 |
c Loop over all trajectories
|
|
|
265 |
do i=1,ntra
|
|
|
266 |
|
|
|
267 |
c Decide whether trajectory comes close to point
|
|
|
268 |
isnear = 0
|
|
|
269 |
near = radius
|
|
|
270 |
do j=1,ntim
|
|
|
271 |
indx = nint( ( tra(i,j,2) - rxmin ) / rdx + 1. )
|
|
|
272 |
indy = nint( ( tra(i,j,3) - rymin ) / rdy + 1. )
|
|
|
273 |
if ( (indx.ge.1).and.(indx.le.rnx).and.
|
|
|
274 |
> (indy.ge.1).and.(indy.le.rny) )
|
|
|
275 |
> then
|
|
|
276 |
if (rmask(indx,indy).lt.near) then
|
|
|
277 |
near = rmask(indx,indy)
|
|
|
278 |
isnear = j
|
|
|
279 |
endif
|
|
|
280 |
endif
|
|
|
281 |
enddo
|
|
|
282 |
|
|
|
283 |
c No close point was found - go to next trajectory
|
|
|
284 |
do k=1,ncol
|
|
|
285 |
proxy(i,k) = tra(i,1,k)
|
|
|
286 |
enddo
|
|
|
287 |
proxy(i,ncol+1) = radius
|
|
|
288 |
if ( isnear.eq.0 ) goto 310
|
|
|
289 |
|
|
|
290 |
c Get the exact position and time with binary search
|
|
|
291 |
j0 = isnear - 1
|
|
|
292 |
if ( j0.eq.0 ) j0 = 1
|
|
|
293 |
rlon0 = tra(i,j0,2)
|
|
|
294 |
rlat0 = tra(i,j0,3)
|
|
|
295 |
near0 = sdis(rlon0,rlat0,0.,0.)
|
|
|
296 |
r0 = real(j0)
|
|
|
297 |
|
|
|
298 |
j1 = isnear + 1
|
|
|
299 |
if ( j1.gt.ntim ) j1 = ntim
|
|
|
300 |
rlon1 = tra(i,j1,2)
|
|
|
301 |
rlat1 = tra(i,j1,3)
|
|
|
302 |
near1 = sdis(tra(i,j1,2),tra(i,j1,3),0.,0.)
|
|
|
303 |
r1 = real(j1)
|
|
|
304 |
|
|
|
305 |
do k=1,nbin
|
|
|
306 |
rm = 0.5 * ( r0 + r1 )
|
|
|
307 |
rlatm = 0.5 * ( rlat0 + rlat1 )
|
|
|
308 |
rlonm = 0.5 * ( rlon0 + rlon1 )
|
|
|
309 |
nearm = sdis(rlonm,rlatm,0.,0.)
|
|
|
310 |
if ( near0.lt.near1 ) then
|
|
|
311 |
r1 = rm
|
|
|
312 |
rlon1 = rlonm
|
|
|
313 |
rlat1 = rlatm
|
|
|
314 |
near1 = nearm
|
|
|
315 |
else
|
|
|
316 |
r0 = rm
|
|
|
317 |
rlon0 = rlonm
|
|
|
318 |
rlat0 = rlatm
|
|
|
319 |
near0 = nearm
|
|
|
320 |
endif
|
|
|
321 |
enddo
|
|
|
322 |
|
|
|
323 |
c Now get the final distance and position
|
|
|
324 |
proxy(i,ncol+1) = nearm
|
|
|
325 |
if ( test.eq.1 ) then
|
|
|
326 |
indx = nint( ( rlonm - rxmin ) / rdx + 1. )
|
|
|
327 |
indy = nint( ( rlatm - rymin ) / rdy + 1. )
|
|
|
328 |
if ( (indx.ge.1).and.(indx.le.rnx).and.
|
|
|
329 |
> (indy.ge.1).and.(indy.le.rny) )
|
|
|
330 |
> then
|
|
|
331 |
rcount(indx,indy) = rcount(indx,indy) + 1.
|
|
|
332 |
endif
|
|
|
333 |
endif
|
|
|
334 |
|
|
|
335 |
c Get all values at exactly this position
|
|
|
336 |
j0 = int(rm)
|
|
|
337 |
frac = rm - real(j0)
|
|
|
338 |
j1 = j0 + 1
|
|
|
339 |
if (j1.gt.ntim ) j1 = ntim
|
|
|
340 |
do k=1,ncol
|
|
|
341 |
proxy(i,k) = (1.-frac)*tra(i,j0,k)+frac*tra(i,j1,k)
|
|
|
342 |
enddo
|
|
|
343 |
|
|
|
344 |
c Next trajectory
|
|
|
345 |
310 continue
|
|
|
346 |
|
|
|
347 |
enddo
|
|
|
348 |
|
|
|
349 |
c Write radius mask to netCDF file for tests
|
|
|
350 |
if ( test.eq.1 ) then
|
|
|
351 |
varname = 'MASK'
|
|
|
352 |
call writecdf2D(testfile,varname,rmask,0.,
|
|
|
353 |
> rdx,rdy,rxmin,rymin,rnx,rny,1,1)
|
|
|
354 |
varname = 'COUNT'
|
|
|
355 |
call writecdf2D(testfile,varname,rcount,0.,
|
|
|
356 |
> rdx,rdy,rxmin,rymin,rnx,rny,0,1)
|
|
|
357 |
endif
|
|
|
358 |
|
|
|
359 |
end
|
|
|
360 |
|
|
|
361 |
|
|
|
362 |
|
|
|
363 |
c --------------------------------------------------------------------------------
|
|
|
364 |
c Forward coordinate transformation (True lon/lat -> Rotated lon/lat)
|
|
|
365 |
c --------------------------------------------------------------------------------
|
|
|
366 |
|
|
|
367 |
SUBROUTINE getenvir_f (clon,clat,rotation,
|
|
|
368 |
> rlon,rlat,lon,lat,n)
|
|
|
369 |
|
|
|
370 |
implicit none
|
|
|
371 |
|
|
|
372 |
c Declaration of input and output parameters
|
|
|
373 |
integer n
|
|
|
374 |
real clon,clat,rotation
|
|
|
375 |
real lon(n), lat(n)
|
|
|
376 |
real rlon(n),rlat(n)
|
|
|
377 |
|
|
|
378 |
c Auxiliary variables
|
|
|
379 |
real pollon,pollat
|
|
|
380 |
integer i
|
|
|
381 |
real rlon1(n),rlat1(n)
|
|
|
382 |
|
|
|
383 |
c Externals
|
|
|
384 |
real lmtolms,phtophs
|
|
|
385 |
external lmtolms,phtophs
|
|
|
386 |
|
|
|
387 |
c First rotation
|
|
|
388 |
pollon=clon-180.
|
|
|
389 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
390 |
pollat=90.-clat
|
|
|
391 |
do i=1,n
|
|
|
392 |
|
|
|
393 |
c First rotation
|
|
|
394 |
pollon=clon-180.
|
|
|
395 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
396 |
pollat=90.-clat
|
|
|
397 |
rlon1(i)=lmtolms(lat(i),lon(i),pollat,pollon)
|
|
|
398 |
rlat1(i)=phtophs(lat(i),lon(i),pollat,pollon)
|
|
|
399 |
|
|
|
400 |
c Second coordinate transformation
|
|
|
401 |
pollon=-180.
|
|
|
402 |
pollat=90.+rotation
|
|
|
403 |
rlon(i)=90.+lmtolms(rlat1(i),rlon1(i)-90.,pollat,pollon)
|
|
|
404 |
rlat(i)=phtophs(rlat1(i),rlon1(i)-90.,pollat,pollon)
|
|
|
405 |
|
|
|
406 |
enddo
|
|
|
407 |
|
|
|
408 |
END
|
|
|
409 |
|
|
|
410 |
c --------------------------------------------------------------------------------
|
|
|
411 |
c Backward coordinate transformation (Rotated lon/lat -> True lon/lat)
|
|
|
412 |
c --------------------------------------------------------------------------------
|
|
|
413 |
|
|
|
414 |
SUBROUTINE getenvir_b (clon,clat,rotation,
|
|
|
415 |
> lon,lat,rlon,rlat,n)
|
|
|
416 |
|
|
|
417 |
implicit none
|
|
|
418 |
|
|
|
419 |
c Declaration of input and output parameters
|
|
|
420 |
integer n
|
|
|
421 |
real clon,clat,rotation
|
|
|
422 |
real lon(n), lat(n)
|
|
|
423 |
real rlon(n),rlat(n)
|
|
|
424 |
|
|
|
425 |
c Auxiliary variables
|
|
|
426 |
real pollon,pollat
|
|
|
427 |
integer i
|
|
|
428 |
real rlon1(n),rlat1(n)
|
|
|
429 |
|
|
|
430 |
c Externals
|
|
|
431 |
real lmstolm,phstoph
|
|
|
432 |
external lmstolm,phstoph
|
|
|
433 |
|
|
|
434 |
c First coordinate transformation (make the local coordinate system parallel to equator)
|
|
|
435 |
pollon=-180.
|
|
|
436 |
pollat=90.+rotation
|
|
|
437 |
do i=1,n
|
|
|
438 |
rlon1(i)=90.+lmstolm(rlat(i),rlon(i)-90.,pollat,pollon)
|
|
|
439 |
rlat1(i)=phstoph(rlat(i),rlon(i)-90.,pollat,pollon)
|
|
|
440 |
enddo
|
|
|
441 |
|
|
|
442 |
c Second coordinate transformation (make the local coordinate system parallel to equator)
|
|
|
443 |
pollon=clon-180.
|
|
|
444 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
445 |
pollat=90.-clat
|
|
|
446 |
do i=1,n
|
|
|
447 |
lon(i)=lmstolm(rlat1(i),rlon1(i),pollat,pollon)
|
|
|
448 |
lat(i)=phstoph(rlat1(i),rlon1(i),pollat,pollon)
|
|
|
449 |
enddo
|
|
|
450 |
|
|
|
451 |
END
|
|
|
452 |
|
|
|
453 |
c --------------------------------------------------------------------------------
|
|
|
454 |
c Transformation routine: LMSTOLM and PHSTOPH from library gm2em
|
|
|
455 |
c --------------------------------------------------------------------------------
|
|
|
456 |
|
|
|
457 |
REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
458 |
C
|
|
|
459 |
C**** LMSTOLM - FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
|
|
|
460 |
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
461 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
462 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
463 |
C** AUFRUF : LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
464 |
C** ENTRIES : KEINE
|
|
|
465 |
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
|
|
|
466 |
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
467 |
C** IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
468 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
469 |
C** VERSIONS-
|
|
|
470 |
C** DATUM : 03.05.90
|
|
|
471 |
C**
|
|
|
472 |
C** EXTERNALS: KEINE
|
|
|
473 |
C** EINGABE-
|
|
|
474 |
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.
|
|
|
475 |
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
|
|
|
476 |
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS
|
|
|
477 |
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS
|
|
|
478 |
C** AUSGABE-
|
|
|
479 |
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
|
|
|
480 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
481 |
C**
|
|
|
482 |
C** COMMON-
|
|
|
483 |
C** BLOECKE : KEINE
|
|
|
484 |
C**
|
|
|
485 |
C** FEHLERBE-
|
|
|
486 |
C** HANDLUNG : KEINE
|
|
|
487 |
C** VERFASSER: D.MAJEWSKI
|
|
|
488 |
|
|
|
489 |
REAL LAMS,PHIS,POLPHI,POLLAM
|
|
|
490 |
|
|
|
491 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
492 |
|
|
|
493 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
494 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
495 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
496 |
ZPHIS = ZPIR18*PHIS
|
|
|
497 |
ZLAMS = LAMS
|
|
|
498 |
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
|
|
|
499 |
ZLAMS = ZPIR18*ZLAMS
|
|
|
500 |
|
|
|
501 |
ZARG1 = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +
|
|
|
502 |
1 ZCOSPOL* SIN(ZPHIS)) -
|
|
|
503 |
2 COS(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)
|
|
|
504 |
ZARG2 = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +
|
|
|
505 |
1 ZCOSPOL* SIN(ZPHIS)) +
|
|
|
506 |
2 SIN(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)
|
|
|
507 |
IF (ABS(ZARG2).LT.1.E-30) THEN
|
|
|
508 |
IF (ABS(ZARG1).LT.1.E-30) THEN
|
|
|
509 |
LMSTOLM = 0.0
|
|
|
510 |
ELSEIF (ZARG1.GT.0.) THEN
|
|
|
511 |
LMSTOLAM = 90.0
|
|
|
512 |
ELSE
|
|
|
513 |
LMSTOLAM = -90.0
|
|
|
514 |
ENDIF
|
|
|
515 |
ELSE
|
|
|
516 |
LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
|
|
|
517 |
ENDIF
|
|
|
518 |
|
|
|
519 |
RETURN
|
|
|
520 |
END
|
|
|
521 |
|
|
|
522 |
REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
523 |
C
|
|
|
524 |
C**** PHSTOPH - FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
|
|
|
525 |
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
526 |
C**** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
527 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
528 |
C** AUFRUF : PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
529 |
C** ENTRIES : KEINE
|
|
|
530 |
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
|
|
|
531 |
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
532 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
533 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
534 |
C** VERSIONS-
|
|
|
535 |
C** DATUM : 03.05.90
|
|
|
536 |
C**
|
|
|
537 |
C** EXTERNALS: KEINE
|
|
|
538 |
C** EINGABE-
|
|
|
539 |
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.
|
|
|
540 |
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
|
|
|
541 |
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS
|
|
|
542 |
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS
|
|
|
543 |
C** AUSGABE-
|
|
|
544 |
C** PARAMETER: WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
|
|
|
545 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
546 |
C**
|
|
|
547 |
C** COMMON-
|
|
|
548 |
C** BLOECKE : KEINE
|
|
|
549 |
C**
|
|
|
550 |
C** FEHLERBE-
|
|
|
551 |
C** HANDLUNG : KEINE
|
|
|
552 |
C** VERFASSER: D.MAJEWSKI
|
|
|
553 |
|
|
|
554 |
REAL LAMS,PHIS,POLPHI,POLLAM
|
|
|
555 |
|
|
|
556 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
557 |
|
|
|
558 |
SINPOL = SIN(ZPIR18*POLPHI)
|
|
|
559 |
COSPOL = COS(ZPIR18*POLPHI)
|
|
|
560 |
ZPHIS = ZPIR18*PHIS
|
|
|
561 |
ZLAMS = LAMS
|
|
|
562 |
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
|
|
|
563 |
ZLAMS = ZPIR18*ZLAMS
|
|
|
564 |
ARG = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
|
|
|
565 |
|
|
|
566 |
PHSTOPH = ZRPI18*ASIN(ARG)
|
|
|
567 |
|
|
|
568 |
RETURN
|
|
|
569 |
END
|
|
|
570 |
|
|
|
571 |
REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
572 |
C
|
|
|
573 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
574 |
C
|
|
|
575 |
C**** LMTOLMS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
|
|
|
576 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
577 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
578 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
579 |
C** AUFRUF : LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
580 |
C** ENTRIES : KEINE
|
|
|
581 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
|
|
|
582 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
583 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
584 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
585 |
C** VERSIONS-
|
|
|
586 |
C** DATUM : 03.05.90
|
|
|
587 |
C**
|
|
|
588 |
C** EXTERNALS: KEINE
|
|
|
589 |
C** EINGABE-
|
|
|
590 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
591 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
592 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
593 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
594 |
C** AUSGABE-
|
|
|
595 |
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
|
|
|
596 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
597 |
C**
|
|
|
598 |
C** COMMON-
|
|
|
599 |
C** BLOECKE : KEINE
|
|
|
600 |
C**
|
|
|
601 |
C** FEHLERBE-
|
|
|
602 |
C** HANDLUNG : KEINE
|
|
|
603 |
C** VERFASSER: G. DE MORSIER
|
|
|
604 |
|
|
|
605 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
606 |
|
|
|
607 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
608 |
|
|
|
609 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
610 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
611 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
612 |
ZPHI = ZPIR18*PHI
|
|
|
613 |
ZLAM = LAM
|
|
|
614 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
615 |
ZLAM = ZPIR18*ZLAM
|
|
|
616 |
|
|
|
617 |
ZARG1 = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
|
|
|
618 |
ZARG2 = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
|
|
|
619 |
IF (ABS(ZARG2).LT.1.E-30) THEN
|
|
|
620 |
IF (ABS(ZARG1).LT.1.E-30) THEN
|
|
|
621 |
LMTOLMS = 0.0
|
|
|
622 |
ELSEIF (ZARG1.GT.0.) THEN
|
|
|
623 |
LMTOLMS = 90.0
|
|
|
624 |
ELSE
|
|
|
625 |
LMTOLMS = -90.0
|
|
|
626 |
ENDIF
|
|
|
627 |
ELSE
|
|
|
628 |
LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
|
|
|
629 |
ENDIF
|
|
|
630 |
|
|
|
631 |
RETURN
|
|
|
632 |
END
|
|
|
633 |
|
|
|
634 |
REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
635 |
C
|
|
|
636 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
637 |
C
|
|
|
638 |
C**** PHTOPHS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
|
|
|
639 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
640 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
641 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
642 |
C** AUFRUF : PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
643 |
C** ENTRIES : KEINE
|
|
|
644 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
|
|
|
645 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
646 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
647 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
648 |
C** VERSIONS-
|
|
|
649 |
C** DATUM : 03.05.90
|
|
|
650 |
C**
|
|
|
651 |
C** EXTERNALS: KEINE
|
|
|
652 |
C** EINGABE-
|
|
|
653 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
654 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
655 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
656 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
657 |
C** AUSGABE-
|
|
|
658 |
C** PARAMETER: ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
|
|
|
659 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
660 |
C**
|
|
|
661 |
C** COMMON-
|
|
|
662 |
C** BLOECKE : KEINE
|
|
|
663 |
C**
|
|
|
664 |
C** FEHLERBE-
|
|
|
665 |
C** HANDLUNG : KEINE
|
|
|
666 |
C** VERFASSER: G. DE MORSIER
|
|
|
667 |
|
|
|
668 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
669 |
|
|
|
670 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
671 |
|
|
|
672 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
673 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
674 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
675 |
ZPHI = ZPIR18*PHI
|
|
|
676 |
ZLAM = LAM
|
|
|
677 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
678 |
ZLAM = ZPIR18*ZLAM
|
|
|
679 |
ZARG = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
|
|
|
680 |
|
|
|
681 |
PHTOPHS = ZRPI18*ASIN(ZARG)
|
|
|
682 |
|
|
|
683 |
RETURN
|
|
|
684 |
END
|
|
|
685 |
|
|
|
686 |
c --------------------------------------------------------------------
|
|
|
687 |
c Subroutines to write the netcdf output file
|
|
|
688 |
c --------------------------------------------------------------------
|
|
|
689 |
|
|
|
690 |
subroutine writecdf2D(cdfname,
|
|
|
691 |
> varname,arr,time,
|
|
|
692 |
> dx,dy,xmin,ymin,nx,ny,
|
|
|
693 |
> crefile,crevar)
|
|
|
694 |
|
|
|
695 |
c Create and write to the netcdf file <cdfname>. The variable
|
|
|
696 |
c with name <varname> and with time <time> is written. The data
|
|
|
697 |
c are in the two-dimensional array <arr>. The list <dx,dy,xmin,
|
|
|
698 |
c ymin,nx,ny> specifies the output grid. The flags <crefile> and
|
|
|
699 |
c <crevar> determine whether the file and/or the variable should
|
|
|
700 |
c be created
|
|
|
701 |
|
|
|
702 |
IMPLICIT NONE
|
|
|
703 |
|
|
|
704 |
c Declaration of input parameters
|
|
|
705 |
character*80 cdfname,varname
|
|
|
706 |
integer nx,ny
|
|
|
707 |
real arr(nx,ny)
|
|
|
708 |
real dx,dy,xmin,ymin
|
|
|
709 |
real time
|
|
|
710 |
integer crefile,crevar
|
|
|
711 |
|
|
|
712 |
c Further variables
|
|
|
713 |
real varmin(4),varmax(4),stag(4)
|
|
|
714 |
integer ierr,cdfid,ndim,vardim(4)
|
|
|
715 |
character*80 cstname
|
|
|
716 |
real mdv
|
|
|
717 |
integer datar(14),stdate(5)
|
|
|
718 |
integer i
|
|
|
719 |
|
|
|
720 |
C Definitions
|
|
|
721 |
varmin(1)=xmin
|
|
|
722 |
varmin(2)=ymin
|
|
|
723 |
varmin(3)=1050.
|
|
|
724 |
varmax(1)=xmin+real(nx-1)*dx
|
|
|
725 |
varmax(2)=ymin+real(ny-1)*dy
|
|
|
726 |
varmax(3)=1050.
|
|
|
727 |
cstname=trim(cdfname)//'_cst'
|
|
|
728 |
ndim=4
|
|
|
729 |
vardim(1)=nx
|
|
|
730 |
vardim(2)=ny
|
|
|
731 |
vardim(3)=1
|
|
|
732 |
stag(1)=0.
|
|
|
733 |
stag(2)=0.
|
|
|
734 |
stag(3)=0.
|
|
|
735 |
mdv=-999.98999
|
|
|
736 |
|
|
|
737 |
C Create the file
|
|
|
738 |
if (crefile.eq.0) then
|
|
|
739 |
call cdfwopn(cdfname,cdfid,ierr)
|
|
|
740 |
if (ierr.ne.0) goto 906
|
|
|
741 |
else if (crefile.eq.1) then
|
|
|
742 |
call crecdf(cdfname,cdfid,
|
|
|
743 |
> varmin,varmax,ndim,cstname,ierr)
|
|
|
744 |
if (ierr.ne.0) goto 903
|
|
|
745 |
|
|
|
746 |
C Write the constants file
|
|
|
747 |
datar(1)=vardim(1)
|
|
|
748 |
datar(2)=vardim(2)
|
|
|
749 |
datar(3)=int(1000.*varmax(2))
|
|
|
750 |
datar(4)=int(1000.*varmin(1))
|
|
|
751 |
datar(5)=int(1000.*varmin(2))
|
|
|
752 |
datar(6)=int(1000.*varmax(1))
|
|
|
753 |
datar(7)=int(1000.*dx)
|
|
|
754 |
datar(8)=int(1000.*dy)
|
|
|
755 |
datar(9)=1
|
|
|
756 |
datar(10)=1
|
|
|
757 |
datar(11)=0 ! data version
|
|
|
758 |
datar(12)=0 ! cstfile version
|
|
|
759 |
datar(13)=0 ! longitude of pole
|
|
|
760 |
datar(14)=90000 ! latitude of pole
|
|
|
761 |
do i=1,5
|
|
|
762 |
stdate(i)=0
|
|
|
763 |
enddo
|
|
|
764 |
c call wricst(cstname,datar,0.,0.,0.,0.,stdate)
|
|
|
765 |
endif
|
|
|
766 |
|
|
|
767 |
c Write the data to the netcdf file, and close the file
|
|
|
768 |
if (crevar.eq.1) then
|
|
|
769 |
call putdef(cdfid,varname,ndim,mdv,
|
|
|
770 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
771 |
if (ierr.ne.0) goto 904
|
|
|
772 |
endif
|
|
|
773 |
call putdat(cdfid,varname,time,0,arr,ierr)
|
|
|
774 |
if (ierr.ne.0) goto 905
|
|
|
775 |
call clscdf(cdfid,ierr)
|
|
|
776 |
|
|
|
777 |
return
|
|
|
778 |
|
|
|
779 |
c Error handling
|
|
|
780 |
903 print*,'*** Problem to create netcdf file ***'
|
|
|
781 |
stop
|
|
|
782 |
904 print*,'*** Problem to write definition ***'
|
|
|
783 |
stop
|
|
|
784 |
905 print*,'*** Problem to write data ***'
|
|
|
785 |
stop
|
|
|
786 |
906 print*,'*** Problem to open netcdf file ***'
|
|
|
787 |
stop
|
|
|
788 |
|
|
|
789 |
END
|
|
|
790 |
|
|
|
791 |
|
|
|
792 |
|