!******************************************************************************* subroutine bias_rmse_cv_mpi(cgrid,itotmodel,imodel,cmodel,rmusecv,rmuseb, & lpadjust,fnamesuffix) !*********************************************************************** ! abstract: compute domain averaged bias and rmse for background, * ! interim analysis, and analysis. Also compute domain * ! averaged cross-validation scores * ! * ! program history log: * ! 2005-10-08 pondeca * !*********************************************************************** use mpi use mpitaskmod, only: mype,npe,ista,iend,jsta,jend use mpitaskmod, only: nx => nlon use mpitaskmod, only: ny => nlat use controlvars, only: igust,ivis,ipblh,idist, & iwspd10m,itd2m,imxtm,imitm,ipmsl,ihowv,itcamt,ilcbas, & iuwnd10m,ivwnd10m use cressanl_common, only: nobsmax, xlocs, ylocs, hgt0s, & hobs, rmuses, & rtvts, cstations, & insubdoms, dups, bckgs, xberrs, & ipointer,uvpointer,nmax,nflds, & kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, & ktds,kgusts,kvis,kpblhs,kdists, & kwspd10ms,ktd2ms,kmxtms,kmitms,kpmsls,khowvs,ktcamts,klcbas, & kuwnd10ms,kvwnd10ms,kaltws,kaltw2s,kaltwds implicit none ! note: consider mpi version ! ! output ! bias(1) , rmse(1) - psfc bias and rmse ! bias(2) , rmse(2) - 2m t bias and rmse ! bias(3) , rmse(3) - 2m q bias and rmse ! bias(4) , rmse(4) - 10m u bias and rmse ! bias(5) , rmse(5) - 10m v bias and rmse ! bias(6) , rmse(6) - 10m wind bias and rmse ! bias(7) , rmse(7) - 10m wind2 bias and rmse ! bias(8) , rmse(8) - 10m wdir bias and rmse ! bias(9) , rmse(9) - 2m td bias and rmse of ! bias(10) , rmse(10) - 10m gust bias and rmse ! bias(11) , rmse(11) - vis bias and rmse ! bias(12) , rmse(12) - pblh bias and rmse ! bias(13) , rmse(13) - dist bias and rmse ! bias(14) , rmse(14) - wspd20m bias and rmse ! bias(15) , rmse(15) - td2m bias and rmse ! bias(16) , rmse(16) - mxtm bias and rmse ! bias(17) , rmse(17) - mitm bias and rmse ! bias(18) , rmse(18) - pmls bias and rmse ! bias(19) , rmse(19) - howv bias and rmse ! bias(20) , rmse(20) - tcamt bias and rmse ! bias(21) , rmse(21) - lcbas bias and rmse ! bias(22) , rmse(22) - uwnd10m bias and rmse ! bias(23) , rmse(23) - vwnd10m bias and rmse ! bias(24) , rmse(24) - altw bias and rmse ! alt stands for "alternate" ! bias(25) , rmse(25) - altw2 bias and rmse ! bias(26) , rmse(26) - altwd bias and rmse ! ! note: order for nflds is psfc,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist, ! wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas,uwnd10m,vwnd10m,altw,altw2,,altwd !------------------------------------------------------------------------------ ! Declare passed variables integer(4),intent(in):: itotmodel integer(4),intent(in)::imodel real(4),intent(in)::rmusecv real(4),intent(in)::rmuseb logical,intent(in)::lpadjust character(60),intent(in)::cgrid character(60),intent(in)::cmodel character(*),intent(in)::fnamesuffix ! Declare local parameters integer(4),parameter::nflds0=26 !same as nflds integer(4),parameter::itotmodel0=5 !if itmodel > itmodel0 must increase itotmodel0 real(4),parameter :: gravity=9.81 real(4),parameter :: epsilon=1.e-03 real(8),parameter :: r360_d=360._8 real(4),parameter::spval=-999. ! Declare local variables integer(4),allocatable,dimension(:):: i1_info integer(4),allocatable,dimension(:):: i2_info integer(4),allocatable,dimension(:):: j1_info integer(4),allocatable,dimension(:):: j2_info integer(4),allocatable,dimension(:):: iaux1,iaux2,iaux3,iaux4 integer(4) i1,i2,j1,j2 integer(4) m1,m2,n1,n2 integer(4) mype1 integer(4) i,j,k,kk,n,nn,iu,iv,iu0,iv0 integer(4) ierror,islab,nobsmax0 integer(4),allocatable,dimension(:,:):: ipointer0 real(4),allocatable,dimension(:):: xlocs0, ylocs0,hobs0,rmuses0,rtvts0,hgt0s0 logical,allocatable,dimension(:):: insubdoms0 integer(4), save:: allnobs(nflds0,3,itotmodel0) real(4), save:: allbias(nflds0,3,itotmodel0) !the 3 accounts for rmuse=rmusecv.,1., and (rmusecv .or. 1.). real(4), save:: allrmse(nflds0,3,itotmodel0) real(4), save:: allabse(nflds0,3,itotmodel0) character(60), save:: allcmodel(itotmodel0) integer(4) nobs0(nflds) real(4) bias0(nflds) real(4) rmse0(nflds) real(4) abse0(nflds) character(60) cmodel0 integer(4) nobs(nflds,3) real(4) bias(nflds,3) real(4) rmse(nflds,3) real(4) abse(nflds,3) real(8) s1(nflds,3) real(8) s2(nflds,3) real(8) s3(nflds,3) real(8) s4(nflds,3) integer(4) nobs99(nflds,3) real(8) s99(nflds,3) integer(4) rtime(6),nlon,nlat,nsig real(4),allocatable,dimension(:,:)::glon real(4),allocatable,dimension(:,:)::fis1,psfc1,t1,q1,u1,v1,w1,wd1,td1 real(4),allocatable,dimension(:,:)::gust1,vis1,pblh1,dist1 real(4),allocatable,dimension(:,:)::wspd10m1,td2m1,mxtm1,mitm1,pmsl1,howv1,tcamt1,lcbas1 real(4),allocatable,dimension(:,:)::uwnd10m1,vwnd10m1 real(4),allocatable,dimension(:,:)::uwnd10me,vwnd10me real(4),allocatable,dimension(:,:)::altw1,altwd1 real(4),allocatable,dimension(:,:)::tv1 real(4),allocatable,dimension(:,:)::ue,ve real(4),allocatable,dimension(:,:)::fis0 real(4),allocatable,dimension(:,:)::tv0 integer(4),allocatable,dimension(:,:):: iauxfield real(4),allocatable,dimension(:,:)::auxfield,auxfield2 real(4),allocatable,dimension(:,:,:)::field0 integer(4) npts integer(4) kobs0(nflds) real(4) amin0,amin,amax0,amax real(4) xx,yy,ob,rmuse real(4) rtvts9,zsges,tges real(4) fint real(4) uob,uob_model,vob,vob_model,waux real(8) diff,diff2,du2,dv2 real(8) zeroone(3) character(60) fname character(10) chvar logical fexist logical lstats1,lstats2,lstats3 logical lanlerr !******************************************************************************* if (mype.eq.0) print*,'in bias_rmse_cv:nflds0,nflds=',nflds0,nflds if (nflds0 .ne. nflds) then print*,'in bias_rmse_cv: trouble!!! nfld0 and nfld must be identical' call mpi_finalize(ierror) stop endif if (itotmodel > itotmodel0) then if (mype.eq.0) then print*,'in bias_rmse_cv:' print*,'itotmodel,itotmodel0=',itotmodel,itotmodel0 print*,' itotmodel must be <= itotmodel0. ... returning' endif return endif if (trim(cgrid)/='conus' .and. trim(cgrid)/='alaska' .and. & trim(cgrid)/='prico' .and. trim(cgrid)/='hawaii' .and. & trim(cgrid)/='cohres' .and. trim(cgrid)/='akhres' .and. & trim(cgrid)/='hrrr' .and. trim(cgrid)/='guam' .and. & trim(cgrid)/='juneau' .and. trim(cgrid)/='cohresext' .and. & trim(cgrid)/='cohreswexp' ) then if (mype.eq.0) print*,& 'in bias_rmse_cv: unknown grid ', cgrid, 'returning ...' return endif if(mype.eq.0) then print*,'in bias_rmse_cv:: cgrid,nx,ny=',trim(cgrid),nx,ny print*,'in bias_rmse_cv:: cmodel,rmusecv,rmuseb=',trim(cmodel),rmusecv,rmuseb print*,'in bias_rmse_cv:: fnamesuffix=',trim(fnamesuffix) endif allocate(i1_info(npe)) allocate(i2_info(npe)) allocate(j1_info(npe)) allocate(j2_info(npe)) allocate(iaux1(npe)) allocate(iaux2(npe)) allocate(iaux3(npe)) allocate(iaux4(npe)) i1_info(:)=0 ; iaux1(:)=0 i2_info(:)=0 ; iaux2(:)=0 j1_info(:)=0 ; iaux3(:)=0 j2_info(:)=0 ; iaux4(:)=0 i1=max(1,ista-1) ; iaux1(mype+1)=i1 i2=min(nx,iend+1) ; iaux2(mype+1)=i2 j1=max(1,jsta-1) ; iaux3(mype+1)=j1 j2=min(ny,jend+1) ; iaux4(mype+1)=j2 call mpi_allreduce(iaux1,i1_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux2,i2_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux3,j1_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux4,j2_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) allocate ( iauxfield (nx,ny) ) allocate ( auxfield (nx,ny) ) allocate ( auxfield2 (nx,ny) ) allocate ( glon (i1:i2,j1:j2) ) allocate ( psfc1 (i1:i2,j1:j2) ) allocate ( fis1 (i1:i2,j1:j2) ) allocate ( t1 (i1:i2,j1:j2) ) allocate ( tv1 (i1:i2,j1:j2) ) allocate ( q1 (i1:i2,j1:j2) ) allocate ( u1 (i1:i2,j1:j2) ) ; allocate ( ue(i1:i2,j1:j2) ) allocate ( v1 (i1:i2,j1:j2) ) ; allocate ( ve(i1:i2,j1:j2) ) allocate ( w1 (i1:i2,j1:j2) ) allocate ( wd1 (i1:i2,j1:j2) ) allocate ( td1 (i1:i2,j1:j2) ) allocate ( gust1 (i1:i2,j1:j2) ) allocate ( vis1 (i1:i2,j1:j2) ) allocate ( pblh1 (i1:i2,j1:j2) ) allocate ( dist1 (i1:i2,j1:j2) ) allocate ( wspd10m1 (i1:i2,j1:j2) ) allocate ( td2m1 (i1:i2,j1:j2) ) allocate ( mxtm1 (i1:i2,j1:j2) ) allocate ( mitm1 (i1:i2,j1:j2) ) allocate ( pmsl1 (i1:i2,j1:j2) ) allocate ( howv1 (i1:i2,j1:j2) ) allocate ( tcamt1(i1:i2,j1:j2) ) allocate ( lcbas1(i1:i2,j1:j2) ) allocate ( uwnd10m1 (i1:i2,j1:j2) ) ; allocate ( uwnd10me (i1:i2,j1:j2) ) allocate ( vwnd10m1 (i1:i2,j1:j2) ) ; allocate ( vwnd10me (i1:i2,j1:j2) ) allocate ( altw1 (i1:i2,j1:j2) ) allocate ( altwd1 (i1:i2,j1:j2) ) fis1=0. gust1=0. vis1=0. pblh1=0. dist1=0. wspd10m1=0. td2m1=0. mxtm1=0. mitm1=0. pmsl1=0. howv1=0. tcamt1=0. lcbas1=0. uwnd10m1=0. vwnd10m1=0. altw1=0. altwd1=0. lanlerr=trim(cmodel)=='errfield.dat'.or.trim(cmodel)=='errfield.dat_precval' open (52,file=trim(cmodel),form='unformatted') if (.not.lanlerr) then read(52) rtime,nlon,nlat,nsig read(52) auxfield,auxfield2 !glat,dx read(52) auxfield,auxfield2 !glon,dy if (mype.eq.0) print*,'in bias_rmse_cv nlon,nlat,nsig=',nlon,nlat,nsig glon (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; psfc1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; fis1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; t1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; q1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; u1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; v1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) nn=12 if (trim(cmodel)=='sigges') nn=11 do n=1,nn !must jump nn records to get to gust1 read(52) enddo read(52) auxfield ; gust1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; vis1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; pblh1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; dist1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; wspd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; td2m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; mxtm1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; mitm1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; pmsl1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; howv1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; tcamt1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; lcbas1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; uwnd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) read(52) auxfield ; vwnd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) do j=j1,j2 do i=i1,i2 tv1(i,j)=t1(i,j)+(1.+0.608*max(0.,q1(i,j))) enddo enddo call wspdwdir(cgrid,glon,u1,v1,ue,ve,w1,wd1,i1,i2,j1,j2,mype) call wspdwdir(cgrid,glon,uwnd10m1,vwnd10m1,uwnd10me,vwnd10me,altw1,altwd1,i1,i2,j1,j2,mype) inquire(file='td.dat_'//trim(cmodel),exist=fexist) if (fexist) then open(57,file='td.dat_'//trim(cmodel),form='unformatted') read (57) auxfield ; td1(i1:i2,j1:j2)=auxfield(i1:i2,j1:j2) close(57) else call get_dewpt(psfc1,q1,t1,td1,i1,i2,j1,j2) endif endif if (lanlerr) then read(52) auxfield ; psfc1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !perr read(52) auxfield ; t1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !terr read(52) auxfield ; td1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !tderr read(52) auxfield ; u1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !uerr read(52) auxfield ; v1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !verr read(52) auxfield ; q1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !qerr read(52) auxfield ; wd1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !wdirerr2 read(52) auxfield ; w1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !wspderr read(52) auxfield ; ue (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !uerr2 read(52) auxfield ; ve (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !verr2 read(52) auxfield !wdirerr if (igust > 0 ) read(52) auxfield ; gust1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !gusterr if (ivis > 0 ) read(52) auxfield ; vis1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !viserr if (ipblh > 0 ) read(52) auxfield ; pblh1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !pblherr if (idist > 0 ) read(52) auxfield ; dist1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !disterr if (iwspd10m > 0 ) read(52) auxfield ; wspd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !wspd10merr if (itd2m > 0 ) read(52) auxfield ; td2m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !td2merr if (imxtm > 0 ) read(52) auxfield ; mxtm1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !mxtmerr if (imitm > 0 ) read(52) auxfield ; mitm1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !mitmerr if (ipmsl > 0 ) read(52) auxfield ; pmsl1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !pmslerr if (ihowv > 0 ) read(52) auxfield ; howv1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !howverr if (itcamt > 0 ) read(52) auxfield ; tcamt1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !tcamterr if (ilcbas > 0 ) read(52) auxfield ; lcbas1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !lcbaserr if (iuwnd10m > 0 ) read(52) auxfield ; uwnd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !uwnd10merr if (ivwnd10m > 0 ) read(52) auxfield ; vwnd10m1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !vwnd10merr if (iuwnd10m > 0 .and. ivwnd10m > 0) then !use this for now /MPondeca 20Oct20116 uwnd10me (i1:i2,j1:j2)=uwnd10m1(i1:i2,j1:j2) vwnd10me (i1:i2,j1:j2)=vwnd10m1(i1:i2,j1:j2) altw1 (i1:i2,j1:j2)=w1(i1:i2,j1:j2) altwd1 (i1:i2,j1:j2)=wd1(i1:i2,j1:j2) endif tv1(i1:i2,j1:j2)=t1(i1:i2,j1:j2) endif close(52) if (mype==0 ) print*,'*************in bias_rmse_cv********************' do nn=1,26 if ( nn==11 .and. igust <= 0 ) cycle if ( nn==12 .and. ivis <= 0 ) cycle if ( nn==13 .and. ipblh <= 0 ) cycle if ( nn==14 .and. idist <= 0 ) cycle if ( nn==15 .and. iwspd10m <= 0 ) cycle if ( nn==16 .and. itd2m <= 0 ) cycle if ( nn==17 .and. imxtm <= 0 ) cycle if ( nn==18 .and. imitm <= 0 ) cycle if ( nn==19 .and. ipmsl <= 0 ) cycle if ( nn==20 .and. ihowv <= 0 ) cycle if ( nn==21 .and. itcamt <= 0 ) cycle if ( nn==22 .and. ilcbas <= 0 ) cycle if ( nn==23 .and. iuwnd10m <= 0 ) cycle if ( nn==24 .and. ivwnd10m <= 0 ) cycle if ( nn==25 .and. (iuwnd10m <= 0 .or. ivwnd10m <= 0 ) ) cycle if ( nn==26 .and. (iuwnd10m <= 0 .or. ivwnd10m <= 0 ) ) cycle if ( nn == 1 ) then ; amin0=minval(psfc1) ; amax0=maxval(psfc1) ; chvar='psfc1' ; endif if ( nn == 2 ) then ; amin0=minval(fis1) ; amax0=maxval(fis1) ; chvar='fis1' ; endif if ( nn == 3 ) then ; amin0=minval(t1) ; amax0=maxval(t1) ; chvar='t1' ; endif if ( nn == 4 ) then ; amin0=minval(tv1) ; amax0=maxval(tv1) ; chvar='tv1' ; endif if ( nn == 5 ) then ; amin0=minval(td1) ; amax0=maxval(td1) ; chvar='td1' ; endif if ( nn == 6 ) then ; amin0=minval(ue) ; amax0=maxval(ue) ; chvar='ue' ; endif if ( nn == 7 ) then ; amin0=minval(ve) ; amax0=maxval(ve) ; chvar='ve' ; endif if ( nn == 8 ) then ; amin0=minval(w1) ; amax0=maxval(w1) ; chvar='w1' ; endif if ( nn == 9 ) then ; amin0=minval(wd1) ; amax0=maxval(wd1) ; chvar='wd1' ; endif if ( nn == 10 ) then ; amin0=minval(q1) ; amax0=maxval(q1) ; chvar='q1' ; endif if ( nn == 11 ) then ; amin0=minval(gust1) ; amax0=maxval(gust1) ; chvar='gust1' ; endif if ( nn == 12 ) then ; amin0=minval(vis1) ; amax0=maxval(vis1) ; chvar='vis1' ; endif if ( nn == 13 ) then ; amin0=minval(pblh1) ; amax0=maxval(pblh1) ; chvar='pblh1' ; endif if ( nn == 14 ) then ; amin0=minval(dist1) ; amax0=maxval(dist1) ; chvar='dist1' ; endif if ( nn == 15 ) then ; amin0=minval(wspd10m1) ; amax0=maxval(wspd10m1);chvar='wspd10m1' ; endif if ( nn == 16 ) then ; amin0=minval(td2m1) ; amax0=maxval(td2m1) ; chvar='td2m1' ; endif if ( nn == 17 ) then ; amin0=minval(mxtm1) ; amax0=maxval(mxtm1) ; chvar='mxtm1' ; endif if ( nn == 18 ) then ; amin0=minval(mitm1) ; amax0=maxval(mitm1) ; chvar='mitm1' ; endif if ( nn == 19 ) then ; amin0=minval(pmsl1) ; amax0=maxval(pmsl1) ; chvar='pmsl1' ; endif if ( nn == 20 ) then ; amin0=minval(howv1) ; amax0=maxval(howv1) ; chvar='howv1' ; endif if ( nn == 21 ) then ; amin0=minval(tcamt1) ; amax0=maxval(tcamt1); chvar='tcamt1'; endif if ( nn == 22 ) then ; amin0=minval(lcbas1) ; amax0=maxval(lcbas1); chvar='lcbas1'; endif if ( nn == 23 ) then ; amin0=minval(uwnd10me) ; amax0=maxval(uwnd10m1);chvar='uwnd10m1' ; endif if ( nn == 24 ) then ; amin0=minval(vwnd10me) ; amax0=maxval(vwnd10m1);chvar='vwnd10m1' ; endif if ( nn == 25 ) then ; amin0=minval(altw1) ; amax0=maxval(altw1) ;chvar='altw1' ; endif if ( nn == 26 ) then ; amin0=minval(altwd1) ; amax0=maxval(altwd1) ;chvar='altwd1' ; endif 150 continue call mpi_allreduce(amin0,amin,1,mpi_real,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax0,amax,1,mpi_real,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'in bias_rmse_cv:', trim(chvar),' min,max=',amin,amax enddo if (mype==0) print*,'************************************************' mype1=mype+1 if (imodel==1) then allnobs(:,:,:)=9999 allbias(:,:,:)=9999.000 allrmse(:,:,:)=9999.000 allabse(:,:,:)=9999.000 endif nobs(:,:)=0 s1(:,:)=0._8 s2(:,:)=0._8 s3(:,:)=0._8 s4(:,:)=0._8 fis1=fis1/gravity !terrain field do 300 islab=1,npe m1=i1_info(islab) !all processors operate on the m2=i2_info(islab) !same subdomain (m1:m2,n1:n2) n1=j1_info(islab) n2=j2_info(islab) allocate(field0 (m1:m2,n1:n2,nflds)) allocate(fis0 (m1:m2,n1:n2)) allocate(tv0 (m1:m2,n1:n2)) call mpi_barrier(mpi_comm_world,ierror) if (mype1==islab) then do j=n1,n2 do i=m1,m2 field0 (i,j,1) = psfc1(i,j) field0 (i,j,2) = t1 (i,j) field0 (i,j,3) = q1 (i,j) field0 (i,j,4) = ue (i,j) field0 (i,j,5) = ve (i,j) field0 (i,j,6) = w1 (i,j) !use same field0 (i,j,7) = w1 (i,j) !use same field0 (i,j,8) = wd1 (i,j) field0 (i,j,9) = td1 (i,j) field0 (i,j,10) = gust1(i,j) field0 (i,j,11) = vis1 (i,j) field0 (i,j,12) = pblh1(i,j) field0 (i,j,13) = dist1(i,j) field0 (i,j,14) = wspd10m1(i,j) field0 (i,j,15) = td2m1(i,j) field0 (i,j,16) = mxtm1(i,j) field0 (i,j,17) = mitm1(i,j) field0 (i,j,18) = pmsl1(i,j) field0 (i,j,19) = howv1(i,j) field0 (i,j,20) = tcamt1(i,j) field0 (i,j,21) = lcbas1(i,j) field0 (i,j,22) = uwnd10me(i,j) field0 (i,j,23) = vwnd10me(i,j) field0 (i,j,24) = altw1(i,j) !use same field0 (i,j,25) = altw1(i,j) !use same field0 (i,j,26) = altwd1(i,j) fis0 (i,j) = fis1(i,j) tv0 (i,j) = tv1 (i,j) enddo enddo nobsmax0=nobsmax kobs0(1)=kps kobs0(2)=kts kobs0(3)=kqs kobs0(4)=kus kobs0(5)=kvs kobs0(6)=kws kobs0(7)=kw2s kobs0(8)=kwds kobs0(9)=ktds kobs0(10)=kgusts kobs0(11)=kvis kobs0(12)=kpblhs kobs0(13)=kdists kobs0(14)=kwspd10ms kobs0(15)=ktd2ms kobs0(16)=kmxtms kobs0(17)=kmitms kobs0(18)=kpmsls kobs0(19)=khowvs kobs0(20)=ktcamts kobs0(21)=klcbas kobs0(22)=kuwnd10ms kobs0(23)=kvwnd10ms kobs0(24)=kaltws kobs0(25)=kaltw2s kobs0(26)=kaltwds endif call mpi_barrier(mpi_comm_world,ierror) npts=(n2-n1+1)*(m2-m1+1)*nflds call mpi_bcast (field0, npts, mpi_real, islab-1, mpi_comm_world, ierror) npts=(n2-n1+1)*(m2-m1+1) call mpi_bcast (fis0, npts, mpi_real, islab-1, mpi_comm_world, ierror) npts=(n2-n1+1)*(m2-m1+1) call mpi_bcast (tv0, npts, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (nobsmax0, 1, mpi_integer, islab-1, mpi_comm_world, ierror) call mpi_bcast (kobs0, nflds, mpi_integer, islab-1, mpi_comm_world, ierror) allocate(xlocs0(nobsmax0)) allocate(ylocs0(nobsmax0)) allocate(hobs0(nobsmax0)) allocate(rmuses0(nobsmax0)) allocate(rtvts0(nobsmax0)) allocate(hgt0s0(nobsmax0)) allocate(insubdoms0(nobsmax0)) allocate(ipointer0(nmax,nflds)) call mpi_barrier(mpi_comm_world,ierror) if (mype1==islab) then xlocs0(:)=xlocs(:) ylocs0(:)=ylocs(:) hobs0(:)=hobs(:) rmuses0(:)=rmuses(:) rtvts0(:)=rtvts(:) hgt0s0(:)=hgt0s(:) insubdoms0(:)=insubdoms(:) ipointer0(:,:)=ipointer(:,:) endif call mpi_barrier(mpi_comm_world,ierror) call mpi_bcast (xlocs0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (ylocs0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (hobs0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (rmuses0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (rtvts0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (hgt0s0, nobsmax0, mpi_real, islab-1, mpi_comm_world, ierror) call mpi_bcast (insubdoms0, nobsmax0, mpi_logical, islab-1, mpi_comm_world, ierror) npts=nmax*nflds call mpi_bcast (ipointer0, npts, mpi_integer, islab-1, mpi_comm_world, ierror) do 200 nn=1,nflds do 100 kk=1,kobs0(nn) if (mype==mod(kk-1,npe)) then n=ipointer0(kk,nn) if (.not.insubdoms0(n)) cycle xx=xlocs0(n) yy=ylocs0(n) ob=hobs0(n) rmuse=float(int(rmuses0(n))) !rmuses0(n) IF (NN==20 .AND. OB < 0.) CYCLE !must fix gsi for tcamt and remove this statement call bilinear_2d0v2(field0(m1:m2,n1:n2,nn),m1,m2,n1,n2,fint,yy,xx) if (nn==1) then if (lpadjust .and. .not.lanlerr) then call bilinear_2d0v2(fis0,m1,m2,n1,n2,zsges,yy,xx) rtvts9=1. !assume this for now if (rtvts9==0.) then call bilinear_2d0v2(tv0,m1,m2,n1,n2,tges,yy,xx) else call bilinear_2d0v2(field0(m1:m2,n1:n2,2),m1,m2,n1,n2,tges,yy,xx) endif call adjust_pressure(zsges,tges,hgt0s0(n),spval,spval,fint) endif fint=fint/100. !convert to hPa ob=ob/100. !convert to hPa endif if (nn==2) then if(rtvts0(n)==0. .and. .not.lanlerr) & call bilinear_2d0v2(tv0,m1,m2,n1,n2,fint,yy,xx) endif if (nn==3) then fint=fint*1000. !use g/kg ob=ob*1000. !use g/kg endif if (nn==8 .or. nn==26) then !Recompute fint for wind direction. Straight bilinear interpolation won't do. if (nn==8) then !Instead, compute wind direction from the interpolated u- and v values iu0=4 ; iv0=5 endif if (nn==26) then iu0=22 ; iv0=23 endif call bilinear_2d0v2(field0(m1:m2,n1:n2,iu0),m1,m2,n1,n2,uob_model,yy,xx) call bilinear_2d0v2(field0(m1:m2,n1:n2,iv0),m1,m2,n1,n2,vob_model,yy,xx) call speeddir(uob_model,vob_model,fint,waux) endif if (lanlerr) then diff=dble(fint) else diff=dble(fint-ob) if ( (nn==8.or.nn==26) .and. abs(diff) > 0._8) then !!! diff2=min(abs(diff),abs(diff+r360_d),abs(diff-r360_d)) !!! diff=diff2*diff/abs(diff) if (abs(diff) <= abs(diff+r360_d) .and. abs(diff) <= abs(diff-r360_d) ) then diff2=diff else if (abs(diff+r360_d) <= abs(diff) .and. abs(diff+r360_d) <= abs(diff-r360_d) ) then diff2=diff+r360_d else if (abs(diff-r360_d) <= abs(diff) .and. abs(diff-r360_d) <= abs(diff+r360_d) ) then diff2=diff-r360_d endif diff=diff2 endif endif lstats1=abs(rmuse-rmusecv)<=epsilon lstats2=abs(rmuse-rmuseb)<=epsilon lstats3=lstats1.or.lstats2 zeroone(:)=0._8 if (lstats1) then ; nobs(nn,1)=nobs(nn,1)+1 ; zeroone(1)=1._8 ; endif if (lstats2) then ; nobs(nn,2)=nobs(nn,2)+1 ; zeroone(2)=1._8 ; endif if (lstats3) then ; nobs(nn,3)=nobs(nn,3)+1 ; zeroone(3)=1._8 ; endif do k=1,3 if (zeroone(k) > 0.1_8) then s1(nn,k)=s1(nn,k)+abs(diff) s2(nn,k)=s2(nn,k)+diff s3(nn,k)=s3(nn,k)+diff*diff endif enddo ! IF (nn==20 .and. lstats1) then !! PRINT*,'TCAMT,cmodel,rmuse-rmusecv,fint,ob,diff=',trim(cmodel),rmuse,rmusecv,fint,ob,diff ! PRINT*,'TCAMT,=',trim(cmodel),xx,yy,fint,ob,diff ! ENDIF if (.not.lanlerr) then if (nn==7 .or. nn==25) then if (nn==7) then iu0=4 ; iv0=5 endif if (nn==25) then iu0=22 ; iv0=23 endif call bilinear_2d0v2(field0(m1:m2,n1:n2,iu0),m1,m2,n1,n2,uob_model,yy,xx) call bilinear_2d0v2(field0(m1:m2,n1:n2,iv0),m1,m2,n1,n2,vob_model,yy,xx) !!!iu=ipointer0(kk,iu0) ; uob=hobs0(iu) !NO!!! iv=ipointer0(kk,iv0) ; vob=hobs0(iv) uob=rtvts0(iv) !Ugly but correct!!! See how rtvts0 is defined for vwnd10m-ob in cressanl_common.f90 du2=dble(uob_model-uob)*dble(uob_model-uob) !(delta u squared) dv2=dble(vob_model-vob)*dble(vob_model-vob) !(delta v squared) do k=1,3 if (zeroone(k) > 0.1_8) s4(nn,k)=s4(nn,k)+zeroone(k)*(du2+dv2) enddo endif endif endif !mype condition 100 continue call mpi_barrier(mpi_comm_world,ierror) 200 continue deallocate (field0) deallocate (fis0) deallocate (tv0) deallocate (xlocs0) deallocate (ylocs0) deallocate (hobs0) deallocate(rmuses0) deallocate(rtvts0) deallocate(hgt0s0) deallocate(insubdoms0) deallocate (ipointer0) 300 continue npts=nflds*3 nobs99=nobs ; call mpi_allreduce(nobs99, nobs, npts, mpi_integer, mpi_sum, mpi_comm_world, ierror) s99=s1 ; call mpi_allreduce(s99, s1, npts, mpi_real8, mpi_sum, mpi_comm_world, ierror) s99=s2 ; call mpi_allreduce(s99, s2, npts, mpi_real8, mpi_sum, mpi_comm_world, ierror) s99=s3 ; call mpi_allreduce(s99, s3, npts, mpi_real8, mpi_sum, mpi_comm_world, ierror) s99=s4 ; call mpi_allreduce(s99, s4, npts, mpi_real8, mpi_sum, mpi_comm_world, ierror) if (.not.lanlerr) s3(7,:)=s4(7,:) if (.not.lanlerr) s3(25,:)=s4(25,:) bias(:,:)=9999.000 rmse(:,:)=9999.000 abse(:,:)=9999.000 do nn=1,nflds do k=1,3 n=nobs(nn,k) if (n > 0) then abse(nn,k)=sngl(s1(nn,k)/dble(n)) bias(nn,k)=sngl(s2(nn,k)/dble(n)) rmse(nn,k)=sngl(sqrt(s3(nn,k)/dble(n))) endif enddo enddo if (lanlerr) then !force this for now abse(8,:)=9999.000 !30. bias(8,:)=9999.000 !10. rmse(8,:)=9999.000 !30. abse(24:26,:)=9999.000 !30. bias(24:26,:)=9999.000 !10. rmse(24:26,:)=9999.000 !30. endif allcmodel(imodel)=cmodel allnobs(:,:,imodel)=nobs(:,:) allabse(:,:,imodel)=abse(:,:) allbias(:,:,imodel)=bias(:,:) allrmse(:,:,imodel)=rmse(:,:) if (imodel==itotmodel) then do 500 k=1,3 !accounts for rmuse=rmusecv.,1., and (rmusecv. & 1.). if (k==1) then if (rmusecv==-2.) fname='for_rmuse_neg2' if (rmusecv==-3.) fname='for_rmuse_neg3' endif if (k==2) fname='for_rmuse_pos1' if (k==3) fname='for_rmuse_neg3pos1' fname=trim(fname)//trim(fnamesuffix) if (mype==0) then open (17,file='stats_'//trim(fname),form='formatted') open (18,file='rawstats.dat_'//trim(fname),form='unformatted') print*,'in bias_rmse_cv: k,fname=',k,trim(fname) do i=1,itotmodel cmodel0=allcmodel(i) nobsmax0=maxval(allnobs(:,k,i)) !no need really. but leave to conform with content of older files nobs0(:)=allnobs(:,k,i) bias0(:)=allbias(:,k,i) abse0(:)=allabse(:,k,i) rmse0(:)=allrmse(:,k,i) print*,' ' print*,'verification statistics for ',trim(cmodel0) write(17,*) trim(cmodel0) do n=1,nflds print 123,'n,nobs,abse,bias,rmse=',n,nobs0(n),abse0(n),bias0(n),rmse0(n) write(17,123) 'n,nobs,abse,bias,rmse=',n,nobs0(n),abse0(n),bias0(n),rmse0(n) enddo write(18) cmodel0 write(18) nflds,nobsmax0 write(18) nobs0,bias0,abse0,rmse0 enddo close(17) close(18) endif call mpi_barrier(mpi_comm_world,ierror) 500 continue endif 123 format(a23,i3,i7,3(f14.6)) deallocate(i1_info) deallocate(i2_info) deallocate(j1_info) deallocate(j2_info) deallocate(iaux1) deallocate(iaux2) deallocate(iaux3) deallocate(iaux4) deallocate(glon) deallocate(psfc1) deallocate(fis1) deallocate(t1) deallocate(tv1) deallocate(q1) deallocate(u1) ; deallocate(ue) deallocate(v1) ; deallocate(ve) deallocate(w1) deallocate(wd1) deallocate(td1) deallocate(gust1) deallocate(vis1) deallocate(pblh1) deallocate(dist1) deallocate(wspd10m1) deallocate(td2m1) deallocate(mxtm1) deallocate(mitm1) deallocate(pmsl1) deallocate(howv1) deallocate(tcamt1) deallocate(lcbas1) deallocate(uwnd10m1) ; deallocate(uwnd10me) deallocate(vwnd10m1) ; deallocate(vwnd10me) deallocate (altw1) deallocate (altwd1) deallocate(iauxfield) deallocate(auxfield) deallocate(auxfield2) end subroutine bias_rmse_cv_mpi !------------------------------------------------------ subroutine adjust_pressure(zsges,tges0,hgt0,tob,spval,fint) implicit none !Declare passed variables real(4),intent(in):: zsges real(4),intent(in):: tges0 real(4),intent(in):: hgt0 real(4),intent(in):: tob real(4),intent(in):: spval real(4),intent(inout):: fint !Declare local variables real(4) rdp,rdelz,pges,tges real(4) half,half_tlapse,g_over_rd half=0.5 half_tlapse=0.00325 g_over_rd=9.81/287.04 pges=fint rdp=0. rdelz=hgt0-zsges tges=tges0 if (abs(tob-spval).gt.1.e-03) then tges = half*(tges0+tob) else if(rdelz < 0.)then tges=tges0-half_tlapse*rdelz endif endif rdp = g_over_rd*rdelz/tges !Adjust hydrostatically fint=1000.*exp(log(pges/1000.) - rdp) !/100. !hPa end subroutine adjust_pressure !------------------------------------------------------ !------------------------------------------------------