Subversion Repositories pvinversion.ecmwf

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
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