! !################################################################## !################################################################## !###### ###### !###### 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(j,i,5) ) 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