Subversion Repositories pvinversion.ecmwf

      PROGRAM rotate_lalo

c     *******************************************************************
c     * Rotate to a geographical latitude/longitude grid                *
c     * Michael Sprenger / Autumn 2006                                  *
c     *******************************************************************

      implicit none

c     -------------------------------------------------------------------
c     Declaration of parameters and variables
c     -------------------------------------------------------------------

c     Specification of input parameters
      real         clon,clat,rotation
      integer      nx,ny
      real         dx,dy,xmin,ymin
      integer      nfield
      character*80 fieldname(100)

c     Numerical and physical parameters
      real                 degrad
      parameter            (degrad=0.0174532925)
      real                 omegaerd
      parameter            (omegaerd=7.292e-5 )
      real                 eps
      parameter            (eps=0.001)

c     Variables for input Z file : height level
      character*80 in_cfn
      real         in_varmin(4),in_varmax(4),in_stag(4)
      integer      in_vardim(4)
      real         in_mdv
      integer      in_ndim
      integer      in_nx,in_ny,in_nz
      real         in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy
      integer      in_ntimes
      real         in_aklev(500),in_bklev(500)
      real         in_aklay(500),in_bklay(500)
      real         in_time
      real         in_pollon,in_pollat
      integer      in_nvars
      character*80 in_vnam(100)
      integer      in_idate(5)
      real,allocatable, dimension (:,:,:) :: in_field3
      real,allocatable, dimension (:,:,:) :: in_vect1
      real,allocatable, dimension (:,:,:) :: in_vect2

c     Variables for output Z file : height level
      character*80 out_cfn
      real         out_varmin(4),out_varmax(4),out_stag(4)
      integer      out_vardim(4)
      real         out_mdv
      integer      out_ndim
      integer      out_nx,out_ny,out_nz
      real         out_xmin,out_xmax,out_ymin,out_ymax,out_dx,out_dy
      integer      out_ntimes
      real         out_aklev(500),out_bklev(500)
      real         out_aklay(500),out_bklay(500)
      real         out_time
      real         out_pollon,out_pollat
      integer      out_idate(5)
      real,allocatable, dimension (:,:,:) :: out_field3
      real,allocatable, dimension (:,:,:) :: out_vect1
      real,allocatable, dimension (:,:,:) :: out_vect2
      real,allocatable, dimension (:,:)   :: out_lat,out_lon
      real,allocatable, dimension (:,:)   :: out_x,out_y
      real,allocatable, dimension (:,:)   :: out_coriol

c     Auxiliary variables
      integer      ifield
      integer      i,j,k
      integer      cdfid,cstid
      character*80 cfn
      integer      stat,ierr,isok
      real         time
      character*80 varname,cdfname,varname1,varname2
      integer      idate(5),stdate(5),datar(14)
      integer      rotmode
      integer      i1,i2,i3,i4
      integer      isvector
      real         lat
      character*80 name

c     Externals
      real         sdis
      external     sdis

c     -------------------------------------------------------------------
c     Preparations
c     -------------------------------------------------------------------
            
      print*,'*********************************************************'
      print*,'* rotate_lalo                                          *'
      print*,'*********************************************************'

c     Read parameter file
      open(10,file='fort.10')

       read(10,*) in_cfn
       read(10,*) out_cfn

       read(10,*) name,nx
       if ( trim(name).ne.'GEO_NX'  ) stop
       read(10,*) name,ny
       if ( trim(name).ne.'GEO_NY'  ) stop
       read(10,*) name,dx
       if ( trim(name).ne.'GEO_DX'  ) stop
       read(10,*) name,dy
       if ( trim(name).ne.'GEO_DY'  ) stop
       read(10,*) name,xmin
       if ( trim(name).ne.'GEO_XMIN') stop
       read(10,*) name,ymin
       if ( trim(name).ne.'GEO_YMIN') stop
       read(10,*) name,clon
       if ( trim(name).ne.'CLON'    ) stop
       read(10,*) name,clat
       if ( trim(name).ne.'CLAT'    ) stop
       read(10,*) name,rotation
       if ( trim(name).ne.'CROT'    ) stop
       
       read(10,*) nfield
       do i=1,nfield
          read(10,*) fieldname(i)
       enddo
      close(10)
      
      print*,clon,clat,rotation
      print*,nx,ny,dx,dy,xmin,ymin
      print*,trim(in_cfn),' -> ',trim(out_cfn)

c     Get grid description for input Z file : height level
      call cdfopn(in_cfn,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getcfn(cdfid,cfn,ierr)
      if (ierr.ne.0) goto 998
      call cdfopn(cfn,cstid,ierr)
      if (ierr.ne.0) goto 998
      call getvars(cdfid,in_nvars,in_vnam,ierr)
      if (ierr.ne.0) goto 998
      isok=0
      varname=fieldname(1)
      call check_varok (isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) goto 998
      call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim,
     >            in_varmin,in_varmax,in_stag,ierr)
      if (ierr.ne.0) goto 998
      in_nx  =in_vardim(1)
      in_ny  =in_vardim(2)
      in_xmin=in_varmin(1)
      in_ymin=in_varmin(2)
      call getlevs(cstid,in_nz,in_aklev,in_bklev,in_aklay,in_bklay,ierr)
      call getgrid(cstid,in_dx,in_dy,ierr)
      in_xmax=in_xmin+real(in_nx-1)*in_dx
      in_ymax=in_ymin+real(in_ny-1)*in_dy
      call gettimes(cdfid,in_time,in_ntimes,ierr)
      call getstart(cstid,in_idate,ierr)
      call getpole(cstid,in_pollon,in_pollat,ierr)
      call clscdf(cstid,ierr)
      call clscdf(cdfid,ierr)

c     Set grid description for output file : height level
      out_vardim(1) = nx
      out_vardim(2) = ny
      out_vardim(3) = in_nz
      out_varmin(1) = xmin
      out_varmin(2) = ymin
      out_varmin(3) = in_varmin(3)
      out_varmax(1) = xmin+real(nx-1)*dx
      out_varmax(2) = ymin+real(ny-1)*dy
      out_varmax(3) = in_varmax(3)
      do i=1,in_nz
         out_aklay(i) = in_aklay(i)
         out_bklay(i) = in_bklay(i)
         out_aklev(i) = in_aklev(i)
         out_bklev(i) = in_bklev(i)
      enddo
      out_dx       = dx
      out_dy       = dy
      out_time     = in_time
      out_ntimes   = in_ntimes
      out_ndim     = 4
      out_mdv      = in_mdv
      out_nx       = out_vardim(1)
      out_ny       = out_vardim(2)
      out_nz       = out_vardim(3)
      out_xmin     = out_varmin(1)
      out_ymin     = out_varmin(2)
      out_pollon   = 0.
      out_pollat   = 90.
      do i=1,5
         out_idate(i) = in_idate(i)
      enddo

c     Allocate memory for all fields
      allocate(in_field3(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_field3 ***'
      allocate(in_vect1(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_vect1 ***'
      allocate(in_vect2(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_vect2 ***'
      allocate(out_field3(out_nx,out_ny,out_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_field3 ***'
      allocate(out_vect1(out_nx,out_ny,out_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_vect1 ***'
      allocate(out_vect2(out_nx,out_ny,out_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_vect2 ***'
      allocate(out_lat(out_nx,out_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_lat ***'
      allocate(out_lon(out_nx,out_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_lon ***'
      allocate(out_x(out_nx,out_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_x ***'
      allocate(out_y(out_nx,out_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_y ***'
      allocate(out_coriol(out_nx,out_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_coriol ***'

c     Create constants file and output file (if requested) 
      datar(1)=out_nx
      datar(2)=out_ny
      datar(3)=int(1000.*out_varmax(2))
      datar(4)=int(1000.*out_varmin(1))
      datar(5)=int(1000.*out_varmin(2))
      datar(6)=int(1000.*out_varmax(1))
      datar(7)=int(1000.*out_dx)
      datar(8)=int(1000.*out_dy)
      datar(9)=out_nz
      datar(10)=1
      datar(11)=1      
      datar(12)=0      
      datar(13)=int(1000.*out_pollon) 
      datar(14)=int(1000.*out_pollat) 
         
      cfn=trim(out_cfn)//'_cst'
      call wricst(cfn,datar,out_aklev,out_bklev,
     >            out_aklay,out_bklay,out_idate)
      
      call crecdf(trim(out_cfn),cdfid,out_varmin,out_varmax,
     >           out_ndim,trim(cfn),ierr)
      if (ierr.ne.0) goto 997
      call clscdf(cdfid,ierr)

c     -----------------------------------------------------------------
c     Loop through all fields - rotate to new grid
c     -----------------------------------------------------------------

      do ifield=1,nfield

c        Check if scalar or vectorial field (X for scalar, U.V for vector)
         varname=fieldname(ifield)
         i1=1
         i2=1
         i3=0
         i4=0
 100     if (varname(i1:i1).eq.' ') then
            i1=i1+1
            goto 100
         endif
         i2=i1
 101     if ((varname(i2:i2).ne.' ').and.(varname(i2:i2).ne.'.')) then
            i2=i2+1
            goto 101
         endif
         if (varname(i2:i2).eq.'.') then
            i3=i2+1
 102        if (varname(i3:i3).eq.' ') then
               i3=i3+1
               goto 102
            endif
            i4=i3
 104        if (varname(i4:i4).ne.' ') then
               i4=i4+1
               goto 104
            endif
         endif
         if ((i3.ne.0).and.(i4.ne.0)) then
            isvector=1
            i2=i2-1
            varname1=varname(i1:i2)
            i4=i4-1         
            varname2=varname(i3:i4)
            print*,'Rotating vector : ',
     >             trim(varname1),' / ',trim(varname2)
         else
            isvector=0
            i2=i2-1
            varname=varname(i1:i2)
            print*,'Rotating scalar : ',trim(varname)
         endif
         
c        Rotate a scalar field
         if (isvector.eq.0) then

c           Read input field
            call cdfopn(in_cfn,cdfid,ierr)
            if (ierr.ne.0) goto 998
            call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim,
     >                  in_varmin,in_varmax,in_stag,ierr)
            in_nz=in_vardim(3)
            if (ierr.ne.0) goto 998
            call getdat(cdfid,varname,in_time,0,in_field3,ierr)
            if (ierr.ne.0) goto 998
            call clscdf(cdfid,ierr) 

c           Rotate to new coordinate system
            out_nz=in_nz
            call getenvir (clon,clat,rotation,0,
     >                     in_field3,in_dx,in_dy,
     >                     in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
     >                     out_field3,out_dx,out_dy,
     >                     out_xmin,out_ymin,out_nx,out_ny,out_nz)

c           Write output scalar field
            call cdfwopn(trim(out_cfn),cdfid,ierr)
            if (ierr.ne.0) goto 997
            out_vardim(3)=out_nz
            call putdef(cdfid,varname,4,out_mdv,out_vardim,
     >                  out_varmin,out_varmax,out_stag,ierr)         
            if (ierr.ne.0) goto 997
            call putdat(cdfid,varname,out_time,0,out_field3,ierr)
            if (ierr.ne.0) goto 997
            call clscdf(cdfid,ierr)

c        Rotate vector field
         else if (isvector.eq.1) then
               
c           Read input vector field
            call cdfopn(in_cfn,cdfid,ierr)
            if (ierr.ne.0) goto 998
            call getdef(cdfid,varname1,in_ndim,in_mdv,in_vardim,
     >                  in_varmin,in_varmax,in_stag,ierr)
            in_nz=in_vardim(3)
            if (ierr.ne.0) goto 998
            call getdat(cdfid,varname1,in_time,0,in_vect1,ierr)
            if (ierr.ne.0) goto 998
            call getdef(cdfid,varname2,in_ndim,in_mdv,in_vardim,
     >                  in_varmin,in_varmax,in_stag,ierr)
            in_nz=in_vardim(3)
            if (ierr.ne.0) goto 998
            call getdat(cdfid,varname2,in_time,0,in_vect2,ierr)
            if (ierr.ne.0) goto 998
            call clscdf(cdfid,ierr) 

c           Get new vector component in x direction
            out_nz=in_nz
            call getenvir (clon,clat,rotation,1,
     >                     in_vect1,in_dx,in_dy,
     >                     in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
     >                     out_field3,out_dx,out_dy,
     >                     out_xmin,out_ymin,out_nx,out_ny,out_nz)
            do i=1,out_nx
               do j=1,out_ny
                  do k=1,out_nz
                     out_vect1(i,j,k)=out_field3(i,j,k)
                  enddo
               enddo
            enddo
            call getenvir (clon,clat,rotation,2,
     >                     in_vect2,in_dx,in_dy,
     >                     in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
     >                     out_field3,out_dx,out_dy,
     >                     out_xmin,out_ymin,out_nx,out_ny,out_nz)
            do i=1,out_nx
               do j=1,out_ny
                  do k=1,out_nz
                     if ( (abs(out_vect1 (i,j,k)-out_mdv).gt.eps).and.
     >                    (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then
                        out_vect1(i,j,k)=out_vect1(i,j,k)-
     >                                   out_field3(i,j,k)
                     endif
                  enddo
               enddo
            enddo           

c           Get new vector component in y direction
            out_nz=in_nz
            call getenvir (clon,clat,rotation,2,
     >                     in_vect1,in_dx,in_dy,
     >                     in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
     >                     out_field3,out_dx,out_dy,
     >                     out_xmin,out_ymin,out_nx,out_ny,out_nz)
            do i=1,out_nx
               do j=1,out_ny
                  do k=1,out_nz
                     out_vect2(i,j,k)=out_field3(i,j,k)
                  enddo
               enddo
            enddo
            call getenvir (clon,clat,rotation,1,
     >                     in_vect2,in_dx,in_dy,
     >                     in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv,
     >                     out_field3,out_dx,out_dy,
     >                     out_xmin,out_ymin,out_nx,out_ny,out_nz)
            do i=1,out_nx
               do j=1,out_ny
                  do k=1,out_nz
                     if ( (abs(out_vect2 (i,j,k)-out_mdv).gt.eps).and.
     >                    (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then
                        out_vect2(i,j,k)=out_vect2(i,j,k)+
     >                                   out_field3(i,j,k)
                     endif
                  enddo
               enddo
            enddo           

c           Write output vector field
            call cdfwopn(trim(out_cfn),cdfid,ierr)
            if (ierr.ne.0) goto 997
            out_vardim(3)=out_nz
            call putdef(cdfid,varname1,4,out_mdv,out_vardim,
     >                  out_varmin,out_varmax,out_stag,ierr)         
            if (ierr.ne.0) goto 997
            call putdat(cdfid,varname1,out_time,0,out_vect1,ierr)
            if (ierr.ne.0) goto 997
            call putdef(cdfid,varname2,4,out_mdv,out_vardim,
     >                  out_varmin,out_varmax,out_stag,ierr)         
            if (ierr.ne.0) goto 997
            call putdat(cdfid,varname2,out_time,0,out_vect2,ierr)
            if (ierr.ne.0) goto 997
            call clscdf(cdfid,ierr)

         endif

      enddo

c     -----------------------------------------------------------------
c     Write additional fields: latitude, longitude, Coriolis parameter
c     -----------------------------------------------------------------

c     Open the output file
      call cdfwopn(trim(out_cfn),cdfid,ierr)
      if (ierr.ne.0) goto 997     

c     Geographical latitude
      varname='RLAT'
      print*,'Rotating scalar : ',trim(varname)
      do i=1,in_nx
         do j=1,in_ny
            in_field3(i,j,1)=in_ymin+real(j-1)*in_dy
         enddo
      enddo
      call getenvir (clon,clat,rotation,0,
     >               in_field3,in_dx,in_dy,
     >               in_xmin,in_ymin,in_nx,in_ny,1,in_mdv,
     >               out_lat,out_dx,out_dy,
     >               out_xmin,out_ymin,out_nx,out_ny,1)      
      out_vardim(3)=1
      call putdef(cdfid,varname,4,out_mdv,out_vardim,
     >            out_varmin,out_varmax,out_stag,ierr)         
      if (ierr.ne.0) goto 997
      call putdat(cdfid,varname,out_time,0,out_lat,ierr)
      if (ierr.ne.0) goto 997

c     Geographical longitude
      varname='RLON'
      print*,'Rotating scalar : ',trim(varname)
      do i=1,in_nx
         do j=1,in_ny
            in_field3(i,j,1)=in_xmin+real(i-1)*in_dx
         enddo
      enddo
      call getenvir (clon,clat,rotation,0,
     >               in_field3,in_dx,in_dy,
     >               in_xmin,in_ymin,in_nx,in_ny,1,in_mdv,
     >               out_lon,out_dx,out_dy,
     >               out_xmin,out_ymin,out_nx,out_ny,1)      
      out_vardim(3)=1
      call putdef(cdfid,varname,4,out_mdv,out_vardim,
     >            out_varmin,out_varmax,out_stag,ierr)         
      if (ierr.ne.0) goto 997
      call putdat(cdfid,varname,out_time,0,out_lon,ierr)
      if (ierr.ne.0) goto 997



c     Close output file
      call clscdf(cdfid,ierr)

c     -----------------------------------------------------------------
c     Exception handling and format specs
c     -----------------------------------------------------------------

      stop

 998  print*,'Z: Problems with input rotated grid'
      stop

 997  print*,'Z: Problems with output to lat/lon grid'
      stop

      end


c     ********************************************************************************
c     * SUBROUTINE: ROTATION TO LATITUDE/LONGITUDE COORDINATE SYSTEM                    *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Subroutine to get environment of strcof
c     --------------------------------------------------------------------------------

      SUBROUTINE getenvir (clon,clat,rotation,rotmode,
     >                     inar, rdx,rdy,rxmin,rymin,rnx,rny,rnz,mdv,
     >                     outar, dx, dy, xmin, ymin, nx, ny, nz)

c     Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate 
c     system.

      implicit none

c     Declaration of input parameters
      integer     rnx,rny,rnz
      real        rdx,rdy,rxmin,rymin
      real        inar(rnx,rny,rnz)
      real        clon,clat,rotation
      real        mdv
      integer     rotmode

c     Declaration of output parameters
      integer     nx,ny,nz
      real        dx,dy,xmin,ymin
      real        outar(nx,ny,nz)

c     Set numerical and physical constants
      real        g2r
      parameter   (g2r=0.0174533)
      real        pi180
      parameter   (pi180=3.14159265359/180.)
      real        eps
      parameter   (eps=0.0001)


c     Flag for test mode
      integer      test
      parameter    (test=0)
      character*80 testfile
      parameter    (testfile='ROTGRID')

c     Auxiliary variables 
      real         pollon,pollat
      integer      i,j,k,l
      real         rlon(nx,ny),rlat(nx,ny)
      real         rlon1(nx,ny),rlat1(nx,ny)
      real         lon(nx,ny),lat(nx,ny)
      real         rotangle1(nx,ny),rotangle2(nx,ny)
      real         rotangle(nx,ny)
      real         sinoutar(nx,ny,nz),cosoutar(nx,ny,nz)
      real         rxmax,rymax
      real         xind,yind,pind
      real         outval
      integer      ix,iy
      real         ax,ay,az,bx,by,bz,zx,zy,zz
      real         clon1,clat1,clon2,clat2
      real         rindx,rindy
      integer      indx,indy,indr,indu
      real         frac0i,frac0j,frac1i,frac1j
      character*80 cdfname,varname,cstname
      real         vx1,vx2,vy1,vy2,angle,vx2min
      integer      s

c     Externals
      real     lmtolms,phtophs
      external lmtolms,phtophs

c     Get geographical coordinates
      do i=1,nx
         do j=1,ny               
            lon(i,j)=xmin+real(i-1)*dx
            lat(i,j)=ymin+real(j-1)*dy
         enddo
      enddo

c     First rotation
      pollon=clon-180.
      if (pollon.lt.-180.) pollon=pollon+360.
      pollat=90.-clat
      do i=1,nx
         do j=1,ny               
            rlon1(i,j)=lmtolms(lat(i,j),lon(i,j),pollat,pollon)
            rlat1(i,j)=phtophs(lat(i,j),lon(i,j),pollat,pollon)            
         enddo
      enddo

c     Second coordinate transformation 
      pollon=-180.
      pollat=90.+rotation
      do i=1,nx
       do j=1,ny
         rlon(i,j)=90.+lmtolms(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)
         rlat(i,j)=phtophs(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)            
       enddo
      enddo

c     Get the angle between the rotated and non-rotated coordinate frame
      if ((rotmode.eq.1).or.(rotmode.eq.2)) then
         do i=2,nx-1
            do j=2,ny-1
               
c              Angle between latitude circles
               vx1=1.
               vy1=0.
               vx2min=(rlon(i+1,j)-rlon(i-1,j))*cos(pi180*rlat(i,j))
               do s=-1,1,2
                  vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)*
     >                 cos(pi180*rlat(i,j))
                  if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
               enddo
               vx2=vx2min
               vy2=rlat(i+1,j)-rlat(i-1,j)
               call getangle(vx1,vy1,vx2,vy2,angle)           
               rotangle1(i,j)=angle

c              Angle between longitude circles
               vx1=0.
               vy1=1.
               vx2min=(rlon(i,j+1)-rlon(i,j-1))*cos(pi180*rlat(i,j))
               do s=-1,1,2
                  vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)*
     >                 cos(pi180*rlat(i,j))
                  if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
               enddo
               vx2=vx2min
               vy2=rlat(i,j+1)-rlat(i,j-1)
               call getangle(vx1,vy1,vx2,vy2,angle)           
               rotangle2(i,j)=angle
               
            enddo
         enddo

c        Set the angle at the boundaries
         do i=1,nx
            rotangle1(i,ny)=2.0*rotangle1(i,ny-1)-rotangle1(i,ny-2)
            rotangle1(i,1) =2.0*rotangle1(i,2)-rotangle1(i,3)
            rotangle2(i,ny)=2.0*rotangle2(i,ny-1)-rotangle2(i,ny-2)
            rotangle2(i,1) =2.0*rotangle2(i,2)-rotangle2(i,3)
         enddo
         do j=1,ny
            rotangle1(nx,j)=2.0*rotangle1(nx-1,j)-rotangle1(nx-2,j)
            rotangle1(1,j) =2.0*rotangle1(2,j)-rotangle1(3,j)
            rotangle2(nx,j)=2.0*rotangle2(nx-1,j)-rotangle2(nx-2,j)
            rotangle2(1,j) =2.0*rotangle2(2,j)-rotangle2(3,j)
         enddo   
         
c        Set the final rotation angle
         do i=1,nx
            do j=1,ny
               rotangle(i,j)=0.5*(rotangle1(i,j)+rotangle2(i,j))
            enddo
         enddo

      endif

c     Bring longitude into domain of geographical grid (shift by 360 deg)
      do i=1,nx
         do j=1,ny
 100        if (rlon(i,j).lt.rxmin) then
               rlon(i,j)=rlon(i,j)+360.
               goto 100
            endif
 102        if (rlon(i,j).gt.(rxmin+real(rnx-1)*rdx)) then
               rlon(i,j)=rlon(i,j)-360.
               goto 102
            endif
         enddo
      enddo

c     Rotate the scalar or the vector component
      do i=1,nx
         do j=1,ny
            do k=1,nz   

c              Get index
               rindx=(rlon(i,j)-rxmin)/rdx+1.
               rindy=(rlat(i,j)-rymin)/rdy+1.
               indx=int(rindx)
               indy=int(rindy)
               if ((indx.lt.1).or.(indx.gt.rnx).or.
     >             (indy.lt.1).or.(indy.gt.rny)) then
                  outar(i,j,k)=mdv
                  goto 103
               endif

c              Get inidices of left and upper neighbours
               indr=indx+1
               if (indr.gt.rnx) indr=1
               indu=indy+1
               if (indu.gt.rny) indu=ny
               
c              Do linear interpolation
               if ( ( abs(inar(indx ,indy, k)-mdv).gt.eps ).and.
     &              ( abs(inar(indx ,indu, k)-mdv).gt.eps ).and.
     &              ( abs(inar(indr ,indy, k)-mdv).gt.eps ).and.
     &              ( abs(inar(indr ,indu, k)-mdv).gt.eps ) ) then
                  frac0i=rindx-float(indx)
                  frac0j=rindy-float(indy)
                  frac1i=1.-frac0i
                  frac1j=1.-frac0j
                  outval = inar(indx ,indy, k ) * frac1i * frac1j
     &                   + inar(indx ,indu, k ) * frac1i * frac0j
     &                   + inar(indr ,indy, k ) * frac0i * frac1j
     &                   + inar(indr ,indu, k ) * frac0i * frac0j
               else
                  outval=mdv
               endif                  

c              Update output array
               outar(i,j,k)=outval

c              Next
 103           continue

            enddo            
         enddo
      enddo      

c     Get components for tests
      if (test.eq.1) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  cosoutar(i,j,k)=outar(i,j,k)*cos(pi180*rotangle(i,j))
                  sinoutar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j))
               enddo
            enddo
         enddo         
      endif

c     Get correct component of rotated field
      do i=1,nx
         do j=1,ny
            do k=1,nz
               if ( abs(outar(i,j,k)-mdv).gt.eps ) then
                  if (rotmode.eq.1) then
                     outar(i,j,k)= outar(i,j,k)*cos(pi180*rotangle(i,j))
                  else if (rotmode.eq.2) then
                     outar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j))
                  else if (rotmode.eq.0) then
                     outar(i,j,k)=outar(i,j,k)
                  endif
               endif
            enddo
         enddo
      enddo

c     Test mode: Write grids to cdf file
      if (test.eq.1) then
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='RLON1'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlon1,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    1,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='RLON'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlon,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='LON'
         call  writecdf2D(cdfname,cstname,
     >                    varname,lon,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='RLAT1'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlat1,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='RLAT'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlat,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='LAT'
         call  writecdf2D(cdfname,cstname,
     >                    varname,lat,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='ANGLE1'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rotangle1,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='ANGLE2'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rotangle2,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='U'
         call  writecdf2D(cdfname,cstname,
     >                    varname,cosoutar,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)
         cdfname=testfile
         cstname=trim(testfile)//'_cst'
         varname='V'
         call  writecdf2D(cdfname,cstname,
     >                    varname,sinoutar,0.,
     >                    rdx,rdy,rxmin,rymin,rnx,rny,
     >                    0,1)

      endif

      END


c     --------------------------------------------------------------------------------
c     Auxiliary routines: angle between two vectors
c     --------------------------------------------------------------------------------

      SUBROUTINE getangle (vx1,vy1,vx2,vy2,angle)

c     Given two vectors <vx1,vy1> and <vx2,vy2>, determine the angle (in deg)
c     between the two vectors.

      implicit none

c     Declaration of subroutine parameters
      real vx1,vy1
      real vx2,vy2
      real angle

c     Auxiliary variables and parameters
      real len1,len2,len3
      real val1,val2,val3
      real pi
      parameter (pi=3.14159265359)

      len1=sqrt(vx1*vx1+vy1*vy1)
      len2=sqrt(vx2*vx2+vy2*vy2)

      if ((len1.gt.0.).and.(len2.gt.0.)) then
         vx1=vx1/len1
         vy1=vy1/len1
         vx2=vx2/len2
         vy2=vy2/len2

         val1=vx1*vx2+vy1*vy2
         val2=-vy1*vx2+vx1*vy2

         len3=sqrt(val1*val1+val2*val2)

         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
            val3=acos(val1/len3)
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
            val3=pi-acos(abs(val1)/len3)
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
            val3=-acos(val1/len3)
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
            val3=-pi+acos(abs(val1)/len3)
         endif
      else
         val3=0.
      endif

      angle=180./pi*val3

      END

c     --------------------------------------------------------------------------------
c     Transformation routine: LMSTOLM and PHSTOPH from library gm2em
c     --------------------------------------------------------------------------------

      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL =     ZPIR18*POLLAM
      ZPHI    =     ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
 
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMTOLMS =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMTOLMS =  90.0
            ELSE
              LMTOLMS = -90.0
            ENDIF
      ELSE
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END


      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHI    = ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 
      RETURN
      END


c     ---------------------------------------------------------
c     Spherical distance between two lat/lon points
c     ---------------------------------------------------------

      real function sdis(xp,yp,xq,yq)
c
c     calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
c
      real      re
      parameter (re=6370.)
      real      xp,yp,xq,yq,arg
 
      arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
      if (arg.lt.-1.) arg=-1.
      if (arg.gt.1.) arg=1.
      sdis=re*acos(arg)
 
      end


c     ********************************************************************************
c     * NETCDF INPUT/OUTPUT                                                          *
c     ********************************************************************************

c     ---------------------------------------------------------
c     Subroutines to write the netcdf output file
c     ---------------------------------------------------------

      subroutine writecdf2D(cdfname,cstname,
     >                      varname,arr,ctime,
     >                      dx,dy,xmin,ymin,nx,ny,
     >                      crefile,crevar)

c     Create and write to the netcdf file <cdfname>. The variable
c     with name <varname> and with time <time> is written. The data
c     are in the two-dimensional array <arr>. The list <dx,dy,xmin,
c     ymin,nx,ny> specifies the output grid. The flags <crefile> and
c     <crevar> determine whether the file and/or the variable should
c     be created. 1: create / 0: not created

      IMPLICIT NONE

c     Declaration of input parameters
      character*80 cdfname,cstname,varname
      integer      nx,ny
      real         arr(nx,ny)
      real         dx,dy,xmin,ymin
      integer      ctime
      integer      crefile,crevar

c     Further variables
      real         varmin(4),varmax(4),stag(4)
      integer      ierr,cdfid,ndim,vardim(4)
      real         mdv
      integer      datar(14),stdate(5)
      integer      i
      real         time
      
C     Definitions 
      varmin(1)=xmin
      varmin(2)=ymin
      varmin(3)=1050.
      varmax(1)=xmin+real(nx-1)*dx
      varmax(2)=ymin+real(ny-1)*dy
      varmax(3)=1050.
      ndim=4
      vardim(1)=nx
      vardim(2)=ny
      vardim(3)=1
      stag(1)=0.
      stag(2)=0.
      stag(3)=0.
      mdv=-999.98999
      time=real(ctime)

C     Create the file
      if (crefile.eq.0) then
         call cdfwopn(cdfname,cdfid,ierr)
         if (ierr.ne.0) goto 906
      else if (crefile.eq.1) then
         call crecdf(cdfname,cdfid,
     >        varmin,varmax,ndim,cstname,ierr)
         if (ierr.ne.0) goto 903 

C        Write the constants file
         datar(1)=vardim(1)
         datar(2)=vardim(2)
         datar(3)=int(1000.*varmax(2))
         datar(4)=int(1000.*varmin(1))
         datar(5)=int(1000.*varmin(2))
         datar(6)=int(1000.*varmax(1))
         datar(7)=int(1000.*dx)
         datar(8)=int(1000.*dy)
         datar(9)=1
         datar(10)=1
         datar(11)=0            ! data version
         datar(12)=0            ! cstfile version
         datar(13)=0            ! longitude of pole
         datar(14)=90000        ! latitude of pole     
         do i=1,5
            stdate(i)=0
         enddo
         call wricst(cstname,datar,0.,0.,0.,0.,stdate)
      endif

c     Write the data to the netcdf file, and close the file
      if (crevar.eq.1) then
         call putdef(cdfid,varname,ndim,mdv,
     >        vardim,varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 904
      endif
      call putdat(cdfid,varname,time,0,arr,ierr)
      if (ierr.ne.0) goto 905
      call clscdf(cdfid,ierr)

      return
       
c     Error handling
 903  print*,'*** Problem to create netcdf file ***'
      stop
 904  print*,'*** Problem to write definition ***'
      stop
 905  print*,'*** Problem to write data ***'
      stop
 906  print*,'*** Problem to open netcdf file ***'
      stop

      END

c     ----------------------------------------------------------------
c     Check whether variable is found on netcdf file
c     ----------------------------------------------------------------

      subroutine check_varok (isok,varname,varlist,nvars)

c     Check whether the variable <varname> is in the list <varlist(nvars)>.
c     If this is the case, <isok> is incremented by 1. Otherwise <isok>
c     keeps its value.

      implicit none

c     Declaraion of subroutine parameters
      integer      isok
      integer      nvars
      character*80 varname
      character*80 varlist(nvars)

c     Auxiliary variables
      integer      i

c     Main
      do i=1,nvars
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
      enddo

      end