copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
*** s/r dc_topo2
*
subroutine dc_topo2 (hh0,meh,coef,minx,maxx,miny,maxy,lni,lnj,hx, 2,2
$ hy,period_x,period_y,maxhh01,maxhh02,prout)
implicit none
*
logical period_x,period_y,prout
integer minx,maxx,miny,maxy,lni,lnj,hx,hy
real hh0(minx-1:maxx,miny-1:maxy,2),meh(*),coef,maxhh01,maxhh02
*
*AUTHOR
* Daniel Leuenberger Apr 2001
*
*OBJECT
* decomposes a given topography field: INPUT/OUTPUT hh0(,,1)
* in a large-scale part: hh0(,,1) - hh0(,,2)
* and a small-scale part: OUTPUT hh0(,,2)
*
* a standard laplace filter (1st order) is applied.
* the boundary values are treated seperately to assure, that the
* also these points are smoothed:
* i.e. at the i=1 boundary: A(0,j) = A(2,j) for all j
* i=ni boundary: A(ni+1,j) = A(ni-1,j) for all j
* j=1 boundary: A(i,0) = A(i,2) for all i
* j=nj boundary: A(i,nj+1) = A(i,nj-1) for all i
* nflt determines, how often the filter is applied
* (nflt has to be an even number)
* Additionally, the maxima of h0, hh0(,,1) and hh0(,,2) are
* computed and written to standard output.
*
#include "partopo.cdk"
*
include 'mpif.h'
integer i,j,n,nflt,err,cnt
real t1(minx-1:maxx+1,miny-1:maxy+1,2),
$ w1(minx-1:maxx+1,miny-1:maxy+1),maxtopo(2),maxtopo_g(2)
real*8 con1
**
*----------------------------------------------------------------------
*
if (coef.gt.0.) then
*
con1 = 1. - coef
*
cnt = 0
do j=-hy,lnj+hy
do i=-hx,lni+hx
cnt = cnt+1
t1 (i,j,1) = meh(cnt)
end do
end do
*
call rpn_comm_xch_halo (t1,minx-1,maxx+1,miny-1,maxy+1,
$ lni,lnj,1,hx+1,hy+1,period_x,period_y,lni,0)
*
do j=-hy,lnj+hy
do i=1-hx,lni+hx-1
w1(i,j) = coef * (t1(i-1,j,1) + t1(i+1,j,1)) / 2.
$ + con1 * t1 (i,j,1)
enddo
if (west_L)
$ w1(-hx,j) = coef * ( t1(-hx,j,1) + t1(1-hx,j,1) ) / 2.
$ + con1 * t1 (-hx,j,1)
if (east_L)
$ w1(lni+hx,j) = coef * (t1(lni+hx-1,j,1)+t1(lni+hx,j,1)) / 2.
$ + con1 * t1 (lni+hx,j,1)
end do
*
do j=1-hy,lnj+hy-1
do i=-hx,lni+hx
t1 (i,j,1) = coef * (w1(i,j-1)+w1(i,j+1)) / 2.
$ + con1 * w1(i,j)
end do
end do
if (south_L) then
do i=-hx,lni+hx
t1 (i,-hy,1) = coef * (w1(i,-hy)+w1(i,1-hy)) / 2.
$ + con1 * w1(i,-hy)
end do
endif
if (north_L) then
do i=-hx,lni+hx
t1 (i,lnj+hy,1) = coef*(w1(i,lnj+hy-1)+w1(i,lnj+hy)) / 2.
$ + con1*w1(i,lnj+hy)
end do
endif
*
call rpn_comm_xch_halo (t1,minx-1,maxx+1,miny-1,maxy+1,
$ lni,lnj,1,hx+1,hy+1,period_x,period_y,lni,0)
*
do j=-hy,lnj+hy
do i=-hx,lni+hx
t1 (i,j,1) = max(0.,t1(i,j,1))
hh0(i,j,1) = t1(i,j,1)
end do
end do
*
else
*
cnt = 0
do j=-hy,lnj+hy
do i=-hx,lni+hx
cnt = cnt+1
t1 (i,j,1) = max(0.,meh(cnt))
hh0(i,j,1) = t1(i,j,1)
end do
end do
*
endif
*
call rpn_comm_xch_halo (t1,minx-1,maxx+1,miny-1,maxy+1,
$ lni,lnj,1,hx+1,hy+1,period_x,period_y,lni,0)
nflt = 100
*
* apply nflt times an ideal 2d filter to compute the
* longwave part hh0(,,1) of the topography
*
do n=1,nflt/2
*
call rpn_comm_xch_halo (t1,minx-1,maxx+1,miny-1,maxy+1,
$ lni,lnj,2,hx+1,hy+1,period_x,period_y,lni,0)
call dc_topop (t1(minx-1,miny-1,2),t1(minx-1,miny-1,1),
$ minx-1,maxx+1,miny-1,maxy+1,lni,lnj,hx,hy)
call dc_topop (t1(minx-1,miny-1,1),t1(minx-1,miny-1,2),
$ minx-1,maxx+1,miny-1,maxy+1,lni,lnj,hx,hy)
*
end do
*
call rpn_comm_xch_halo (t1,minx-1,maxx+1,miny-1,maxy+1,
$ lni,lnj,2,hx+1,hy+1,period_x,period_y,lni,0)
*
* compute the shortwave part hh0(,,2) of the topo with the relation
* hh0(,,2) = h0(,) - hh0(,,1)
*
maxtopo = 0.
do j=-hy,lnj+hy
do i=-hx,lni+hx
hh0(i,j,2) = hh0(i,j,1) - t1(i,j,1)
maxtopo(1) = max(maxtopo(1), t1(i,j,1))
maxtopo(2) = max(maxtopo(2),hh0(i,j,2))
end do
end do
call MPI_ALLREDUCE (maxtopo,maxtopo_g,2,
$ MPI_REAL,MPI_MAX,MPI_COMM_WORLD,err)
maxhh01 = maxtopo_g(1)
maxhh02 = maxtopo_g(2)
*
if (prout) then
print *, 'MAXIMA OF TOPOGRAPHY:'
print *, 'MAX TOPO(large scale): ',maxhh01
print *, 'MAX TOPO(small scale): ',maxhh02
endif
*
*----------------------------------------------------------------------
return
end
*
subroutine dc_topop (d,s,lminx,lmaxx,lminy,lmaxy,ldni,ldnj,hx,hy) 2
implicit none
*
integer lminx,lmaxx,lminy,lmaxy,ldni,ldnj,hx,hy
real d(lminx:lmaxx,lminy:lmaxy), s(lminx:lmaxx,lminy:lmaxy)
*
#include "partopo.cdk"
*
integer i,j
**
*----------------------------------------------------------------------
*
* inner points
*
do j=1-hy,ldnj+hy-1
do i=1-hx,ldni+hx-1
d(i,j) = s(i,j)/2.+(s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i,j+1))/8.
end do
end do
*
* corner points
*
if (south_L.and.west_L)
$ d(-hx,-hy) = s(-hx,-hy)/2. + s(1-hx,-hy)/4. + s(-hx,1-hy)/4.
if (north_L.and.west_L)
$ d(-hx,ldnj+hy) = s(-hx,ldnj+hy)/2. + s(1-hx,ldnj+hy)/4. + s(-hx,ldnj+hy-1)/4.
if (south_L.and.east_L)
$ d(ldni+hx,-hy) = s(ldni+hx,-hy)/2. + s(ldni+hx-1,-hy)/4. + s(ldni+hx,1-hy)/4.
if (north_L.and.east_L)
$ d(ldni+hx,ldnj+hy) = s(ldni+hx,ldnj+hy)/2. + s(ldni+hx-1,ldnj+hy)/4. + s(ldni+hx,ldnj+hy-1)/4.
*
* edge points
*
if (west_L) then
do j=1-hy,ldnj+hy-1
d(-hx,j) = s(-hx,j)/2. + s(1-hx,j)/4. + (s(-hx,j-1) + s(-hx,j+1))/8.
end do
endif
if (east_L) then
do j=1-hy,ldnj+hy-1
d(ldni+hx,j)= s(ldni+hx,j)/2.+s(ldni+hx-1,j)/4.+(s(ldni+hx,j-1)+s(ldni+hx,j+1))/8.
end do
endif
*
if (south_L) then
do i=1-hx,ldni+hx-1
d(i,-hy) = s(i,-hy)/2. + s(i,1-hy)/4. + (s(i-1,-hy) + s(i+1,-hy))/8.
end do
endif
if (north_L) then
do i=1-hx,ldni+hx-1
d(i,ldnj+hy) = s(i,ldnj+hy)/2.+s(i,ldnj+hy-1)/4.+(s(i-1,ldnj+hy)+s(i+1,ldnj+hy))/8.
end do
endif
*
*----------------------------------------------------------------------
return
end