copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r e_geopfld
subroutine e_geopfld ( topo_h) 1,19
implicit none
*
real topo_h(*)
*
*AUTHOR M. Desgagne May 2001
*
*IMPLICIT
#include "lesbus.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "grd.cdk"
#include "filename.cdk"
#include "lcldim.cdk"
#include "cdate.cdk"
#include "rec.cdk"
#include "hinterpo.cdk"
#include "nesting.cdk"
#include "lun.cdk"
*
**
integer fnom,fstouv,fclos,fstfrm,s_rdhint,longueur,pilotf,
$ newdate,vrtstc
external fnom,fstouv,fclos,fstfrm,s_rdhint,longueur,pilotf,
$ newdate,vrtstc
*
character*2 typvar
character*4 nomvar
character*8 etk,typ,interp,seq,dum
logical flag
integer ip1,ip2,ip3,iun2,iun3,unf,iun,
$ i,j,indx,ier,im,ivar,ic,offk,offj,idatev,ni,nj,
$ jj, month,day,mi,mf,fstdate,dat2,dat3
integer , dimension (:) , allocatable :: err
real f2d(Grd_ni,Grd_nj),wk2(Grd_ni,Grd_nj),poi,pof
real , dimension (:) , allocatable :: geobus
real*8 dayfrac,sec_in_day
parameter (sec_in_day=86400.0d0)
*
integer nkfin
parameter (nkfin = 1000)
real topo_l(Grd_ni*Grd_nj*2)
integer nilvls,err_mx,err_im,err_v1,err_v2,n,refip1(nkfin*3)
*---------------------------------------------------------------------
*
un_geo = 33
open(un_geo,file='process/geophy.bin',access='SEQUENTIAL',
$ form='UNFORMATTED')
write (un_geo) gni,gnj,hx,hy
call wpilpar (un_geo)
*
if (sfc2d) then
ni = gni-2*hx
nj = gnj-2*hy
etk = ' '
typvar = ' '
*
* Interpolating Target Topography (ME)
*
call e_gettopo_h ( topo_h )
*
* Open climato and geophy files
*
iun2 = 0
iun3 = 0
ier = fnom (iun2,climato ,'RND+OLD+R/O',0)
ier = fstouv (iun2,'RND')
if (ier.lt.0) then
write (6,1002) climato(1:longueur(climato))
stop
endif
ier = fnom (iun3,geophy ,'RND+OLD+R/O',0)
ier = fstouv (iun3,'RND')
if (ier.lt.0) then
write (6,1002) geophy(1:longueur(geophy))
stop
endif
*
* Open proper input analysis file for sfc fields
*
call datp2f (idatev,gcrunstrt)
ier = pilotf (idatev,'UU',etk,typvar,-1,-1,-1)
flag = .false.
if (ier.lt.0) then
write (6,900)
flag = .true.
goto 8877
endif
*
* Interpolating Topography (MX) of the analysis
*
nilvls = vrtstc (refip1,idatev,nkfin,un_pil,.false.)
if (nilvls.le.0) then
if (pilotf (idatev,'UU',' ',' ',-1,-1,-1).lt.0) then
write (6,1006) gcrunstrt
stop
endif
nilvls = vrtstc (refip1,idatev,nkfin,un_pil,.false.)
endif
*
if (gngalsig.eq.2) then
err_mx = s_rdhint (topo_l,xpx,ypx,Grd_ni,Grd_nj,'GZ',-1,
$ refip1(1),-1,-1,etk,typvar,.false.,hint_ntr,un_pil,6)
if (err_mx.ge.0) then
do i=1,Grd_ni*Grd_nj
topo_l(i) = topo_l(i) * 10.
end do
call dc_topo (topo_l,maxhh01_l,maxhh02_l,Grd_ni,Grd_nj )
else
stop
endif
endif
if (gngalsig.eq.0) then
do i=1,Grd_ni*Grd_nj*2
topo_l(i) = topo_h(i)
end do
err_mx = 0
endif
if (gngalsig.eq.1) then
call convip (ip1,1.,3,2,dum,.false.)
err_mx = s_rdhint (topo_l,xpx,ypx,Grd_ni,Grd_nj,'MX',idatev,
$ ip1,-1,-1,etk,typvar,.false.,hint_ntr,un_pil,6)
if (err_mx.ge.0) then
err_mx = -1
call convip (ip1,2.,3,2,dum,.false.)
err_mx = s_rdhint (topo_l(Grd_ni*Grd_nj+1),xpx,ypx,
$ Grd_ni,Grd_nj,'MX',idatev,ip1,-1,-1,
$ etk,typvar,.false.,hint_ntr,un_pil,6)
endif
maxhh01_l = 0.
maxhh02_l = 0.
n = Grd_ni*Grd_nj
do i=1,n
maxhh01_l = max(maxhh01_l,topo_l(i)-topo_l(n+i))
maxhh02_l = max(maxhh02_l,topo_l(n+i))
end do
endif
*
* Filling geobus
*
8877 print*
call geopini (ni,nj,.true.)
allocate (geobus(max(1,geospc)),err(geotop))
geobus = 0.
err = -1
*
call datp2f ( fstdate,gcrunstrt )
ier = newdate ( fstdate,dat2,dat3,-3 )
dat3 = dat2 / 100
day = dat2 - dat3*100
dat2 = dat3 / 100
month = dat3 - dat2*100
*
if ( day .ge. 15 ) then
jj = day
mi = month
else
jj = day + 30
mi = month - 1
if ( mi .lt. 1 ) mi = 12
endif
mf = mi + 1
if ( mf .gt. 12 ) mf = 1
pof = (real(jj-15)/30.)
poi = 1. - pof
*
do 100 ivar=1,geotop
870 nomvar = geonm(ivar,2)
if (nomvar.eq.'00') then
err(ivar) = 0
goto 100
endif
etk = geonm(ivar,3)
typ = geonm(ivar,4)
interp= geonm(ivar,5)
seq(1:8) = ' '
seq(1:longueur(geonm(ivar,6))) = geonm(ivar,6)
write(6,1001) nomvar,geonm(ivar,1),etk,typ,interp,seq
670 if (seq(1:1).eq." ") goto 100
ip1 = -1
ip2 = -1
unf = un_pil
if ((seq(1:1).eq."C").or.(seq(1:1).eq."V")) then
unf = iun2
ip2 = month
endif
if (seq(1:1).eq."G") unf = iun3
do 110 im=1,geopar(ivar,3)
if (geopar(ivar,3).gt.1) then
if ((seq(1:1).eq."A").and.(gngalsig.eq.1)) then
call convip ( ip1, real(im), 3, 2, dum, .false. )
else
call convip ( ip1, real(im), 3, 1, dum, .false. )
endif
endif
indx = geopar(ivar,1) + (im-1)*ni*nj
err_im = -1
if (seq(1:1).ne."V") then
if (unf.gt.0)
$ err_im= s_rdhint (f2d,xpx,ypx,Grd_ni,Grd_nj,nomvar,-1,
$ ip1,ip2,-1,etk,typ,.false.,interp,unf,6)
else
write (6,1000) day, month,poi, mi, pof, mf
if (unf.gt.0) then
err_v1= s_rdhint (f2d,xpx,ypx,Grd_ni,Grd_nj,nomvar,-1,
$ ip1,mi,-1,etk,typ,.false.,interp,unf,6)
err_v2= s_rdhint (wk2,xpx,ypx,Grd_ni,Grd_nj,nomvar,-1,
$ ip1,mf,-1,etk,typ,.false.,interp,unf,6)
if ((err_v1.ge.0).and.(err_v2.ge.0)) then
err_im = 0
do j=1,Grd_nj
do i=1,Grd_ni
f2d(i,j) = poi*f2d(i,j) + pof*wk2(i,j)
end do
end do
endif
endif
endif
if (err_im.lt.0) then
if ((nomvar(1:2).eq.'LG').and.(seq(2:2).eq.' '))then
geonm(ivar,2) = 'GL'
goto 870
else if (nomvar(1:2).eq.'Z0') then
geonm(ivar,2) = 'ZP'
goto 870
else
if (seq(2:2).eq.' ') then
write(6,1004) nomvar,seq(1:1)
else
write(6,1003) nomvar,seq(1:1),seq(2:2)
endif
seq = seq(2:longueur(seq))
goto 670
endif
else
if (nomvar(1:2).eq.'ZP') then
do j=1,Grd_nj
do i=1,Grd_ni
f2d(i,j) = exp(f2d(i,j))
end do
end do
endif
err(ivar) = 0
endif
if (err(ivar).ge.0) then
offk = geopar(ivar,1)+(ni*nj)*(im-1)
do j=1,nj
offj=ni*(j-1)
do i=1,ni
geobus(offk+offj+i-1)= f2d(i+hx+1,j+hy+1)
end do
end do
endif
*
110 continue
100 continue
*
if (.not.flag) then
flag=.false.
do ivar = 1, geotop
if (err(ivar).lt.0) then
flag=.true.
write (6,1005) geonm(ivar,2)
endif
end do
else
flag=.false.
endif
if (flag) stop
*
if (gnmaphy.gt.0) call geopfld (geobus,geospc,topo_h,
$ ni,nj,Grd_ni,Grd_nj)
*
do i =1,geotop
do ic=1,geopar(i,3)
call statfld ( geobus(geopar(i,1)+(ic-1)*ni*nj), geonm(i,1),
$ ic, "geopfld", .false.,1,ni,1,nj,1,
$ 1,1,1,ni,nj,1)
end do
end do
*
* Write geobus and topography ME,MX on file ungeo
*
call wrgeop4 ( topo_l,topo_h,geobus,geospc,Grd_ni,Grd_nj )
*
deallocate (geobus,err)
ier = fstfrm (iun2)
ier = fclos (iun2)
ier = fstfrm (iun3)
ier = fclos (iun3)
*
endif
*
close (un_geo)
*
write (6,600)
*
600 format (' ',' rdgeoanl complete')
900 format (' TOPOGRAPHY MX NOT PROCESSED')
1000 format ('AVERAGING CLIMATOLOGY for (dd/mm): ',i2,'/',i2,
$ ' (month/weight): ',2('(',f7.5,'/',i2,') '))
1001 format (/'====>',' Processing ',a2,' ',a8,4a8,' <====')
1002 format ( ' FILE: ',a,' NOT AVAILABLE --- ABORT ---'/)
1003 format ('Could not find ',a,' in ',a,' - switching to ',a)
1004 format ('Could not find ',a,' in ',a,' - quitting')
1005 format (' GEOPHYSICAL FIELD: 'a,' NOT AVAILABLE --- ABORT ---')
1006 format (/' MISSING DATA VALID AT: 'a/,10x,'----- ABORT -----'/)
*
*---------------------------------------------------------------------
*
return
end
*