!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COSTF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Define the total costfunction for analysis
!
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! Jidong Gao, CAPS, July, 2000
!
!
SUBROUTINE costf(numctr,ctrv, cfun_single, & 3,21
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,qobsng2,qualsng,isrcsng,icatsng,nobsng, &
xua,yua,hgtua,theua, &
obsua,odifua,qobsua2,qualua,isrcua,nlevsua,nobsua, &
elvrad,xradc,yradc, &
distrad,uazmrad,vazmrad,hgtradc,theradc,dsdr,dhdr, &
obsrad,odifrad,qobsrad2,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,nradius_z,turn_div,cntl_var
REAL :: hradius
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 :: qobsng2(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 :: qobsua2(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 :: qobsrad2(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
!
!
!
! output:
! ------
! cfun: value of cost function
!
!-----------------------------------------------------------------------
!
! argument
! --------
INTEGER :: numctr,iflag,num
REAL :: ctrv(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 :: 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)
!additional memory in model space
!
REAL, DIMENSION (:,:,:), allocatable :: u,v,w,wcont
REAL :: zuu,zvv,zpp,ztt,zqq,zww
!
!
REAL :: cfun,f_b,f_bu,f_bv,f_bp,f_bt,f_bq,f_bw,ugrid,vgrid,vr
REAL :: wgrid,sum1
REAL :: f_osng,f_ousng,f_ovsng,f_opsng,f_otsng,f_oqsng,f_owsng
REAL :: f_oua,f_ouua,f_ovua,f_opua,f_otua,f_oqua,f_owua
REAL :: f_orad,f_ourad,f_ovrad,f_oprad,f_otrad,f_oqrad,f_owrad
REAL :: f_div
REAL :: div3(nx,ny,nz)
real :: maxdiv3, idiv,jdiv,kdiv
real :: tdiv1, tdiv2,tdiv3,tdiv4
! array save coast function of each variable, the nvar+1 is divergence and radar
REAL :: cfun_single(nvar+1), cfun_total
!
!
! allocate control variable arrays
! -------------------------------------------------------
!
allocate ( u(nx,ny,nz) )
allocate ( v(nx,ny,nz) )
allocate ( w(nx,ny,nz) )
allocate ( wcont(nx,ny,nz) )
! print*,'sngsw====',sngsw,uasw,radsw
!
DO i = 1, nvar+1
cfun_single(i) = 0.0
END DO
!
DO i = 1, nvar
DO j = 1, mxsng
obsng(i,j) = 0.0
END DO
END DO
!
DO i = 1, nvar
DO j = 1, nzua
DO k = 1, mxua
obsua(i,j,k) = 0.0
END DO
END DO
END DO
!
!
DO i = 1, nvarradin
DO j = 1, nzrdr
DO k = 1, mxcolrad
obsrad(i,j,k) = 0.0
END DO
END DO
END DO
!
!
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, psi, phi, p_ctr, &
t_ctr,q_ctr,w_ctr,ctrv)
!
!
f_b = 0.
!
f_bu = 0.
f_bv = 0.
f_bp = 0.
f_bt = 0.
f_bq = 0.
f_bw = 0.
!
DO k = 1,nz
DO j = 1,ny
DO i = 1,nx
zuu = psi(i,j,k) * psi(i,j,k)
zvv = phi(i,j,k) * phi(i,j,k)
zpp = p_ctr(i,j,k) * p_ctr(i,j,k)
ztt = t_ctr(i,j,k) * t_ctr(i,j,k)
zqq = q_ctr(i,j,k) * q_ctr(i,j,k)
zww = w_ctr(i,j,k) * w_ctr(i,j,k)
!
f_bu = f_bu + zuu
f_bv = f_bv + zvv
f_bp = f_bp + zpp
f_bt = f_bt + ztt
f_bq = f_bq + zqq
f_bw = f_bw + zww
END DO
END DO
END DO
!
f_b = f_bu+f_bv+f_bp+f_bt+f_bq+f_bw
IF( 1 == 2 ) THEN
f_b=0.0
cfun_single(1)=0.0
cfun_single(2)=0.0
cfun_single(3)=0.0
cfun_single(4)=0.0
cfun_single(5)=0.0
cfun_single(6)=0.0
ELSE
cfun_single(1)=f_bu
cfun_single(2)=f_bv
cfun_single(3)=f_bp
cfun_single(4)=f_bt
cfun_single(5)=f_bq
cfun_single(6)=f_bw
ENDIF
!
! print*,'f_b===',f_b
!
! --------------------------------------------------------
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdu_err, &
gdscal, psi)
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdv_err, &
gdscal, phi)
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdp_err, &
gdscal,p_ctr)
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdt_err, &
gdscal,t_ctr)
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdq_err, &
gdscal,q_ctr)
CALL ctr_to_vbl(ipass_filt,hradius,nradius_z,nx,ny,nz,gdw_err, &
gdscal,w_ctr)
!
!
! if (cntl_var ==0) option, u, v as control variables
!
!
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
!
! loading single level data
! --------------------------------------------------------
!
!
f_osng = 0.
f_ousng = 0.
f_ovsng = 0.
f_opsng = 0.
f_otsng = 0.
f_oqsng = 0.
f_owsng = 0.
!
!
! sum1=0
! do i = 1, nobsng
! sum1=sum1+odifsng(3,i)**2
! end do
! print*,'in Gcostsum1_odifsng==',sum1
! stop
!
print*,'sngsw===========',sngsw
!
IF(sngsw > 0) THEN
!
!
CALL linearint_3d(nx,ny,nz,u_ctr(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 linearint_3d(nx,ny,nz,v_ctr(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 linearint_3d(nx,ny,nz,p_ctr(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 linearint_3d(nx,ny,nz,t_ctr(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 linearint_3d(nx,ny,nz,q_ctr(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 linearint_3d(nx,ny,nz,w_ctr(1,1,1),xs(1),ys(1),zs(1,1,1),
!!! 1, mxsng, icatsng, nobsng, &
!!! 6,xsng_p,ysng_p,zsng_1,zsng_2,hgtsng,ihgtsng,corsng(1,6) )
!!!
!
!
num=0
DO i = 1, nobsng
iflag = 1
IF(qualsng(1,i) <= 0 .or. ihgtsng(i)<0 ) iflag = 0
zuu = corsng(i,1) - odifsng(1,i)
IF(iflag /= 0) THEN
corsng(i,1) = iflag*zuu / ( qobsng2(1,i) )
ELSE
corsng(i,1) = 0.0
END IF
zuu = zuu * corsng(i,1)
f_ousng = f_ousng + zuu
!
!
! if(iflag.eq.1) num=num+1
! print*,' num==',num,nobsng
! print*,'qobsng2==', qobsng2(1,i),odifsng(1,i),iflag
!
!
iflag = 1
IF(qualsng(2,i) <= 0 .or. ihgtsng(i)<0 ) iflag = 0
zvv = corsng(i,2) - odifsng(2,i)
IF(iflag /= 0) THEN
corsng(i,2) = iflag*zvv / ( qobsng2(2,i) )
ELSE
corsng(i,2) = 0.0
END IF
zvv = zvv * corsng(i,2)
f_ovsng = f_ovsng + zvv
!
! if(iflag.eq.1) num=num+1
! print*,' num==',num,nobsng
! print*,'corsng==', corsng(i,2),odifsng(2,i)
!
!
iflag = 1
IF(qualsng(3,i) <= 0 .or. ihgtsng(i)<0) iflag = 0
zpp = corsng(i,3) - odifsng(3,i)
IF(iflag /= 0) THEN
corsng(i,3) = iflag*zpp / ( qobsng2(3,i) )
ELSE
corsng(i,3) = 0.0
END IF
zpp = zpp * corsng(i,3)
f_opsng = f_opsng + zpp
!
! if(iflag.eq.1) num=num+1
! print*,' num==',num,nobsng
! print*,'corsng==', corsng(i,3),odifsng(3,i)
! print*,'qobsng2==', qobsng2(3,i),odifsng(3,i),iflag
! print*,'f_opsng====',qobsng2(1,i),qobsng2(2,i),
! : qobsng2(3,i),
! : qobsng2(4,i),qobsng2(5,i)
!
!
!
iflag = 1
IF(qualsng(4,i) <= 0 .or. ihgtsng(i)<0) iflag = 0
ztt = corsng(i,4) - odifsng(4,i)
IF(iflag /= 0) THEN
corsng(i,4) = iflag*ztt / ( qobsng2(4,i) )
ELSE
corsng(i,4) = 0.0
END IF
ztt = ztt * corsng(i,4)
f_otsng = f_otsng + ztt
!
! if(iflag.eq.1) num=num+1
! print*,' num==',num,nobsng
! print*,'corsng==', corsng(i,4),odifsng(4,i)
!
!
iflag = 1
IF(qualsng(5,i) <= 0 .or. ihgtsng(i)<0) iflag = 0
zqq = corsng(i,5) - odifsng(5,i)
IF(iflag /= 0) THEN
corsng(i,5) = iflag*zqq / ( qobsng2(5,i) )
ELSE
corsng(i,5) = 0.0
END IF
zqq = zqq * corsng(i,5)
f_oqsng = f_oqsng + zqq
!
!
!! if(iflag.eq.1) num=num+1
!! print*,'corsng==', corsng(i,5),odifsng(5,i)
!! print*,' qobsng2 odifsng=',qobsng2(5,i),odifsng(5,i)
! : ,qualsng(5,i),iflag
!
!
!
! if(1.eq.0) THEN
!
! iflag = 1
! if(qualsng(6,i).le.0 .and. ihgtsng<0) iflag = 0
! zww = corsng(i,6) - odifsng(6,i)
! if(iflag.eq.1) num=num+1
! print*,' num==',num,nobsng
! if(iflag.ne.0) then
! corsng(i,6) = iflag*zww / ( qobsng2(6,i) )
! else
! corsng(i,6) = 0.0
! end if
! zww = zww * corsng(i,6)
! f_owsng = f_owsng + zww
! print*,'obsng==', corsng(i,6),odifsng(6,i)
! end if
!
! observation cost function for single level data
! --------------------------------------------------------------
!
END DO
!
f_osng=f_ousng+f_ovsng+f_opsng+f_otsng+f_oqsng+f_owsng
cfun_single(1)=cfun_single(1)+f_ousng
cfun_single(2)=cfun_single(2)+f_ovsng
cfun_single(3)=cfun_single(3)+f_opsng
cfun_single(4)=cfun_single(4)+f_otsng
cfun_single(5)=cfun_single(5)+f_oqsng
cfun_single(6)=cfun_single(6)+f_owsng
!
print*,' npass=',npass,' f_ousng=', f_ousng
print*,' npass=',npass,' f_ovsng=', f_ovsng
print*,' npass=',npass,' f_opsng=', f_opsng
print*,' npass=',npass,' f_otsng=', f_otsng
print*,' npass=',npass,' f_oqsng=', f_oqsng
print*,' npass=',npass,' f_owsng=', f_owsng
print*,' npass=',npass,' f_osng =', f_osng
print*,' num==',num,nobsng
!
END IF
! stop
!
!
!
!
! loading upper level data
! --------------------------------------------------------
!
!
f_oua = 0.
f_ouua = 0.
f_ovua = 0.
f_opua = 0.
f_otua = 0.
f_oqua = 0.
f_owua = 0.
!
IF(uasw > 0) THEN
! print*,'nobsua nlevsua==========',nobsua,nlevsua(1)
! print*,'nobsua nlevsua==========',nobsua,nlevsua(2)
! print*,'nobsua nlevsua==========',nobsua,nlevsua(3)
! print*,'nobsua nlevsua==========',nobsua,nlevsua(4)
! print*,'nobsua nlevsua==========',nobsua,nlevsua(5)
! print*,'nobsua nlevsua==========',nobsua,nlevsua(6)
! stop
!
!
CALL linearint_3d(nx,ny,nz,u_ctr(1,1,1),xs(1),ys(1),zs(1,1,1), &
! CALL linearint_3d(nx,ny,nz,zp(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) )
if(1 == 1 ) THEN
CALL linearint_3d(nx,ny,nz,v_ctr(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 linearint_3d(nx,ny,nz,p_ctr(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 linearint_3d(nx,ny,nz,t_ctr(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 linearint_3d(nx,ny,nz,q_ctr(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) )
end if
!
! CALL linearint_3d(nx,ny,nz,w_ctr(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) )
!
!
do i = 1, nobsua
DO j = 1, nlevsua(i)
!
iflag = 1
IF(qualua(1,j,i) <= 0 .or. ihgtua(j,i)<0) iflag = 0
zuu = corua(j,i,1) - odifua(1,j,i)
! print*,'iflag=',iflag,corua(j,i,1),odifua(1,j,i),zuu
IF(iflag /= 0) THEN
corua(j,i,1) = iflag*zuu /( qobsua2(1,j,i) )
ELSE
corua(j,i,1) = 0.0
END IF
zuu = zuu * corua(j,i,1)
f_ouua = f_ouua + zuu
! goto 20000
!
! print*,'qobsua2=',qobsua2(1,j,i),qobsua2(2,j,i),
! : qobsua2(3,j,i),
! : qobsua2(4,j,i),
! : qobsua2(5,j,i)
!
! print*,'odifua==',odifua(1,j,i),odifua(2,j,i),
! : odifua(3,j,i),
! : odifua(4,j,i),
! : odifua(5,j,i)
!
!
!
iflag = 1
IF(qualua(2,j,i) <= 0 .or. ihgtua(j,i)<0) iflag = 0
zvv = corua(j,i,2) - odifua(2,j,i)
IF(iflag /= 0) THEN
corua(j,i,2) = iflag*zvv/( qobsua2(2,j,i) )
ELSE
corua(j,i,2) = 0.0
END IF
zvv = zvv * corua(j,i,2)
f_ovua = f_ovua + zvv
!
!
iflag = 1
IF(qualua(3,j,i) <= 0 .or. ihgtua(j,i)<0) iflag = 0
zpp = corua(j,i,3) - odifua(3,j,i)
IF(iflag /= 0) THEN
corua(j,i,3) = iflag*zpp /( qobsua2(3,j,i) )
ELSE
corua(j,i,3) = 0.0
END IF
zpp = zpp * corua(j,i,3)
f_opua = f_opua + zpp
!
!
iflag = 1
IF(qualua(4,j,i) <= 0 .or. ihgtua(j,i)<0) iflag = 0
ztt = corua(j,i,4) - odifua(4,j,i)
IF(iflag /= 0) THEN
corua(j,i,4) = iflag*ztt /( qobsua2(4,j,i) )
ELSE
corua(j,i,4) = 0.0
END IF
ztt = ztt * corua(j,i,4)
f_otua = f_otua + ztt
!
iflag = 1
IF(qualua(5,j,i) <= 0 .or. ihgtua(j,i)<0) iflag = 0
zqq = corua(j,i,5) - odifua(5,j,i)
IF(iflag /= 0) THEN
corua(j,i,5) = iflag*zqq / ( qobsua2(5,j,i) )
ELSE
corua(j,i,5) = 0.0
END IF
zqq = zqq * corua(j,i,5)
f_oqua = f_oqua + zqq
!
! print*,' qobsua2 odifua=',qobsua2(5,k,i),odifua(5,k,i)
! : ,qualua(5,k,i),iflag
!
! iflag = 1
! if(qualua(6,k,i).le.0 .or. ihgtua(j,i)<0) iflag = 0
! zww = corua(j,i,6) - odifua(6,j,i)
! if(iflag.ne.0) then
! corua(j,i,6) = iflag*zww / ( qobsua2(6,j,i) )
! else
! corua(j,i,6) = 0.0
! end if
! zww = zww * corua(j,i,6)
! f_owua = f_owua + zww
!
! 20000 continue
END DO
END DO
f_oua = f_ouua+f_ovua+f_opua+f_otua+f_oqua+f_owua
cfun_single(1)=cfun_single(1)+f_ouua
cfun_single(2)=cfun_single(2)+f_ovua
cfun_single(3)=cfun_single(3)+f_opua
cfun_single(4)=cfun_single(4)+f_otua
cfun_single(5)=cfun_single(5)+f_oqua
cfun_single(6)=cfun_single(6)+f_owua
!
print*,'f_oua = ',f_oua
print*,'f_ouua = ',f_ouua
print*,'f_ovua = ',f_ovua
print*,'f_opua = ',f_opua
print*,'f_otua = ',f_otua
print*,'f_oqua = ',f_oqua
print*,'f_owua = ',f_owua
print*,' num==',nobsua
! stop
!
END IF
!
print*,'f_oua=',f_ouua,f_ovua,f_opua,f_otua,f_oqua,f_oua
!
! print*, 'odifsng============'
! write(*,'(16f5.1)') ( odifsng(4,i),i = 1, nobsng )
! print*, 'obsua============'
! print*,'nobsua ==========',nobsua,f_oua
! stop
!
! do i = 1, 5
! write(*,*)
! write(*,'(I5,2x,a)') nlevsua(i),'u'
! write(*,'(15f7.1)') (odifua(1,j,i),j=1,nlevsua(i) )
! end do
! do i = 1, 5
! write(*,*)
! write(*,'(I5,2x,a)') nlevsua(i),'v'
! write(*,'(15f7.1)') (odifua(2,j,i),j=1,nlevsua(i) )
! end do
! do i = 1, 5
! write(*,*)
! write(*,'(I5,2x,a)') nlevsua(i),'p'
! write(*,'(15f7.1)') (odifua(3,j,i),j=1,nlevsua(i) )
! end do
! do i = 1, 5
! write(*,*)
! write(*,'(I5,2x,a)') nlevsua(i),'t'
! write(*,'(15f7.1)') (odifua(4,j,i),j=1,nlevsua(i) )
! end do
! do i = 1, 5
! write(*,*)
! write(*,'(I5,2x,a)') nlevsua(i),'q'
! write(*,'(15f7.1)') (odifua(5,j,i),j=1,nlevsua(i) )
! end do
!
! DO 401 i=1,nobsua
! write(*,*)
! write(*,*)'in cost obsua====',i
! DO 301 j=1,nlevsua(i)
! write(*,'(6f7.1)') hgtua(j,i),(odifua(k,j,i),k=1,5)
!301 continue
!401 continue
!
! loading radar data
! --------------------------------------------------------
!
!
f_orad = 0.
f_ovrad = 0.
f_oqrad = 0.
!
IF(radsw > 0) THEN
!
CALL linearint_3d(nx,ny,nz,u_ctr(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 linearint_3d(nx,ny,nz,v_ctr(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 linearint_3d(nx,ny,nz,w_ctr(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) )
!
! write(*,*) iuserad(1),iuserad(2)
!
DO i = 1, ncolrad
if(iuserad(isrcrad(irad(i))).GT.0) THEN
!
DO k=1,nlevrad(i)
!
! IF(qualrad(2,k,i) > 0 ) THEN
IF(qualrad(2,k,i) > 0 .and.ihgtradc(k,i)>=0) THEN
!
vr= ( uazmrad(i)*corrad(k,i,1) &
+vazmrad(i)*corrad(k,i,2) ) * dsdr(k,i) &
+corrad(k,i,3) * dhdr(k,i)
!
zvv = vr - odifrad(2,k,i)
!
! corrad(k,i,2)=zvv /( qobsrad2(2,k,i) )
! zvv = zvv * corrad(k,i,2)
!
obsrad(2,k,i)=zvv /( qobsrad2(2,k,i) )
zvv = zvv * obsrad(2,k,i)
f_ovrad = f_ovrad + zvv
!
!!
!! print*,'sqrt=',qobsrad2(1,k,i),qobsrad2(2,k,i)
!!
ENDIF ! (qualrad(2,k,i) > 0 )
!!
!! IF(1 == 0) THEN
!!
!! CALL linearint_3d(nx,ny,nz,q_ctr(1,1,1),xs(1),ys(1),zs(1,1,1), &
!! 5,xradc(i),yradc(i),hgtradc(k,i),obsrad(1,k,i),iflag )
!! IF(qualrad(3,k,i) <= 0) iflag = 0
!!_gao zqq = obsrad(1,k,i) - odifrad(1,k,i)
!! zqq = obsrad(1,k,i) - odifrad(1,k,i)
!! IF(iflag /= 0) THEN
!!_gao obsrad(1,k,i)=iflag*zqq /(qobsrad2(1,k,i)*qobsrad2(1,k,i))
!! obsrad(1,k,i)=iflag*zqq /( qobsrad2(1,k,i) )
!! ELSE
!! obsrad(1,k,i)=0.0
!! END IF
!! zqq = zqq * obsrad(1,k,i)
!! f_oqrad = f_oqrad + zqq
!! END IF
!!
END DO ! nlevrad(i)
!!
ENDIF
END DO
END IF
print*,'f_ovrad==',f_ovrad
! stop
!-----------------------------------------------------------------------
!
! Compute the momentum divergence term, defined as
!
! div = d(u*rhostr)/dx + d(v*rhostr)/dy + d(wcont*rhostr)/dz.
!
!-----------------------------------------------------------------------
!
!
f_div = 0.0
if( turn_div == 1 ) then
!
DO k= 1,nz
DO j= 1,ny
DO i= 1,nx
u (i,j,k) = anx(i,j,k,1)+u_ctr(i,j,k)
v (i,j,k) = anx(i,j,k,2)+v_ctr(i,j,k)
w (i,j,k) = anx(i,j,k,6)+w_ctr(i,j,k)
END DO
END DO
END DO
!
CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
!
IF(wei_div_h > 0.0 ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6)
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
tem3(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k)
END DO
END DO
END DO
ENDIF
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
div3(i,j,k) = 1./wei_div_h * j3inv(i,j,k) &
* ( mapfct(i,j,7) &
* ((tem1(i+1,j,k)-tem1(i,j,k))*dxinv &
+(tem2(i,j+1,k)-tem2(i,j,k))*dyinv)) &
+1./wei_div_v * j3inv(i,j,k) &
* ((tem3(i,j,k+1)-tem3(i,j,k))*dzinv )
f_div = f_div + div3(i,j,k) * div3(i,j,k)
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
div3(i,j,k) = 1./wei_div_h * j3inv(i,j,k) &
* ( mapfct(i,j,7) &
* ((tem1(i+1,j,k)-tem1(i,j,k))*dxinv &
+(tem2(i,j+1,k)-tem2(i,j,k))*dyinv))
f_div = f_div + div3(i,j,k) * div3(i,j,k)
END DO
END DO
END DO
ENDIF
ENDIF
f_orad = f_ovrad+f_oqrad+f_div
cfun_single(nvar+1)=f_orad
write(*,*) 'f_orad =', f_ovrad,f_oqrad,f_div
!
! print*,'f_orad==',f_orad
!
!mhu END IF
!
!
cfun = (f_b+f_osng+f_oua+f_orad)/2.0
DO i=1,nvar+1
cfun_single(i)=cfun_single(i)/2.0
ENDDO
cfun_total=0
DO i=1,nvar+1
cfun_total=cfun_total+cfun_single(i)
ENDDO
! write(*,*) 'check cost function==',cfun_total, cfun
!
! PRINT*,' cfun b o===',f_b,f_osng,f_oua,f_orad,f_div
! PRINT*,' cfun b o===',f_b,2*cfun-f_b, f_div
!
print*,' cfun===',cfun,' fb=',f_b,' f_osng=',f_osng &
,' f_oua==',f_oua,' f_orad=',f_orad
!
! if (icall .eq. 0) then
! write (911, '(/,a,a,a,a)')
! $ ' call u+v+t+q+p u v',
! $ ' t q',
! $ ' p'
! endif
!
! write (911, '(2x,i3,2x,1pe15.8,2x,a,2x,5(1pe15.8,2x))')
! $ icall, f_b/2., 'f_b', f_bu/2., f_bv/2., f_bt/2.,
! $ f_bq/2., f_bp/2.
!
! write (911, '(2x,i3,2x,1pe15.8,2x,a,2x,5(1pe15.8,2x))')
! $ icall, f_o/2., 'f_o', f_ou/2., f_ov/2., f_ot/2.,
! $ f_oq/2., f_op/2.
!
! write (911, '(2x,i3,2x,1pe15.8,2x,a,/))')
! $ icall, cfun,'f_cost'
!
!
deallocate( u )
deallocate( v )
deallocate( w )
deallocate( wcont )
!
!
! ------
!
RETURN
END SUBROUTINE costf