subroutine assiwght(a,b,j,kz,tlapse,ifield) c c====================================================================== c === c This routine computes the melding weights A and B for the === c observations and model field IFIELD for the vertical slab === c J during OI assimilation. === c === c====================================================================== c c---------------------------------------------------------------------- c Define global data. c---------------------------------------------------------------------- c #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,ifield,j,k,kz,m #if defined fcsterr FLOAT * rseterr #elif !defined wgterr FLOAT * delta,errfct,errmax,errmin,range,ratio,wght #endif FLOAT * diag1,diag2,err,errobs,fac,tlapse FLOAT * a(imt,km),b(imt,km) c c======================================================================= c Begin executable code. c======================================================================= c #if defined fcsterr c----------------------------------------------------------------------- c Set forecast error increment. c----------------------------------------------------------------------- c fac = c2*( c1 - exp(-(tlapse/tsat)**2) ) c #endif c----------------------------------------------------------------------- c Get assimilation errors according to field type. c----------------------------------------------------------------------- c do 20 k=1,kz do 20 i=1,imt c c Transport stream function c if(ifield.eq.1) then if(iobserr.eq.0) then errobs=errp(icoi) else errobs=pobserr(i) endif c c Zonal component of internal velocity. c elseif(ifield.eq.2) then if(iobserr.eq.0) then errobs=erruv(icoi,k) else errobs=uobserr(i,k) endif c c Meridional component of internal velocity. c elseif(ifield.eq.3) then if(iobserr.eq.0) then errobs=erruv(icoi,k) else errobs=vobserr(i,k) endif c c Tracers. c elseif((4.le.ifield).and.(ifield.le.nt+3)) then m=ifield-3 if(iobserr.eq.0) then errobs=errts(icoi,k,m) else errobs=tobserr(i,k,m) endif endif #if defined fcsterr c c Transport stream function c if(ifield.eq.1) then err=pfcterr(i,j) c c Zonal component of internal velocity. c elseif(ifield.eq.2) then err=ufcterr(i,k,j) c c Meridional component of internal velocity. c elseif(ifield.eq.3) then err=vfcterr(i,k,j) c c Tracers. c elseif((4.le.ifield).and.(ifield.le.nt+3)) then err=tfcterr(i,k,m,j) endif c c----------------------------------------------------------------------- c Reset forecast error for model drift. c----------------------------------------------------------------------- c err = err + fac c c----------------------------------------------------------------------- c Reset observation error for ramping. c----------------------------------------------------------------------- c errobs = rseterr (errobs,err,cor,obswgt(icoi,ifield)) c c----------------------------------------------------------------------- c Compute assimilation coefficient. c----------------------------------------------------------------------- c diag1=cor*sqrt(errobs*err) diag2=errobs+err-c2*diag1 if(diag2.eq.c0) then a(i,k)=p5 else a(i,k)=(err-diag1)/diag2 endif #elif defined wgterr a(i,k)=(c1-min(c1,errobs))*obswgt(icoi,ifield) #else c c----------------------------------------------------------------------- c Compute forecast error from observation error using provided c observation weights. If variable observation errors, assign c that weigth to minimum error value and linearly interpolate c for other error values. Plus a great deal of boloni. c----------------------------------------------------------------------- c errmin=oerrmin(ifield,k) errmax=oerrmax(ifield,k) range=errmax-errmin delta=errmax-errobs ratio=c1 if(range.ne.c0) ratio=delta/range wght=ratio*obswgt(icoi,ifield) fac=(cor*(c1-c2*wght)+sqrt(c1-((c2*wght-c1)**2)*(c1-cor*cor))) * /(c2*(c1-wght)) errfct=fac*fac*errobs c c----------------------------------------------------------------------- c Compute weights for the melding of forecast and observation fields c----------------------------------------------------------------------- c diag1=cor*sqrt(errobs*errfct) diag2=errobs+errfct-2.*diag1 if(diag2.eq.c0) diag2=c1em3 a(i,k)=(errfct-diag1)/diag2 #endif if(a(i,k).gt.c1) a(i,k)=c1 if(a(i,k).lt.c0) a(i,k)=c0 b(i,k)=c1-a(i,k) c c----------------------------------------------------------------------- c Compute expected mean square error associated with assimilation c----------------------------------------------------------------------- #ifndef fcsterr c c Back out effective pre-assimilation forecast errors based on c observation error and assimilation weight. c diag1=c2*(c1-a(i,k)) diag2=(c1-c2*a(i,k))*cor if (diag1.ne.c0) then fac = (diag2 - sqrt(diag2*diag2+c2*diag1)) / diag1 err = errobs*fac*fac else err = errobs / (cor*cor) end if c c Compute expected mean square error associated with assimilation c diag1=cor*sqrt(errobs*err) diag2=errobs+err-c2*diag1 #endif c if(diag2.eq.c0) then err=c1em12 else err=(c1-cor*cor)*err*errobs/diag2 endif c c----------------------------------------------------------------------- c Update assimilation expected mean square error c----------------------------------------------------------------------- c if(iobserr.eq.0) then errfld(icoi+1,k,ifield)=err else #ifdef fcsterr if(ifield.eq.1) then pfcterr(i,j)=err elseif(ifield.eq.2) then ufcterr(i,k,j)=err elseif(ifield.eq.3) then vfcterr(i,k,j)=err else tfcterr(i,k,m,j)=err #else if(ifield.eq.1) then pfcterr(i)=err elseif(ifield.eq.2) then ufcterr(i,k)=err elseif(ifield.eq.3) then vfcterr(i,k)=err else tfcterr(i,k,m)=err #endif endif endif c c----------------------------------------------------------------------- c If error fields are uniform, and if requested, write out c assimilation information c----------------------------------------------------------------------- c if((ioiwrt.eq.1).and.(iobserr.eq.0)) then if((i.eq.1).and.(j.eq.1)) then write(ioiout,10) ifield,i,j,k,a(i,k),b(i,k),cor,errobs, #if defined fcsterr | defined wgterr * err #else * errfct,err #endif #if defined fcsterr | defined wgterr 10 format(i2,3i4,1x,2f8.5,3(1pe12.4)) #else 10 format(i2,3i4,1x,2f8.5,4(1pe12.4)) #endif endif endif 20 continue c return end