Subversion Repositories lagranto.um

Rev

Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 

C     ********************************************************************

      program ptos

C     ********************************************************************

C     Purpose
C     -------
C
C       Calculates secondary data files from primary data files
C       (based upon IVE-routines).
C       The program is invoked by the shell-script p2s.
C
C     Author
C     ------
C
C       H. Wernli       April 96
C
C     Modifications
C     -------------
C
CSue Adpated to calculate f*geostrophic wind components and geostrophic 
C momentum components. 
C lat, lon, rlat, rlon modified to arrays so that can be precalculated
C lat in radians, rlat in degrees (similarly for lon)
C     ********************************************************************


      include "um_dims.inc"

      real      time(ntmax)
      real      sp(nxmax*nymax),cl(nxmax*nymax),f(nxmax*nymax)
      real      oro(nxmax*nymax),ps(nxmax*nymax)
      real      var(nxmax*nymax*nzmax),th(nxmax*nymax*nzmax),
     >          pv(nxmax*nymax*nzmax),the(nxmax*nymax*nzmax),
     >          thw(nxmax*nymax*nzmax),
     >          rh(nxmax*nymax*nzmax),dhr(nxmax*nymax*nzmax),
     >          tt(nxmax*nymax*nzmax),qq(nxmax*nymax*nzmax),
     >          uu(nxmax*nymax*nzmax),vv(nxmax*nymax*nzmax),
     >          ww(nxmax*nymax*nzmax),fug(nxmax*nymax*nzmax),
     >          fvg(nxmax*nymax*nzmax),dspdx(nxmax*nymax),
     >          dspdy(nxmax*nymax),dvardx(nxmax*nymax*nzmax),
     >          dvardy(nxmax*nymax*nzmax),Mg(nxmax*nymax*nzmax),
     >          Ng(nxmax*nymax*nzmax)
      character*80 cdfnam,cstnam,cdf_file
      integer   cdfid,cdfid1,cstid,ierr,ndim,vardim(4)
      real      dx,dy,mdv,varmin(4),varmax(4),stag(4)
      real      aklev(nzmax),bklev(nzmax),aklay(nzmax),bklay(nzmax),
     >          ak(nzmax),bk(nzmax)
      real lon_eq,lat_eq,u_ll,v_ll,Mg_ll,Ng_ll,lon1,lat1
      integer   nx,ny,nz,ntimes,i,j,k,n,kk,jj
      integer   stdate(5)
      integer   mode,qmode,zdef
      real      rlat(nxmax*nymax),rlon(nxmax*nymax),
     >          lat(nxmax*nymax),lon(nxmax*nymax)
      real      pollon,pollat,yphys,xphys
      real      phstoph,lmstolm
      logical   prelev

      real      pi
      data      pi /3.141592654/

      integer iw, jw, kkw, w
      parameter (w=7)
      real weight(w,w), sumweight, tmpu(nxmax*nymax*nzmax), 
     > tmpv(nxmax*nymax*nzmax), ug(nxmax*nymax*nzmax), 
     > vg(nxmax*nymax*nzmax)
      write(*,*)'*** start of program ptos ***'

C     Read filename

      read(9,10)cdfnam
      write(*,*) 'cdfnam is ',cdfnam
   10 format(a13)

C     Read mode and qmode
  
      read(9,*)mode
      read(9,*)qmode
      if (mode.eq.10) read(9,*)zdef

C     Open files and get infos about data domain

      if (mode.eq.10) then
        call cdfwopn('P'//trim(cdfnam),cdfid1,ierr)
      else
        call cdfopn('P'//trim(cdfnam),cdfid1,ierr)
      endif
      call getcfn(cdfid1,cstnam,ierr)
      call cdfopn(trim(cstnam),cstid,ierr)

      call getdef(cdfid1,'T',ndim,mdv,vardim,varmin,varmax,stag,ierr)
      if (ierr.ne.0) goto 920
      mdv=-999.98999

C     Get the levels, pole, etc.
 
      nx=vardim(1)
      ny=vardim(2)
      nz=vardim(3)
 
      call getgrid(cstid,dx,dy,ierr)
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
      call getpole(cstid,pollon,pollat,ierr)
      call getstart(cstid,stdate,ierr)
C      write(*,*) 'dx ',dx,' dy ',dy
C      write(*,*) 'pollon ', pollon,' pollat ',pollat
 
C     Determine if data is on pressure or model levels
 
      prelev=.true.
      do k=1,nz
        if (bklev(k).ne.0.) prelev=.false.
      enddo
C      write(*,*) ' prelev ',prelev,'aklev ',aklev,' bklev ',bklev

CSue Calculate real lats and lons
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
        do j=1,ny
          do i=1,nx
            jj=i+(j-1)*nx
            rlat(jj)=varmin(2)+(j-1)*dy
            rlon(jj)=varmin(1)+(i-1)*dx      
            yphys=phstoph(rlat(jj),rlon(jj),pollat,pollon)
C           if I use sind(lat in deg): troubles at the N-pole
            lat(jj)=2.*pi*yphys/360.
            xphys=lmstolm(rlat(jj),rlon(jj),pollat,pollon)
            lon(jj)=2.*pi*xphys/360.
          enddo
        enddo
        else
        do j=1,ny
          do i=1,nx
            jj=i+(j-1)*nx
            lon(jj)=varmin(1)+(i-1)*dx 
            lat(jj)=varmin(2)+(j-1)*dy
          enddo
        enddo
      endif

C ---------------------------------

C     Calculate cos(latitude) array and the coriolis parameter

      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
        do j=1,ny
          do i=1,nx
            jj=i+(j-1)*nx
            cl(i+(j-1)*nx)=cos(pi*rlat(jj)/180.)
            f(i+(j-1)*nx)=0.000145444*sin(lat(jj))
          enddo
        enddo
      else
        do j=1,ny
          do i=1,nx
            cl(i+(j-1)*nx)=cos(lat(jj))
            f(i+(j-1)*nx)=0.000145444*sin(lat(jj))
          enddo
        enddo
      endif

C     Determine if data is on levels or layers

      if (stag(3).eq.-0.5) then
        do k=1,nz
          ak(k)=aklay(k)
          bk(k)=bklay(k)
        enddo
      else
        do k=1,nz
          ak(k)=aklev(k)
          bk(k)=bklev(k)
        enddo
      endif

C     Get all the fields
 
      call gettimes(cdfid1,time,ntimes,ierr)

C     Loop over all times
      write(*,*) 'ntimes = ',ntimes

      do n=1,ntimes
 
        if (.not.prelev) then
          call getdat(cdfid1,'PS',time(n),0,sp,ierr)
          if (ierr.ne.0) goto 921
        else
          call getdat(cdfid1,'PS',time(n),0,ps,ierr)
          if (ierr.ne.0) goto 921
        endif
        call getdat(cdfid1,'T',time(n),0,tt,ierr)
        if (ierr.ne.0) goto 920
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
     >      (mode.eq.5).or.(mode.eq.6).or.(mode.eq.1).or.
     >      (mode.eq.9).or.(mode.eq.10)) then
          if (qmode.eq.1) then
            call getdat(cdfid1,'Q',time(n),0,qq,ierr)
          else
            call getdat(cdfid1,'QD',time(n),0,qq,ierr)
          endif
          if (ierr.ne.0) goto 922
        endif
        if (mode.lt.9 .or. mode.eq.10) then
          call getdat(cdfid1,'U',time(n),0,uu,ierr)
          if (ierr.ne.0) goto 923
          call getdat(cdfid1,'V',time(n),0,vv,ierr)
          if (ierr.ne.0) goto 924
          if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.8)) then
            call getdat(cdfid1,'OMEGA',time(n),0,ww,ierr)
            if (ierr.ne.0) goto 925
          endif
        endif

        if (mode.eq.10) then 

C         Calculation of the geopotential
 
C         first get the orography
          call getoro(oro,dx,dy,stdate(1),
     >                varmin(1),varmin(2),nx,ny,ierr)
          if (ierr.eq.2) then
            stop '*** error in subrountine getoro ***'
          endif
 
          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)
 
          if (zdef.eq.0) then
            if (n.eq.1) then
              call putdef(cdfid1,'Z',4,mdv,
     >                    vardim,varmin,varmax,stag,ierr)
              write(*,*)
              write(*,*)'*** variable Z created on P-file'
            endif
          endif

          call putdat(cdfid1,'Z',time(n),0,var,ierr)
c          goto 900
        endif

C       Create the secondary data file

        if (n.eq.1) then
          cdf_file = 'S'//trim(cdfnam)
          call crecdf(cdf_file,cdfid,varmin,varmax,
     >                3,trim(cstnam),ierr)
          if (ierr.ne.0) goto 996
          write(*,*)
          write(*,*)'*** NetCDF file S',trim(cdfnam),
     >              ' created'
        endif

C       Put surface pressure on S-file. 

        if (.not.prelev) then
          vardim(3)=1
          if (n.eq.1) then
            call putdef(cdfid,'PS',4,mdv,vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable PS created on S-file'
          endif
          call putdat(cdfid,'PS',time(n),0,sp,ierr)
          vardim(3)=nz
        else
          vardim(3)=1
          if (n.eq.1) then
            call putdef(cdfid,'PS',4,mdv,vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable PS created on S-file'
          endif
          call putdat(cdfid,'PS',time(n),0,ps,ierr)
          vardim(3)=nz
        endif

C       Calculate the secondary data variables

C       Calculation of potential temperature

        if (mode.lt.9 .or. mode.eq.10) then
          call pottemp(th,tt,sp,nx,ny,nz,ak,bk)
  
          if (n.eq.1) then
            call putdef(cdfid,'TH',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable TH created on S-file'
          endif

          call putdat(cdfid,'TH',time(n),0,th,ierr)
        endif

C       Calculation of relative humidity

        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4).or.
     >      (mode.eq.5).or.(mode.eq.6).or.(mode.eq.1).or.
     >      (mode.eq.10)) then

          call relhum(rh,qq,tt,sp,nx,ny,nz,ak,bk)
  
          if (n.eq.1) then
            call putdef(cdfid,'RH',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable RH created on S-file'
          endif
 
          call putdat(cdfid,'RH',time(n),0,rh,ierr)
        endif

C       Calculation of relative vorticity

        if (mode.eq.999) then
          call relvort(var,uu,vv,sp,cl,f,nx,ny,nz,ak,bk,
     >                 varmin,varmax)

          if (n.eq.1) then
            call putdef(cdfid,'VO',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable VO created on S-file'
          endif

          call putdat(cdfid,'VO',time(n),0,var,ierr)
        endif

C       Calculation of potential vorticity
 
        if (mode.eq.5) then
          call potvort(pv,uu,vv,th,sp,cl,f,nx,ny,nz,ak,bk,
     >               varmin,varmax)
 
          if (n.eq.1) then
            call putdef(cdfid,'PV',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable PV created on S-file'
          endif
 
          call putdat(cdfid,'PV',time(n),0,pv,ierr)
        endif

C       Calculation of equivalent potential temperature

        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.1)
     >        .or.(mode.eq.10)) then
          call equpot(var,tt,qq,sp,nx,ny,nz,ak,bk)

          if (n.eq.1) then
            call putdef(cdfid,'THE',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable THE created on S-file'
          endif
           call putdat(cdfid,'THE',time(n),0,var,ierr)
        endif

C       Calculation of wet-bulb potential temperature
 
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.4)
     >       .or.(mode.eq.1).or.(mode.eq.5)) then
          call equpot(the,tt,qq,sp,nx,ny,nz,ak,bk)
          call wetbpt(var,the,sp,nx,ny,nz,ak,bk)
 
          if (n.eq.1) then
            call putdef(cdfid,'THW',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable THW created on S-file'
          endif
 
          call putdat(cdfid,'THW',time(n),0,var,ierr)
        endif

C       Calculation of the vertical gradient of the wet-bulb potential
C       temperature

        if ((mode.eq.4).or.(mode.eq.5)) then
c          if (mode.eq.5) then
c must calculate thw
c          call equpot(the,tt,qq,sp,nx,ny,nz,ak,bk)
c          call wetbpt(thw,the,sp,nx,ny,nz,ak,bk)
c          endif
        
          call dwetbptdp(var,thw,sp,nx,ny,nz,ak,bk)

          if (n.eq.1) then
            call putdef(cdfid,'DTHWDP',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable DTHWDP created on S-file'
          endif

          call putdat(cdfid,'DTHWDP',time(n),0,var,ierr)
        endif

C       Calculation of the vertical gradient of potential temperature

        if ((mode.eq.999).or.(mode.eq.5)) then
          call dpottdp(var,th,sp,nx,ny,nz,ak,bk)

          if (n.eq.1) then
            call putdef(cdfid,'DTHDP',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable DTHDP created on S-file'
          endif

          call putdat(cdfid,'DTHDP',time(n),0,var,ierr)
        endif

C       Calculation of the Richardson number

        if (mode.eq.4) then
          call richardson(var,uu,vv,th,sp,nx,ny,nz,ak,bk,mdv)

          if (n.eq.1) then
            call putdef(cdfid,'RI',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable RI created on S-file'
          endif

          call putdat(cdfid,'RI',time(n),0,var,ierr)
        endif

C       Calculation of diabatic heating rate

        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.8)) then
          call diabheat(dhr,th,ww,rh,sp,nx,ny,nz,ak,bk)

          if (n.eq.1) then
            call putdef(cdfid,'CH',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable CH created on S-file'
          endif
 
          call putdat(cdfid,'CH',time(n),0,dhr,ierr)
        endif

C       Calculation of diabatic PV rate
 
        if ((mode.eq.0).or.(mode.eq.2).or.(mode.eq.7)) then
          call diabpvr(var,uu,vv,dhr,sp,cl,f,nx,ny,nz,ak,bk,
     >                 varmin,varmax)
 
          if (n.eq.1) then
            call putdef(cdfid,'PVR',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable PVR created on S-file'
          endif
 
          call putdat(cdfid,'PVR',time(n),0,var,ierr)
        endif

C       Calculation of the geopotential

        if (mode.eq.9) then

C         first get the orography
          call getoro(oro,dx,dy,stdate(1),
     >                varmin(1),varmin(2),nx,ny,ierr)
          if (ierr.eq.2) then
            stop '*** error in subrountine getoro ***'
          endif

          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)

          if (n.eq.1) then
            call putdef(cdfid,'Z',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable Z created on S-file'
          endif
 
          call putdat(cdfid,'Z',time(n),0,var,ierr)
        endif

C       Calculation of the theta gradient
 
        if (mode.eq.11) then
 
          call pottemp(th,tt,sp,nx,ny,nz,ak,bk)
          call gradth(var,th,sp,prelev,cl,nx,ny,nz,ak,bk,varmin,varmax)
 
          if (n.eq.1) then
            call putdef(cdfid,'GRADTH',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable GRADTH created on S-file'
          endif

          call putdat(cdfid,'GRADTH',time(n),0,var,ierr)
        endif
 
C     Calculate Real W-E and N-S winds

      if (mode.eq.4) then
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
      do k=1,nz
        do j=1,ny
          lat_eq=varmin(2)+(j-1)*dy
          do i=1,nx
            lon_eq=varmin(1)+(i-1)*dx
            kk=i+(j-1)*nx+(k-1)*nx*ny
            u_ll=uvtougm(uu(kk),vv(kk),lon_eq,lat_eq,pollon,pollat)
            v_ll=uvtovgm(uu(kk),vv(kk),lon_eq,lat_eq,pollon,pollat)
            uu(kk)=u_ll
            vv(kk)=v_ll
          enddo
        enddo
      enddo
          if (n.eq.1) then
           call putdef(cdfid,'UREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable UREAL created on S-file'
          endif
 
          call putdat(cdfid,'UREAL',time(n),0,uu,ierr)

          if (n.eq.1) then
           call putdef(cdfid,'VREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable VREAL created on S-file'
          endif
 
          call putdat(cdfid,'VREAL',time(n),0,vv,ierr)
      endif
      endif

CSue  Calculate f*geostrophic winds

      if (mode.eq.999) then

C         first get the orography
          call getoro(oro,dx,dy,stdate(1),
     >                varmin(1),varmin(2),nx,ny,ierr)
          if (ierr.eq.2) then
            stop '*** error in subroutine getoro ***'
          endif

          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)


C         Then calculate the horizontal derivatives 

      if (prelev) then
        call ddh2(var,dvardx,cl,'X',nx,ny,1,varmin,varmax)
        call ddh2(var,dvardy,cl,'Y',nx,ny,1,varmin,varmax)
      else
        call ddh2(sp,dspdx,cl,'X',nx,ny,1,varmin,varmax)
        call ddh2(sp,dspdy,cl,'Y',nx,ny,1,varmin,varmax)
        call ddh3(var,dvardx,sp,dspdx,cl,'X',nx,ny,nz,
     >       varmin,varmax,ak,bk)
        call ddh3(var,dvardy,sp,dspdy,cl,'Y',nx,ny,nz,
     >       varmin,varmax,ak,bk)
      endif


C          Finally calculate Real W-E and N-S geostrophic winds

      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
      do k=1,nz
        do j=1,ny
          lat_eq=varmin(2)+(j-1)*dy
          do i=1,nx
            lon_eq=varmin(1)+(i-1)*dx
            kk=i+(j-1)*nx+(k-1)*nx*ny
            fug(kk) = -9.806665*dvardy(kk)
            fvg(kk) = 9.806665*dvardx(kk)
            u_ll=uvtougm(fug(kk),fvg(kk),
     >           lon_eq,lat_eq,pollon,pollat)
            v_ll=uvtovgm(fug(kk),fvg(kk),
     >           lon_eq,lat_eq,pollon,pollat)
            fug(kk)=u_ll
            fvg(kk)=v_ll
          enddo
        enddo
      enddo
      else
      do k=1,nz
        do j=1,ny
          do i=1,nx          
            kk=i+(j-1)*nx+(k-1)*nx*ny
            fug(kk) = -9.806665*dvardy(kk)
            fvg(kk) = 9.806665*dvardx(kk)
          enddo
        enddo
       enddo
       endif
c       write(21) fug
c       write(22) fvg
c       write(23) f

          if (n.eq.1) then
           call putdef(cdfid,'FUGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable FUGREAL created on S-file'
          endif
 
          call putdat(cdfid,'FUGREAL',time(n),0,fug,ierr)

          if (n.eq.1) then
           call putdef(cdfid,'FVGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable FVGREAL created on S-file'
          endif
 
          call putdat(cdfid,'FVGREAL',time(n),0,fvg,ierr)

      endif ! if mode
 
CSue  Calculate smoothed geostrophic winds

      if (mode.eq.4) then

C         first get the orography
          call getoro(oro,dx,dy,stdate(1),
     >                varmin(1),varmin(2),nx,ny,ierr)
          if (ierr.eq.2) then
            stop '*** error in subroutine getoro ***'
          endif

          call geopot(var,qq,tt,oro,sp,nx,ny,nz,ak,bk)


C         Then calculate the horizontal derivatives 

      if (prelev) then
        call ddh2(var,dvardx,cl,'X',nx,ny,1,varmin,varmax)
        call ddh2(var,dvardy,cl,'Y',nx,ny,1,varmin,varmax)
      else
        call ddh2(sp,dspdx,cl,'X',nx,ny,1,varmin,varmax)
        call ddh2(sp,dspdy,cl,'Y',nx,ny,1,varmin,varmax)
        call ddh3(var,dvardx,sp,dspdx,cl,'X',nx,ny,nz,
     >       varmin,varmax,ak,bk)
        call ddh3(var,dvardy,sp,dspdy,cl,'Y',nx,ny,nz,
     >       varmin,varmax,ak,bk)
      endif


c Calculated model grid orientated ug and vg
      do k=1,nz
        do j=1,ny
          do i=1,nx          
            kk=i+(j-1)*nx+(k-1)*nx*ny
            jj=i+(j-1)*nx
            ug(kk) = -9.806665*dvardy(kk)/f(jj)
            vg(kk) = 9.806665*dvardx(kk)/f(jj)
          enddo
        enddo
       enddo


c Smooth winds
c      data weight/1, 2, 1, 2, 3, 2, 1, 2, 1/ 
      data weight/49*1/
      do k=1,nz
        do j=1,ny
          do i=1,nx
            kk=i+(j-1)*nx+(k-1)*nx*ny             
            sumweight = 0
            do jw = 1, w
               do iw = 1, w
                  iiw = i+iw-(int(w/2.)+1)
                  jjw = j+jw-(int(w/2.)+1)
                  if (iiw.lt.1) iiw = 1
                  if (jjw.lt.1) jjw = 1 

c            if (iiw.gt.nx) then 
c               write(*,*) 'nx ',nx,' ny ',ny,' i ',i,' j ',j
c               write(*,*) 'jw ',jw,' iw ',iw
c               write(*,*) 'ug ', ug(iiw+(jjw-1)*nx+(k-1)*nx*ny),
c     >                    'vg ', vg(iiw+(jjw-1)*nx+(k-1)*nx*ny)
c            do tmpjw = 1, w
c               do tmpiw = 1, w
c                  tmpiiw = i+tmpiw-(int(w/2.)+1)
c                  tmpjjw = j+tmpjw-(int(w/2.)+1)
c                  if (tmpiiw.lt.1) tmpiiw = 1
c                  if (tmpjjw.lt.1) tmpjjw = 1                
c                  if (tmpiiw.gt.nx) tmpiiw = nx
c                  if (tmpjjw.gt.ny) tmpjjw = ny
c                  kkw = tmpiiw+(tmpjjw-1)*nx+(k-1)*nx*ny   
c                  write(*,*) 'tmpiiw ', tmpiiw,' tmpjjw ', tmpjjw,
c     >                       ' vg ',vg(kkw),' vg ',vg(kkw)
c                enddo
c            enddo
c            stop
c            endif

                  if (iiw.gt.nx) iiw = nx
                  if (jjw.gt.ny) jjw = ny                  
                  kkw = iiw+(jjw-1)*nx+(k-1)*nx*ny   
                  tmpu(kk) = tmpu(kk) + ug(kkw)*weight(iw,jw)
                  tmpv(kk) = tmpv(kk) + vg(kkw)*weight(iw,jw) 
                  sumweight = sumweight+weight(iw,jw)
                enddo
            enddo
            tmpu(kk) = tmpu(kk)/sumweight
            tmpv(kk) = tmpv(kk)/sumweight
         enddo
        enddo
      enddo

      do k=1,nz
        do j=1,ny
          do i=1,nx
            kk=i+(j-1)*nx+(k-1)*nx*ny
            ug(kk) = tmpu(kk)
            vg(kk) = tmpv(kk)
         enddo
        enddo
      enddo

C          Finally calculate Real W-E and N-S geostrophic winds

      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
      do k=1,nz
        do j=1,ny
          lat_eq=varmin(2)+(j-1)*dy
          do i=1,nx
            lon_eq=varmin(1)+(i-1)*dx
            kk=i+(j-1)*nx+(k-1)*nx*ny
            jj=i+(j-1)*nx
            u_ll=uvtougm(ug(kk),vg(kk),
     >           lon_eq,lat_eq,pollon,pollat)
            v_ll=uvtovgm(ug(kk),vg(kk),
     >           lon_eq,lat_eq,pollon,pollat)
            ug(kk)=u_ll
            vg(kk)=v_ll
          enddo
        enddo
      enddo
      endif
c       write(21) fug
c       write(22) fvg
c       write(23) f
        
c      open(111,file='ug.dat',status = 'unknown',form = 'unformatted')
c      write(111) ug      
c      open(112,file='vg.dat',status = 'unknown',form = 'unformatted')
c      write(112) vg
c      open(113,file='geo.dat',status = 'unknown',form = 'unformatted')
c      write(113) var
c      stop

          if (n.eq.1) then
           call putdef(cdfid,'UGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable UGREAL created on S-file'
          endif
 
          call putdat(cdfid,'UGREAL',time(n),0,ug,ierr)

          if (n.eq.1) then
           call putdef(cdfid,'VGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable VGREAL created on S-file'
          endif
 
          call putdat(cdfid,'VGREAL',time(n),0,vg,ierr)


cSue Calculate geostropic momentum components

      do k=1,nz
        do j=1,ny
          do i=1,nx    
            jj=i+(j-1)*nx
            kk=i+(j-1)*nx+(k-1)*nx*ny
            Mg(kk) = 0.000145444*6.378e6*sin(lat(jj))*
     >           (lon(jj) - lon(1))*cos(lat(jj)) 
     >           + vg(kk)
            Ng(kk) = 0.000145444*6.378e6*
     >           (cos(lat(1)) - cos(lat(jj))) - 
     >           ug(kk)
          enddo
        enddo
      enddo
 
          if (n.eq.1) then
           call putdef(cdfid,'MGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable MGREAL created on S-file'
          endif
 
          call putdat(cdfid,'MGREAL',time(n),0,Mg,ierr)

          if (n.eq.1) then
           call putdef(cdfid,'NGREAL',4,mdv,
     >                  vardim,varmin,varmax,stag,ierr)
            write(*,*)
            write(*,*)'*** variable NGREAL created on S-file'
          endif
 
          call putdat(cdfid,'NGREAL',time(n),0,Ng,ierr)

       endif

      enddo !loop over n


C     Close the NetCDF files

      call clscdf(cdfid,ierr)
 900  continue
      call clscdf(cdfid1,ierr)
      call clscdf(cstid,ierr)

      goto 999

  920 stop '*** error: variable T not found on P-file ***'
  921 stop '*** error: variable PS not found on P-file ***'
  922 stop '*** error: variable Q not found on P-file ***'
  923 stop '*** error: variable U not found on P-file ***'
  924 stop '*** error: variable V not found on P-file ***'
  925 stop '*** error: variable OMEGA not found on P-file ***'

  996 stop '*** error: could not create S-file ***'
  999 continue
      end

      subroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)
c     ==========================================
 
c     argument declaration
      integer   ie,je,ke
      real      pt(ie,je,ke),t(ie,je,ke),sp(ie,je),
     >          ak(ke),bk(ke)
 
c     variable declaration
      integer   i,j,k
      real      rdcp,tzero
      data      rdcp,tzero /0.286,273.15/
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
c     computation of potential temperature
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
c         distinction of temperature in K and deg C
          if (t(i,j,k).lt.100.) then
            pt(i,j,k)=(t(i,j,k)+tzero)*( (1000./pp(k))**rdcp )
          else
            pt(i,j,k)=t(i,j,k)*( (1000./pp(k))**rdcp )
          endif
        enddo
      enddo
      enddo
      end

      subroutine gradth(gth,th,sp,prelev,cl,ie,je,ke,ak,bk,vmin,vmax)
C     ===============================================================
 
c     argument declaration
      integer   ie,je,ke
      real      gth(ie,je,ke),th(ie,je,ke),sp(ie,je),cl(ie,je)
      real      ak(ke),bk(ke),vmin(4),vmax(4)
      logical   prelev
 
c     variable declaration
      include "um_dims.inc"
      real      dthdx(nxmax*nymax*nzmax),dthdy(nxmax*nymax*nzmax)
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
      integer   i,j,k,ind,ind2
 
      if (prelev) then
        call ddh2(th,dthdx,cl,'X',ie,je,1,vmin,vmax)
        call ddh2(th,dthdy,cl,'Y',ie,je,1,vmin,vmax)
      else
        call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
        call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
        call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
        call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
      endif
 
      do j=1,je
      do i=1,ie
        ind2=i+(j-1)*ie
        do k=1,ke
          ind=ind2+(k-1)*ie*je
          gth(i,j,k)=sqrt(dthdx(ind)**2.+dthdy(ind)**2.)
        enddo
      enddo
      enddo
      end

      subroutine geopot(psi,q,t,oro,sp,ie,je,ke,ak,bk)
c     ================================================
 
c     argument declaration
      integer   ie,je,ke
      real      psi(ie,je,ke),t(ie,je,ke),q(ie,je,ke),oro(ie,je),
     >          sp(ie,je),ak(ke),bk(ke)
 
c     variable declaration
      integer   i,j,k
      real      r,c,g
      data      r,c,g /287.,0.608,9.8/
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
c     integration of geopotential height(special for first layer)
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        psi(i,j,1)=oro(i,j)+r/g*
     >                 ((t(i,j,1)+273.15)*(1.+c*q(i,j,1))*
     >                (psrf-pp(1))/(0.5*(psrf+pp(1))))
c        psi(i,j,1)=1./g*(oro(i,j)
c     >                 +r*(t(i,j,1)+273.15)*(1.+c*q(i,j,1))*
c     >                (psrf-pp(1))/(0.5*(psrf+pp(1))))
      enddo
      enddo
      do j=1,je
      do i=1,ie
        psrf=sp(i,j)
        do k=2,ke
          psi(i,j,k)=psi(i,j,k-1)+r/g*
     >                  ((t(i,j,k-1)+273.15)*(1.+c*q(i,j,k-1))+
     >                   (t(i,j,k)+273.15)*(1.+c*q(i,j,k)))*
     >                  (pp(k-1)-pp(k))/(pp(k-1)+pp(k))
        enddo
      enddo
      enddo
      end

      subroutine getoro(oro,dx,dy,starty,lonmin,latmin,ie,je,ierr)
c     ===========================================================
c     reads the orography for the actual data domain from file
 
      include 'netcdf.inc'
      include "um_dims.inc"
 
c     argument declaration
      integer  ie,je
      real     or(nxmax*nymax)
      real     oro(ie,je), fld(nxmax,nymax)
 
c     variable declaration
      integer look(45)
      real rook(19)
      character*52 name
      character*52 filename
      integer  starty
 
c     open file with orography values and get them
c      write(*,*) 'Inset name of orography file'
c      read(*,*) name
c      filename = '/export/wombat/wombat-01/sws98slg/'//name
      write(*,*) ' Using orography file MESO17_orog_high_res.pp ' 
      open(10,file = 
     >'/export/cloud/stingjet/rb904381/umdata/orog_NAE.pp32',
     >  form = 'unformatted', status = 'old')
cNH     >'/export/wombat/wombat-01/sws98slg/MESO17_orog_high_res.pp',
cNH     >  form = 'unformatted', status = 'old')
c      open(10,file = '/export/wombat/wombat-01/sws98slg/'//name,
c     >  form = 'unformatted', status = 'old')
      read(10) look, rook
      nx = look(19)
      ny = look(18)
      read(10) (or(i), i = 1, nx*ny)
      close(10)

c data written starting at N-W corner, following loops assign data from S-W 
c corner and write out.
      do j = ny-1,0,-1
         i = j*nx
         do n = 1, nx
            fld(n,ny-j) = or(n+i)
         enddo
      enddo
        do i=1,ie
          do j=1,je
            oro(i,j)=fld(i+1,j+1)
          enddo
        enddo
      
 
      end

      subroutine wetbpt(thw,the,sp,ie,je,ke,ak,bk)
c     ============================================

c     argument declaration
      integer   ie,je,ke
      real      thw(ie,je,ke),the(ie,je,ke),
     >          sp(ie,je),ak(ke),bk(ke)

c     variable declaration
      real      tsa
      integer   i,j,k

      do k=1,ke
      do j=1,je
      do i=1,ie
        thw(i,j,k)=tsa(the(i,j,k)-273.15,1000.)+273.15
      enddo
      enddo
      enddo
      end

      subroutine dwetbptdp(dthwdp,thw,sp,ie,je,ke,ak,bk)
c     ==================================================

c     argument declaration
      integer   ie,je,ke
      real      dthwdp(ie,je,ke),thw(ie,je,ke),
     >          sp(ie,je),ak(ke),bk(ke)

      call ddp(thw,dthwdp,sp,ie,je,ke,ak,bk)

      end

      subroutine dpottdp(dthdp,th,sp,ie,je,ke,ak,bk)
c     ==============================================

c     argument declaration
      integer   ie,je,ke
      real      dthdp(ie,je,ke),th(ie,je,ke),
     >          sp(ie,je),ak(ke),bk(ke)

      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)

      end

      subroutine richardson(ri,uu,vv,th,sp,ie,je,ke,ak,bk,mdv)
c     ========================================================

c     argument declaration
      integer   ie,je,ke
      real      ri(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
     >          th(ie,je,ke),sp(ie,je),ak(ke),bk(ke),mdv
c     variable declaration
      include "um_dims.inc"
      real      vel(nxmax,nymax,nzmax),dthdp(nxmax,nymax,nzmax),
     &          dveldp(nxmax,nymax,nzmax)

c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf

      do k=1,ke
      do j=1,je
      do i=1,ie
        vel(i,j,k)=(uu(i,j,k)**2.+vv(i,j,k)**2.)**0.5
      enddo
      enddo
      enddo

      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
      call ddp(vel,dveldp,sp,ie,je,ke,ak,bk)

c     computation of the richardson number
c     use: kappa-1 = -0.714
c          p0**(-kappa) = 0.0372
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
          if (abs(dveldp(i,j,k)).lt.0.0001) then
            ri(i,j,k)=1000.
          else
            ri(i,j,k)=-287.*((100.*pp(k))**(-.714))*.0372*
     >                dthdp(i,j,k)/(dveldp(i,j,k)**2.)
          endif
          if (ri(i,j,k).gt.1000.) ri(i,j,k)=1000.
        enddo
      enddo
      enddo

      end

      real function tsa(os,p)
c     =======================
 
C     This function returns the temperature tsa (celsius) on a saturation
C     adiabat at pressure p (millibars). os is the equivalent potential
C     temperature of the parcel (celsius). sign(a,b) replaces the
C     algebraic sign of a with that of b.
C     b is an empirical constant approximately equal to 0.001 of the latent
C     heat of vaporization for water divided by the specific heat at constant
C     pressure for dry air.
 
      real      a,b,os,p,tq,d,tqk,x,w
      integer   i
 
      data b/2.6518986/
      a=os+273.16
 
C     tq is the first guess for tsa
 
      tq=253.16
 
C     d is an initial value used in the iteration below
 
      d=120.
 
C     Iterate to obtain sufficient accuracy....see table 1, p.8
C     of Stipanuk (1973) for equation used in iteration
      do 1 i=1,12
        tqk=tq-273.16
        d=d/2.
        x=a*exp(-b*w(tqk,p)/tq)-tq*((1000./p)**.286)
        if (abs(x).lt.1e-7) go to 2
        tq=tq+sign(d,x)
 1    continue
 2    tsa=tq-273.16
      end

      real function w(t,p)
c     ====================
 
C     This function returns the mixing ratio (grams of water vapor per
C     kilogram of dry air) given the dew point (celsius) and pressure
C     (millibars). If the temperature is input instead of the
C     dew point, then saturation mixing ratio (same units) is returned.
C     The formula is found in most meteorological texts.
 
      real      t,p,tkel,x,esat
 
      tkel=t+273.16
      x=esat(tkel)      ! our function esat requires t in Kelvin
      w=622.*x/(p-x)
      end

      real function esat(t)
c     =====================
 
C     This function returns the saturation vapor pressure over water (mb)
C     given the temperature (Kelvin).
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
C     Approximations of Selected Meteorological Parameters for Cloud
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory, U. S.
C     Army Electronics Command, White Sands Missile Range, New Mexico 88002.
 
      real p1,p2,c1,t
 
      p1=11.344-0.0303998*t
      p2=3.49149-1302.8844/t
      c1=23.832241-5.02808*log10(t)
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
      end
 
      subroutine equpot(ap,t,q,sp,ie,je,ke,ak,bk)
c     ===========================================
 
c     argument declaration
      integer  ie,je,ke
      real     ap(ie,je,ke),t(ie,je,ke),sp(ie,je)
      real     q(ie,je,ke),ak(ke),bk(ke)
 
c     variable declaration
      integer  i,j,k
      real     rdcp,tzero
      data     rdcp,tzero /0.286,273.15/
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
c     computation of potential temperature
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
          ap(i,j,k) = (t(i,j,k)+tzero)*(1000./pp(k))
     +       **(0.2854*(1.0-0.28*q(i,j,k)))*exp(
     +       (3.376/(2840.0/(3.5*alog(t(i,j,k)+tzero)-alog(
     +       100.*pp(k)*max(1.0E-10,q(i,j,k))/(0.622+0.378*
     +       q(i,j,k)))-0.1998)+55.0)-0.00254)*1.0E3*
     +       max(1.0E-10,q(i,j,k))*(1.0+0.81*q(i,j,k)))
        enddo
      enddo
      enddo
      end

      subroutine relhum(rh,q,t,sp,ie,je,ke,ak,bk)
c     ===========================================
 
c     argument declaration
      integer   ie,je,ke
      real      rh(ie,je,ke),t(ie,je,ke),q(ie,je,ke),
     >          sp(ie,je),ak(ke),bk(ke)
 
c     variable declaration
      integer   i,j,k
      real      rdcp,tzero
      real      b1,b2w,b3,b4w,r,rd,gqd,ge
      data      rdcp,tzero /0.286,273.15/
      data      b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,
     &          287.05, 461.51/
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
          ge = b1*exp(b2w*(t(i,j,k))/(t(i,j,k)+b3-b4w))
          gqd= r/rd*ge/(pp(k)-(1.-r/rd)*ge)
          rh(i,j,k)=100.*q(i,j,k)/gqd
        enddo
      enddo
      enddo
      end

      subroutine diabheat(dhr,t,w,rh,sp,ie,je,ke,ak,bk)
c     =================================================
 
c     argument declaration
      integer   ie,je,ke
      real      dhr(ie,je,ke),t(ie,je,ke),w(ie,je,ke),
     &          rh(ie,je,ke),sp(ie,je),ak(ke),bk(ke)
 
c     variable declaration
      integer   i,j,k
      real      p0,kappa,tzero
      data      p0,kappa,tzero /1000.,0.286,273.15/
      real      blog10,cp,r,lw,eps
      data      blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/
      real      esat,c,tt
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
c     computation of diabatic heating rate
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
          if (rh(i,j,k).lt.80.) then    ! only moist air of interest
            dhr(i,j,k)=0.               ! cond. heating rate set to zero
          else if (w(i,j,k).gt.0.) then ! cond. heating only for ascent
            dhr(i,j,k)=0.
          else
            tt=t(i,j,k)*((pp(k)/p0)**kappa)   ! temp. from pot.temp.
            c=lw/cp*eps*blog10*esat(tt)/pp(k)
            dhr(i,j,k)=21600.*                ! in units K per 6 hours
     >        (1.-exp(.2*(80.-rh(i,j,k))))    ! weighting fun. for 80<RH<100
     >        *(-c*kappa*t(i,j,k)*w(i,j,k)/(100.*pp(k)))/(1.+c)
          endif
        enddo
      enddo
      enddo
      end

      subroutine diabpvr(dpvr,uu,vv,dhr,sp,cl,f,ie,je,ke,ak,bk,
     >           vmin,vmax)
C     =========================================================
 
c     argument declaration
      integer   ie,je,ke
      real      dpvr(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
     >          dhr(ie,je,ke),sp(ie,je),cl(ie,je),f(ie,je)
      real      ak(ke),bk(ke),vmin(4),vmax(4)
 
c     variable declaration
      include "um_dims.inc"
      real      dhrdp(nxmax*nymax*nzmax),dudp(nxmax*nymax*nzmax),
     >          dvdp(nxmax*nymax*nzmax),dvdx(nxmax*nymax*nzmax),
     >          dhrdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax),
     >          dhrdy(nxmax*nymax*nzmax)
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
      integer   i,j,k,ind,ind2
 
      call ddp(dhr,dhrdp,sp,ie,je,ke,ak,bk)
      call ddp(uu,dudp,sp,ie,je,ke,ak,bk)
      call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)
      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
      call ddh3(dhr,dhrdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(dhr,dhrdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
 
      do j=1,je
      do i=1,ie
        ind2=i+(j-1)*ie
        do k=1,ke
          ind=ind2+(k-1)*ie*je
          dpvr(i,j,k)=-1.E6*9.80665*(
     >              -dvdp(ind)*dhrdx(ind)+dudp(ind)*dhrdy(ind)
     >              +(-dudy(ind)+dvdx(ind)+f(i,j))*dhrdp(ind))
        enddo
      enddo
      enddo
      end

      subroutine potvort(pv,uu,vv,th,sp,cl,f,ie,je,ke,ak,bk,
     >           vmin,vmax)
C     ======================================================

c     argument declaration
      integer   ie,je,ke
      real      pv(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
     >          th(ie,je,ke),sp(ie,je),
     >          cl(ie,je),f(ie,je)
      real      ak(ke),bk(ke),vmin(4),vmax(4)
 
c     variable declaration
      include "um_dims.inc"
      real      dthdp(nxmax*nymax*nzmax),dudp(nxmax*nymax*nzmax),
     >          dvdp(nxmax*nymax*nzmax),dvdx(nxmax*nymax*nzmax),
     >          dthdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax),
     >          dthdy(nxmax*nymax*nzmax)
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
      integer   i,j,k,ind,ind2

      call ddp(th,dthdp,sp,ie,je,ke,ak,bk)
      call ddp(uu,dudp,sp,ie,je,ke,ak,bk)
      call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)
      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
      call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)

      do j=1,je
      do i=1,ie
        ind2=i+(j-1)*ie
        do k=1,ke
          ind=ind2+(k-1)*ie*je
          pv(i,j,k)=1.E6*9.80665*(
     >              -(-dudy(ind)+dvdx(ind)+f(i,j))*dthdp(ind)
     >              -(dudp(ind)*dthdy(ind)-dvdp(ind)*dthdx(ind)))
        enddo
      enddo
      enddo
      end

      subroutine relvort(vo,uu,vv,sp,cl,f,ie,je,ke,ak,bk,
     >           vmin,vmax)
C     ===================================================

c     argument declaration
      integer   ie,je,ke
      real      vo(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),
     >          sp(ie,je),cl(ie,je),f(ie,je)
      real      ak(ke),bk(ke),vmin(4),vmax(4)

c     variable declaration
      include "um_dims.inc"
      real      dvdx(nxmax*nymax*nzmax),dudy(nxmax*nymax*nzmax)
      real      dspdx(nxmax*nymax),dspdy(nxmax*nymax)
      integer   i,j,k,ind,ind2

      call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)
      call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)
      call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)
      call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)

      do j=1,je
      do i=1,ie
        ind2=i+(j-1)*ie
        do k=1,ke
          ind=ind2+(k-1)*ie*je
          vo(i,j,k)=1.E4*(-dudy(ind)+dvdx(ind))
        enddo
      enddo
      enddo
      end


      subroutine ddp(a,d,sp,ie,je,ke,ak,bk)
c--------------------------------------------------------------------------
c     Purpose: VERTICAL DERIVATIVE
c       Compute the vertical derivative without missing data checking.
c       The derivative is taken from array 'a' in the direction of 'P'.
c       The result is stored in array 'd'.
c       3 point weighted centered differencing is used.
c       The vertical level-structure of the data is of the form
c       p=ak(k)+bk(k)*ps.
c--------------------------------------------------------------------------
 
c     declaration of arguments
      integer       ie,je,ke
      real          a(ie,je,ke),d(ie,je,ke),sp(ie,je)
 
c     variable declaration
      integer       i,j,k
      real          dpu,dpl,quot,fac,psrf
      real          ak(ke),bk(ke)
 
c     specify scaling factor associated with derivation with respect
c     to pressure
      fac=0.01
 
c     compute vertical 3 point derivative
c     ---------------------------
c     3-point vertical derivative
c     ---------------------------
      do j=1,je
      do i=1,ie
c       get surface pressure at current grid-point
        psrf=sp(i,j)
c       points at k=1
        dpu=(ak(1)+bk(1)*psrf)-(ak(2)+bk(2)*psrf)
        d(i,j,1)=(a(i,j,1)-a(i,j,2))*fac/dpu
c       points at 1<k<ke
        do k=2,ke-1
          dpu=(ak(k)+bk(k)*psrf)-(ak(k+1)+bk(k+1)*psrf)
          dpl=(ak(k-1)+bk(k-1)*psrf)-(ak(k)+bk(k)*psrf)
          quot=dpu/dpl
          d(i,j,k)=(quot*(a(i,j,k-1)-a(i,j,k))
     &              +1./quot*(a(i,j,k)-a(i,j,k+1)))*fac/(dpu+dpl)
        enddo
c       points at k=ke
        dpl=(ak(ke-1)+bk(ke-1)*psrf)-(ak(ke)+bk(ke)*psrf)
        d(i,j,ke)=(a(i,j,ke-1)-a(i,j,ke))*fac/dpl
      enddo
      enddo

      end

      subroutine ddh3(a,d,ps,dps,cl,dir,ie,je,ke,datmin,datmax,ak,bk)
c-----------------------------------------------------------------------
c     Purpose: HORIZONTAL DERIVATIVE ON PRESSURE-SURFACES WITHOUT MISSING DATA
c       The derivative is taken from array 'a' in the direction of 'dir',
c       where 'dir' is either 'X','Y'. The result is stored in array 'd'.
c       The routine accounts for derivatives at the pole and periodic
c       boundaries in the longitudinal direction (depending on
c       the value of datmin, datmax). If the data-set does not reach to
c       the pole, a one-sided derivative is taken. Pole-treatment is only
c       carried out if the data-set covers 360 deg in longitude, and it
c       requires that ie=4*ii+1, where ii is an integer.
c     History:
c       Daniel Luethi
c-----------------------------------------------------------------------
 
c     declaration of arguments
      integer       ie,je,ke
      real          a(ie,je,ke),d(ie,je,ke),cl(ie,je)
      real          ps(ie,je),dps(ie,je)
      real          datmin(4),datmax(4)
      character*(*) dir
 
c     variable declaration
      integer       i,j,k
      real          ak(ke),bk(ke),as(ke),bs(ke)

c     compute vertical derivatives of ak's and bk's
      do k=2,ke-1
        as(k)=(ak(k-1)-ak(k+1))/2.
        bs(k)=(bk(k-1)-bk(k+1))/2.
      enddo
      as(1 )=ak(1)-ak(2)
      bs(1 )=bk(1)-bk(2)
      as(ke)=ak(ke-1)-ak(ke)
      bs(ke)=bk(ke-1)-bk(ke)
 
c     compute horizontal derivatives on sigma surfaces
      call ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)
 
c     apply correction for horizontal derivative on p-surfaces
      do j=1,je
      do i=1,ie
        do k=2,ke-1
          d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/2./(as(k)+
     &               bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k-1))
        enddo
        k=1
        d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+
     &             bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k))
        k=ke
        d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+
     &             bs(k)*ps(i,j))*(a(i,j,k)-a(i,j,k-1))
      enddo
      enddo
      end

      subroutine ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)
c-----------------------------------------------------------------------
c     Purpose: HORIZONTAL DERIVATIVE ON DATA-SURFACES WITHOUT MISSING DATA
c       Compute the horizontal derivative without missing data checking.
c       The derivative is taken from array 'a' in the direction of 'dir',
c       where 'dir' is either 'X','Y'. The result is stored in array 'd'.
c          The routine accounts for derivatives at the pole and periodic
c       boundaries in the longitudinal direction (depending on
c       the value of datmin, datmax). If the data-set does not reach to
c       the pole, a one-sided derivative is taken. Pole-treatment is only
c       carried out if the data-set covers 360 deg in longitude, and it
c       requires that ie=4*ii+1, where ii is an integer.
c-----------------------------------------------------------------------
 
c     declaration of arguments
      integer       ie,je,ke
      real          a(ie,je,ke),d(ie,je,ke),cl(ie,je)
      real          datmin(4),datmax(4)
      character*(*) dir
 
c     local variable declaration
      integer       i,j,k,ip1,im1,jp1,jm1,ip,im,j1,j2
      real          dlat,dlon,coslat,dx,dy,dxr,dyr
      integer       northpl,southpl,lonper
 
c     rerd and circ are the mean radius and diameter of the earth in meter
      real          rerd,circ,pi
      data          rerd,circ,pi /6.37e6,4.e7,3.141592654/
 
c     compute flags for pole and periodic treatment
      southpl=0
      northpl=0
      lonper =0
      j1=1
      j2=je
      if (abs(datmax(1)-datmin(1)-360.).lt.1.e-3) then
        lonper=1
        if (abs(datmin(2)+90.).lt.1.e-3) then
          southpl=1
          j1=2
        endif
        if (abs(datmax(2)-90.).lt.1.e-3) then
          northpl=1
          j2=je-1
        endif
      endif
 
      dlon=((datmax(1)-datmin(1))/float(ie-1)) *pi/180.
      dlat=((datmax(2)-datmin(2))/float(je-1)) *pi/180.
 
c      print *,'Computing derivative ',dir(1:1),
c     &        ' of an array dimensioned ',ie,je,ke
 
      if (dir(1:1).eq.'X') then
c       -----------------------------
c       derivation in the x-direction
c       -----------------------------
        do k=1,ke
 
c         do gridpoints at j1<=j<=j2
          do j=j1,j2
           coslat=cl(1,j)
 
c          do regular gridpoints at 1<i<ie, 1<j<je
           dx =rerd*coslat*dlon
           dxr=1./(2.*dx)
           do i=2,ie-1
            ip1=i+1
            im1=i-1
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
           enddo ! i-loop
c          completed regular gridpoints at 1<i<ie, 1<j<je
 
c          do gridpoints at i=1, i=ie, 1<j<je
           if (lonper.eq.1) then
c           use periodic boundaries
            i=1
            ip1=2
            im1=ie-1
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
            d(ie,j,k)=d(1,j,k)
           else
c           use one-sided derivatives
            dxr=1./dx
            i=1
            ip1=2
            im1=1
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
            i=ie
            ip1=ie
            im1=ie-1
            d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))
           endif
c          completed gridpoints at i=1, i=ie, j1<=j<=j2
 
          enddo ! j-loop
c         completed gridpoints at 1<j<je
 
c         do gridpoints at j=je
          if (northpl.eq.1) then
c          for these gridpoints, the derivative in the x-direction is a
c          derivative in the y-direction at another pole-gridpoint
           dy =rerd*dlat
           dyr=1./(2.*dy)
           j=je
           jp1=je-1
           jm1=je-1
           do i=1,ie
            ip=mod(i-1+  (ie-1)/4,ie)+1
            im=mod(i-1+3*(ie-1)/4,ie)+1
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
           enddo    ! i-loop
c          completed gridpoints at j=je
          endif
c         do gridpoints at j=1
          if (southpl.eq.1) then
           dy =rerd*dlat
           dyr=1./(2.*dy)
           j=1
           jp1=2
           jm1=2
           do i=1,ie
            ip=mod(i-1+  (ie-1)/4,ie)+1
            im=mod(i-1+3*(ie-1)/4,ie)+1
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
           enddo    ! i-loop
          endif
c         completed gridpoints at j=1
 
        enddo    ! k-loop
 
      else if (dir(1:1).eq.'Y') then
c       -----------------------------
c       derivation in the y-direction
c       -----------------------------
        dy =dlat*rerd
        dyr=1./(2.*dy)
        do k=1,ke
         do i=1,ie
 
c          do regular gridpoints
           do j=2,je-1
            jp1=j+1
            jm1=j-1
            d(i,j,k)=dyr*(a(i,jp1,k)-a(i,jm1,k))
           enddo
 
c          do gridpoints at j=je
           if (northpl.eq.1) then
c           pole-treatment
            j=je
            jm1=j-1
            jp1=j-1
            ip=mod(i-1+(ie-1)/2,ie)+1
            im=i
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
           else
c           one-sided derivative
            j=je
            jm1=j-1
            jp1=j
            d(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))
           endif
c          completed gridpoints at j=je
 
c          do gridpoints at j=1
           if (southpl.eq.1) then
c           pole-treatment
            j=1
            jm1=2
            jp1=2
            ip=i
            im=mod(i-1+(ie-1)/2,ie)+1
            d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))
           else
c           one-sided derivative
            j=1
            jm1=1
            jp1=2
            d(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))
           endif
c          completed gridpoints at j=1
 
         enddo
        enddo
 
 
      endif
      end
c--------------------------------------------------------------------------
c     transformation of coordinates with a rotated pole
c--------------------------------------------------------------------------
c     the following routines allow for the transformation of
c     spherical coordinates to a spherical coordinate system with a rotated
c     pole, and vice versa.
c
c     the arguments of the routines use the following conventions
c        PHI,LAM:       latitude and longitude in regular spherical coordinates
c        lat,lon:       dito
c        PHIS,LAMS:     latitude and longitude in rotated spherical coordinates
c        latp,lonp:     dito
c        polphi,pollat: latitude and longitude of the rotated pole as
c                       expressed in regular spherical coordinates
c        latpol,lonpol: dito
c        u,v:           windvector in regular spherical coordinates
c        up,vp:         windvector in rotated spherical coordinates
c     all angles are expressed in degrees.
c
c     For a given rotated pole location (pollat,polphi), the routines   
c     provide the following transformations:
c         REAL FUNCTION PHTOPHS: (lam,phi) -> phis
c         REAL FUNCTION PHSTOPH: (lams,phis) -> phi
c         REAL FUNCTION LMTOLMS: (lam,phi) -> lams
c         REAL FUNCTION LMSTOLM: (lams,phis) -> lams
c         real function uvtouem: (u,v,lon,lat) -> up
c         real function uvtovem: (u,v,lon,lat) -> vp
c         real function uvtougm: (up,vp,lonp,latp) -> u
c         real function uvtovgm: (up,vp,lonp,latp) -> v
c--------------------------------------------------------------------------

      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
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
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
      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 LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
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
c-----------------------------------------------------------------
      real function uvtouem(u,v,lon,lat,lonpol,latpol)
c-----------------------------------------------------------------
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
c     INPUT:
c     - winds in geographical coordinated system
c     - geographical  coordinates
c     - coordinates of the rotated pole
c     OUTPUT
c     - u-component of the rotated wind
c-----------------------------------------------------------------

      real               u,v,lon,lat,lonpol,latpol
      real               lampol,phipol,lam,phi
      real               cosdlam,cosphipol,cosphi
      real               sindlam,sinphipol,sinphi

      data               pid180 / 0.0174532925199 /

c-----------------------------------------------------------------

      lampol=lonpol*pid180
      phipol=latpol*pid180
      lam=lon*pid180
      phi=lat*pid180
      cosdlam=cos(lampol-lam)
      cosphipol=cos(phipol)
      cosphi=cos(phi)
      sinphipol=sin(phipol)
      sinphi=sin(phi)
      cosdlam=cosdlam
      sindlam=sin(lampol-lam)

      uvtouem=-(sqrt(1.-(cosdlam*cosphipol*cosphi+ 
     >     sinphipol*sinphi)**2.)*(-(u*cosdlam/
     >     ((cosdlam*cosphi*sinphipol- 
     >     cosphipol*sinphi))) + sindlam*
     >     (-(v*sinphi/
     >     ((cosdlam*cosphi*sinphipol - 
     >     cosphipol*sinphi))) - 
     >     cosphi*(-(v*cosphipol*cosphi) + 
     >     sinphipol*(u*sindlam - v*cosdlam*sinphi))/
     >     (cos(lampol - lam)*cosphi*sinphipol - 
     >     cosphipol*sinphi)**2))/
     >     (1. + cosphi**2. *sindlam**2./
     >     (cosdlam*cosphi*sinphipol - 
     >     cosphipol*sinphi)**2.))
    
      return
      end
c-----------------------------------------------------------------
      real function uvtovem(u,v,lon,lat,lonpol,latpol)
c-----------------------------------------------------------------
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
c     INPUT:
c     - winds in geographical coordinated system
c     - geographical  coordinates
c     - coordinates of the rotated pole
c     OUTPUT
c     - v-component of the rotated wind
c-----------------------------------------------------------------

      real               u,v,lon,lat,lonpol,latpol
      real               lampol,phipol,lam,phi

      real               cosdlam,cosphipol,cosphi
      real               sindlam,sinphipol,sinphi

      data               pid180 / 0.0174532925199 /

c-----------------------------------------------------------------

      lampol=lonpol*pid180
      phipol=latpol*pid180
      lam=lon*pid180
      phi=lat*pid180
      cosdlam=cos(lampol-lam)
      cosphipol=cos(phipol)
      cosphi=cos(phi)
      sinphipol=sin(phipol)
      sinphi=sin(phi)
      cosdlam=cosdlam
      sindlam=sin(lampol-lam)

      uvtovem=(v*cosphi*sinphipol + cosphipol*
     >     (u*sindlam-v*cosdlam*sinphi))/
     >     sqrt(1.-(cosdlam*cosphipol*cosphi + 
     >     sinphipol*sinphi)**2.)

      return
      end
c-----------------------------------------------------------------
      real function uvtougm(up,vp,lonp,latp,lonpol,latpol)
c-----------------------------------------------------------------
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
c     INPUT:
c     - rotated winds (i.e. wrt equatorial grid)
c     - rotated coordinates (i.e. on equatorial grid)
c     - coordinates of the rotated pole
c     OUTPUT
c     - u-component of the wind (i.e. real W-E wind)
c-----------------------------------------------------------------

      real               up,vp,lonp,latp,lonpol,latpol
      real               lampol,phipol,lamp,phip,pid180

      data               pid180 / 0.0174532925199 /

c-----------------------------------------------------------------

      phip=latp*pid180
      lamp=lonp*pid180
      phipol=latpol*pid180
      lampol=lonpol*pid180

      uvtougm=8*Cos(ASin(Cos(lamp)*Cos(phipol)*Cos(phip) +
     -     Sin(phipol)*Sin(phip)))*
     -   (-2*vp*Sin(lamp - phipol) - 2*vp*Sin(lamp + phipol) - 
     -     up*Sin(lamp - phipol - phip) - 
     >     2*up*Sin(phipol - phip) - 
     -     up*Sin(lamp + phipol - phip) + up*Sin(lamp - 
     >     phipol + phip) - 
     -     2*up*Sin(phipol + phip) + up*Sin(lamp + phipol + phip))/
     -  (-20 + 4*Cos(2*lamp) + 2*Cos(2*(lamp - phipol)) - 
     >     4*Cos(2*phipol) + 
     -    2*Cos(2*(lamp + phipol)) - 4*Cos(lamp - 2*phipol - 2*phip) + 
     -    4*Cos(lamp + 2*phipol - 2*phip) + 2*Cos(2*(lamp - phip)) + 
     -    Cos(2*(lamp - phipol - phip)) + 6*Cos(2*(phipol - phip)) + 
     -    Cos(2*(lamp + phipol - phip)) - 4*Cos(2*phip) + 
     -    2*Cos(2*(lamp + phip)) + Cos(2*(lamp - phipol + phip)) + 
     -    6*Cos(2*(phipol + phip)) + Cos(2*(lamp + phipol + phip)) + 
     -    4*Cos(lamp - 2*phipol + 2*phip) - 
     -    4*Cos(lamp + 2*phipol + 2*phip))

      return
      end
c-----------------------------------------------------------------
      real function uvtovgm(up,vp,lonp,latp,lonpol,latpol)
c-----------------------------------------------------------------
c     Rene Fehlmann, Wed Nov  1 13:31:28 MET 1995
c     INPUT:
c     - rotated winds (i.e. wrt equatorial grid)
c     - rotated coordinates (i.e. on equatorial grid)
c     - coordinates of the rotated pole
c     OUTPUT
c     - v-component of the wind (i.e. real N-S wind)
c-----------------------------------------------------------------

      real               up,vp,lonp,latp,lonpol,latpol
      real               lampol,phipol,lamp,phip

      data               pid180 / 0.0174532925199 /

c-----------------------------------------------------------------

      phip=latp*pid180
      lamp=lonp*pid180
      phipol=latpol*pid180
      lampol=lonpol*pid180

      uvtovgm=(vp*Cos(phip)*Sin(phipol) + 
     -     Cos(phipol)*(-(up*Sin(lamp)) -
     -     vp*Cos(lamp)*Sin(phip)))/
     -     Sqrt(1 - (Cos(lamp)*Cos(phipol)*Cos(phip) + 
     -       Sin(phipol)*Sin(phip))**2)

      return
      end