subroutine bilinb(cob,inb,jnb,kb,nx,ny,nz,ain,aout)
c
c *** Bilinear interpolation routine.
c
implicit none
c
integer*4 i00,i10,i01,i11
. ,j00,j10,j01,j11
. ,nx,ny,nz,kb,k,l
c
real*4 cob(3,kb),inb(4,kb),jnb(4,kb),ain(nx,ny,nz),aout(kb,nz)
. ,p,q,pq
c_______________________________________________________________________________
c
do l=1,nz
do k=1,kb
c
i00=inb(1,k)
i10=inb(2,k)
i01=inb(3,k)
i11=inb(4,k)
c
j00=jnb(1,k)
j10=jnb(2,k)
j01=jnb(3,k)
j11=jnb(4,k)
c
p =cob(2,k)
q =cob(3,k)
pq=cob(1,k)
c
aout(k,l)=ain(i00,j00,l)
. +p*(ain(i10,j10,l)-ain(i00,j00,l))
. +q*(ain(i01,j01,l)-ain(i00,j00,l))
. +pq*(ain(i00,j00,l)-ain(i10,j10,l)
. -ain(i01,j01,l)+ain(i11,j11,l))
c
enddo
enddo
c
return
end
c
c===============================================================================
Cmp
subroutine bilinb_sfc(cob,inb,jnb,kb,nx,ny,ain,aout)
c
c *** Bilinear interpolation routine.
c
implicit none
c
integer*4 i00,i10,i01,i11
. ,j00,j10,j01,j11
. ,nx,ny,nz,kb,k,l
c
real*4 cob(3,kb),inb(4,kb),jnb(4,kb),ain(nx,ny),
. aout(kb)
. ,p,q,pq
c_______________________________________________________________________________
c
do k=1,kb
c
i00=inb(1,k)
i10=inb(2,k)
i01=inb(3,k)
i11=inb(4,k)
c
j00=jnb(1,k)
j10=jnb(2,k)
j01=jnb(3,k)
j11=jnb(4,k)
c
p =cob(2,k)
q =cob(3,k)
pq=cob(1,k)
c
aout(k)=ain(i00,j00)
. +p*(ain(i10,j10)-ain(i00,j00))
. +q*(ain(i01,j01)-ain(i00,j00))
. +pq*(ain(i00,j00)-ain(i10,j10)
. -ain(i01,j01)+ain(i11,j11))
c
enddo
c
return
end
c
Cmp
c===============================================================================
c
subroutine ltlwin(almd,aphd,imjm1,tlm0d,ctph0,stph0,pus,pvs)
c
c *** ll to tll transformation of velocity.
c *** Programer: Z. Janjic, Yugoslav Fed. Hydromet. Inst.,
c Beograd, 1982
c
implicit none
c
integer*4 imjm1,l
c
real*4 pus(imjm1),pvs(imjm1),lpus,lpvs
. ,almd(imjm1),aphd(imjm1)
. ,tlm0d,ctph0,stph0
. ,relmd,srlm,crlm,sph,cph,cc,tphd,rctph,cray,dray
c_______________________________________________________________________________
c
do l=1,imjm1
c
relmd=almd(l)-tlm0d
srlm=sind(relmd)
crlm=cosd(relmd)
c
sph=sind(aphd(l))
cph=cosd(aphd(l))
c
cc=cph*crlm
tphd=asin(ctph0*sph-stph0*cc)
rctph=1./cosd(tphd)
c
cray=stph0*srlm*rctph
dray=(ctph0*cph+stph0*sph*crlm)*rctph
c
lpus=pus(l)
lpvs=pvs(l)
pus(l)=dray*lpus-cray*lpvs
pvs(l)=cray*lpus+dray*lpvs
enddo
c
return
end
c
c===============================================================================
c
subroutine blwts(coh,inh,jnh,cov,inv,jnv,ald,apd
. ,gproj)
c
c *** Compute bilinear interpolation weights.
c
implicit none
c
include 'ecommons.h'
c
real*4 coh(3,imjm),inh(4,imjm),jnh(4,imjm)
. ,cov(3,imjm1),inv(4,imjm1),jnv(4,imjm1)
. ,ald(imjm1),apd(imjm1)
. ,khl0(jm),kvl0(jm),khh0(jm),kvh0(jm)
. ,tdlmd
. ,tphd,tlmd,aphd,almd,x,y
c
integer*4 khl,khh,kvl,kvh
. ,i,j,k
c
character*2 gproj
c
integer*4 nxsec,nysec,nzsec
common /sectorsize/nxsec,nysec,nzsec
c_______________________________________________________________________________
c
C write(6,*) 'in blwts:'
C print *,'im=',im,'jm=',jm
C print *,'Input data dimensions:',nxsec,nysec,nzsec
c
c *** Hibu model conts.
c
do j=1,jm
khl0(j)=im*(j-1)-(j-1)/2+1
kvl0(j)=im*(j-1) - j/2+1
khh0(j)=im* j - j/2
kvh0(j)=im* j -(j+1)/2
enddo
c
c *** Define some constants.
c
tdlmd=dlmd+dlmd
c
c *** Compute weights for hibu height pts.
c
tphd=sbd-dphd
c
do j=1,jm
khl=khl0(j)
khh=khh0(j)
tphd=tphd+dphd
c
do k=khl,khh
tlmd=wbd+tdlmd*(k-khl)*mod(j,2)
. +mod(j+1,2)*(dlmd+(k-khl)*tdlmd)
c
c ********* tll to ll conversion.
c
call rtll(tlmd,tphd,tlm0d,ctph0,stph0,almd,aphd)
c
c ********* Determine location of eta grid point in native data i,j space.
c
if (gproj .eq. 'PS') then
Cmp call latlon_2_psij(1,aphd,almd,x,y)
Cmp
Cmp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Cmp Below routine is hardwired for grid104 GDS!!!!!!
Cmp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Cmp
call str_ij(aphd,almd,x,y)
Cmp
elseif (gproj .eq. 'LC') then
call latlon_2_lcij(1,aphd,almd,x,y)
elseif (gproj .eq. 'CE') then
call latlon_2_coneqij(1,aphd,almd,x,y)
elseif (gproj .eq. 'LL') then
call latlon_2_llij(1,aphd,almd,x,y)
endif
C if (aphd .gt. 43 .and. aphd .lt. 46) then
if (almd .gt. -105.3 .and. almd .lt. -104.7) then
C write(6,*) 'aphd,almd,x,y', aphd,almd,x,y
endif
c
c ********* Check if hibu pt is out of native data domain.
c
c print *,'k,lat,lon,x,y ',k,aphd,almd,x,y
if (x .lt. 1. .or. x .gt. float(nxsec) .or.
. y .lt. 1. .or. y. gt. float(nysec)) then
Cmp print *,'hibupt outside domain k=',k,aphd,almd,x,y
Cmp if 360 X 181 grid assume periodicity and adjust the x index accordingly
if (nxsec .eq. 360) then
if (x .lt. 1) x=x+360.
if (x .gt. nxsec) x=x-360.
else
print *,'hibupt outside domain k=',k,aphd,almd,x,y
stop
endif
endif
c
coh(2,k)=x-float(int(x))
coh(3,k)=y-float(int(y))
coh(1,k)=coh(2,k)*coh(3,k)
c
inh(1,k)=int(x)
inh(2,k)=inh(1,k)+1
inh(3,k)=inh(1,k)
inh(4,k)=inh(1,k)+1
c
jnh(1,k)=int(y)
jnh(2,k)=jnh(1,k)
jnh(3,k)=jnh(1,k)+1
jnh(4,k)=jnh(1,k)+1
c
enddo
enddo
c
c *** Compute weights for hibu wind points
c
tphd=sbd-dphd
c
do j=1,jm
kvl=kvl0(j)
kvh=kvh0(j)
tphd=tphd+dphd
c
do k=kvl,kvh
tlmd=wbd+dlmd*mod(j,2)+(k-kvl)*tdlmd
c
c ********* tll to ll conversion.
c
call rtll(tlmd,tphd,tlm0d,ctph0,stph0,almd,aphd)
c
ald(k)=almd
apd(k)=aphd
c
c ********* Determine location of eta grid point in native data i,j space.
c
if (gproj .eq. 'PS') then
Cnot call latlon_2_psij(1,aphd,almd,x,y)
Cmp
call str_ij(aphd,almd,x,y)
Cmp
elseif (gproj .eq. 'LC') then
call latlon_2_lcij(1,aphd,almd,x,y)
elseif (gproj .eq. 'CE') then
call latlon_2_coneqij(1,aphd,almd,x,y)
elseif (gproj .eq. 'LL') then
call latlon_2_llij(1,aphd,almd,x,y)
endif
c
c *** check if hibu pt is out of native data domain.
c
if (x .lt. 1. .or. x .gt. float(nxsec) .or.
. y .lt. 1. .or. y. gt. float(nysec)) then
cmp print *,'wind point outside domain k=',k,aphd,almd,x,y
Cmp if 360 X 181 grid assume periodicity and adjust the x index accordingly
if (nxsec .eq. 360) then
if (x .lt. 1) x=x+360.
if (x .gt. nxsec) x=x-360.
else
print *,'wind point outside domain k=',k,aphd,almd,x,y
stop
endif
endif
c
cov(2,k)=x-float(int(x))
cov(3,k)=y-float(int(y))
cov(1,k)=cov(2,k)*cov(3,k)
c
inv(1,k)=int(x)
inv(2,k)=inv(1,k)+1
inv(3,k)=inv(1,k)
inv(4,k)=inv(1,k)+1
c
jnv(1,k)=int(y)
jnv(2,k)=jnv(1,k)
jnv(3,k)=jnv(1,k)+1
jnv(4,k)=jnv(1,k)+1
c
enddo
enddo
c
return
end
c
c===============================================================================
c
subroutine get_sector_size(fname)
c
implicit none
c
integer*4 nx,ny,nz,nxsec,nysec,ip,jp,nsfcfld,l
character*(*) fname
character*2 gproj
c
CTO common/sectorsize/nxsec,nysec,nz,ip,jp
common/sectorsize/nxsec,nysec,nz
c
c *** Polar stereographic common block.
c
integer*4 psnx,psny,psnz !no. of ps domain grid points
real*4 pslat0,pslon0,psrota, !pol ste. std lat, lon and rotation
. pssw(2),psne(2) !sw lat, lon, ne lat, lon
common /psgrid/psnx,psny,psnz,pslat0,pslon0,psrota,pssw,psne
c
c *** Lambert-conformal common block.
c
integer*4 lcnx,lcny,lcnz !no. of lc domain grid points
real*4 lclat1,lclat2,lclon0, !lambert-conformal std lat1, lat, lon
. lcsw(2),lcne(2) !sw lat, lon, ne lat, lon
common /lcgrid/lcnx,lcny,lcnz,lclat1,lclat2,lclon0,lcsw,lcne
c
c *** Conical equidistant common block.
c
integer*4 cenx,ceny,cenz
real*4 celat0,celon0,dphi,dlam
common /coneqgrid/cenx,ceny,cenz,celat0,celon0,dphi,dlam
c
c *** Lat-lon common block.
c
integer*4 llnx,llny,llnz !no. of ll domain grid points
real*4 lllat0,lllon0,dlat,dlon
common /llgrid/llnx,llny,llnz,lllat0,lllon0,dlat,dlon
Cmp
Cmp
Cmp lat0,lon0 is SW corner....dlat,dlon is grid spacing
c_______________________________________________________________________________
c
print *,' '
l=index(fname//' ',' ')-1
print *,'Ingesting data - file = ',fname(1:l)
open(1,file=fname,status='old',form='unformatted')
c
c *** Read sector grid information.
c
read(1) nx,ny,nz,nxsec,nysec,ip,jp,nsfcfld,gproj
c
c *** Read grid projection specific information.
c
if (gproj .eq. 'PS') then
read(1) psnx,psny,psnz,pslat0,pslon0,psrota,pssw,psne
elseif (gproj .eq. 'LC') then
read(1) lcnx,lcny,lcnz,lclat1,lclat2,lclon0,lcsw,lcne
elseif (gproj .eq. 'CE') then
read(1) cenx,ceny,cenz,celat0,celon0,dphi,dlam
elseif (gproj .eq. 'LL') then
read(1) llnx,llny,llnz,lllat0,lllon0,dlat,dlon
else
print *,'Unrecognizable grid projection - ',gproj
print *,' abort.'
stop
endif
c
close(1)
return
end
c
c===============================================================================
c
Cmp subroutine getdata(fname,pr,uw,vw,ht,mr,gproj)
subroutine getdata(fname,pr,uw,vw,ht,mr,th,sfcgrid,gproj,
+ nsfcfld)
c
c *** Read input upper air data.
c
implicit none
c
real*4 cpi,kappi
parameter (cpi=1./1004.,kappi=1004./287.)
c
integer*4 nx,ny,nz,dummy(7),nsfcfld,i,j,k,pres,ii,jj
c
Cmp real*4 sfcgrid(nx,ny,nsfcfld)
real*4 sfcgrid(nx,ny,12),tmp
c
c *** Output data variables.
c
real*4 uw(nx,ny,nz) !u-wind (m/s)
. ,vw(nx,ny,nz) !v-wind (m/s)
. ,th(nx,ny,nz) !theta (k)
. ,ht(nx,ny,nz) !height (m)
. ,mr(nx,ny,nz) !mixing ratio (kg/kg)
. ,ex(nx,ny,nz) !exner
. ,pr(nz) !pressure (Pa)
c
character*(*) fname
character*2 gproj
c
common /sectorsize/nx,ny,nz
c_______________________________________________________________________________
c
print *,'In getdata'
print *,' - data sector size: ',nx,ny,nz
c print *,' '
c
c *** Open data file.
c
open(1,file=fname,status='old',form='unformatted',err=100)
c
c *** Skip through header stuff and surface fields.
c
rewind(1)
read(1) dummy,nsfcfld,gproj
read(1) dummy(1)
C do i=1,nsfcfld
C read(1) sfcgrid
C enddo
c
c *** Read upper air data.
c (Read in upside down since that is how ETA wants it.)
c
read(1) (((uw(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
read(1) (((vw(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
read(1) (((th(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
read(1) (((ht(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
read(1) (((mr(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
read(1) (((ex(i,j,k),i=1,nx),j=1,ny),k=nz,1,-1)
cfhs
c print *," ex in dutil "
c call maxmn3(ex,nx,ny,nz,"ex ")
c
do k=1,nz
C write(6,*) (ht(i,20,k),I=20,30)
enddo
Cmp
write(6,*) 'nsfcfld= ', nsfcfld
C do k=1,nsfcfld
C write(6,*) 'reading k= ', k
Cmp read(1) ((sfcgrid(i,j,k),i=1,nx),j=1,ny)
read(1) sfcgrid
C enddo
C write(6,*) 'thta values after getdata read (s-->n, 180W)'
do j=1,ny
C write(6,*) th(nx/2,j,nz)
enddo
close(1)
c
c *** Computer pressure from Exner.
c
write(6,*) 'in dutil...computing pr levels'
do k=1,nz
Cmp pr(k)=100000.*(ex(1,1,k)*cpi)**kappi
Cmp changed line to avoid problems with 221 data (no data at 1,1)
pr(k)=100000.*(ex(nx/2,ny/2,k)*cpi)**kappi
cfhs
cfhs print *,"pr inside dutil.f",pr(k),k
enddo
c
return
100 print*,'Error reading upper data.'
c
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine str_ij(RLAT,RLON,XPTS,YPTS)
INTEGER KGDS(12)
REAL XPTS,YPTS,RLON,RLAT
PARAMETER(RERTH=6.3712E6)
PARAMETER(PI=3.14159265358979,DPR=180./PI)
C
C Hardwired for Grid 104
data kgds/005,147,110,-268,-139475,8,-105000,90755,90755,
+0,64,0/
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IM=KGDS(2)
JM=KGDS(3)
RLAT1=KGDS(4)*1.E-3
RLON1=KGDS(5)*1.E-3
IROT=MOD(KGDS(6)/8,2)
ORIENT=KGDS(7)*1.E-3
DX=KGDS(8)
DY=KGDS(9)
IPROJ=MOD(KGDS(10)/128,2)
ISCAN=MOD(KGDS(11)/128,2)
JSCAN=MOD(KGDS(11)/64,2)
NSCAN=MOD(KGDS(11)/32,2)
H=(-1.)**IPROJ
HI=(-1.)**ISCAN
HJ=(-1.)**(1-JSCAN)
DXS=DX*HI
DYS=DY*HJ
DE=(1.+SIN(60./DPR))*RERTH
DR=DE*COS(RLAT1/DPR)/(1+H*SIN(RLAT1/DPR))
XP=1-H*SIN((RLON1-ORIENT)/DPR)*DR/DXS
YP=1+COS((RLON1-ORIENT)/DPR)*DR/DYS
DE2=DE**2
XMIN=0
XMAX=IM+1
YMIN=0
YMAX=JM+1
NRET=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C TRANSLATE EARTH COORDINATES TO GRID COORDINATES
IOPT=-1
IF(IOPT.EQ.-1) THEN
IF(ABS(RLON).LE.360.AND.ABS(RLAT).LE.90.AND.
& H*RLAT.NE.-90) THEN
DR=DE*TAN((90-H*RLAT)/2/DPR)
C write(6,*) 'DE,DR ', de,dr
XPTS=XP+H*SIN((RLON-ORIENT)/DPR)*DR/DXS
YPTS=YP-COS((RLON-ORIENT)/DPR)*DR/DYS
C write(6,*) 'xpts,ypts ', xpts,ypts
IF(XPTS.GE.XMIN.AND.XPTS.LE.XMAX.AND.
& YPTS.GE.YMIN.AND.YPTS.LE.YMAX) THEN
NRET=NRET+1
ELSE
XPTS=FILL
YPTS=FILL
ENDIF
ELSE
XPTS=FILL
YPTS=FILL
ENDIF
ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
return
END