!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRADT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculating the gradient of costfunction.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! Jidong Gao, CAPS, July, 2000
!
!-----------------------------------------------------------------------
!
SUBROUTINE gradt(numctr,ctrv,grad, & 2,22
gdu_err,gdv_err,gdp_err,gdt_err,gdq_err,gdw_err, &
u_ctr,v_ctr,p_ctr,t_ctr,q_ctr,w_ctr, psi, phi, &
gdscal, &
nx,ny,nz, &
nvar,nvarradin,nvarrad,nzua,nzrdr,nzret, &
mapfct,j1,j2,j3,aj3x,aj3y,aj3z,j3inv,rhostr, &
rhostru, rhostrv, rhostrw, div3, &
mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret, &
nsrcsng,nsrcua,nsrcrad,nsrcret,ncat, &
mxpass,npass,iwstat,xs,ys,zs,x,y,z,zp,hterain, &
icatg,xcor,nam_var,xsng,ysng,hgtsng,thesng, &
obsng,odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng, &
xua,yua,hgtua,theua, &
obsua,odifua,qobsua,qualua,isrcua,nlevsua,nobsua, &
elvrad,xradc,yradc, &
distrad,uazmrad,vazmrad,hgtradc,theradc,dsdr,dhdr, &
obsrad,odifrad,qobsrad,qualrad, &
irad,isrcrad,nlevrad,ncolrad, &
xretc,yretc,hgtretc,theretc, &
obsret,odifret,qobsret,qualret, &
iret,isrcret,nlevret,ncolret, &
srcsng,srcua,srcrad,srcret, &
ianxtyp,iusesng,iuseua,iuserad,iuseret, &
xyrange,kpvrsq,wlim,zrange,zwlim, &
thrng,rngsqi,knt,wgtsum,zsum, &
corsng,corua,corrad,corret, &
xsng_p,ysng_p,ihgtsng,xua_p,yua_p,ihgtua, &
xradc_p,yradc_p,ihgtradc,zsng_1,zsng_2, &
zua_1,zua_2,zradc_1,zradc_2, &
oanxsng,oanxua,oanxrad,oanxret, &
sngsw,uasw,radsw,retsw, &
ipass_filt,hradius,nradius_z,turn_div,cntl_var, &
wei_div_h,wei_div_v, &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INCLUDE 'varpara.inc'
!
INCLUDE 'grid.inc'
!-----------------------------------------------------------------------
!
! Input Sizing Arguments
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz,sngsw,uasw,radsw,retsw
INTEGER :: nvar,nvarradin,nvarrad
INTEGER :: nzua,nzrdr,nzret
INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
INTEGER :: mxpass,npass
INTEGER :: ipass_filt,hradius,nradius_z,turn_div,cntl_var
REAL :: wei_div_h, wei_div_v
!
!-----------------------------------------------------------------------
!
! input grid arguments
!
!-----------------------------------------------------------------------
!
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z (nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL :: hterain(nx,ny) ! The height of the terrain.
!
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( x ).
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( y ).
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as d( zp )/d( z ).
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: j3inv (nx,ny,nz) ! Inverse of j3
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
!
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zs(nx,ny,nz)
INTEGER :: icatg(nx,ny)
REAL :: xcor(ncat,ncat)
!
!-----------------------------------------------------------------------
!
! Input Observation Arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=6) :: nam_var(nvar)
REAL :: xsng(mxsng)
REAL :: ysng(mxsng)
REAL :: hgtsng(mxsng)
REAL :: thesng(mxsng)
REAL :: obsng(nvar,mxsng)
REAL :: odifsng(nvar,mxsng)
REAL :: qobsng(nvar,mxsng)
INTEGER :: qualsng(nvar,mxsng)
INTEGER :: isrcsng(mxsng)
INTEGER :: icatsng(mxsng)
INTEGER :: nobsng
REAL :: xsng_p(mxsng),ysng_p(mxsng)
REAL :: zsng_1(mxsng),zsng_2(mxsng)
INTEGER :: ihgtsng(mxsng)
!
REAL :: xua(mxua)
REAL :: yua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: theua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
REAL :: odifua(nvar,nzua,mxua)
REAL :: qobsua(nvar,nzua,mxua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: nlevsua(mxua)
INTEGER :: isrcua(mxua)
INTEGER :: nobsua
!
REAL :: xua_p(mxua),yua_p(mxua)
REAL :: zua_1(nzua,mxua),zua_2(nzua,mxua)
INTEGER :: ihgtua(nzua,mxua)
!
REAL :: elvrad(mxrad)
REAL :: xradc(mxcolrad)
REAL :: yradc(mxcolrad)
REAL :: distrad(mxcolrad)
REAL :: uazmrad(mxcolrad)
REAL :: vazmrad(mxcolrad)
REAL :: hgtradc(nzrdr,mxcolrad)
REAL :: theradc(nzrdr,mxcolrad)
REAL :: dsdr(nzrdr,mxcolrad)
REAL :: dhdr(nzrdr,mxcolrad)
REAL :: obsrad(nvarradin,nzrdr,mxcolrad)
REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
REAL :: qobsrad(nvarrad,nzrdr,mxcolrad)
INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
INTEGER :: nlevrad(mxcolrad)
INTEGER :: irad(mxcolrad)
INTEGER :: isrcrad(0:mxrad)
INTEGER :: ncolrad
!
REAL :: xradc_p(mxcolrad),yradc_p(mxcolrad)
REAL :: zradc_1(nzrdr,mxcolrad),zradc_2(nzrdr,mxcolrad)
INTEGER :: ihgtradc(nzrdr,mxcolrad)
!
REAL :: xretc(mxcolret)
REAL :: yretc(mxcolret)
REAL :: hgtretc(nzret,mxcolret)
REAL :: theretc(nzret,mxcolret)
REAL :: obsret(nvar,nzret,mxcolret)
REAL :: odifret(nvar,nzret,mxcolret)
REAL :: qobsret(nvar,nzret,mxcolret)
INTEGER :: qualret(nvar,nzret,mxcolret)
INTEGER :: nlevret(mxcolret)
INTEGER :: iret(mxcolret)
INTEGER :: isrcret(0:mxret)
INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
! Input Analysis Control Variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: srcsng(nsrcsng)
CHARACTER (LEN=8) :: srcua (nsrcua )
CHARACTER (LEN=8) :: srcrad(nsrcrad)
CHARACTER (LEN=8) :: srcret(nsrcret)
INTEGER :: ianxtyp(mxpass)
INTEGER :: iusesng(0:nsrcsng)
INTEGER :: iuseua(0:nsrcua)
INTEGER :: iuserad(0:nsrcrad)
INTEGER :: iuseret(0:nsrcret)
REAL :: xyrange(mxpass)
REAL :: kpvrsq(nvar)
REAL :: wlim
REAL :: zrange(mxpass)
REAL :: zwlim
REAL :: thrng(mxpass)
INTEGER :: iwstat
!
!-----------------------------------------------------------------------
!
! Scratch Space
!
!-----------------------------------------------------------------------
!
REAL :: rngsqi(nvar)
INTEGER :: knt(nvar,nz)
REAL :: wgtsum(nvar,nz)
REAL :: zsum(nvar,nz)
!
!-----------------------------------------------------------------------
!
! Output Variables at Observation Locations
!
!-----------------------------------------------------------------------
!
REAL :: corsng(mxsng,nvar)
REAL :: corua(nzua,mxua,nvar)
REAL :: corrad(nzrdr,mxcolrad,nvarrad)
REAL :: corret(nzret,mxcolret,nvar)
REAL :: oanxsng(nvar,mxsng)
REAL :: oanxua(nvar,nzua,mxua)
REAL :: oanxrad(nvarrad,nzrdr,mxcolrad)
REAL :: oanxret(nvar,nzret,mxcolret)
!
!-----------------------------------------------------------------------
!
! Output Grid
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
!
!-----------------------------------------------------------------------
!
! Work arrays
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny,nz)
REAL :: tem2(nx,ny,nz)
REAL :: tem3(nx,ny,nz)
INTEGER :: ibegin(ny)
INTEGER :: iend(ny)
!
!-----------------------------------------------------------------------
!
! Return status
!
!-----------------------------------------------------------------------
!
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Misc.local variables
!
!-----------------------------------------------------------------------
!
REAL :: ftabinv,setexp
INTEGER :: i,j,k,isrc
REAL :: rpass,zrngsq,thrngsq
!
!
!
! input:
! -----
! numctr: dimension of control variables
! ctrv: current perturbation (1d arrays)
!
! output:
! ------
! gradient of costfunction
!
!------------------------------------------------------------------------------c
!
!
INTEGER :: numctr,iflag,num
REAL :: ctrv (numctr), grad (numctr)
REAL :: gdu_err(nx,ny,nz)
REAL :: gdv_err(nx,ny,nz)
REAL :: gdp_err(nx,ny,nz)
REAL :: gdt_err(nx,ny,nz)
REAL :: gdq_err(nx,ny,nz)
REAL :: gdw_err(nx,ny,nz)
REAL :: gdscal(nx,ny,nz)
REAL :: ugrid,vgrid,wgrid,vr,zv
REAL :: u_ctr(nx,ny,nz)
REAL :: v_ctr(nx,ny,nz)
REAL :: p_ctr(nx,ny,nz)
REAL :: t_ctr(nx,ny,nz)
REAL :: q_ctr(nx,ny,nz)
REAL :: w_ctr(nx,ny,nz)
REAL :: psi(nx,ny,nz)
REAL :: phi(nx,ny,nz)
!
REAL,DIMENSION (:,:,:), allocatable :: &
u,v,w,wcont
!
REAL,DIMENSION (:,:,:), allocatable :: &
u_grd,v_grd,p_grd,t_grd,q_grd,w_grd
!
REAL :: sum1,sum2
integer:: isum
REAL :: f_div
REAL :: div3(nx,ny,nz)
!
INTEGER :: ivar
REAL :: temp1,temp2
!
!
allocate ( u(nx,ny,nz) )
allocate ( v(nx,ny,nz) )
allocate ( w(nx,ny,nz) )
allocate ( wcont(nx,ny,nz) )
!
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
u_ctr(i,j,k) = 0.
v_ctr(i,j,k) = 0.
psi(i,j,k) = 0.
phi(i,j,k) = 0.
p_ctr(i,j,k) = 0.
t_ctr(i,j,k) = 0.
q_ctr(i,j,k) = 0.
w_ctr(i,j,k) = 0.
END DO
END DO
END DO
!
!
CALL transi
(numctr,nx,ny,nz,u_ctr,v_ctr,p_ctr, &
t_ctr,q_ctr,w_ctr,ctrv)
IF( 2 == 1 ) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
u_ctr(i,j,k) = 0.
v_ctr(i,j,k) = 0.
psi(i,j,k) = 0.
phi(i,j,k) = 0.
p_ctr(i,j,k) = 0.
t_ctr(i,j,k) = 0.
q_ctr(i,j,k) = 0.
w_ctr(i,j,k) = 0.
END DO
END DO
END DO
ENDIF
!
!
!
!
allocate (u_grd(nx,ny,nz))
allocate (v_grd(nx,ny,nz))
allocate (p_grd(nx,ny,nz))
allocate (t_grd(nx,ny,nz))
allocate (q_grd(nx,ny,nz))
allocate (w_grd(nx,ny,nz))
!
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
u_grd(i,j,k) = 0.
v_grd(i,j,k) = 0.
p_grd(i,j,k) = 0.
t_grd(i,j,k) = 0.
q_grd(i,j,k) = 0.
w_grd(i,j,k) = 0.
END DO
END DO
END DO
!
!
!
! observation for single level data
! --------------------------------------------------------------
!
IF(sngsw > 0) THEN
!
num = 0
DO i = 1, nobsng
!
!
CALL alinearint_3d(nx,ny,nz,u_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
1, mxsng, icatsng, nobsng, &
1,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,1) )
!
CALL alinearint_3d(nx,ny,nz,v_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
1, mxsng, icatsng, nobsng, &
2,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,2) )
!
CALL alinearint_3d(nx,ny,nz,p_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
1, mxsng, icatsng, nobsng, &
3,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,3) )
!
CALL alinearint_3d(nx,ny,nz,t_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
1, mxsng, icatsng, nobsng, &
4,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,4) )
!
CALL alinearint_3d(nx,ny,nz,q_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
1, mxsng, icatsng, nobsng, &
5,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,5) )
!
!
!
! call alinearint_3d(nx,ny,nz,w_grd(1,1,1),xs(1),ys(1),zs(1,1,1),
! : 6,xsng(i),ysng(i),hgtsng(i),obsng(6,i),iflag )
!
!
IF(iflag == 1) num = num+1
END DO
!
! print*, 'num of valid data in grad cal===',num
!
END IF
!
! sum1=0.0
! do i =1, nx
! do j =1, ny
! do k =1, nz
! sum1 = sum1 + u_grd(i,j,k)*u_grd(i,j,k)
! sum1 = sum1 + v_grd(i,j,k)*v_grd(i,j,k)
! sum1 = sum1 + p_grd(i,j,k)*p_grd(i,j,k)
! sum1 = sum1 + t_grd(i,j,k)*t_grd(i,j,k)
! sum1 = sum1 + q_grd(i,j,k)*q_grd(i,j,k)
! end do
! end do
! end do
! print*,'sum1===',sum1
! stop
!
!
!
! observation cost function for upper level data
! --------------------------------------------------------------
!
IF(uasw > 0) THEN
!
!
do i = 1, nobsua
DO j = 1, nlevsua(i)
!
!
CALL alinearint_3d(nx,ny,nz,u_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzua, mxua, nlevsua, nobsua, &
1,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,1) )
CALL alinearint_3d(nx,ny,nz,v_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzua, mxua, nlevsua, nobsua, &
2,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,2) )
CALL alinearint_3d(nx,ny,nz,p_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzua, mxua, nlevsua, nobsua, &
3,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,3) )
CALL alinearint_3d(nx,ny,nz,t_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzua, mxua, nlevsua, nobsua, &
4,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,4) )
CALL alinearint_3d(nx,ny,nz,q_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzua, mxua, nlevsua, nobsua, &
5,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,5) )
!
! CALL alinearint_3d(nx,ny,nz,w_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
! nzua, mxua, nlevsua, nobsua, &
! 6,xua_p,yua_p,zua_1,zua_2,hgtua,ihgtua,corua(1,1,6) )
!
!
END DO
END DO
!
END IF
!
!
! sum1=0.0
! do i =1, nx
! do j =1, ny
! do k =1, nz
! sum1 = sum1 + u_grd(i,j,k)*u_grd(i,j,k)
! sum1 = sum1 + v_grd(i,j,k)*v_grd(i,j,k)
! sum1 = sum1 + p_grd(i,j,k)*p_grd(i,j,k)
! sum1 = sum1 + t_grd(i,j,k)*t_grd(i,j,k)
! sum1 = sum1 + q_grd(i,j,k)*q_grd(i,j,k)
! end do
! end do
! end do
! print*,'sum1===',sum1
! stop
!
! loading radar data
! --------------------------------------------------------
!
sum1 = 0.0
sum2 = 0.0
!
IF(radsw > 0) THEN
!
corrad=0.0
sum1= 0.
sum2= 0.
isum = 0
DO i = 1, ncolrad
! if(iuserad(isrcrad(irad(i))).GT.0) THEN
!
DO k=1,nlevrad(i)
!
IF(qualrad(2,k,i) > 0 .and. ihgtradc(k,i)>=0 ) THEN
! corrad(k,i,1) = uazmrad(i)*corrad(k,i,2)*dsdr(k,i)
! corrad(k,i,3) = corrad(k,i,2)*dhdr(k,i)
! corrad(k,i,2) = vazmrad(i)*corrad(k,i,2)*dsdr(k,i)
!
! sum2 = sum2 + obsrad(2,k,i)*obsrad(2,k,i)
corrad(k,i,1) = uazmrad(i)*obsrad(2,k,i)*dsdr(k,i)
corrad(k,i,3) = obsrad(2,k,i)*dhdr(k,i)
corrad(k,i,2) = vazmrad(i)*obsrad(2,k,i)*dsdr(k,i)
isum = isum +1
sum1=sum1+corrad(k,i,1)*corrad(k,i,1)
sum1=sum1+corrad(k,i,2)*corrad(k,i,2)
sum1=sum1+corrad(k,i,3)*corrad(k,i,3)
! ugrid = uazmrad(i)*obsrad(2,k,i)*dsdr(k,i)
! wgrid = obsrad(2,k,i)*dhdr(k,i)
! vgrid = vazmrad(i)*obsrad(2,k,i)*dsdr(k,i)
! sum1=sum1+ugrid*ugrid
! sum1=sum1+wgrid*wgrid
! sum1=sum1+vgrid*vgrid
! sum1=sum1+dsdr(k,i)*dsdr(k,i)+dhdr(k,i)*dhdr(k,i)
! sum2=sum2+uazmrad(i)*uazmrad(i)+vazmrad(i)*vazmrad(i)
ENDIF
!
END DO
! END IF
!
END DO
print*,'sum1=========================',sum1,sum2
!
CALL alinearint_3d(nx,ny,nz,u_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzrdr, mxcolrad, nlevrad, ncolrad, &
1,xradc_p,yradc_p,zradc_1,zradc_2,hgtradc, &
ihgtradc,corrad(1,1,1) )
!
CALL alinearint_3d(nx,ny,nz,v_grd(1,1,1),xs(1),ys(1),zs(1,1,1), &
nzrdr, mxcolrad, nlevrad, ncolrad, &
2,xradc_p,yradc_p,zradc_1,zradc_2,hgtradc, &
ihgtradc,corrad(1,1,2) )
!
CALL alinearint_3d(nx,ny,nz,w_grd(1,1,1),xs(1),ys(1),zp(1,1,1), &
nzrdr, mxcolrad, nlevrad, ncolrad, &
6,xradc_p,yradc_p,zradc_1,zradc_2,hgtradc, &
ihgtradc,corrad(1,1,3) )
!
!
END IF
!
! sum1=0.0
! do i =1, nx
! do j =1, ny
! do k =1, nz
! sum1 = sum1 + u_grd(i,j,k)*u_grd(i,j,k)
! sum1 = sum1 + v_grd(i,j,k)*v_grd(i,j,k)
! sum1 = sum1 + w_grd(i,j,k)*w_grd(i,j,k)
! end do
! end do
! end do
! print*,'sum1===',sum1,isum
! stop
IF( turn_div == 1 ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k) = 0.0
tem2(i,j,k) = 0.0
tem3(i,j,k) = 0.0
u(i,j,k) = 0.0
v(i,j,k) = 0.0
w(i,j,k) = 0.0
wcont(i,j,k) = 0.0
END DO
END DO
END DO
IF( (wei_div_h > 0.0) .and. (wei_div_v > 0.0)) THEN
DO k=3,nz-2
DO j=2,ny-2
DO i=2,nx-2
tem1(i+1,j,k)=tem1(i+1,j,k) + 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dxinv
tem1(i ,j,k)=tem1(i ,j,k) - 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dxinv
tem2(i,j+1,k)=tem2(i,j+1,k) + 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dyinv
tem2(i ,j,k)=tem2(i ,j,k) - 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dyinv
tem3(i,j,k+1)=tem3(i,j,k+1) + 1./wei_div_v*j3inv(i,j,k)* &
div3(i,j,k)*dzinv
tem3(i ,j,k)=tem3(i ,j,k) - 1./wei_div_v*j3inv(i,j,k)* &
div3(i,j,k)*dzinv
END DO
END DO
END DO
ELSEIF( (wei_div_h > 0.0) .and. (wei_div_v < 0.0)) THEN
DO k=3,nz-2
DO j=2,ny-2
DO i=2,nx-2
tem1(i+1,j,k)=tem1(i+1,j,k) + 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dxinv
tem1(i ,j,k)=tem1(i ,j,k) - 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dxinv
tem2(i,j+1,k)=tem2(i,j+1,k) + 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dyinv
tem2(i ,j,k)=tem2(i ,j,k) - 1./wei_div_h*j3inv(i,j,k)* &
mapfct(i,j,7)*div3(i,j,k)*dyinv
END DO
END DO
END DO
ENDIF
IF(wei_div_v > 0.0 ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
wcont(i,j,k)=tem3(i,j,k)*rhostrw(i,j,k)
tem3(i,j,k)=0.0
END DO
END DO
END DO
ENDIF
IF(wei_div_h > 0.0 ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
v(i,j,k)=tem2(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6)
tem2(i,j,k)=0.0
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
u(i,j,k)=tem1(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5)
tem1(i,j,k) = 0.0
END DO
END DO
END DO
ENDIF
!
CALL adwcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
!
DO k= 1,nz
DO j= 1,ny
DO i= 1,nx
u_grd(i,j,k) = u_grd(i,j,k) + u(i,j,k)
v_grd(i,j,k) = v_grd(i,j,k) + v(i,j,k)
w_grd(i,j,k) = w_grd(i,j,k) + w(i,j,k)
u(i,j,k) = 0.0
v(i,j,k) = 0.0
w(i,j,k) = 0.0
END DO
END DO
END DO
END IF
!
!
!
!--IF cntl_var == 0 u & v are control varibles---------------------
!
!
IF( cntl_var == 0 ) THEN
DO k = 1, nz
DO i = 1, nx
DO j = 1, ny
!
psi(i,j,k)= u_grd(i,j,k)
phi(i,j,k)= v_grd(i,j,k)
!
END DO
END DO
END DO
ELSE
DO k = 1, nz
v_grd(nx-1,ny,k)=v_grd(nx-1,ny,k)+0.5*v_grd(nx,ny,k)
v_grd(nx,ny-1,k)=v_grd(nx,ny-1,k)+0.5*v_grd(nx,ny,k)
v_grd(nx,ny,k)=0.0
u_grd(nx-1,ny,k)=u_grd(nx-1,ny,k)+0.5*u_grd(nx,ny,k)
u_grd(nx,ny-1,k)=u_grd(nx,ny-1,k)+0.5*u_grd(nx,ny,k)
u_grd(nx,ny,k)=0.0
v_grd(nx,2,k)=v_grd(nx,2,k)+0.5*v_grd(nx,1,k)
v_grd(nx-1,1,k)=v_grd(nx-1,1,k)+0.5*v_grd(nx,1,k)
v_grd(nx,1,k)=0.0
u_grd(nx,2,k)=u_grd(nx,2,k)+0.5*u_grd(nx,1,k)
u_grd(nx-1,1,k)=u_grd(nx-1,1,k)+0.5*u_grd(nx,1,k)
u_grd(nx,1,k)=0.0
v_grd(2,ny,k)=v_grd(2,ny,k)+0.5*v_grd(1,ny,k)
v_grd(1,ny-1,k)=v_grd(1,ny-1,k)+0.5*v_grd(1,ny,k)
v_grd(1,ny,k)=0.0
u_grd(2,ny,k)=u_grd(2,ny,k)+0.5*u_grd(1,ny,k)
u_grd(1,ny-1,k)=u_grd(1,ny-1,k)+0.5*u_grd(1,ny,k)
u_grd(1,ny,k)=0.0
v_grd(1,2,k)=v_grd(1,2,k)+0.5*v_grd(1,1,k)
v_grd(2,1,k)=v_grd(2,1,k)+0.5*v_grd(1,1,k)
v_grd(1,1,k)=0.0
u_grd(1,2,k)=u_grd(1,2,k)+0.5*u_grd(1,1,k)
u_grd(2,1,k)=u_grd(2,1,k)+0.5*u_grd(1,1,k)
u_grd(1,1,k)=0.0
DO i=2,nx-1
v_grd(i,ny-2,k)=v_grd(i,ny-2,k)-v_grd(i,ny,k)
v_grd(i,ny-1,k)=v_grd(i,ny-1,k)+v_grd(i,ny,k)
v_grd(i,ny-1,k)=v_grd(i,ny-1,k)+v_grd(i,ny,k)
v_grd(i,ny,k)=0.0
u_grd(i,ny-2,k)=u_grd(i,ny-2,k)-u_grd(i,ny,k)
u_grd(i,ny-1,k)=u_grd(i,ny-1,k)+u_grd(i,ny,k)
u_grd(i,ny-1,k)=u_grd(i,ny-1,k)+u_grd(i,ny,k)
u_grd(i,ny,k)=0.0
v_grd(i, 3,k)=v_grd(i, 3,k)-v_grd(i, 1,k)
v_grd(i, 2,k)=v_grd(i, 2,k)+v_grd(i, 1,k)
v_grd(i, 2,k)=v_grd(i, 2,k)+v_grd(i, 1,k)
v_grd(i, 1,k)=0.0
u_grd(i, 3,k)=u_grd(i, 3,k)-u_grd(i, 1,k)
u_grd(i, 2,k)=u_grd(i, 2,k)+u_grd(i, 1,k)
u_grd(i, 2,k)=u_grd(i, 2,k)+u_grd(i, 1,k)
u_grd(i, 1,k)=0.0
END DO
DO j=2,ny-1
v_grd(nx-2,j,k)=v_grd(nx-2,j,k)-v_grd(nx,j,k)
v_grd(nx-1,j,k)=v_grd(nx-1,j,k)+v_grd(nx,j,k)
v_grd(nx-1,j,k)=v_grd(nx-1,j,k)+v_grd(nx,j,k)
v_grd(nx,j,k)=0.0
u_grd(nx-2,j,k)=u_grd(nx-2,j,k)-u_grd(nx,j,k)
u_grd(nx-1,j,k)=u_grd(nx-1,j,k)+u_grd(nx,j,k)
u_grd(nx-1,j,k)=u_grd(nx-1,j,k)+u_grd(nx,j,k)
u_grd(nx,j,k)=0.0
v_grd( 3,j,k)=v_grd( 3,j,k)-v_grd( 1,j,k)
v_grd( 2,j,k)=v_grd( 2,j,k)+v_grd( 1,j,k)
v_grd( 2,j,k)=v_grd( 2,j,k)+v_grd( 1,j,k)
v_grd( 1,j,k)=0.0
u_grd( 3,j,k)=u_grd( 3,j,k)-u_grd( 1,j,k)
u_grd( 2,j,k)=u_grd( 2,j,k)+u_grd( 1,j,k)
u_grd( 2,j,k)=u_grd( 2,j,k)+u_grd( 1,j,k)
u_grd( 1,j,k)=0.0
END DO
DO i = 2, nx-1
DO j = 2, ny-1
!
u_grd(i,j,k) = u_grd(i,j,k)*mapfct(i,j,2)
psi(i-1,j+1,k)= psi(i-1,j+1,k)+ u_grd(i,j,k)/dy/4.
psi(i,j+1,k) = psi(i,j+1,k) + u_grd(i,j,k)/dy/4.
psi(i-1,j-1,k)= psi(i-1,j-1,k)- u_grd(i,j,k)/dy/4.
psi(i,j-1,k) = psi(i,j-1,k) - u_grd(i,j,k)/dy/4.
phi(i, j, k)= phi(i, j, k)+ u_grd(i,j,k)/dx
phi(i-1,j,k) = phi(i-1,j,k) - u_grd(i,j,k)/dx
u_grd(i,j,k) = 0.0
!
v_grd(i,j,k) = v_grd(i,j,k)*mapfct(i,j,3)
psi(i+1,j-1,k)= psi(i+1,j-1,k)+ v_grd(i,j,k)/dx/4.
psi(i+1,j,k) = psi(i+1,j,k) + v_grd(i,j,k)/dx/4.
psi(i-1,j-1,k)= psi(i-1,j-1,k)- v_grd(i,j,k)/dx/4.
psi(i-1,j,k) = psi(i-1,j,k) - v_grd(i,j,k)/dx/4.
phi(i, j, k)= phi(i, j, k)+ v_grd(i,j,k)/dy
phi(i,j-1,k) = phi(i,j-1,k) - v_grd(i,j,k)/dy
v_grd(i,j,k) = 0.0
END DO
END DO
END DO
!
END IF
!
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdu_err, &
gdscal, psi)
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdv_err, &
gdscal, phi)
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdp_err, &
gdscal,p_grd)
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdt_err, &
gdscal,t_grd)
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdq_err, &
gdscal,q_grd)
!
CALL vbl_to_ctr(ipass_filt,hradius,nradius_z,nx,ny,nz,gdw_err, &
gdscal,w_grd)
!
!
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
u_grd(i,j,k) = 0.0
v_grd(i,j,k) = 0.0
u_grd(i,j,k) = psi(i,j,k) + u_ctr(i,j,k)
v_grd(i,j,k) = phi(i,j,k) + v_ctr(i,j,k)
p_grd(i,j,k) = p_grd(i,j,k) + p_ctr(i,j,k)
t_grd(i,j,k) = t_grd(i,j,k) + t_ctr(i,j,k)
q_grd(i,j,k) = q_grd(i,j,k) + q_ctr(i,j,k)
w_grd(i,j,k) = w_grd(i,j,k) + w_ctr(i,j,k)
END DO
END DO
END DO
!
!
CALL trans(numctr,nx,ny,nz,u_grd,v_grd,p_grd, &
t_grd,q_grd,w_grd,grad)
!
!
deallocate( u )
deallocate( v )
deallocate( w )
deallocate( wcont )
!
deallocate ( u_grd )
deallocate ( v_grd )
deallocate ( p_grd )
deallocate ( t_grd )
deallocate ( q_grd )
deallocate ( w_grd )
!
!
RETURN
END SUBROUTINE gradt