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