WebSVN – lagranto.wrf – Diff – /trunk/goodies/wrfmap.f

       character*80 radunit
       real         rdummy
       real         rid,rjd,rkd
       real         xpos,ypos,zpos,ppos,lonpos,latpos
       real         lon1,lat1,lon2,lat2
+      integer      nfilter
 c     Externals
       real         int_index3
       external     int_index3
       real         sdis
        if ( mode.eq.'-z2p.single'  ) then     ! convert single x/y/z -> x/y/p
           read(10,*) xpos,ypos,zpos
        endif
+       if ( mode.eq.'-mapscale'  ) then       ! calculate mapping scale factors
+          read(10,*) nfilter
+       endif
       close(10)
 c     Read the mapping file - except for mode '-create'
       if ( mode.ne.'-create' ) then
       if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 c     Get format of input and output file
       call mode_tra(inpmode,inpfile)
       call mode_tra(outmode,outfile)
-      if ( outmode.eq.-1 ) outmode = 1
+      if ( outmode.eq.-1 ) outmode = inpmode
 c     Read coordinates from file - Format (lon,lat,lev)
       if (inpmode.eq.-1) then
          fid = 10
          npoints = 0
       if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 c     Get format of input and output file
       call mode_tra(inpmode,inpfile)
       call mode_tra(outmode,outfile)
-      if ( outmode.eq.-1) outmode=1
+      if ( outmode.eq.-1) outmode=inpmode
 c     Read coordinates from file - Format (x,y,lev)
       if (inpmode.eq.-1) then
          fid = 10
          npoints = 0
       mpx(nx,ny) = mpx(nx-1,ny-1)
       mpy(nx,ny) = mpy(nx-1,ny-1)
       mpx( 1,ny) = mpx(   2,ny-1)
       mpy( 1,ny) = mpy(   2,ny-1)
+c     Apply a diffuive filter
+      do i=1,nfilter
+         call filt2d (mpx,mpx,nx,ny,1.0,-999.,0,0,0,0)
+         call filt2d (mpy,mpy,nx,ny,1.0,-999.,0,0,0,0)
+      enddo
 c     Write map factors to mapping file
       cdfname = mapfile
       varname = 'MAPSCALE_X'
       call writecdf2D(cdfname,varname,mpx,0.,
      >                dx,dy,xmin,ymin,nx,ny,1,0,1)
       if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 c     Get format of input and output file
       call mode_tra(inpmode,inpfile)
       call mode_tra(outmode,outfile)
-      if ( outmode.eq.-1) outmode=1
+      if ( outmode.eq.-1) outmode=inpmode
 c     Read coordinates from file - Format (x,y,lev)
       if (inpmode.eq.-1) then
          fid = 10
          npoints = 0
       if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 c     Get format of input and output file
       call mode_tra(inpmode,inpfile)
       call mode_tra(outmode,outfile)
-      if ( outmode.eq.-1) outmode=1
+      if ( outmode.eq.-1) outmode=inpmode
 c     Read coordinates from file - Format (x,y,lev)
       if (inpmode.eq.-1) then
          fid = 10
          npoints = 0
       real varmin(4),varmax(4),stag(4)
       integer ierr,cdfid,ndim,vardim(4)
       character*80 cstname
       real mdv
       integer i
+      integer nvars
+      character*80 vars(300)
 C     Definitions
       varmin(1)=xmin
       varmin(2)=ymin
       varmin(3)=1050.
          call crecdf(cdfname,cdfid,
      >        varmin,varmax,ndim,cstname,ierr)
          if (ierr.ne.0) goto 903
       endif
+c     Check whether the variable is already on the file
+      call getvars(cdfid,nvars,vars,ierr)
+      if (ierr.ne.0) goto 906
+      do i=1,nvars
+         if ( vars(i).eq.varname ) crevar = 0
+      enddo
 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
       stop
       end
+c     -----------------------------------------------
+c     Diffusive filter
+c     -----------------------------------------------
+      subroutine filt2d (a,af,nx,ny,fil,misdat,
+     &                   iperx,ipery,ispol,inpol)
+c     Apply a conservative diffusion operator onto the 2d field a,
+c     with full missing data checking.
+c
+c     a     real   inp  array to be filtered, dimensioned (nx,ny)
+c     af    real   out  filtered array, dimensioned (nx,ny), can be
+c                       equivalenced with array a in the calling routine
+c     f1    real        workarray, dimensioned (nx+1,ny)
+c     f2    real        workarray, dimensioned (nx,ny+1)
+c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
+c                       corresponds to one application of the linear filter.
+c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
+c                       the corresponding value is not available. The
+c                       misdat-checking can be switched off with with misdat=0.
+c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
+c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
+c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
+c     ispol int    inp  southpole at j=1   (1=yes,0=no)
+c
+c     Christoph Schaer, 1993
+c     argument declaration
+      integer     nx,ny
+      real        a(nx,ny),af(nx,ny),fil,misdat
+      integer     iperx,ipery,inpol,ispol
+c     local variable declaration
+      integer     i,j,is
+      real        fh
+      real        f1(nx+1,ny),f2(nx,ny+1)
+c     compute constant fh
+      fh=0.125*fil
+c     compute fluxes in x-direction
+      if (misdat.eq.0.) then
+        do j=1,ny
+        do i=2,nx
+          f1(i,j)=a(i-1,j)-a(i,j)
+        enddo
+        enddo
+      else
+        do j=1,ny
+        do i=2,nx
+          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
+            f1(i,j)=0.
+          else
+            f1(i,j)=a(i-1,j)-a(i,j)
+          endif
+        enddo
+        enddo
+      endif
+      if (iperx.eq.1) then
+c       do periodic boundaries in the x-direction
+        do j=1,ny
+          f1(1,j)=f1(nx,j)
+          f1(nx+1,j)=f1(2,j)
+        enddo
+      else
+c       set boundary-fluxes to zero
+        do j=1,ny
+          f1(1,j)=0.
+          f1(nx+1,j)=0.
+        enddo
+      endif
+c     compute fluxes in y-direction
+      if (misdat.eq.0.) then
+        do j=2,ny
+        do i=1,nx
+          f2(i,j)=a(i,j-1)-a(i,j)
+        enddo
+        enddo
+      else
+        do j=2,ny
+        do i=1,nx
+          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
+            f2(i,j)=0.
+          else
+            f2(i,j)=a(i,j-1)-a(i,j)
+          endif
+        enddo
+        enddo
+      endif
+c     set boundary-fluxes to zero
+      do i=1,nx
+        f2(i,1)=0.
+        f2(i,ny+1)=0.
+      enddo
+      if (ipery.eq.1) then
+c       do periodic boundaries in the x-direction
+        do i=1,nx
+          f2(i,1)=f2(i,ny)
+          f2(i,ny+1)=f2(i,2)
+        enddo
+      endif
+      if (iperx.eq.1) then
+        if (ispol.eq.1) then
+c         do south-pole
+          is=(nx-1)/2
+          do i=1,nx
+            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
+          enddo
+        endif
+        if (inpol.eq.1) then
+c         do north-pole
+          is=(nx-1)/2
+          do i=1,nx
+            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
+          enddo
+        endif
+      endif
+c     compute flux-convergence -> filter
+      if (misdat.eq.0.) then
+        do j=1,ny
+        do i=1,nx
+            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
+        enddo
+        enddo
+      else
+        do j=1,ny
+        do i=1,nx
+          if (a(i,j).eq.misdat) then
+            af(i,j)=misdat
+          else
+            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
+          endif
+        enddo
+        enddo
+      endif
+      end

Subversion Repositories lagranto.wrf

(root)/trunk/goodies/wrfmap.f – Rev 8 → 11