| 3 |
michaesp |
1 |
PROGRAM rotate_lalo
|
|
|
2 |
|
|
|
3 |
c *******************************************************************
|
|
|
4 |
c * Rotate to a geographical latitude/longitude grid *
|
|
|
5 |
c * Michael Sprenger / Autumn 2006 *
|
|
|
6 |
c *******************************************************************
|
|
|
7 |
|
|
|
8 |
implicit none
|
|
|
9 |
|
|
|
10 |
c -------------------------------------------------------------------
|
|
|
11 |
c Declaration of parameters and variables
|
|
|
12 |
c -------------------------------------------------------------------
|
|
|
13 |
|
|
|
14 |
c Specification of input parameters
|
|
|
15 |
real clon,clat,rotation
|
|
|
16 |
integer nx,ny
|
|
|
17 |
real dx,dy,xmin,ymin
|
|
|
18 |
integer nfield
|
|
|
19 |
character*80 fieldname(100)
|
|
|
20 |
|
|
|
21 |
c Numerical and physical parameters
|
|
|
22 |
real degrad
|
|
|
23 |
parameter (degrad=0.0174532925)
|
|
|
24 |
real omegaerd
|
|
|
25 |
parameter (omegaerd=7.292e-5 )
|
|
|
26 |
real eps
|
|
|
27 |
parameter (eps=0.001)
|
|
|
28 |
|
|
|
29 |
c Variables for input Z file : height level
|
|
|
30 |
character*80 in_cfn
|
|
|
31 |
real in_varmin(4),in_varmax(4),in_stag(4)
|
|
|
32 |
integer in_vardim(4)
|
|
|
33 |
real in_mdv
|
|
|
34 |
integer in_ndim
|
|
|
35 |
integer in_nx,in_ny,in_nz
|
|
|
36 |
real in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy
|
|
|
37 |
integer in_ntimes
|
|
|
38 |
real in_aklev(500),in_bklev(500)
|
|
|
39 |
real in_aklay(500),in_bklay(500)
|
|
|
40 |
real in_time
|
|
|
41 |
real in_pollon,in_pollat
|
|
|
42 |
integer in_nvars
|
|
|
43 |
character*80 in_vnam(100)
|
|
|
44 |
integer in_idate(5)
|
|
|
45 |
real,allocatable, dimension (:,:,:) :: in_field3
|
|
|
46 |
real,allocatable, dimension (:,:,:) :: in_vect1
|
|
|
47 |
real,allocatable, dimension (:,:,:) :: in_vect2
|
|
|
48 |
|
|
|
49 |
c Variables for output Z file : height level
|
|
|
50 |
character*80 out_cfn
|
|
|
51 |
real out_varmin(4),out_varmax(4),out_stag(4)
|
|
|
52 |
integer out_vardim(4)
|
|
|
53 |
real out_mdv
|
|
|
54 |
integer out_ndim
|
|
|
55 |
integer out_nx,out_ny,out_nz
|
|
|
56 |
real out_xmin,out_xmax,out_ymin,out_ymax,out_dx,out_dy
|
|
|
57 |
integer out_ntimes
|
|
|
58 |
real out_aklev(500),out_bklev(500)
|
|
|
59 |
real out_aklay(500),out_bklay(500)
|
|
|
60 |
real out_time
|
|
|
61 |
real out_pollon,out_pollat
|
|
|
62 |
integer out_idate(5)
|
|
|
63 |
real,allocatable, dimension (:,:,:) :: out_field3
|
|
|
64 |
real,allocatable, dimension (:,:,:) :: out_vect1
|
|
|
65 |
real,allocatable, dimension (:,:,:) :: out_vect2
|
|
|
66 |
real,allocatable, dimension (:,:) :: out_lat,out_lon
|
|
|
67 |
real,allocatable, dimension (:,:) :: out_x,out_y
|
|
|
68 |
real,allocatable, dimension (:,:) :: out_coriol
|
|
|
69 |
|
|
|
70 |
c Auxiliary variables
|
|
|
71 |
integer ifield
|
|
|
72 |
integer i,j,k
|
|
|
73 |
integer cdfid,cstid
|
|
|
74 |
character*80 cfn
|
|
|
75 |
integer stat,ierr,isok
|
|
|
76 |
real time
|
|
|
77 |
character*80 varname,cdfname,varname1,varname2
|
|
|
78 |
integer idate(5),stdate(5),datar(14)
|
|
|
79 |
integer rotmode
|
|
|
80 |
integer i1,i2,i3,i4
|
|
|
81 |
integer isvector
|
|
|
82 |
real lat
|
|
|
83 |
character*80 name
|
|
|
84 |
|
|
|
85 |
c Externals
|
|
|
86 |
real sdis
|
|
|
87 |
external sdis
|
|
|
88 |
|
|
|
89 |
c -------------------------------------------------------------------
|
|
|
90 |
c Preparations
|
|
|
91 |
c -------------------------------------------------------------------
|
|
|
92 |
|
|
|
93 |
print*,'*********************************************************'
|
|
|
94 |
print*,'* rotate_lalo *'
|
|
|
95 |
print*,'*********************************************************'
|
|
|
96 |
|
|
|
97 |
c Read parameter file
|
|
|
98 |
open(10,file='fort.10')
|
|
|
99 |
|
|
|
100 |
read(10,*) in_cfn
|
|
|
101 |
read(10,*) out_cfn
|
|
|
102 |
|
|
|
103 |
read(10,*) name,nx
|
|
|
104 |
if ( trim(name).ne.'GEO_NX' ) stop
|
|
|
105 |
read(10,*) name,ny
|
|
|
106 |
if ( trim(name).ne.'GEO_NY' ) stop
|
|
|
107 |
read(10,*) name,dx
|
|
|
108 |
if ( trim(name).ne.'GEO_DX' ) stop
|
|
|
109 |
read(10,*) name,dy
|
|
|
110 |
if ( trim(name).ne.'GEO_DY' ) stop
|
|
|
111 |
read(10,*) name,xmin
|
|
|
112 |
if ( trim(name).ne.'GEO_XMIN') stop
|
|
|
113 |
read(10,*) name,ymin
|
|
|
114 |
if ( trim(name).ne.'GEO_YMIN') stop
|
|
|
115 |
read(10,*) name,clon
|
|
|
116 |
if ( trim(name).ne.'CLON' ) stop
|
|
|
117 |
read(10,*) name,clat
|
|
|
118 |
if ( trim(name).ne.'CLAT' ) stop
|
|
|
119 |
read(10,*) name,rotation
|
|
|
120 |
if ( trim(name).ne.'CROT' ) stop
|
|
|
121 |
|
|
|
122 |
read(10,*) nfield
|
|
|
123 |
do i=1,nfield
|
|
|
124 |
read(10,*) fieldname(i)
|
|
|
125 |
enddo
|
|
|
126 |
close(10)
|
|
|
127 |
|
|
|
128 |
print*,clon,clat,rotation
|
|
|
129 |
print*,nx,ny,dx,dy,xmin,ymin
|
|
|
130 |
print*,trim(in_cfn),' -> ',trim(out_cfn)
|
|
|
131 |
|
|
|
132 |
c Get grid description for input Z file : height level
|
|
|
133 |
call cdfopn(in_cfn,cdfid,ierr)
|
|
|
134 |
if (ierr.ne.0) goto 998
|
|
|
135 |
call getcfn(cdfid,cfn,ierr)
|
|
|
136 |
if (ierr.ne.0) goto 998
|
|
|
137 |
call cdfopn(cfn,cstid,ierr)
|
|
|
138 |
if (ierr.ne.0) goto 998
|
|
|
139 |
call getvars(cdfid,in_nvars,in_vnam,ierr)
|
|
|
140 |
if (ierr.ne.0) goto 998
|
|
|
141 |
isok=0
|
|
|
142 |
varname=fieldname(1)
|
|
|
143 |
call check_varok (isok,varname,in_vnam,in_nvars)
|
|
|
144 |
if (isok.eq.0) goto 998
|
|
|
145 |
call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim,
|
|
|
146 |
> in_varmin,in_varmax,in_stag,ierr)
|
|
|
147 |
if (ierr.ne.0) goto 998
|
|
|
148 |
in_nx =in_vardim(1)
|
|
|
149 |
in_ny =in_vardim(2)
|
|
|
150 |
in_xmin=in_varmin(1)
|
|
|
151 |
in_ymin=in_varmin(2)
|
|
|
152 |
call getlevs(cstid,in_nz,in_aklev,in_bklev,in_aklay,in_bklay,ierr)
|
|
|
153 |
call getgrid(cstid,in_dx,in_dy,ierr)
|
|
|
154 |
in_xmax=in_xmin+real(in_nx-1)*in_dx
|
|
|
155 |
in_ymax=in_ymin+real(in_ny-1)*in_dy
|
|
|
156 |
call gettimes(cdfid,in_time,in_ntimes,ierr)
|
|
|
157 |
call getstart(cstid,in_idate,ierr)
|
|
|
158 |
call getpole(cstid,in_pollon,in_pollat,ierr)
|
|
|
159 |
call clscdf(cstid,ierr)
|
|
|
160 |
call clscdf(cdfid,ierr)
|
|
|
161 |
|
|
|
162 |
c Set grid description for output file : height level
|
|
|
163 |
out_vardim(1) = nx
|
|
|
164 |
out_vardim(2) = ny
|
|
|
165 |
out_vardim(3) = in_nz
|
|
|
166 |
out_varmin(1) = xmin
|
|
|
167 |
out_varmin(2) = ymin
|
|
|
168 |
out_varmin(3) = in_varmin(3)
|
|
|
169 |
out_varmax(1) = xmin+real(nx-1)*dx
|
|
|
170 |
out_varmax(2) = ymin+real(ny-1)*dy
|
|
|
171 |
out_varmax(3) = in_varmax(3)
|
|
|
172 |
do i=1,in_nz
|
|
|
173 |
out_aklay(i) = in_aklay(i)
|
|
|
174 |
out_bklay(i) = in_bklay(i)
|
|
|
175 |
out_aklev(i) = in_aklev(i)
|
|
|
176 |
out_bklev(i) = in_bklev(i)
|
|
|
177 |
enddo
|
|
|
178 |
out_dx = dx
|
|
|
179 |
out_dy = dy
|
|
|
180 |
out_time = in_time
|
|
|
181 |
out_ntimes = in_ntimes
|
|
|
182 |
out_ndim = 4
|
|
|
183 |
out_mdv = in_mdv
|
|
|
184 |
out_nx = out_vardim(1)
|
|
|
185 |
out_ny = out_vardim(2)
|
|
|
186 |
out_nz = out_vardim(3)
|
|
|
187 |
out_xmin = out_varmin(1)
|
|
|
188 |
out_ymin = out_varmin(2)
|
|
|
189 |
out_pollon = 0.
|
|
|
190 |
out_pollat = 90.
|
|
|
191 |
do i=1,5
|
|
|
192 |
out_idate(i) = in_idate(i)
|
|
|
193 |
enddo
|
|
|
194 |
|
|
|
195 |
c Allocate memory for all fields
|
|
|
196 |
allocate(in_field3(in_nx,in_ny,in_nz),stat=stat)
|
|
|
197 |
if (stat.ne.0) print*,'*** error allocating array in_field3 ***'
|
|
|
198 |
allocate(in_vect1(in_nx,in_ny,in_nz),stat=stat)
|
|
|
199 |
if (stat.ne.0) print*,'*** error allocating array in_vect1 ***'
|
|
|
200 |
allocate(in_vect2(in_nx,in_ny,in_nz),stat=stat)
|
|
|
201 |
if (stat.ne.0) print*,'*** error allocating array in_vect2 ***'
|
|
|
202 |
allocate(out_field3(out_nx,out_ny,out_nz),stat=stat)
|
|
|
203 |
if (stat.ne.0) print*,'*** error allocating array out_field3 ***'
|
|
|
204 |
allocate(out_vect1(out_nx,out_ny,out_nz),stat=stat)
|
|
|
205 |
if (stat.ne.0) print*,'*** error allocating array out_vect1 ***'
|
|
|
206 |
allocate(out_vect2(out_nx,out_ny,out_nz),stat=stat)
|
|
|
207 |
if (stat.ne.0) print*,'*** error allocating array out_vect2 ***'
|
|
|
208 |
allocate(out_lat(out_nx,out_ny),stat=stat)
|
|
|
209 |
if (stat.ne.0) print*,'*** error allocating array out_lat ***'
|
|
|
210 |
allocate(out_lon(out_nx,out_ny),stat=stat)
|
|
|
211 |
if (stat.ne.0) print*,'*** error allocating array out_lon ***'
|
|
|
212 |
allocate(out_x(out_nx,out_ny),stat=stat)
|
|
|
213 |
if (stat.ne.0) print*,'*** error allocating array out_x ***'
|
|
|
214 |
allocate(out_y(out_nx,out_ny),stat=stat)
|
|
|
215 |
if (stat.ne.0) print*,'*** error allocating array out_y ***'
|
|
|
216 |
allocate(out_coriol(out_nx,out_ny),stat=stat)
|
|
|
217 |
if (stat.ne.0) print*,'*** error allocating array out_coriol ***'
|
|
|
218 |
|
|
|
219 |
c Create constants file and output file (if requested)
|
|
|
220 |
datar(1)=out_nx
|
|
|
221 |
datar(2)=out_ny
|
|
|
222 |
datar(3)=int(1000.*out_varmax(2))
|
|
|
223 |
datar(4)=int(1000.*out_varmin(1))
|
|
|
224 |
datar(5)=int(1000.*out_varmin(2))
|
|
|
225 |
datar(6)=int(1000.*out_varmax(1))
|
|
|
226 |
datar(7)=int(1000.*out_dx)
|
|
|
227 |
datar(8)=int(1000.*out_dy)
|
|
|
228 |
datar(9)=out_nz
|
|
|
229 |
datar(10)=1
|
|
|
230 |
datar(11)=1
|
|
|
231 |
datar(12)=0
|
|
|
232 |
datar(13)=int(1000.*out_pollon)
|
|
|
233 |
datar(14)=int(1000.*out_pollat)
|
|
|
234 |
|
|
|
235 |
cfn=trim(out_cfn)//'_cst'
|
|
|
236 |
call wricst(cfn,datar,out_aklev,out_bklev,
|
|
|
237 |
> out_aklay,out_bklay,out_idate)
|
|
|
238 |
|
|
|
239 |
call crecdf(trim(out_cfn),cdfid,out_varmin,out_varmax,
|
|
|
240 |
> out_ndim,trim(cfn),ierr)
|
|
|
241 |
if (ierr.ne.0) goto 997
|
|
|
242 |
call clscdf(cdfid,ierr)
|
|
|
243 |
|
|
|
244 |
c -----------------------------------------------------------------
|
|
|
245 |
c Loop through all fields - rotate to new grid
|
|
|
246 |
c -----------------------------------------------------------------
|
|
|
247 |
|
|
|
248 |
do ifield=1,nfield
|
|
|
249 |
|
|
|
250 |
c Check if scalar or vectorial field (X for scalar, U.V for vector)
|
|
|
251 |
varname=fieldname(ifield)
|
|
|
252 |
i1=1
|
|
|
253 |
i2=1
|
|
|
254 |
i3=0
|
|
|
255 |
i4=0
|
|
|
256 |
100 if (varname(i1:i1).eq.' ') then
|
|
|
257 |
i1=i1+1
|
|
|
258 |
goto 100
|
|
|
259 |
endif
|
|
|
260 |
i2=i1
|
|
|
261 |
101 if ((varname(i2:i2).ne.' ').and.(varname(i2:i2).ne.'.')) then
|
|
|
262 |
i2=i2+1
|
|
|
263 |
goto 101
|
|
|
264 |
endif
|
|
|
265 |
if (varname(i2:i2).eq.'.') then
|
|
|
266 |
i3=i2+1
|
|
|
267 |
102 if (varname(i3:i3).eq.' ') then
|
|
|
268 |
i3=i3+1
|
|
|
269 |
goto 102
|
|
|
270 |
endif
|
|
|
271 |
i4=i3
|
|
|
272 |
104 if (varname(i4:i4).ne.' ') then
|
|
|
273 |
i4=i4+1
|
|
|
274 |
goto 104
|
|
|
275 |
endif
|
|
|
276 |
endif
|
|
|
277 |
if ((i3.ne.0).and.(i4.ne.0)) then
|
|
|
278 |
isvector=1
|
|
|
279 |
i2=i2-1
|
|
|
280 |
varname1=varname(i1:i2)
|
|
|
281 |
i4=i4-1
|
|
|
282 |
varname2=varname(i3:i4)
|
|
|
283 |
print*,'Rotating vector : ',
|
|
|
284 |
> trim(varname1),' / ',trim(varname2)
|
|
|
285 |
else
|
|
|
286 |
isvector=0
|
|
|
287 |
i2=i2-1
|
|
|
288 |
varname=varname(i1:i2)
|
|
|
289 |
print*,'Rotating scalar : ',trim(varname)
|
|
|
290 |
endif
|
|
|
291 |
|
|
|
292 |
c Rotate a scalar field
|
|
|
293 |
if (isvector.eq.0) then
|
|
|
294 |
|
|
|
295 |
c Read input field
|
|
|
296 |
call cdfopn(in_cfn,cdfid,ierr)
|
|
|
297 |
if (ierr.ne.0) goto 998
|
|
|
298 |
call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim,
|
|
|
299 |
> in_varmin,in_varmax,in_stag,ierr)
|
|
|
300 |
in_nz=in_vardim(3)
|
|
|
301 |
if (ierr.ne.0) goto 998
|
|
|
302 |
call getdat(cdfid,varname,in_time,0,in_field3,ierr)
|
|
|
303 |
if (ierr.ne.0) goto 998
|
|
|
304 |
call clscdf(cdfid,ierr)
|
|
|
305 |
|
|
|
306 |
c Rotate to new coordinate system
|
|
|
307 |
out_nz=in_nz
|
|
|
308 |
call getenvir (clon,clat,rotation,0,
|
|
|
309 |
> in_field3,in_dx,in_dy,
|
|
|
310 |
> in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
|
|
|
311 |
> out_field3,out_dx,out_dy,
|
|
|
312 |
> out_xmin,out_ymin,out_nx,out_ny,out_nz)
|
|
|
313 |
|
|
|
314 |
c Write output scalar field
|
|
|
315 |
call cdfwopn(trim(out_cfn),cdfid,ierr)
|
|
|
316 |
if (ierr.ne.0) goto 997
|
|
|
317 |
out_vardim(3)=out_nz
|
|
|
318 |
call putdef(cdfid,varname,4,out_mdv,out_vardim,
|
|
|
319 |
> out_varmin,out_varmax,out_stag,ierr)
|
|
|
320 |
if (ierr.ne.0) goto 997
|
|
|
321 |
call putdat(cdfid,varname,out_time,0,out_field3,ierr)
|
|
|
322 |
if (ierr.ne.0) goto 997
|
|
|
323 |
call clscdf(cdfid,ierr)
|
|
|
324 |
|
|
|
325 |
c Rotate vector field
|
|
|
326 |
else if (isvector.eq.1) then
|
|
|
327 |
|
|
|
328 |
c Read input vector field
|
|
|
329 |
call cdfopn(in_cfn,cdfid,ierr)
|
|
|
330 |
if (ierr.ne.0) goto 998
|
|
|
331 |
call getdef(cdfid,varname1,in_ndim,in_mdv,in_vardim,
|
|
|
332 |
> in_varmin,in_varmax,in_stag,ierr)
|
|
|
333 |
in_nz=in_vardim(3)
|
|
|
334 |
if (ierr.ne.0) goto 998
|
|
|
335 |
call getdat(cdfid,varname1,in_time,0,in_vect1,ierr)
|
|
|
336 |
if (ierr.ne.0) goto 998
|
|
|
337 |
call getdef(cdfid,varname2,in_ndim,in_mdv,in_vardim,
|
|
|
338 |
> in_varmin,in_varmax,in_stag,ierr)
|
|
|
339 |
in_nz=in_vardim(3)
|
|
|
340 |
if (ierr.ne.0) goto 998
|
|
|
341 |
call getdat(cdfid,varname2,in_time,0,in_vect2,ierr)
|
|
|
342 |
if (ierr.ne.0) goto 998
|
|
|
343 |
call clscdf(cdfid,ierr)
|
|
|
344 |
|
|
|
345 |
c Get new vector component in x direction
|
|
|
346 |
out_nz=in_nz
|
|
|
347 |
call getenvir (clon,clat,rotation,1,
|
|
|
348 |
> in_vect1,in_dx,in_dy,
|
|
|
349 |
> in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
|
|
|
350 |
> out_field3,out_dx,out_dy,
|
|
|
351 |
> out_xmin,out_ymin,out_nx,out_ny,out_nz)
|
|
|
352 |
do i=1,out_nx
|
|
|
353 |
do j=1,out_ny
|
|
|
354 |
do k=1,out_nz
|
|
|
355 |
out_vect1(i,j,k)=out_field3(i,j,k)
|
|
|
356 |
enddo
|
|
|
357 |
enddo
|
|
|
358 |
enddo
|
|
|
359 |
call getenvir (clon,clat,rotation,2,
|
|
|
360 |
> in_vect2,in_dx,in_dy,
|
|
|
361 |
> in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
|
|
|
362 |
> out_field3,out_dx,out_dy,
|
|
|
363 |
> out_xmin,out_ymin,out_nx,out_ny,out_nz)
|
|
|
364 |
do i=1,out_nx
|
|
|
365 |
do j=1,out_ny
|
|
|
366 |
do k=1,out_nz
|
|
|
367 |
if ( (abs(out_vect1 (i,j,k)-out_mdv).gt.eps).and.
|
|
|
368 |
> (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then
|
|
|
369 |
out_vect1(i,j,k)=out_vect1(i,j,k)-
|
|
|
370 |
> out_field3(i,j,k)
|
|
|
371 |
endif
|
|
|
372 |
enddo
|
|
|
373 |
enddo
|
|
|
374 |
enddo
|
|
|
375 |
|
|
|
376 |
c Get new vector component in y direction
|
|
|
377 |
out_nz=in_nz
|
|
|
378 |
call getenvir (clon,clat,rotation,2,
|
|
|
379 |
> in_vect1,in_dx,in_dy,
|
|
|
380 |
> in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
|
|
|
381 |
> out_field3,out_dx,out_dy,
|
|
|
382 |
> out_xmin,out_ymin,out_nx,out_ny,out_nz)
|
|
|
383 |
do i=1,out_nx
|
|
|
384 |
do j=1,out_ny
|
|
|
385 |
do k=1,out_nz
|
|
|
386 |
out_vect2(i,j,k)=out_field3(i,j,k)
|
|
|
387 |
enddo
|
|
|
388 |
enddo
|
|
|
389 |
enddo
|
|
|
390 |
call getenvir (clon,clat,rotation,1,
|
|
|
391 |
> in_vect2,in_dx,in_dy,
|
|
|
392 |
> in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
|
|
|
393 |
> out_field3,out_dx,out_dy,
|
|
|
394 |
> out_xmin,out_ymin,out_nx,out_ny,out_nz)
|
|
|
395 |
do i=1,out_nx
|
|
|
396 |
do j=1,out_ny
|
|
|
397 |
do k=1,out_nz
|
|
|
398 |
if ( (abs(out_vect2 (i,j,k)-out_mdv).gt.eps).and.
|
|
|
399 |
> (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then
|
|
|
400 |
out_vect2(i,j,k)=out_vect2(i,j,k)+
|
|
|
401 |
> out_field3(i,j,k)
|
|
|
402 |
endif
|
|
|
403 |
enddo
|
|
|
404 |
enddo
|
|
|
405 |
enddo
|
|
|
406 |
|
|
|
407 |
c Write output vector field
|
|
|
408 |
call cdfwopn(trim(out_cfn),cdfid,ierr)
|
|
|
409 |
if (ierr.ne.0) goto 997
|
|
|
410 |
out_vardim(3)=out_nz
|
|
|
411 |
call putdef(cdfid,varname1,4,out_mdv,out_vardim,
|
|
|
412 |
> out_varmin,out_varmax,out_stag,ierr)
|
|
|
413 |
if (ierr.ne.0) goto 997
|
|
|
414 |
call putdat(cdfid,varname1,out_time,0,out_vect1,ierr)
|
|
|
415 |
if (ierr.ne.0) goto 997
|
|
|
416 |
call putdef(cdfid,varname2,4,out_mdv,out_vardim,
|
|
|
417 |
> out_varmin,out_varmax,out_stag,ierr)
|
|
|
418 |
if (ierr.ne.0) goto 997
|
|
|
419 |
call putdat(cdfid,varname2,out_time,0,out_vect2,ierr)
|
|
|
420 |
if (ierr.ne.0) goto 997
|
|
|
421 |
call clscdf(cdfid,ierr)
|
|
|
422 |
|
|
|
423 |
endif
|
|
|
424 |
|
|
|
425 |
enddo
|
|
|
426 |
|
|
|
427 |
c -----------------------------------------------------------------
|
|
|
428 |
c Write additional fields: latitude, longitude, Coriolis parameter
|
|
|
429 |
c -----------------------------------------------------------------
|
|
|
430 |
|
|
|
431 |
c Open the output file
|
|
|
432 |
call cdfwopn(trim(out_cfn),cdfid,ierr)
|
|
|
433 |
if (ierr.ne.0) goto 997
|
|
|
434 |
|
|
|
435 |
c Geographical latitude
|
|
|
436 |
varname='RLAT'
|
|
|
437 |
print*,'Rotating scalar : ',trim(varname)
|
|
|
438 |
do i=1,in_nx
|
|
|
439 |
do j=1,in_ny
|
|
|
440 |
in_field3(i,j,1)=in_ymin+real(j-1)*in_dy
|
|
|
441 |
enddo
|
|
|
442 |
enddo
|
|
|
443 |
call getenvir (clon,clat,rotation,0,
|
|
|
444 |
> in_field3,in_dx,in_dy,
|
|
|
445 |
> in_xmin,in_ymin,in_nx,in_ny,1,in_mdv,
|
|
|
446 |
> out_lat,out_dx,out_dy,
|
|
|
447 |
> out_xmin,out_ymin,out_nx,out_ny,1)
|
|
|
448 |
out_vardim(3)=1
|
|
|
449 |
call putdef(cdfid,varname,4,out_mdv,out_vardim,
|
|
|
450 |
> out_varmin,out_varmax,out_stag,ierr)
|
|
|
451 |
if (ierr.ne.0) goto 997
|
|
|
452 |
call putdat(cdfid,varname,out_time,0,out_lat,ierr)
|
|
|
453 |
if (ierr.ne.0) goto 997
|
|
|
454 |
|
|
|
455 |
c Geographical longitude
|
|
|
456 |
varname='RLON'
|
|
|
457 |
print*,'Rotating scalar : ',trim(varname)
|
|
|
458 |
do i=1,in_nx
|
|
|
459 |
do j=1,in_ny
|
|
|
460 |
in_field3(i,j,1)=in_xmin+real(i-1)*in_dx
|
|
|
461 |
enddo
|
|
|
462 |
enddo
|
|
|
463 |
call getenvir (clon,clat,rotation,0,
|
|
|
464 |
> in_field3,in_dx,in_dy,
|
|
|
465 |
> in_xmin,in_ymin,in_nx,in_ny,1,in_mdv,
|
|
|
466 |
> out_lon,out_dx,out_dy,
|
|
|
467 |
> out_xmin,out_ymin,out_nx,out_ny,1)
|
|
|
468 |
out_vardim(3)=1
|
|
|
469 |
call putdef(cdfid,varname,4,out_mdv,out_vardim,
|
|
|
470 |
> out_varmin,out_varmax,out_stag,ierr)
|
|
|
471 |
if (ierr.ne.0) goto 997
|
|
|
472 |
call putdat(cdfid,varname,out_time,0,out_lon,ierr)
|
|
|
473 |
if (ierr.ne.0) goto 997
|
|
|
474 |
|
|
|
475 |
|
|
|
476 |
|
|
|
477 |
c Close output file
|
|
|
478 |
call clscdf(cdfid,ierr)
|
|
|
479 |
|
|
|
480 |
c -----------------------------------------------------------------
|
|
|
481 |
c Exception handling and format specs
|
|
|
482 |
c -----------------------------------------------------------------
|
|
|
483 |
|
|
|
484 |
stop
|
|
|
485 |
|
|
|
486 |
998 print*,'Z: Problems with input rotated grid'
|
|
|
487 |
stop
|
|
|
488 |
|
|
|
489 |
997 print*,'Z: Problems with output to lat/lon grid'
|
|
|
490 |
stop
|
|
|
491 |
|
|
|
492 |
end
|
|
|
493 |
|
|
|
494 |
|
|
|
495 |
c ********************************************************************************
|
|
|
496 |
c * SUBROUTINE: ROTATION TO LATITUDE/LONGITUDE COORDINATE SYSTEM *
|
|
|
497 |
c ********************************************************************************
|
|
|
498 |
|
|
|
499 |
c --------------------------------------------------------------------------------
|
|
|
500 |
c Subroutine to get environment of strcof
|
|
|
501 |
c --------------------------------------------------------------------------------
|
|
|
502 |
|
|
|
503 |
SUBROUTINE getenvir (clon,clat,rotation,rotmode,
|
|
|
504 |
> inar, rdx,rdy,rxmin,rymin,rnx,rny,rnz,mdv,
|
|
|
505 |
> outar, dx, dy, xmin, ymin, nx, ny, nz)
|
|
|
506 |
|
|
|
507 |
c Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate
|
|
|
508 |
c system.
|
|
|
509 |
|
|
|
510 |
implicit none
|
|
|
511 |
|
|
|
512 |
c Declaration of input parameters
|
|
|
513 |
integer rnx,rny,rnz
|
|
|
514 |
real rdx,rdy,rxmin,rymin
|
|
|
515 |
real inar(rnx,rny,rnz)
|
|
|
516 |
real clon,clat,rotation
|
|
|
517 |
real mdv
|
|
|
518 |
integer rotmode
|
|
|
519 |
|
|
|
520 |
c Declaration of output parameters
|
|
|
521 |
integer nx,ny,nz
|
|
|
522 |
real dx,dy,xmin,ymin
|
|
|
523 |
real outar(nx,ny,nz)
|
|
|
524 |
|
|
|
525 |
c Set numerical and physical constants
|
|
|
526 |
real g2r
|
|
|
527 |
parameter (g2r=0.0174533)
|
|
|
528 |
real pi180
|
|
|
529 |
parameter (pi180=3.14159265359/180.)
|
|
|
530 |
real eps
|
|
|
531 |
parameter (eps=0.0001)
|
|
|
532 |
|
|
|
533 |
|
|
|
534 |
c Flag for test mode
|
|
|
535 |
integer test
|
|
|
536 |
parameter (test=0)
|
|
|
537 |
character*80 testfile
|
|
|
538 |
parameter (testfile='ROTGRID')
|
|
|
539 |
|
|
|
540 |
c Auxiliary variables
|
|
|
541 |
real pollon,pollat
|
|
|
542 |
integer i,j,k,l
|
|
|
543 |
real rlon(nx,ny),rlat(nx,ny)
|
|
|
544 |
real rlon1(nx,ny),rlat1(nx,ny)
|
|
|
545 |
real lon(nx,ny),lat(nx,ny)
|
|
|
546 |
real rotangle1(nx,ny),rotangle2(nx,ny)
|
|
|
547 |
real rotangle(nx,ny)
|
|
|
548 |
real sinoutar(nx,ny,nz),cosoutar(nx,ny,nz)
|
|
|
549 |
real rxmax,rymax
|
|
|
550 |
real xind,yind,pind
|
|
|
551 |
real outval
|
|
|
552 |
integer ix,iy
|
|
|
553 |
real ax,ay,az,bx,by,bz,zx,zy,zz
|
|
|
554 |
real clon1,clat1,clon2,clat2
|
|
|
555 |
real rindx,rindy
|
|
|
556 |
integer indx,indy,indr,indu
|
|
|
557 |
real frac0i,frac0j,frac1i,frac1j
|
|
|
558 |
character*80 cdfname,varname,cstname
|
|
|
559 |
real vx1,vx2,vy1,vy2,angle,vx2min
|
|
|
560 |
integer s
|
|
|
561 |
|
|
|
562 |
c Externals
|
|
|
563 |
real lmtolms,phtophs
|
|
|
564 |
external lmtolms,phtophs
|
|
|
565 |
|
|
|
566 |
c Get geographical coordinates
|
|
|
567 |
do i=1,nx
|
|
|
568 |
do j=1,ny
|
|
|
569 |
lon(i,j)=xmin+real(i-1)*dx
|
|
|
570 |
lat(i,j)=ymin+real(j-1)*dy
|
|
|
571 |
enddo
|
|
|
572 |
enddo
|
|
|
573 |
|
|
|
574 |
c First rotation
|
|
|
575 |
pollon=clon-180.
|
|
|
576 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
577 |
pollat=90.-clat
|
|
|
578 |
do i=1,nx
|
|
|
579 |
do j=1,ny
|
|
|
580 |
rlon1(i,j)=lmtolms(lat(i,j),lon(i,j),pollat,pollon)
|
|
|
581 |
rlat1(i,j)=phtophs(lat(i,j),lon(i,j),pollat,pollon)
|
|
|
582 |
enddo
|
|
|
583 |
enddo
|
|
|
584 |
|
|
|
585 |
c Second coordinate transformation
|
|
|
586 |
pollon=-180.
|
|
|
587 |
pollat=90.+rotation
|
|
|
588 |
do i=1,nx
|
|
|
589 |
do j=1,ny
|
|
|
590 |
rlon(i,j)=90.+lmtolms(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)
|
|
|
591 |
rlat(i,j)=phtophs(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)
|
|
|
592 |
enddo
|
|
|
593 |
enddo
|
|
|
594 |
|
|
|
595 |
c Get the angle between the rotated and non-rotated coordinate frame
|
|
|
596 |
if ((rotmode.eq.1).or.(rotmode.eq.2)) then
|
|
|
597 |
do i=2,nx-1
|
|
|
598 |
do j=2,ny-1
|
|
|
599 |
|
|
|
600 |
c Angle between latitude circles
|
|
|
601 |
vx1=1.
|
|
|
602 |
vy1=0.
|
|
|
603 |
vx2min=(rlon(i+1,j)-rlon(i-1,j))*cos(pi180*rlat(i,j))
|
|
|
604 |
do s=-1,1,2
|
|
|
605 |
vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)*
|
|
|
606 |
> cos(pi180*rlat(i,j))
|
|
|
607 |
if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
|
|
|
608 |
enddo
|
|
|
609 |
vx2=vx2min
|
|
|
610 |
vy2=rlat(i+1,j)-rlat(i-1,j)
|
|
|
611 |
call getangle(vx1,vy1,vx2,vy2,angle)
|
|
|
612 |
rotangle1(i,j)=angle
|
|
|
613 |
|
|
|
614 |
c Angle between longitude circles
|
|
|
615 |
vx1=0.
|
|
|
616 |
vy1=1.
|
|
|
617 |
vx2min=(rlon(i,j+1)-rlon(i,j-1))*cos(pi180*rlat(i,j))
|
|
|
618 |
do s=-1,1,2
|
|
|
619 |
vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)*
|
|
|
620 |
> cos(pi180*rlat(i,j))
|
|
|
621 |
if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
|
|
|
622 |
enddo
|
|
|
623 |
vx2=vx2min
|
|
|
624 |
vy2=rlat(i,j+1)-rlat(i,j-1)
|
|
|
625 |
call getangle(vx1,vy1,vx2,vy2,angle)
|
|
|
626 |
rotangle2(i,j)=angle
|
|
|
627 |
|
|
|
628 |
enddo
|
|
|
629 |
enddo
|
|
|
630 |
|
|
|
631 |
c Set the angle at the boundaries
|
|
|
632 |
do i=1,nx
|
|
|
633 |
rotangle1(i,ny)=2.0*rotangle1(i,ny-1)-rotangle1(i,ny-2)
|
|
|
634 |
rotangle1(i,1) =2.0*rotangle1(i,2)-rotangle1(i,3)
|
|
|
635 |
rotangle2(i,ny)=2.0*rotangle2(i,ny-1)-rotangle2(i,ny-2)
|
|
|
636 |
rotangle2(i,1) =2.0*rotangle2(i,2)-rotangle2(i,3)
|
|
|
637 |
enddo
|
|
|
638 |
do j=1,ny
|
|
|
639 |
rotangle1(nx,j)=2.0*rotangle1(nx-1,j)-rotangle1(nx-2,j)
|
|
|
640 |
rotangle1(1,j) =2.0*rotangle1(2,j)-rotangle1(3,j)
|
|
|
641 |
rotangle2(nx,j)=2.0*rotangle2(nx-1,j)-rotangle2(nx-2,j)
|
|
|
642 |
rotangle2(1,j) =2.0*rotangle2(2,j)-rotangle2(3,j)
|
|
|
643 |
enddo
|
|
|
644 |
|
|
|
645 |
c Set the final rotation angle
|
|
|
646 |
do i=1,nx
|
|
|
647 |
do j=1,ny
|
|
|
648 |
rotangle(i,j)=0.5*(rotangle1(i,j)+rotangle2(i,j))
|
|
|
649 |
enddo
|
|
|
650 |
enddo
|
|
|
651 |
|
|
|
652 |
endif
|
|
|
653 |
|
|
|
654 |
c Bring longitude into domain of geographical grid (shift by 360 deg)
|
|
|
655 |
do i=1,nx
|
|
|
656 |
do j=1,ny
|
|
|
657 |
100 if (rlon(i,j).lt.rxmin) then
|
|
|
658 |
rlon(i,j)=rlon(i,j)+360.
|
|
|
659 |
goto 100
|
|
|
660 |
endif
|
|
|
661 |
102 if (rlon(i,j).gt.(rxmin+real(rnx-1)*rdx)) then
|
|
|
662 |
rlon(i,j)=rlon(i,j)-360.
|
|
|
663 |
goto 102
|
|
|
664 |
endif
|
|
|
665 |
enddo
|
|
|
666 |
enddo
|
|
|
667 |
|
|
|
668 |
c Rotate the scalar or the vector component
|
|
|
669 |
do i=1,nx
|
|
|
670 |
do j=1,ny
|
|
|
671 |
do k=1,nz
|
|
|
672 |
|
|
|
673 |
c Get index
|
|
|
674 |
rindx=(rlon(i,j)-rxmin)/rdx+1.
|
|
|
675 |
rindy=(rlat(i,j)-rymin)/rdy+1.
|
|
|
676 |
indx=int(rindx)
|
|
|
677 |
indy=int(rindy)
|
|
|
678 |
if ((indx.lt.1).or.(indx.gt.rnx).or.
|
|
|
679 |
> (indy.lt.1).or.(indy.gt.rny)) then
|
|
|
680 |
outar(i,j,k)=mdv
|
|
|
681 |
goto 103
|
|
|
682 |
endif
|
|
|
683 |
|
|
|
684 |
c Get inidices of left and upper neighbours
|
|
|
685 |
indr=indx+1
|
|
|
686 |
if (indr.gt.rnx) indr=1
|
|
|
687 |
indu=indy+1
|
|
|
688 |
if (indu.gt.rny) indu=ny
|
|
|
689 |
|
|
|
690 |
c Do linear interpolation
|
|
|
691 |
if ( ( abs(inar(indx ,indy, k)-mdv).gt.eps ).and.
|
|
|
692 |
& ( abs(inar(indx ,indu, k)-mdv).gt.eps ).and.
|
|
|
693 |
& ( abs(inar(indr ,indy, k)-mdv).gt.eps ).and.
|
|
|
694 |
& ( abs(inar(indr ,indu, k)-mdv).gt.eps ) ) then
|
|
|
695 |
frac0i=rindx-float(indx)
|
|
|
696 |
frac0j=rindy-float(indy)
|
|
|
697 |
frac1i=1.-frac0i
|
|
|
698 |
frac1j=1.-frac0j
|
|
|
699 |
outval = inar(indx ,indy, k ) * frac1i * frac1j
|
|
|
700 |
& + inar(indx ,indu, k ) * frac1i * frac0j
|
|
|
701 |
& + inar(indr ,indy, k ) * frac0i * frac1j
|
|
|
702 |
& + inar(indr ,indu, k ) * frac0i * frac0j
|
|
|
703 |
else
|
|
|
704 |
outval=mdv
|
|
|
705 |
endif
|
|
|
706 |
|
|
|
707 |
c Update output array
|
|
|
708 |
outar(i,j,k)=outval
|
|
|
709 |
|
|
|
710 |
c Next
|
|
|
711 |
103 continue
|
|
|
712 |
|
|
|
713 |
enddo
|
|
|
714 |
enddo
|
|
|
715 |
enddo
|
|
|
716 |
|
|
|
717 |
c Get components for tests
|
|
|
718 |
if (test.eq.1) then
|
|
|
719 |
do i=1,nx
|
|
|
720 |
do j=1,ny
|
|
|
721 |
do k=1,nz
|
|
|
722 |
cosoutar(i,j,k)=outar(i,j,k)*cos(pi180*rotangle(i,j))
|
|
|
723 |
sinoutar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j))
|
|
|
724 |
enddo
|
|
|
725 |
enddo
|
|
|
726 |
enddo
|
|
|
727 |
endif
|
|
|
728 |
|
|
|
729 |
c Get correct component of rotated field
|
|
|
730 |
do i=1,nx
|
|
|
731 |
do j=1,ny
|
|
|
732 |
do k=1,nz
|
|
|
733 |
if ( abs(outar(i,j,k)-mdv).gt.eps ) then
|
|
|
734 |
if (rotmode.eq.1) then
|
|
|
735 |
outar(i,j,k)= outar(i,j,k)*cos(pi180*rotangle(i,j))
|
|
|
736 |
else if (rotmode.eq.2) then
|
|
|
737 |
outar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j))
|
|
|
738 |
else if (rotmode.eq.0) then
|
|
|
739 |
outar(i,j,k)=outar(i,j,k)
|
|
|
740 |
endif
|
|
|
741 |
endif
|
|
|
742 |
enddo
|
|
|
743 |
enddo
|
|
|
744 |
enddo
|
|
|
745 |
|
|
|
746 |
c Test mode: Write grids to cdf file
|
|
|
747 |
if (test.eq.1) then
|
|
|
748 |
cdfname=testfile
|
|
|
749 |
cstname=trim(testfile)//'_cst'
|
|
|
750 |
varname='RLON1'
|
|
|
751 |
call writecdf2D(cdfname,cstname,
|
|
|
752 |
> varname,rlon1,0.,
|
|
|
753 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
754 |
> 1,1)
|
|
|
755 |
cdfname=testfile
|
|
|
756 |
cstname=trim(testfile)//'_cst'
|
|
|
757 |
varname='RLON'
|
|
|
758 |
call writecdf2D(cdfname,cstname,
|
|
|
759 |
> varname,rlon,0.,
|
|
|
760 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
761 |
> 0,1)
|
|
|
762 |
cdfname=testfile
|
|
|
763 |
cstname=trim(testfile)//'_cst'
|
|
|
764 |
varname='LON'
|
|
|
765 |
call writecdf2D(cdfname,cstname,
|
|
|
766 |
> varname,lon,0.,
|
|
|
767 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
768 |
> 0,1)
|
|
|
769 |
cdfname=testfile
|
|
|
770 |
cstname=trim(testfile)//'_cst'
|
|
|
771 |
varname='RLAT1'
|
|
|
772 |
call writecdf2D(cdfname,cstname,
|
|
|
773 |
> varname,rlat1,0.,
|
|
|
774 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
775 |
> 0,1)
|
|
|
776 |
cdfname=testfile
|
|
|
777 |
cstname=trim(testfile)//'_cst'
|
|
|
778 |
varname='RLAT'
|
|
|
779 |
call writecdf2D(cdfname,cstname,
|
|
|
780 |
> varname,rlat,0.,
|
|
|
781 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
782 |
> 0,1)
|
|
|
783 |
cdfname=testfile
|
|
|
784 |
cstname=trim(testfile)//'_cst'
|
|
|
785 |
varname='LAT'
|
|
|
786 |
call writecdf2D(cdfname,cstname,
|
|
|
787 |
> varname,lat,0.,
|
|
|
788 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
789 |
> 0,1)
|
|
|
790 |
cdfname=testfile
|
|
|
791 |
cstname=trim(testfile)//'_cst'
|
|
|
792 |
varname='ANGLE1'
|
|
|
793 |
call writecdf2D(cdfname,cstname,
|
|
|
794 |
> varname,rotangle1,0.,
|
|
|
795 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
796 |
> 0,1)
|
|
|
797 |
cdfname=testfile
|
|
|
798 |
cstname=trim(testfile)//'_cst'
|
|
|
799 |
varname='ANGLE2'
|
|
|
800 |
call writecdf2D(cdfname,cstname,
|
|
|
801 |
> varname,rotangle2,0.,
|
|
|
802 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
803 |
> 0,1)
|
|
|
804 |
cdfname=testfile
|
|
|
805 |
cstname=trim(testfile)//'_cst'
|
|
|
806 |
varname='U'
|
|
|
807 |
call writecdf2D(cdfname,cstname,
|
|
|
808 |
> varname,cosoutar,0.,
|
|
|
809 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
810 |
> 0,1)
|
|
|
811 |
cdfname=testfile
|
|
|
812 |
cstname=trim(testfile)//'_cst'
|
|
|
813 |
varname='V'
|
|
|
814 |
call writecdf2D(cdfname,cstname,
|
|
|
815 |
> varname,sinoutar,0.,
|
|
|
816 |
> rdx,rdy,rxmin,rymin,rnx,rny,
|
|
|
817 |
> 0,1)
|
|
|
818 |
|
|
|
819 |
endif
|
|
|
820 |
|
|
|
821 |
END
|
|
|
822 |
|
|
|
823 |
|
|
|
824 |
c --------------------------------------------------------------------------------
|
|
|
825 |
c Auxiliary routines: angle between two vectors
|
|
|
826 |
c --------------------------------------------------------------------------------
|
|
|
827 |
|
|
|
828 |
SUBROUTINE getangle (vx1,vy1,vx2,vy2,angle)
|
|
|
829 |
|
|
|
830 |
c Given two vectors <vx1,vy1> and <vx2,vy2>, determine the angle (in deg)
|
|
|
831 |
c between the two vectors.
|
|
|
832 |
|
|
|
833 |
implicit none
|
|
|
834 |
|
|
|
835 |
c Declaration of subroutine parameters
|
|
|
836 |
real vx1,vy1
|
|
|
837 |
real vx2,vy2
|
|
|
838 |
real angle
|
|
|
839 |
|
|
|
840 |
c Auxiliary variables and parameters
|
|
|
841 |
real len1,len2,len3
|
|
|
842 |
real val1,val2,val3
|
|
|
843 |
real pi
|
|
|
844 |
parameter (pi=3.14159265359)
|
|
|
845 |
|
|
|
846 |
len1=sqrt(vx1*vx1+vy1*vy1)
|
|
|
847 |
len2=sqrt(vx2*vx2+vy2*vy2)
|
|
|
848 |
|
|
|
849 |
if ((len1.gt.0.).and.(len2.gt.0.)) then
|
|
|
850 |
vx1=vx1/len1
|
|
|
851 |
vy1=vy1/len1
|
|
|
852 |
vx2=vx2/len2
|
|
|
853 |
vy2=vy2/len2
|
|
|
854 |
|
|
|
855 |
val1=vx1*vx2+vy1*vy2
|
|
|
856 |
val2=-vy1*vx2+vx1*vy2
|
|
|
857 |
|
|
|
858 |
len3=sqrt(val1*val1+val2*val2)
|
|
|
859 |
|
|
|
860 |
if ( (val1.ge.0.).and.(val2.ge.0.) ) then
|
|
|
861 |
val3=acos(val1/len3)
|
|
|
862 |
else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
|
|
|
863 |
val3=pi-acos(abs(val1)/len3)
|
|
|
864 |
else if ( (val1.ge.0.).and.(val2.le.0.) ) then
|
|
|
865 |
val3=-acos(val1/len3)
|
|
|
866 |
else if ( (val1.lt.0.).and.(val2.le.0.) ) then
|
|
|
867 |
val3=-pi+acos(abs(val1)/len3)
|
|
|
868 |
endif
|
|
|
869 |
else
|
|
|
870 |
val3=0.
|
|
|
871 |
endif
|
|
|
872 |
|
|
|
873 |
angle=180./pi*val3
|
|
|
874 |
|
|
|
875 |
END
|
|
|
876 |
|
|
|
877 |
c --------------------------------------------------------------------------------
|
|
|
878 |
c Transformation routine: LMSTOLM and PHSTOPH from library gm2em
|
|
|
879 |
c --------------------------------------------------------------------------------
|
|
|
880 |
|
|
|
881 |
REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
882 |
C
|
|
|
883 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
884 |
C
|
|
|
885 |
C**** LMTOLMS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
|
|
|
886 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
887 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
888 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
889 |
C** AUFRUF : LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
|
|
|
890 |
C** ENTRIES : KEINE
|
|
|
891 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
|
|
|
892 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
893 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
894 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
895 |
C** VERSIONS-
|
|
|
896 |
C** DATUM : 03.05.90
|
|
|
897 |
C**
|
|
|
898 |
C** EXTERNALS: KEINE
|
|
|
899 |
C** EINGABE-
|
|
|
900 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
901 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
902 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
903 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
904 |
C** AUSGABE-
|
|
|
905 |
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
|
|
|
906 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
907 |
C**
|
|
|
908 |
C** COMMON-
|
|
|
909 |
C** BLOECKE : KEINE
|
|
|
910 |
C**
|
|
|
911 |
C** FEHLERBE-
|
|
|
912 |
C** HANDLUNG : KEINE
|
|
|
913 |
C** VERFASSER: G. DE MORSIER
|
|
|
914 |
|
|
|
915 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
916 |
|
|
|
917 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
918 |
|
|
|
919 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
920 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
921 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
922 |
ZPHI = ZPIR18*PHI
|
|
|
923 |
ZLAM = LAM
|
|
|
924 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
925 |
ZLAM = ZPIR18*ZLAM
|
|
|
926 |
|
|
|
927 |
ZARG1 = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
|
|
|
928 |
ZARG2 = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
|
|
|
929 |
IF (ABS(ZARG2).LT.1.E-30) THEN
|
|
|
930 |
IF (ABS(ZARG1).LT.1.E-30) THEN
|
|
|
931 |
LMTOLMS = 0.0
|
|
|
932 |
ELSEIF (ZARG1.GT.0.) THEN
|
|
|
933 |
LMTOLMS = 90.0
|
|
|
934 |
ELSE
|
|
|
935 |
LMTOLMS = -90.0
|
|
|
936 |
ENDIF
|
|
|
937 |
ELSE
|
|
|
938 |
LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
|
|
|
939 |
ENDIF
|
|
|
940 |
|
|
|
941 |
RETURN
|
|
|
942 |
END
|
|
|
943 |
|
|
|
944 |
|
|
|
945 |
REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
946 |
C
|
|
|
947 |
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
|
|
|
948 |
C
|
|
|
949 |
C**** PHTOPHS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
|
|
|
950 |
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
|
|
|
951 |
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
|
|
|
952 |
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
953 |
C** AUFRUF : PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
|
|
|
954 |
C** ENTRIES : KEINE
|
|
|
955 |
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
|
|
|
956 |
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
|
|
|
957 |
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
|
|
|
958 |
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
|
|
|
959 |
C** VERSIONS-
|
|
|
960 |
C** DATUM : 03.05.90
|
|
|
961 |
C**
|
|
|
962 |
C** EXTERNALS: KEINE
|
|
|
963 |
C** EINGABE-
|
|
|
964 |
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
965 |
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
|
|
|
966 |
C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
|
|
|
967 |
C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
|
|
|
968 |
C** AUSGABE-
|
|
|
969 |
C** PARAMETER: ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
|
|
|
970 |
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
|
|
|
971 |
C**
|
|
|
972 |
C** COMMON-
|
|
|
973 |
C** BLOECKE : KEINE
|
|
|
974 |
C**
|
|
|
975 |
C** FEHLERBE-
|
|
|
976 |
C** HANDLUNG : KEINE
|
|
|
977 |
C** VERFASSER: G. DE MORSIER
|
|
|
978 |
|
|
|
979 |
REAL LAM,PHI,POLPHI,POLLAM
|
|
|
980 |
|
|
|
981 |
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /
|
|
|
982 |
|
|
|
983 |
ZSINPOL = SIN(ZPIR18*POLPHI)
|
|
|
984 |
ZCOSPOL = COS(ZPIR18*POLPHI)
|
|
|
985 |
ZLAMPOL = ZPIR18*POLLAM
|
|
|
986 |
ZPHI = ZPIR18*PHI
|
|
|
987 |
ZLAM = LAM
|
|
|
988 |
IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
|
|
|
989 |
ZLAM = ZPIR18*ZLAM
|
|
|
990 |
ZARG = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
|
|
|
991 |
|
|
|
992 |
PHTOPHS = ZRPI18*ASIN(ZARG)
|
|
|
993 |
|
|
|
994 |
RETURN
|
|
|
995 |
END
|
|
|
996 |
|
|
|
997 |
|
|
|
998 |
c ---------------------------------------------------------
|
|
|
999 |
c Spherical distance between two lat/lon points
|
|
|
1000 |
c ---------------------------------------------------------
|
|
|
1001 |
|
|
|
1002 |
real function sdis(xp,yp,xq,yq)
|
|
|
1003 |
c
|
|
|
1004 |
c calculates spherical distance (in km) between two points given
|
|
|
1005 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
1006 |
c
|
|
|
1007 |
real re
|
|
|
1008 |
parameter (re=6370.)
|
|
|
1009 |
real xp,yp,xq,yq,arg
|
|
|
1010 |
|
|
|
1011 |
arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
|
|
|
1012 |
if (arg.lt.-1.) arg=-1.
|
|
|
1013 |
if (arg.gt.1.) arg=1.
|
|
|
1014 |
sdis=re*acos(arg)
|
|
|
1015 |
|
|
|
1016 |
end
|
|
|
1017 |
|
|
|
1018 |
|
|
|
1019 |
c ********************************************************************************
|
|
|
1020 |
c * NETCDF INPUT/OUTPUT *
|
|
|
1021 |
c ********************************************************************************
|
|
|
1022 |
|
|
|
1023 |
c ---------------------------------------------------------
|
|
|
1024 |
c Subroutines to write the netcdf output file
|
|
|
1025 |
c ---------------------------------------------------------
|
|
|
1026 |
|
|
|
1027 |
subroutine writecdf2D(cdfname,cstname,
|
|
|
1028 |
> varname,arr,ctime,
|
|
|
1029 |
> dx,dy,xmin,ymin,nx,ny,
|
|
|
1030 |
> crefile,crevar)
|
|
|
1031 |
|
|
|
1032 |
c Create and write to the netcdf file <cdfname>. The variable
|
|
|
1033 |
c with name <varname> and with time <time> is written. The data
|
|
|
1034 |
c are in the two-dimensional array <arr>. The list <dx,dy,xmin,
|
|
|
1035 |
c ymin,nx,ny> specifies the output grid. The flags <crefile> and
|
|
|
1036 |
c <crevar> determine whether the file and/or the variable should
|
|
|
1037 |
c be created. 1: create / 0: not created
|
|
|
1038 |
|
|
|
1039 |
IMPLICIT NONE
|
|
|
1040 |
|
|
|
1041 |
c Declaration of input parameters
|
|
|
1042 |
character*80 cdfname,cstname,varname
|
|
|
1043 |
integer nx,ny
|
|
|
1044 |
real arr(nx,ny)
|
|
|
1045 |
real dx,dy,xmin,ymin
|
|
|
1046 |
integer ctime
|
|
|
1047 |
integer crefile,crevar
|
|
|
1048 |
|
|
|
1049 |
c Further variables
|
|
|
1050 |
real varmin(4),varmax(4),stag(4)
|
|
|
1051 |
integer ierr,cdfid,ndim,vardim(4)
|
|
|
1052 |
real mdv
|
|
|
1053 |
integer datar(14),stdate(5)
|
|
|
1054 |
integer i
|
|
|
1055 |
real time
|
|
|
1056 |
|
|
|
1057 |
C Definitions
|
|
|
1058 |
varmin(1)=xmin
|
|
|
1059 |
varmin(2)=ymin
|
|
|
1060 |
varmin(3)=1050.
|
|
|
1061 |
varmax(1)=xmin+real(nx-1)*dx
|
|
|
1062 |
varmax(2)=ymin+real(ny-1)*dy
|
|
|
1063 |
varmax(3)=1050.
|
|
|
1064 |
ndim=4
|
|
|
1065 |
vardim(1)=nx
|
|
|
1066 |
vardim(2)=ny
|
|
|
1067 |
vardim(3)=1
|
|
|
1068 |
stag(1)=0.
|
|
|
1069 |
stag(2)=0.
|
|
|
1070 |
stag(3)=0.
|
|
|
1071 |
mdv=-999.98999
|
|
|
1072 |
time=real(ctime)
|
|
|
1073 |
|
|
|
1074 |
C Create the file
|
|
|
1075 |
if (crefile.eq.0) then
|
|
|
1076 |
call cdfwopn(cdfname,cdfid,ierr)
|
|
|
1077 |
if (ierr.ne.0) goto 906
|
|
|
1078 |
else if (crefile.eq.1) then
|
|
|
1079 |
call crecdf(cdfname,cdfid,
|
|
|
1080 |
> varmin,varmax,ndim,cstname,ierr)
|
|
|
1081 |
if (ierr.ne.0) goto 903
|
|
|
1082 |
|
|
|
1083 |
C Write the constants file
|
|
|
1084 |
datar(1)=vardim(1)
|
|
|
1085 |
datar(2)=vardim(2)
|
|
|
1086 |
datar(3)=int(1000.*varmax(2))
|
|
|
1087 |
datar(4)=int(1000.*varmin(1))
|
|
|
1088 |
datar(5)=int(1000.*varmin(2))
|
|
|
1089 |
datar(6)=int(1000.*varmax(1))
|
|
|
1090 |
datar(7)=int(1000.*dx)
|
|
|
1091 |
datar(8)=int(1000.*dy)
|
|
|
1092 |
datar(9)=1
|
|
|
1093 |
datar(10)=1
|
|
|
1094 |
datar(11)=0 ! data version
|
|
|
1095 |
datar(12)=0 ! cstfile version
|
|
|
1096 |
datar(13)=0 ! longitude of pole
|
|
|
1097 |
datar(14)=90000 ! latitude of pole
|
|
|
1098 |
do i=1,5
|
|
|
1099 |
stdate(i)=0
|
|
|
1100 |
enddo
|
|
|
1101 |
call wricst(cstname,datar,0.,0.,0.,0.,stdate)
|
|
|
1102 |
endif
|
|
|
1103 |
|
|
|
1104 |
c Write the data to the netcdf file, and close the file
|
|
|
1105 |
if (crevar.eq.1) then
|
|
|
1106 |
call putdef(cdfid,varname,ndim,mdv,
|
|
|
1107 |
> vardim,varmin,varmax,stag,ierr)
|
|
|
1108 |
if (ierr.ne.0) goto 904
|
|
|
1109 |
endif
|
|
|
1110 |
call putdat(cdfid,varname,time,0,arr,ierr)
|
|
|
1111 |
if (ierr.ne.0) goto 905
|
|
|
1112 |
call clscdf(cdfid,ierr)
|
|
|
1113 |
|
|
|
1114 |
return
|
|
|
1115 |
|
|
|
1116 |
c Error handling
|
|
|
1117 |
903 print*,'*** Problem to create netcdf file ***'
|
|
|
1118 |
stop
|
|
|
1119 |
904 print*,'*** Problem to write definition ***'
|
|
|
1120 |
stop
|
|
|
1121 |
905 print*,'*** Problem to write data ***'
|
|
|
1122 |
stop
|
|
|
1123 |
906 print*,'*** Problem to open netcdf file ***'
|
|
|
1124 |
stop
|
|
|
1125 |
|
|
|
1126 |
END
|
|
|
1127 |
|
|
|
1128 |
c ----------------------------------------------------------------
|
|
|
1129 |
c Check whether variable is found on netcdf file
|
|
|
1130 |
c ----------------------------------------------------------------
|
|
|
1131 |
|
|
|
1132 |
subroutine check_varok (isok,varname,varlist,nvars)
|
|
|
1133 |
|
|
|
1134 |
c Check whether the variable <varname> is in the list <varlist(nvars)>.
|
|
|
1135 |
c If this is the case, <isok> is incremented by 1. Otherwise <isok>
|
|
|
1136 |
c keeps its value.
|
|
|
1137 |
|
|
|
1138 |
implicit none
|
|
|
1139 |
|
|
|
1140 |
c Declaraion of subroutine parameters
|
|
|
1141 |
integer isok
|
|
|
1142 |
integer nvars
|
|
|
1143 |
character*80 varname
|
|
|
1144 |
character*80 varlist(nvars)
|
|
|
1145 |
|
|
|
1146 |
c Auxiliary variables
|
|
|
1147 |
integer i
|
|
|
1148 |
|
|
|
1149 |
c Main
|
|
|
1150 |
do i=1,nvars
|
|
|
1151 |
if (trim(varname).eq.trim(varlist(i))) isok=isok+1
|
|
|
1152 |
enddo
|
|
|
1153 |
|
|
|
1154 |
end
|