subroutine geodata ( geobus,topo,lminx,lmaxx,lminy,lmaxy, 1,20
$ maxhh01,maxhh02,dimgx,dimgy,unf)
implicit none
*
integer dimgx,dimgy,lminx,lmaxx,lminy,lmaxy,unf
real geobus(*),topo(lminx-1:lmaxx,lminy-1:lmaxy,2),maxhh01,maxhh02
*
#include "lcldim.cdk"
#include "grd.cdk"
#include "rec.cdk"
#include "ifd.cdk"
#include "partopo.cdk"
#include "hinterpo.cdk"
#include "lesbus.cdk"
#include "filename.cdk"
#include "yomdyn.cdk"
#include "topo.cdk"
*
integer fnom,fstouv,fstinf,fstprm,fstluk,fstfrm,fclos,
$ fstopc,nav_3df
external fnom,fstouv,fstinf,fstprm,fstluk,fstfrm,fclos,
$ fstopc,nav_3df
*
character*1 typ, grd
character*2 var
character*8 lab, inttyp
character*256 fn
integer dte, det, ipas, p1, p2, p3, g1, g2, g3, g4, bit,
$ dty, swa, lng, dlf, ubc, ex1, ex2, ex3, longueur
*
integer unf2,errop,ier,ni1,nj1,nk1,key,ni,nj,n,nia,nja,
$ i,j,mgi,z0i,lhi,y7i,y8i,y9i,gai,vgi,vfi,gid,
$ lai,loi,mti,mfi,niext,njext,errrot
integer, dimension (:), allocatable :: idx,idy
real xlat1,xlon1,xlat2,xlon2,epsil
real , dimension (: ) , allocatable :: xps,yps
real , dimension (:,:) , allocatable :: meh
real*8, dimension (: ) , allocatable :: xpaq,ypaq,
$ cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd
real*8 xpxext(0:dimgx), ypxext(0:dimgy)
data epsil /1.0e-5/
*-----------------------------------------------------------------------
*
inttyp = 'LINEAR'
if (myproc.ne.0) ier = fstopc('MSGLVL','SYSTEM',.false.)
*
do i=1,dimgx
xpxext(i) = xpx(i)
end do
do i=1,dimgy
ypxext(i) = ypx(i)
end do
xpxext(0) = xpxext(1) - (xpxext(2)-xpxext(1))
ypxext(0) = ypxext(1) - (ypxext(2)-ypxext(1))
*
unf2 = 0
*
errrot = 0
c ier = nav_3df (unf,0,0,0.)
ier = nav_3df (unf,hx,hy,1.2)
nia = ifd_niaf - ifd_niad + 1
nja = ifd_njaf - ifd_njad + 1
close (unf)
*
do n=1,ifd_nf
if (ifd_needit(n)) then
fn =
$ '../geophy/'//prefgeo(1:longueur(prefgeo))//'_'//ifd_fnext(n)
ier = fnom (unf2,fn,'RND+OLD+R/O',0)
ier = fstouv (unf2,'RND')
*
* Use first file to establish geo-references
*
allocate ( xpaq(nia), ypaq(nja) )
key = fstinf(unf2,ni1,nj1,nk1,-1,' ',-1,-1,-1,' ','>>')
allocate (xps(ni1))
ier = fstluk ( xps, key, ni,nj1,nk1 )
key = fstinf(unf2,ni1,nj1,nk1,-1,' ',-1,-1,-1,' ','^^')
allocate (yps(nj1))
ier = fstluk ( yps, key, ni1,nj,nk1 )
do i=1,nia
xpaq(i) = xps(ifd_niad+i-1)
end do
do j=1,nja
ypaq(j) = yps(ifd_njad+j-1)
end do
deallocate (xps,yps)
*
ier= fstprm (key, DTE, DET, IPAS, ni1, nj1, nk1, BIT, DTY,
$ P1, P2, P3, TYP, VAR, LAB, GRD, G1, G2, G3, G4,
$ SWA, LNG, DLF, UBC, EX1, EX2, EX3)
call cigaxg (GRD,xlat1,xlon1,xlat2,xlon2,g1,g2,g3,g4)
if ((abs(Grd_xlat1-xlat1)/(Grd_xlat1+91.).gt.epsil) .or.
$ (abs(Grd_xlat2-xlat2)/(Grd_xlat2+91.).gt.epsil) .or.
$ (abs(Grd_xlon1-xlon1)/(Grd_xlon1+ 1.).gt.epsil) .or.
$ (abs(Grd_xlon2-xlon2)/(Grd_xlon2+ 1.).gt.epsil) )
$ errrot = -1
ier = fstfrm (unf2)
ier = fclos (unf2)
goto 57
endif
enddo
*
57 if ((errrot.lt.0).and.(myproc.eq.0)) write (6,1001)
call mc2stop(errrot)
*
niext = ldni+2*hx+1
njext = ldnj+2*hy+1
allocate (idx(niext), idy(njext))
allocate (cxa(niext),cxb(niext),cxc(niext),cxd(niext),
$ cya(njext),cyb(njext),cyc(njext),cyd(njext),
$ meh(niext,njext))
*
call grid_to_grid_coef (xpxext(gc_ld(1,myproc)-1),niext,xpaq,
$ nia,idx,cxa,cxb,cxc,cxd,inttyp)
call grid_to_grid_coef (ypxext(gc_ld(3,myproc)-1),njext,ypaq,
$ nja,idy,cya,cyb,cyc,cyd,inttyp)
*
call fstrhint (meh,'ME',niext,njext,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
*
call dc_topo2 ( topo,meh,topo_flt_coef,minx,maxx,miny,maxy,
$ ldni,ldnj,hx,hy,period_x,period_y,
$ maxhh01,maxhh02,myproc.eq.0 )
*
deallocate (meh)
*
if (gnmaphy.gt.0) then
do gid=1,geotop
if (geonm(gid,2).eq.'MG') mgi = geopar(gid,1)
if (geonm(gid,2).eq.'Z0') z0i = geopar(gid,1)
if (geonm(gid,2).eq.'LH') lhi = geopar(gid,1)
if (geonm(gid,2).eq.'Y7') y7i = geopar(gid,1)
if (geonm(gid,2).eq.'Y8') y8i = geopar(gid,1)
if (geonm(gid,2).eq.'Y9') y9i = geopar(gid,1)
if (geonm(gid,2).eq.'GA') gai = geopar(gid,1)
if (geonm(gid,2).eq.'VG') vgi = geopar(gid,1)
if (geonm(gid,2).eq.'VF') vfi = geopar(gid,1)
if (geonm(gid,1).eq.'DLATEN') lai = geopar(gid,1)
if (geonm(gid,1).eq.'DLONEN') loi = geopar(gid,1)
if (geonm(gid,1).eq.'MT' ) mti = geopar(gid,1)
if (geonm(gid,1).eq.'MF' ) mfi = geopar(gid,1)
end do
*
do j=1,ldnj
do i=1,ldni
geobus(mti +(j-1)*ldni+i-1) = topo(i,j,1)
geobus(mfi +(j-1)*ldni+i-1) = topo(i,j,1)
end do
end do
*
call grid_to_grid_coef (xpx(gc_ld(1,myproc)+hx),ldni,
$ xpaq,nia,idx,cxa,cxb,cxc,cxd,inttyp)
call grid_to_grid_coef (ypx(gc_ld(3,myproc)+hy),ldnj,
$ ypaq,nja,idy,cya,cyb,cyc,cyd,inttyp)
*
call fstrhint (geobus(mgi),'MG',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(z0i),'Z0',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(lhi),'LH',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(y7i),'Y7',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(y8i),'Y8',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(y9i),'Y9',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(gai),'GA',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
*
call grid_to_grid_coef (xpx(gc_ld(1,myproc)+hx),ldni,
$ xpaq,nia,idx,cxa,cxb,cxc,cxd,'NEAREST')
call grid_to_grid_coef (ypx(gc_ld(3,myproc)+hy),ldnj,
$ ypaq,nja,idy,cya,cyb,cyc,cyd,'NEAREST')
*
call fstrhint (geobus(vgi),'VG',ldni,ldnj,1,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
call fstrhint (geobus(vfi),'VF',ldni,ldnj,26,nia,nja,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
*
call latlon2 (geobus(lai),geobus(loi),
$ xpx(gc_ld(1,myproc)+hx),ypx(gc_ld(3,myproc)+hy),ldni,ldnj)
*
do i=1,ldni*ldnj
geobus(mgi +i-1) = min(max(0.,geobus(mgi +i-1)),1.)
geobus(gai +i-1) = min(max(0.,geobus(gai +i-1)),1.)
geobus(lhi +i-1) = max(0.,geobus(lhi+i-1))
geobus(z0i +i-1) = max(1.0e-30,geobus(z0i+i-1))
end do
*
endif
*
deallocate (idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,xpaq,ypaq)
if (myproc.ne.0) ier = fstopc('MSGLVL','INFORM',.false.)
*
1001 format (/' INPUT DATA NOT ON SAME GRID ROTATION AS MODEL '
$ '- ABORT IN geodata -'/)
*-----------------------------------------------------------------------
return
end
*
subroutine fstrhint (f,nv,ni,nj,nk,nia,nja,idx,idy, 10,2
$ cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
implicit none
*
character* (*) nv,inttyp
integer ni,nj,nk,nia,nja,idx(*),idy(*)
real f(ni*nj,*)
real*8 cxa(*),cxb(*),cxc(*),cxd(*),cya(*),cyb(*),cyc(*),cyd(*)
*
#include "filename.cdk"
#include "ifd.cdk"
*
integer fnom,fstouv,fstfrm,fclos,longueur
character*256 fn
integer k,n,ier,unf,ofi,ofj
real wk1(nia*nja)
*
unf = 0
*
do k=1,nk
*
do n=1,ifd_nf
if (ifd_needit(n)) then
fn =
$ '../geophy/'//prefgeo(1:longueur(prefgeo))//'_'//ifd_fnext(n)
ier = fnom (unf,fn,'RND+OLD+R/O',0)
ier = fstouv (unf,'RND')
ofi = ifd_minx(n)-1
ofj = ifd_miny(n)-1
call filsfc (wk1,nv,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,k,nk,
$ unf,ofi,ofj)
ier = fstfrm (unf)
ier = fclos (unf)
endif
end do
*
call hinterpo (f(1,k),ni,nj,wk1,nia,nja,1,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,inttyp)
*
end do
*
return
end
*
subroutine filsfc ( f,nv,n1,n2,n3,n4,ki,nk,unf,ofi,ofj ) 1
implicit none
*
character* (*) nv
integer n1,n2,n3,n4,ki,nk,unf,ofi,ofj,cumerr
real f(n1:n2,n3:n4)
*
integer fstinf,fstluk
external fstinf,fstluk
*
character*8 dum
integer key,ni1,nj1,nk1,ip1,ier,i,j
real, dimension (:,:), allocatable :: tr1
*
*-----------------------------------------------------------------------
*
ip1 = 0
if (nk.gt.1) call convip ( ip1, real(ki), 3, 1, dum, .false.)
key = fstinf(unf,ni1,nj1,nk1,-1,' ',ip1,-1,-1,' ',nv)
allocate (tr1(ni1,nj1))
*
ier = fstluk ( tr1, key, ni1,nj1,nk1 )
*
do j=1,nj1
do i=1,ni1
f(ofi+i,ofj+j) = tr1(i,j)
end do
end do
*
deallocate (tr1)
*-----------------------------------------------------------------------
return
end