copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r casc_3df
*
subroutine casc_3df (geobus,trname,dimgx,dimgy,unf) 1,46
implicit none
*
character*8 trname(*)
integer dimgx,dimgy,unf
real geobus(*)
*
#include "dynmem.cdk"
#include "cdate.cdk"
#include "grd.cdk"
#include "rec.cdk"
#include "partopo.cdk"
#include "ifd.cdk"
#include "lesbus.cdk"
#include "vinterpo.cdk"
#include "hinterpo.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "bcsdim.cdk"
#include "bcsgrds.cdk"
#include "topo.cdk"
#include "nestpnt.cdk"
*
include 'mpif.h'
integer longueur,sid3df
external longueur,sid3df
character*8 dynophy,nomvar
character*8, dimension (:), allocatable :: trname_a,ponm
character*15 datev
character*256 fn
logical done,topo_local,dyn_done,phy_done
integer i,j,k,nia,nja,nka,ntra,ni1,nj1,nk1,ioerr,n,err,ng,
$ iscala(2),errop,iun1,ofi,ofj,l_i0,l_j0,l_in,l_jn,
$ iuv,itt,iwz,cnt,nkphy,errdyn,errphy,cumerr,pid,gid,nfphy,
$ niext,njext,id,ic,midx,midy,mnmtx,mnmty,mnmx,mnmy
integer, dimension (: ), allocatable :: ktop,kbot,idx,idu,idy
integer, dimension (:,:), allocatable :: popar
real xi,xf,yi,yf,htopa,maxtopo(2),maxtopo_g(2)
real , dimension (: ), allocatable :: phybr
real*8, dimension (: ), allocatable ::
$ xpaq,ypaq,xpau,ypav,xpuu,ypvv,
$ cxa,cxb,cxc,cxd,cua,cub,cuc,cud,cya,cyb,cyc,cyd
real, dimension (:,:), allocatable :: hgeow,hgeot,hgeom,
$ hgeow_anal,hgeot_anal,hgeom_anal,mxn,uun,vvn,wwn,
$ bbn,ppn,hun,phybn, wwr,bbr,ppr,hur
real, dimension (:,:,:), allocatable :: mxr,mxr_ext,trn,trr,w1,w2,
$ uur,vvr,topo_l,topo_h,topo_lg,topo_hg
real*8 xpxext(0:dimgx+1), ypxext(0:dimgy+1)
data nfphy,nkphy /0,0/
*-----------------------------------------------------------------------
*
do i=1,dimgx
xpxext(i) = xpx(i)
end do
xpxext(0) = xpxext(1) - (xpxext(2)-xpxext(1))
xpxext(dimgx+1) = xpxext(dimgx) + (xpxext(dimgx)-xpxext(dimgx-1))
*
do i=1,dimgy
ypxext(i) = ypx(i)
end do
ypxext(0) = ypxext(1) - (ypxext(2)-ypxext(1))
ypxext(dimgy+1) = ypxext(dimgy) + (ypxext(dimgy)-ypxext(dimgy-1))
*
* Read all needed files and construct the source domain for
* the horozontal interpolation
*
bcs_nia = ifd_niaf - ifd_niad + 1
bcs_nja = ifd_njaf - ifd_njad + 1
nia = bcs_nia
nja = bcs_nja
*
datev= gcrunstrt
done = .false.
iun1 = unf
*
ntra = 0
err = 0
*
do n=1,ifd_nf
if (ifd_needit(n)) then
errdyn = -1
errphy = -1
dyn_done = .false.
phy_done = .false.
fn ='../casc/3df_'//datev//'_'//ifd_fnext(n)
open (iun1,file=fn(1:longueur(fn)),access='SEQUENTIAL',
$ form='UNFORMATTED',status='OLD',iostat=errop)
if (errop.ne.0) goto 33
*
* Use first file to establish 3D grid dimensions and geo-references
* of all input stagerred grids (xpaq, ypaq, xpau and ypva).
*
if (.not.done)
$ allocate (xpaq(nia), ypaq(nja), xpau(nia), ypav(nja))
err = sid3df (xpaq,ypaq,xpau,ypav,unf,done,nia,nja,
$ bcs_nka,bcs_ntra)
nka = bcs_nka
ntra = bcs_ntra
if (.not.done) then
allocate ( mxn(nia*nja, 2), uun(nia*nja,nka ),
$ vvn(nia*nja,nka), wwn(nia*nja,nka ),
$ bbn(nia*nja,nka), ppn(nia*nja,nka+1),
$ hun(nia*nja,nka), trn(nia*nja,nka,ntra),
$ trname_a(ntra) )
trname_a = '!@@NOT@@'
endif
*
ofi = ifd_minx(n)-1
ofj = ifd_miny(n)-1
read (iun1,end=33) dynophy,cnt
*
23 if (dynophy.eq.'PHYSICSS') then
nfphy = cnt
if (.not.done) allocate (ponm(nfphy), popar(2,nfphy))
read (iun1,end=33) (ponm(i),popar(1,i),popar(2,i),i=1,nfphy)
if (.not.done) then
nkphy=0
do i=1,nfphy
nkphy = nkphy + max(popar(1,i),popar(2,i))
end do
allocate ( phybn(nia*nja,nkphy) )
endif
cumerr = 0
nkphy = 1
do i=1,nfphy
k = max(popar(1,i),popar(2,i))
call filmup (phybn(1,nkphy),ifd_niad,ifd_niaf,ifd_njad,
$ ifd_njaf,k,iun1,ofi,ofj,cumerr)
nkphy = nkphy + k
end do
nkphy = nkphy - 1
errphy = cumerr
phy_done = .true.
if (dyn_done) goto 33
read (iun1,end=33) dynophy,cnt
if (dynophy.ne.'DYNAMICS') goto 33
endif
*
cumerr=0
call filmup ( mxn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf, 2,
$ iun1,ofi,ofj,cumerr )
call filmup ( uun,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,cumerr )
call filmup ( vvn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,cumerr )
call filmup ( wwn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,cumerr )
call filmup ( bbn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,cumerr )
call filmup ( ppn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka+1,
$ iun1,ofi,ofj,cumerr )
call filmup ( hun,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,cumerr )
if (ntra.gt.0) then
call filuptr ( trn,ifd_niad,ifd_niaf,ifd_njad,ifd_njaf,nka,
$ iun1,ofi,ofj,trname,ntr,trname_a,ntra,cumerr )
endif
errdyn = cumerr
dyn_done = .true.
if ((gnmaphy.lt.1).or.(phy_done)) goto 33
read (iun1,end=33) dynophy,cnt
goto 23
*
33 done = .true.
*
if (gnmaphy.lt.1) errphy = 0
err = err + errdyn + errphy
if (err.lt.0) then
write (6,203) fn(1:longueur(fn)),myproc
goto 999
endif
close (iun1)
endif
end do
999 call mc2stop(err)
*
* Establish geo-references of model target horizontal grids
* (xp1, yp1, xpuu and ypvv).
l_i0 = 1 - west *hx
l_j0 = 1 - south*hy
l_in = ldni + east *hx
l_jn = ldnj + north*hy
ni1 = l_in - l_i0 + 1
nj1 = l_jn - l_j0 + 1
niext = ldni+2*hx+1
njext = ldnj+2*hy+1
allocate (mxr (l_i0:l_in ,l_j0:l_jn ,2),
$ mxr_ext(-hx:ldni+hx,-hy:ldnj+hy,2),
$ uur(ni1,nj1,nka+1),vvr(ni1,nj1,nka+1),wwr(ni1*nj1,nka),
$ bbr(ni1*nj1,nka ),ppr(ni1*nj1,nka+1),hur(ni1*nj1,nka),
$ trr(ni1*nj1,nka,ntra),phybr(ldni*ldnj*nkphy))
allocate (xpuu(l_i0:l_in),ypvv(l_j0:l_jn))
*
ofi = gc_ld(1,myproc) + hx - 1
ofj = gc_ld(3,myproc) + hy - 1
*
do i=l_i0,l_in
xpuu(i) = 0.5d0 * (xpxext(ofi+i-1)+xpxext(ofi+i))
end do
do j=l_j0,l_jn
ypvv(j) = 0.5d0 * (ypxext(ofj+j-1)+ypxext(ofj+j))
end do
*
* Horizontal interpolation (xpaq,ypaq) ===> (xp1,yp1)
*
allocate (idx(niext), idu(max(niext,njext)),idy(njext))
allocate (cxa(niext),cxb(niext),cxc(niext),cxd(niext),
$ cua(max(niext,njext)),cub(max(niext,njext)),
$ cuc(max(niext,njext)),cud(max(niext,njext)),
$ cya(njext),cyb(njext),cyc(njext),cyd(njext))
*
ofi = 1 - hx*west
ofj = 1 - hy*south
*
call grid_to_grid_coef (xpxext(gc_ld(1,myproc)-1),niext,xpaq,
$ nia,idx,cxa,cxb,cxc,cxd,hint_model)
call grid_to_grid_coef (ypxext(gc_ld(3,myproc)-1),njext,ypaq,
$ nja,idy,cya,cyb,cyc,cyd,hint_model)
call hinterpo (mxr_ext,niext,njext,mxn,nia,nja,2,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
*
call grid_to_grid_coef (xpxext(gc_ld(1,myproc)+(1- west)*hx),ni1,
$ xpaq,nia,idx,cxa,cxb,cxc,cxd,hint_model)
call grid_to_grid_coef (ypxext(gc_ld(3,myproc)+(1-south)*hy),nj1,
$ ypaq,nja,idy,cya,cyb,cyc,cyd,hint_model)
*
call hinterpo (wwr,ni1,nj1,wwn,nia,nja,nka,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
call hinterpo (bbr,ni1,nj1,bbn,nia,nja,nka,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
call hinterpo (ppr,ni1,nj1,ppn,nia,nja,nka+1,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
call hinterpo (hur,ni1,nj1,hun,nia,nja,nka,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
*
do k=1,ntra
if (trname_a(k).ne.'!@@NOT@@')
$ call hinterpo (trr(1,1,k),ni1,nj1,trn(1,1,k),nia,nja,nka,
$ idx,idy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint_model)
end do
*
* Horizontal interpolation (xpaq,ypaq) ===> (xp1,yp1) (physics grid)
*
i = 1 + west *hx
j = 1 + south*hy
call hinterpo (phybr,ldni,ldnj,phybn,nia,nja,nkphy,
$ idx(i),idy(j),cxa(i),cxb(i),cxc(i),cxd(i),
$ cya(j),cyb(j),cyc(j),cyd(j),hint_model)
*
* Horizontal interpolation (xpau,ypaq) ===> (xpuu,yp1)
*
call grid_to_grid_coef (xpuu,ni1,xpau,nia,idu,cua,cub,cuc,cud,
$ hint_model)
call hinterpo (uur(1,1,2),ni1,nj1,uun,nia,nja,nka,
$ idu,idy,cua,cub,cuc,cud,cya,cyb,cyc,cyd,hint_model)
*
* Horizontal interpolation (xpaq,ypav) ===> (xp1,ypvv)
*
call grid_to_grid_coef (ypvv,nj1,ypav,nja,idu,cua,cub,cuc,cud,
$ hint_model)
call hinterpo (vvr(1,1,2),ni1,nj1,vvn,nia,nja,nka,
$ idx,idu,cxa,cxb,cxc,cxd,cua,cub,cuc,cud,hint_model)
*
deallocate (xpaq,ypaq,xpau,ypav,xpuu,ypvv,
$ idx,idy,idu,cxa,cxb,cxc,cxd,cua,cub,cuc,cud,
$ cya,cyb,cyc,cyd,mxn,uun,vvn,wwn,bbn,ppn,hun)
if (nkphy.gt.0) deallocate (phybn)
*
do j=l_j0,l_jn
do i=l_i0,l_in
mxr(i,j,:) = mxr_ext(i,j,:)
end do
end do
*
maxtopo = 0.
do j=-hy,ldnj+hy
do i=-hx,ldni+hx
hh0i(i,j,1)= mxr_ext(i,j,1)
hh0i(i,j,2)= mxr_ext(i,j,2)
maxtopo(1) = max(maxtopo(1),hh0i(i,j,1)-hh0i(i,j,2))
maxtopo(2) = max(maxtopo(2),hh0i(i,j,2))
end do
end do
do j = 1-hy, ldnj+hy
do i = 1-hx, ldni+hx
hh0 (i,j,1) = hh0i(i,j,1)
hh0 (i,j,2) = hh0i(i,j,2)
end do
end do
call MPI_ALLREDUCE (maxtopo,maxtopo_g,2,
$ MPI_REAL,MPI_MAX,MPI_COMM_WORLD,err)
maxhh01_l = maxtopo_g(1)
maxhh02_l = maxtopo_g(2)
*
midx = gni/2+mod(gni,2)
midy = gnj/2+mod(gnj,2)
mnmtx = min(max(nesmt_bgx,0),midx-mod(gni,2))
mnmty = min(max(nesmt_bgy,0),midy-mod(gnj,2))
mnmx = min(max(nesmt_ndx,0),(gni-2*mnmtx-4)/2)
mnmy = min(max(nesmt_ndy,0),(gnj-2*mnmty-4)/2)
topo_local = .true.
if ((ldni.le.9+mnmtx+mnmx).or.(ldnj.le.9+mnmty+mnmy))
$ topo_local=.false.
*
if (topo_local) then
call nesmt_loc ( hh0f,hh0i,minx-1,maxx,miny-1,maxy,2,
$ mnmtx,mnmty,mnmx,mnmy )
else
if (myproc.eq.0) write(6,1001)
allocate (topo_l(dimgx+1,dimgy+1,2),topo_h(dimgx+1,dimgy+1,2),
$ topo_lg(dimgx+1,dimgy+1,2),topo_hg(dimgx+1,dimgy+1,2))
topo_l = 0.
topo_h = 0.
do j=-hy,ldnj+hy
do i=-hx,ldni+hx
topo_l(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,1)=hh0i(i,j,1)
topo_l(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,2)=hh0i(i,j,2)
topo_h(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,1)=hh0f(i,j,1)
topo_h(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,2)=hh0f(i,j,2)
end do
end do
call MPI_REDUCE (topo_l,topo_lg,(dimgx+1)*(dimgy+1) * 2,
$ MPI_INTEGER,MPI_BOR,0,MPI_COMM_WORLD,err)
call MPI_REDUCE (topo_h,topo_hg,(dimgx+1)*(dimgy+1) * 2,
$ MPI_INTEGER,MPI_BOR,0,MPI_COMM_WORLD,err)
if (myproc.eq.0) call nesmt (topo_hg,topo_lg,dimgx+1,dimgy+1,2,
$ hx,hy,nesmt_bgx,nesmt_bgy,nesmt_ndx,nesmt_ndy)
call MPI_bcast(topo_hg,(dimgx+1)*(dimgy+1) * 2,MPI_REAL,0,
$ MPI_COMM_WORLD,err)
*
do j=-hy,ldnj+hy
do i=-hx,ldni+hx
hh0f(i,j,1)=topo_hg(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,1)
hh0f(i,j,2)=topo_hg(gc_ld(1,myproc)+i+hx,gc_ld(3,myproc)+j+hy,2)
end do
end do
deallocate (topo_l,topo_h,topo_lg,topo_hg)
endif
*
uur(:,:,1) = uur(:,:,2)
vvr(:,:,1) = vvr(:,:,2)
*
ofj = 1
do pid=1,nfphy
do gid=1,geotop
nomvar = geonm(gid,2)
c next line must be commented for phy41
c if (geonm(gid,1).eq.'GLSEAEN') nomvar = 'GL'
if (ponm(pid).eq.nomvar) then
ofi = geopar(gid,1) - 1
if ((nomvar.eq.'LG').or.(nomvar.eq.'AL')
$ .or.(nomvar.eq.'HS')) then
do i=1,ldni*ldnj*geopar(gid,3)
geobus(ofi+i) = min(max(0.,phybr(ofj+i-1)),1.)
end do
else
do i=1,ldni*ldnj*geopar(gid,3)
geobus(ofi+i) = phybr(ofj+i-1)
end do
endif
goto 30
endif
end do
30 ofj = ofj + ldni*ldnj*popar(2,pid)
end do
if (nfphy.gt.0) deallocate (phybr, ponm, popar )
*
call statf_dm (hh0i,'MI',
$ 0, "geob", .false.,minx-1,maxx,miny-1,maxy,2,
$ -3,-3,1,gni+hx,gnj+hy,2)
call statf_dm (hh0f,'MF',
$ 0, "geob", .false.,minx-1,maxx,miny-1,maxy,2,
$ -3,-3,1,gni+hx,gnj+hy,2)
do id=1,geotop
do ic=1,geopar(id,3)
call statf_dm (geobus(geopar(id,1)+(ic-1)*ldni*ldnj),geonm(id,1),
$ ic, "geob", .false.,1,ldni,1,ldnj,1,
$ 1,1,1,gni,gnj,1)
end do
end do
*
* Preparation for vertical interpolation
*
ng = ni1*nj1
iuv = 1
itt = ng*gnk*2 + 1
iwz = ng*gnk*4 + 1
call hpalloc (paposit,ng*gnk*6 , err,1)
call hpalloc (pahuv ,ng*(nka+4) , err,1)
call hpalloc (pahtt ,ng*(nka+4) , err,1)
call hpalloc (pahww ,ng*(nka+4) , err,1)
*
allocate (hgeow(ng,gnk),hgeot(ng,gnk),hgeom(ng,0:gnk),
$ hgeow_anal(ng,nka),hgeot_anal(ng,nka),hgeom_anal(ng, nka+1))
*
* Compute target heights in hgeow, hgeot and hgeom
*
call hauteur (hgeow ,'HW',mxr,ng,1,gnk)
call hauteur (hgeot ,'HT',mxr,ng,1,gnk)
call hauteur (hgeom(1,1),'HM',mxr,ng,1,gnk)
cnt = 0
do j=l_j0,l_jn
do i=l_i0,l_in
cnt = cnt + 1
hgeom(cnt,0) = mxr(i,j,1)
end do
end do
*
* Compute reference (input data) heights in hgeow_anal,
* hgeot_anal and hgeom_anal.
*
call ref_h2 ( hgeow_anal,hgeot_anal,hgeom_anal,mxr,zta,ng,nka )
*
* Precompute vertical interpolation positions
*
allocate (ktop(ng*gnk), kbot(ng*gnk))
call posiz3 (posit(iuv),huv,hgeom_anal,hgeom(1,1),ktop,kbot,
$ ng,gnk,nka+1)
call posiz3 (posit(itt),htt,hgeot_anal,hgeot ,ktop,kbot,
$ ng,gnk,nka)
call posiz3 (posit(iwz),hww,hgeow_anal,hgeow ,ktop,kbot,
$ ng,gnk,nka)
deallocate ( ktop,kbot,hgeow,hgeow_anal, hgeom_anal )
*
* Perform vertical interpolations and transfer into model space
*
allocate ( w1(l_i0:l_in,l_j0:l_jn,gnk+1),
$ w2(l_i0:l_in,l_j0:l_jn,gnk+1) )
call vertint3 ( w1,uur,posit(iuv),huv,ng,gnk,nka+1)
call vertint3 ( w2,vvr,posit(iuv),huv,ng,gnk,nka+1)
do k=1,gnk-1
do j=l_j0,l_jn
do i=l_i0,l_in
uup(i,j,k) = w1(i,j,k) / sqrt(sby(i,j))
vvp(i,j,k) = w2(i,j,k) / sqrt(sbx(i,j))
end do
end do
end do
call vertint3 (w1,wwr,posit(iwz),hww,ng,gnk,nka )
call vertint3 (w2,hur,posit(itt),htt,ng,gnk,nka )
do k=1,gnk
do j=l_j0,l_jn
do i=l_i0,l_in
wwp(i,j,k) = w1(i,j,k)
hup(i,j,k) = w2(i,j,k)
end do
end do
end do
call vertint3 (w1,bbr,posit(itt),htt,ng,gnk,nka )
call vertint3 (w2,ppr,posit(iuv),huv,ng,gnk,nka+1)
call corbusfc (w1,hgeot,hgeom,hgeot_anal,ng,gnk)
do k=1,gnk
do j=l_j0,l_jn
do i=l_i0,l_in
bbp(i,j,k) = w1(i,j,k)
ppp(i,j,k) = w2(i,j,k)
end do
end do
end do
*
do n=1,ntr
j=-1
do k=1,ntra
if (trname(n).eq.trname_a(k)) j=k
end do
if (j.gt.0) then
call vertint3 (w1,trr(1,1,j),posit(itt),htt,ng,gnk,nka)
do k=1,gnk
do j=l_j0,l_jn
do i=l_i0,l_in
trp(i,j,k,n) = w1(i,j,k)
end do
end do
end do
else
trp(:,:,:,n) = 0.
endif
end do
*
cnt = 0
do j=l_j0,l_jn
do i=l_i0,l_in
cnt = cnt + 1
uup(i,j,gnk) = uup(i,j,gnk-1)
vvp(i,j,gnk) = vvp(i,j,gnk-1)
ppp(i,j,0 ) = ppr(cnt,1)
end do
end do
*
deallocate (mxr,mxr_ext,uur,vvr,wwr,bbr,ppr,hur,trr,
$ hgeot,hgeom,hgeot_anal,w1,w2)
call hpdeallc (paposit,err,1)
call hpdeallc (pahuv ,err,1)
call hpdeallc (pahtt ,err,1)
call hpdeallc (pahww ,err,1)
paposit = 0
pahuv = 0
pahtt = 0
pahww = 0
*
if (myproc.eq.0) then
write(6,100)
write(6,101) datev
write(6,100)
endif
*
100 format (' ',65('*'))
101 format (' (CASC_3DF) JUST READ INPUT DATA FOR DATE: ',a15)
201 format (/' CAN NOT OPEN FILE: ',a,', PROC#:',i4,' --ABORT--'/)
202 format (/' INSUFFICIENT INPUT DATA COVERAGE: ',4L,', PROC#:',i4,' --ABORT--'/)
203 format (/' PROBLEM WITH FILE: ',a,', PROC#:',i4,' --ABORT--'/)
303 format (a,4e15.7,2i10)
1001 format (/' Global treatment of topography - local domain too small')
*-----------------------------------------------------------------------
return
end
*
subroutine filmup ( f,n1,n2,n3,n4,nk,unf,ofi,ofj,cumerr ) 15
implicit none
*
integer n1,n2,n3,n4,nk,unf,ofi,ofj,cumerr
real f(n1:n2,n3:n4,nk)
*
character*4 nomvar
integer i,j,k,ni1,nj1,nk1,err,n,nbits,nb
real, dimension (:), allocatable :: wkc
real, dimension (:,:,:), allocatable :: tr1
*
*-----------------------------------------------------------------------
nb = 0
err = -1
read (unf,end=44) nomvar,ni1,nj1,nk1,nbits
allocate (tr1(ni1,nj1,nk1))
if (nbits.ge.32) then
read (unf,end=45) tr1
else
n = (ni1*nj1*nbits+120+32-1)/32
allocate (wkc(n))
do k=1,nk1
read (unf) wkc
call xxpak (tr1(1,1,k), wkc, ni1, nj1, -nbits, nb, 2)
end do
deallocate (wkc)
endif
do k=1,nk1
do j=1,nj1
do i=1,ni1
f(ofi+i,ofj+j,k) = tr1(i,j,k)
end do
end do
end do
err = 0
45 deallocate (tr1)
44 cumerr = cumerr + err
*-----------------------------------------------------------------------
return
end
*
subroutine filuptr ( f,n1,n2,n3,n4,nk,unf,ofi,ofj,trname,ntr, 2
$ trname_a,ntra,cumerr )
implicit none
*
integer n1,n2,n3,n4,nk,unf,ofi,ofj,ntr,ntra,cumerr
character*8 trname(ntr),trname_a(ntra)
real f(n1:n2,n3:n4,nk,*)
*
character*4 nomvar
integer i,j,k,n,m,ni1,nj1,nk1,takeit,err,nw,nbits,nb
real, dimension (:), allocatable :: wkc
real, dimension (:,:,:), allocatable :: tr1
*-----------------------------------------------------------------------
nb = 0
do n=1,ntra
err = -1
read (unf,end=54) nomvar,ni1,nj1,nk1,nbits
if (n.eq.1) allocate (tr1(ni1,nj1,nk1))
takeit=-1
do m=1,ntr
if (trname(m)(1:4).eq.nomvar) takeit=m
end do
if (takeit.gt.0) then
trname_a(n) = trname(takeit)
if (nbits.ge.32) then
read (unf,end=55) tr1
else
nw = (ni1*nj1*nbits+120+32-1)/32
allocate (wkc(nw))
do k=1,nk1
read (unf) wkc
call xxpak (tr1(1,1,k), wkc, ni1, nj1, -nbits, nb, 2)
end do
deallocate (wkc)
endif
do k=1,nk1
do j=1,nj1
do i=1,ni1
f(ofi+i,ofj+j,k,n) = tr1(i,j,k)
end do
end do
end do
err = 0
else
trname_a(n) = '!@@NOT@@'
read (unf,end=55)
err = 0
endif
end do
55 deallocate (tr1)
54 cumerr = cumerr + err
*-----------------------------------------------------------------------
return
end
*
subroutine hinterpo ( dst,nid,njd, src,nis,njs, nk, indxx,indxy, 18,1
$ cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint)
implicit none
*
character* (*) hint
integer nid,njd,nis,njs,nk
integer indxx(nid),indxy(njd)
real dst(nid,njd,nk), src(nis,njs,nk),
$ cxa(*),cxb(*),cxc(*),cxd(*),cya(*),cyb(*),cyc(*),cyd(*)
*
integer k,err,ezsint
external ezsint
*-----------------------------------------------------------------------
do k=1,nk
call grid_to_grid_interp (dst(1,1,k),nid,njd,src(1,1,k),
$ nis,njs,indxx,indxy,cxa,cxb,cxc,cxd,cya,cyb,cyc,cyd,hint)
end do
*-----------------------------------------------------------------------
return
end
*
subroutine corbusfc (bu,hgeot,hgeom,hgeot_anal,ng,lnk) 2,1
implicit none
*
integer ng,lnk
real bu(ng,lnk),hgeot(ng,lnk),hgeom(ng,lnk),hgeot_anal(ng)
*
#include "consdyn_8.cdk"
integer i,k
real , dimension (:,:), allocatable :: gotsr,ortsr,qtsr,nstr2
real*8 rtt,zh
*-----------------------------------------------------------------------
*
allocate ( gotsr(ng,lnk ), ortsr(ng,0:lnk),
$ qtsr (ng,0:lnk), nstr2(ng, lnk) )
call qntstar (qtsr,nstr2,gotsr,ortsr,hgeot,hgeom,ng,0,lnk)
do k=1,lnk
do i=1,ng
zh = hgeot_anal(i)-hgeot(i,k)
if (zh.gt.0.d0) then
rtt = (bu(i,k)+grav_8)/gotsr(i,k) - tcdk_8
rtt = rtt + 0.0065*zh
bu(i,k) = (rtt+tcdk_8)*gotsr(i,k) - grav_8
endif
end do
end do
deallocate ( gotsr, ortsr, qtsr, nstr2 )
*
*-----------------------------------------------------------------------
return
end
*