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



      real,allocatable,dimension (:,:)     :: ps    ! surface pressure
      real,allocatable,dimension (:,:,:)   :: z3    ! 3D geopotential height
      real,allocatable,dimension (:,:)     :: zb    ! surface geopotential height
 
c     Transformation between numerical and geographical grid
      integer,allocatable,dimension (:,:)  :: connect1 
      integer      connectval1 
      integer,allocatable,dimension (:,:)  :: connect2 
      integer      connectval2
      real,allocatable,dimension (:,:)  :: xc,yc
      real         radius
      real         xval,yval
 

c     Write grid information to screen
      print*,' xmin,xmax  = ',xmin,xmax
      print*,' ymin,ymax  = ',ymin,ymax
      print*,' dx,dy      = ',dx,dy
      print*,' nx,ny,nz   = ',nx,ny,nz 
      print*,' anglemode  = ',trim(anglemode)
 
c     Transform grid specifier into grid coordinates
      xmin = 1.
      ymin = 1.
      dx   = 1.

 
      cdfname = mapfile
      varname = 'P'
      call writecdf2D(cdfname,varname,p3,0.,
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
 
      cdfname = mapfile
      varname = 'TOPO'
      call writecdf2D(cdfname,varname,zb,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 

 
 
c     ----------------------------------------------------------------------
c     Get spherical distance between lat/lon points
c     ----------------------------------------------------------------------
        
      real function sdis(xp,yp,xq,yq,unit)
 
c     Calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
 

      character*80 unit
      real         dlon,dlat,alpha,rlat1,ralt2
 
      if ( unit.eq.'km' ) then
 
         arg=sin(yp*deg2rad)*sin(yq*deg2rad)+
     >           cos(yp*deg2rad)*cos(yq*deg2rad)*cos((xp-xq)*deg2rad)
         if (arg.lt.-1.) arg=-1.
         if (arg.gt.1.) arg=1.
         sdis=re*acos(arg)
         
      elseif ( unit.eq.'deg' ) then

Rev 11	Rev 15
Line 43...	Line 43...
43	`real,allocatable,dimension (:,:) :: ps ! surface pressure`	43	`real,allocatable,dimension (:,:) :: ps ! surface pressure`
44	`real,allocatable,dimension (:,:,:) :: z3 ! 3D geopotential height`	44	`real,allocatable,dimension (:,:,:) :: z3 ! 3D geopotential height`
45	`real,allocatable,dimension (:,:) :: zb ! surface geopotential height`	45	`real,allocatable,dimension (:,:) :: zb ! surface geopotential height`
46		46
47	`c Transformation between numerical and geographical grid`	47	`c Transformation between numerical and geographical grid`
48	`real,allocatable,dimension (:,:) :: connect1`	48	`integer,allocatable,dimension (:,:) :: connect1`
49	`integer connectval1`	49	`integer connectval1`
50	`real,allocatable,dimension (:,:) :: connect2`	50	`integer,allocatable,dimension (:,:) :: connect2`
51	`integer connectval2`	51	`integer connectval2`
52	`real,allocatable,dimension (:,:) :: xc,yc`	52	`real,allocatable,dimension (:,:) :: xc,yc`
53	`real radius`	53	`real radius`
54	`real xval,yval`	54	`real xval,yval`
55		55
Line 233...	Line 233...
233	`c Write grid information to screen`	233	`c Write grid information to screen`
234	`print*,' xmin,xmax = ',xmin,xmax`	234	`print*,' xmin,xmax = ',xmin,xmax`
235	`print*,' ymin,ymax = ',ymin,ymax`	235	`print*,' ymin,ymax = ',ymin,ymax`
236	`print*,' dx,dy = ',dx,dy`	236	`print*,' dx,dy = ',dx,dy`
237	`print*,' nx,ny,nz = ',nx,ny,nz`	237	`print*,' nx,ny,nz = ',nx,ny,nz`
-		238	`print*,' anglemode = ',trim(anglemode)`
238		239
239	`c Transform grid specifier into grid coordinates`	240	`c Transform grid specifier into grid coordinates`
240	`xmin = 1.`	241	`xmin = 1.`
241	`ymin = 1.`	242	`ymin = 1.`
242	`dx = 1.`	243	`dx = 1.`
Line 450...	Line 451...
450		451
451	`cdfname = mapfile`	452	`cdfname = mapfile`
452	`varname = 'P'`	453	`varname = 'P'`
453	`call writecdf2D(cdfname,varname,p3,0.,`	454	`call writecdf2D(cdfname,varname,p3,0.,`
454	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`	455	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`
455		-
456	`cdfname = mapfile`	456	`cdfname = mapfile`
457	`varname = 'TOPO'`	457	`varname = 'TOPO'`
458	`call writecdf2D(cdfname,varname,zb,0.,`	458	`call writecdf2D(cdfname,varname,zb,0.,`
459	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`	459	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`
460		460
Line 1177...	Line 1177...
1177		1177
1178		1178
1179	`c ----------------------------------------------------------------------`	1179	`c ----------------------------------------------------------------------`
1180	`c Get spherical distance between lat/lon points`	1180	`c Get spherical distance between lat/lon points`
1181	`c ----------------------------------------------------------------------`	1181	`c ----------------------------------------------------------------------`
1182		1182
1183	`real function sdis(xp,yp,xq,yq,unit)`	1183	`real function sdis(xp,yp,xq,yq,unit)`
1184		1184
1185	`c Calculates spherical distance (in km) between two points given`	1185	`c Calculates spherical distance (in km) between two points given`
1186	`c by their spherical coordinates (xp,yp) and (xq,yq), respectively.`	1186	`c by their spherical coordinates (xp,yp) and (xq,yq), respectively.`
1187		1187
Line 1193...	Line 1193...
1193	`character*80 unit`	1193	`character*80 unit`
1194	`real dlon,dlat,alpha,rlat1,ralt2`	1194	`real dlon,dlat,alpha,rlat1,ralt2`
1195		1195
1196	`if ( unit.eq.'km' ) then`	1196	`if ( unit.eq.'km' ) then`
1197		1197
1198	`arg=sind(yp)sind(yq)+cosd(yp)cosd(yq)*cosd(xp-xq)`	1198	`arg=sin(ypdeg2rad)sin(yq*deg2rad)+`
-		1199	`> cos(ypdeg2rad)cos(yqdeg2rad)cos((xp-xq)*deg2rad)`
1199	`if (arg.lt.-1.) arg=-1.`	1200	`if (arg.lt.-1.) arg=-1.`
1200	`if (arg.gt.1.) arg=1.`	1201	`if (arg.gt.1.) arg=1.`
1201	`sdis=re*acos(arg)`	1202	`sdis=re*acos(arg)`
1202		1203
1203	`elseif ( unit.eq.'deg' ) then`	1204	`elseif ( unit.eq.'deg' ) then`

Subversion Repositories lagranto.wrf

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