!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MINIMIZATION ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Main driver for gradient check and minimization.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! Jidong Gao, CAPS, July, 2000
!
!-----------------------------------------------------------------------
!
SUBROUTINE minimization(nx,ny,nz, & 1,33
nvar,nvarradin,nvarrad,nzua,nzrdr,nzret, &
mapfct,j1,j2,j3,aj3x,aj3y,aj3z,j3inv,rhostr, &
mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret, &
nsrcsng,nsrcua,nsrcrad,nsrcret,ncat, &
mxpass,ipass,iwstat, &
xs,ys,zs,x,y,z,zp,hterain, &
icatg,xcor,nam_var, &
xsng,ysng,hgtsng,thesng, &
odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng, &
xua,yua,hgtua,theua, &
odifua,qobsua,qualua,isrcua,nlevsua,nobsua, &
elvrad,xradc,yradc, &
distrad,uazmrad,vazmrad,hgtradc,theradc,dsdr,dhdr, &
odifrad,qobsrad,qualrad, &
irad,isrcrad,nlevrad,ncolrad, &
xretc,yretc,hgtretc,theretc, &
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, &
oanxsng,oanxua,oanxrad,oanxret, &
nz_tab,hqback,qback, &
nsngfil,nuafil,nradfil,nretfil, &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nvar Number of analysis variables
! nvarradin Number of variables in the obsrad array
! nvarrad Number of variables in the odifrad array
! nzua Maximumnumber of levels in UA observations
! nzrdr Maximum number of levels in a radar column
! nzret Maximum number of levels in a retrieval column
! mxsng Maximum number of single level data
! mxua Maximum number of upper air observations
! mxrad Maximum number of radars
! mxcolrad Maximum number of radar columns
! mxcolret Maximum number of retrieval columns
! mxpass Maximum number of iterations
! npass Number of iterations
! iwstat Status indicator for writing statistics
!
! xs x location of scalar pts (m)
! ys y location of scalar pts (m)
! zs z location of scalar pts (m)
!
! xsng x location of single-level data
! ysng y location of single-level data
! hgtsng elevation of single-level data
! thesng theta (potential temperature) of single-level data
!
! obsng single-level observations
! odifsng difference between single-level obs and analysis
! qobsng normalized observation error
! qualsng single-level data quality indicator
! isrcsng index of single-level data source
! nobsng number of single-level observations
!
! xua x location of upper air data
! yua y location of upper air data
! hgtua elevation of upper air data
! theua theta (potential temperature) of upper air data
!
! obsua upper air observations
! odifua difference between upper air obs and analysis
! qobsua normalized observation error
! qualua upper air data quality indicator
! isrcua index of upper air data source
! nlevsua number of levels of data for each upper air location
! nobsua number of upper air observations
!
! anx Analyzed variables (or first guess)
! qback Background (first guess) error
!
! nradfil number of radar files
! fradname file name for radar dataset
! srcrad name of radar sources
!
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of feed horn of radar (m MSL)
! xradc x location of radar column
! yradc y location of radar column
! irad radar number
! nlevrad number of levels of radar data in each column
! distrad distance of radar column from source radar
! uazmrad u-component of radar beam for each column
! vazmrad v-component of radar beam for each column
! hgtradc height (m MSL) of radar observations
! obsrad radar observations
! oanxrad analysis (first guess) value at radar data location
! odifrad difference between radar observation and analysis
! qobsrad normalized observation error
! qualrad radar data quality indicator
! ncolrad number of radar columns read-in
! istatus status indicator
!
! latret latitude of retrieval radar (degrees N)
! lonret longitude of retrieval radar (degrees E)
! xretc x location of retrieval column
! yretc y location of retrieval column
! iret retrieval number
! nlevret number of levels of retrieval data in each column
! hgtretc height (m MSL) of retrieval observations
! obsret retrieval observations
! odifret difference between retrieval observation and analysis
! qobsret normalized observation error
! qualret retrieval data quality indicator
! ncolret number of retr columns read-in
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'grid.inc'
INCLUDE 'varpara.inc'
INCLUDE '3dvarcputime.inc'
!-----------------------------------------------------------------------
!
! Input Sizing Arguments
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
INTEGER :: nvar,nvarradin,nvarrad
INTEGER :: nzua,nzrdr,nzret
INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
INTEGER :: mxpass,ipass
INTEGER :: nsngfil,nuafil,nradfil,nretfil
!
!-----------------------------------------------------------------------
!
! 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 :: 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, DIMENSION(:), allocatable:: xsng_p, ysng_p
REAL, DIMENSION(:), allocatable:: zsng_1, zsng_2
INTEGER, DIMENSION(:), allocatable:: ihgtsng
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 :: xua(mxua)
REAL :: yua(mxua)
REAL, DIMENSION(:), allocatable:: xua_p, yua_p
REAL, DIMENSION(:,:), allocatable:: zua_1, zua_2
INTEGER, DIMENSION(:,:), allocatable:: ihgtua
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 :: elvrad(mxrad)
REAL :: xradc(mxcolrad)
REAL :: yradc(mxcolrad)
REAL, DIMENSION(:), allocatable:: xradc_p, yradc_p
REAL, DIMENSION(:,:), allocatable:: zradc_1, zradc_2
INTEGER, DIMENSION(:,:), allocatable:: ihgtradc
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 :: 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,mxpass)
INTEGER :: iuseua(0:nsrcua,mxpass)
INTEGER :: iuserad(0:nsrcrad,mxpass)
INTEGER :: iuseret(0:nsrcret,mxpass)
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)
!
INTEGER :: nz_tab
REAL :: hqback(nz_tab)
REAL :: qback(nvar,nz_tab)
!
REAL,DIMENSION(nz_tab) :: qback_wk
REAL,DIMENSION(nz_tab) :: hqback_wk
!
!-----------------------------------------------------------------------
!
! Output Grid
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
!
!
!
!-----------------------------------------------------------------------
!
! variable declarations for varirational analysis: by gao
!
!-----------------------------------------------------------------------
!
!
INTEGER :: numctr
!
INTEGER :: ndim,mgra,nxm,nwork
!
REAL :: cfun,eps,gtol,ftol
REAL :: cfun_single(nvar+1),cfun_total
INTEGER :: mp,lp,iprint(2),iflag,point
LOGICAL :: diagco
COMMON /va15dd/mp,lp, gtol
!
INTEGER :: icall
!
REAL, DIMENSION(:), allocatable:: ctrv,grad,xgus
REAL, DIMENSION(:), allocatable:: swork,ywork,diag,work
!
!
REAL :: cfun1,rchek,gxnn
REAL,DIMENSION(:), allocatable:: fa, fr
REAL,DIMENSION(:,:,:), allocatable:: gdscal
REAL,DIMENSION(:,:,:), allocatable:: gdu_err,gdv_err, &
gdp_err,gdt_err,gdq_err,gdw_err
REAL,DIMENSION (:,:,:), allocatable :: &
u_ctr,v_ctr,p_ctr,t_ctr,q_ctr,w_ctr,psi,phi
REAL,DIMENSION (:,:,:), allocatable :: &
rhostru, rhostrv,rhostrw,div3
!
!
!-----------------------------------------------------------------------
!
! end of variable declarations for variational anaylsis:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! 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,ierr
!
!-----------------------------------------------------------------------
!
! Misc.local variables
!
!-----------------------------------------------------------------------
!
REAL :: ftabinv,setexp
INTEGER :: i,j,k,isrc,sngsw,uasw,radsw,retsw
INTEGER :: num, numsingle
REAL :: rpass,zrngsq,thrngsq
REAL :: sum1,sum2,sum3,sum4,sum5,sum6
REAL :: f_cputime
REAL :: cpuin,cpuout
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! sum1=0
! do i = 1, nobsng
! sum1=sum1+odifsng(3,i)**2
! end do
! print*,'sum1_odifsng==',sum1
!
!
! allocate control variable array
! -------------------------------------------------------
!
numctr=6*nx*ny*nz
numsingle=nx*ny*nz
if( allocated(xgus) ) deallocate (xgus)
allocate ( xgus(numctr), stat=ierr )
if( allocated(ctrv) ) deallocate (ctrv)
allocate ( ctrv(numctr) )
if( allocated(grad) ) deallocate (grad)
allocate ( grad(numctr) )
if( allocated(gdu_err) ) deallocate (gdu_err)
allocate ( gdu_err(nx,ny,nz) )
if( allocated(gdv_err) ) deallocate (gdv_err)
allocate ( gdv_err(nx,ny,nz) )
if( allocated(gdp_err) ) deallocate (gdp_err)
allocate ( gdp_err(nx,ny,nz) )
if( allocated(gdt_err) ) deallocate (gdt_err)
allocate ( gdt_err(nx,ny,nz) )
if( allocated(gdq_err) ) deallocate (gdq_err)
allocate ( gdq_err(nx,ny,nz) )
if( allocated(gdw_err) ) deallocate (gdw_err)
allocate ( gdw_err(nx,ny,nz) )
!------------------------------------------------------------------------
if( allocated( u_ctr) ) deallocate ( u_ctr)
allocate ( u_ctr(nx,ny,nz) )
if( allocated( v_ctr) ) deallocate ( v_ctr)
allocate ( v_ctr(nx,ny,nz) )
if( allocated( psi) ) deallocate ( psi)
allocate ( psi(nx,ny,nz) )
if( allocated( phi) ) deallocate ( phi)
allocate ( phi(nx,ny,nz) )
if( allocated( p_ctr) ) deallocate ( p_ctr)
allocate ( p_ctr(nx,ny,nz) )
if( allocated( t_ctr) ) deallocate ( t_ctr)
allocate ( t_ctr(nx,ny,nz) )
if( allocated( q_ctr) ) deallocate ( q_ctr)
allocate ( q_ctr(nx,ny,nz) )
if( allocated( w_ctr) ) deallocate ( w_ctr)
allocate ( w_ctr(nx,ny,nz) )
!-----------------------------------------------------------------------
if( allocated(gdscal) ) deallocate (gdscal)
allocate ( gdscal(nx,ny,nz) )
if( allocated(rhostru) ) deallocate (rhostru)
allocate ( rhostru(nx,ny,nz) )
if( allocated(rhostrv) ) deallocate (rhostrv)
allocate ( rhostrv(nx,ny,nz) )
if( allocated(rhostrw) ) deallocate (rhostrw)
allocate ( rhostrw(nx,ny,nz) )
if( allocated(div3 ) ) deallocate (div3 )
allocate ( div3 (nx,ny,nz) )
!
CALL rhouvw
(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw)
!
! PRINT*,'qbackin=',(qback(1,i),i=1,nz_tab)
!
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
gdw_err(i,j,k) = 3.0
end do
end do
end do
num = 0
DO isrc = 1, nz_tab
if( qback(1,isrc) > -900.0) then
num = num+1
qback_wk(num) = qback(1,isrc)
hqback_wk(num) = hqback(isrc)
end if
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
CALL linear(num,hqback_wk,qback_wk,zs(i,j,k),gdu_err(i,j,k) )
END DO
END DO
END DO
!
!
num = 0
DO isrc = 1, nz_tab
if( qback(2,isrc) > -900.0) then
num = num+1
qback_wk(num) = qback(2,isrc)
hqback_wk(num) = hqback(isrc)
end if
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
CALL linear(num,hqback_wk,qback_wk,zs(i,j,k),gdv_err(i,j,k) )
END DO
END DO
END DO
!
!
num = 0
DO isrc = 1, nz_tab
if( qback(3,isrc) > -900.0) then
num = num+1
qback_wk(num) = qback(3,isrc)
hqback_wk(num) = hqback(isrc)
end if
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
CALL linear(num,hqback_wk,qback_wk,zs(i,j,k),gdp_err(i,j,k) )
END DO
END DO
END DO
!
!
num = 0
DO isrc = 1, nz_tab
if( qback(4,isrc) > -900.0) then
num = num+1
qback_wk(num) = qback(4,isrc)
hqback_wk(num) = hqback(isrc)
end if
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
CALL linear(num,hqback_wk,qback_wk,zs(i,j,k),gdt_err(i,j,k) )
END DO
END DO
END DO
!
num = 0
DO isrc = 1, nz_tab
if( qback(5,isrc) > -900.0) then
num = num+1
qback_wk(num) = qback(5,isrc)
hqback_wk(num) = hqback(isrc)
end if
END DO
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
CALL linear(num,hqback_wk,qback_wk,zs(i,j,k),gdq_err(i,j,k) )
END DO
END DO
END DO
!
!
! num = 0
! DO isrc = 1, nz_tab
! if( qback(1,isrc) > -900.0) then
! num = num+1
! qback_wk(num) = qback(1,isrc)
! hqback_wk(num) = hqback(isrc)
! end if
! END DO
!
! DO k = 1,nz
! DO j = 1,ny
! DO i = 1,nx
! call linear(num,hqback_wk,qback_wk,zs(i,j,k),gdw_err(i,j,k))
! END DO
! END DO
! END DO
!
!
!----------------if ( cntl_var==0) u & v as control variable-------------
!
IF(cntl_var == 0) THEN
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
gdu_err(i,j,k) = SQRT( gdu_err(i,j,k) )
gdv_err(i,j,k) = SQRT( gdv_err(i,j,k) )
gdp_err(i,j,k) = SQRT( gdp_err(i,j,k) )
gdt_err(i,j,k) = SQRT( gdt_err(i,j,k) )
gdq_err(i,j,k) = SQRT( gdq_err(i,j,k) )
gdw_err(i,j,k) = sqrt( gdw_err(i,j,k) )
IF( 2 == 1 ) THEN
gdu_err(i,j,k) = 1.0
gdv_err(i,j,k) = 1.0
gdp_err(i,j,k) = 1.0
gdt_err(i,j,k) = 1.0
gdq_err(i,j,k) = 1.0
gdw_err(i,j,k) = 1.0
ENDIF
END DO
END DO
END DO
ELSE
!
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
gdu_err(i,j,k) = SQRT( gdu_err(i,j,k) )*(dx+dy)/2.
gdv_err(i,j,k) = SQRT( gdv_err(i,j,k) )*(dx+dy)/2.
gdp_err(i,j,k) = SQRT( gdp_err(i,j,k) )
gdt_err(i,j,k) = SQRT( gdt_err(i,j,k) )
gdq_err(i,j,k) = SQRT( gdq_err(i,j,k) )
gdw_err(i,j,k) = sqrt( gdw_err(i,j,k) )
END DO
END DO
END DO
!
END IF
!
DO i=1,numctr
ctrv(i)=0.0
xgus(i)=0.0
grad(i)=0.0
END DO
!
IF( nsngfil == 0 ) THEN
sngsw = 0
ELSE
sngsw = 1
END IF
IF( nuafil == 0 ) THEN
uasw = 0
ELSE
uasw = 1
ENDIF
IF( nradfil == 0 ) THEN
radsw = 0
ELSE
radsw = 1
END IF
IF( nretfil == 0 ) THEN
retsw = 0
ELSE
retsw = 1
END IF
!
! caculating some parameters for interpolation ! added by Jidong Gao
!
do i = 1, nobsng
icatsng(i) = 1
end do
if( allocated( xsng_p) ) deallocate ( xsng_p)
allocate ( xsng_p(mxsng) )
if( allocated( ysng_p) ) deallocate ( ysng_p)
allocate ( ysng_p(mxsng) )
if( allocated( zsng_1) ) deallocate ( zsng_1)
allocate ( zsng_1(mxsng) )
if( allocated( zsng_2) ) deallocate ( zsng_2)
allocate ( zsng_2(mxsng) )
if( allocated( ihgtsng) ) deallocate ( ihgtsng)
allocate ( ihgtsng(mxsng) )
call map_to_mod2(nx,ny,mxsng,nobsng,xs,ys,xsng,ysng,xsng_p,ysng_p)
call map_to_modz(1, mxsng, icatsng, nobsng,nx,ny,nz,zs, &
xsng_p,ysng_p,hgtsng,ihgtsng,zsng_1,zsng_2)
!do i=1,nobsng
! print*,'xsng==',xsng_p(i),ysng_p(i)
!end do
!do i=1,nobsng
! print*,'hgtsng==',zsng_1(i),hgtsng(i),zsng_2(i),ihgtsng(i)
!end do
!stop
!
if( allocated( xua_p) ) deallocate ( xua_p)
allocate ( xua_p(mxua) )
if( allocated( yua_p) ) deallocate ( yua_p)
allocate ( yua_p(mxua) )
if( allocated( zua_1) ) deallocate ( zua_1)
allocate ( zua_1(nzua,mxua) )
if( allocated( zua_2) ) deallocate ( zua_2)
allocate ( zua_2(nzua,mxua) )
if( allocated( ihgtua) ) deallocate ( ihgtua)
allocate ( ihgtua(nzua,mxua) )
call map_to_mod2(nx,ny,mxua,nobsua,xs,ys,xua,yua,xua_p,yua_p)
call map_to_modz(nzua, mxua, nlevsua, nobsua, nx,ny,nz,zs, &
xua_p,yua_p, hgtua, ihgtua,zua_1,zua_2)
!do i=1,nobsua
! print*,'xua==',xua_p(i),yua_p(i)
!end do
!do i=1,nobsua
! do j = 1, nlevsua(i)
! print*,'hgtua=',zua_1(j,i),hgtua(j,i),zua_2(j,i),ihgtua(j,i)
!end do
!end do
!stop
!
if( allocated( xradc_p) ) deallocate ( xradc_p)
allocate ( xradc_p(mxcolrad) )
if( allocated( yradc_p) ) deallocate ( yradc_p)
allocate ( yradc_p(mxcolrad) )
if( allocated( zradc_1) ) deallocate ( zradc_1)
allocate ( zradc_1(nzrdr,mxcolrad) )
if( allocated( zradc_2) ) deallocate ( zradc_2)
allocate ( zradc_2(nzrdr,mxcolrad) )
if( allocated( ihgtradc) ) deallocate ( ihgtradc)
allocate ( ihgtradc(nzrdr,mxcolrad) )
call map_to_mod2(nx,ny,mxcolrad,ncolrad,xs,ys,xradc,yradc, &
xradc_p,yradc_p)
call map_to_modz(nzrdr, mxcolrad, nlevrad, ncolrad,nx,ny,nz, &
zs,xradc_p,yradc_p,hgtradc,ihgtradc,zradc_1,zradc_2)
!
!do i=1,ncolrad
! print*,'xradc==',xradc_p(i),yradc_p(i)
!end do
! print*,'ncolrad===',ncolrad
!do i=1,ncolrad
! do j = 1, nlevrad(i)
! print*,'hgtrad=',i,j,zradc_1(j,i),hgtradc(j,i),zradc_2(j,i),ihgtradc(j,i)
!end do
!end do
!stop
!
!
!
!
! end of caculating some parameters for interpolation ! added by Jidong
!
! print*,'maxin=',maxin
! print*,'ipass_filt=',ipass_filt
! print*,' hradius=', hradius
! print*,'nradius_z=',nradius_z
! print*,'turn_chk=',turn_chk
! print*,'turn_3dda=',turn_3dda
! print*,'sngsw=',sngsw
! print*,'uasw=',uasw
! print*,'radsw=',radsw
! print*,'retsw=',retsw
! print*,'turn_div=',turn_div
! print*,'wei_div=',wei_div
! stop
!
IF(turn_chk == 1) THEN
ipass = 1
!
CALL scale_factor(nx,ny,nz,gdscal,ipass_filt(1),hradius(1),nradius_z(1))
rchek=1E-1
allocate (fa(20))
allocate (fr(20))
!
PRINT*,'numctr====',numctr
!
!-----------------------------------------------------------------------
!
! OBTAIN THE COST FUNCTION AND THE GRADIENT.
!
!-----------------------------------------------------------------------
!
CALL costf(numctr,ctrv,cfun_single, &
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,ipass,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, &
ianxtyp,iusesng(0,ipass),iuseua(0,ipass), &
iuserad(0,ipass),iuseret(0,ipass), &
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(ipass),hradius(ipass),nradius_z(ipass), &
turn_div,cntl_var,wei_div_h(ipass),wei_div_v(ipass), &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
cfun_total=0
DO i=1,nvar+1
cfun_total=cfun_total+cfun_single(i)
ENDDO
cfun=cfun_total
PRINT*,'cfun====',cfun
!
! do i = 1, ncolrad
! if(iuserad(isrcrad(irad(i))).GT.0) THEN
!
! do k=1,nlevrad(i)
!
! print*,'obsrad in minim_drive==',k,i,obsrad(5,k,i)
! end do
! end if
! end do
! stop
!
CALL gradt(numctr,ctrv,grad, &
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,ipass,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, &
ianxtyp,iusesng(0,ipass),iuseua(0,ipass), &
iuserad(0,ipass),iuseret(0,ipass), &
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(ipass),hradius(ipass),nradius_z(ipass), &
turn_div,cntl_var,wei_div_h(ipass),wei_div_v(ipass), &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
gxnn=0.
DO i=1,numctr
gxnn=gxnn+grad(i)*grad(i)
END DO
PRINT*,'gxnn====',gxnn
! stop
!
DO j=1,20
rchek=rchek*1E-1
DO i=1,numctr
xgus(i)=ctrv(i)+rchek*grad(i)
END DO
!
!-----------------------------------------------------------------------
!
! OBTAIN THE COST FUNCTION AND THE GRADIENT.
!
!-----------------------------------------------------------------------
!
!
CALL costf(numctr,xgus,cfun_single, &
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,ipass,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(0,ipass),iuseua(0,ipass), &
iuserad(0,ipass),iuseret(0,ipass), &
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(ipass),hradius(ipass),nradius_z(ipass), &
turn_div,cntl_var,wei_div_h(ipass),wei_div_v(ipass), &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!
cfun_total=0
DO i=1,nvar+1
cfun_total=cfun_total+cfun_single(i)
ENDDO
cfun1=cfun_total
!
PRINT*,'cfun1====',cfun1,cfun1-cfun
! stop
!
fa(j)=rchek
fr(j)=(cfun1-cfun)/(gxnn*rchek)
END DO
!
OPEN(12,FILE='grad.chk',STATUS='unknown')
!
DO j=1,16
WRITE(*,*) ' Rchek=',fa(j),' fr==',fr(j)
WRITE(12,*) fa(j),fr(j)-1.0
END DO
!
CLOSE(12)
!
deallocate (fa)
deallocate (fr)
STOP
!
END IF
!
IF(turn_3dda == 1) THEN
!
!-----------------------------------------------------------------------
!
! Data assimilation begins here.
! include file in data assimilation
!
!-----------------------------------------------------------------------
!
!
mgra=2
ndim=numctr
nxm=ndim*mgra
nwork=ndim+2*mgra
allocate ( swork(nxm) )
allocate ( ywork(nxm) )
allocate ( diag(ndim) )
allocate ( work(nwork) )
!
!
DO i=1,numctr
ctrv(i)=0.0
END DO
!
!-----------------------------------------------------------------------
!
! Loop through ipass analysis iterations
!
!-----------------------------------------------------------------------
!
! DO ipass=1,npass
icall = 0
eps = 1.0E-12
ftol=1.0E-4
iprint = 0
iflag = 0
diagco = .false.
grad = 0.0
!
!-----------------------------------------------------------------------
CALL scale_factor(nx,ny,nz,gdscal,ipass_filt(ipass), &
hradius(ipass),nradius_z(ipass))
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Set single-level usage switch based in iusesng
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a,i4)') ' Source usage switches for pass ',ipass
WRITE(6,'(/a,i4)') ' Single level data'
sngsw=0
!
IF( nsngfil /= 0 ) THEN
DO isrc=1,nsrcsng
IF(iusesng(isrc,ipass) > 0) THEN
sngsw=1
WRITE(6,'(a,a)') ' Using ',srcsng(isrc)
END IF
END DO
END IF
IF(sngsw == 0) WRITE(6,'(a)') ' none'
!
!-----------------------------------------------------------------------
!
! Set multiple-level usage switch based in iuseua
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a,i4)') ' Multiple level data'
uasw=0
IF( nuafil /= 0 ) THEN
DO isrc=1,nsrcua
IF(iuseua(isrc,ipass) > 0) THEN
uasw=1
WRITE(6,'(a,a)') ' Using ',srcua(isrc)
END IF
END DO
ENDIF
IF(uasw == 0) WRITE(6,'(a)') ' none'
!
!-----------------------------------------------------------------------
!
! Set radar-data usage switch based in iuserad
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a,i4)') ' Raw radar data'
radsw=0
IF( nradfil /= 0 ) THEN
DO isrc=1,nsrcrad
IF(iuserad(isrc,ipass) > 0) THEN
radsw=1
WRITE(6,'(a,a)') ' Using ',srcrad(isrc)
END IF
END DO
END IF
IF(radsw == 0) WRITE(6,'(a)') ' none'
!
!-----------------------------------------------------------------------
!
! Set ret usage switch based in iuseret
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a,i4)') ' Retrieval radar data'
retsw=0
IF( nretfil /= 0 ) THEN
DO isrc=1,nsrcret
IF(iuseret(isrc,ipass) > 0) THEN
retsw=1
WRITE(6,'(a,a)') ' Using ',srcret(isrc)
END IF
END DO
END IF
IF(retsw == 0) WRITE(6,'(a)') ' none'
!
20 CONTINUE
!-----------------------------------------------------------------------
!
! OBTAIN THE COST FUNCTION AND THE GRADIENT.
!
!-----------------------------------------------------------------------
!
print*,'icall===========',icall,sngsw,uasw,radsw
!
cpuin=f_cputime()
CALL costf(numctr,ctrv,cfun_single, &
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,ipass,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(0,ipass),iuseua(0,ipass), &
iuserad(0,ipass),iuseret(0,ipass), &
! 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(ipass),hradius(ipass),nradius_z(ipass), &
turn_div,cntl_var,wei_div_h(ipass),wei_div_v(ipass), &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
cpuout=f_cputime()
cputimearray(1)=cputimearray(1)+(cpuout-cpuin)
!
cfun_total=0
DO i=1,nvar+1
cfun_total=cfun_total+cfun_single(i)
ENDDO
! cfun=cfun_single(5)
cfun=cfun_total
!
cpuin=f_cputime()
CALL gradt(numctr,ctrv,grad, &
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,ipass,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, &
ianxtyp,iusesng(0,ipass),iuseua(0,ipass), &
iuserad(0,ipass),iuseret(0,ipass), &
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(ipass),hradius(ipass),nradius_z(ipass), &
turn_div,cntl_var,wei_div_h(ipass),wei_div_v(ipass), &
anx,tem1,tem2,tem3,ibegin,iend,istatus)
cpuout=f_cputime()
cputimearray(2)=cputimearray(2)+(cpuout-cpuin)
!
!
!-----------------------------------------------------------------------
!
! Call the LBFGS minimization algorithm.
!
!-----------------------------------------------------------------------
!
!
! IF(2 == 1) THEN
! DO i=1,4*numsingle
! grad(i)=0.0
! ENDDO
! DO i=5*numsingle+1,6*numsingle
! grad(i)=0.0
! ENDDO
! ENDIF
cpuin=f_cputime()
CALL va15ad(numctr,mgra,ctrv,cfun,grad,diagco,diag,iprint, &
eps,swork,ywork,point,work,iflag,ftol)
cpuout=f_cputime()
cputimearray(3)=cputimearray(3)+(cpuout-cpuin)
!
!
IF(iflag <= 0) GO TO 50
icall=icall + 1
IF(icall > maxin(ipass) ) GO TO 50
GO TO 20
50 CONTINUE
!
! END DO ! end of do loop over ipass
!
deallocate ( swork )
deallocate ( ywork )
deallocate ( diag )
deallocate ( work )
!
PRINT*,'-----------------after minimization!---------------------'
!
!-----------------------------------------------------------------------
!
!c---------------------allocate the arrays for ----------------------'
!c-----------------storing the optimal perturbation--------------------'
!c
!c
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
u_ctr(i,j,k) = 0.0
v_ctr(i,j,k) = 0.0
psi(i,j,k) = 0.0
phi(i,j,k) = 0.0
p_ctr(i,j,k) = 0.0
t_ctr(i,j,k) = 0.0
q_ctr(i,j,k) = 0.0
w_ctr(i,j,k) = 0.0
END DO
END DO
END DO
!
!
CALL transi(numctr,nx,ny,nz, psi, phi, p_ctr, &
t_ctr,q_ctr,w_ctr,ctrv)
!
!
! --------------------------------------------------------
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdu_err, gdscal, psi)
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdv_err, gdscal, phi)
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdp_err, gdscal,p_ctr)
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdt_err, gdscal,t_ctr)
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdq_err, gdscal,q_ctr)
CALL ctr_to_vbl(ipass_filt(ipass),hradius(ipass),nradius_z(ipass), &
nx,ny,nz,gdw_err, gdscal,w_ctr)
!
!------if cntl_var =0 option if U,V as control variables --------------------
!------else if cntl_var =1 option, then psi, phi as control variable----------
!
!
IF(cntl_var == 0 ) THEN
!
DO k = 1, nz
DO i = 1, nx
DO j = 1, ny
u_ctr(i,j,k) = psi(i,j,k)
v_ctr(i,j,k) = phi(i,j,k)
END DO
END DO
END DO
ELSE
!
DO k = 1, nz
DO i = 2, nx-1
DO j = 2, ny-1
u_ctr(i,j,k) = ( psi(i-1,j+1,k)+psi(i,j+1,k) &
-psi(i-1,j-1,k)-psi(i,j-1,k) )/dy/4. &
+( phi(i, j, k)-phi(i-1,j,k) )/dx
u_ctr(i,j,k) = u_ctr(i,j,k)*mapfct(i,j,2)
v_ctr(i,j,k) = ( psi(i+1,j-1,k)+psi(i+1,j,k) &
-psi(i-1,j-1,k)-psi(i-1,j,k) )/dx/4. &
+( phi(i, j, k)-phi(i,j-1,k) )/dy
v_ctr(i,j,k) = v_ctr(i,j,k)*mapfct(i,j,3)
END DO
END DO
!
DO j=2,ny-1
u_ctr( 1,j,k)=u_ctr( 2,j,k)+u_ctr( 2,j,k)-u_ctr( 3,j,k)
v_ctr( 1,j,k)=v_ctr( 2,j,k)+v_ctr( 2,j,k)-v_ctr( 3,j,k)
u_ctr(nx,j,k)=u_ctr(nx-1,j,k)+u_ctr(nx-1,j,k)-u_ctr(nx-2,j,k)
v_ctr(nx,j,k)=v_ctr(nx-1,j,k)+v_ctr(nx-1,j,k)-v_ctr(nx-2,j,k)
END DO
DO i=2,nx-1
u_ctr(i, 1,k)=u_ctr(i, 2,k)+u_ctr(i, 2,k)-u_ctr(i, 3,k)
v_ctr(i, 1,k)=v_ctr(i, 2,k)+v_ctr(i, 2,k)-v_ctr(i, 3,k)
u_ctr(i,ny,k)=u_ctr(i,ny-1,k)+u_ctr(i,ny-1,k)-u_ctr(i,ny-2,k)
v_ctr(i,ny,k)=v_ctr(i,ny-1,k)+v_ctr(i,ny-1,k)-v_ctr(i,ny-2,k)
END DO
u_ctr(1,1 ,k)=0.5*( u_ctr(2,1,k)+u_ctr(1,2,k) )
v_ctr(1,1 ,k)=0.5*( v_ctr(2,1,k)+v_ctr(1,2,k) )
u_ctr(1,ny,k)=0.5*( u_ctr(1,ny-1,k)+u_ctr(2,ny,k) )
v_ctr(1,ny,k)=0.5*( v_ctr(1,ny-1,k)+v_ctr(2,ny,k) )
u_ctr(nx,1,k)=0.5*( u_ctr(nx-1,1,k)+u_ctr(nx,2,k) )
v_ctr(nx,1,k)=0.5*( v_ctr(nx-1,1,k)+v_ctr(nx,2,k) )
u_ctr(nx,ny,k)=0.5*( u_ctr(nx,ny-1,k)+u_ctr(nx-1,ny,k) )
v_ctr(nx,ny,k)=0.5*( v_ctr(nx,ny-1,k)+v_ctr(nx-1,ny,k) )
END DO
!
END IF
!
!
! call smooth ( u_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 1)
! call smooth ( v_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 1)
! call smooth ( p_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 0)
! call smooth ( t_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 0)
! call smooth ( q_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 0)
! call smooth ( w_ctr(1,1,1), NX, NY, NZ, NUM_SMTH, 0)
!
!
!
!-----------------------------------------------------------------------
!
!c add the optimal perturbation to backgd
!
!-----------------------------------------------------------------------
!
!
!----------------temparily setting by Jidong Gao------------------------
!
sum1=0.0
sum2=0.0
sum3=0.0
sum4=0.0
sum5=0.0
sum6=0.0
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
sum1=sum1+u_ctr(i,j,k)
sum2=sum2+v_ctr(i,j,k)
sum3=sum3+p_ctr(i,j,k)
sum4=sum4+t_ctr(i,j,k)
sum5=sum5+q_ctr(i,j,k)
sum6=sum6+w_ctr(i,j,k)
!
anx(i,j,k,1) = anx(i,j,k,1) + u_ctr(i,j,k)
anx(i,j,k,2) = anx(i,j,k,2) + v_ctr(i,j,k)
anx(i,j,k,3) = anx(i,j,k,3) + p_ctr(i,j,k)
anx(i,j,k,4) = anx(i,j,k,4) + t_ctr(i,j,k)
anx(i,j,k,5) = anx(i,j,k,5) + q_ctr(i,j,k)
anx(i,j,k,6) = anx(i,j,k,6) + w_ctr(i,j,k)
!
!
IF ( anx(i,j,k,5) <= 0.) anx(i,j,k,5) = 0.0
!
END DO
END DO
END DO
!
IF(ipass == 1 ) THEN
OPEN( 13, file='radialwind.dat',form='unformatted')
write(13) u_ctr
write(13) v_ctr
write(13) w_ctr
CLOSE(13)
OPEN( 13, file='wind3d.dat',form='unformatted')
write(13) (((anx(i,j,k,1),i=1,nx),j=1,ny),k=1,nz)
write(13) (((anx(i,j,k,2),i=1,nx),j=1,ny),k=1,nz)
write(13) (((anx(i,j,k,6),i=1,nx),j=1,ny),k=1,nz)
CLOSE(13)
ENDIF
!
! print*,'GDSCAL====',GDSCAL
IF(1==1) THEN
PRINT*,'u_ctr_sum=',sum1
DO j=1,8
! write(*,'(10f10.4)') (u_ctr(i,j,nz/2),i=nx/2-4,nx/2+4)
write(*,'(10f10.4)') (u_ctr(104+i,84+j,26),i=1,8)
ENDDO
PRINT*,'v_ctr_sum=',sum2
! DO j=ny/2-4,ny/2+4
DO j=1,8
write(*,'(10f10.4)') (v_ctr(104+i,84+j,26),i=1,8)
ENDDO
PRINT*,'p_ctr_sum=',sum3
DO j=ny/2-4,ny/2+4
! write(*,'(10f10.4)') (p_ctr(i,j,nz/2),i=nx/2-4,nx/2+4)
ENDDO
PRINT*,'t_ctr_sum=',sum4
DO j=ny/2-4,ny/2+4
! write(*,'(10f10.4)') (t_ctr(i,j,nz/2),i=nx/2-4,nx/2+4)
ENDDO
PRINT*,'q_ctr_sum=',sum5
DO j=ny/2-4,ny/2+4
! write(*,'(10f10.4)') (q_ctr(i,j,nz/2),i=nx/2-4,nx/2+4)
ENDDO
PRINT*,'w_ctr_sum=',sum6
ENDIF
!
! K=20
! DO I=1,NX
! write(*,'(16f5.1)') (u_ctr(I,J,K),J=1,NY)
! END DO
! DO I=1,NX
! write(*,'(16f5.1)') (v_ctr(I,J,K),J=1,NY)
! END DO
! DO K = 1, NZ
! DO I=1,NX
! write(*,'(16f5.1)') (w_ctr(I,J,K),J=1,NY)
! END DO
! END DO
! print*,'stop here? 2'
! stop
!c
!
deallocate (u_ctr)
deallocate (v_ctr)
deallocate (p_ctr)
deallocate (t_ctr)
deallocate (q_ctr)
deallocate (w_ctr)
!
END IF
!
deallocate ( xgus )
deallocate ( ctrv )
deallocate ( grad )
deallocate ( gdu_err )
deallocate ( gdv_err )
deallocate ( gdp_err )
deallocate ( gdt_err )
deallocate ( gdq_err )
deallocate ( gdw_err )
deallocate ( gdscal )
deallocate ( rhostru )
deallocate ( rhostrv )
deallocate ( rhostrw )
deallocate ( div3 )
!
!
! IF(1 == 0) THEN
!
!
!-----------------------------------------------------------------------
!
! Note in the following the analysis scratch arrays
! are used to hold the output statistics
!
!-----------------------------------------------------------------------
!
! CALL diffnst3d(nvar, &
! nzua,nzret,mxsng,mxua,mxcolret, &
! obsng,odifsng,oanxsng,isrcsng,qualsng,nobsng, &
! obsua,odifua,oanxua,isrcua,qualua,nlevsua,nobsua, &
! obsret,odifret,oanxret,iret, &
! qualret,nlevret,ncolret, &
! knt,wgtsum,zsum)
!
! 820 FORMAT(' Statistics for analysis pass ',i4,/ &
! ' ivar name knt bias rms')
! DO k=1,nvar
! WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
! END DO
! IF(iwstat > 0) THEN
! DO k=1,nvar
! WRITE(iwstat,812) k,nam_var(k), &
! knt(k,1),wgtsum(k,1),zsum(k,1)
! END DO
! WRITE(iwstat,806)
! END IF
! CALL diffnstra(nvar,nvarrad,nvarradin,nzrdr,mxrad,mxcolrad, &
! elvrad,distrad,hgtradc, &
! uazmrad,vazmrad,qualrad,irad, &
! obsrad,odifrad,oanxrad, &
! nlevrad,ncolrad, &
! knt,wgtsum,zsum)
! WRITE(6,815)
! DO k=1,nvar
! WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
! END DO
! WRITE(6,806)
! END IF
!
! 200 CONTINUE
!
istatus=0
!
RETURN
END SUBROUTINE minimization
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE LINEAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!
! Linear interplation in one direction.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! Jidong Gao, CAPS, July, 2000
!
!-----------------------------------------------------------------------
!
SUBROUTINE linear(n,xx,yy,u,f) 5
INTEGER :: n
REAL :: xx(n),yy(n),u,f
IF( u <= xx(1)) THEN
a1=( u-xx(2))/(xx(1)-xx(2))
a2=( u-xx(1))/(xx(2)-xx(1))
f=a1*yy(1)+a2*yy(2)
ELSE
DO i=2,n
IF( u <= xx(i) ) EXIT
END DO
i = min(i,n)
a1=( u-xx(i) )/(xx(i-1)-xx(i))
a2=( u-xx(i-1))/(xx(i)-xx(i-1))
f=a1*yy(i-1)+a2*yy(i)
END IF
! print*,"hqback==",xx
! print*,"qback==",yy
! stop
!
RETURN
END SUBROUTINE linear