subroutine create_rjlist(cgrid,mype,npe,mpi_comm_new) !******************************************************************** ! abstract: use output from gsi diagnostic files to generate * ! reject lists of observations * ! * ! program history log: * ! 2007-03-13 pondeca * ! * !******************************************************************** use mpi implicit none character(60),intent(in):: cgrid integer(4),intent(in):: mype,npe,mpi_comm_new integer(4),parameter::nobsmax=250000 integer(4),parameter::nflds=5 integer(4),parameter::lun=10 real(4),parameter::rmuse_c=1. !-5000. integer(4) nx,ny integer(4) ierror,nobs integer(4) j,m,n,it,nt,ifld,itype0,ll integer(4) ntw,ntt,ntp,ntq,ntspd,ntmax,ntfld(nflds) integer(4) itype(nobsmax) integer(4) iflag(nobsmax),jflag(nobsmax,12) integer(4) nrj(nflds),nrjsum,iaux(nflds) integer(4),allocatable,dimension(:)::ikx,kflag integer(4) ithwindqc real(4) bmax(nflds,0:19) !maximum bias values real(4) wmax280, tmax180, pmax180, qmax180, spdmax280, & wmax281, tmax181, pmax181, qmax181, spdmax281, & wmax282, tmax182, pmax182, qmax182, spdmax282, & wmax283, tmax183, pmax183, qmax183, spdmax283, & wmax284, tmax184, pmax184, qmax184, spdmax284, & wmax285, tmax185, pmax185, qmax185, spdmax285, & wmax286, tmax186, pmax186, qmax186, spdmax286, & wmax287, tmax187, pmax187, qmax187, spdmax287, & wmax288, tmax188, pmax188, qmax188, spdmax288, & wmax289, tmax189, pmax189, qmax189, spdmax289, & wmax290, tmax190, pmax190, qmax190, spdmax290, & wmax291, tmax191, pmax191, qmax191, spdmax291, & wmax292, tmax192, pmax192, qmax192, spdmax292, & wmax293, tmax193, pmax193, qmax193, spdmax293, & wmax294, tmax194, pmax194, qmax194, spdmax294, & wmax295, tmax195, pmax195, qmax195, spdmax295, & wmax296, tmax196, pmax196, qmax196, spdmax296, & wmax297, tmax197, pmax197, qmax197, spdmax297, & wmax298, tmax198, pmax198, qmax198, spdmax298, & wmax299, tmax199, pmax199, qmax199, spdmax299 real(4) epsw,epst,epsp,epsq,epsspd real(4) fepsw,fepst,fepsp,fepsq,fepsspd real(4) wcutoff,rcutoff real(4) whigh,wlow real(4) ds real(4) valleyfact,valleyfact0 logical valleygcheck real(4) rlat0,rlon0,xx0,yy0,oberr0,ob0,ob_model0,dtime0, & rmuse0,a0,b0 real(4) bias1,f,g,wo,wg real(4) rlat(nobsmax),rlon(nobsmax),& xx(nobsmax),yy(nobsmax),& oberr(nobsmax),ob(nobsmax),ob_model(nobsmax), & dtime(nobsmax),rmuse(nobsmax), & a(nobsmax),b(nobsmax) real(4) eps(nflds),feps(nflds) real(4),allocatable,dimension(:)::alat,alon real(4),allocatable,dimension(:,:)::valleys character*60 statsfile(nflds),rjfile(nflds) character(60) cvar(nflds) character*10 cycall(500),cyc(500),tcyc character(8) cstation(nobsmax),cstation0 character*80 char1,char3 character*80 char2(nflds) character(8),allocatable,dimension(:)::station_id character(2) cloc0 character(7) corigin0 character(1),parameter:: star='*' logical fexist logical lres1,lres2 logical first(nobsmax) logical sfctype logical windstormqc logical lcondition1,lcondition2,lcondition3,mesonetwind namelist/rj_param/tcyc, & wmax280, tmax180, pmax180, qmax180, spdmax280, & wmax281, tmax181, pmax181, qmax181, spdmax281, & wmax282, tmax182, pmax182, qmax182, spdmax282, & wmax283, tmax183, pmax183, qmax183, spdmax283, & wmax284, tmax184, pmax184, qmax184, spdmax284, & wmax285, tmax185, pmax185, qmax185, spdmax285, & wmax286, tmax186, pmax186, qmax186, spdmax286, & wmax287, tmax187, pmax187, qmax187, spdmax287, & wmax288, tmax188, pmax188, qmax188, spdmax288, & wmax289, tmax189, pmax189, qmax189, spdmax289, & wmax290, tmax190, pmax190, qmax190, spdmax290, & wmax291, tmax191, pmax191, qmax191, spdmax291, & wmax292, tmax192, pmax192, qmax192, spdmax292, & wmax293, tmax193, pmax193, qmax193, spdmax293, & wmax294, tmax194, pmax194, qmax194, spdmax294, & wmax295, tmax195, pmax195, qmax195, spdmax295, & wmax296, tmax196, pmax196, qmax196, spdmax296, & wmax297, tmax197, pmax197, qmax197, spdmax297, & wmax298, tmax198, pmax198, qmax198, spdmax298, & wmax299, tmax199, pmax199, qmax199, spdmax299, & ntw,ntt,ntp,ntq,ntspd,cycall, & epsw,epst,epsp,epsq,epsspd, & fepsw,fepst,fepsp,fepsq,fepsspd, & windstormqc,wcutoff,rcutoff, & ithwindqc,whigh,wlow, & valleygcheck,valleyfact0 data statsfile/'stats_w_rj.out', & 'stats_t_rj.out', & 'stats_ps_rj.out',& 'stats_q_rj.out',& 'stats_spd_rj.out'/ data rjfile/'w_rjlist.txt_dynamic', & 't_rjlist.txt_dynamic', & 'p_rjlist.txt_dynamic',& 'q_rjlist.txt_dynamic',& 'spd_rjlist.txt_dynamic'/ data cvar/'wind','temperature', & 'surface pressure',& 'moisture','wind speed'/ data epsw/0.01/ ; data fepsw/10./ data epst/0.01/ ; data fepst/10./ data epsp/10./ ; data fepsp/10./ !note: input pressure data to this routine are in Pa data epsq/1.e-06/ ; data fepsq/5./ data epsspd/0.01/ ; data fepsspd/10./ data windstormqc/.false./ data wcutoff/7.732/ !this is 15./1.94 ,that is 15 knots converted to SI units data rcutoff/0.75/ data ithwindqc/0/ !if <0 do not apply tyndall/horel wind qc. If 0 apply to mesonets only, if > 0 apply to all ob types data whigh/10./ !values for tyndall/horel data wlow/1./ !wind QC data valleygcheck/.false./ data valleyfact0/3./ ! data char1/'***********************************& ! &*******************'/ data char2/'reject list for wind obs', & 'reject list for temperature obs', & 'reject list for surface pressure obs', & 'reject list for specific humidity obs', & 'reject list for satellite wind speed obs'/ wmax280=5.0; tmax180=10.; pmax180=3000.; qmax180=5.e-03; spdmax280=15. wmax281=5.0; tmax181=10.; pmax181=3000.; qmax181=5.e-03; spdmax281=15. wmax282=5.0; tmax182=10.; pmax182=3000.; qmax182=5.e-03; spdmax282=15. wmax283=5.0; tmax183=10.; pmax183=3000.; qmax183=5.e-03; spdmax283=15. wmax284=5.0; tmax184=10.; pmax184=3000.; qmax184=5.e-03; spdmax284=15. wmax285=5.0; tmax185=10.; pmax185=3000.; qmax185=5.e-03; spdmax285=15. wmax286=5.0; tmax186=10.; pmax186=3000.; qmax186=5.e-03; spdmax286=15. wmax287=100.; tmax187=10.; pmax187=3000.; qmax187=5.e-03; spdmax287=15. wmax288=3.0; tmax188=10.; pmax188=3000.; qmax188=5.e-03; spdmax288=15. wmax289=5.0; tmax189=10.; pmax189=3000.; qmax189=5.e-03; spdmax289=15. wmax290=5.0; tmax190=10.; pmax190=3000.; qmax190=5.e-03; spdmax290=15. wmax291=5.0; tmax191=10.; pmax191=3000.; qmax191=5.e-03; spdmax291=15. wmax292=5.0; tmax192=10.; pmax192=3000.; qmax192=5.e-03; spdmax292=15. wmax293=5.0; tmax193=10.; pmax193=3000.; qmax193=5.e-03; spdmax293=15. wmax294=5.0; tmax194=10.; pmax194=3000.; qmax194=5.e-03; spdmax294=15. wmax295=5.0; tmax195=10.; pmax195=3000.; qmax195=5.e-03; spdmax295=15. wmax296=5.0; tmax196=10.; pmax196=3000.; qmax196=5.e-03; spdmax296=15. wmax297=5.0; tmax197=10.; pmax197=3000.; qmax197=5.e-03; spdmax297=15. wmax298=5.0; tmax198=10.; pmax198=3000.; qmax198=5.e-03; spdmax298=15. wmax299=5.0; tmax199=10.; pmax199=3000.; qmax199=5.e-03; spdmax299=15. ntw=1 ntt=3 ntp=3 ntq=3 ntspd=3 do n=1,80 char1(n:n)=star enddo cloc0(1:2)='US' corigin0='dynamic' ntfld(:)=1 cycall(:)=' ' bmax(:,:)=1.e+20 print*,'in create_rjlist: npe,mype=',npe,mype if (mype.eq.0) then print*,'in create_rjlist:: cgrid=',trim(cgrid) call streamline_for_rjlist(cgrid) endif call mpi_barrier(mpi_comm_new,ierror) open (55,file='rj_param_input',form='formatted') read(55,rj_param) close(55) cyc(:)=cycall(:) !Note: use of intermediary field cycall(:) is last minute fix ! to avoid variable naming conflict with NCO RTMA post script/30May2007 if (mype == 0 ) then write(6,100) 100 format(' calling create_rjlist with following input parameters:',//) write(6,*) 'tcyc=',tcyc ntmax=max(ntw,ntt,ntp,ntq,ntspd) do n=1,ntmax write(6,*) 'cyc(',n,')=',cyc(n) enddo write(6,1001) wmax280, tmax180, pmax180, qmax180, spdmax280 write(6,1002) wmax281, tmax181, pmax181, qmax181, spdmax281 write(6,1003) wmax282, tmax182, pmax182, qmax182, spdmax282 write(6,1004) wmax283, tmax183, pmax183, qmax183, spdmax283 write(6,1005) wmax284, tmax184, pmax184, qmax184, spdmax284 write(6,1006) wmax285, tmax185, pmax185, qmax185, spdmax285 write(6,1007) wmax286, tmax186, pmax186, qmax186, spdmax286 write(6,1008) wmax287, tmax187, pmax187, qmax187, spdmax287 write(6,1009) wmax288, tmax188, pmax188, qmax188, spdmax288 write(6,1010) wmax289, tmax189, pmax189, qmax189, spdmax289 write(6,1011) wmax290, tmax190, pmax190, qmax190, spdmax290 write(6,1012) wmax291, tmax191, pmax191, qmax191, spdmax291 write(6,1013) wmax292, tmax192, pmax192, qmax192, spdmax292 write(6,1014) wmax293, tmax193, pmax193, qmax193, spdmax293 write(6,1015) wmax294, tmax194, pmax194, qmax194, spdmax294 write(6,1016) wmax295, tmax195, pmax195, qmax195, spdmax295 write(6,1016) wmax296, tmax196, pmax196, qmax196, spdmax296 write(6,1016) wmax297, tmax197, pmax197, qmax197, spdmax297 write(6,1016) wmax298, tmax198, pmax198, qmax198, spdmax298 write(6,1016) wmax299, tmax199, pmax199, qmax199, spdmax299 write(6,*) 'ntw,ntt,ntp,ntq,ntspd=',ntw,ntt,ntp,ntq,ntspd write(6,*) 'epsw,epst,epsp,epsq,epsspd=',epsw,epst,epsp,epsq,epsspd write(6,*) 'fepsw,fepst,fepsp,fepsq,fepsspd=',fepsw,fepst,fepsp,fepsq,fepsspd write(6,*) 'windstormqc,wcutoff,rcutoff=',windstormqc,wcutoff,rcutoff write(6,*) 'ithwindqc,whigh,wlow=',ithwindqc,whigh,wlow write(6,*) 'valleygcheck=',valleygcheck endif bmax(1,0)=wmax280; bmax(2,0)=tmax180; bmax(3,0)=pmax180 bmax(1,1)=wmax281; bmax(2,1)=tmax181; bmax(3,1)=pmax181 bmax(1,2)=wmax282; bmax(2,2)=tmax182; bmax(3,2)=pmax182 bmax(1,3)=wmax283; bmax(2,3)=tmax183; bmax(3,3)=pmax183 bmax(1,4)=wmax284; bmax(2,4)=tmax184; bmax(3,4)=pmax184 bmax(1,5)=wmax285; bmax(2,5)=tmax185; bmax(3,5)=pmax185 bmax(1,6)=wmax286; bmax(2,6)=tmax186; bmax(3,6)=pmax186 bmax(1,7)=wmax287; bmax(2,7)=tmax187; bmax(3,7)=pmax187 bmax(1,8)=wmax288; bmax(2,8)=tmax188; bmax(3,8)=pmax188 bmax(1,9)=wmax289; bmax(2,9)=tmax189; bmax(3,9)=pmax189 bmax(1,10)=wmax290; bmax(2,10)=tmax190; bmax(3,10)=pmax190 bmax(1,11)=wmax291; bmax(2,11)=tmax191; bmax(3,11)=pmax191 bmax(1,12)=wmax292; bmax(2,12)=tmax192; bmax(3,12)=pmax192 bmax(1,13)=wmax293; bmax(2,13)=tmax193; bmax(3,13)=pmax193 bmax(1,14)=wmax294; bmax(2,14)=tmax194; bmax(3,14)=pmax194 bmax(1,15)=wmax295; bmax(2,15)=tmax195; bmax(3,15)=pmax195 bmax(1,16)=wmax296; bmax(2,16)=tmax196; bmax(3,16)=pmax196 bmax(1,17)=wmax297; bmax(2,17)=tmax197; bmax(3,17)=pmax197 bmax(1,18)=wmax298; bmax(2,18)=tmax198; bmax(3,18)=pmax198 bmax(1,19)=wmax299; bmax(2,19)=tmax199; bmax(3,19)=pmax199 bmax(4,0)=qmax180; bmax(5,0)=spdmax280 bmax(4,1)=qmax181; bmax(5,1)=spdmax281 bmax(4,2)=qmax182; bmax(5,2)=spdmax282 bmax(4,3)=qmax183; bmax(5,3)=spdmax283 bmax(4,4)=qmax184; bmax(5,4)=spdmax284 bmax(4,5)=qmax185; bmax(5,5)=spdmax285 bmax(4,6)=qmax186; bmax(5,6)=spdmax286 bmax(4,7)=qmax187; bmax(5,7)=spdmax287 bmax(4,8)=qmax188; bmax(5,8)=spdmax288 bmax(4,9)=qmax189; bmax(5,9)=spdmax289 bmax(4,10)=qmax190; bmax(5,10)=spdmax290 bmax(4,11)=qmax191; bmax(5,11)=spdmax291 bmax(4,12)=qmax192; bmax(5,12)=spdmax292 bmax(4,13)=qmax193; bmax(5,13)=spdmax293 bmax(4,14)=qmax194; bmax(5,14)=spdmax294 bmax(4,15)=qmax195; bmax(5,15)=spdmax295 bmax(4,16)=qmax196; bmax(5,16)=spdmax296 bmax(4,17)=qmax197; bmax(5,17)=spdmax297 bmax(4,18)=qmax198; bmax(5,18)=spdmax298 bmax(4,19)=qmax199; bmax(5,19)=spdmax299 ntfld(1)=ntw ; ntfld(2)=ntt ; ntfld(3)=ntp ntfld(4)=ntq ; ntfld(5)=ntspd eps(1)=epsw ; feps(1)=fepsw eps(2)=epst ; feps(2)=fepst eps(3)=epsp ; feps(3)=fepsp eps(4)=epsq ; feps(4)=fepsq eps(5)=epsspd ; feps(5)=fepsspd if (valleygcheck) then call domain_dims(cgrid,nx,ny,ds) if (mype==0) print*,'in create_rjlist: cgrid,nx,ny=',cgrid,nx,ny allocate(valleys(nx,ny)) valleys(:,:)=1. inquire(file='valley_map.dat',exist=fexist) if (fexist) then open (55,file='valley_map.dat',form='unformatted') read(55) valleys close(55) endif print*,'in create_rjlist: valleys,min,max=',minval(valleys),maxval(valleys) endif nrj(:)=0 do 5000 ifld=1,nflds if (mod(ifld-1,npe) == mype ) then print*,'in create_rjlist: mype,ifld,statsfile=',mype,ifld,trim(statsfile(ifld)) inquire(file=trim(statsfile(ifld)),exist=fexist) if (.not.fexist) cycle open (lun,file=trim(statsfile(ifld)),form='formatted') nobs=0 do n=1,nobsmax read(lun,'(a8)',end=200) cstation0 read(lun,*,end=200) itype0,rlat0,rlon0,xx0,yy0,dtime0, & oberr0,ob0,ob_model0,a0,b0,rmuse0 sfctype=(itype0>=180.and.itype0<190).or.(itype0>=280.and.itype0<=290).or. & (itype0>=192.and.itype0<=199).or.(itype0>=292.and.itype0<=299) if (sfctype) then nobs=nobs+1 cstation(nobs)(1:8)=cstation0(1:8) itype(nobs)=itype0 rlat(nobs)=rlat0 rlon(nobs)=rlon0 xx(nobs)=xx0 yy(nobs)=yy0 dtime(nobs)=dtime0 oberr(nobs)=oberr0 ob(nobs)=ob0 ob_model(nobs)=ob_model0 rmuse(nobs)=rmuse0 a(nobs)=a0 b(nobs)=b0 endif enddo 200 continue close(lun) print*,'in create_rjlist: mype,cvar(ifld),nobs=', & mype,trim(cvar(ifld)),nobs iflag(:)=-1 jflag(:,:)=+1 do n=1,nobs bias1=ob_model(n)-ob(n) if (ifld == 1 ) bias1=sqrt( (ob_model(n)-ob(n))*(ob_model(n)-ob(n)) + & (b(n)-a(n))*(b(n)-a(n)) ) if (itype(n) < 200) ll=itype(n)-180 if (itype(n) >= 200) ll=itype(n)-280 !bias-based QC valleyfact=1. ! if ((ifld==2 .or. ifld==4) .and. valleygcheck) & ! call valley_probe(valleys,nx,ny,xx(n),yy(n),valleyfact0,valleyfact) if (valleygcheck) call valley_probe(valleys,nx,ny,xx(n),yy(n),valleyfact0,valleyfact) if (ifld.ne.2 .and. ifld.ne.4) valleyfact=1. !apply valleygchect to T and Q only if ( abs(bias1) .le. valleyfact*bmax(ifld,ll)) iflag(n)=+1 !added QC for winds wo=sqrt(ob(n)*ob(n)+a(n)*a(n)) wg=sqrt(ob_model(n)*ob_model(n)+b(n)*b(n)) mesonetwind=itype(n)==288.or.itype(n)==295 !mesonet winds can be unrealistically light during high impact wind events if ( mesonetwind .and. ( windstormqc .and. (wg>wcutoff .and. (wo/wg)0.and.wg>whigh.and.wowhigh.and.wo=180.and.itype0<190).or.(itype0>=280.and.itype0<=290).or. & (itype0>=192.and.itype0<=199).or.(itype0>=292.and.itype0<=299) if (.not.sfctype) cycle bias1=ob_model0-ob0 if (ifld == 1 ) bias1=sqrt( (ob_model0-ob0)*(ob_model0-ob0) + & (b0-a0)*(b0-a0) ) wo=sqrt(ob0*ob0+a0*a0) wg=sqrt(ob_model0*ob_model0+b0*b0) mesonetwind=itype0==288.or.itype0==295 lcondition1=ithwindqc>0.and.wg>whigh.and.wowhigh.and.wo= 200) ll=itype0-280 do n=1,nobs if (cstation(n)(1:8)==cstation0(1:8) .and. itype(n)==itype0) then if (itype0>=200) then lcondition3=(abs(bias1)>bmax(ifld,ll)).or. & (windstormqc.and.(wg>wcutoff .and. (wo/wg) 0 .and. iflag(m) < 0) iflag(n)=-1 iflag(m)=2 endif enddo enddo write(lun,'(a80)') char1 write(lun,'(a80)') char2(ifld) write(lun,'(a80)') char1 do n=1,nobs if (iflag(n) .eq. -1 ) then nrj(ifld)=nrj(ifld)+1 if (rlon(n) > 180.) rlon(n)=rlon(n)-360. write(lun,1234) cstation(n),itype(n),rlat(n),rlon(n),cloc0,corigin0 endif enddo print*,'in create_rjlist: ifld,nrj=',ifld,nrj(ifld) close(lun) endif 5000 continue !1234 format("'",a8,'| itype=',i3,3x,'lat=',f10.4,3x,'lon=',f10.4,"'") 1234 format("'",a8,'| itype=',i3,2x,'lat=',f10.4,2x,'lon=',f10.4,2x, & 'loc=',a2,2x,'origin:',a7,"'") if (valleygcheck) deallocate(valleys) !==>now processors '0' and '1' generate mass and-wind reject lists. iaux(:)=0 call mpi_allreduce(nrj,iaux,nflds,mpi_integer,mpi_sum, & mpi_comm_new,ierror) nrj(:)=iaux(:) if (mype==mod(0,npe)) then !see mark-0 open (20,file='mass_rjlist.txt_dynamic'//'_'//tcyc, & form='formatted') char3='reject list for all mass obs' write(20,'(a80)') char1 write(20,'(a80)') char3 write(20,'(a80)') char1 nrjsum=nrj(2)+nrj(4) !+nrj(3) !put back on to include sfcp if (nrjsum .gt. 0) then !see mark-1 allocate(station_id(nrjsum)) allocate(ikx(nrjsum)) allocate(alat(nrjsum)) allocate(alon(nrjsum)) allocate(kflag(nrjsum)) if (nrj(2).gt.0) then open (10,file=trim(rjfile(2))//'_'//tcyc, & form='formatted') read(10,*) read(10,*) read(10,*) do n=1,nrj(2) read(10,1234) cstation0,itype0,rlat0,rlon0,cloc0,corigin0 station_id(n)=cstation0 ikx(n)=itype0 alat(n)=rlat0 alon(n)=rlon0 enddo close(10) endif if (nrj(4).gt.0) then open (10,file=trim(rjfile(4))//'_'//tcyc, & form='formatted') read(10,*) read(10,*) read(10,*) do n=nrj(2)+1,nrj(2)+nrj(4) read(10,1234) cstation0,itype0,rlat0,rlon0,cloc0,corigin0 station_id(n)=cstation0 ikx(n)=itype0 alat(n)=rlat0 alon(n)=rlon0 enddo close(10) endif ! if (nrj(3).gt.0) then !put back on to include sfcp ! open (10,file=trim(rjfile(3))//'_'//tcyc, & ! form='formatted') ! read(10,*) ! read(10,*) ! read(10,*) ! do n=nrj(2)+nrj(4)+1,nrjsum ! read(10,1234) cstation0,itype0,rlat0,rlon0,cloc0,corigin0 ! station_id(n)=cstation0 ! ikx(n)=itype0 ! alat(n)=rlat0 ! alon(n)=rlon0 ! enddo ! close(10) ! endif kflag(:)=-1 do n=1,nrjsum do m=n+1,nrjsum if (station_id(m)==station_id(n)) kflag(m)=+1 enddo enddo j=0 do n=1,nrjsum if (kflag(n) .eq. -1 ) then j=j+1 write(20,1234) station_id(n),ikx(n),alat(n),alon(n),cloc0,corigin0 endif enddo print*,'in create_rjlist: mass,nrj=',j close(20) deallocate(station_id) deallocate(ikx) deallocate(alat) deallocate(alon) deallocate(kflag) endif !mark-1 endif !mark-0 if (mype==mod(1,npe)) then !see mark-00 open (20,file='wind_rjlist.txt_dynamic'//'_'//tcyc, & form='formatted') char3='reject list for all wind obs' write(20,'(a80)') char1 write(20,'(a80)') char3 write(20,'(a80)') char1 nrjsum=nrj(1)+nrj(5) if (nrjsum .gt. 0) then !see mark-2 allocate(station_id(nrjsum)) allocate(ikx(nrjsum)) allocate(alat(nrjsum)) allocate(alon(nrjsum)) allocate(kflag(nrjsum)) if (nrj(1) .gt.0) then open (10,file=trim(rjfile(1))//'_'//tcyc, & form='formatted') read(10,*) read(10,*) read(10,*) do n=1,nrj(1) read(10,1234) cstation0,itype0,rlat0,rlon0,cloc0,corigin0 station_id(n)=cstation0 ikx(n)=itype0 alat(n)=rlat0 alon(n)=rlon0 enddo close(10) endif if (nrj(5) .gt.0) then open (10,file=trim(rjfile(5))//'_'//tcyc, & form='formatted') read(10,*) read(10,*) read(10,*) do n=nrj(1)+1,nrjsum read(10,1234) cstation0,itype0,rlat0,rlon0,cloc0,corigin0 station_id(n)=cstation0 ikx(n)=itype0 alat(n)=rlat0 alon(n)=rlon0 enddo close(10) endif kflag(:)=-1 do n=1,nrjsum do m=n+1,nrjsum if (station_id(m)==station_id(n)) kflag(m)=+1 enddo enddo j=0 do n=1,nrjsum if (kflag(n) .eq. -1 ) then j=j+1 write(20,1234) station_id(n),ikx(n),alat(n),alon(n),cloc0,corigin0 endif enddo print*,'in create_rjlist: wind,nrj=',j close(20) deallocate(station_id) deallocate(ikx) deallocate(alat) deallocate(alon) deallocate(kflag) endif !mark-2 endif !mark-00 1001 format('wmax280=',f5.1,1x,'tmax180=',f5.1,1x,'pmax180=',f7.0,1x,'qmax180=',e9.2,1x,'spdmax280=',f5.1) 1002 format('wmax281=',f5.1,1x,'tmax181=',f5.1,1x,'pmax181=',f7.0,1x,'qmax181=',e9.2,1x,'spdmax281=',f5.1) 1003 format('wmax282=',f5.1,1x,'tmax182=',f5.1,1x,'pmax182=',f7.0,1x,'qmax182=',e9.2,1x,'spdmax282=',f5.1) 1004 format('wmax283=',f5.1,1x,'tmax183=',f5.1,1x,'pmax183=',f7.0,1x,'qmax183=',e9.2,1x,'spdmax283=',f5.1) 1005 format('wmax284=',f5.1,1x,'tmax184=',f5.1,1x,'pmax184=',f7.0,1x,'qmax184=',e9.2,1x,'spdmax284=',f5.1) 1006 format('wmax285=',f5.1,1x,'tmax185=',f5.1,1x,'pmax185=',f7.0,1x,'qmax185=',e9.2,1x,'spdmax285=',f5.1) 1007 format('wmax286=',f5.1,1x,'tmax186=',f5.1,1x,'pmax186=',f7.0,1x,'qmax186=',e9.2,1x,'spdmax286=',f5.1) 1008 format('wmax287=',f5.1,1x,'tmax187=',f5.1,1x,'pmax187=',f7.0,1x,'qmax187=',e9.2,1x,'spdmax287=',f5.1) 1009 format('wmax288=',f5.1,1x,'tmax188=',f5.1,1x,'pmax188=',f7.0,1x,'qmax188=',e9.2,1x,'spdmax288=',f5.1) 1010 format('wmax289=',f5.1,1x,'tmax189=',f5.1,1x,'pmax189=',f7.0,1x,'qmax189=',e9.2,1x,'spdmax289=',f5.1) 1011 format('wmax290=',f5.1,1x,'tmax190=',f5.1,1x,'pmax190=',f7.0,1x,'qmax190=',e9.2,1x,'spdmax290=',f5.1) 1012 format('wmax291=',f5.1,1x,'tmax191=',f5.1,1x,'pmax191=',f7.0,1x,'qmax191=',e9.2,1x,'spdmax291=',f5.1) 1013 format('wmax292=',f5.1,1x,'tmax192=',f5.1,1x,'pmax192=',f7.0,1x,'qmax192=',e9.2,1x,'spdmax292=',f5.1) 1014 format('wmax293=',f5.1,1x,'tmax193=',f5.1,1x,'pmax193=',f7.0,1x,'qmax193=',e9.2,1x,'spdmax293=',f5.1) 1015 format('wmax294=',f5.1,1x,'tmax194=',f5.1,1x,'pmax194=',f7.0,1x,'qmax194=',e9.2,1x,'spdmax294=',f5.1) 1016 format('wmax295=',f5.1,1x,'tmax195=',f5.1,1x,'pmax195=',f7.0,1x,'qmax195=',e9.2,1x,'spdmax295=',f5.1) call mpi_barrier(mpi_comm_new,ierror) end subroutine create_rjlist subroutine valley_probe(valleys,nx,ny,xob,yob,valleyfact0,valleyfact) !================================================================ ! (im,jp) (i,jp) ip,jp) ! -------------------------- ! | | | ! | | | ! | | | ! | | | ! (im,j)----------(i,j)----------- (ip.j) ! | | | ! | | | ! | | | ! | | | ! -------------------------- ! (im,jm) (i,jm) ip,jm) !================================================================ implicit none !Declare passed variables integer(4),intent(in):: nx,ny real(4),intent(in):: valleys(nx,ny) real(4),intent(in):: xob,yob,valleyfact0 real(4),intent(inout):: valleyfact !Declare local variables integer(4) i,j,im,jm,ip,jp i=int(xob) ; im=max(1,i-1) ; ip=min(nx,i+1) j=int(yob) ; jm=max(1,j-1) ; jp=min(ny,j+1) valleyfact=1. if (valleys(i,j )<1. .or. valleys(im,j)<1. .or. valleys(ip,j)<1. .or. & valleys(i,jm)<1. .or. valleys(i,jp)<1. .or. & valleys(im,jm)<1. .or. valleys(ip,jm)<1. .or. & valleys(im,jp)<1. .or. valleys(ip,jp)<1. ) valleyfact=valleyfact0 end subroutine valley_probe