Subversion Repositories pvinversion.ecmwf

      PROGRAM rotate_grid

c     *******************************************************************
c     * Rotate to a grid with the PV anomaly at its center              *
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      rnx,rny
      real         rdx,rdy
      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_grid                                           *'
      print*,'*********************************************************'

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

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

       read(10,*) name,rnx
       if ( trim(name).ne.'ROT_NX') stop
       read(10,*) name,rny
       if ( trim(name).ne.'ROT_NY') stop
       read(10,*) name,rdx
       if ( trim(name).ne.'ROT_DX') stop
       read(10,*) name,rdy
       if ( trim(name).ne.'ROT_DY') 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*,rnx,rny,rdx,rdy
      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) = rnx
      out_vardim(2) = rny
      out_vardim(3) = in_nz
      out_varmin(1) = -real(rnx/2)*rdx
      out_varmin(2) = -real(rny/2)*rdy
      out_varmin(3) = in_varmin(3)
      out_varmax(1) = real(rnx/2)*rdx
      out_varmax(2) = real(rny/2)*rdy
      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       = rdx
      out_dy       = rdy
      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   = clon
      out_pollat   = clat
      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
      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='LAT'
      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='LON'
      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     Coriolis parameter
      varname='CORIOL'
      print*,'Rotating scalar : ',trim(varname)
      do i=1,in_nx
         do j=1,in_ny
            lat=in_ymin+real(j-1)*in_dy
            in_field3(i,j,1)=2.*omegaerd*sin(lat*degrad)
         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_coriol,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_coriol,ierr)
      if (ierr.ne.0) goto 997

c     X coordinate
      varname='X'
      print*,'Getting X grid : ',trim(varname)
      do j=1,out_ny
         out_x(out_nx/2,j)=0.
         do i=out_nx/2+1,out_nx
            out_x(i,j)=out_x(i-1,j)+
     >                 sdis(out_lon(i-1,j),out_lat(i-1,j),
     >                      out_lon(i  ,j),out_lat(i  ,j))
         enddo
         do i=out_nx/2-1,1,-1
            out_x(i,j)=out_x(i+1,j)-
     >                 sdis(out_lon(i+1,j),out_lat(i+1,j),
     >                      out_lon(i  ,j),out_lat(i  ,j))
         enddo
      enddo
      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_x,ierr)
      if (ierr.ne.0) goto 997

c     Y coordinate
      varname='Y'
      print*,'Getting Y grid : ',trim(varname)
      do i=1,out_nx
         out_y(i,out_ny/2)=0.
         do j=out_ny/2+1,out_ny
            out_y(i,j)=out_y(i,j-1)+
     >                 sdis(out_lon(i,j-1),out_lat(i,j-1),
     >                      out_lon(i  ,j),out_lat(i  ,j))
         enddo
         do j=out_ny/2-1,1,-1
            out_y(i,j)=out_y(i,j+1)-
     >                 sdis(out_lon(i,j+1),out_lat(i,j+1),
     >                      out_lon(i  ,j),out_lat(i  ,j))
         enddo
      enddo
      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_y,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 lat/lon grid'
      stop

 997  print*,'Z: Problems with output to rotated grid'
      stop

      end


c     ********************************************************************************
c     * SUBROUTINE: ROTATION TO LOCAL CARTESIAN COORDINATE SYSTEM                    *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Subroutine to get environment (vector+scalar)
c     --------------------------------------------------------------------------------

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

c     Get an input field <inar> and the center <clon,clat> of a structure. The input 
c     array is defined on a lat/lon grid with grid specification <dx,dy,xmin,ymin,nx,ny>. 
c     The input field is then rotated to an output grid <rdx,rdy,rxmin,rymin,rnx,rny>, 
c     whereby the structure <clon,clat> is at the center of the new grid and the output 
c     grid is rotated with rotation angle <rotation>. This routine rotates the components
c     of a vector field for <rotmode=1|2>: <rotmode=1> refers to the x xomponent of the 
c     field, <rotmode=2> to the y component. If <rotmode=0>, scalar field is assumed.

      implicit none

c     Declaration of input parameters
      integer     nx,ny,nz
      real        dx,dy,xmin,ymin
      real        inar(nx,ny,nz)
      real        clon,clat,rotation
      real        mdv
      integer     rotmode

c     Declaration of output parameters
      integer     rnx,rny,rnz
      real        rdx,rdy,rxmin,rymin
      real        outar(rnx,rny,rnz)

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         rlon1(rnx,rny),rlat1(rnx,rny)
      real         rlon2(rnx,rny),rlat2(rnx,rny)
      real         lon(rnx,rny),lat(rnx,rny)
      real         rotangle1(rnx,rny),rotangle2(rnx,rny)
      real         rotangle(rnx,rny)
      real         sinoutar(rnx,rny,rnz),cosoutar(rnx,rny,rnz)
      real         xmax,ymax
      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
      real         lonmean,latmean
      character*80 cdfname,varname,cstname
      real         vx1,vx2,vy1,vy2,angle,vx2min
      integer      s

c     Externals
      real         lmstolm,phstoph
      external     lmstolm,phstoph

c     Get rotated coordinates
      do i=1,rnx
         do j=1,rny               
            rlon1(i,j)=rxmin+real(i-1)*rdx
            rlat1(i,j)=rymin+real(j-1)*rdy
         enddo
      enddo

c     First coordinate transformation (make the local coordinate system parallel to equator)
      pollon=-180.
      pollat=90.+rotation
      do i=1,rnx
       do j=1,rny
         rlon2(i,j)=90.+lmstolm(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)
         rlat2(i,j)=phstoph(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon)            
       enddo
      enddo

c     Second coordinate transformation (make the local coordinate system parallel to equator)
      pollon=clon-180.
      if (pollon.lt.-180.) pollon=pollon+360.
      pollat=90.-clat
      do i=1,rnx
         do j=1,rny               
            lon(i,j)=lmstolm(rlat2(i,j),rlon2(i,j),pollat,pollon)
            lat(i,j)=phstoph(rlat2(i,j),rlon2(i,j),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,rnx-1
            do j=2,rny-1
               
c              Angle between latitude circles
               vx1=1.
               vy1=0.
               vx2min=(lon(i+1,j)-lon(i-1,j))*cos(pi180*lat(i,j))
               do s=-1,1,2
                  vx2=(lon(i+1,j)-lon(i-1,j)+real(s)*360.)*
     >                 cos(pi180*lat(i,j))
                  if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
               enddo
               vx2=vx2min
               vy2=lat(i+1,j)-lat(i-1,j)
               call getangle(vx1,vy1,vx2,vy2,angle)           
               rotangle1(i,j)=angle

c              Angle between longitude circles
               vx1=0.
               vy1=1.
               vx2min=(lon(i,j+1)-lon(i,j-1))*cos(pi180*lat(i,j))
               do s=-1,1,2
                  vx2=(lon(i+1,j)-lon(i-1,j)+real(s)*360.)*
     >                 cos(pi180*lat(i,j))
                  if (abs(vx2).lt.abs(vx2min)) vx2min=vx2
               enddo
               vx2=vx2min
               vy2=lat(i,j+1)-lat(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,rnx
            rotangle1(i,rny)=2.0*rotangle1(i,rny-1)-rotangle1(i,rny-2)
            rotangle1(i,1)  =2.0*rotangle1(i,2)-rotangle1(i,3)
            rotangle2(i,rny)=2.0*rotangle2(i,rny-1)-rotangle2(i,rny-2)
            rotangle2(i,1)  =2.0*rotangle2(i,2)-rotangle2(i,3)
         enddo
         do j=1,rny
            rotangle1(rnx,j)=2.0*rotangle1(rnx-1,j)-rotangle1(rnx-2,j)
            rotangle1(1,j)  =2.0*rotangle1(2,j)-rotangle1(3,j)
            rotangle2(rnx,j)=2.0*rotangle2(rnx-1,j)-rotangle2(rnx-2,j)
            rotangle2(1,j)  =2.0*rotangle2(2,j)-rotangle2(3,j)
         enddo   
         
c        Set the final rotation angle
         do i=1,rnx
            do j=1,rny
               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,rnx
         do j=1,rny
 100        if (lon(i,j).lt.xmin) then
               lon(i,j)=lon(i,j)+360.
               goto 100
            endif
 102        if (lon(i,j).gt.(xmin+real(nx-1)*dx)) then
               lon(i,j)=lon(i,j)-360.
               goto 102
            endif
         enddo
      enddo

c     Rotate the scalar or the vector component
      do i=1,rnx
         do j=1,rny
            do k=1,rnz   

c              Get index
               rindx=(lon(i,j)-xmin)/dx+1.
               rindy=(lat(i,j)-ymin)/dy+1.
               indx=int(rindx)
               indy=int(rindy)
               if ((indx.lt.1).or.(indx.gt.nx).or.
     >             (indy.lt.1).or.(indy.gt.ny)) then
                  print*,'lat/lon interpolation not possible....'
                  print*,i,j,k,indx,indy,lon(i,j),lat(i,j)
                  stop
               endif

c              Get inidices of left and upper neighbours
               indr=indx+1
               if (indr.gt.nx) indr=1
               indu=indy+1
               if (indu.gt.ny) 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

            enddo            
         enddo
      enddo      

c     Get components for tests
      if (test.eq.1) then
         do i=1,rnx
            do j=1,rny
               do k=1,rnz
                  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,rnx
         do j=1,rny
            do k=1,rnz
               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='RLON2'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlon2,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='RLAT2'
         call  writecdf2D(cdfname,cstname,
     >                    varname,rlat2,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 LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C**                IM 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:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
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:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHIS   = ZPIR18*PHIS
      ZLAMS   = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS   = ZPIR18*ZLAMS
 
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) -
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) +
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMSTOLM =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMSTOLAM =  90.0
            ELSE
              LMSTOLAM = -90.0
            ENDIF
      ELSE
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END


      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C**                EINEN 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:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE 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:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      SINPOL = SIN(ZPIR18*POLPHI)
      COSPOL = COS(ZPIR18*POLPHI)
      ZPHIS  = ZPIR18*PHIS
      ZLAMS  = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS  = ZPIR18*ZLAMS
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
 
      PHSTOPH = ZRPI18*ASIN(ARG)
 
      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