| 3 |
michaesp |
1 |
PROGRAM density
|
|
|
2 |
|
|
|
3 |
use netcdf
|
|
|
4 |
|
|
|
5 |
implicit none
|
|
|
6 |
|
|
|
7 |
c ---------------------------------------------------------------------
|
|
|
8 |
c Declaration of variables
|
|
|
9 |
c ---------------------------------------------------------------------
|
|
|
10 |
|
|
|
11 |
c Parameter and working arrays
|
|
|
12 |
real radius
|
|
|
13 |
character*80 runit
|
|
|
14 |
integer nx,ny
|
|
|
15 |
integer nlonlat
|
|
|
16 |
real dlonlat
|
|
|
17 |
real xmin,ymin,dx,dy
|
|
|
18 |
real clon,clat
|
|
|
19 |
integer ntime,nfield,ntra
|
|
|
20 |
character*80 inpfile
|
|
|
21 |
character*80 outfile
|
|
|
22 |
character*80 mode
|
|
|
23 |
real param
|
|
|
24 |
integer opts,npts
|
|
|
25 |
integer step
|
|
|
26 |
character*80 gridtype
|
|
|
27 |
character*80 field
|
|
|
28 |
integer crefile,crevar
|
|
|
29 |
real,allocatable, dimension (:,:) :: cnt,res,fld,area
|
|
|
30 |
real,allocatable, dimension (:) :: traj
|
|
|
31 |
real,allocatable, dimension (:) :: olon,olat,otim,ofld
|
|
|
32 |
real,allocatable, dimension (:) :: nlon,nlat,ntim,nfld
|
|
|
33 |
|
|
|
34 |
c Output format
|
|
|
35 |
character*80 outformat
|
|
|
36 |
|
|
|
37 |
c Physical and mathematical constants
|
|
|
38 |
real pi180
|
|
|
39 |
parameter (pi180=3.14159/180.)
|
|
|
40 |
real deltay
|
|
|
41 |
parameter (deltay=111.)
|
|
|
42 |
real eps
|
|
|
43 |
parameter (eps=0.001)
|
|
|
44 |
|
|
|
45 |
c Input trajectories (see iotra.f)
|
|
|
46 |
integer inpmode
|
|
|
47 |
real,allocatable, dimension (:,:,:) :: trainp
|
|
|
48 |
integer reftime(6)
|
|
|
49 |
character*80 varsinp(100)
|
|
|
50 |
integer,allocatable, dimension (:) :: sel_flag
|
|
|
51 |
character*80 sel_file
|
|
|
52 |
character*80 sel_format
|
|
|
53 |
|
|
|
54 |
c Auxiliary variables
|
|
|
55 |
character*80 cdfname,varname
|
|
|
56 |
integer i,j,k
|
|
|
57 |
integer stat
|
|
|
58 |
integer,allocatable, dimension (:,:) :: connect0
|
|
|
59 |
integer connectval0
|
|
|
60 |
integer,allocatable, dimension (:,:) :: connect1
|
|
|
61 |
integer connectval1
|
|
|
62 |
integer,allocatable, dimension (:,:) :: connect2
|
|
|
63 |
integer connectval2
|
|
|
64 |
real slat
|
|
|
65 |
integer ipre
|
|
|
66 |
real addvalue
|
|
|
67 |
real xmax,ymax
|
|
|
68 |
real ,allocatable, dimension (:) :: odist,ndist
|
|
|
69 |
real dt
|
|
|
70 |
integer fid
|
|
|
71 |
integer dynamic_grid
|
|
|
72 |
real ycen,xcen
|
|
|
73 |
integer indx,indy
|
|
|
74 |
character*80 unit
|
|
|
75 |
real pollon,pollat
|
|
|
76 |
real rlon0,rlat0,rlon,rlat
|
|
|
77 |
real lon,lat
|
|
|
78 |
real crot
|
|
|
79 |
integer count
|
|
|
80 |
character*80 longname, varunit
|
|
|
81 |
real time
|
|
|
82 |
integer ind
|
|
|
83 |
integer ifield
|
|
|
84 |
real hhmm,frac
|
|
|
85 |
integer ierr,ncID
|
|
|
86 |
|
|
|
87 |
c External functions
|
|
|
88 |
real lmstolm,lmtolms
|
|
|
89 |
real phstoph,phtophs
|
|
|
90 |
external lmstolm,lmtolms,phstoph,phtophs
|
|
|
91 |
|
|
|
92 |
real sdis
|
|
|
93 |
external sdis
|
|
|
94 |
|
|
|
95 |
c ---------------------------------------------------------------------
|
|
|
96 |
c Preparations
|
|
|
97 |
c ---------------------------------------------------------------------
|
|
|
98 |
|
|
|
99 |
c Write start message
|
|
|
100 |
print*,'========================================================='
|
|
|
101 |
print*,' *** START OF PROGRAM DENSITY ***'
|
|
|
102 |
print*
|
|
|
103 |
|
|
|
104 |
c Read input parameters
|
|
|
105 |
open(10,file='density.param')
|
|
|
106 |
read(10,*) inpfile
|
|
|
107 |
read(10,*) outfile
|
|
|
108 |
read(10,*) field
|
|
|
109 |
read(10,*) ntime,nfield,ntra
|
|
|
110 |
read(10,*) gridtype
|
|
|
111 |
if ( gridtype.eq.'latlon' ) then
|
|
|
112 |
read(10,*) nx,ny,xmin,ymin,dx,dy
|
|
|
113 |
elseif ( gridtype.eq.'rotated') then
|
|
|
114 |
read(10,*) clon,clat,nlonlat,dlonlat
|
|
|
115 |
else
|
|
|
116 |
print*,' ERROR: unsupported grid type ',trim(gridtype)
|
|
|
117 |
stop
|
|
|
118 |
endif
|
|
|
119 |
read(10,*) radius,runit
|
|
|
120 |
read(10,*) mode
|
|
|
121 |
read(10,*) param
|
|
|
122 |
read(10,*) step
|
|
|
123 |
read(10,*) sel_file
|
|
|
124 |
read(10,*) sel_format
|
|
|
125 |
read(10,*) crefile
|
|
|
126 |
read(10,*) crevar
|
|
|
127 |
close(10)
|
|
|
128 |
|
|
|
129 |
c Get the grid parameters if <crefile=0>
|
|
|
130 |
if ( crefile.eq.0 ) then
|
|
|
131 |
|
|
|
132 |
ierr = nf90_open (trim(outfile), NF90_NOWRITE , ncID)
|
|
|
133 |
|
| 39 |
michaesp |
134 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'grid' ,gridtype )
|
|
|
135 |
if ( gridtype.eq.'rotated' ) then
|
|
|
136 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon' ,clon )
|
|
|
137 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat' ,clat )
|
|
|
138 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )
|
|
|
139 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )
|
|
|
140 |
endif
|
| 3 |
michaesp |
141 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx' ,nx )
|
|
|
142 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny' ,ny )
|
|
|
143 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx' ,dx )
|
|
|
144 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy' ,dy )
|
|
|
145 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin' ,xmin )
|
|
|
146 |
ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin' ,ymin )
|
|
|
147 |
|
|
|
148 |
ierr = nf90_close(ncID)
|
|
|
149 |
|
|
|
150 |
print*,'**** GRID PARAMETERS IMPORTED ',
|
|
|
151 |
> 'FROM NETCDF FILE!!!! ****'
|
|
|
152 |
print*
|
|
|
153 |
|
|
|
154 |
endif
|
|
|
155 |
|
|
|
156 |
c Check for consistency
|
|
|
157 |
if ( (step.ne.0).and.(mode.ne.'keep') ) then
|
|
|
158 |
print*," ERROR: interpolation is only possible for all",
|
|
|
159 |
> ' time steps... Stop'
|
|
|
160 |
stop
|
|
|
161 |
endif
|
|
|
162 |
|
|
|
163 |
c Set the number of times (just code aesthetics)
|
|
|
164 |
opts=ntime
|
|
|
165 |
|
|
|
166 |
c Set grid parameters for rotated grid
|
|
|
167 |
if ( gridtype.eq.'rotated' ) then
|
|
|
168 |
nx = nlonlat
|
|
|
169 |
ny = nlonlat
|
|
|
170 |
dx = dlonlat
|
|
|
171 |
dy = dlonlat
|
|
|
172 |
xmin = - real(nlonlat-1)/2. * dx
|
|
|
173 |
xmax = + real(nlonlat-1)/2. * dx
|
|
|
174 |
ymin = - real(nlonlat-1)/2. * dy
|
|
|
175 |
ymax = + real(nlonlat-1)/2. * dy
|
|
|
176 |
endif
|
|
|
177 |
|
|
|
178 |
c Set the flag for dynamic grid adjustment
|
|
|
179 |
if ( (nx.eq.0).or.(ny.eq.0) ) then
|
|
|
180 |
dynamic_grid = 1
|
|
|
181 |
else
|
|
|
182 |
dynamic_grid = 0
|
|
|
183 |
endif
|
|
|
184 |
|
|
|
185 |
c Print status information
|
|
|
186 |
print*,'---- INPUT PARAMETERS -----------------------------------'
|
|
|
187 |
print*
|
|
|
188 |
print*,'Input : ',trim(inpfile)
|
|
|
189 |
print*,'Output : ',trim(outfile)
|
|
|
190 |
print*,'Field : ',trim(field)
|
|
|
191 |
print*,'Trajectory : ',ntime,nfield,ntra
|
|
|
192 |
print*,'Grid type : ',trim(gridtype)
|
|
|
193 |
if ( dynamic_grid.eq.1 ) then
|
|
|
194 |
print*,'Grid : dynamic (see below)'
|
|
|
195 |
elseif ( gridtype.eq.'latlon' ) then
|
|
|
196 |
print*,'Grid nlon,nlat : ',nx,ny
|
|
|
197 |
print*,' lonmin,latmin : ',xmin,ymin
|
|
|
198 |
print*,' dlon,dlat : ',dx,dy
|
|
|
199 |
elseif ( gridtype.eq.'rotated' ) then
|
|
|
200 |
print*,'Grid clon,clat : ',clon,clat
|
|
|
201 |
print*,' nlonlat : ',nlonlat
|
|
|
202 |
print*,' dlonlat : ',dlonlat
|
|
|
203 |
endif
|
|
|
204 |
print*,'Filter radius : ',radius,' ',trim(runit)
|
|
|
205 |
print*,'Mode : ',trim(mode)
|
|
|
206 |
if ( ( mode.eq.'time' ).or.
|
|
|
207 |
> ( mode.eq.'space' ).or.
|
|
|
208 |
> (mode.eq.'grid' ) )
|
|
|
209 |
>then
|
|
|
210 |
print*,'Parameter : ',param
|
|
|
211 |
endif
|
|
|
212 |
if ( step.eq.0 ) then
|
|
|
213 |
print*,'Time step : all'
|
|
|
214 |
elseif (step.gt.0) then
|
|
|
215 |
print*,'Time step : ',step
|
|
|
216 |
endif
|
|
|
217 |
print*,'Selection file : ',trim(sel_file)
|
|
|
218 |
print*,'Selection format : ',trim(sel_file)
|
|
|
219 |
print*,'Flag <crefile> : ',crefile
|
|
|
220 |
print*,'Flag <crevar> : ',crevar
|
|
|
221 |
|
|
|
222 |
c Check whether mode is valid
|
|
|
223 |
if ((mode.ne.'keep' ).and.
|
|
|
224 |
> (mode.ne.'time' ).and.
|
|
|
225 |
> (mode.ne.'space' ).and.
|
|
|
226 |
> (mode.ne.'grid' ))
|
|
|
227 |
>then
|
|
|
228 |
print*,' ERROR: Invalid mode ',trim(mode)
|
|
|
229 |
stop
|
|
|
230 |
endif
|
|
|
231 |
|
|
|
232 |
c Allocate memory for old and new (reparameterised) trajectory
|
|
|
233 |
allocate(olon(ntime),stat=stat)
|
|
|
234 |
if (stat.ne.0) print*,'*** error allocating array olon ***'
|
|
|
235 |
allocate(olat(ntime),stat=stat)
|
|
|
236 |
if (stat.ne.0) print*,'*** error allocating array olat ***'
|
|
|
237 |
allocate(otim(ntime),stat=stat)
|
|
|
238 |
if (stat.ne.0) print*,'*** error allocating array otim ***'
|
|
|
239 |
allocate(nlon(1000*ntime),stat=stat)
|
|
|
240 |
if (stat.ne.0) print*,'*** error allocating array nlon ***'
|
|
|
241 |
allocate(nlat(1000*ntime),stat=stat)
|
|
|
242 |
if (stat.ne.0) print*,'*** error allocating array nlat ***'
|
|
|
243 |
allocate(ntim(1000*ntime),stat=stat)
|
|
|
244 |
if (stat.ne.0) print*,'*** error allocating array ntim ***'
|
|
|
245 |
allocate(odist(ntime),stat=stat)
|
|
|
246 |
if (stat.ne.0) print*,'*** error allocating array odist ***'
|
|
|
247 |
allocate(ndist(1000*ntime),stat=stat)
|
|
|
248 |
if (stat.ne.0) print*,'*** error allocating array ndist ***'
|
|
|
249 |
allocate(ofld(ntime),stat=stat)
|
|
|
250 |
if (stat.ne.0) print*,'*** error allocating array ofld ***'
|
|
|
251 |
allocate(nfld(1000*ntime),stat=stat)
|
|
|
252 |
if (stat.ne.0) print*,'*** error allocating array nfld ***'
|
|
|
253 |
|
|
|
254 |
c Allocate memory for complete trajectory set
|
|
|
255 |
allocate(trainp(ntra,ntime,nfield),stat=stat)
|
|
|
256 |
if (stat.ne.0) print*,'*** error allocating array trainp ***'
|
|
|
257 |
allocate(sel_flag(ntra),stat=stat)
|
|
|
258 |
if (stat.ne.0) print*,'*** error allocating array sel_flag ***'
|
|
|
259 |
|
|
|
260 |
c Allocate memory for auxiliary fields
|
|
|
261 |
allocate(traj(nfield),stat=stat)
|
|
|
262 |
if (stat.ne.0) print*,'*** error allocating array traj ***'
|
|
|
263 |
|
|
|
264 |
c Set the format of the input file
|
|
|
265 |
call mode_tra(inpmode,inpfile)
|
|
|
266 |
if (inpmode.eq.-1) inpmode=1
|
|
|
267 |
|
|
|
268 |
c Read the input trajectory file
|
|
|
269 |
call ropen_tra(fid,inpfile,ntra,ntime,nfield,
|
|
|
270 |
> reftime,varsinp,inpmode)
|
|
|
271 |
call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)
|
|
|
272 |
call close_tra(fid,inpmode)
|
|
|
273 |
|
| 15 |
michaesp |
274 |
c Bring trajectories into grid interval
|
|
|
275 |
print*,xmin,ymin,dx,dy,nx,ny
|
|
|
276 |
if ( ( abs(xmin+180.).lt.eps).and.
|
|
|
277 |
> ( abs(ymin+90. ).lt.eps).and.
|
|
|
278 |
> ( abs(dx-1. ).lt.eps).and.
|
|
|
279 |
> ( abs(dy-1. ).lt.eps).and.
|
|
|
280 |
> ( nx.eq.360 ).and.
|
|
|
281 |
> ( ny.eq.180 ) )
|
|
|
282 |
>then
|
|
|
283 |
print*,'HALLO'
|
|
|
284 |
do i=1,ntra
|
|
|
285 |
do j=1,ntime
|
|
|
286 |
if ( trainp(i,j,2).lt.-180. ) then
|
|
|
287 |
trainp(i,j,2) = trainp(i,j,2) + 360.
|
|
|
288 |
endif
|
|
|
289 |
if ( trainp(i,j,2).gt.180. ) then
|
|
|
290 |
trainp(i,j,2) = trainp(i,j,2) - 360.
|
|
|
291 |
endif
|
|
|
292 |
enddo
|
|
|
293 |
enddo
|
|
|
294 |
endif
|
|
|
295 |
|
| 3 |
michaesp |
296 |
c Check that first four columns correspond to time,lon,lat,p
|
|
|
297 |
if ( (varsinp(1).ne.'time' ).or.
|
|
|
298 |
> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
|
|
|
299 |
> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
|
|
|
300 |
> (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p' ) )
|
|
|
301 |
>then
|
|
|
302 |
print*,' ERROR: problem with input trajectories ...'
|
|
|
303 |
stop
|
|
|
304 |
endif
|
|
|
305 |
varsinp(1) = 'TIME'
|
|
|
306 |
varsinp(2) = 'lon'
|
|
|
307 |
varsinp(3) = 'lat'
|
|
|
308 |
varsinp(4) = 'p'
|
|
|
309 |
|
|
|
310 |
c Get the index of the field (if needed)
|
|
|
311 |
if ( field.ne.'nil' ) then
|
|
|
312 |
ifield = 0
|
|
|
313 |
do i=1,nfield
|
|
|
314 |
if ( varsinp(i).eq.field ) ifield = i
|
|
|
315 |
enddo
|
|
|
316 |
if ( ifield.eq.0 ) then
|
|
|
317 |
print*,' ERROR: field ',trim(field),' not found... Stop'
|
|
|
318 |
stop
|
|
|
319 |
endif
|
|
|
320 |
endif
|
|
|
321 |
|
|
|
322 |
c Write some status information of the input trajectories
|
|
|
323 |
print*
|
|
|
324 |
print*,'---- INPUT TRAJECTORIES ---------------------------------'
|
|
|
325 |
print*
|
|
|
326 |
print*,' Reference time (year) : ',reftime(1)
|
|
|
327 |
print*,' (month) : ',reftime(2)
|
|
|
328 |
print*,' (day) : ',reftime(3)
|
|
|
329 |
print*,' (hour) : ',reftime(4)
|
|
|
330 |
print*,' (min) : ',reftime(5)
|
|
|
331 |
print*,' Time range (min) : ',reftime(6)
|
|
|
332 |
do i=1,nfield
|
|
|
333 |
if ( i.ne.ifield ) then
|
|
|
334 |
print*,' Var :',i,trim(varsinp(i))
|
|
|
335 |
else
|
|
|
336 |
print*,' Var :',i,trim(varsinp(i)),
|
|
|
337 |
> ' [ gridding ]'
|
|
|
338 |
endif
|
|
|
339 |
enddo
|
|
|
340 |
print*,' List of selected times'
|
|
|
341 |
do i=1,ntime
|
|
|
342 |
if ( (step.eq.0).or.(step.eq.i) ) then
|
|
|
343 |
print*,' ',i,' -> ',trainp(1,i,1)
|
|
|
344 |
endif
|
|
|
345 |
enddo
|
|
|
346 |
print*
|
|
|
347 |
|
|
|
348 |
c Select flag: all trajectories are selected
|
|
|
349 |
if ( sel_file.eq.'nil' ) then
|
|
|
350 |
|
|
|
351 |
do i=1,ntra
|
|
|
352 |
sel_flag(i) = 1
|
|
|
353 |
enddo
|
|
|
354 |
|
|
|
355 |
c Select flag: index file
|
|
|
356 |
elseif ( sel_format.eq.'index' ) then
|
|
|
357 |
|
|
|
358 |
do i=1,ntra
|
|
|
359 |
sel_flag(i) = 0
|
|
|
360 |
enddo
|
|
|
361 |
|
|
|
362 |
open(10,file=sel_file)
|
|
|
363 |
142 read(10,*,end=141) ind
|
|
|
364 |
sel_flag(ind) = 1
|
|
|
365 |
goto 142
|
|
|
366 |
141 continue
|
|
|
367 |
close(10)
|
|
|
368 |
|
|
|
369 |
c Select flag: boolean file
|
|
|
370 |
elseif ( sel_format.eq.'boolean' ) then
|
|
|
371 |
|
|
|
372 |
open(10,file=sel_file)
|
|
|
373 |
do i=1,ntra
|
|
|
374 |
read(10,*) ind
|
|
|
375 |
if ( ind.eq.1 ) sel_flag(i) = ind
|
|
|
376 |
enddo
|
|
|
377 |
close(10)
|
|
|
378 |
|
|
|
379 |
endif
|
|
|
380 |
|
|
|
381 |
c Write status information
|
|
|
382 |
if ( sel_file.eq.'nil' ) then
|
|
|
383 |
print*,' Selected trajectories : all ',ntra
|
|
|
384 |
else
|
|
|
385 |
count = 0
|
|
|
386 |
do i=1,ntra
|
|
|
387 |
if ( sel_flag(i).eq.1 ) count = count + 1
|
|
|
388 |
enddo
|
|
|
389 |
print*,' #selected trajectories : ',count,
|
|
|
390 |
> ' [ ',real(count)/real(ntra) * 100.,' % ] '
|
|
|
391 |
endif
|
|
|
392 |
print*
|
|
|
393 |
|
|
|
394 |
c ---------------------------------------------------------------------
|
|
|
395 |
c Coordinate transformations and grid adjustment
|
|
|
396 |
c ---------------------------------------------------------------------
|
|
|
397 |
|
|
|
398 |
c Transform from lat/lon to rotated lat/lon, if requested
|
|
|
399 |
if ( gridtype.eq.'rotated') then
|
|
|
400 |
|
|
|
401 |
crot = 0.
|
|
|
402 |
|
|
|
403 |
pollon=clon-180.
|
|
|
404 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
405 |
pollat=90.-clat
|
|
|
406 |
do i=1,ntra
|
|
|
407 |
do j=1,ntime
|
|
|
408 |
|
|
|
409 |
if ( sel_flag(i).eq.1 ) then
|
|
|
410 |
|
|
|
411 |
c Get lat/lon coordinates for trajectory point
|
|
|
412 |
lon = trainp(i,j,2)
|
|
|
413 |
lat = trainp(i,j,3)
|
|
|
414 |
|
|
|
415 |
c First Rotation
|
|
|
416 |
pollon=clon-180.
|
|
|
417 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
418 |
pollat=90.-clat
|
|
|
419 |
rlon0=lmtolms(lat,lon,pollat,pollon)
|
|
|
420 |
rlat0=phtophs(lat,lon,pollat,pollon)
|
|
|
421 |
|
|
|
422 |
c Second rotation
|
|
|
423 |
pollon=-180.
|
|
|
424 |
pollat=90.+crot
|
|
|
425 |
rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)
|
|
|
426 |
rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)
|
|
|
427 |
|
|
|
428 |
c Get rotated latitude and longitude
|
|
|
429 |
100 if (rlon.lt.xmin) then
|
|
|
430 |
rlon=rlon+360.
|
|
|
431 |
goto 100
|
|
|
432 |
endif
|
|
|
433 |
102 if (rlon.gt.(xmin+real(nx-1)*dx)) then
|
|
|
434 |
rlon=rlon-360.
|
|
|
435 |
goto 102
|
|
|
436 |
endif
|
|
|
437 |
|
|
|
438 |
c Set the new trajectory coordinates
|
|
|
439 |
trainp(i,j,2) = rlon
|
|
|
440 |
trainp(i,j,3) = rlat
|
|
|
441 |
|
|
|
442 |
endif
|
|
|
443 |
|
|
|
444 |
enddo
|
|
|
445 |
enddo
|
|
|
446 |
|
|
|
447 |
endif
|
|
|
448 |
|
|
|
449 |
c Dynamic grid adjustment
|
|
|
450 |
if ( dynamic_grid.eq.1 ) then
|
|
|
451 |
|
|
|
452 |
c Get the grid parameters
|
|
|
453 |
xmin = 180.
|
|
|
454 |
ymin = 90.
|
|
|
455 |
xmax = -180.
|
|
|
456 |
ymax = -90.
|
|
|
457 |
|
|
|
458 |
do i=1,ntra
|
|
|
459 |
|
|
|
460 |
if ( sel_flag(i).eq.1 ) then
|
|
|
461 |
|
|
|
462 |
if ( step.eq.0 ) then
|
|
|
463 |
do j=1,ntime
|
|
|
464 |
if ( trainp(i,j,2).lt.xmin) xmin = trainp(i,j,2)
|
|
|
465 |
if ( trainp(i,j,2).gt.xmax) xmax = trainp(i,j,2)
|
|
|
466 |
if ( trainp(i,j,3).lt.ymin) ymin = trainp(i,j,3)
|
|
|
467 |
if ( trainp(i,j,3).gt.ymax) ymax = trainp(i,j,3)
|
|
|
468 |
enddo
|
|
|
469 |
else
|
|
|
470 |
if ( trainp(i,step,2).lt.xmin) xmin = trainp(i,step,2)
|
|
|
471 |
if ( trainp(i,step,2).gt.xmax) xmax = trainp(i,step,2)
|
|
|
472 |
if ( trainp(i,step,3).lt.ymin) ymin = trainp(i,step,3)
|
|
|
473 |
if ( trainp(i,step,3).gt.ymax) ymax = trainp(i,step,3)
|
|
|
474 |
endif
|
|
|
475 |
|
|
|
476 |
endif
|
|
|
477 |
|
|
|
478 |
enddo
|
|
|
479 |
|
|
|
480 |
c Get first guess for "optimal" grid
|
|
|
481 |
nx = 400
|
|
|
482 |
ny = 400
|
|
|
483 |
dx = (xmax - xmin)/real(nx-1)
|
|
|
484 |
dy = (ymax - ymin)/real(ny-1)
|
|
|
485 |
|
|
|
486 |
c Make the grid spacing equal in zonal and meridional direction
|
|
|
487 |
if ( dx.gt.dy ) then
|
|
|
488 |
|
|
|
489 |
dy = dx
|
|
|
490 |
ny = (ymax - ymin)/dy + 1
|
|
|
491 |
if (ny.lt.nx/2) ny = nx / 2
|
|
|
492 |
if ( real(ny)*dy .ge. 180. ) ny = 180./dy + 1
|
|
|
493 |
ycen = 0.5* (ymin+ymax)
|
|
|
494 |
ymin = ycen - 0.5 * real(ny/2) * dy
|
|
|
495 |
if (ymin.le.-90.) ymin = -90.
|
|
|
496 |
|
|
|
497 |
else
|
|
|
498 |
|
|
|
499 |
dx = dy
|
|
|
500 |
nx = (xmax - xmin)/dx + 1
|
|
|
501 |
if (nx.lt.ny/2) nx = ny / 2
|
|
|
502 |
if ( real(nx)*dx .ge. 360. ) nx = 360./dx + 1
|
|
|
503 |
xcen = 0.5* (xmin+xmax)
|
|
|
504 |
xmin = xcen - 0.5 * real(nx/2) * dx
|
|
|
505 |
if (xmin.le.-180.) xmin = -180.
|
|
|
506 |
|
|
|
507 |
endif
|
|
|
508 |
|
|
|
509 |
c Write information
|
|
|
510 |
print*
|
|
|
511 |
print*,'---- DYNAMIC GRID ADJUSTMENT',
|
|
|
512 |
> ' ----------------------------'
|
|
|
513 |
print*
|
|
|
514 |
print*,'Grid nlon,nlat : ',nx,ny
|
|
|
515 |
print*,' lonmin,latmin : ',xmin,ymin
|
|
|
516 |
print*,' dlon,dlat : ',dx,dy
|
|
|
517 |
print*
|
|
|
518 |
|
|
|
519 |
c Write grid information for rotated grid (if not already done
|
|
|
520 |
elseif ( gridtype.eq.'rotated') then
|
|
|
521 |
|
|
|
522 |
print*
|
|
|
523 |
print*,'---- GRID PARAMETERS -------',
|
|
|
524 |
> ' ----------------------------'
|
|
|
525 |
print*
|
|
|
526 |
print*,'Grid nlon,nlat : ',nx,ny
|
|
|
527 |
print*,' lonmin,latmin : ',xmin,ymin
|
|
|
528 |
print*,' dlon,dlat : ',dx,dy
|
|
|
529 |
print*
|
|
|
530 |
|
|
|
531 |
|
|
|
532 |
endif
|
|
|
533 |
|
|
|
534 |
c Set the grid boundaries
|
|
|
535 |
xmax=xmin+real(nx-1)*dx
|
|
|
536 |
ymax=ymin+real(ny-1)*dy
|
|
|
537 |
|
|
|
538 |
c Allocate memory for output array and auxiliary gridding array
|
|
|
539 |
allocate(cnt(nx,ny),stat=stat)
|
|
|
540 |
if (stat.ne.0) print*,'*** error allocating array cnt ***'
|
|
|
541 |
allocate(res(nx,ny),stat=stat)
|
|
|
542 |
if (stat.ne.0) print*,'*** error allocating array res ***'
|
|
|
543 |
allocate(fld(nx,ny),stat=stat)
|
|
|
544 |
if (stat.ne.0) print*,'*** error allocating array fld ***'
|
|
|
545 |
allocate(area(nx,ny),stat=stat)
|
|
|
546 |
if (stat.ne.0) print*,'*** error allocating array area ***'
|
|
|
547 |
|
|
|
548 |
allocate(connect0(nx,ny),stat=stat)
|
|
|
549 |
if (stat.ne.0) print*,'*** error allocating array connect0 ***'
|
|
|
550 |
allocate(connect1(nx,ny),stat=stat)
|
|
|
551 |
if (stat.ne.0) print*,'*** error allocating array connect1 ***'
|
|
|
552 |
allocate(connect2(nx,ny),stat=stat)
|
|
|
553 |
if (stat.ne.0) print*,'*** error allocating array connect2 ***'
|
|
|
554 |
|
|
|
555 |
|
|
|
556 |
c Init the output array
|
|
|
557 |
do i=1,nx
|
|
|
558 |
do j=1,ny
|
|
|
559 |
connect0(i,j) = 0
|
|
|
560 |
connect1(i,j) = 0
|
|
|
561 |
connect2(i,j) = 0
|
|
|
562 |
cnt(i,j) = 0.
|
|
|
563 |
res(i,j) = 0.
|
|
|
564 |
fld(i,j) = 0.
|
|
|
565 |
enddo
|
|
|
566 |
enddo
|
|
|
567 |
|
|
|
568 |
c ---------------------------------------------------------------------
|
|
|
569 |
c Gridding
|
|
|
570 |
c ---------------------------------------------------------------------
|
|
|
571 |
|
|
|
572 |
c Write some status information
|
|
|
573 |
print*,'---- GRIDDING -------------------------------------------'
|
|
|
574 |
print*
|
|
|
575 |
|
|
|
576 |
c Loop over all entries of sampling table
|
|
|
577 |
connectval0 = 0
|
|
|
578 |
connectval1 = 0
|
|
|
579 |
connectval2 = 0
|
|
|
580 |
count = 0
|
|
|
581 |
|
|
|
582 |
do i=1,ntra
|
|
|
583 |
|
|
|
584 |
if (mod(i,100).eq.0) print*,i,' of ',ntra
|
|
|
585 |
|
|
|
586 |
c Skip all trajectories which are not selected
|
|
|
587 |
if ( sel_flag(i).eq.0 ) goto 300
|
|
|
588 |
|
|
|
589 |
c ------- Read a complete trajectory ---------------------------
|
|
|
590 |
do j=1,ntime
|
|
|
591 |
otim(j) = trainp(i,j,1)
|
|
|
592 |
olon(j) = trainp(i,j,2)
|
|
|
593 |
olat(j) = trainp(i,j,3)
|
|
|
594 |
if ( field.ne.'nil' ) then
|
|
|
595 |
ofld(j) =trainp(i,j,ifield)
|
|
|
596 |
endif
|
|
|
597 |
enddo
|
|
|
598 |
|
|
|
599 |
c -------- Convert hh.m time into fractional time --------------
|
|
|
600 |
do j=1,ntime
|
|
|
601 |
hhmm = otim(j)
|
|
|
602 |
call hhmm2frac (hhmm,frac)
|
|
|
603 |
otim(j) = frac
|
|
|
604 |
enddo
|
|
|
605 |
|
|
|
606 |
c -------- Interpolation ---------------------------------------
|
|
|
607 |
|
|
|
608 |
c Keep the trajectory points as they are
|
|
|
609 |
if ( ( mode.eq.'keep').and.(step.eq.0) ) then
|
|
|
610 |
npts=opts
|
|
|
611 |
do j=1,opts
|
|
|
612 |
ntim(j)=otim(j)
|
|
|
613 |
nlon(j)=olon(j)
|
|
|
614 |
nlat(j)=olat(j)
|
|
|
615 |
if ( field.ne.'nil' ) then
|
|
|
616 |
nfld(j)=ofld(j)
|
|
|
617 |
endif
|
|
|
618 |
enddo
|
|
|
619 |
|
|
|
620 |
c Select a single time step
|
|
|
621 |
elseif ( ( mode.eq.'keep').and.(step.gt.0) ) then
|
|
|
622 |
npts = 1
|
|
|
623 |
ntim(1) = otim(step)
|
|
|
624 |
nlon(1) = olon(step)
|
|
|
625 |
nlat(1) = olat(step)
|
|
|
626 |
if ( field.ne.'nil' ) then
|
|
|
627 |
nfld(1) = ofld(step)
|
|
|
628 |
endif
|
|
|
629 |
|
|
|
630 |
c Perform a reparameterisation in time
|
|
|
631 |
else if ( (mode.eq.'time').and.(step.eq.0) ) then
|
|
|
632 |
|
|
|
633 |
c Get the new number of trajectory points
|
|
|
634 |
npts=nint(abs(otim(opts)-otim(1))/param)+1
|
|
|
635 |
|
|
|
636 |
c Handle date line problem
|
|
|
637 |
do j=2,opts
|
|
|
638 |
if ( (olon(j-1)-olon(j)).gt.180. ) then
|
|
|
639 |
olon(j) = olon(j) + 360.
|
|
|
640 |
else if ( (olon(j-1)-olon(j)).lt.-180. ) then
|
|
|
641 |
olon(j) = olon(j) - 360.
|
|
|
642 |
endif
|
|
|
643 |
enddo
|
|
|
644 |
|
|
|
645 |
c Cubic spline fitting
|
|
|
646 |
call curvefit(otim,olon,opts,ntim,nlon,npts)
|
|
|
647 |
call curvefit(otim,olat,opts,ntim,nlat,npts)
|
|
|
648 |
if ( field.ne.'nil' ) then
|
|
|
649 |
call curvefit(otim,ofld,opts,ntim,nfld,npts)
|
|
|
650 |
endif
|
|
|
651 |
|
|
|
652 |
c Reverse date line handling
|
|
|
653 |
do j=1,npts
|
|
|
654 |
if ( nlon(j).gt.xmax ) then
|
|
|
655 |
nlon(j) = nlon(j) -360.
|
|
|
656 |
else if ( nlon(j).lt.xmin ) then
|
|
|
657 |
nlon(j) = nlon(j) +360.
|
|
|
658 |
endif
|
|
|
659 |
enddo
|
|
|
660 |
|
|
|
661 |
c Perform a reparameterisation with equally spaced gridpoint
|
|
|
662 |
elseif ( (mode.eq.'space').and.(step.eq.0) ) then
|
|
|
663 |
|
|
|
664 |
c Calculate the distance and spacing
|
|
|
665 |
odist(1) = 0.
|
|
|
666 |
unit = 'km'
|
|
|
667 |
do j=2,ntime
|
|
|
668 |
odist(j)=odist(j-1) +
|
|
|
669 |
> sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
|
|
|
670 |
enddo
|
|
|
671 |
|
|
|
672 |
c Determine the new number of trajectory points
|
|
|
673 |
npts=nint(odist(ntime)/param)+1
|
|
|
674 |
if (npts.eq.0) then
|
|
|
675 |
npts=1.
|
|
|
676 |
endif
|
|
|
677 |
|
|
|
678 |
c Handle date line problem
|
|
|
679 |
do j=2,opts
|
|
|
680 |
if ( (olon(j-1)-olon(j)).gt.180. ) then
|
|
|
681 |
olon(j) = olon(j) + 360.
|
|
|
682 |
else if ( (olon(j-1)-olon(j)).lt.-180. ) then
|
|
|
683 |
olon(j) = olon(j) - 360.
|
|
|
684 |
endif
|
|
|
685 |
enddo
|
|
|
686 |
|
|
|
687 |
c Cubic spline fitting
|
|
|
688 |
call curvefit(odist,olon,opts,ndist,nlon,npts)
|
|
|
689 |
call curvefit(odist,olat,opts,ndist,nlat,npts)
|
|
|
690 |
call curvefit(odist,otim,opts,ndist,ntim,npts)
|
|
|
691 |
if ( field.ne.'nil' ) then
|
|
|
692 |
call curvefit(odist,ofld,opts,ndist,nfld,npts)
|
|
|
693 |
endif
|
|
|
694 |
|
|
|
695 |
c Reverse date line handling
|
|
|
696 |
do j=1,npts
|
|
|
697 |
if ( nlon(j).gt.xmax ) then
|
|
|
698 |
nlon(j) = nlon(j) -360.
|
|
|
699 |
else if ( nlon(j).lt.xmin ) then
|
|
|
700 |
nlon(j) = nlon(j) +360.
|
|
|
701 |
endif
|
|
|
702 |
enddo
|
|
|
703 |
|
|
|
704 |
c Perform a reparameterisation with equally spaced gridpoint
|
|
|
705 |
elseif ( (mode.eq.'grid').and.(step.eq.0) ) then
|
|
|
706 |
|
|
|
707 |
c Calculate the distance and spacing
|
|
|
708 |
odist(1) = 0.
|
|
|
709 |
unit = 'deg'
|
|
|
710 |
do j=2,ntime
|
|
|
711 |
odist(j)=odist(j-1) +
|
|
|
712 |
> sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
|
|
|
713 |
enddo
|
|
|
714 |
|
|
|
715 |
c Determine the new number of trajectory points
|
|
|
716 |
npts=nint(odist(ntime)/param)+1
|
|
|
717 |
if (npts.eq.0) then
|
|
|
718 |
npts=1.
|
|
|
719 |
endif
|
|
|
720 |
|
|
|
721 |
c Handle date line problem
|
|
|
722 |
do j=2,opts
|
|
|
723 |
if ( (olon(j-1)-olon(j)).gt.180. ) then
|
|
|
724 |
olon(j) = olon(j) + 360.
|
|
|
725 |
else if ( (olon(j-1)-olon(j)).lt.-180. ) then
|
|
|
726 |
olon(j) = olon(j) - 360.
|
|
|
727 |
endif
|
|
|
728 |
enddo
|
|
|
729 |
|
|
|
730 |
c Cubic spline fitting
|
|
|
731 |
call curvefit(odist,olon,opts,ndist,nlon,npts)
|
|
|
732 |
call curvefit(odist,olat,opts,ndist,nlat,npts)
|
|
|
733 |
call curvefit(odist,otim,opts,ndist,ntim,npts)
|
|
|
734 |
if ( field.ne.'nil' ) then
|
|
|
735 |
call curvefit(odist,ofld,opts,ndist,nfld,npts)
|
|
|
736 |
endif
|
|
|
737 |
|
|
|
738 |
c Reverse date line handling
|
|
|
739 |
do j=1,npts
|
|
|
740 |
if ( nlon(j).gt.xmax ) then
|
|
|
741 |
nlon(j) = nlon(j) -360.
|
|
|
742 |
else if ( nlon(j).lt.xmin ) then
|
|
|
743 |
nlon(j) = nlon(j) +360.
|
|
|
744 |
endif
|
|
|
745 |
enddo
|
|
|
746 |
|
|
|
747 |
endif
|
|
|
748 |
|
|
|
749 |
c -------- Do the gridding -------------------------------------
|
|
|
750 |
|
|
|
751 |
c Gridding of trajectory
|
|
|
752 |
do j=1,npts
|
|
|
753 |
|
|
|
754 |
c Check whether point is in data domain
|
|
|
755 |
if ( (nlon(j).gt.xmin).and.(nlon(j).lt.xmax).and.
|
| 15 |
michaesp |
756 |
> (nlat(j).gt.ymin).and.(nlat(j).lt.ymax))
|
| 3 |
michaesp |
757 |
> then
|
|
|
758 |
|
|
|
759 |
c Increase counter for gridded points
|
|
|
760 |
count = count + 1
|
|
|
761 |
|
|
|
762 |
c ----------------- Gridding: simple count -----------------
|
|
|
763 |
connectval0 = connectval0+1
|
|
|
764 |
|
|
|
765 |
addvalue = 1.
|
|
|
766 |
|
|
|
767 |
call gridding1
|
|
|
768 |
> (nlat(j),nlon(j),addvalue,
|
|
|
769 |
> radius,runit,connect0,connectval0,
|
|
|
770 |
> cnt,nx,ny,xmin,ymin,dx,dy)
|
|
|
771 |
|
|
|
772 |
c ----------------- Gridding: residence time ---------------
|
|
|
773 |
connectval1 = connectval1+1
|
|
|
774 |
|
|
|
775 |
if ( ntime.eq.1 ) then
|
|
|
776 |
addvalue = 0.
|
|
|
777 |
elseif ( j.eq.1 ) then
|
|
|
778 |
addvalue=abs(ntim(2)-ntim(1))
|
|
|
779 |
else
|
|
|
780 |
addvalue=abs(ntim(j)-ntim(j-1))
|
|
|
781 |
endif
|
|
|
782 |
|
|
|
783 |
call gridding1
|
|
|
784 |
> (nlat(j),nlon(j),addvalue,
|
|
|
785 |
> radius,runit,connect1,connectval1,
|
|
|
786 |
> res,nx,ny,xmin,ymin,dx,dy)
|
|
|
787 |
|
|
|
788 |
|
|
|
789 |
c --------------- Gridding: field -------------------------
|
|
|
790 |
if ( field.ne.'nil' ) then
|
|
|
791 |
|
|
|
792 |
connectval2 = connectval2+1
|
|
|
793 |
|
|
|
794 |
addvalue = nfld(j)
|
|
|
795 |
|
|
|
796 |
call gridding1
|
|
|
797 |
> (nlat(j),nlon(j),addvalue,
|
|
|
798 |
> radius,runit,connect2,connectval2,
|
|
|
799 |
> fld,nx,ny,xmin,ymin,dx,dy)
|
|
|
800 |
|
|
|
801 |
endif
|
|
|
802 |
|
|
|
803 |
endif
|
|
|
804 |
|
|
|
805 |
enddo
|
|
|
806 |
|
|
|
807 |
c Exit point for loop over all trajectories
|
|
|
808 |
300 continue
|
|
|
809 |
|
|
|
810 |
enddo
|
|
|
811 |
|
|
|
812 |
c Write status information
|
|
|
813 |
print*
|
|
|
814 |
print*,' # gridded points : ',count
|
|
|
815 |
|
|
|
816 |
c ---------------------------------------------------------------------
|
|
|
817 |
c Unit conversions and output to netCDF file
|
|
|
818 |
c ---------------------------------------------------------------------
|
|
|
819 |
|
|
|
820 |
c Write some status information
|
|
|
821 |
print*
|
|
|
822 |
print*,'---- WRITE OUTPUT ---------------------------------------'
|
|
|
823 |
print*
|
|
|
824 |
|
|
|
825 |
c Area (in km^2)
|
|
|
826 |
do i=1,nx
|
|
|
827 |
do j=1,ny
|
|
|
828 |
slat=ymin+real(j-1)*dy
|
|
|
829 |
if (abs(abs(slat)-90.).gt.eps) then
|
|
|
830 |
area(i,j) = dy*dx*cos(pi180*slat)*deltay**2
|
|
|
831 |
else
|
|
|
832 |
area(i,j) = 0.
|
|
|
833 |
endif
|
|
|
834 |
enddo
|
|
|
835 |
enddo
|
|
|
836 |
|
|
|
837 |
c Normalise gridded field
|
|
|
838 |
if ( field.ne.'nil' ) then
|
|
|
839 |
do i=1,nx
|
|
|
840 |
do j=1,ny
|
|
|
841 |
if ( cnt(i,j).gt.0. ) then
|
|
|
842 |
fld(i,j) = fld(i,j) / cnt(i,j)
|
|
|
843 |
endif
|
|
|
844 |
enddo
|
|
|
845 |
enddo
|
|
|
846 |
endif
|
|
|
847 |
|
|
|
848 |
c Set the time for the output netCDF files - if a composite is
|
|
|
849 |
c calculatd, then the time is set to
|
|
|
850 |
if ( step.eq.0 ) then
|
|
|
851 |
time = -999.
|
|
|
852 |
print*,' ... COMPOSITE OVER ALL TRAJECTORY TIMES (-999)'
|
|
|
853 |
print*
|
|
|
854 |
else
|
|
|
855 |
time = trainp(1,step,1)
|
|
|
856 |
endif
|
|
|
857 |
|
|
|
858 |
c Write output to CF netCDF
|
|
|
859 |
cdfname = outfile
|
|
|
860 |
|
|
|
861 |
varname = 'COUNT'
|
|
|
862 |
longname = 'trajectory counts'
|
|
|
863 |
varunit = 'counts per grid point'
|
|
|
864 |
call writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
|
|
|
865 |
> clon,clat,nlonlat,dlonlat,cnt,time,dx,dy,xmin,ymin,nx,
|
|
|
866 |
> ny,crefile,crefile,1)
|
|
|
867 |
write(*,'(a8,a10,a5,a10,a10,f7.2,a2)')
|
|
|
868 |
> ' ... ',trim(varname),' -> ',trim(cdfname),
|
|
|
869 |
> ' [ time = ',time,' ]'
|
|
|
870 |
|
|
|
871 |
varname = 'RESIDENCE'
|
|
|
872 |
longname = 'residence time'
|
|
|
873 |
varunit = 'hours per grid point'
|
| 13 |
michaesp |
874 |
|
| 3 |
michaesp |
875 |
call writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
|
|
|
876 |
> clon,clat,nlonlat,dlonlat,res,time,dx,dy,xmin,ymin,nx,
|
|
|
877 |
> ny,0,crefile,1)
|
|
|
878 |
write(*,'(a8,a10,a5,a10,a10,f7.2,a2)')
|
|
|
879 |
> ' ... ',trim(varname),' -> ',trim(cdfname),
|
|
|
880 |
> ' [ time = ',time,' ]'
|
|
|
881 |
|
|
|
882 |
varname = 'AREA'
|
|
|
883 |
longname = 'area corresponding to grid points'
|
|
|
884 |
varunit = 'square kilometers'
|
|
|
885 |
call writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
|
|
|
886 |
> clon,clat,nlonlat,dlonlat,area,time,dx,dy,xmin,ymin,nx,
|
|
|
887 |
> ny,0,crefile,1)
|
|
|
888 |
write(*,'(a8,a10,a5,a10,a10,f7.2,a2)')
|
|
|
889 |
> ' ... ',trim(varname),' -> ',trim(cdfname),
|
|
|
890 |
> ' [ time = ',time,' ]'
|
|
|
891 |
|
|
|
892 |
if ( field.ne.'nil' ) then
|
|
|
893 |
varname = field
|
|
|
894 |
longname = field
|
|
|
895 |
varunit = 'as on trajectory file'
|
|
|
896 |
call writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
|
|
|
897 |
> clon,clat,nlonlat,dlonlat,fld,time,dx,dy,xmin,ymin,nx,
|
|
|
898 |
> ny,0,crevar,1)
|
|
|
899 |
|
|
|
900 |
write(*,'(a8,a10,a5,a10,a10,f7.2,a2)')
|
|
|
901 |
> ' ... ',trim(varname),' -> ',trim(cdfname),
|
|
|
902 |
> ' [ time = ',time,' ]'
|
|
|
903 |
endif
|
|
|
904 |
|
|
|
905 |
c Write status information
|
|
|
906 |
print*
|
|
|
907 |
print*,' *** END OF PROGRAM DENSITY **'
|
|
|
908 |
print*,'========================================================='
|
|
|
909 |
|
|
|
910 |
end
|
|
|
911 |
|
|
|
912 |
c ********************************************************************
|
|
|
913 |
c * GRIDDING SUBROUTINES *
|
|
|
914 |
c ********************************************************************
|
|
|
915 |
|
|
|
916 |
c ---------------------------------------------------------------------
|
|
|
917 |
c Gridding of one single data point (smoothing in km, deg, gridp)
|
|
|
918 |
c ---------------------------------------------------------------------
|
|
|
919 |
|
|
|
920 |
subroutine gridding1 (lat,lon,addval,radius,unit,
|
|
|
921 |
> connect,connectval,
|
|
|
922 |
> out,nx,ny,xmin,ymin,dx,dy)
|
|
|
923 |
|
|
|
924 |
implicit none
|
|
|
925 |
|
|
|
926 |
c Declaration of subroutine parameters
|
|
|
927 |
real lat,lon
|
|
|
928 |
integer nx,ny
|
|
|
929 |
real xmin,ymin,dx,dy
|
|
|
930 |
real out(nx,ny)
|
|
|
931 |
real radius
|
|
|
932 |
character*80 unit
|
|
|
933 |
integer connectval
|
|
|
934 |
integer connect(nx,ny)
|
|
|
935 |
real addval
|
|
|
936 |
|
|
|
937 |
c Auxiliary variables
|
|
|
938 |
integer i,j,k
|
|
|
939 |
integer mu,md,nr,nl,n,m
|
|
|
940 |
integer stackx(nx*ny),stacky(nx*ny)
|
|
|
941 |
integer tab_x(nx*ny),tab_y(nx*ny)
|
|
|
942 |
real tab_r(nx*ny)
|
|
|
943 |
integer sp
|
|
|
944 |
real lat2,lon2
|
|
|
945 |
real dist,sum
|
|
|
946 |
real xmax
|
|
|
947 |
integer periodic
|
|
|
948 |
integer test
|
|
|
949 |
|
|
|
950 |
c Numerical epsilon
|
|
|
951 |
real eps
|
|
|
952 |
parameter (eps=0.01)
|
|
|
953 |
|
|
|
954 |
c Externals
|
|
|
955 |
real sdis,weight
|
|
|
956 |
external sdis,weight
|
|
|
957 |
|
|
|
958 |
c Check whether lat/lon point is valid
|
|
|
959 |
xmax=xmin+real(nx-1)*dx
|
|
|
960 |
if (lon.lt.xmin-eps) lon=lon+360.
|
|
|
961 |
if (lon.gt.xmax+eps) lon=lon-360.
|
|
|
962 |
if (abs(lat-90).lt.eps) lat=90.
|
|
|
963 |
if (abs(lat+90).lt.eps) lat=-90.
|
|
|
964 |
if ((abs(lat).gt.(90.+eps)).or.
|
|
|
965 |
> (lon.lt.xmin-eps).or.(lon.gt.xmax+eps)) then
|
|
|
966 |
print*,'Invalid lat/lon point ',lat,lon
|
|
|
967 |
return
|
|
|
968 |
endif
|
|
|
969 |
|
|
|
970 |
c Set flag for periodic domain
|
|
|
971 |
if (abs(xmax-xmin-360.).lt.eps) then
|
|
|
972 |
periodic=1
|
|
|
973 |
else if (abs(xmax-xmin-360+dx).lt.eps) then
|
|
|
974 |
periodic=2
|
|
|
975 |
else
|
|
|
976 |
periodic=0
|
|
|
977 |
endif
|
|
|
978 |
|
|
|
979 |
c Get indices of one coarse grid point within search radius
|
|
|
980 |
i=nint((lon-xmin)/dx)+1
|
|
|
981 |
if ((i.eq.nx).and.(periodic.eq.1)) i=1
|
|
|
982 |
j=nint((lat-ymin)/dy)+1
|
|
|
983 |
lat2=ymin+real(j-1)*dy
|
|
|
984 |
lon2=xmin+real(i-1)*dx
|
|
|
985 |
dist=sdis(lon,lat,lon2,lat2,unit)
|
|
|
986 |
if (dist.gt.radius) then
|
|
|
987 |
print*,'1: Search radius is too small...'
|
|
|
988 |
stop
|
|
|
989 |
endif
|
|
|
990 |
|
|
|
991 |
c Get connected points
|
|
|
992 |
k=0
|
|
|
993 |
stackx(1)=i
|
|
|
994 |
stacky(1)=j
|
|
|
995 |
sp=1
|
|
|
996 |
do while (sp.ne.0)
|
|
|
997 |
|
|
|
998 |
c Get an element from stack
|
|
|
999 |
n=stackx(sp)
|
|
|
1000 |
m=stacky(sp)
|
|
|
1001 |
sp=sp-1
|
|
|
1002 |
|
|
|
1003 |
c Get distance from reference point
|
|
|
1004 |
lat2=ymin+real(m-1)*dy
|
|
|
1005 |
lon2=xmin+real(n-1)*dx
|
|
|
1006 |
dist=sdis(lon,lat,lon2,lat2,unit)
|
|
|
1007 |
|
|
|
1008 |
c Check whether distance is smaller than search radius: connected
|
|
|
1009 |
if (dist.lt.radius) then
|
|
|
1010 |
|
|
|
1011 |
c Make entry in filter mask
|
|
|
1012 |
k=k+1
|
|
|
1013 |
tab_x(k)=n
|
|
|
1014 |
tab_y(k)=m
|
|
|
1015 |
tab_r(k)=weight(dist,radius)
|
|
|
1016 |
|
|
|
1017 |
c Mark this point as visited
|
|
|
1018 |
connect(n,m)=connectval
|
|
|
1019 |
|
|
|
1020 |
c Get coordinates of neighbouring points
|
|
|
1021 |
nr=n+1
|
|
|
1022 |
if ((nr.gt.nx) .and.(periodic.eq.0)) nr=nx
|
|
|
1023 |
if ((nr.gt.nx-1).and.(periodic.eq.1)) nr=1
|
|
|
1024 |
if ((nr.gt.nx) .and.(periodic.eq.2)) nr=1
|
|
|
1025 |
nl=n-1
|
|
|
1026 |
if ((nl.lt.1).and.(periodic.eq.0)) nl=1
|
|
|
1027 |
if ((nl.lt.1).and.(periodic.eq.1)) nl=nx-1
|
|
|
1028 |
if ((nl.lt.1).and.(periodic.eq.2)) nl=nx
|
|
|
1029 |
mu=m+1
|
|
|
1030 |
if (mu.gt.ny) mu=ny
|
|
|
1031 |
md=m-1
|
|
|
1032 |
if (md.lt.1) md=1
|
|
|
1033 |
|
|
|
1034 |
c Update stack
|
|
|
1035 |
if (connect(nr,m).ne.connectval) then
|
|
|
1036 |
connect(nr,m)=connectval
|
|
|
1037 |
sp=sp+1
|
|
|
1038 |
stackx(sp)=nr
|
|
|
1039 |
stacky(sp)=m
|
|
|
1040 |
endif
|
|
|
1041 |
if (connect(nl,m).ne.connectval) then
|
|
|
1042 |
connect(nl,m)=connectval
|
|
|
1043 |
sp=sp+1
|
|
|
1044 |
stackx(sp)=nl
|
|
|
1045 |
stacky(sp)=m
|
|
|
1046 |
endif
|
|
|
1047 |
if (connect(n,mu).ne.connectval) then
|
|
|
1048 |
connect(n,mu)=connectval
|
|
|
1049 |
sp=sp+1
|
|
|
1050 |
stackx(sp)=n
|
|
|
1051 |
stacky(sp)=mu
|
|
|
1052 |
endif
|
|
|
1053 |
if (connect(n,md).ne.connectval) then
|
|
|
1054 |
connect(n,md)=connectval
|
|
|
1055 |
sp=sp+1
|
|
|
1056 |
stackx(sp)=n
|
|
|
1057 |
stacky(sp)=md
|
|
|
1058 |
endif
|
|
|
1059 |
endif
|
|
|
1060 |
|
|
|
1061 |
end do
|
|
|
1062 |
|
|
|
1063 |
if (k.ge.1) then
|
|
|
1064 |
sum=0.
|
|
|
1065 |
do i=1,k
|
|
|
1066 |
sum=sum+tab_r(i)
|
|
|
1067 |
enddo
|
|
|
1068 |
do i=1,k
|
|
|
1069 |
out(tab_x(i),tab_y(i))=out(tab_x(i),tab_y(i))+
|
|
|
1070 |
> addval*tab_r(i)/sum
|
|
|
1071 |
|
|
|
1072 |
if ((tab_x(i).eq.1).and.(periodic.eq.1)) then
|
|
|
1073 |
out(nx,tab_y(i))=out(nx,tab_y(i))+
|
|
|
1074 |
> addval*tab_r(i)/sum
|
|
|
1075 |
endif
|
|
|
1076 |
enddo
|
|
|
1077 |
else
|
|
|
1078 |
print*,'2: Search radius is too small...'
|
|
|
1079 |
stop
|
|
|
1080 |
endif
|
|
|
1081 |
|
|
|
1082 |
end
|
|
|
1083 |
|
|
|
1084 |
|
|
|
1085 |
c ----------------------------------------------------------------------
|
|
|
1086 |
c Get spherical distance between lat/lon points
|
|
|
1087 |
c ----------------------------------------------------------------------
|
|
|
1088 |
|
|
|
1089 |
real function sdis(xp,yp,xq,yq,unit)
|
|
|
1090 |
|
|
|
1091 |
c Calculates spherical distance (in km) between two points given
|
|
|
1092 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
1093 |
|
|
|
1094 |
real re
|
|
|
1095 |
parameter (re=6370.)
|
|
|
1096 |
real xp,yp,xq,yq,arg
|
|
|
1097 |
character*80 unit
|
|
|
1098 |
real dlon
|
|
|
1099 |
|
|
|
1100 |
if ( unit.eq.'km' ) then
|
|
|
1101 |
|
|
|
1102 |
arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
|
|
|
1103 |
if (arg.lt.-1.) arg=-1.
|
|
|
1104 |
if (arg.gt.1.) arg=1.
|
|
|
1105 |
sdis=re*acos(arg)
|
|
|
1106 |
|
|
|
1107 |
elseif ( unit.eq.'deg' ) then
|
|
|
1108 |
|
|
|
1109 |
dlon = xp-xq
|
|
|
1110 |
if ( dlon.gt. 180. ) dlon = dlon - 360.
|
|
|
1111 |
if ( dlon.lt.-180. ) dlon = dlon + 360.
|
|
|
1112 |
sdis = sqrt( dlon**2 + (yp-yq)**2 )
|
|
|
1113 |
|
|
|
1114 |
endif
|
|
|
1115 |
|
|
|
1116 |
|
|
|
1117 |
c Quick and dirty trick to avoid zero distances
|
|
|
1118 |
if (sdis.eq.0.) sdis=0.1
|
|
|
1119 |
|
|
|
1120 |
end
|
|
|
1121 |
|
|
|
1122 |
c ----------------------------------------------------------------------
|
|
|
1123 |
c Weight function for the filter mask
|
|
|
1124 |
c ----------------------------------------------------------------------
|
|
|
1125 |
|
|
|
1126 |
real function weight (r,radius)
|
|
|
1127 |
|
|
|
1128 |
c Attribute to each distanc r its corresponding weight in the filter mask
|
|
|
1129 |
|
|
|
1130 |
implicit none
|
|
|
1131 |
|
|
|
1132 |
c Declaration of subroutine parameters
|
|
|
1133 |
real r
|
|
|
1134 |
real radius
|
|
|
1135 |
|
|
|
1136 |
c Simple 0/1 mask
|
|
|
1137 |
if (r.lt.radius) then
|
|
|
1138 |
weight=exp(-r/radius)
|
|
|
1139 |
else
|
|
|
1140 |
weight=0.
|
|
|
1141 |
endif
|
|
|
1142 |
|
|
|
1143 |
end
|
|
|
1144 |
|
|
|
1145 |
|
|
|
1146 |
c ********************************************************************
|
|
|
1147 |
c * REPARAMETERIZATION SUBROUTINES *
|
|
|
1148 |
c ********************************************************************
|
|
|
1149 |
|
|
|
1150 |
c -------------------------------------------------------------
|
|
|
1151 |
c Interpolation of the trajectory with a natural cubic spline
|
|
|
1152 |
c -------------------------------------------------------------
|
|
|
1153 |
|
|
|
1154 |
SUBROUTINE curvefit (time,lon,n,
|
|
|
1155 |
> sptime,splon,spn)
|
|
|
1156 |
|
|
|
1157 |
c Given the curve <time,lon> with <n> data points, fit a
|
|
|
1158 |
c cubic spline to this curve. The new curve is returned in
|
|
|
1159 |
c <sptime,splon,spn> with <spn> data points. The parameter
|
|
|
1160 |
c <spn> specifies on entry the number of spline interpolated points
|
|
|
1161 |
c along the curve.
|
|
|
1162 |
|
|
|
1163 |
implicit none
|
|
|
1164 |
|
|
|
1165 |
c Declaration of subroutine parameters
|
|
|
1166 |
integer n
|
|
|
1167 |
real time(n),lon(n)
|
|
|
1168 |
integer spn
|
|
|
1169 |
real sptime(spn),splon(spn)
|
|
|
1170 |
|
|
|
1171 |
c Auxiliary variables
|
|
|
1172 |
real y2ax(n)
|
|
|
1173 |
real dt
|
|
|
1174 |
real s
|
|
|
1175 |
integer i
|
|
|
1176 |
real order
|
|
|
1177 |
|
|
|
1178 |
c Determine whether the input array is ascending or descending
|
|
|
1179 |
if (time(1).gt.time(n)) then
|
|
|
1180 |
order=-1.
|
|
|
1181 |
else
|
|
|
1182 |
order= 1.
|
|
|
1183 |
endif
|
|
|
1184 |
|
|
|
1185 |
c Bring the time array into ascending order
|
|
|
1186 |
do i=1,n
|
|
|
1187 |
time(i)=order*time(i)
|
|
|
1188 |
enddo
|
|
|
1189 |
|
|
|
1190 |
c Prepare the (natural) cubic spline interpolation
|
|
|
1191 |
call spline (time,lon,n,1.e30,1.e30,y2ax)
|
|
|
1192 |
dt=(time(n)-time(1))/real(spn-1)
|
|
|
1193 |
do i=1,spn
|
|
|
1194 |
sptime(i)=time(1)+real(i-1)*dt
|
|
|
1195 |
enddo
|
|
|
1196 |
|
|
|
1197 |
c Do the spline interpolation
|
|
|
1198 |
do i=1,spn
|
|
|
1199 |
call splint(time,lon,y2ax,n,sptime(i),s)
|
|
|
1200 |
splon(i)=s
|
|
|
1201 |
enddo
|
|
|
1202 |
|
|
|
1203 |
c Change the time arrays back
|
|
|
1204 |
do i=1,spn
|
|
|
1205 |
sptime(i)=order*sptime(i)
|
|
|
1206 |
enddo
|
|
|
1207 |
do i=1,n
|
|
|
1208 |
time(i)=order*time(i)
|
|
|
1209 |
enddo
|
|
|
1210 |
|
|
|
1211 |
return
|
|
|
1212 |
end
|
|
|
1213 |
|
|
|
1214 |
c -------------------------------------------------------------
|
|
|
1215 |
c Basic routines for spline interpolation (Numerical Recipes)
|
|
|
1216 |
c -------------------------------------------------------------
|
|
|
1217 |
|
|
|
1218 |
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
|
|
|
1219 |
INTEGER n,NMAX
|
|
|
1220 |
REAL yp1,ypn,x(n),y(n),y2(n)
|
|
|
1221 |
PARAMETER (NMAX=500)
|
|
|
1222 |
INTEGER i,k
|
|
|
1223 |
REAL p,qn,sig,un,u(NMAX)
|
|
|
1224 |
if (yp1.gt..99e30) then
|
|
|
1225 |
y2(1)=0.
|
|
|
1226 |
u(1)=0.
|
|
|
1227 |
else
|
|
|
1228 |
y2(1)=-0.5
|
|
|
1229 |
u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
|
|
|
1230 |
endif
|
|
|
1231 |
do 11 i=2,n-1
|
|
|
1232 |
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
|
|
|
1233 |
p=sig*y2(i-1)+2.
|
|
|
1234 |
y2(i)=(sig-1.)/p
|
|
|
1235 |
u(i)=(6.*((y(i+1)-y(i))/(x(i+
|
|
|
1236 |
*1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
|
|
|
1237 |
*u(i-1))/p
|
|
|
1238 |
11 continue
|
|
|
1239 |
if (ypn.gt..99e30) then
|
|
|
1240 |
qn=0.
|
|
|
1241 |
un=0.
|
|
|
1242 |
else
|
|
|
1243 |
qn=0.5
|
|
|
1244 |
un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
|
|
|
1245 |
endif
|
|
|
1246 |
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
|
|
|
1247 |
do 12 k=n-1,1,-1
|
|
|
1248 |
y2(k)=y2(k)*y2(k+1)+u(k)
|
|
|
1249 |
12 continue
|
|
|
1250 |
return
|
|
|
1251 |
END
|
|
|
1252 |
|
|
|
1253 |
SUBROUTINE splint(xa,ya,y2a,n,x,y)
|
|
|
1254 |
INTEGER n
|
|
|
1255 |
REAL x,y,xa(n),y2a(n),ya(n)
|
|
|
1256 |
INTEGER k,khi,klo
|
|
|
1257 |
REAL a,b,h
|
|
|
1258 |
klo=1
|
|
|
1259 |
khi=n
|
|
|
1260 |
1 if (khi-klo.gt.1) then
|
|
|
1261 |
k=(khi+klo)/2
|
|
|
1262 |
if(xa(k).gt.x)then
|
|
|
1263 |
khi=k
|
|
|
1264 |
else
|
|
|
1265 |
klo=k
|
|
|
1266 |
endif
|
|
|
1267 |
goto 1
|
|
|
1268 |
endif
|
|
|
1269 |
h=xa(khi)-xa(klo)
|
| 11 |
michaesp |
1270 |
if (h.eq.0.) then
|
|
|
1271 |
print*,'bad xa input in splint'
|
|
|
1272 |
stop
|
|
|
1273 |
endif
|
| 3 |
michaesp |
1274 |
a=(xa(khi)-x)/h
|
|
|
1275 |
b=(x-xa(klo))/h
|
|
|
1276 |
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
|
|
|
1277 |
*2)/6.
|
|
|
1278 |
return
|
|
|
1279 |
END
|
|
|
1280 |
|
|
|
1281 |
c ********************************************************************
|
|
|
1282 |
c * INPUT / OUTPUT SUBROUTINES *
|
|
|
1283 |
c ********************************************************************
|
|
|
1284 |
|
|
|
1285 |
|
|
|
1286 |
c --------------------------------------------------------------------
|
|
|
1287 |
c Subroutines to write the CF netcdf output file
|
|
|
1288 |
c --------------------------------------------------------------------
|
|
|
1289 |
|
|
|
1290 |
subroutine writecdf2D_cf
|
|
|
1291 |
> (cdfname,varname,longname,unit,gridtype,clon,clat,
|
|
|
1292 |
> nlonlat,dlonlat,arr,time,dx,dy,xmin,ymin,nx,ny,
|
|
|
1293 |
> crefile,crevar,cretime)
|
|
|
1294 |
|
|
|
1295 |
c Create and write to the CF netcdf file <cdfname>. The variable
|
|
|
1296 |
c with name <varname> and with time <time> is written. The data
|
|
|
1297 |
c are in the two-dimensional array <arr>. The list <dx,dy,xmin,
|
|
|
1298 |
c ymin,nx,ny> specifies the output grid. The flags <crefile> and
|
|
|
1299 |
c <crevar> determine whether the file and/or the variable should
|
|
|
1300 |
c be created; correspondingly for the unlimited dimension <time>
|
|
|
1301 |
c with the flag <cretime>.
|
|
|
1302 |
|
|
|
1303 |
USE netcdf
|
|
|
1304 |
|
|
|
1305 |
IMPLICIT NONE
|
|
|
1306 |
|
|
|
1307 |
c Declaration of input parameters
|
|
|
1308 |
character*80 cdfname
|
|
|
1309 |
character*80 varname,longname,unit
|
|
|
1310 |
integer nx,ny
|
|
|
1311 |
real arr(nx,ny)
|
|
|
1312 |
real dx,dy,xmin,ymin
|
|
|
1313 |
real time
|
|
|
1314 |
integer crefile,crevar,cretime
|
|
|
1315 |
character*80 gridtype
|
|
|
1316 |
real clon,clat
|
|
|
1317 |
integer nlonlat
|
|
|
1318 |
real dlonlat
|
|
|
1319 |
|
|
|
1320 |
c Local variables
|
|
|
1321 |
integer ierr
|
|
|
1322 |
integer ncID
|
|
|
1323 |
integer LonDimId, varLonID
|
|
|
1324 |
integer LatDimID, varLatID
|
|
|
1325 |
integer TimeDimID, varTimeID
|
|
|
1326 |
real longitude(nx)
|
|
|
1327 |
real latitude (ny)
|
|
|
1328 |
real timeindex
|
|
|
1329 |
integer i
|
|
|
1330 |
integer nvars,varids(100)
|
|
|
1331 |
integer ndims,dimids(100)
|
|
|
1332 |
real timelist(1000)
|
|
|
1333 |
integer ntimes
|
|
|
1334 |
integer ind
|
|
|
1335 |
integer varID
|
|
|
1336 |
|
|
|
1337 |
c Quick an dirty solution for fieldname conflict
|
|
|
1338 |
if ( varname.eq.'time' ) varname = 'TIME'
|
|
|
1339 |
|
|
|
1340 |
c Initially set error to indicate no errors.
|
|
|
1341 |
ierr = 0
|
|
|
1342 |
|
|
|
1343 |
c ---- Create the netCDF - skip if <crefile=0> ----------------------
|
|
|
1344 |
if ( crefile.ne.1 ) goto 100
|
|
|
1345 |
|
|
|
1346 |
c Create the file
|
|
|
1347 |
ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
|
|
|
1348 |
|
|
|
1349 |
c Define dimensions
|
|
|
1350 |
ierr=nf90_def_dim(ncID,'longitude',nx , LonDimID )
|
|
|
1351 |
ierr=nf90_def_dim(ncID,'latitude' ,ny , LatDimID )
|
|
|
1352 |
ierr=nf90_def_dim(ncID,'time' ,nf90_unlimited, TimeDimID)
|
|
|
1353 |
|
|
|
1354 |
c Define coordinate Variables
|
|
|
1355 |
ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
|
|
|
1356 |
> (/ LonDimID /),varLonID)
|
|
|
1357 |
ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
|
|
|
1358 |
ierr = nf90_put_att(ncID, varLonID, "units" ,"degree_east")
|
|
|
1359 |
|
|
|
1360 |
ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
|
|
|
1361 |
> (/ LatDimID /),varLatID)
|
|
|
1362 |
ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
|
|
|
1363 |
ierr = nf90_put_att(ncID, varLatID, "units" ,"degree_north")
|
|
|
1364 |
|
|
|
1365 |
ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
|
|
|
1366 |
> (/ TimeDimID /), varTimeID)
|
|
|
1367 |
ierr = nf90_put_att(ncID, varTimeID, "axis", "T")
|
|
|
1368 |
ierr = nf90_put_att(ncID, varTimeID, "calendar", "standard")
|
|
|
1369 |
ierr = nf90_put_att(ncID, varTimeID, "long_name", "time")
|
|
|
1370 |
ierr = nf90_put_att(ncID, varTimeID, "units", "hours")
|
|
|
1371 |
|
|
|
1372 |
c Write global attributes
|
|
|
1373 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
|
|
|
1374 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',
|
|
|
1375 |
> 'Trajectory Densities')
|
|
|
1376 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source',
|
|
|
1377 |
> 'Lagranto Trajectories')
|
|
|
1378 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution',
|
|
|
1379 |
> 'ETH Zurich, IACETH')
|
|
|
1380 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'grid',trim(gridtype) )
|
| 39 |
michaesp |
1381 |
if ( gridtype.eq.'rotated' ) then
|
|
|
1382 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clon',clon )
|
|
|
1383 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clat',clat )
|
|
|
1384 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )
|
|
|
1385 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )
|
|
|
1386 |
endif
|
| 3 |
michaesp |
1387 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
|
|
|
1388 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
|
|
|
1389 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dx',dx )
|
|
|
1390 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dy',dy )
|
|
|
1391 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
|
|
|
1392 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
|
|
|
1393 |
|
|
|
1394 |
c Check whether the definition was successful
|
|
|
1395 |
ierr = nf90_enddef(ncID)
|
|
|
1396 |
if (ierr.gt.0) then
|
|
|
1397 |
print*, 'An error occurred while attempting to ',
|
|
|
1398 |
> 'finish definition mode.'
|
|
|
1399 |
stop
|
|
|
1400 |
endif
|
|
|
1401 |
|
|
|
1402 |
c Write coordinate data
|
| 39 |
michaesp |
1403 |
do i = 1,nx
|
|
|
1404 |
longitude(i) = xmin + real(i-1) * dx
|
|
|
1405 |
enddo
|
|
|
1406 |
do i = 1,ny
|
|
|
1407 |
latitude(i) = ymin + real(i-1) * dy
|
|
|
1408 |
enddo
|
| 3 |
michaesp |
1409 |
ierr = nf90_put_var(ncID,varLonID ,longitude)
|
|
|
1410 |
ierr = nf90_put_var(ncID,varLatID ,latitude )
|
|
|
1411 |
|
|
|
1412 |
c Close netCDF file
|
|
|
1413 |
ierr = nf90_close(ncID)
|
|
|
1414 |
|
|
|
1415 |
100 continue
|
|
|
1416 |
|
|
|
1417 |
c ---- Define a new variable - skip if <crevar=0> -----------------------
|
|
|
1418 |
|
|
|
1419 |
if ( crevar.ne.1 ) goto 110
|
|
|
1420 |
|
|
|
1421 |
c Open the file for read(write access
|
|
|
1422 |
ierr = nf90_open (trim(cdfname), NF90_WRITE , ncID)
|
|
|
1423 |
|
|
|
1424 |
c Get the IDs for dimensions
|
|
|
1425 |
ierr = nf90_inq_dimid(ncID,'longitude', LonDimID )
|
|
|
1426 |
ierr = nf90_inq_dimid(ncID,'latitude' , LatDimID )
|
|
|
1427 |
ierr = nf90_inq_dimid(ncID,'time' , TimeDimID)
|
|
|
1428 |
|
|
|
1429 |
c Enter define mode
|
|
|
1430 |
ierr = nf90_redef(ncID)
|
|
|
1431 |
|
|
|
1432 |
c Write definition and add attributes
|
|
|
1433 |
ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
|
| 13 |
michaesp |
1434 |
> (/ LonDimID, LatDimID, TimeDimID /),varID)
|
| 39 |
michaesp |
1435 |
|
| 3 |
michaesp |
1436 |
ierr = nf90_put_att(ncID, varID, "long_name" , longname )
|
|
|
1437 |
ierr = nf90_put_att(ncID, varID, "units" , unit )
|
|
|
1438 |
ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99 )
|
| 39 |
michaesp |
1439 |
|
| 3 |
michaesp |
1440 |
c Check whether definition was successful
|
|
|
1441 |
ierr = nf90_enddef(ncID)
|
|
|
1442 |
if (ierr.gt.0) then
|
|
|
1443 |
print*, 'An error occurred while attempting to ',
|
|
|
1444 |
> 'finish definition mode.'
|
|
|
1445 |
stop
|
|
|
1446 |
endif
|
| 13 |
michaesp |
1447 |
print*,trim(varname),' defined on ',trim(cdfname)
|
| 3 |
michaesp |
1448 |
|
|
|
1449 |
c Close netCDF file
|
|
|
1450 |
ierr = nf90_close(ncID)
|
|
|
1451 |
|
|
|
1452 |
110 continue
|
|
|
1453 |
|
|
|
1454 |
c ---- Create a new time (unlimited dimension) - skip if <cretime=0> ------
|
|
|
1455 |
|
|
|
1456 |
if ( cretime.ne.1 ) goto 120
|
|
|
1457 |
|
|
|
1458 |
c Open the file for read/write access
|
|
|
1459 |
ierr = nf90_open (trim(cdfname), NF90_WRITE, ncID)
|
|
|
1460 |
|
|
|
1461 |
c Get the list of times on the netCDF file
|
|
|
1462 |
ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
|
|
|
1463 |
if ( ierr.ne.0 ) then
|
|
|
1464 |
print*,'Time dimension is not defined on ',trim(cdfname),
|
|
|
1465 |
> ' .... Stop'
|
|
|
1466 |
stop
|
|
|
1467 |
endif
|
|
|
1468 |
ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
|
|
|
1469 |
ierr = nf90_inq_varid(ncID,'time', varTimeID)
|
|
|
1470 |
if ( ierr.ne.0 ) then
|
|
|
1471 |
print*,'Variable time is not defined on ',trim(cdfname),
|
|
|
1472 |
> ' ... Stop'
|
|
|
1473 |
stop
|
|
|
1474 |
endif
|
|
|
1475 |
ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
1476 |
|
|
|
1477 |
c Decide whether a new time must be written
|
|
|
1478 |
ind = 0
|
|
|
1479 |
do i=1,ntimes
|
|
|
1480 |
if ( time.eq.timelist(i) ) ind = i
|
|
|
1481 |
enddo
|
|
|
1482 |
|
|
|
1483 |
c Extend the time list if required
|
|
|
1484 |
if ( ind.eq.0 ) then
|
|
|
1485 |
ntimes = ntimes + 1
|
|
|
1486 |
timelist(ntimes) = time
|
|
|
1487 |
ierr = nf90_put_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
1488 |
endif
|
|
|
1489 |
|
|
|
1490 |
c Close netCDF file
|
|
|
1491 |
ierr = nf90_close(ncID)
|
|
|
1492 |
|
|
|
1493 |
120 continue
|
|
|
1494 |
|
|
|
1495 |
c ---- Write data --------------------------------------------------
|
|
|
1496 |
|
|
|
1497 |
c Open the file for read/write access
|
|
|
1498 |
ierr = nf90_open (trim(cdfname), NF90_WRITE , ncID)
|
|
|
1499 |
|
|
|
1500 |
c Get the varID
|
|
|
1501 |
ierr = nf90_inq_varid(ncID,varname, varID )
|
|
|
1502 |
if (ierr.ne.0) then
|
|
|
1503 |
print*,'Variable ',trim(varname),' is not defined on ',
|
|
|
1504 |
> trim(cdfname)
|
|
|
1505 |
stop
|
|
|
1506 |
endif
|
|
|
1507 |
|
|
|
1508 |
c Get the time index
|
|
|
1509 |
ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
|
|
|
1510 |
if ( ierr.ne.0 ) then
|
|
|
1511 |
print*,'Time dimension is not defined on ',trim(cdfname),
|
|
|
1512 |
> ' .... Stop'
|
|
|
1513 |
stop
|
|
|
1514 |
endif
|
|
|
1515 |
ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
|
|
|
1516 |
ierr = nf90_inq_varid(ncID,'time', varTimeID)
|
|
|
1517 |
if ( ierr.ne.0 ) then
|
|
|
1518 |
print*,'Variable time is not defined on ',trim(cdfname),
|
|
|
1519 |
> ' ... Stop'
|
|
|
1520 |
stop
|
|
|
1521 |
endif
|
|
|
1522 |
ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
1523 |
ind = 0
|
|
|
1524 |
do i=1,ntimes
|
|
|
1525 |
if ( time.eq.timelist(i) ) ind = i
|
|
|
1526 |
enddo
|
|
|
1527 |
if (ind.eq.0) then
|
|
|
1528 |
print*,'Time',time,' is not defined on the netCDF file',
|
|
|
1529 |
> trim(cdfname),' ... Stop'
|
|
|
1530 |
stop
|
|
|
1531 |
endif
|
|
|
1532 |
|
|
|
1533 |
c Write data block
|
|
|
1534 |
ierr = nf90_put_var(ncID,varID,arr,
|
|
|
1535 |
> start = (/ 1, 1, ind /),
|
|
|
1536 |
> count = (/ nx, ny, 1 /) )
|
|
|
1537 |
|
|
|
1538 |
c Check whether writing was successful
|
|
|
1539 |
ierr = nf90_close(ncID)
|
|
|
1540 |
if (ierr.ne.0) then
|
|
|
1541 |
write(*,*) trim(nf90_strerror(ierr))
|
|
|
1542 |
write(*,*) 'An error occurred while attempting to ',
|
|
|
1543 |
> 'close the netcdf file.'
|
|
|
1544 |
write(*,*) 'in clscdf_CF'
|
|
|
1545 |
endif
|
|
|
1546 |
|
|
|
1547 |
end
|
|
|
1548 |
|
|
|
1549 |
|
|
|
1550 |
c ********************************************************************************
|
|
|
1551 |
c * Transformation routine: LMSTOLM and PHSTOPH from library gm2em *
|
|
|
1552 |
c ********************************************************************************
|
|
|
1553 |
|
|
|
1554 |
REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
1555 |
C
|
|
|
1556 |
C**** LMSTOLM - FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
|
|
|
1557 |
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
1558 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
1559 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1560 |
C** AUFRUF : LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
1561 |
C** ENTRIES : KEINE
|
|
|
1562 |
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
|
|
|
1563 |
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
1564 |
C** IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
1565 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1566 |
C** VERSIONS-
|
|
|
1567 |
C** DATUM : 03.05.90
|
|
|
1568 |
C**
|
|
|
1569 |
C** EXTERNALS: KEINE
|
|
|
1570 |
C** EINGABE-
|
|
|
1571 |
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.
|
|
|
1572 |
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
|
|
|
1573 |
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS
|
|
|
1574 |
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS
|
|
|
1575 |
C** AUSGABE-
|
|
|
1576 |
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
|
|
|
1577 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
1578 |
C**
|
|
|
1579 |
C** COMMON-
|
|
|
1580 |
C** BLOECKE : KEINE
|
|
|
1581 |
C**
|
|
|
1582 |
C** FEHLERBE-
|
|
|
1583 |
C** HANDLUNG : KEINE
|
|
|
1584 |
C** VERFASSER: D.MAJEWSKI
|
|
|
1585 |
|
|
|
1586 |
REAL LAMS,PHIS,POLPHI,POLLAM
|
|
|
1587 |
|
|
|
1588 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
1589 |
|
|
|
1590 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
1591 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
1592 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
1593 |
ZPHIS = ZPIR18*PHIS
|
|
|
1594 |
ZLAMS = LAMS
|
|
|
1595 |
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
|
|
|
1596 |
ZLAMS = ZPIR18*ZLAMS
|
|
|
1597 |
|
|
|
1598 |
ZARG1 = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +
|
|
|
1599 |
1 ZCOSPOL* SIN(ZPHIS)) -
|
|
|
1600 |
2 COS(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)
|
|
|
1601 |
ZARG2 = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +
|
|
|
1602 |
1 ZCOSPOL* SIN(ZPHIS)) +
|
|
|
1603 |
2 SIN(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)
|
|
|
1604 |
IF (ABS(ZARG2).LT.1.E-30) THEN
|
|
|
1605 |
IF (ABS(ZARG1).LT.1.E-30) THEN
|
|
|
1606 |
LMSTOLM = 0.0
|
|
|
1607 |
ELSEIF (ZARG1.GT.0.) THEN
|
|
|
1608 |
LMSTOLAM = 90.0
|
|
|
1609 |
ELSE
|
|
|
1610 |
LMSTOLAM = -90.0
|
|
|
1611 |
ENDIF
|
|
|
1612 |
ELSE
|
|
|
1613 |
LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
|
|
|
1614 |
ENDIF
|
|
|
1615 |
|
|
|
1616 |
RETURN
|
|
|
1617 |
END
|
|
|
1618 |
|
|
|
1619 |
|
|
|
1620 |
REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
1621 |
C
|
|
|
1622 |
C**** PHSTOPH - FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
|
|
|
1623 |
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
1624 |
C**** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
1625 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1626 |
C** AUFRUF : PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
|
|
|
1627 |
C** ENTRIES : KEINE
|
|
|
1628 |
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
|
|
|
1629 |
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
1630 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
1631 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1632 |
C** VERSIONS-
|
|
|
1633 |
C** DATUM : 03.05.90
|
|
|
1634 |
C**
|
|
|
1635 |
C** EXTERNALS: KEINE
|
|
|
1636 |
C** EINGABE-
|
|
|
1637 |
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.
|
|
|
1638 |
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
|
|
|
1639 |
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS
|
|
|
1640 |
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS
|
|
|
1641 |
C** AUSGABE-
|
|
|
1642 |
C** PARAMETER: WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
|
|
|
1643 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
1644 |
C**
|
|
|
1645 |
C** COMMON-
|
|
|
1646 |
C** BLOECKE : KEINE
|
|
|
1647 |
C**
|
|
|
1648 |
C** FEHLERBE-
|
|
|
1649 |
C** HANDLUNG : KEINE
|
|
|
1650 |
C** VERFASSER: D.MAJEWSKI
|
|
|
1651 |
|
|
|
1652 |
REAL LAMS,PHIS,POLPHI,POLLAM
|
|
|
1653 |
|
|
|
1654 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
1655 |
|
|
|
1656 |
SINPOL = SIN(ZPIR18*POLPHI)
|
|
|
1657 |
COSPOL = COS(ZPIR18*POLPHI)
|
|
|
1658 |
ZPHIS = ZPIR18*PHIS
|
|
|
1659 |
ZLAMS = LAMS
|
|
|
1660 |
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
|
|
|
1661 |
ZLAMS = ZPIR18*ZLAMS
|
|
|
1662 |
ARG = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
|
|
|
1663 |
|
|
|
1664 |
PHSTOPH = ZRPI18*ASIN(ARG)
|
|
|
1665 |
|
|
|
1666 |
RETURN
|
|
|
1667 |
END
|
|
|
1668 |
|
|
|
1669 |
|
|
|
1670 |
REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
1671 |
C
|
|
|
1672 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
1673 |
C
|
|
|
1674 |
C**** LMTOLMS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
|
|
|
1675 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
1676 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
1677 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1678 |
C** AUFRUF : LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
1679 |
C** ENTRIES : KEINE
|
|
|
1680 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
|
|
|
1681 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
1682 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
1683 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1684 |
C** VERSIONS-
|
|
|
1685 |
C** DATUM : 03.05.90
|
|
|
1686 |
C**
|
|
|
1687 |
C** EXTERNALS: KEINE
|
|
|
1688 |
C** EINGABE-
|
|
|
1689 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
1690 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
1691 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
1692 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
1693 |
C** AUSGABE-
|
|
|
1694 |
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
|
|
|
1695 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
1696 |
C**
|
|
|
1697 |
C** COMMON-
|
|
|
1698 |
C** BLOECKE : KEINE
|
|
|
1699 |
C**
|
|
|
1700 |
C** FEHLERBE-
|
|
|
1701 |
C** HANDLUNG : KEINE
|
|
|
1702 |
C** VERFASSER: G. DE MORSIER
|
|
|
1703 |
|
|
|
1704 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
1705 |
|
|
|
1706 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
1707 |
|
|
|
1708 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
1709 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
1710 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
1711 |
ZPHI = ZPIR18*PHI
|
|
|
1712 |
ZLAM = LAM
|
|
|
1713 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
1714 |
ZLAM = ZPIR18*ZLAM
|
|
|
1715 |
|
|
|
1716 |
ZARG1 = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
|
|
|
1717 |
ZARG2 = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
|
|
|
1718 |
IF (ABS(ZARG2).LT.1.E-30) THEN
|
|
|
1719 |
IF (ABS(ZARG1).LT.1.E-30) THEN
|
|
|
1720 |
LMTOLMS = 0.0
|
|
|
1721 |
ELSEIF (ZARG1.GT.0.) THEN
|
|
|
1722 |
LMTOLMS = 90.0
|
|
|
1723 |
ELSE
|
|
|
1724 |
LMTOLMS = -90.0
|
|
|
1725 |
ENDIF
|
|
|
1726 |
ELSE
|
|
|
1727 |
LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
|
|
|
1728 |
ENDIF
|
|
|
1729 |
|
|
|
1730 |
RETURN
|
|
|
1731 |
END
|
|
|
1732 |
|
|
|
1733 |
|
|
|
1734 |
REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
1735 |
C
|
|
|
1736 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
1737 |
C
|
|
|
1738 |
C**** PHTOPHS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
|
|
|
1739 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
1740 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
1741 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1742 |
C** AUFRUF : PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
1743 |
C** ENTRIES : KEINE
|
|
|
1744 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
|
|
|
1745 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
1746 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
1747 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
1748 |
C** VERSIONS-
|
|
|
1749 |
C** DATUM : 03.05.90
|
|
|
1750 |
C**
|
|
|
1751 |
C** EXTERNALS: KEINE
|
|
|
1752 |
C** EINGABE-
|
|
|
1753 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
1754 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
1755 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
1756 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
1757 |
C** AUSGABE-
|
|
|
1758 |
C** PARAMETER: ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
|
|
|
1759 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
1760 |
C**
|
|
|
1761 |
C** COMMON-
|
|
|
1762 |
C** BLOECKE : KEINE
|
|
|
1763 |
C**
|
|
|
1764 |
C** FEHLERBE-
|
|
|
1765 |
C** HANDLUNG : KEINE
|
|
|
1766 |
C** VERFASSER: G. DE MORSIER
|
|
|
1767 |
|
|
|
1768 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
1769 |
|
|
|
1770 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
1771 |
|
|
|
1772 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
1773 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
1774 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
1775 |
ZPHI = ZPIR18*PHI
|
|
|
1776 |
ZLAM = LAM
|
|
|
1777 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
1778 |
ZLAM = ZPIR18*ZLAM
|
|
|
1779 |
ZARG = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
|
|
|
1780 |
|
|
|
1781 |
PHTOPHS = ZRPI18*ASIN(ZARG)
|
|
|
1782 |
|
|
|
1783 |
RETURN
|
|
|
1784 |
END
|
| 11 |
michaesp |
1785 |
|
|
|
1786 |
c ------------------------------------------------------------------
|
|
|
1787 |
c Compute Cos/Sin of an argument in Degree instead of Radian
|
|
|
1788 |
c ------------------------------------------------------------------
|
|
|
1789 |
|
|
|
1790 |
real function cosd(arg)
|
|
|
1791 |
|
|
|
1792 |
real,intent(IN) :: arg
|
| 29 |
michaesp |
1793 |
real,parameter :: grad2rad=3.1415926/180.
|
| 11 |
michaesp |
1794 |
cosd=cos(arg*grad2rad)
|
|
|
1795 |
return
|
|
|
1796 |
end
|
|
|
1797 |
|
|
|
1798 |
real function sind(arg)
|
|
|
1799 |
|
|
|
1800 |
real,intent(IN) :: arg
|
| 29 |
michaesp |
1801 |
real,parameter :: grad2rad=3.1415926/180.
|
| 11 |
michaesp |
1802 |
sind=sin(arg*grad2rad)
|
|
|
1803 |
return
|
|
|
1804 |
end
|
|
|
1805 |
|
|
|
1806 |
|