Subversion Repositories lagranto.wrf

Rev

Rev 20 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 20 Rev 21
Line 23... Line 23...
23 
      parameter    (mapfile ='wrfmap.nc')
 23 
      parameter    (mapfile ='wrfmap.nc')
24 
      real         mdv                              ! missing data value
 24 
      real         mdv                              ! missing data value
25 
      parameter    (mdv=-999.)
 25 
      parameter    (mdv=-999.)
26 
      real         eps                              ! Numerical epsilon
 26 
      real         eps                              ! Numerical epsilon
27 
      parameter    (eps= 0.001)
 27 
      parameter    (eps= 0.001)
- 
 
 28 
      real         rad2deg                          ! from radian to deg
- 
 
 29 
      parameter    (rad2deg=180./3.141592653)
28 
 
 30 
 
29 
c     Numerical grid
 31 
c     Numerical grid
30 
      real         xmin,xmax,ymin,ymax              ! Domain size
 32 
      real         xmin,xmax,ymin,ymax              ! Domain size
31 
      real         dx,dy                            ! Horizontal resolution
 33 
      real         dx,dy                            ! Horizontal resolution
32 
      integer      nx,ny,nz                         ! Grid dimensions
 34 
      integer      nx,ny,nz                         ! Grid dimensions
33 
      real,allocatable,dimension (:,:)  :: lon,lat  ! lon/lat on numerical grid
 35 
      real,allocatable,dimension (:,:)  :: lon,lat  ! lon/lat on numerical grid
34 
      real,allocatable,dimension (:,:)  :: mpx,mpy  ! map scale factors in x/y direction
 36 
      real,allocatable,dimension (:,:)  :: mpx,mpy  ! map scale factors in x/y direction
35 
 
 - 
 
- 
 
 37 
      real,allocatable,dimension (:,:)  :: cosalpha,sinalpha  ! cos and sin of grid rotation
- 
 
 38 
      real,allocatable,dimension (:,:)  :: alpha    ! grid rotation
36 
 
 39 
 
37 
c     Geographical grid
 40 
c     Geographical grid
38 
      real         latmin,latmax,lonmin,lonmax      ! Geographical domain
 41 
      real         latmin,latmax,lonmin,lonmax      ! Geographical domain
39 
      real         dlon,dlat                        ! Minimum spacing in geographical space
 42 
      real         dlon,dlat                        ! Minimum spacing in geographical space
40 
      integer      nlon,nlat,nlev                        ! Geographical grid dimension
 43 
      integer      nlon,nlat,nlev                        ! Geographical grid dimension
Line 73... Line 76...
73 
      real         rdummy
 76 
      real         rdummy
74 
      real         rid,rjd,rkd
 77 
      real         rid,rjd,rkd
75 
      real         xpos,ypos,zpos,ppos,lonpos,latpos
 78 
      real         xpos,ypos,zpos,ppos,lonpos,latpos
76 
      real         lon1,lat1,lon2,lat2
 79 
      real         lon1,lat1,lon2,lat2
77 
      integer      nfilter
 80 
      integer      nfilter
- 
 
 81 
      real         uu0,vv0,ca0,sa0
- 
 
 82 
      character*80 vectorx,vectory
- 
 
 83 
      integer      indu,indv,len
- 
 
 84 
      character    char0,char1,char2
- 
 
 85 
      character*80 string
78 
 
 86 
 
79 
c     Externals
 87 
c     Externals
80 
      real         int_index3
 88 
      real         int_index3
81 
      external     int_index3
 89 
      external     int_index3
82 
      real         sdis
 90 
      real         sdis
Line 109... Line 117...
109 
       if ( mode.eq.'-xy2ll'  ) then         ! convert x/y -> lon/lat
 117 
       if ( mode.eq.'-xy2ll'  ) then         ! convert x/y -> lon/lat
110 
          read(10,*) inpfile
 118 
          read(10,*) inpfile
111 
          read(10,*) outfile
 119 
          read(10,*) outfile
112 
          read(10,*) ntra,ntim,ncol
 120 
          read(10,*) ntra,ntim,ncol
113 
       endif
 121 
       endif
- 
 
 122
- 
 
 123 
       if ( mode.eq.'-vec2ll'  ) then        ! convert vector to lon/lat components
- 
 
 124 
          read(10,*) inpfile
- 
 
 125 
          read(10,*) outfile
- 
 
 126 
          read(10,*) ntra,ntim,ncol
- 
 
 127 
          read(10,*) vectorx
- 
 
 128 
          read(10,*) vectory
- 
 
 129 
       endif
114 
 
 130
- 
 
 131 
       if ( mode.eq.'-vec2xy'  ) then        ! convert vector to x/y components
- 
 
 132 
          read(10,*) inpfile
- 
 
 133 
          read(10,*) outfile
- 
 
 134 
          read(10,*) ntra,ntim,ncol
- 
 
 135 
          read(10,*) vectorx
- 
 
 136 
          read(10,*) vectory
- 
 
 137 
       endif
- 
 
 138
115 
       if ( mode.eq.'-xy2ll.single'  ) then   ! convert single x/y -> lon/lat
 139 
       if ( mode.eq.'-xy2ll.single'  ) then   ! convert single x/y -> lon/lat
116 
          read(10,*) xpos,ypos
 140 
          read(10,*) xpos,ypos
117 
       endif
 141 
       endif
118 
 
 142 
 
119 
       if ( mode.eq.'-p2z'  ) then            ! convert x/y/p -> x/y/z
 143 
       if ( mode.eq.'-p2z'  ) then            ! convert x/y/p -> x/y/z
Line 162... Line 186...
162 
         call readcdf2D(cdfname,varname,rdummy,0.,
 186 
         call readcdf2D(cdfname,varname,rdummy,0.,
163 
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,nlev)
 187 
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,nlev)
164 
 
 188 
 
165 
c        Allocate memory
 189 
c        Allocate memory
166 
         allocate(lon(nx,ny),stat=stat)
 190 
         allocate(lon(nx,ny),stat=stat)
167 
         if (stat.ne.0) print*,'*** error allocating array lon ***'
 191 
         if (stat.ne.0) print*,'*** error allocating array lon      ***'
168 
         allocate(lat(nx,ny),stat=stat)
 192 
         allocate(lat(nx,ny),stat=stat)
169 
         if (stat.ne.0) print*,'*** error allocating array lat ***'
 193 
         if (stat.ne.0) print*,'*** error allocating array lat      ***'
170 
         allocate(x(nlon,nlat),stat=stat)
 194 
         allocate(x(nlon,nlat),stat=stat)
171 
         if (stat.ne.0) print*,'*** error allocating array x   ***'
 195 
         if (stat.ne.0) print*,'*** error allocating array x        ***'
172 
         allocate(y(nlon,nlat),stat=stat)
 196 
         allocate(y(nlon,nlat),stat=stat)
173 
         if (stat.ne.0) print*,'*** error allocating array y   ***'
 197 
         if (stat.ne.0) print*,'*** error allocating array y        ***'
174 
         allocate(p3(nx,ny,nz),stat=stat)
 198 
         allocate(p3(nx,ny,nz),stat=stat)
175 
         if (stat.ne.0) print*,'*** error allocating array p3  ***'
 199 
         if (stat.ne.0) print*,'*** error allocating array p3       ***'
176 
         allocate(z3(nx,ny,nz),stat=stat)
 200 
         allocate(z3(nx,ny,nz),stat=stat)
177 
         if (stat.ne.0) print*,'*** error allocating array z3  ***'
 201 
         if (stat.ne.0) print*,'*** error allocating array z3       ***'
178 
         allocate(zb(nx,ny),stat=stat)
 202 
         allocate(zb(nx,ny),stat=stat)
179 
         if (stat.ne.0) print*,'*** error allocating array zb  ***'
 203 
         if (stat.ne.0) print*,'*** error allocating array zb       ***'
180 
         allocate(ps(nx,ny),stat=stat)
 204 
         allocate(ps(nx,ny),stat=stat)
181 
         if (stat.ne.0) print*,'*** error allocating array ps  ***'
 205 
         if (stat.ne.0) print*,'*** error allocating array ps       ***'
- 
 
 206 
         allocate(cosalpha(nx,ny),stat=stat)
- 
 
 207 
         if (stat.ne.0) print*,'*** error allocating array cosalpha ***'
- 
 
 208 
         allocate(sinalpha(nx,ny),stat=stat)
- 
 
 209 
         if (stat.ne.0) print*,'*** error allocating array cosalpha ***'
- 
 
 210 
         allocate(alpha(nx,ny),stat=stat)
- 
 
 211 
         if (stat.ne.0) print*,'*** error allocating array alpha    ***'
182 
 
 212 
 
183 
c        Read data
 213 
c        Read data
184 
         cdfname = mapfile
 214 
         cdfname = mapfile
185 
         varname = 'LON'
 215 
         varname = 'LON'
186 
         call readcdf2D(cdfname,varname,lon,0.,
 216 
         call readcdf2D(cdfname,varname,lon,0.,
Line 207... Line 237...
207 
     >                  dx,dy,xmin,ymin,nx,ny,nz)
 237 
     >                  dx,dy,xmin,ymin,nx,ny,nz)
208 
         cdfname = mapfile
 238 
         cdfname = mapfile
209 
         varname = 'TOPO'
 239 
         varname = 'TOPO'
210 
         call readcdf2D(cdfname,varname,zb,0.,
 240 
         call readcdf2D(cdfname,varname,zb,0.,
211 
     >                  dx,dy,xmin,ymin,nx,ny,1)
 241 
     >                  dx,dy,xmin,ymin,nx,ny,1)
- 
 
 242 
         cdfname = mapfile
- 
 
 243 
         varname = 'COSALPHA'
- 
 
 244 
         call readcdf2D(cdfname,varname,cosalpha,0.,
- 
 
 245 
     >                  dx,dy,xmin,ymin,nx,ny,1)
- 
 
 246 
         cdfname = mapfile
- 
 
 247 
         varname = 'SINALPHA'
- 
 
 248 
         call readcdf2D(cdfname,varname,sinalpha,0.,
- 
 
 249 
     >                  dx,dy,xmin,ymin,nx,ny,1)
212 
 
 250
- 
 
 251 
         cdfname = mapfile
- 
 
 252 
         varname = 'ALPHA'
- 
 
 253 
         call readcdf2D(cdfname,varname,alpha,0.,
- 
 
 254 
     >                  dx,dy,xmin,ymin,nx,ny,1)
- 
 
 255
213 
c	     Set surface prtessure to lowest 3d pressure level
 256 
c	     Set surface prtessure to lowest 3d pressure level
214 
         do i=1,nx
 257 
         do i=1,nx
215 
            do j=1,ny
 258 
            do j=1,ny
216 
               ps(i,j) = p3(i,j,1)
 259 
               ps(i,j) = p3(i,j,1)
217 
            enddo
 260 
            enddo
Line 245... Line 288...
245 
      dx   = 1.
 288 
      dx   = 1.
246 
      dy   = 1.
 289 
      dy   = 1.
247 
 
 290 
 
248 
c     Allocate memory for lon/lat on numerical grid
 291 
c     Allocate memory for lon/lat on numerical grid
249 
      allocate(lon(nx,ny),stat=stat)
 292 
      allocate(lon(nx,ny),stat=stat)
250 
      if (stat.ne.0) print*,'*** error allocating array lon ***'
 293 
      if (stat.ne.0) print*,'*** error allocating array lon       ***'
251 
      allocate(lat(nx,ny),stat=stat)
 294 
      allocate(lat(nx,ny),stat=stat)
252 
      if (stat.ne.0) print*,'*** error allocating array lat ***'
 295 
      if (stat.ne.0) print*,'*** error allocating array lat       ***'
253 
      allocate(p3(nx,ny,nz),stat=stat)
 296 
      allocate(p3(nx,ny,nz),stat=stat)
254 
      if (stat.ne.0) print*,'*** error allocating array p3  ***'
 297 
      if (stat.ne.0) print*,'*** error allocating array p3        ***'
255 
      allocate(z3(nx,ny,nz),stat=stat)
 298 
      allocate(z3(nx,ny,nz),stat=stat)
256 
      if (stat.ne.0) print*,'*** error allocating array z3  ***'
 299 
      if (stat.ne.0) print*,'*** error allocating array z3        ***'
257 
      allocate(zb(nx,ny),stat=stat)
 300 
      allocate(zb(nx,ny),stat=stat)
258 
      if (stat.ne.0) print*,'*** error allocating array zb  ***'
 301 
      if (stat.ne.0) print*,'*** error allocating array zb        ***'
- 
 
 302 
      allocate(cosalpha(nx,ny),stat=stat)
- 
 
 303 
      if (stat.ne.0) print*,'*** error allocating array cosalpha  ***'
- 
 
 304 
      allocate(sinalpha(nx,ny),stat=stat)
- 
 
 305 
      if (stat.ne.0) print*,'*** error allocating array sinalpha  ***'
- 
 
 306 
      allocate(alpha(nx,ny),stat=stat)
- 
 
 307 
      if (stat.ne.0) print*,'*** error allocating array alpha     ***'
259 
 
 308 
 
260 
c     Read lon/lat on the model grid
 309 
c     Read lon/lat on the model grid
261 
      varname='XLONG'
 310 
      varname='XLONG'
262 
      call input_var_wrf(fid,varname,lon,nx,ny,1)
 311 
      call input_var_wrf(fid,varname,lon      ,nx,ny,1)
263 
      varname='XLAT'
 312 
      varname='XLAT'
264 
      call input_var_wrf(fid,varname,lat,nx,ny,1)
 313 
      call input_var_wrf(fid,varname,lat      ,nx,ny,1)
265 
      varname='pressure'
 314 
      varname='pressure'
266 
      call input_var_wrf(fid,varname,p3 ,nx,ny,nz)
 315 
      call input_var_wrf(fid,varname,p3       ,nx,ny,nz)
267 
      varname='geopot'
 316 
      varname='geopot'
268 
      call input_var_wrf(fid,varname,z3 ,nx,ny,nz)
 317 
      call input_var_wrf(fid,varname,z3       ,nx,ny,nz)
269 
      varname='geopots'
 318 
      varname='geopots'
270 
      call input_var_wrf(fid,varname,zb ,nx,ny,1)
 319 
      call input_var_wrf(fid,varname,zb       ,nx,ny,1)
- 
 
 320 
      varname='COSALPHA'
- 
 
 321 
      call input_var_wrf(fid,varname,cosalpha ,nx,ny,1)
- 
 
 322 
      varname='SINALPHA'
- 
 
 323 
      call input_var_wrf(fid,varname,sinalpha ,nx,ny,1)
- 
 
 324 
 
- 
 
 325 
c     Determine the grid rotation <alpha> from <cosalpha> and <sinalpha>
- 
 
 326 
      do i=1,nx
- 
 
 327 
        do j=1,ny
- 
 
 328 
           alpha(i,j) = rad2deg * atan2(sinalpha(i,j),cosalpha(i,j))
- 
 
 329 
        enddo
- 
 
 330 
      enddo
271 
 
 331 
 
272 
c     Transform longitude depending on <anglemode>
 332 
c     Transform longitude depending on <anglemode>
273 
      if ( anglemode.eq.'greenwich' ) then
 333 
      if ( anglemode.eq.'greenwich' ) then
274 
         do i=1,nx
 334 
         do i=1,nx
275 
           do j=1,ny
 335 
           do j=1,ny
Line 407... Line 467...
407 
 
 467 
 
408 
 
 468 
 
409 
	enddo
 469 
	enddo
410 
      enddo
 470 
      enddo
411 
 
 471 
 
412 
c     Normalize output and set miising data flag
 472 
c     Normalize output and set missing data flag
413 
      do i=1,nlon
 473 
      do i=1,nlon
414 
         do j=1,nlat
 474 
         do j=1,nlat
415 
            if ( xc(i,j).gt.0. ) then
 475 
            if ( xc(i,j).gt.0. ) then
416 
               x(i,j) = x(i,j)/xc(i,j)
 476 
               x(i,j) = x(i,j)/xc(i,j)
417 
            else
 477 
            else
Line 457... Line 517...
457 
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
 517 
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
458 
      cdfname = mapfile
 518 
      cdfname = mapfile
459 
      varname = 'TOPO'
 519 
      varname = 'TOPO'
460 
      call writecdf2D(cdfname,varname,zb,0.,
 520 
      call writecdf2D(cdfname,varname,zb,0.,
461 
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 521 
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
- 
 
 522 
      cdfname = mapfile
- 
 
 523 
      varname = 'COSALPHA'
- 
 
 524 
      call writecdf2D(cdfname,varname,cosalpha,0.,
- 
 
 525 
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
- 
 
 526 
      cdfname = mapfile
- 
 
 527 
      varname = 'SINALPHA'
462 
 
 - 
 
- 
 
 528 
      call writecdf2D(cdfname,varname,sinalpha,0.,
- 
 
 529 
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
- 
 
 530 
      cdfname = mapfile
- 
 
 531 
      varname = 'ALPHA'
- 
 
 532 
      call writecdf2D(cdfname,varname,alpha,0.,
- 
 
 533 
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
463 
 
 534 
 
464 
      endif
 535 
      endif
465 
 
 536 
 
466 
 
 - 
 
467 
c     ------------------------------------------------------------
 537 
c     ------------------------------------------------------------
468 
c     Convert a lat/lon/z list to a x/y/z list
 538 
c     Convert a lat/lon/z list to a x/y/z list
469 
c     ------------------------------------------------------------
 539 
c     ------------------------------------------------------------
470
540
471 
      if ( mode.eq.'-ll2xy' ) then
 541 
      if ( mode.eq.'-ll2xy' ) then
Line 580... Line 650...
580 
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
 650 
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
581 
      allocate(yy0(ntra*ntim),stat=stat)
 651 
      allocate(yy0(ntra*ntim),stat=stat)
582 
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
 652 
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
583 
      allocate(zz0(ntra*ntim),stat=stat)
 653 
      allocate(zz0(ntra*ntim),stat=stat)
584 
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 654 
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
585 
 
 655
586 
c     Get format of input and output file
 656 
c     Get format of input and output file
587 
      call mode_tra(inpmode,inpfile)
 657 
      call mode_tra(inpmode,inpfile)
588 
      call mode_tra(outmode,outfile)
 658 
      call mode_tra(outmode,outfile)
589 
      if ( outmode.eq.-1) outmode=inpmode
 659 
      if ( outmode.eq.-1) outmode=inpmode
590 
 
 660 
 
Line 667... Line 737...
667 
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
 737 
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
668 
         call close_tra(fid,outmode)
 738 
         call close_tra(fid,outmode)
669 
      endif
 739 
      endif
670 
 
 740 
 
671 
      endif
 741 
      endif
- 
 
 742
- 
 
 743 
c     ------------------------------------------------------------
- 
 
 744 
c     Convert a evctor from grid-relative to lon/lat-relative components
- 
 
 745 
c     ------------------------------------------------------------
- 
 
 746
- 
 
 747 
      if ( mode.eq.'-vec2ll' ) then
- 
 
 748 
 
- 
 
 749 
c     Allocate memory for input and output lists
- 
 
 750 
      allocate(tra(ntra,ntim,ncol),stat=stat)
- 
 
 751 
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
- 
 
 752
- 
 
 753 
c     Get format of input and output file
- 
 
 754 
      call mode_tra(inpmode,inpfile)
- 
 
 755 
      call mode_tra(outmode,outfile)
- 
 
 756 
      if ( outmode.eq.-1) outmode=inpmode
- 
 
 757 
 
- 
 
 758 
c     Read coordinates from file - Format (x,y,lev)
- 
 
 759 
      if (inpmode.eq.-1) then
- 
 
 760 
         print*,' ERROR: vector transformation only for trajectory file'
- 
 
 761 
         stop
- 
 
 762 
      endif
- 
 
 763 
 
- 
 
 764 
c     Read coordinates from trajectory file
- 
 
 765 
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
- 
 
 766 
      call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
- 
 
 767 
      call close_tra(fid,inpmode)
- 
 
 768 
      if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
- 
 
 769 
            print*,' ERROR: input must be in x/y format... Stop'
- 
 
 770 
            stop
- 
 
 771 
      endif
- 
 
 772 
 
- 
 
 773 
c     Check whether vector components are on trajectory file
- 
 
 774 
      indu   = -1
- 
 
 775 
      indv   = -1
- 
 
 776 
      do i=1,ncol
- 
 
 777 
         if ( vars(i).eq.vectorx ) indu = i
- 
 
 778 
         if ( vars(i).eq.vectory ) indv = i
- 
 
 779 
      enddo
- 
 
 780 
      if ( (indu.lt.1).or.(indv.lt.1) ) then
- 
 
 781 
        print*,' WARNING: VECTOR cannot be transformed... '
- 
 
 782 
        stop
- 
 
 783 
      endif
- 
 
 784 
 
- 
 
 785 
c     Transform vectors
- 
 
 786 
      do i=1,ntra
- 
 
 787 
        do j=1,ntim
- 
 
 788 
 
- 
 
 789 
c         Check whether the position is valid
- 
 
 790 
          if ( ( abs(tra(i,j,2)-mdv).lt.eps ).or.
- 
 
 791 
     >         ( abs(tra(i,j,3)-mdv).lt.eps ) ) goto 100
- 
 
 792
- 
 
 793 
c         Get cosalpa and sinalpha at position
- 
 
 794 
          rid = (tra(i,j,2)-xmin)/dx+1.
- 
 
 795 
          rjd = (tra(i,j,3)-ymin)/dy+1.
- 
 
 796 
          ca0 = int_index3(cosalpha,nx,ny,1,rid,rjd,1.,mdv)
- 
 
 797 
          sa0 = int_index3(sinalpha,nx,ny,1,rid,rjd,1.,mdv)
- 
 
 798 
 
- 
 
 799 
c         Get grid-relative vector at position
- 
 
 800 
          uu0 = tra(i,j,indu)
- 
 
 801 
          vv0 = tra(i,j,indv)
- 
 
 802
- 
 
 803 
c         Get the new vector components
- 
 
 804 
          tra(i,j,indu) = ca0 * uu0 - sa0 * vv0
- 
 
 805 
          tra(i,j,indv) = sa0 * uu0 + ca0 * vv0
- 
 
 806 
 
- 
 
 807 
c         Exit point for loop over all tra's and tim's
- 
 
 808 
 100      continue
- 
 
 809 
 
- 
 
 810 
        enddo
- 
 
 811 
      enddo
- 
 
 812 
 
- 
 
 813 
c     Write output to trajectory
- 
 
 814 
      string = vars(indu)
- 
 
 815 
      len    = len_trim(string)
- 
 
 816 
      char0  = string((len-2):(len-2))
- 
 
 817 
      char1  = string((len-1):(len-1))
- 
 
 818 
      char2  = string(len:len)
- 
 
 819 
      if ( (char0.eq.'.').and.(char1.eq.'x').and.(char2.eq.'y') ) then
- 
 
 820 
         vars(indu) = trim(string(1:len-3))//'ll'
- 
 
 821 
      else
- 
 
 822 
         vars(indu) = trim(string)//'.ll'
- 
 
 823 
      endif
- 
 
 824 
      string = vars(indv)
- 
 
 825 
      len    = len_trim(string)
- 
 
 826 
      char0  = string((len-2):(len-2))
- 
 
 827 
      char1  = string((len-1):(len-1))
- 
 
 828 
      char2  = string(len:len)
- 
 
 829 
      if ( (char0.eq.'.').and.(char1.eq.'x').and.(char2.eq.'y') ) then
- 
 
 830 
         vars(indv) = trim(string(1:len-3))//'ll'
- 
 
 831 
      else
- 
 
 832 
         vars(indv) = trim(string)//'.ll'
- 
 
 833 
      endif
- 
 
 834 
      call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
- 
 
 835 
      call write_tra(fid,tra,ntra,ntim,ncol,outmode)
- 
 
 836 
      call close_tra(fid,outmode)
- 
 
 837 
 
- 
 
 838 
      endif
- 
 
 839
- 
 
 840 
c     ------------------------------------------------------------
- 
 
 841 
c     Convert a evctor from lat/lon-relative to grid-relative components
- 
 
 842 
c     ------------------------------------------------------------
- 
 
 843
- 
 
 844 
      if ( mode.eq.'-vec2xy' ) then
- 
 
 845 
 
- 
 
 846 
c     Allocate memory for input and output lists
- 
 
 847 
      allocate(tra(ntra,ntim,ncol),stat=stat)
- 
 
 848 
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
- 
 
 849
- 
 
 850 
c     Get format of input and output file
- 
 
 851 
      call mode_tra(inpmode,inpfile)
- 
 
 852 
      call mode_tra(outmode,outfile)
- 
 
 853 
      if ( outmode.eq.-1) outmode=inpmode
- 
 
 854 
 
- 
 
 855 
c     Read coordinates from file - Format (x,y,lev)
- 
 
 856 
      if (inpmode.eq.-1) then
- 
 
 857 
         print*,' ERROR: vector transformation only for trajectory file'
- 
 
 858 
         stop
- 
 
 859 
      endif
- 
 
 860 
 
- 
 
 861 
c     Read coordinates from trajectory file
- 
 
 862 
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
- 
 
 863 
      call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
- 
 
 864 
      call close_tra(fid,inpmode)
- 
 
 865 
      if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
- 
 
 866 
            print*,' ERROR: input must be in x/y format... Stop'
- 
 
 867 
            stop
- 
 
 868 
      endif
- 
 
 869 
 
- 
 
 870 
c     Check whether vector components are on trajectory file
- 
 
 871 
      indu   = -1
- 
 
 872 
      indv   = -1
- 
 
 873 
      do i=1,ncol
- 
 
 874 
         if ( vars(i).eq.vectorx ) indu = i
- 
 
 875 
         if ( vars(i).eq.vectory ) indv = i
- 
 
 876 
      enddo
- 
 
 877 
      if ( (indu.lt.1).or.(indv.lt.1) ) then
- 
 
 878 
        print*,' WARNING: VECTOR cannot be transformed... '
- 
 
 879 
        stop
- 
 
 880 
      endif
- 
 
 881 
 
- 
 
 882 
c     Transform vectors
- 
 
 883 
      do i=1,ntra
- 
 
 884 
        do j=1,ntim
- 
 
 885 
 
- 
 
 886 
c         Check whether the position is valid
- 
 
 887 
          if ( ( abs(tra(i,j,2)-mdv).lt.eps ).or.
- 
 
 888 
     >         ( abs(tra(i,j,3)-mdv).lt.eps ) ) goto 110
- 
 
 889
- 
 
 890 
c         Get cosalpa and sinalpha at position
- 
 
 891 
          rid = (tra(i,j,2)-xmin)/dx+1.
- 
 
 892 
          rjd = (tra(i,j,3)-ymin)/dy+1.
- 
 
 893 
          ca0 = int_index3(cosalpha,nx,ny,1,rid,rjd,1.,mdv)
- 
 
 894 
          sa0 = int_index3(sinalpha,nx,ny,1,rid,rjd,1.,mdv)
- 
 
 895 
 
- 
 
 896 
c         Get grid-relative vector at position
- 
 
 897 
          uu0 = tra(i,j,indu)
- 
 
 898 
          vv0 = tra(i,j,indv)
- 
 
 899
- 
 
 900 
c         Get the new vector components
- 
 
 901 
          tra(i,j,indu) =  ca0 * uu0 + sa0 * vv0
- 
 
 902 
          tra(i,j,indv) = -sa0 * uu0 + ca0 * vv0
- 
 
 903 
 
- 
 
 904 
c         Exit point for loop over all tra's and tim's
- 
 
 905 
 110      continue
- 
 
 906 
 
- 
 
 907 
        enddo
- 
 
 908 
      enddo
- 
 
 909 
 
- 
 
 910 
c     Write output to trajectory
- 
 
 911 
      string = vars(indu)
- 
 
 912 
      len    = len_trim(string)
- 
 
 913 
      char0  = string((len-2):(len-2))
- 
 
 914 
      char1  = string((len-1):(len-1))
- 
 
 915 
      char2  = string(len:len)
- 
 
 916 
      if ( (char0.eq.'.').and.(char1.eq.'l').and.(char2.eq.'l') ) then
- 
 
 917 
         vars(indu) = trim(string(1:len-3))//'.xy'
- 
 
 918 
      else
- 
 
 919 
         vars(indu) = trim(string)//'.xy'
- 
 
 920 
      endif
- 
 
 921 
      string = vars(indv)
- 
 
 922 
      len    = len_trim(string)
- 
 
 923 
      char0  = string((len-2):(len-2))
- 
 
 924 
      char1  = string((len-1):(len-1))
- 
 
 925 
      char2  = string(len:len)
- 
 
 926 
      if ( (char0.eq.'.').and.(char1.eq.'l').and.(char2.eq.'l') ) then
- 
 
 927 
         vars(indv) = trim(string(1:len-3))//'.xy'
- 
 
 928 
      else
- 
 
 929 
         vars(indv) = trim(string)//'.xy'
- 
 
 930 
      endif
- 
 
 931 
      call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
- 
 
 932 
      call write_tra(fid,tra,ntra,ntim,ncol,outmode)
- 
 
 933 
      call close_tra(fid,outmode)
- 
 
 934 
 
- 
 
 935 
      endif
672 
 
 936 
 
673 
c     ------------------------------------------------------------
 937 
c     ------------------------------------------------------------
674 
c     Convert a single lat/lon/z list to a x/y/z list
 938 
c     Convert a single lat/lon/z list to a x/y/z list
675 
c     ------------------------------------------------------------
 939 
c     ------------------------------------------------------------
676 
 
 940 
 
Line 1025... Line 1289...
1025 
      print*,xpos,ypos,ppos
 1289 
      print*,xpos,ypos,ppos
1026 
 
 1290 
 
1027 
      endif
 1291 
      endif
1028 
 
 1292 
 
1029 
      end
 1293 
      end
1030 
 
 - 
 
1031 
 
 1294
1032 
c     ********************************************************************
 1295 
c     ********************************************************************
1033 
c     * GRIDDING SUBROUTINES                                             *
 1296 
c     * GRIDDING SUBROUTINES                                             *
1034 
c     ********************************************************************
 1297 
c     ********************************************************************
1035 
 
 1298 
 
1036 
c     ---------------------------------------------------------------------
 1299 
c     ---------------------------------------------------------------------