program anl_error !*********************************************************************** ! abstract: use barnes or cressman analysis to obtain gridded field * ! of analysis uncertainty * ! * ! program history log: * ! 2005-10-08 pondeca * ! 2018-02-06 carley/yang - add option to skip calculation of anlerr * ! 2018-05-18 yang/pondeca - modified C & V error conversion in NLTR * ! * !*********************************************************************** use mpi use kinds, only: i_kind use mpitaskmod, only: mype,npe implicit none ! Declare local variables integer(i_kind) ierror logical projway,skipanlerr logical anlqlty logical lobjanl logical climatology logical cvbasedrecalibration logical lanczos_writediag character(60) cvbasedcmodel0 integer(4) miter0 namelist/anlerrmethod/projway,skipanlerr,anlqlty,lobjanl,climatology,cvbasedrecalibration, & cvbasedcmodel0,miter0,lanczos_writediag ! MPI setup call mpi_init(ierror) call mpi_comm_size(mpi_comm_world,npe,ierror) call mpi_comm_rank(mpi_comm_world,mype,ierror) call run_biascor(mype) call mpi_barrier(mpi_comm_world,ierror) projway=.true. skipanlerr=.false. climatology=.false. anlqlty=.false. lobjanl=.false. cvbasedrecalibration=.false. cvbasedcmodel0='siganl' miter0=2 !default number of outer loops in gsi lanczos_writediag=.false. open (55,file='lanczosparm.anl',form='formatted') read(55,anlerrmethod) close(55) projway=.true. !must force this after reading namelist, since !projway=.false. is no longer a valid option /18Oct2012 if (mype==0) then print*,'projway=',projway print*,'skipanlerr=',skipanlerr print*,'anlqlty=',anlqlty print*,'lobjanl=',lobjanl print*,'climatology=',climatology print*,'cvbasedrecalibration=',cvbasedrecalibration print*,'cvbasedcmodel0=',cvbasedcmodel0 print*,'miter0=',miter0 print*,'lanczos_writediag=',lanczos_writediag endif if (projway) call lanczos_anlerr(anlqlty,skipanlerr,lobjanl,climatology,cvbasedrecalibration, & cvbasedcmodel0,miter0,lanczos_writediag) ! if (.not.projway) call simplified_anlerr(mype,npe) call mpi_finalize(ierror) end !================================================================================= !================================================================================= subroutine lanczos_anlerr(anlqlty,skipanlerr,lobjanl,climatology,cvbasedrecalibration, & cvbasedcmodel0,miter0,lanczos_writediag) !******************************************************************** ! abstract: use barnes or cressman analysis to obtain gridded field * ! of analysis uncertainty * ! Must run with a minimum of 5 tasks * ! Change this later so that less tasks also possible * ! * ! program history log: * ! 2006-11-20 pondeca * ! 2018-02-06 carley/yang - add option to skip calculation * ! of anlerr * ! * !******************************************************************** use mpi use kinds, only: r_single,r_kind,r_double,i_kind use constants, only: zero,one,tiny_r_kind,tiny_single,half use controlvars, only: nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw,& nrf3_sfwter,nrf3_vpwter,nrf2_ps,nrf2_sst,nrf2_gust,& nrf2_vis,nrf2_pblh,nrf2_dist,nrf2_wspd10m,nrf2_td2m,& nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas,& nrf2_uwnd10m,nrf2_vwnd10m,nrf1_rad,nrf1_prec use controlvars, only: aka_controlvars use mpitaskmod, only: mype,npe use mpitaskmod, only: nlon,nlat,nsig,ista,iend,jsta,jend use cressanl_common, only: make_cressanl_common, & destroy_cressanl_common use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, & wspderr,wdirerr,wdirerr2,terr,tderr, & qerr,perr,gusterr,viserr,pblherr,disterr, & wspd10merr,td2merr,mxtmerr,mitmerr,pmslerr, & howverr,tcamterr,lcbaserr,uwnd10merr,vwnd10merr use errs_common, only: create_errs_common,destroy_errs_common implicit none ! Declare passed variables logical,intent(in):: anlqlty,skipanlerr logical,intent(in):: lobjanl logical,intent(in):: climatology logical,intent(in):: cvbasedrecalibration logical,intent(in):: lanczos_writediag character(60),intent(in):: cvbasedcmodel0 integer(4),intent(in):: miter0 ! Declare local variables character(60) cgrid integer(i_kind) ierror,ierror2 integer(i_kind) jpch,npred,jtype,npredp,jiter,nv,nvused integer(i_kind) nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1 integer(i_kind) i,j,k,kall,lun1,lun2,m,n,ifield,nsym,kk integer(i_kind) num_slabs,num_pad,islab,islab_prev integer(i_kind) isub,jsub,ijsub integer(i_kind) nx,ny,nk integer(i_kind) npointer(24) !total number of control variables integer(i_kind),allocatable,dimension(:)::iscnt,isdisp,ircnt,irdisp integer(i_kind),allocatable,dimension(:)::ijpts_in_sub integer(i_kind),allocatable,dimension(:)::ista_info,iend_info,jsta_info,jend_info integer(i_kind),allocatable,dimension(:)::iaux0,iaux1,iaux2,iaux3,iaux4 real(r_double) innerprod,dotp0,dotp real(r_kind),allocatable,dimension(:,:)::h_slab,g_slab real(r_kind),allocatable,dimension(:)::h_rad,h_prec,g_rad,g_prec real(r_kind),allocatable,dimension(:)::hfield,gfield real(r_kind),allocatable,dimension(:,:,:,:)::hall_loc,gall_loc real(r_kind),allocatable,dimension(:,:,:,:)::hall_loc2,gall_loc2 real(r_double),allocatable,dimension(:)::cnorm character(2) clun1 character(3) clun2 character(3) clun3 character(20) varname1,varname2 logical twodvar_regional logical fexist1 logical trouble logical,parameter:: leigenvector_write=.false. real(r_kind),allocatable,dimension(:,:,:)::psi real(r_kind),allocatable,dimension(:,:,:)::chi real(r_kind),allocatable,dimension(:,:,:)::t real(r_kind),allocatable,dimension(:,:,:)::q real(r_kind),allocatable,dimension(:,:,:)::oz real(r_kind),allocatable,dimension(:,:,:)::cwmr real(r_kind),allocatable,dimension(:,:,:)::sfwter real(r_kind),allocatable,dimension(:,:,:)::vpwter real(r_kind),allocatable,dimension(:,:)::sfcp real(r_kind),allocatable,dimension(:,:)::sfct real(r_kind),allocatable,dimension(:,:)::gust,vis,pblh,dist real(r_kind),allocatable,dimension(:,:)::wspd10m,td2m,mxtm,mitm real(r_kind),allocatable,dimension(:,:)::pmsl,howv,tcamt,lcbas real(r_kind),allocatable,dimension(:,:)::uwnd10m,vwnd10m real(r_kind),allocatable,dimension(:,:,:,:)::u,u2 real(r_kind),allocatable,dimension(:,:,:,:)::v,v2 real(r_kind),allocatable,dimension(:,:,:,:)::wdir,wdir2 real(r_kind),allocatable,dimension(:,:,:,:)::uwter,u2wter real(r_kind),allocatable,dimension(:,:,:,:)::vwter,v2wter real(r_kind),allocatable,dimension(:,:,:,:)::wdirwter,wdir2wter real(r_kind),allocatable,dimension(:,:,:)::usum,usum2 real(r_kind),allocatable,dimension(:,:,:)::vsum,vsum2 real(r_kind),allocatable,dimension(:,:,:)::wdirsum,wdirsum2 real(r_kind),allocatable,dimension(:,:)::uaux1 real(r_kind),allocatable,dimension(:,:)::vaux1 real(r_kind),allocatable,dimension(:,:)::uaux2 real(r_kind),allocatable,dimension(:,:)::vaux2 real(r_kind),allocatable,dimension(:,:)::coeffx real(r_kind),allocatable,dimension(:,:)::coeffy real(r_kind),allocatable,dimension(:,:)::grad,hrad real(r_kind),allocatable,dimension(:,:)::gprec,hprec real(r_kind),allocatable,dimension(:,:)::grad2,hrad2 real(r_kind),allocatable,dimension(:,:)::gprec2,hprec2 real(r_double),allocatable,dimension(:,:)::hgdotaux,hgdot real(r_kind),allocatable,dimension(:,:)::tmatrix real(r_kind),allocatable,dimension(:,:)::slab1,slab2 real(r_kind),allocatable,dimension(:)::alpha,gamma real(r_kind),allocatable,dimension(:)::beta,delta,theta real(r_kind),allocatable,dimension(:)::atest,eval real(r_kind),allocatable,dimension(:,:)::rtest real(r_kind),allocatable,dimension(:,:)::tempb8,tempc8 real(r_single),allocatable,dimension(:,:)::tempa,tempb real(r_single),allocatable,dimension(:,:)::glon real(r_single),allocatable,dimension(:,:)::rlons,rlats real(r_single),allocatable,dimension(:,:)::rmask real(r_single),allocatable,dimension(:,:):: climerr real(r_single) ds0 real(r_kind),parameter::ubar=10. real(r_kind),parameter::vbar=10. real(r_kind),parameter::r25=25._r_kind real(r_kind),parameter::r210=210._r_kind real(r_kind),parameter::r265=265._r_kind real(r_kind),parameter::r38_5=38.5_r_kind real(r_kind),parameter::r262_5=262.5_r_kind real(r_kind),parameter::r225=225._r_kind real(r_kind) dg2rad,xn,elonv,angle2,sinx2,cosx2 real(r_kind) alpha0,gamma0 integer(i_kind),parameter::nvmin=4 real(r_single) pbiascor,tbiascor,qbiascor,ubiascor,vbiascor,tdbiascor, & gustbiascor,visbiascor,pblhbiascor,tcamtbiascor,lcbasbiascor, & wspdbiascor,mxtmbiascor,mitmbiascor,pmslbiascor,howvbiascor, & uwndbiascor,vwndbiascor real(r_single) sfcr0 logical lbiascor logical mkrjlists logical hwrfblend logical use_sfcr0 logical neutral_stability_windfact_2dvar logical use_similarity_2dvar integer(i_kind) icolor, ikey, mpi_comm_new, mypenew, npenew integer(i_kind) mypecutoff real(r_single),parameter::ribuffer_km=250. real(r_single),parameter::rjbuffer_km=250. character(60) tdfname character(60) filename integer(i_kind) ista2,iend2,jsta2,jend2 !memory dimensions real(8) amin8,amax8,amin82,amax82 namelist/gridname/cgrid,lbiascor,pbiascor,tbiascor,qbiascor,ubiascor,vbiascor, & gustbiascor,visbiascor,pblhbiascor, & tdbiascor,wspdbiascor,mxtmbiascor,mitmbiascor,pmslbiascor, & howvbiascor,tcamtbiascor,lcbasbiascor,uwndbiascor,vwndbiascor, & mkrjlists,hwrfblend,use_sfcr0,sfcr0, & neutral_stability_windfact_2dvar,use_similarity_2dvar !************************************************************************* !************************************************************************* data cgrid/'conus'/ open (55,file='gridname_input',form='formatted') read(55,gridname) close(55) if (mype==0) print*,'in lanczos_anlerr:, cgrid=',cgrid if (trim(cgrid) /= 'conus' .and. & trim(cgrid) /= 'alaska' .and. & trim(cgrid) /= 'hawaii' .and. & trim(cgrid) /= 'guam' .and. & trim(cgrid) /= 'prico' .and. & trim(cgrid) /= 'cohres' .and. & trim(cgrid) /= 'akhres' .and. & trim(cgrid) /= 'hrrr' .and. & trim(cgrid) /= 'juneau' .and. & trim(cgrid) /= 'cohresext' .and. & trim(cgrid) /= 'cohreswexp') then print*,'in anl_error: unknown grid ',trim(cgrid),'...aborting' call mpi_abort(mpi_comm_world,ierror2,ierror) endif call domain_dims(cgrid,nx,ny,ds0) if (mype==0) print*,'in lanczos_anlerr:, ds0=',ds0 lun1=19 lun2=29 open (lun1,file='grady.dat',form='unformatted') read (lun1) nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv,alpha0,gamma0, & nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw,nrf3_sfwter,nrf3_vpwter, & nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,nrf2_wspd10m, & nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas, & nrf2_uwnd10m,nrf2_vwnd10m if (mype.eq.0) then print*,' -----------in lanczos_anlerr-----------------' print*,'nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv=', & nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv print*,'nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw=', & nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw print*,'nrf3_sfwter,nrf3_vpwter=',nrf3_sfwter,nrf3_vpwter print*,'nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist=',& nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist print*,'nrf2_wspd10m,nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv=', & nrf2_wspd10m,nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv print*,'nrf2_tcamt,nrf2_lcbas=',nrf2_tcamt,nrf2_lcbas print*,'nrf2_uwnd10m,nrf2_vwnd10m=',nrf2_uwnd10m,nrf2_vwnd10m nrf1_rad=-1 if ( jpch*npred > 0 ) nrf1_rad=10 !any positive integer nrf1_prec=-1 if (jtype*npredp > 0 ) nrf1_prec=11 !any positive integer print*,'nrf1_rad,nrf1_prec=',nrf1_rad,nrf1_prec endif close(lun1) !may have used less iterations than original 'nv=niter(1)'. address that possibility: if (.not.skipanlerr) then open (lun1,file='used_iterations.dat',form='formatted') read (lun1,*) nvused nvused=nvused-1 close(lun1) else nvused=0 end if if (mype==0) print*,'nv,nvused,nvmin=',nv,nvused,nvmin nv=min(nvused,30) !added 19Nov2007 call aka_controlvars !assign alternative names to control variables !eg., igust becomes alternative name for nrf2_gust !*********************************************************** !*********************************************************** !==>Divide up full horizontal domain into subdomains call horiz_domain_partition(nlon,nlat,mype,npe, & ista,iend,jsta,jend) ! print*,'ista,iend,jsta,jend=',ista,iend,jsta,jend !*********************************************************** !*********************************************************** if (anlqlty .or. cvbasedrecalibration) then tdfname='sigges' ; call td_flds(nlon,nlat,tdfname,mype,npe) tdfname='siganl' ; call td_flds(nlon,nlat,tdfname,mype,npe) if (cvbasedrecalibration) call td_flds(nlon,nlat,cvbasedcmodel0,mype,npe) call make_cressanl_common(cgrid,ista,iend,jsta,jend,ribuffer_km,mype,npe) call mpi_barrier(mpi_comm_world,ierror) endif !*********************************************************** !*********************************************************** !==>Let every processor know about the number of horizontal ! points and the starting and ending i and j values on each ! subdomain allocate(ijpts_in_sub(npe)) allocate(ista_info(npe)) allocate(iend_info(npe)) allocate(jsta_info(npe)) allocate(jend_info(npe)) allocate(iaux0(npe)) allocate(iaux1(npe)) allocate(iaux2(npe)) allocate(iaux3(npe)) allocate(iaux4(npe)) ijpts_in_sub(:)=0 ; iaux0(:)=0 ista_info(:)=0 ; iaux1(:)=0 iend_info(:)=0 ; iaux2(:)=0 jsta_info(:)=0 ; iaux3(:)=0 jend_info(:)=0 ; iaux4(:)=0 iaux0(mype+1)=(iend-ista+1)*(jend-jsta+1) iaux1(mype+1)=ista iaux2(mype+1)=iend iaux3(mype+1)=jsta iaux4(mype+1)=jend call mpi_allreduce(iaux0,ijpts_in_sub,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux1,ista_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux2,iend_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux3,jsta_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(iaux4,jend_info,npe,mpi_integer, & mpi_sum,mpi_comm_world,ierror) if (mype.eq.0) then print*,'=====================================================================' print*,' task ijpts_in_sub ista_info iend_info jsta_info jend_info' do i=1,npe print'(i5,9x,i5,4(7x,i5))',i-1,ijpts_in_sub(i),ista_info(i),iend_info(i), & jsta_info(i),jend_info(i) enddo print*,'=====================================================================' endif !*********************************************************** !==>Total number of horizontal slabs in one gradient file !excluding gradient with respect to rad or prec num_slabs=0 if (nrf3_sf > 0 ) num_slabs = num_slabs + nsig if (nrf3_vp > 0 ) num_slabs = num_slabs + nsig if (nrf3_t > 0 ) num_slabs = num_slabs + nsig if (nrf3_q > 0 ) num_slabs = num_slabs + nsig if (nrf3_oz > 0 ) num_slabs = num_slabs + nsig if (nrf3_cw > 0 ) num_slabs = num_slabs + nsig if (nrf3_sfwter > 0 ) num_slabs = num_slabs + nsig if (nrf3_vpwter > 0 ) num_slabs = num_slabs + nsig if (nrf2_ps > 0 ) num_slabs = num_slabs + 1 if (nrf2_sst > 0 ) num_slabs = num_slabs + 1 if (nrf2_gust > 0 ) num_slabs = num_slabs + 1 if (nrf2_vis > 0 ) num_slabs = num_slabs + 1 if (nrf2_pblh > 0 ) num_slabs = num_slabs + 1 if (nrf2_dist > 0 ) num_slabs = num_slabs + 1 if (nrf2_wspd10m > 0 ) num_slabs = num_slabs + 1 if (nrf2_td2m > 0 ) num_slabs = num_slabs + 1 if (nrf2_mxtm > 0 ) num_slabs = num_slabs + 1 if (nrf2_mitm > 0 ) num_slabs = num_slabs + 1 if (nrf2_pmsl > 0 ) num_slabs = num_slabs + 1 if (nrf2_howv > 0 ) num_slabs = num_slabs + 1 if (nrf2_tcamt > 0 ) num_slabs = num_slabs + 1 if (nrf2_lcbas > 0 ) num_slabs = num_slabs + 1 if (nrf2_uwnd10m > 0 ) num_slabs = num_slabs + 1 if (nrf2_vwnd10m > 0 ) num_slabs = num_slabs + 1 if (mype .eq. 0) & print*,'# of horizontal slabs in gradient file=',num_slabs !*********************************************************** !*********************************************************** skip_analysis_error: if (.not. skipanlerr) then Lanczos_Method: if (nv >= nvmin .and. (.not. climatology)) then !*********************************************************** !*********************************************************** !*********************************************************** !==>prepare for alltoallv if (mod(num_slabs,npe) == 0) then num_pad=num_slabs else num_pad=(num_slabs/npe +1)*npe endif if (mype .eq. 0) print*,'num_pad=',num_pad allocate(hall_loc(ista:iend,jsta:jend,num_pad,0:nv)) allocate(gall_loc(ista:iend,jsta:jend,num_pad,0:nv)) if (nrf1_rad > 0) then allocate(hrad(jpch*npred,0:nv)) allocate(grad(jpch*npred,0:nv)) endif if (nrf1_prec > 0) then allocate(hprec(jtype*npredp,0:nv)) allocate(gprec(jtype*npredp,0:nv)) endif allocate(iscnt(npe)) allocate(isdisp(npe)) allocate(ircnt(npe)) allocate(irdisp(npe)) do i=1,npe iscnt(i)=ijpts_in_sub(i) enddo isdisp(1)=0 do i=1,npe if (i /=1 ) isdisp(i)=isdisp(i-1)+iscnt(i-1) enddo do i=1,npe ircnt(i)=ijpts_in_sub(mype+1) irdisp(i)=(i-1)*ijpts_in_sub(mype+1) enddo !*********************************************************** !*********************************************************** !==>for each grad field, have all processors read in the ! data slabs. evoke alltoallv to send fields to subdomains ! everytime a group of npe slabs have been read in or when ! end of slab count has been reached allocate(h_slab(nlon,nlat)) allocate(hfield(nlon*nlat)) allocate(h_rad(jpch*npred)) allocate(h_prec(jtype*npredp)) allocate(g_slab(nlon,nlat)) allocate(gfield(nlon*nlat)) allocate(g_rad(jpch*npred)) allocate(g_prec(jtype*npredp)) allocate(alpha(0:nv)) allocate(gamma(-1:nv-1)) write (clun1(1:2),'(i2.2)') jiter do 400 n=0,nv ! print*,'gradient field number=',n write (clun2(1:3),'(i3.3)') n open (lun1,file='grady.dat_'//clun1//'_'//clun2,form='unformatted') open (lun2,file='gradx.dat_'//clun1//'_'//clun2,form='unformatted') rewind (lun1) rewind (lun2) read (lun1) nlon1,nlat1,nsig1,jpch1,npred1, & jtype1,npredp1,jiter1,nv1,alpha(n),gamma(n-1) read (lun2) nlon1,nlat1,nsig1,jpch1,npred1, & jtype1,npredp1,jiter1,nv1,alpha(n),gamma(n-1) ! print*,'nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1=',& ! nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1 islab_prev=1 do 200 islab=1,num_slabs if(mype==mod(islab-1,npe)) then read(lun1) h_slab call vectorform(h_slab,hfield,nlon,nlat, & ista_info,iend_info,jsta_info,jend_info,npe) read(lun2) g_slab call vectorform(g_slab,gfield,nlon,nlat, & ista_info,iend_info,jsta_info,jend_info,npe) ! print*,'min,max,g_slab,gfield=', & ! minval(g_slab),minval(gfield),maxval(g_slab),maxval(gfield) ! print*,'min,max,h_slab,hfield=', & ! minval(h_slab),minval(hfield),maxval(h_slab),maxval(hfield) else read(lun1) read(lun2) end if if(mod(islab,npe) == 0 .or. islab == num_slabs) then call mp_flush(1) call mpi_alltoallv(hfield,iscnt,isdisp,mpi_real8, & hall_loc(ista,jsta,islab_prev,n),ircnt,irdisp, & mpi_real8,mpi_comm_world,ierror) call mpi_alltoallv(gfield,iscnt,isdisp,mpi_real8, & gall_loc(ista,jsta,islab_prev,n),ircnt,irdisp, & mpi_real8,mpi_comm_world,ierror) islab_prev=islab+1 end if 200 continue read(lun1) h_rad read(lun1) h_prec read(lun2) g_rad read(lun2) g_prec if (nrf1_rad > 0) then hrad(:,n)=h_rad(:) grad(:,n)=g_rad(:) endif if (nrf1_prec > 0) then hprec(:,n)=h_prec(:) gprec(:,n)=g_prec(:) endif close (lun1) close (lun2) 400 continue !*********************************************************** !*********************************************************** !==>compute renormalization coefficients cnorm isub=(iend-ista+1) jsub=(jend-jsta+1) ijsub=isub*jsub allocate(slab1(1:isub,1:jsub)) allocate(slab2(1:isub,1:jsub)) allocate(hgdotaux(0:nv,0:nv)) allocate(hgdot(0:nv,0:nv)) hgdotaux(0:nv,0:nv)=0._r_double do n=0,nv do m=0,nv do islab=1,num_slabs do j=1,jsub do i=1,isub slab1(i,j)=hall_loc(i+ista-1,j+jsta-1,islab,n) slab2(i,j)=gall_loc(i+ista-1,j+jsta-1,islab,m) enddo enddo hgdotaux(n,m)=hgdotaux(n,m)+innerprod(slab1,slab2,ijsub) ! print*,'n,m,islab,min,max,slab1,slab2=', & ! n,m,islab,minval(slab1),minval(slab2),maxval(slab1),maxval(slab2) enddo enddo enddo hgdot(0:nv,0:nv)=0._r_double call mpi_allreduce(hgdotaux,hgdot,(nv+1)*(nv+1),mpi_real8, & mpi_sum,mpi_comm_world,ierror) if (nrf1_rad > 0) then do n=0,nv do m=0,nv hgdot(n,m)=hgdot(n,m)+innerprod(hrad(1,n),grad(1,m),jpch*npred) enddo enddo endif if (nrf1_prec > 0) then do n=0,nv do m=0,nv hgdot(n,m)=hgdot(n,m)+innerprod(hprec(1,n),gprec(1,m),jtype*npredp) enddo enddo endif if (lanczos_writediag .and. mype.eq.0) then do n=0,nv do m=0,nv print*,'n,m,hgdot=',n,m,hgdot(n,m) enddo enddo endif trouble=.false. allocate(cnorm(0:nv)) do n=0,nv if (hgdot(n,n) .le. tiny_r_kind*1._r_double) then trouble=.true. if (mype==0) then print*,'n,hgdot(n,n)=',n,hgdot(n,n) print*,'trouble with normalization--must abort' endif else cnorm(n)=1._r_double/sqrt(hgdot(n,n)) endif enddo if (trouble) then call mpi_finalize(ierror) stop endif if (mype .eq. 0) print*,'min,max,cnorm=',minval(cnorm),maxval(cnorm) deallocate(ijpts_in_sub) deallocate(ista_info) deallocate(iend_info) deallocate(jsta_info) deallocate(jend_info) deallocate(iaux0) deallocate(iaux1) deallocate(iaux2) deallocate(iaux3) deallocate(iaux4) deallocate(iscnt) deallocate(isdisp) deallocate(ircnt) deallocate(irdisp) deallocate(h_slab) deallocate(g_slab) deallocate(hfield) deallocate(gfield) deallocate(h_rad) deallocate(g_rad) deallocate(h_prec) deallocate(g_prec) deallocate(hgdot) deallocate(hgdotaux) !****************************************************************************** ! step-II : build the T matrix * !****************************************************************************** allocate(beta(0:nv)) allocate(delta(0:nv)) allocate(theta(0:nv)) beta(0:nv)=zero delta(0:nv)=zero theta(0:nv)=zero delta(0)=1./alpha(0) do n=0,nv if (n.lt.nv) then delta(n+1)=one/alpha(n+1) + gamma(n)/alpha(n) beta(n+1)=-cnorm(n)/cnorm(n+1)/alpha(n) endif if (n.gt.0) theta(n)=-gamma(n-1)/alpha(n-1)*cnorm(n)/cnorm(n-1) enddo if (mype.eq.0) then print*,'delta(0)=',delta(0) do n=1,nv print*,'n,delta(n),beta(n),theta(n)=',n,delta(n),beta(n),theta(n) enddo endif allocate(tmatrix(1:nv,1:nv)) tmatrix(:,:)=zero do n=1,nv tmatrix(n,n)=delta(n-1) if (n.lt.nv) tmatrix(n,n+1)=theta(n) if (n.gt.1) tmatrix(n,n-1)=beta(n-1) enddo do n=1,nv do m=n+1,nv tmatrix(n,m)=tmatrix(m,n) enddo enddo !****************************************************************************** ! step-III : compute eigenvectors and eigenvalues of matrix !****************************************************************************** nsym=(nv*nv-nv)/2+nv allocate(atest(nsym)) j=0 do n=1,nv do i=1,n j=j+1 atest(j)=tmatrix(i,n) enddo enddo if (mype.eq.0) print*,'nsym,j=',nsym,j allocate(rtest(1:nv,1:nv)) rtest(:,:)=zero call eigen_3(atest,rtest,nv,0) allocate(eval(nv)) j=0 do n=1,nv do i=1,n j=j+1 if (i.eq.n) eval(n)=atest(j) enddo if (mype==0) print*,'n,eval=',n,eval(n) enddo ! check orthogonality of eigenvectors do n=1,nv do m=1,nv dotp=zero do k=1,nv dotp=dotp+rtest(k,m)*rtest(k,n) enddo if (mype==0) print*,'orthogonality of T eigvectors, m,n,dotp=',m,n,dotp enddo enddo !****************************************************************************** ! step-V : compute the eigenvectors of the Hessian. Store in gall_loc ! grad, and gprec !****************************************************************************** allocate(hall_loc2(ista:iend,jsta:jend,num_slabs,1:nv)) allocate(gall_loc2(ista:iend,jsta:jend,num_slabs,1:nv)) hall_loc2(ista:iend,jsta:jend,1:num_slabs,1:nv)=zero gall_loc2(ista:iend,jsta:jend,1:num_slabs,1:nv)=zero do n=1,nv do islab=1,num_slabs do j=jsta,jend do i=ista,iend do k=1,nv hall_loc2(i,j,islab,n)=hall_loc2(i,j,islab,n)+ & hall_loc(i,j,islab,k-1)*cnorm(k-1)*rtest(k,n) gall_loc2(i,j,islab,n)=gall_loc2(i,j,islab,n)+ & gall_loc(i,j,islab,k-1)*cnorm(k-1)*rtest(k,n) enddo enddo enddo enddo enddo if (nrf1_rad > 0) then allocate(hrad2(jpch*npred,1:nv)) allocate(grad2(jpch*npred,1:nv)) hrad2(:,:)=zero grad2(:,:)=zero do n=1,nv do j=1,jpch*npred do k=1,nv hrad2(j,n)=hrad2(j,n)+hrad(j,k-1)*cnorm(k-1)*rtest(k,n) grad2(j,n)=grad2(j,n)+grad(j,k-1)*cnorm(k-1)*rtest(k,n) enddo enddo enddo endif if (nrf1_prec > 0) then allocate(hprec2(jtype*npredp,1:nv)) allocate(gprec2(jtype*npredp,1:nv)) hprec2(:,:)=zero gprec2(:,:)=zero do n=1,nv do j=1,jtype*npredp do k=1,nv hprec2(j,n)=hprec2(j,n)+hprec(j,k-1)*cnorm(k-1)*rtest(k,n) gprec2(j,n)=gprec2(j,n)+gprec(j,k-1)*cnorm(k-1)*rtest(k,n) enddo enddo enddo endif !****************************************************************************** ! step-VI : construct vector of error variance. Store in ! hall_loc(:,:,:,0), hrad(:,0), and hprec(:,0) !****************************************************************************** do n=1,nv dotp0=0._r_double do islab=1,num_slabs do j=1,jsub do i=1,isub slab1(i,j)=hall_loc2(i+ista-1,j+jsta-1,islab,n) slab2(i,j)=gall_loc2(i+ista-1,j+jsta-1,islab,n) enddo enddo dotp0=dotp0+innerprod(slab1,slab2,ijsub) enddo dotp=0._r_double call mpi_allreduce(dotp0,dotp,1,mpi_real8, & mpi_sum,mpi_comm_world,ierror) if (nrf1_rad > 0) then dotp=dotp+innerprod(hrad2(1,n),grad2(1,n),jpch*npred) endif if (nrf1_prec > 0) then dotp=dotp+innerprod(hprec2(1,n),gprec2(1,n),jtype*npredp) endif if (dotp .gt. tiny_r_kind*1._r_double) then hall_loc2(ista:iend,jsta:jend,1:num_slabs,n)= & hall_loc2(ista:iend,jsta:jend,1:num_slabs,n)/sqrt(dotp) gall_loc2(ista:iend,jsta:jend,1:num_slabs,n)= & gall_loc2(ista:iend,jsta:jend,1:num_slabs,n)/sqrt(dotp) if (nrf1_rad > 0) then hrad2(:,n)=hrad2(:,n)/sqrt(dotp) grad2(:,n)=grad2(:,n)/sqrt(dotp) endif if (nrf1_prec > 0) then hprec2(:,n)=hprec2(:,n)/sqrt(dotp) gprec2(:,n)=gprec2(:,n)/sqrt(dotp) endif else if (mype==0) print*,'WARNING: for Hessian, getting null eigevector #',n endif enddo !test L2-orthogonality of eigenvectors of the hessian do n=1,nv do m=1,nv dotp0=0._r_double do islab=1,num_slabs do j=1,jsub do i=1,isub slab1(i,j)=hall_loc2(i+ista-1,j+jsta-1,islab,n) slab2(i,j)=gall_loc2(i+ista-1,j+jsta-1,islab,m) enddo enddo dotp0=dotp0+innerprod(slab1,slab2,ijsub) enddo dotp=0._r_double call mpi_allreduce(dotp0,dotp,1,mpi_real8, & mpi_sum,mpi_comm_world,ierror) if (mype.eq.0) print*,'n,m,dotprod=',n,m,dotp enddo enddo !==> Perform variable transformation to satisfy specific needs (RTMA, NMM, GFS). For ! example, RTMA needs u and v variances. Hence, convert eigenvectors of hessian from ! (psi,chi) to (u,v). allocate(psi(ista:iend,jsta:jend,nsig)) allocate(chi(ista:iend,jsta:jend,nsig)) allocate(t(ista:iend,jsta:jend,nsig)) allocate(q(ista:iend,jsta:jend,nsig)) allocate(oz(ista:iend,jsta:jend,nsig)) allocate(cwmr(ista:iend,jsta:jend,nsig)) allocate(sfcp(ista:iend,jsta:jend)) allocate(sfct(ista:iend,jsta:jend)) allocate(u(ista:iend,jsta:jend,nsig,nv)) allocate(v(ista:iend,jsta:jend,nsig,nv)) allocate(u2(ista:iend,jsta:jend,nsig,nv)) allocate(v2(ista:iend,jsta:jend,nsig,nv)) allocate(wdir(ista:iend,jsta:jend,nsig,nv)) allocate(wdir2(ista:iend,jsta:jend,nsig,nv)) allocate(uwter(ista:iend,jsta:jend,nsig,nv)) allocate(vwter(ista:iend,jsta:jend,nsig,nv)) allocate(u2wter(ista:iend,jsta:jend,nsig,nv)) allocate(v2wter(ista:iend,jsta:jend,nsig,nv)) allocate(wdirwter(ista:iend,jsta:jend,nsig,nv)) allocate(wdir2wter(ista:iend,jsta:jend,nsig,nv)) if (nrf3_sfwter > 0) allocate(sfwter(ista:iend,jsta:jend,nsig)) if (nrf3_vpwter > 0) allocate(vpwter(ista:iend,jsta:jend,nsig)) if (nrf2_gust > 0) allocate(gust(ista:iend,jsta:jend)) if (nrf2_vis > 0) allocate(vis(ista:iend,jsta:jend)) if (nrf2_pblh > 0) allocate(pblh(ista:iend,jsta:jend)) if (nrf2_dist > 0) allocate(dist(ista:iend,jsta:jend)) if (nrf2_wspd10m > 0) allocate(wspd10m(ista:iend,jsta:jend)) if (nrf2_td2m > 0) allocate(td2m(ista:iend,jsta:jend)) if (nrf2_mxtm > 0) allocate(mxtm(ista:iend,jsta:jend)) if (nrf2_mitm > 0) allocate(mitm(ista:iend,jsta:jend)) if (nrf2_pmsl > 0) allocate(pmsl(ista:iend,jsta:jend)) if (nrf2_howv > 0) allocate(howv(ista:iend,jsta:jend)) if (nrf2_tcamt > 0) allocate(tcamt(ista:iend,jsta:jend)) if (nrf2_lcbas > 0) allocate(lcbas(ista:iend,jsta:jend)) if (nrf2_uwnd10m > 0) allocate(uwnd10m(ista:iend,jsta:jend)) if (nrf2_vwnd10m > 0) allocate(vwnd10m(ista:iend,jsta:jend)) ista2=max(1,ista-1) iend2=min(nx,iend+1) jsta2=max(1,jsta-1) jend2=min(ny,jend+1) allocate(tempa(nlon,nlat)) allocate ( coeffx (jsta2:jend2,ista2:iend2) ) !transposed allocate ( coeffy (jsta2:jend2,ista2:iend2) ) !transposed allocate ( glon (ista2:iend2,jsta2:jend2) ) !non-transposed allocate ( rlats (ista:iend,jsta:jend) ) !non-transposed allocate ( rlons (ista:iend,jsta:jend) ) !non-transposed open (33,file='rtma_latlon_mpfactor.dat',form='unformatted') read(33) tempa !glat if (mype==0) print*,'in lanczos_anlerr: glat,min,max=',minval(tempa),maxval(tempa) rlats(ista:iend,jsta:jend)=tempa(ista:iend,jsta:jend) read(33) tempa !glon if (mype==0) print*,'in lanczos_anlerr: glon,min,max=',minval(tempa),maxval(tempa) rlons(ista:iend,jsta:jend)=tempa(ista:iend,jsta:jend) do j=jsta2,jend2 do i=ista2,iend2 glon(i,j)=tempa(i,j) if (glon(i,j) < 0._r_single) & glon(i,j)=glon(i,j)+360._r_single! convert to eastern longitude enddo enddo read(33) tempa !mapfact do j=jsta2,jend2 do i=ista2,iend2 coeffx(j,i)=half/(ds0*tempa(i,j)*one) coeffy(j,i)=coeffx(j,i) enddo enddo close(33) amin8=minval(one/coeffx) ; amax8=maxval(one/coeffx) call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror) if (mype==0) then print*,'in lanczos_anlerr: 1/coeffx,min,max=',amin82,amax82 print*,'in lanczos_anlerr: 1/coeffy,min,max=',amin82,amax82 endif dg2rad=atan(one)/45._r_kind if (trim(cgrid)=='alaska') xn=one if (trim(cgrid)=='akhres') xn=one if (trim(cgrid)=='juneau') xn=one if (trim(cgrid)=='conus') xn=sin(r25*dg2rad) if (trim(cgrid)=='cohres') xn=sin(r25*dg2rad) if (trim(cgrid)=='hrrr') xn=sin(r38_5*dg2rad) if (trim(cgrid)=='cohresext') xn=sin(r25*dg2rad) if (trim(cgrid)=='cohreswexp') xn=sin(r25*dg2rad) if (trim(cgrid)=='alaska') elonv=r210 if (trim(cgrid)=='akhres') elonv=r210 if (trim(cgrid)=='juneau') elonv=r225 if (trim(cgrid)=='conus') elonv=r265 if (trim(cgrid)=='cohres') elonv=r265 if (trim(cgrid)=='hrrr') elonv=r262_5 if (trim(cgrid)=='cohresext') elonv=r265 if (trim(cgrid)=='cohreswexp') elonv=r265 if (leigenvector_write) then open (14,file='hessian_eigenvectors.dat',form='unformatted') write(14) nv write(14) eval allocate(tempb8(nlon,nlat)) allocate(tempc8(nlon,nlat)) endif do n=1,nv kall=0 do ifield=1,8 if (ifield==5 .and. nrf3_oz <= 0) cycle if (ifield==6 .and. nrf3_cw <= 0) cycle if (ifield==7 .and. nrf3_sfwter <= 0) cycle if (ifield==8 .and. nrf3_vpwter <= 0) cycle if (ifield.eq.1) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend psi(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.2) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend chi(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.3) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend t(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.4) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend q(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.5) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend oz(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.6) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend cwmr(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.7) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend sfwter(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif if (ifield.eq.8) then do k=1,nsig kall=kall+1 do j=jsta,jend do i=ista,iend vpwter(i,j,k)=hall_loc2(i,j,kall,n) enddo enddo enddo endif enddo kall=kall+1 do j=jsta,jend do i=ista,iend sfcp(i,j)=hall_loc2(i,j,kall,n) enddo enddo if (nrf2_sst > 0 ) then kall=kall+1 sfct(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_gust > 0 ) then kall=kall+1 gust(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_vis > 0 ) then kall=kall+1 vis(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_pblh > 0 ) then kall=kall+1 pblh(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_dist > 0 ) then kall=kall+1 dist(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_wspd10m > 0 ) then kall=kall+1 wspd10m(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_td2m > 0 ) then kall=kall+1 td2m(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_mxtm > 0 ) then kall=kall+1 mxtm(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_mitm > 0 ) then kall=kall+1 mitm(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_pmsl > 0 ) then kall=kall+1 pmsl(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_howv > 0 ) then kall=kall+1 howv(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_tcamt > 0 ) then kall=kall+1 tcamt(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_lcbas > 0 ) then kall=kall+1 lcbas(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_uwnd10m > 0 ) then kall=kall+1 uwnd10m(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif if (nrf2_vwnd10m > 0 ) then kall=kall+1 vwnd10m(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n) endif ! note: rad and prec parts are hrad2(:,n) and hprec2(:,n) allocate(uaux1(ista:iend,jsta:jend)) allocate(vaux1(ista:iend,jsta:jend)) allocate(uaux2(ista:iend,jsta:jend)) allocate(vaux2(ista:iend,jsta:jend)) do k=1,nsig ! (psi,chi) ==> (u,v) conversion call psichi2uv_eigenvec(psi(:,:,k),chi(:,:,k), & coeffx,coeffy,glon,uaux1,vaux1,uaux2,vaux2, & dg2rad,xn,elonv,k,n,ista,iend,jsta,jend, & ista2,iend2,jsta2,jend2,nlon,nlat,cgrid,mype,npe,lanczos_writediag) do j=jsta,jend do i=ista,iend u(i,j,k,n)=uaux1(i,j) v(i,j,k,n)=vaux1(i,j) wdir(i,j,k,n)=(u(i,j,k,n)/vbar-ubar/vbar**2*v(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad u2(i,j,k,n)=uaux2(i,j) v2(i,j,k,n)=vaux2(i,j) wdir2(i,j,k,n)=(u2(i,j,k,n)/vbar-ubar/vbar**2*v2(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad enddo enddo enddo if (nrf3_sfwter > 0 .and. nrf3_vpwter > 0) then allocate ( rmask (ista:iend,jsta:jend) ) !non-transposed open (33,file='rtma_slmask.dat',form='unformatted') read(33) tempa rmask(ista:iend,jsta:jend)=tempa(ista:iend,jsta:jend) close(33) do k=1,nsig call psichi2uv_eigenvec(sfwter(:,:,k),vpwter(:,:,k), & coeffx,coeffy,glon,uaux1,vaux1,uaux2,vaux2, & dg2rad,xn,elonv,k,n,ista,iend,jsta,jend, & ista2,iend2,jsta2,jend2,nlon,nlat,cgrid,mype,npe,lanczos_writediag) do j=jsta,jend do i=ista,iend uwter(i,j,k,n)=uaux1(i,j) vwter(i,j,k,n)=vaux1(i,j) wdirwter(i,j,k,n)=(uwter(i,j,k,n)/vbar-ubar/vbar**2*vwter(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad u2wter(i,j,k,n)=uaux2(i,j) v2wter(i,j,k,n)=vaux2(i,j) wdir2wter(i,j,k,n)=(u2wter(i,j,k,n)/vbar-ubar/vbar**2*v2wter(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad enddo enddo call landlake_uvmerge_V2(uaux1,vaux1,u(:,:,k,n),v(:,:,k,n),uwter(:,:,k,n),vwter(:,:,k,n), & rmask,rlats,rlons,ista,iend,jsta,jend,nsig,1) u(:,:,k,n)=uaux1(:,:) v(:,:,k,n)=vaux1(:,:) call landlake_uvmerge_V2(uaux1,vaux1,u2(:,:,k,n),v2(:,:,k,n),u2wter(:,:,k,n),v2wter(:,:,k,n), & rmask,rlats,rlons,ista,iend,jsta,jend,nsig,1) u2(:,:,k,n)=uaux1(:,:) v2(:,:,k,n)=vaux1(:,:) enddo deallocate (rmask) endif deallocate(uaux1) deallocate(vaux1) deallocate(uaux2) deallocate(vaux2) if (leigenvector_write) then do ifield=1,30 if ( ifield == 5 .and. nrf3_oz <= 0 ) cycle if ( ifield == 6 .and. nrf3_cw <= 0 ) cycle if ( ifield == 7 .and. nrf3_sfwter <= 0 ) cycle if ( ifield == 8 .and. nrf3_vpwter <= 0 ) cycle if ( ifield == 10 .and. nrf2_sst <= 0 ) cycle if ( ifield == 17 .and. nrf2_gust <= 0 ) cycle if ( ifield == 18 .and. nrf2_vis <= 0 ) cycle if ( ifield == 19 .and. nrf2_pblh <= 0 ) cycle if ( ifield == 20 .and. nrf2_dist <= 0 ) cycle if ( ifield == 21 .and. nrf2_wspd10m<= 0 ) cycle if ( ifield == 22 .and. nrf2_td2m <= 0 ) cycle if ( ifield == 23 .and. nrf2_mxtm <= 0 ) cycle if ( ifield == 24 .and. nrf2_mitm <= 0 ) cycle if ( ifield == 25 .and. nrf2_pmsl <= 0 ) cycle if ( ifield == 26 .and. nrf2_howv <= 0 ) cycle if ( ifield == 27 .and. nrf2_tcamt <= 0 ) cycle if ( ifield == 28 .and. nrf2_lcbas <= 0 ) cycle if ( ifield == 29 .and. nrf2_uwnd10m<= 0 ) cycle if ( ifield == 30 .and. nrf2_vwnd10m<= 0 ) cycle nk=nsig if (ifield==9 .or. ifield==10 .or. ifield > 16) nk=1 do k=1,nk tempb8=zero if (ifield == 1) tempb8(ista:iend,jsta:jend) = psi(ista:iend,jsta:jend,k) if (ifield == 2) tempb8(ista:iend,jsta:jend) = chi(ista:iend,jsta:jend,k) if (ifield == 3) tempb8(ista:iend,jsta:jend) = t(ista:iend,jsta:jend,k) if (ifield == 4) tempb8(ista:iend,jsta:jend) = q(ista:iend,jsta:jend,k) if (ifield == 5) tempb8(ista:iend,jsta:jend) = oz(ista:iend,jsta:jend,k) if (ifield == 6) tempb8(ista:iend,jsta:jend) = cwmr(ista:iend,jsta:jend,k) if (ifield == 7) tempb8(ista:iend,jsta:jend) = sfwter(ista:iend,jsta:jend,k) if (ifield == 8) tempb8(ista:iend,jsta:jend) = vpwter(ista:iend,jsta:jend,k) if (ifield == 9) tempb8(ista:iend,jsta:jend) = sfcp(ista:iend,jsta:jend) if (ifield == 10) tempb8(ista:iend,jsta:jend) = sfct(ista:iend,jsta:jend) if (ifield == 11) tempb8(ista:iend,jsta:jend) = u(ista:iend,jsta:jend,k,n) if (ifield == 12) tempb8(ista:iend,jsta:jend) = v(ista:iend,jsta:jend,k,n) if (ifield == 13) tempb8(ista:iend,jsta:jend) = wdir(ista:iend,jsta:jend,k,n) if (ifield == 14) tempb8(ista:iend,jsta:jend) = u2(ista:iend,jsta:jend,k,n) if (ifield == 15) tempb8(ista:iend,jsta:jend) = v2(ista:iend,jsta:jend,k,n) if (ifield == 16) tempb8(ista:iend,jsta:jend) = wdir2(ista:iend,jsta:jend,k,n) if (ifield == 17) tempb8(ista:iend,jsta:jend) = gust(ista:iend,jsta:jend) if (ifield == 18) tempb8(ista:iend,jsta:jend) = vis(ista:iend,jsta:jend) if (ifield == 19) tempb8(ista:iend,jsta:jend) = pblh(ista:iend,jsta:jend) if (ifield == 20) tempb8(ista:iend,jsta:jend) = dist(ista:iend,jsta:jend) if (ifield == 21) tempb8(ista:iend,jsta:jend) = wspd10m(ista:iend,jsta:jend) if (ifield == 22) tempb8(ista:iend,jsta:jend) = td2m(ista:iend,jsta:jend) if (ifield == 23) tempb8(ista:iend,jsta:jend) = mxtm(ista:iend,jsta:jend) if (ifield == 24) tempb8(ista:iend,jsta:jend) = mitm(ista:iend,jsta:jend) if (ifield == 25) tempb8(ista:iend,jsta:jend) = pmsl(ista:iend,jsta:jend) if (ifield == 26) tempb8(ista:iend,jsta:jend) = howv(ista:iend,jsta:jend) if (ifield == 27) tempb8(ista:iend,jsta:jend) = tcamt(ista:iend,jsta:jend) if (ifield == 28) tempb8(ista:iend,jsta:jend) = lcbas(ista:iend,jsta:jend) if (ifield == 29) tempb8(ista:iend,jsta:jend) = uwnd10m(ista:iend,jsta:jend) if (ifield == 30) tempb8(ista:iend,jsta:jend) = vwnd10m(ista:iend,jsta:jend) tempc8=zero call mpi_reduce(tempb8,tempc8,nlon*nlat,mpi_real8,& mpi_sum,0,mpi_comm_world,ierror) if (mype==0) write(14) tempc8 enddo enddo endif enddo if (leigenvector_write) then deallocate(tempb8) deallocate(tempc8) close(14) endif !==>Perform final matrix multiplication to get variances hall_loc(ista:iend,jsta:jend,1:num_pad,0)=zero do islab=1,num_slabs do j=jsta,jend do i=ista,iend do k=1,nv if (eval(k).gt.tiny_r_kind) then hall_loc(i,j,islab,0)=hall_loc(i,j,islab,0)+ & hall_loc2(i,j,islab,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hall_loc2(i,j,islab,k)*float(k*k) ! hall_loc2(i,j,islab,k)*(one/eval(k)-one)*hall_loc2(i,j,islab,k) !exact endif enddo enddo enddo enddo if (nrf1_rad > 0) then hrad(:,0)=zero do j=1,jpch*npred do k=1,nv hrad(j,0)=hrad(j,0)+hrad2(j,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hrad2(j,k)*float(k*k) ! hrad(j,0)=hrad(j,0)+hrad2(j,k)*(one/eval(k)-one)*hrad2(j,k) !exact enddo enddo endif if (nrf1_prec > 0) then hprec(:,0)=zero do j=1,jtype*npredp do k=1,nv hprec(j,0)=hprec(j,0)+hprec2(j,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hprec2(j,k)*float(k*k) ! hprec(j,0)=hprec(j,0)+hprec2(j,k)*(one/eval(k)-one)*hprec2(j,k) !exact enddo enddo endif allocate(usum(ista:iend,jsta:jend,1:nsig)) allocate(vsum(ista:iend,jsta:jend,1:nsig)) allocate(wdirsum(ista:iend,jsta:jend,1:nsig)) allocate(usum2(ista:iend,jsta:jend,1:nsig)) allocate(vsum2(ista:iend,jsta:jend,1:nsig)) allocate(wdirsum2(ista:iend,jsta:jend,1:nsig)) usum=zero vsum=zero wdirsum=zero usum2=zero vsum2=zero wdirsum2=zero do kk=1,nsig do j=jsta,jend do i=ista,iend do k=1,nv if (eval(k).gt.tiny_r_kind) then usum(i,j,kk)=usum(i,j,kk)+ & u(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*u(i,j,kk,k)*float(k*k) ! u(i,j,kk,k)*(one/eval(k)-one)*u(i,j,kk,k) !exact vsum(i,j,kk)=vsum(i,j,kk)+ & v(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*v(i,j,kk,k)*float(k*k) ! v(i,j,kk,k)*(one/eval(k)-one)*v(i,j,kk,k) !exact wdirsum(i,j,kk)=wdirsum(i,j,kk)+ & wdir(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*wdir(i,j,kk,k)*float(k*k) ! wdir(i,j,kk,k)*(one/eval(k)-one)*wdir(i,j,kk,k) !exact usum2(i,j,kk)=usum2(i,j,kk)+ & u2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*u2(i,j,kk,k)*float(k*k) ! u2(i,j,kk,k)*(one/eval(k)-one)*u2(i,j,kk,k) !exact vsum2(i,j,kk)=vsum2(i,j,kk)+ & v2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*v2(i,j,kk,k)*float(k*k) ! v2(i,j,kk,k)*(one/eval(k)-one)*v2(i,j,kk,k) !exact wdirsum2(i,j,kk)=wdirsum2(i,j,kk)+ & wdir2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*wdir2(i,j,kk,k)*float(k*k) ! wdir2(i,j,kk,k)*(one/eval(k)-one)*wdir2(i,j,kk,k) !exact endif enddo enddo enddo enddo !****************************************************************************** ! step-VIII : write out solution (valid only for twodvar_regional) !****************************************************************************** call create_errs_common(ista,iend,jsta,jend) !allocates error fields, eg: psierr,qerr do j=jsta,jend do i=ista,iend uerr(i,j)=usum(i,j,1) verr(i,j)=vsum(i,j,1) wdirerr(i,j)=wdirsum(i,j,1) uerr2(i,j)=usum2(i,j,1) verr2(i,j)=vsum2(i,j,1) wdirerr2(i,j)=wdirsum2(i,j,1) enddo enddo uerr=sqrt(max(tiny_single,-uerr)) verr=sqrt(max(tiny_single,-verr)) wdirerr=sqrt(max(tiny_single,-wdirerr)) uerr2=sqrt(max(tiny_single,-uerr2)) verr2=sqrt(max(tiny_single,-verr2)) wdirerr2=sqrt(max(tiny_single,-wdirerr2)) npointer(:)=-1 n=0 if (nrf3_sf > 0 ) then ; n=n+1 ; npointer(1) = n ; endif if (nrf3_vp > 0 ) then ; n=n+1 ; npointer(2) = n ; endif if (nrf3_t > 0 ) then ; n=n+1 ; npointer(3) = n ; endif if (nrf3_q > 0 ) then ; n=n+1 ; npointer(4) = n ; endif if (nrf3_oz > 0 ) then ; n=n+1 ; npointer(5) = n ; endif if (nrf3_cw > 0 ) then ; n=n+1 ; npointer(6) = n ; endif if (nrf3_sfwter > 0 ) then ; n=n+1 ; npointer(7) = n ; endif if (nrf3_vpwter > 0 ) then ; n=n+1 ; npointer(8) = n ; endif if (nrf2_ps > 0 ) then ; n=n+1 ; npointer(9) = n ; endif if (nrf2_sst > 0 ) then ; n=n+1 ; npointer(10) = n ; endif if (nrf2_gust > 0 ) then ; n=n+1 ; npointer(11) = n ; endif if (nrf2_vis > 0 ) then ; n=n+1 ; npointer(12) = n ; endif if (nrf2_pblh > 0 ) then ; n=n+1 ; npointer(13) = n ; endif if (nrf2_dist > 0 ) then ; n=n+1 ; npointer(14) = n ; endif if (nrf2_wspd10m > 0 ) then ; n=n+1 ; npointer(15) = n ; endif if (nrf2_td2m > 0 ) then ; n=n+1 ; npointer(16) = n ; endif if (nrf2_mxtm > 0 ) then ; n=n+1 ; npointer(17) = n ; endif if (nrf2_mitm > 0 ) then ; n=n+1 ; npointer(18) = n ; endif if (nrf2_pmsl > 0 ) then ; n=n+1 ; npointer(19) = n ; endif if (nrf2_howv > 0 ) then ; n=n+1 ; npointer(20) = n ; endif if (nrf2_tcamt > 0 ) then ; n=n+1 ; npointer(21) = n ; endif if (nrf2_lcbas > 0 ) then ; n=n+1 ; npointer(22) = n ; endif if (nrf2_uwnd10m > 0 ) then ; n=n+1 ; npointer(23) = n ; endif if (nrf2_vwnd10m > 0 ) then ; n=n+1 ; npointer(24) = n ; endif if (mype==0) open (15,file='errvar.dat',form='unformatted') allocate(tempb(nlon,nlat)) do n=1,num_slabs tempa=0._r_single tempb=0._r_single do j=jsta,jend do i=ista,iend tempb(i,j)=hall_loc(i,j,n,0) enddo enddo call mpi_allreduce(tempb,tempa,nlon*nlat,mpi_real4,& mpi_sum,mpi_comm_world,ierror) if (mype==0) write(15) tempa tempa=max(tiny_single,-tempa) tempa=sqrt(tempa) if (n==1) psierr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (n==2) chierr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (n==3) terr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (n==4) qerr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (n==npointer(9)) perr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_gust > 0 .and. n==npointer(11)) gusterr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_vis > 0 .and. n==npointer(12)) viserr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_pblh > 0 .and. n==npointer(13)) pblherr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_dist > 0 .and. n==npointer(14)) disterr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_wspd10m > 0 .and. n==npointer(15)) wspd10merr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_td2m > 0 .and. n==npointer(16)) td2merr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_mxtm > 0 .and. n==npointer(17)) mxtmerr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_mitm > 0 .and. n==npointer(18)) mitmerr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_pmsl > 0 .and. n==npointer(19)) pmslerr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_howv > 0 .and. n==npointer(20)) howverr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_tcamt > 0 .and. n==npointer(21)) tcamterr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_lcbas > 0 .and. n==npointer(22)) lcbaserr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_uwnd10m > 0 .and. n==npointer(23)) uwnd10merr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) if (nrf2_vwnd10m > 0 .and. n==npointer(24)) vwnd10merr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend) enddo if (mype==0) close(15) call error_conversion(cgrid,cvbasedrecalibration,cvbasedcmodel0) filename='errfield.dat' call writeout_errfields(filename) deallocate(tempa) deallocate(tempb) deallocate(glon) deallocate(rlats) deallocate(rlons) call destroy_errs_common !deallocates error fields deallocate(hall_loc) deallocate(gall_loc) deallocate(hall_loc2) deallocate(gall_loc2) if (nrf1_rad > 0) then deallocate(hrad) deallocate(hrad2) deallocate(grad) deallocate(grad2) endif if (nrf1_prec > 0) then deallocate(hprec) deallocate(hprec2) deallocate(gprec) deallocate(gprec2) endif deallocate(alpha) deallocate(gamma) deallocate(beta) deallocate(delta) deallocate(theta) deallocate(cnorm) deallocate(tmatrix) deallocate(atest) deallocate(rtest) deallocate(eval) deallocate(slab1) deallocate(slab2) deallocate(psi) deallocate(chi) deallocate(t) deallocate(q) deallocate(oz) deallocate(cwmr) deallocate(sfcp) deallocate(sfct) deallocate(u) deallocate(v) deallocate(u2) deallocate(v2) deallocate(coeffx) deallocate(coeffy) deallocate(usum) deallocate(vsum) deallocate(usum2) deallocate(vsum2) deallocate(wdir) deallocate(wdir2) deallocate(wdirsum) deallocate(wdirsum2) if (nrf2_gust > 0) deallocate(gust) if (nrf2_vis > 0) deallocate(vis) if (nrf2_pblh > 0) deallocate(pblh) if (nrf2_dist > 0) deallocate(dist) if (nrf2_wspd10m > 0) deallocate(wspd10m) if (nrf2_td2m > 0) deallocate(td2m) if (nrf2_mxtm > 0) deallocate(mxtm) if (nrf2_mitm > 0) deallocate(mitm) if (nrf2_pmsl > 0) deallocate(pmsl) if (nrf2_howv > 0) deallocate(howv) if (nrf2_tcamt > 0) deallocate(tcamt) if (nrf2_lcbas > 0) deallocate(lcbas) if (nrf2_uwnd10m > 0) deallocate(uwnd10m) if (nrf2_vwnd10m > 0) deallocate(vwnd10m) !****************************************************************************** !==> force use of climatological error when explicitly so specified or when ! number of gradient vectors is less than a chosen nvmin value !****************************************************************************** else Lanczos_Method !use climatological error allocate(climerr(nlon,nlat)) if (mype==0) then print*,' either number of iterations used in minimization is less than ', nvmin, & 'or climatology is set to .true.' print*,'skip lanczos procedure and use climatological error instead' open (16,file='errfield.dat_clim',form='unformatted') open (17,file='errfield.dat',form='unformatted') do n=1,11 read(16) climerr write(17) climerr enddo if (nrf2_gust > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_vis > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_pblh > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_dist > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_wspd10m > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_td2m > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_mxtm > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_mitm > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_pmsl > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_howv > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_tcamt > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_lcbas > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_uwnd10m > 0) then ; read(16) climerr ; write(17) climerr ; endif if (nrf2_vwnd10m > 0) then ; read(16) climerr ; write(17) climerr ; endif close(16) close(17) endif deallocate(climerr) call mpi_barrier(mpi_comm_world,ierror) call derive_xbvar_mpi(cgrid) end if Lanczos_Method else if (mype==0) write(6,'(A,L)') 'In lanczos_anlerr: Skipping creating analysis error since skipanlerr=',skipanlerr end if skip_analysis_error !****************************************************************************** ! step-IX : create grib2 files for ges, anl, and anlerr !****************************************************************************** call cnv_to_grib2_mpi(cgrid,skipanlerr) call mpi_barrier(mpi_comm_world,ierror) !****************************************************************************** ! step-X : generate observation text files + create dynamic reject lists !****************************************************************************** if (npe==1) then call get_ob_lists(miter0,cgrid,mype,npe,mpi_comm_world) call create_rjlist(cgrid,mype,npe,mpi_comm_world) else if ((npe-(num_slabs-1)) >= 5) then mypecutoff=num_slabs-1 else mypecutoff=npe/2-1 endif if (mype <= mypecutoff) then icolor=1 ikey=mype else icolor=2 ikey=(mype-(mypecutoff+1)) endif call mpi_comm_split(mpi_comm_world,icolor,ikey,mpi_comm_new,ierror) call mpi_comm_size(mpi_comm_new,npenew,ierror) call mpi_comm_rank(mpi_comm_new,mypenew,ierror) print*,'mpi_comm_new,npenew,mypenew=',mpi_comm_new,npenew,mypenew call mpi_barrier(mpi_comm_world,ierror) if (icolor==1) call get_ob_lists(miter0,cgrid,mypenew,npenew,mpi_comm_new) if (icolor==2) call create_rjlist(cgrid,mypenew,npenew,mpi_comm_new) call mpi_barrier(mpi_comm_world,ierror) endif !****************************************************************************** ! step-XI : perform optional cressman analysis and/or cross-validation !****************************************************************************** if (anlqlty) call anl_quality(rjbuffer_km,ds0,lobjanl) if (anlqlty .or. cvbasedrecalibration) call destroy_cressanl_common call stn_interpolator(cgrid) end subroutine lanczos_anlerr !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ subroutine psichi2uv_eigenvec(psi,chi,coeffx,coeffy,glon, & uaux1,vaux1,uaux2,vaux2, & dg2rad,xn,elonv,k,n,ista,iend,jsta,jend, & ista2,iend2,jsta2,jend2,nlon,nlat,cgrid,mype,npe, & lanczos_writediag) !**************************************************************************** ! abstract: convert (psi,chi) eigenvectors into (u,v) eigenvectors * ! * ! program history log: * ! 2005-10-08 pondeca * ! * !**************************************************************************** use mpi use kinds, only: i_kind,r_single,r_kind use constants, only: zero,one implicit none ! Declare passed variables character(60),intent(in):: cgrid integer(i_kind),intent(in):: mype,npe integer(i_kind),intent(in):: nlon,nlat !global dimensions integer(i_kind),intent(in):: ista,iend,jsta,jend !tile dimensions integer(i_kind),intent(in):: ista2,iend2,jsta2,jend2 !memory dimensions integer(i_kind),intent(in):: k !vertical level # integer(i_kind),intent(in):: n !eigenvector number real(r_kind),intent(in)::dg2rad,xn,elonv real(r_single),intent(in):: glon(ista2:iend2,jsta2:jend2) real(r_kind),intent(in):: psi(ista:iend,jsta:jend) real(r_kind),intent(in):: chi(ista:iend,jsta:jend) real(r_kind),intent(in):: coeffx(jsta2:jend2,ista2:iend2) !transposed real(r_kind),intent(in):: coeffy(jsta2:jend2,ista2:iend2) !transposed logical, intent(in):: lanczos_writediag real(r_kind),intent(out):: uaux1(ista:iend,jsta:jend) real(r_kind),intent(out):: vaux1(ista:iend,jsta:jend) real(r_kind),intent(out):: uaux2(ista:iend,jsta:jend) real(r_kind),intent(out):: vaux2(ista:iend,jsta:jend) ! Declare local variables integer(i_kind) i,j,ierror real(r_kind) angle2,sinx2,cosx2 real(r_kind) amin8,amax8,amin82,amax82 real(r_kind),allocatable,dimension(:,:):: dtemp1,dtemp2 real(r_kind),allocatable,dimension(:,:):: daux1,daux2,daux3,daux4 ! !***************************************************************************** if (n==1) then open (10,file='psi_slabs.dat',form='unformatted') open (20,file='chi_slabs.dat',form='unformatted') open (11,file='u_slabs.dat',form='unformatted') open (21,file='v_slabs.dat',form='unformatted') else open (10,file='psi_slabs.dat',form='unformatted',position='append') open (20,file='chi_slabs.dat',form='unformatted',position='append') open (11,file='u_slabs.dat',form='unformatted',position='append') open (21,file='v_slabs.dat',form='unformatted',position='append') endif allocate(dtemp1(nlon,nlat)) allocate(dtemp2(nlon,nlat)) allocate(daux1(jsta2:jend2,ista2:iend2)) !transposed allocate(daux2(jsta2:jend2,ista2:iend2)) !transposed allocate(daux3(jsta2:jend2,ista2:iend2)) !transposed allocate(daux4(jsta2:jend2,ista2:iend2)) !transposed dtemp1(:,:)=zero do j=jsta,jend do i=ista,iend dtemp1(i,j)=psi(i,j) enddo enddo call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, & mpi_sum,mpi_comm_world,ierror) do j=jsta2,jend2 do i=ista2,iend2 daux1(j,i)=dtemp2(i,j) !holds transposed psi in memory subdomains enddo enddo !optional !------------------------------------------------------------------------ if (lanczos_writediag .and. mype==0) then print*,'n,k,psi,min,max=',n,k,minval(dtemp2),maxval(dtemp2) write(10) sngl(dtemp2) endif !------------------------------------------------------------------------ dtemp1(:,:)=zero do j=jsta,jend do i=ista,iend dtemp1(i,j)=chi(i,j) enddo enddo call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, & mpi_sum,mpi_comm_world,ierror) do j=jsta2,jend2 do i=ista2,iend2 daux2(j,i)=dtemp2(i,j) !holds transposed chi in memory subdomains enddo enddo !optional !------------------------------------------------------------------------ if (lanczos_writediag .and. mype==0) then print*,'n,k,chi,min,max=',n,k,minval(dtemp2),maxval(dtemp2) write(20) sngl(dtemp2) endif !------------------------------------------------------------------------ call psichi2uv_reg_V2(daux1,daux2,daux3,daux4,coeffx,coeffy, & 1, nlat, 1, nlon, 1, 1, & jsta2, jend2, ista2, iend2, 1, 1, & jsta , jend , ista , iend , 1, 1) deallocate(daux1) deallocate(daux2) allocate(daux1(ista2:iend2,jsta2:jend2)) !non-transposed allocate(daux2(ista2:iend2,jsta2:jend2)) !non-transposed daux1=transpose(daux3) daux2=transpose(daux4) !optional !------------------------------------------------------------------------ if (lanczos_writediag) then dtemp1(:,:)=zero do j=jsta,jend do i=ista,iend dtemp1(i,j)=daux1(i,j) enddo enddo call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, & mpi_sum,mpi_comm_world,ierror) if (mype==0) print*,'n,k,u,min,max=',n,k,minval(dtemp2),maxval(dtemp2) if (mype==0) write(11) sngl(dtemp2) dtemp1(:,:)=zero do j=jsta,jend do i=ista,iend dtemp1(i,j)=daux2(i,j) enddo enddo call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, & mpi_sum,mpi_comm_world,ierror) if (mype==0) print*,'n,k,v,min,max=',n,k,minval(dtemp2),maxval(dtemp2) if (mype==0) write(21) sngl(dtemp2) endif !------------------------------------------------------------------------ !get earth relative u,v components of eigenvectors and store in daux3,daux4 deallocate(daux3) deallocate(daux4) allocate(daux3(ista2:iend2,jsta2:jend2)) !non-transposed allocate(daux4(ista2:iend2,jsta2:jend2)) !non-transposed daux3=daux1 daux4=daux2 if (trim(cgrid)=='conus' .or. trim(cgrid)=='alaska' .or. & trim(cgrid)=='cohres' .or. trim(cgrid)=='akhres' .or. trim(cgrid)=='hrrr' .or. & trim(cgrid)=='juneau' .or. trim(cgrid)=='cohresext' .or. trim(cgrid)=='cohreswexp') then do j=jsta2,jend2 do i=ista2,iend2 angle2=xn*(glon(i,j)*one-elonv)*dg2rad sinx2=sin(angle2) cosx2 = cos(angle2) daux3(i,j)=cosx2*daux1(i,j)-sinx2*daux2(i,j) daux4(i,j)=sinx2*daux1(i,j)+cosx2*daux2(i,j) enddo enddo endif !optional !------------------------------------------------------------------------ if (lanczos_writediag) then amin8=minval(daux3) ; amax8=maxval(daux3) call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'n,k,uearth,min,max=',n,k,amin82,amax82 amin8=minval(daux4) ; amax8=maxval(daux4) call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'n,k,vearth,min,max=',n,k,amin82,amax82 !------------------------------------------------------------------------ endif do j=jsta,jend do i=ista,iend uaux1(i,j)=daux1(i,j) vaux1(i,j)=daux2(i,j) uaux2(i,j)=daux3(i,j) vaux2(i,j)=daux4(i,j) enddo enddo close(10) close(20) close(11) close(21) deallocate(daux1) deallocate(daux2) deallocate(daux3) deallocate(daux4) deallocate(dtemp1) deallocate(dtemp2) end subroutine psichi2uv_eigenvec !****************************************************************************** !****************************************************************************** !=========================================================================== subroutine landlake_uvmerge_V2(u,v,uland,vland,uwter,vwter, & isli2,rlats,rlons,is,ie,js,je,nsig,iflg) !$$$ subprogram documentation block ! . . . ! subprogram: landlake_uvmerge merge ! prgmmr: pondeca org: np22 date: 2013-07-15 ! ! abstract: merge land and lake u,v increments together (iflg=1) or perform ! the adjoint operation (iflg=0) ! ! program history log: ! 2013-07-15 pondeca ! ! input argument list: ! u - u grid values ! v - v grid values ! uland - land u grid values ! vland - land v grid values ! uwter - lake u grid values ! vwter - lake v grid values ! isli2 - land/sea mask ! rlats - latitude values ! rlons - longitude values ! is,ie,js,je,nsig - tile dimensions ! iflg - flag for u,v manipulation ! 0: u,v --> uland,vland,uwter,vwter ! 1: uland,vland,uwter,vwter --> u,v ! ! output argument list: ! u - u grid values ! v - v grid values ! uland - land u grid values ! vland - land v grid values ! uwter - lake u grid values ! vwter - lake v grid values ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_single,i_kind implicit none ! Declare passed variables integer(i_kind) ,intent(in ) :: is,ie,js,je,nsig integer(i_kind) ,intent(in ) :: iflg real(r_single),dimension(is:ie,js:je),intent(in) :: isli2,rlats,rlons real(r_kind),dimension(is:ie,js:je,nsig),intent(inout) :: u,v,uland,vland,uwter,vwter ! Declare local variables integer(i_kind) i,j,k real(r_single) flon,flat logical glerlarea ! Declare local parameters ! Great Lakes real(r_single),parameter::flon1=-93. real(r_single),parameter::flon2=-75. real(r_single),parameter::flat1=40.5 real(r_single),parameter::flat2=49.5 ! Great Salt Lake real(r_single),parameter::slon1=-113. real(r_single),parameter::slon2=-112. real(r_single),parameter::slat1=40.6 real(r_single),parameter::slat2=41.7 if (iflg == 0) then do j=js,je do i=is,ie flat=rlats(i,j) flon=rlons(i,j) glerlarea=(flat>=flat1.and.flat<=flat2).and.(flon>=flon1.and.flon<=flon2) glerlarea=glerlarea.or.((flat>=slat1.and.flat<=slat2).and.(flon>=slon1.and.flon<=slon2)) if(isli2(i,j) >=0.5_r_single .or. .not.glerlarea) then do k=1,nsig !note: this subroutine is always called with nsig=1 uland(i,j,k)=+u(i,j,k) vland(i,j,k)=+v(i,j,k) enddo else do k=1,nsig uwter(i,j,k)=+u(i,j,k) vwter(i,j,k)=+v(i,j,k) enddo end if enddo enddo else do j=js,je do i=is,ie flat=rlats(i,j) flon=rlons(i,j) glerlarea=(flat>=flat1.and.flat<=flat2).and.(flon>=flon1.and.flon<=flon2) glerlarea=glerlarea.or.((flat>=slat1.and.flat<=slat2).and.(flon>=slon1.and.flon<=slon2)) if(isli2(i,j) >=0.5_r_single .or. .not.glerlarea) then do k=1,nsig u(i,j,k)=uland(i,j,k) v(i,j,k)=vland(i,j,k) enddo else do k=1,nsig u(i,j,k)=uwter(i,j,k) v(i,j,k)=vwter(i,j,k) enddo end if enddo enddo endif end subroutine landlake_uvmerge_V2 !****************************************************************************** !****************************************************************************** subroutine horiz_domain_partition(nx,ny,mype,npe, & ista,iend,jsta,jend) use kinds, only:i_kind implicit none integer(i_kind) nx,ny,mype,npe,ista,iend,jsta,jend jsta=1 jend=ny call para_range2(1,nx,npe,mype,ista,iend) return end subroutine horiz_domain_partition !****************************************************************************** subroutine para_range2(n1,n2,nprocs,irank,ista,iend) use kinds, only:i_kind implicit none integer(i_kind) n1,n2,nprocs,irank,ista,iend integer(i_kind) iwork1,iwork2 iwork1 = (n2-n1+1) / nprocs iwork2 = mod(n2 - n1 + 1 , nprocs) ista = irank * iwork1 + n1 + min(irank, iwork2) iend = ista + iwork1 -1 if (iwork2 > irank) iend = iend+1 return end subroutine para_range2 !****************************************************************************** subroutine para_range3(n1,n2,nparts,ista,iend) use kinds, only:i_kind implicit none integer(i_kind),intent(in) :: n1,n2,nparts integer(i_kind),intent(out) :: ista(nparts),iend(nparts) integer(i_kind) iwork1,iwork2,n iwork1 = (n2-n1+1) / nparts iwork2 = mod(n2 - n1 + 1 , nparts) do n=0,nparts-1 ista(n+1) = n * iwork1 + n1 + min(n, iwork2) iend(n+1) = ista(n+1) + iwork1 -1 if (iwork2 > n) iend(n+1) = iend(n+1)+1 enddo return end subroutine para_range3 !****************************************************************************** subroutine vectorform(g,v,nx,ny,istart,iend,jstart,jend,nprocs) use kinds, only:r_kind,i_kind implicit none ! Declare passed variables integer(i_kind) nx,ny,nprocs integer(i_kind) istart(nprocs) integer(i_kind) iend(nprocs) integer(i_kind) jstart(nprocs) integer(i_kind) jend(nprocs) real(r_kind) g(nx,ny), v(nx*ny) ! Declare local variables integer(i_kind) i,j,ij,n ij=0 do n=1,nprocs do j=jstart(n),jend(n) do i=istart(n),iend(n) ij=ij+1 v(ij)=g(i,j) enddo enddo enddo ! print*,'in vectorform, ij=',ij return end subroutine vectorform !****************************************************************************** real(r_double) function innerprod(dx,dy,n) !$$$ subprogram documentation block ! . . . . ! subprogram: innerprod calculates dot product for data on subdomain ! prgmmr: derber org: np23 date: 2004-05-13 ! ! abstract: calculates dot product for data on subdomain. Note loops over ! streamfunction, velocity potential, temperature, etc. Also, only ! over interior of subdomain. ! ! program history log: ! 2004-05-13 derber, document ! 2004-07-28 treadon - add only on use declarations; add intent in/out ! ! input argument list: ! dx - input vector 1 ! dy - input vector 2 ! ! output argument list ! dplev - dot product ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_double,i_kind use constants, only: zero implicit none ! Declare passed variables integer(i_kind) n real(r_kind),dimension(n),intent(in)::dx,dy ! Declare local variables integer(i_kind) i innerprod=0._r_double do i=1,n innerprod=innerprod+dx(i)*dy(i) end do return end function innerprod !****************************************************************************** subroutine ax_equal_lambda_bx(a,b,evec,eval,n) use kinds, only:r_kind,r_double,i_kind use constants, only:tiny_r_kind implicit none integer(i_kind) n,i,j real(r_kind) a(n,n),b(n,n),evec(n,n),eval(n) real(r_kind) beta(n),aux(3*n) !CHECK DIMENSION OF AUX complex(r_kind) alpha(n) complex(r_kind),allocatable,dimension(:,:)::z allocate(z(n,n)) call dgegv(1,a,n,b,n,alpha,beta,z,n,n,aux,3*n) do j=1,n if (abs(beta(j)) .gt. tiny_r_kind) then eval(j)=dreal(alpha(j))/beta(j) !CHECK DREAL else print*,'trouble in dgegv. eval=infinity for j=',j eval(j)=-9999._r_kind endif do i=1,n evec(i,j)=dreal(z(i,j)) !CHECK DREAL enddo enddo deallocate(z) return end subroutine ax_equal_lambda_bx !****************************************************************************** !****************************************************************************** subroutine error_conversion(cgrid,cvbasedrecalibration,cvbasedcmodel0) !================================================================================== !Note: The number of error fields has increased so much with the inclusion of ! optional control variables that a bundle-type approach would make the ! code easier to read. consider for future work /MPondeca, 18Jul2014 !================================================================================== use mpi use mpitaskmod, only: mype,npe use mpitaskmod, only: i1 => ista use mpitaskmod, only: i2 => iend use mpitaskmod, only: j1 => jsta use mpitaskmod, only: j2 => jend use mpitaskmod, only: nx => nlon use mpitaskmod, only: ny => nlat use constants, only: tiny_single use controlvars, only: igust,ivis,ipblh,idist, & iwspd10m,itd2m,imxtm,imitm, & ipmsl,ihowv,itcamt,ilcbas,iuwnd10m,ivwnd10m use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, & wspderr,wdirerr,wdirerr2,terr,tderr, & qerr,perr,gusterr,viserr,pblherr,disterr, & wspd10merr,td2merr,mxtmerr,mitmerr,pmslerr, & howverr,tcamterr,lcbaserr,uwnd10merr,vwnd10merr implicit none !declare passed variables character(60),intent(in):: cgrid logical,intent(in):: cvbasedrecalibration character(60),intent(in):: cvbasedcmodel0 !declare local parameters integer(4),parameter::nbckgfls=19 !# of bckg error fields. it's equal to # of ctl variables integer(4),parameter::nerrfls=25 !# of error fields for which eigenvectors !where computed. therefore, excludes wspderr,tderr real(4),parameter:: wdirberr0=15. !(dg) real(4),parameter::a=243.5 real(4),parameter::alpha=440.8 real(4),parameter::b=19.48 real(4),parameter::eps=0.62197 real(4),parameter::perrm=300. real(4),parameter::qerrm=6.e-3 !declare local variables integer(4) i,j,im,ip,jm,jp integer(4) n,ierror integer(4) rtime(6),nlon1,nlat1,nsig1 real(4) amin,amax,aave real(4) amin2,amax2,aave2 real(4) amin3,amax3,aave3 real(4) tempamin,tempamax,fieldmin real(4) wspderrmin,wspderrmax real(4) uberrmin,vberrmin real(4) tderrmin,tderrmax character(60) bckgfname(nbckgfls) real(4),allocatable,dimension(:,:)::psiberr real(4),allocatable,dimension(:,:)::chiberr real(4),allocatable,dimension(:,:)::tberr real(4),allocatable,dimension(:,:)::qberr real(4),allocatable,dimension(:,:)::pberr real(4),allocatable,dimension(:,:)::uberr real(4),allocatable,dimension(:,:)::vberr real(4),allocatable,dimension(:,:)::wdirberr real(4),allocatable,dimension(:,:)::gustberr real(4),allocatable,dimension(:,:)::visberr real(4),allocatable,dimension(:,:)::pblhberr real(4),allocatable,dimension(:,:)::distberr real(4),allocatable,dimension(:,:)::wspd10mberr real(4),allocatable,dimension(:,:)::td2mberr real(4),allocatable,dimension(:,:)::mxtmberr real(4),allocatable,dimension(:,:)::mitmberr real(4),allocatable,dimension(:,:)::pmslberr real(4),allocatable,dimension(:,:)::howvberr real(4),allocatable,dimension(:,:)::tcamtberr real(4),allocatable,dimension(:,:)::lcbasberr real(4),allocatable,dimension(:,:)::uwnd10mberr real(4),allocatable,dimension(:,:)::vwnd10mberr real(4),allocatable,dimension(:,:)::field real(4),allocatable,dimension(:,:)::tempa real(4) qs0,qanl0,panl0 real(4) qv,dqv,e,de,dlne,terrmax real(4) vis0,dist0 logical fexist real(4) psicverr,chicverr,tcverr,pcverr,qcverr,ucverr,vcverr, & wdircverr,u2cverr,v2cverr,wdir2cverr, & wspdcverr,tdcverr, & gustcverr,viscverr,pblhcverr,distcverr, & wspd10mcverr,td2mcverr,mxtmcverr,mitmcverr, & pmslcverr,howvcverr,tcamtcverr,lcbascverr, & uwnd10mcverr,vwnd10mcverr character(7) cfldname(nerrfls) !error field names real(4) hoberr(nerrfls) !roughly half the observation error real(4) avgcverr(nerrfls) !sqrt of the cross-validation rmse real(4) errupper(nerrfls) !the maximum allowable error real(4) hoberrwspd,errupperwspd real(4) hoberrtd,erruppertd real(4) errmax character(60) filename namelist/cverrorupdate/psicverr,chicverr,tcverr,pcverr, & qcverr,ucverr,vcverr,wdircverr,u2cverr,v2cverr,wdir2cverr, & wspdcverr,tdcverr, & gustcverr,viscverr,pblhcverr,distcverr, & wspd10mcverr,td2mcverr,mxtmcverr,mitmcverr, & pmslcverr,howvcverr,tcamtcverr,lcbascverr, & uwnd10mcverr,vwnd10mcverr data cfldname(1) /'psi'/ ; data psicverr /-999./ data cfldname(2) /'chi'/ ; data chicverr /-999./ data cfldname(3) /'t'/ ; data tcverr /2.48/ data cfldname(4) /'p'/ ; data pcverr /2.12/ !ps in hPa data cfldname(5) /'q'/ ; data qcverr /0.0011/ !this is for true sphum (in kg/kg) data cfldname(6) /'u'/ ; data ucverr /2.26/ data cfldname(7) /'v'/ ; data vcverr /3.25/ data cfldname(8) /'wdir'/ ; data wdircverr /20./ data cfldname(9) /'u2'/ ; data u2cverr /2.26/ data cfldname(10) /'v2'/ ; data v2cverr /3.25/ data cfldname(11) /'wdir2'/ ; data wdir2cverr /20./ data cfldname(12) /'gust'/ ; data gustcverr /3.25/ data cfldname(13) /'vis'/ ; data viscverr /3000./ data cfldname(14) /'pblh'/ ; data pblhcverr /-999./ data cfldname(15) /'dist'/ ; data distcverr /3000./ data cfldname(16) /'wspd10m'/; data wspd10mcverr /2.60/ data cfldname(17) /'td2m'/ ; data td2mcverr /5.0/ data cfldname(18) /'mxtm'/ ; data mxtmcverr /2.48/ data cfldname(19) /'mitm'/ ; data mitmcverr /2.48/ data cfldname(20) /'pmsl'/ ; data pmslcverr /3.0/ !pmsl in hPa data cfldname(21) /'howv'/ ; data howvcverr /-999./ data cfldname(22) /'tcamt'/ ; data tcamtcverr /30./ data cfldname(23) /'lcbas'/ ; data lcbascverr /3000./ data cfldname(24) /'uwnd10m'/; data uwnd10mcverr /2.60/ data cfldname(25) /'vwnd10m'/; data vwnd10mcverr /2.60/ data wspdcverr /2.60/ data tdcverr /5.0/ data hoberr(1) /-999./ ; data errupper(1) /-999./ data hoberr(2) /-999./ ; data errupper(2) /-999./ data hoberr(3) /0.5/ ; data errupper(3) /4.0/ data hoberr(4) /0.5/ ; data errupper(4) /3./ !ps in hPa data hoberr(5) /0.0005/ ; data errupper(5) /0.004/ !this is for true sphum in (kg/kg) data hoberr(6) /0.5/ ; data errupper(6) /3.0/ data hoberr(7) /0.5/ ; data errupper(7) /3.0/ data hoberr(8) /5./ ; data errupper(8) /30./ data hoberr(9) /0.5/ ; data errupper(9) /3.0/ data hoberr(10) /0.5/ ; data errupper(10) /3.0/ data hoberr(11) /0.5/ ; data errupper(11) /30./ data hoberr(12) /1.0/ ; data errupper(12) /5.0/ data hoberr(13) /2000./ ; data errupper(13) /4500./ data hoberr(14) /-999./ ; data errupper(14) /-999./ data hoberr(15) /2000./ ; data errupper(15) /4000./ data hoberr(16) /0.5 / ; data errupper(16) /3./ data hoberr(17) /1./ ; data errupper(17) /6./ data hoberr(18) /0.5/ ; data errupper(18) /3./ data hoberr(19) /0.5/ ; data errupper(19) /3./ data hoberr(20) /1.5/ ; data errupper(20) /5./ !pmsl in hPa data hoberr(21) /-999./ ; data errupper(21) /-999./ data hoberr(22) /25./ ; data errupper(22) /50./ data hoberr(23) /2000./ ; data errupper(23) /4000./ data hoberr(24) /0.5 / ; data errupper(24) /3./ data hoberr(25) /0.5 / ; data errupper(25) /3./ data hoberrwspd /0.5/ ; data errupperwspd /3.0/ data hoberrtd /1.0/ ; data erruppertd /6.0/ data bckgfname(1)/'bckgvar.dat_ps'/ data bckgfname(2)/'bckgvar.dat_t'/ data bckgfname(3)/'bckgvar.dat_q'/ data bckgfname(4)/'bckgvar.dat_u'/ data bckgfname(5)/'bckgvar.dat_v'/ data bckgfname(6)/'bckgvar.dat_gust'/ data bckgfname(7)/'bckgvar.dat_vis'/ data bckgfname(8)/'bckgvar.dat_pblh'/ data bckgfname(9)/'bckgvar.dat_dist'/ data bckgfname(10)/'bckgvar.dat_wspd10m'/ data bckgfname(11)/'bckgvar.dat_td2m'/ data bckgfname(12)/'bckgvar.dat_mxtm'/ data bckgfname(13)/'bckgvar.dat_mitm'/ data bckgfname(14)/'bckgvar.dat_pmsl'/ data bckgfname(15)/'bckgvar.dat_howv'/ data bckgfname(16)/'bckgvar.dat_tcamt'/ data bckgfname(17)/'bckgvar.dat_lcbas'/ data bckgfname(18)/'bckgvar.dat_uwnd10m'/ data bckgfname(19)/'bckgvar.dat_vwnd10m'/ !========================================================================================== !==> read in the background error standard deviations. they either come from derive_xbar !========================================================================================== allocate ( psiberr (i1:i2,j1:j2) ) allocate ( chiberr (i1:i2,j1:j2) ) allocate ( tberr (i1:i2,j1:j2) ) allocate ( qberr (i1:i2,j1:j2) ) allocate ( pberr (i1:i2,j1:j2) ) allocate ( uberr (i1:i2,j1:j2) ) allocate ( vberr (i1:i2,j1:j2) ) allocate ( wdirberr (i1:i2,j1:j2) ) if (igust > 0 ) allocate ( gustberr (i1:i2,j1:j2) ) if (ivis > 0 ) allocate ( visberr (i1:i2,j1:j2) ) if (ipblh > 0 ) allocate ( pblhberr (i1:i2,j1:j2) ) if (idist > 0 ) allocate ( distberr (i1:i2,j1:j2) ) if (iwspd10m > 0 ) allocate ( wspd10mberr (i1:i2,j1:j2) ) if (itd2m > 0 ) allocate ( td2mberr (i1:i2,j1:j2) ) if (imxtm > 0 ) allocate ( mxtmberr (i1:i2,j1:j2) ) if (imitm > 0 ) allocate ( mitmberr (i1:i2,j1:j2) ) if (ipmsl > 0 ) allocate ( pmslberr (i1:i2,j1:j2) ) if (ihowv > 0 ) allocate ( howvberr (i1:i2,j1:j2) ) if (itcamt > 0 ) allocate ( tcamtberr (i1:i2,j1:j2) ) if (ilcbas > 0 ) allocate ( lcbasberr (i1:i2,j1:j2) ) if (iuwnd10m > 0 ) allocate ( uwnd10mberr (i1:i2,j1:j2) ) if (ivwnd10m > 0 ) allocate ( vwnd10mberr (i1:i2,j1:j2) ) allocate(field(ny,nx)) !transposed !==> psi open (94,file='bckgvar.dat_psi',form='unformatted') read(94) field do j=j1,j2 do i=i1,i2 psiberr(i,j)=field(j,i) enddo enddo close(94) amin2=minval(field) amax2=maxval(field) aave2=sum(field)/float(nx*ny) !==> chi open (94,file='bckgvar.dat_chi',form='unformatted') read(94) field do j=j1,j2 do i=i1,i2 chiberr(i,j)=field(j,i) enddo enddo close(94) amin3=minval(field) amax3=maxval(field) aave3=sum(field)/float(nx*ny) deallocate(field) !deallocate transposed field allocate(field(nx,ny)) !allocate non-transposed field !============================================================================ !==>dump bckg errors for u,v,t,q,gust,vis,pblh,vis on disc. also dump ! background satutation specific humidity qs !============================================================================ ! call derive_xbvar_mpi(cgrid) ! !============================================================================ !==>read in the background satutation specific humidity qs and compute the ! domain averaged value qs0. it is used to convert the bckg error for ! specific humidity into the background error for pseudo-relative humidty, ! which is the true analysis variable !============================================================================ open (94,file='bckg_qsat.dat',form='unformatted') read(94) field close(94) qs0=sum(field(:,:))/float(nx*ny) !0.01 if (mype==0) print*,'in error_conversion: qs0=',qs0 !============================================================================ !==>read in the bckg errors for u,v,t,q,gust,vis,pblh, and vis. also, ! assume constant value of backg error for wind direction !============================================================================ if (mype==0) print*,'============= in error_conversion ================' if (mype==0) print*,'psiberr,min,max,avg=' , amin2,amax2,aave2 if (mype==0) print*,'chiberr,min,max,avg=' , amin3,amax3,aave3 wdirberr=wdirberr0 !assume const background error do n=1,nbckgfls if (n == 6 .and. igust <= 0 ) cycle if (n == 7 .and. ivis <= 0 ) cycle if (n == 8 .and. ipblh <= 0 ) cycle if (n == 9 .and. idist <= 0 ) cycle if (n == 10 .and. iwspd10m <= 0 ) cycle if (n == 11 .and. itd2m <= 0 ) cycle if (n == 12 .and. imxtm <= 0 ) cycle if (n == 13 .and. imitm <= 0 ) cycle if (n == 14 .and. ipmsl <= 0 ) cycle if (n == 15 .and. ihowv <= 0 ) cycle if (n == 16 .and. itcamt <= 0 ) cycle if (n == 17 .and. ilcbas <= 0 ) cycle if (n == 18 .and. iuwnd10m <= 0 ) cycle if (n == 19 .and. ivwnd10m <= 0 ) cycle if (mype==0) print*,'in error_conversion: n,bckgfname=', & n,trim(bckgfname(n)) open (94,file=trim(bckgfname(n)),form='unformatted') read(94) field close(94) if (n==1) field=field/1000. !get in cbar if (n==3) field=field/qs0 if (n==14) field=field/1000. !get in cbar if (n==1) pberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==2) tberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==3) qberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==4) uberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==5) vberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==6) gustberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==7) visberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==8) pblhberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==9) distberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==10) wspd10mberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==11) td2mberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==12) mxtmberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==13) mitmberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==14) pmslberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==15) howvberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==16) tcamtberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==17) lcbasberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==18) uwnd10mberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) if (n==19) vwnd10mberr (i1:i2,j1:j2) = field (i1:i2,j1:j2) amin=minval(field) amax=maxval(field) aave=sum(field)/float(nx*ny) if (mype==0) then if (n==1) print*,'pberr,min,max,avg=' , amin,amax,aave if (n==2) print*,'tberr,min,max,avg=' , amin,amax,aave if (n==3) print*,'qberr,min,max,avg=' , amin,amax,aave if (n==4) print*,'uberr,min,max,avg=' , amin,amax,aave if (n==5) then print*,'vberr,min,max,avg=' , amin,amax,aave print*,'wdirberr,min,max,avg=',wdirberr0,wdirberr0,wdirberr0 endif if (n==6) print*,'gustberr,min,max,avg=' , amin,amax,aave if (n==7) print*,'visberr,min,max,avg=' , amin,amax,aave if (n==8) print*,'pblhberr,min,max,avg=' , amin,amax,aave if (n==9) print*,'distberr,min,max,avg=' , amin,amax,aave if (n==10) print*,'wspd10mberr,min,max,avg=' , amin,amax,aave if (n==11) print*,'td2mberr,min,max,avg=' , amin,amax,aave if (n==12) print*,'mxtmberr,min,max,avg=' , amin,amax,aave if (n==13) print*,'mitmberr,min,max,avg=' , amin,amax,aave if (n==14) print*,'pmslberr,min,max,avg=' , amin,amax,aave if (n==15) print*,'howvberr,min,max,avg=' , amin,amax,aave if (n==16) print*,'tcamtberr,min,max,avg=' , amin,amax,aave if (n==17) print*,'lcabsberr,min,max,avg=' , amin,amax,aave if (n==18) print*,'uwnd10mberr,min,max,avg=' , amin,amax,aave if (n==19) print*,'vwnd10mberr,min,max,avg=' , amin,amax,aave endif enddo deallocate(field) if (mype==0) print*,'===================================================' !========================================================================================== !==> read in the domain-averaged cross-validation standard deviation error info through ! namelist if available. otherwise, use 'climatological' info from data statements above. ! it's used to rescale the analysis uncertainty. !========================================================================================== inquire(file='cverrorupdate_input',exist=fexist) if (fexist) then if (mype==0) print*,'in error_conversion: read-in renormalizing cross-validation stats from input file' open (55,file='cverrorupdate_input',form='formatted') read(55,cverrorupdate) close(55) endif avgcverr(1) = psicverr avgcverr(2) = chicverr avgcverr(3) = tcverr avgcverr(4) = pcverr avgcverr(5) = qcverr avgcverr(6) = ucverr avgcverr(7) = vcverr avgcverr(8) = wdircverr avgcverr(9) = u2cverr avgcverr(10) = v2cverr avgcverr(11) = wdir2cverr avgcverr(12) = gustcverr avgcverr(13) = viscverr avgcverr(14) = pblhcverr avgcverr(15) = distcverr avgcverr(16) = wspd10mcverr avgcverr(17) = td2mcverr avgcverr(18) = mxtmcverr avgcverr(19) = mitmcverr avgcverr(20) = pmslcverr avgcverr(21) = howvcverr avgcverr(22) = tcamtcverr avgcverr(23) = lcbascverr avgcverr(24) = uwnd10mcverr avgcverr(25) = vwnd10mcverr !convert from mb to cbar. do it for ps and pmsl hoberr(4)=hoberr(4)/10. ; hoberr(20)=hoberr(20)/10. errupper(4)=errupper(4)/10. ; errupper(20)=errupper(20)/10. avgcverr(4)=avgcverr(4)/10. ; avgcverr(20)=avgcverr(20)/10. !convert from sphm error to pseudo-rh error hoberr(5)=hoberr(5)/qs0 errupper(5)=errupper(5)/qs0 avgcverr(5)=avgcverr(5)/qs0 !----------------------------------------------------- ! modification: vis error from vis to g space ! dg/dvis at x = x0 , let us denote it by dgdvis ! if p \= 0.0, dg/dvis = x^(p-1) ! if p = 0.0, dg/dvis = 1/x ! since p \= 0.0, so, hardwired the error ! delta(G) = dgdvis(evaluated at x ref) *delta(VIS) !----------------------------------------------------- !convert of vis error to g-space, p=0.2 !******************************************************* vis0=1500. if (hoberr(13) > 0.) hoberr(13)=hoberr(13)*(vis0**(0.2-1.0)) if (errupper(13) > 0.) errupper(13)=errupper(13)*(vis0**(0.2-1.0)) if (avgcverr(13) > 0.) avgcverr(13)=avgcverr(13)*(vis0**(0.2-1.0)) !convert from dist error from dist space to g-space p=0.1 dist0=1500. if (hoberr(15) > 0.) hoberr(15)=hoberr(15)*(dist0**(0.1-1.0)) if (errupper(15) > 0.) errupper(15)=errupper(15)*(dist0**(0.1-1.0)) if (avgcverr(15) > 0.) avgcverr(15)=avgcverr(15)*(dist0**(0.1-1.0)) do n=1,nerrfls if (n == 12 .and. igust <= 0 ) cycle if (n == 13 .and. ivis <= 0 ) cycle if (n == 14 .and. ipblh <= 0 ) cycle if (n == 15 .and. idist <= 0 ) cycle if (n == 16 .and. iwspd10m <= 0 ) cycle if (n == 17 .and. itd2m <= 0 ) cycle if (n == 18 .and. imxtm <= 0 ) cycle if (n == 19 .and. imitm <= 0 ) cycle if (n == 20 .and. ipmsl <= 0 ) cycle if (n == 21 .and. ihowv <= 0 ) cycle if (n == 22 .and. itcamt <= 0 ) cycle if (n == 23 .and. ilcbas <= 0 ) cycle if (n == 24 .and. iuwnd10m <= 0 ) cycle if (n == 25 .and. ivwnd10m <= 0 ) cycle if (mype==0) & print*,' in error_conversion: hoberr,avgcverr,errupper for ',cfldname(n),':', & hoberr(n),avgcverr(n),errupper(n) enddo !========================================================================================== !==> compute the analysis uncertainty !========================================================================================== allocate(tempa(i1:i2,j1:j2)) allocate(field(i1:i2,j1:j2)) !allocate as a subdomain field do n=1,nerrfls if (n == 12 .and. igust <= 0 ) cycle if (n == 13 .and. ivis <= 0 ) cycle if (n == 14 .and. ipblh <= 0 ) cycle if (n == 15 .and. idist <= 0 ) cycle if (n == 16 .and. iwspd10m <= 0 ) cycle if (n == 17 .and. itd2m <= 0 ) cycle if (n == 18 .and. imxtm <= 0 ) cycle if (n == 19 .and. imitm <= 0 ) cycle if (n == 20 .and. ipmsl <= 0 ) cycle if (n == 21 .and. ihowv <= 0 ) cycle if (n == 22 .and. itcamt <= 0 ) cycle if (n == 23 .and. ilcbas <= 0 ) cycle if (n == 24 .and. iuwnd10m <= 0 ) cycle if (n == 25 .and. ivwnd10m <= 0 ) cycle if (n==1) then ; field=psiberr ; tempa=psierr ; endif if (n==2) then ; field=chiberr ; tempa=chierr ; endif if (n==3) then ; field=tberr ; tempa=terr ; endif if (n==4) then ; field=pberr ; tempa=perr ; endif if (n==5) then ; field=qberr ; tempa=qerr ; endif if (n==6) then ; field=uberr ; tempa=uerr ; endif if (n==7) then ; field=vberr ; tempa=verr ; endif if (n==8) then ; field=wdirberr ; tempa=wdirerr ; endif if (n==9) then ; field=uberr ; tempa=uerr2 ; endif !assume uberr2=uberr if (n==10) then ; field=vberr ; tempa=verr2 ; endif !assume vberr2=vberr if (n==11) then ; field=wdirberr ; tempa=wdirerr2 ; endif !assume wdirberr2=wdirberr if (n==12) then ; field=gustberr ; tempa=gusterr ; endif if (n==13) then ; field=visberr ; tempa=viserr ; endif if (n==14) then ; field=pblhberr ; tempa=pblherr ; endif if (n==15) then ; field=distberr ; tempa=disterr ; endif if (n==16) then ; field=wspd10mberr ; tempa=wspd10merr ; endif if (n==17) then ; field=td2mberr ; tempa=td2merr ; endif if (n==18) then ; field=mxtmberr ; tempa=mxtmerr ; endif if (n==19) then ; field=mitmberr ; tempa=mitmerr ; endif if (n==20) then ; field=pmslberr ; tempa=pmslerr ; endif if (n==21) then ; field=howvberr ; tempa=howverr ; endif if (n==22) then ; field=tcamtberr ; tempa=tcamterr ; endif if (n==23) then ; field=lcbasberr ; tempa=lcbaserr ; endif if (n==24) then ; field=uwnd10mberr ; tempa=uwnd10merr ; endif if (n==25) then ; field=vwnd10mberr ; tempa=vwnd10merr ; endif if (n==13) then !visibility exception. assume anlerr is ! tempa=sqrt(sqrt(sqrt(tempa))) !this function of the 'std deviation' reduction tempa=tempa*1. !this function of the 'std deviation' reduction amax=maxval(tempa) call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror) field=tempamax endif amin=minval(tempa) ; amax=maxval(tempa) call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'in error_conversion: n,tempa,min,max=', & n,tempamin,tempamax tempa=field-tempa !preferred way ! tempa=field-sqrt(tempa) !alternate way ! tempa=sqrt(max(0.000001,(field*field-tempa*tempa))) !most correct way amin=minval(tempa) call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror) if (hoberr(n) > 0.) then tempa=tempa-tempamin+hoberr(n) else tempa(:,:)=max(0.,tempa(:,:)) endif amin=minval(field) call mpi_allreduce(amin,fieldmin,1,mpi_real,mpi_min,mpi_comm_world,ierror) if (avgcverr(n) > 0.) then errmax=max(min(avgcverr(n),errupper(n)),fieldmin) !why use fieldmin and not fieldmax? else errmax=1. endif if (mype==0) print*,'in error_conversion: n, errmax=',n, errmax amax=maxval(tempa) call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror) tempa=tempa/tempamax*errmax if (n==1) psierr = tempa if (n==2) chierr = tempa if (n==3) terr = tempa if (n==4) perr = tempa if (n==5) qerr = tempa if (n==6) uerr = tempa if (n==7) verr = tempa if (n==8) wdirerr = tempa if (n==9) uerr2 = tempa if (n==10) verr2 = tempa if (n==11) wdirerr2 = tempa if (n==12) gusterr = tempa if (n==13) viserr = tempa if (n==14) pblherr = tempa if (n==15) disterr = tempa if (n==16) wspd10merr = tempa if (n==17) td2merr = tempa if (n==18) mxtmerr = tempa if (n==19) mitmerr = tempa if (n==20) pmslerr = tempa if (n==21) howverr = tempa if (n==22) tcamterr = tempa if (n==23) lcbaserr = tempa if (n==24) uwnd10merr = tempa if (n==25) vwnd10merr = tempa enddo amin=minval(uberr) call mpi_allreduce(amin,uberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror) amin=minval(vberr) call mpi_allreduce(amin,vberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror) wspderr=sqrt(uerr*uerr+verr*verr) !a very wild approach amin=minval(wspderr) call mpi_allreduce(amin,wspderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror) wspderr=wspderr-wspderrmin+hoberrwspd errmax=max(min(wspdcverr,errupperwspd),0.5*(uberrmin+vberrmin)) !kinda reasonable if (mype==0) print*,'in error_conversion: errmax for wspd=',errmax amax=maxval(wspderr) call mpi_allreduce(amax,wspderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror) wspderr=wspderr/wspderrmax*errmax perr=perr*1000. !convert to Pa qerr=qerr*qs0 !from pseudo-rh to sphum if ( ipmsl > 0) pmslerr=pmslerr*1000. !convert to Pa if ( ivis > 0 ) viserr=viserr/(vis0**(0.2-1.0)) !from g space to vis space if ( idist > 0 ) disterr=disterr/(dist0**(0.1-1.0)) !from g-space to dist space deallocate(field) !========================================================================================== !==> derive dewpoint temperature error allocate(field(nx,ny)) !allocate as a global field open (53,file='siganl',form='unformatted') read(53) rtime,nlon1,nlat1,nsig1 read(53) field !full record contains glat,dx read(53) field !full record contains glon,dy read(53) field !panl panl0=.5*(maxval(field)+minval(field)) read(53) field !terrain read(53) field !temperture read(53) field !qanl close(53) qanl0=.5*(maxval(abs(field))+minval(abs(field))) if (mype==0) print*,'in error_conversion: qanl0,panl0=',qanl0,panl0 do j=j1,j2 do i=i1,i2 qv=qanl0/(1.-qanl0) dqv=qerr(i,j)/(1.-qanl0)**2 e=panl0/100.*qv/(eps+qv) de=e/panl0*perr(i,j) + & eps*panl0/100./(eps+qv)**2*dqv dlne=de/e tderr(i,j)=(a*b-alpha)/(b-alog(e))**2*dlne tderr(i,j)=abs(tderr(i,j)) enddo enddo amin=minval(tderr) call mpi_allreduce(amin,tderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror) tderr=tderr-tderrmin+hoberrtd errmax=max(min(tdcverr,erruppertd),0.) !replace 0. with mininum of bckg error for td once available if (mype==0) print*,'in error_conversion: errmax for td=',errmax amax=maxval(tderr) call mpi_allreduce(amax,tderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror) tderr=tderr/tderrmax*errmax !just in case ... perr=max(tiny_single,perr) terr=max(tiny_single,terr) tderr=max(tiny_single,tderr) uerr=max(tiny_single,uerr) verr=max(tiny_single,verr) qerr=max(tiny_single,qerr) uerr2=max(tiny_single,uerr2) verr2=max(tiny_single,verr2) if (iwspd10m > 0) wspd10merr=max(tiny_single,wspd10merr) if (itd2m > 0) td2merr=max(tiny_single,td2merr) if (imxtm > 0) mxtmerr=max(tiny_single,mxtmerr) if (imitm > 0) mitmerr=max(tiny_single,mitmerr) if (ipmsl > 0) pmslerr=max(tiny_single,pmslerr) if (ihowv > 0) howverr=max(tiny_single,howverr) if (itcamt > 0) tcamterr=max(tiny_single,tcamterr) if (ilcbas > 0) lcbaserr=max(tiny_single,lcbaserr) if (iuwnd10m > 0) uwnd10merr=max(tiny_single,uwnd10merr) if (ivwnd10m > 0) vwnd10merr=max(tiny_single,vwnd10merr) if (cvbasedrecalibration) then filename='errfield.dat_precval' call writeout_errfields(filename) !write out full-domain anl error fields call cvbasedrecal(cgrid,cvbasedcmodel0) endif deallocate(psiberr) deallocate(chiberr) deallocate(tberr) deallocate(qberr) deallocate(pberr) deallocate(uberr) deallocate(vberr) deallocate(wdirberr) if (igust > 0 ) deallocate(gustberr) if (ivis > 0 ) deallocate(visberr) if (ipblh > 0 ) deallocate(pblhberr) if (idist > 0 ) deallocate(distberr) if (iwspd10m > 0 ) deallocate(wspd10mberr) if (itd2m > 0 ) deallocate(td2mberr) if (imxtm > 0 ) deallocate(mxtmberr) if (imitm > 0 ) deallocate(mitmberr) if (ipmsl > 0 ) deallocate(pmslberr) if (ihowv > 0 ) deallocate(howvberr) if (itcamt > 0 ) deallocate(tcamtberr) if (ilcbas > 0 ) deallocate(lcbasberr) if (iuwnd10m > 0 ) deallocate(uwnd10mberr) if (ivwnd10m > 0 ) deallocate(vwnd10mberr) deallocate(field) deallocate(tempa) end subroutine error_conversion !********************************************************************* !********************************************************************* subroutine writeout_errfields(filename) use mpi use mpitaskmod, only: mype,npe use mpitaskmod, only: i1 => ista use mpitaskmod, only: i2 => iend use mpitaskmod, only: j1 => jsta use mpitaskmod, only: j2 => 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 errs_common, only: uerr,verr,uerr2,verr2, & wspderr,wdirerr,wdirerr2,terr,tderr, & qerr,perr,gusterr,viserr,pblherr,disterr, & wspd10merr,td2merr,mxtmerr,mitmerr,pmslerr, & howverr,tcamterr,lcbaserr, & uwnd10merr,vwnd10merr implicit none !declare passed variables character(60),intent(in):: filename !declare local variables integer(4) nn,ierror real(4),allocatable,dimension(:,:):: slab1,slab2 !--------------------------------------------------------------------- ! if (mype==0)open (17,file=trim(filename),form='unformatted') allocate(slab1(nx,ny)) allocate(slab2(nx,ny)) do nn=1,25 if (nn==12 .and. igust <=0) cycle if (nn==13 .and. ivis <=0) cycle if (nn==14 .and. ipblh <=0) cycle if (nn==15 .and. idist <=0) cycle if (nn==16 .and. iwspd10m <=0) cycle if (nn==17 .and. itd2m <=0) cycle if (nn==18 .and. imxtm <=0) cycle if (nn==19 .and. imitm <=0) cycle if (nn==20 .and. ipmsl <=0) cycle if (nn==21 .and. ihowv <=0) cycle if (nn==22 .and. itcamt <=0) cycle if (nn==23 .and. ilcbas <=0) cycle if (nn==24 .and. iuwnd10m <=0) cycle if (nn==25 .and. ivwnd10m <=0) cycle slab1(:,:)=0. if (nn==1) slab1 (i1:i2,j1:j2) = perr (i1:i2,j1:j2) if (nn==2) slab1 (i1:i2,j1:j2) = terr (i1:i2,j1:j2) if (nn==3) slab1 (i1:i2,j1:j2) = tderr (i1:i2,j1:j2) if (nn==4) slab1 (i1:i2,j1:j2) = uerr (i1:i2,j1:j2) if (nn==5) slab1 (i1:i2,j1:j2) = verr (i1:i2,j1:j2) if (nn==6) slab1 (i1:i2,j1:j2) = qerr (i1:i2,j1:j2) if (nn==7) slab1 (i1:i2,j1:j2) = wdirerr2 (i1:i2,j1:j2) if (nn==8) slab1 (i1:i2,j1:j2) = wspderr (i1:i2,j1:j2) if (nn==9) slab1 (i1:i2,j1:j2) = uerr2 (i1:i2,j1:j2) if (nn==10) slab1 (i1:i2,j1:j2) = verr2 (i1:i2,j1:j2) if (nn==11) slab1 (i1:i2,j1:j2) = wdirerr (i1:i2,j1:j2) if (nn==12) slab1 (i1:i2,j1:j2) = gusterr (i1:i2,j1:j2) if (nn==13) slab1 (i1:i2,j1:j2) = viserr (i1:i2,j1:j2) if (nn==14) slab1 (i1:i2,j1:j2) = pblherr (i1:i2,j1:j2) if (nn==15) slab1( i1:i2,j1:j2) = disterr (i1:i2,j1:j2) if (nn==16) slab1( i1:i2,j1:j2) = wspd10merr (i1:i2,j1:j2) if (nn==17) slab1( i1:i2,j1:j2) = td2merr (i1:i2,j1:j2) if (nn==18) slab1( i1:i2,j1:j2) = mxtmerr (i1:i2,j1:j2) if (nn==19) slab1( i1:i2,j1:j2) = mitmerr (i1:i2,j1:j2) if (nn==20) slab1( i1:i2,j1:j2) = pmslerr (i1:i2,j1:j2) if (nn==21) slab1( i1:i2,j1:j2) = howverr (i1:i2,j1:j2) if (nn==22) slab1( i1:i2,j1:j2) = tcamterr (i1:i2,j1:j2) if (nn==23) slab1( i1:i2,j1:j2) = lcbaserr (i1:i2,j1:j2) if (nn==24) slab1( i1:i2,j1:j2) = uwnd10merr (i1:i2,j1:j2) if (nn==25) slab1( i1:i2,j1:j2) = vwnd10merr (i1:i2,j1:j2) call mpi_allreduce(slab1,slab2,nx*ny, & mpi_real,mpi_sum,mpi_comm_world,ierror) if (mype==0) write(17) slab2 enddo if (mype==0) close(17) call mpi_barrier(mpi_comm_world,ierror) deallocate(slab1) deallocate(slab2) end subroutine writeout_errfields !********************************************************************* !********************************************************************* !======================================================================= !------------------------------------------------------------------------------------- ! SUBROUTINE 'EIGEN' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE. ! ! NOTE: TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS ! (1) HAVE BEEN CHANGED TO (*). ! ! .................................................................. ! ! SUBROUTINE EIGEN_3 ! ! PURPOSE ! COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC ! MATRIX ! ! USAGE ! CALL EIGEN_3(A,R,N,MV) ! ! DESCRIPTION OF PARAMETERS ! A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. ! RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF ! MATRIX A IN DESCENDING ORDER. ! R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, ! IN SAME SEQUENCE AS EIGENVALUES) ! N - ORDER OF MATRICES A AND R ! MV- INPUT CODE ! 0 COMPUTE EIGENVALUES AND EIGENVECTORS ! 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE ! DIMENSIONED BUT MUST STILL APPEAR IN CALLING ! SEQUENCE) ! ! REMARKS ! ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1) ! MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R ! ! SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED ! NONE ! ! METHOD ! DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED ! BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL ! METHODS FOR DIGITAL COMPUTERS', EDITED BY A. RALSTON AND ! H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 ! ! .................................................................. ! SUBROUTINE EIGEN_3(A,R,N,MV) use kinds, only: r_kind,i_kind,r_single,r_double,i_long integer(i_kind) N,MV real(r_kind) A(*),R(*),ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,& COSX2,SINCS,RANGE ! ! ............................................................... ! ! IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE ! C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION ! STATEMENT WHICH FOLLOWS. ! ! DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, ! 1 COSX2,SINCS,RANGE ! ! THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS ! APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS ! ROUTINE. ! ! THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO ! CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. SQRT IN STATEMENTS ! 40, 68, 75, AND 78 MUST BE CHANGED TO DSQRT. ABS IN STATEMENT ! 62 MUST BE CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD ! BE CHANGED TO 1.0D-12. ! ! ............................................................... ! ! GENERATE IDENTITY MATRIX ! 5 RANGE=1.0E-12_r_kind IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0 IF(I-J) 20,15,20 15 R(IJ)=1.0 20 CONTINUE ! ! COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) ! 25 ANORM=0.0 DO 35 I=1,N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM) 165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) ! ! INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR ! IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 ! ! COMPUTE SIN AND COS ! 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF( ABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=0.5*(A(LL)-A(MM)) 68 Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/ SQRT(2.0*(1.0+( SQRT(1.0-Y*Y)))) SINX2=SINX*SINX 78 COSX= SQRT(1.0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX ! ! ROTATE L AND M COLUMNS ! ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X ! ! TESTS FOR COMPLETION ! ! TEST FOR M = LAST COLUMN ! 130 IF(M-N) 135,140,135 135 M=M+1 GO TO 60 ! ! TEST FOR L = SECOND FROM LAST COLUMN ! 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 ! ! COMPARE THRESHOLD WITH FINAL NORM ! 160 IF(THR-ANRMX) 165,165,45 ! ! SORT EIGENVALUES AND EIGENVECTORS ! 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE RETURN END ! !======================================================================= subroutine para_range(n1,n2,nprocs,irank,ista,iend) implicit none integer(4) n1,n2,nprocs,irank,ista,iend integer(4) iwork1,iwork2 iwork1 = (n2-n1+1) / nprocs iwork2 = mod(n2 - n1 + 1 , nprocs) ista = irank * iwork1 + n1 + min(irank, iwork2) iend = ista + iwork1 -1 if (iwork2 > irank) iend = iend+1 return end !------------------------------------------------------ subroutine ob_loc(xob,yob,nobs,nvar) implicit none integer(4) nobs,nvar real(4) xob(nobs) real(4) yob(nobs) integer(4) n,lun character*8 filename character*5 clun1,clun2 character*7 clun11,clun22 if (nvar.eq.1) filename='u' if (nvar.eq.2) filename='v' if (nvar.eq.3) filename='t' if (nvar.eq.4) filename='p' if (nvar.eq.5) filename='q' lun=88 open (lun,file='plot_ob_loc_'//trim(filename)//'.gs', & form='formatted') write(lun,*) "'set display color white'" write(lun,*) "'clear'" write(lun,*) "'open anlerr.des'" write(lun,*) "'set mproj off'" if (trim(filename).eq.'u') write(lun,*) "'display uerr'" if (trim(filename).eq.'v') write(lun,*) "'display verr'" if (trim(filename).eq.'t') write(lun,*) "'display terr'" if (trim(filename).eq.'p') write(lun,*) "'display perr'" if (trim(filename).eq.'q') write(lun,*) "'display qerr'" write(lun,*) write(lun,*) do 100 n=1,nobs if (n.le.9) then write (clun1,"(i1.1)") n write (clun2,"(i1.1)") n endif if (n.ge.10.and.n.le.99) then write (clun1,"(i2.2)") n write (clun2,"(i2.2)") n endif if (n.ge.100.and.n.le.999) then write (clun1,"(i3.3)") n write (clun2,"(i3.3)") n endif if (n.ge.1000.and.n.le.9999) then write (clun1,"(i4.4)") n write (clun2,"(i4.4)") n endif if (n.ge.10000.and.n.le.99999) then write (clun1,"(i5.5)") n write (clun2,"(i5.5)") n endif if (n.ge.100000.and.n.le.999999) then write (clun1,"(i6.6)") n write (clun2,"(i6.6)") n endif if (n.ge.1000000.and.n.le.9999999) then write (clun1,"(i7.7)") n write (clun2,"(i7.7)") n endif clun11='a.'//trim(clun1) clun22='b.'//trim(clun2) write(lun,"(1x,a7,3x,'=',3x,f12.6)") clun11,xob(n) write(lun,"(1x,a7,3x,'=',3x,f12.6)") clun22,yob(n) 100 continue write(lun,*) write(lun,*) write(lun,*) 'll=1' write(lun,*) 'while ( ll <',nobs+1,')' write(lun,*) " 'q gr2xy 'a.ll' 'b.ll" write(lun,*) ' line=sublin(result,1)' write(lun,*) ' xval=subwrd(line,3)' write(lun,*) ' yval=subwrd(line,6)' write(lun,*) " 'draw mark 3 'xval' 'yval' 0.06'" write(lun,*) ' ll=ll+1' write(lun,*) ' endwhile' return end !------------------------------------------------------ subroutine smther_one(g1,is,ie,js,je,ns) ! apply 1-2-1 smoother in each direction of data slab ! use kinds,only: r_single,i_long implicit none integer(4) is, ie, js, je integer(4) i,j,l,ip,im,jp,jm integer(4), intent(in) :: ns real(4), dimension(is:ie, js:je), intent(inout) :: g1 ! on input: original data slab ! on ouput: filtered data slab real(4), allocatable:: g2(:,:) allocate(g2(is:ie,js:je)) do l=1,ns do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=.25*(g1(ip,j)+g1(im,j))+.5*g1(i,j) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=.25*(g2(i,jp)+g2(i,jm))+.5*g2(i,j) end do end do end do deallocate(g2) return end subroutine smther_one !------------------------------------------------------ subroutine fetch_anlincs(itime,pinc,tinc,qinc, & uinc,vinc,tdinc, & is,ie,js,je,nx,ny,mype) use mpi implicit none !Declare passed variables integer(4),intent(in):: is,ie,js,je integer(4),intent(in):: nx,ny,mype integer(4),intent(out):: itime(6) real(4),dimension(is:ie,js:je),intent(out):: & pinc,tinc,qinc,uinc,vinc,tdinc !Declare local variables integer(4) nlon,nlat,nsig integer(4) i,j real(4),allocatable,dimension(:,:)::field1,field2 real(4),allocatable,dimension(:,:)::t1,p1,q1,td1 real(4),allocatable,dimension(:,:)::t2,p2,q2,td2 !-------------------------------------------------------------------- if (mype==0) print*,'in fetch_anlincs: nx,ny=',nx,ny allocate(field1(nx,ny)) allocate(field2(nx,ny)) allocate ( t1 (is:ie,js:je) ) allocate ( p1 (is:ie,js:je) ) allocate ( q1 (is:ie,js:je) ) allocate ( td1 (is:ie,js:je) ) allocate ( t2 (is:ie,js:je) ) allocate ( p2 (is:ie,js:je) ) allocate ( q2 (is:ie,js:je) ) allocate ( td2 (is:ie,js:je) ) open (52,file='sigges',form='unformatted') open (53,file='siganl',form='unformatted') read(52) itime,nlon,nlat,nsig read(52) field1,field2! glat,dx read(52) field1,field2! glon,dy print*,'in fetch_anlincs / from first guess' print*,'in fetch_anlincs / itime=',itime print*,'in fetch_anlincs / nlon,nlat,nsig=',nlon,nlat,nsig print*,'**********************************************' read(53) itime,nlon,nlat,nsig read(53) field1,field2! glat,dx read(53) field1,field2! glon,dy print*,'in fetch_anlincs / from analysis' print*,'in fetch_anlincs / itime=',itime print*,'in fetch_anlincs / nlon,nlat,nsig=',nlon,nlat,nsig print*,'**********************************************' read(52) field1 !psfc read(53) field2 !psfc do j=js,je do i=is,ie pinc(i,j)=field2(i,j)-field1(i,j) p1(i,j)=field1(i,j) p2(i,j)=field2(i,j) enddo enddo read(52) field1 !fis read(53) field2 !fis read(52) field1 !t read(53) field2 !t do j=js,je do i=is,ie tinc(i,j)=field2(i,j)-field1(i,j) t1(i,j)=field1(i,j) t2(i,j)=field2(i,j) enddo enddo read(52) field1 !q read(53) field2 !q do j=js,je do i=is,ie qinc(i,j)=field2(i,j)-field1(i,j) q1(i,j)=field1(i,j) q2(i,j)=field2(i,j) enddo enddo read(52) field1 !u read(53) field2 !u do j=js,je do i=is,ie uinc(i,j)=field2(i,j)-field1(i,j) enddo enddo read(52) field1 !v read(53) field2 !v do j=js,je do i=is,ie vinc(i,j)=field2(i,j)-field1(i,j) enddo enddo call get_dewpt(p1,q1,t1,td1,is,ie,js,je) call get_dewpt(p2,q2,t2,td2,is,ie,js,je) !dew point increment do j=js,je do i=is,ie tdinc(i,j) = max ( -20., td2(i,j)-td1(i,j) ) enddo enddo close(52) close(53) deallocate(field1) deallocate(field2) deallocate(t1) deallocate(p1) deallocate(q1) deallocate(td1) deallocate(t2) deallocate(p2) deallocate(q2) deallocate(td2) end subroutine fetch_anlincs !------------------------------------------------------ !------------------------------------------------------ subroutine bilinear_2d0v2(rffcst,ix1,ix2,jx1,jx2,rfobs,xx,yy) !$$$ subprogram documentation block ! . . . . ! subprogram: bilinear_2d0v2 ! prgmmr: ! ! abstract: ! ! ! input argument list: ! rffcst - model grid value ! ix,jx ! xx,yy - define coordinates in grid units ! of point for which interpolation is ! performed ! ! output argument list: ! rfobs - interpolated value ! ! notes: ! ! i+1,j | | i+1,j+1 ! --+----------+--- ! | | dym ! | * + - ! | x,y | dy ! | | ! --+----+-----+--- ! i,j||| i,j+1 ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block implicit none !declare passed variables integer(4),intent(in ) :: ix1,ix2,jx1,jx2 real(4) ,intent(in ) :: rffcst(ix1:ix2,jx1:jx2) real(4) ,intent(in ) :: xx,yy real(4) ,intent( out) :: rfobs !declare local variables integer(4) i,j,ip,jp real(4) dx,dy,dxm,dym i = ifix(yy) j = ifix(xx) dx = xx - float(j) dy = yy - float(i) dxm= 1.0-dx dym= 1.0-dy i=min(max(ix1,i),ix2) ; j=min(max(jx1,j),jx2) ip=min(ix2,i+1) ; jp=min(jx2,j+1) rfobs=dxm*(dym*rffcst(i,j)+dy*rffcst(ip,j)) & + dx *(dym*rffcst(i,jp)+dy*rffcst(ip,jp)) return end subroutine bilinear_2d0v2 !------------------------------------------------------ !------------------------------------------------------ subroutine anl_quality(rjbuffer_km,ds0,lobjanl) !******************************************************************** ! abstract: use barnes or cressman analysis to obtain gridded field * ! of analysis uncertainty * ! * ! program history log: * ! 2005-10-08 pondeca * ! * !******************************************************************** use mpi use kinds, only: i_kind,r_single use cressanl_common, only: nobsmax, xlocs, ylocs, hgt0s, hgts, & hobs, rmuses, oberrs, & dtimes, cstations, & obstypes, dups, bckgs, xberrs, & terrain, jpointer, & kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, & ktds,kgusts,kvis,kpblhs,kdists use cressanl_common, only: kflds !perform cressman analysis for 9 fields !(u,v,t,ps,q,gust,vis,pblh,dist) use controlvars, only: nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist use mpitaskmod, only: mype,npe use mpitaskmod, only: ista,iend,jsta,jend implicit none include 'param.incl' ! Declare passed variables real(r_single),intent(in):: rjbuffer_km real(r_single),intent(in):: ds0 logical,intent(in)::lobjanl ! Declare local parameters real(r_single),parameter :: gravity=9.81 real(r_single),parameter::spval=-99999. character(10),parameter::chstar10='**********' character(1),parameter::cblank1=' ' ! Declare local variables character(60) cgrid integer(i_kind) npass integer(i_kind) nx,ny integer(i_kind) i,j,m,n,nn,nt,np,k integer(i_kind) mlen,nlen integer(i_kind) lun1 integer(i_kind) ierror,itype,jflag integer(i_kind) rtime(6),nlon,nlat,nsig integer(i_kind) ii,jj,nn0 real(r_single) rlat,rlon,xx,yy,oberr,ob,ob_model,rinov real(r_single) qtflg,dtime,shgt0,hgt0,hgt,slm,rmuse real(r_single) rinovmax(kflds),radii(kflds) real(r_single) ds real(r_single) bb !dummy variable real(r_single) w0 real(r_single) radius_u,rinovmax_u,radius_v,rinovmax_v, & radius_t,rinovmax_t,radius_p,rinovmax_p, & radius_q,rinovmax_q, & radius_gust,rinovmax_gust,radius_vis,rinovmax_vis, & radius_pblh,rinovmax_pblh,radius_dist,rinovmax_dist real(r_single),allocatable,dimension(:,:):: & glon,glat,dx,dy,aux real(r_single),allocatable,dimension(:,:)::auxfield integer(i_kind),allocatable,dimension(:,:)::iauxfield character(8) cstation character(8) cprovider,csubprovider character(60) tdfname logical fexist1 logical lpadjust logical usebckg character(2) clun1 character(60) cvmodel(20) !names of verification models real(r_single) rmusecv,rmuseb integer(i_kind) itotmodel,itotrmuse namelist/cress_barnes_anlqlty/cgrid,npass,& radius_u,rinovmax_u,radius_v,rinovmax_v, & radius_t,rinovmax_t,radius_p,rinovmax_p, & radius_q,rinovmax_q, & radius_gust,rinovmax_gust,radius_vis,rinovmax_vis, & radius_pblh,rinovmax_pblh,radius_dist,rinovmax_dist, & lpadjust,usebckg, & cvmodel,rmusecv,rmuseb !===================================================================== ! MPI setup ! call mpi_init(ierror) ! call mpi_comm_size(mpi_comm_world,npe,ierror) ! call mpi_comm_rank(mpi_comm_world,mype,ierror) data cgrid/'conus'/ data npass/1/ data radius_u/18./ ; data rinovmax_u/10./ data radius_v/18./ ; data rinovmax_v/10./ data radius_t/12./ ; data rinovmax_t/7.5/ data radius_p/11./ ; data rinovmax_p/1000./ data radius_q/12./ ; data rinovmax_q/5.e-03/ data radius_gust/18./ ; data rinovmax_gust/10./ data radius_vis /18./ ; data rinovmax_vis/9000./ data radius_pblh/18./ ; data rinovmax_pblh/1000./ data radius_dist/18./ ; data rinovmax_dist/1000./ data lpadjust/.false./ data usebckg/.false./ ! !===================================================================== data rmuseb /+1./ do n=1,20 cvmodel(n)=chstar10 enddo open (55,file='cress_barnesparm.anl',form='formatted') read(55,cress_barnes_anlqlty) close(55) if (mype==0) then print*,'========in anl_quality ===========' print*,'cgrid=',trim(cgrid) print*,'npass=',npass print*,'radius_u=',radius_u print*,'radius_v=',radius_v print*,'radius_t=',radius_t print*,'radius_p=',radius_p print*,'radius_q=',radius_q if (nrf2_gust > 0) print*,'radius_gust =',radius_gust if (nrf2_vis > 0) print*,'radius_vis =',radius_vis if (nrf2_pblh > 0) print*,'radius_pblh =',radius_pblh if (nrf2_dist > 0) print*,'radius_dist =',radius_dist print*,'in anl_quality:' print*,'rinovmax_u=',rinovmax_u print*,'rinovmax_v=',rinovmax_v print*,'rinovmax_t=',rinovmax_t print*,'rinovmax_p=',rinovmax_p print*,'rinovmax_q=',rinovmax_q if (nrf2_gust > 0) print*,'rinovmax_gust =',rinovmax_gust if (nrf2_vis > 0) print*,'rinovmax_vis =',rinovmax_vis if (nrf2_pblh > 0) print*,'rinovmax_pblh =',rinovmax_pblh if (nrf2_dist > 0) print*,'rinovmax_dist =',rinovmax_dist print*,'rmusecv=',rmusecv print*,'rmuseb=',rmuseb print*,'lpadjust=',lpadjust print*,'usebckg=',usebckg endif radii(1)=radius_u radii(2)=radius_v radii(3)=radius_t radii(4)=radius_p radii(5)=radius_q radii(6)=radius_gust radii(7)=radius_vis radii(8)=radius_pblh radii(9)=radius_dist if (usebckg) then rinovmax(1)=rinovmax_u rinovmax(2)=rinovmax_v rinovmax(3)=rinovmax_t rinovmax(4)=rinovmax_p rinovmax(5)=rinovmax_q rinovmax(6)=rinovmax_gust rinovmax(7)=rinovmax_vis rinovmax(8)=rinovmax_pblh rinovmax(9)=rinovmax_dist else if (mype==0) print*,'in anl_quality: reset rinovmax and lpadjust & & to reflect non-use of first guess fields' ! rinovmax(:)=1.e+20 rinovmax(1)=rinovmax_u rinovmax(2)=rinovmax_v rinovmax(3)=rinovmax_t rinovmax(4)=1.e+20! rinovmax_p rinovmax(5)=rinovmax_q rinovmax(6)=rinovmax_gust rinovmax(7)=rinovmax_vis rinovmax(8)=rinovmax_pblh rinovmax(9)=rinovmax_dist lpadjust=.false. if (mype==0) print*,'in anl_quality: lpadjust reset to: ',lpadjust endif call domain_dims(cgrid,nx,ny,ds) if (mype==0) print*,'in anl_quality: nx,ny=', & nx,ny allocate(aux(nx,ny)) allocate(auxfield(nx,ny)) allocate(iauxfield(nx,ny)) ! !===================================================================== !==>all tasks read in the first guess and the bckg error variances !===================================================================== open (53,file='sigges',form='unformatted') read(53) rtime,nlon,nlat,nsig read(53) auxfield, aux !glat,dx read(53) auxfield, aux !glon,dy read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,4)=auxfield(i,j) !psfc1 enddo enddo read(53) auxfield terrain=auxfield/gravity !no need for subdomain distribution read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,3)=auxfield(i,j) !t1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,5)=auxfield(i,j) !q1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,1)=auxfield(i,j) !u1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,2)=auxfield(i,j) !v1 enddo enddo do n=1,11 !must jump 11 records to get to gust read(53) enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,6)=auxfield(i,j) !gust1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,7)=auxfield(i,j) !vis1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,8)=auxfield(i,j) !pblh1 enddo enddo read(53) auxfield do j=jsta,jend do i=ista,iend bckgs(i,j,9)=auxfield(i,j) !dist enddo enddo close(53) xberrs=0. do nn=1,kflds if (nn==6 .and. nrf2_gust <= 0) cycle if (nn==7 .and. nrf2_vis <= 0) cycle if (nn==8 .and. nrf2_pblh <= 0) cycle if (nn==9 .and. nrf2_dist <= 0) cycle if (nn==1) open (94,file='bckgvar.dat_u', form='unformatted') if (nn==2) open (94,file='bckgvar.dat_v', form='unformatted') if (nn==3) open (94,file='bckgvar.dat_t', form='unformatted') if (nn==4) open (94,file='bckgvar.dat_ps',form='unformatted') if (nn==5) open (94,file='bckgvar.dat_q', form='unformatted') if (nn==6) open (94,file='bckgvar.dat_gust',form='unformatted') if (nn==7) open (94,file='bckgvar.dat_vis' ,form='unformatted') if (nn==8) open (94,file='bckgvar.dat_pblh',form='unformatted') if (nn==9) open (94,file='bckgvar.dat_dist',form='unformatted') read(94) auxfield do j=jsta,jend do i=ista,iend xberrs(i,j,nn)=auxfield(i,j) enddo enddo close(94) if (mype==0) print*,'in anl_quality: mype,nn,xberrs,min,max=', & mype,nn,minval(xberrs(:,:,nn)),maxval(xberrs(:,:,nn)) enddo ! !===================================================================== !==>load array dups, which is used to remove duplicate obs in ! the cressman analysis !===================================================================== dups=.false. do nn=1,kflds if (nn==1) nn0=kugrds if (nn==2) nn0=kvgrds if (nn==3) nn0=kts if (nn==4) nn0=kps if (nn==5) nn0=kqs if (nn==6) nn0=kgusts if (nn==7) nn0=kvis if (nn==8) nn0=kpblhs if (nn==9) nn0=kdists do jj=1,nn0 j=jpointer(jj,nn) nlen=0 do k=1,8 if (cstations(j)(k:k) /= cblank1) then nlen=nlen+1 else exit endif enddo do ii=jj+1,nn0 i=jpointer(ii,nn) mlen=0 do k=1,8 if (cstations(i)(k:k) /= cblank1) then mlen=mlen+1 else exit endif enddo if ( mlen==nlen ) then if (cstations(i)(1:mlen)==cstations(j)(1:nlen)) then if (abs(dtimes(i)) >= abs(dtimes(j))) dups(i)=.true. if (abs(dtimes(i)) < abs(dtimes(j))) dups(j)=.true. endif endif enddo enddo enddo !===================================================================== !==> perform cressman (or barnes) analysis !===================================================================== if (lobjanl) & call obj_anl4(rinovmax,lpadjust,usebckg,radii,spval, & nx,ny,kflds,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist, & npass,ista,iend,jsta,jend,rjbuffer_km,ds0,npe,mype) if (mype==0) then open (54,file='sigcress',form='unformatted') ; rewind (54) open (53,file='sigges',form='unformatted') ; rewind (53) read(53) rtime,nlon,nlat,nsig write(54) rtime,nlon,nlat,nsig read(53) auxfield, aux !glat,dx write(54) auxfield, aux read(53) auxfield, aux !glon,dy write(54) auxfield, aux endif call mpi_barrier(mpi_comm_world,ierror) do n=1,6 if (n==1) nn=4 !pfsc if (n==3) nn=3 !temp if (n==4) nn=5 !moisture if (n==5) nn=1 !uwind if (n==6) nn=2 !vwind if (n==2) then auxfield=terrain*gravity else aux(:,:)=0. do j=jsta,jend do i=ista,iend aux(i,j)=bckgs(i,j,nn) enddo enddo call mpi_allreduce(aux,auxfield,nx*ny, & mpi_real4,mpi_sum,mpi_comm_world,ierror) endif if (mype==0) then write(54) auxfield read(53) aux !advance sigges endif enddo call mpi_barrier(mpi_comm_world,ierror) !(see subroutine wr2d_binary) if (mype==0) then auxfield=0. ! mimics increment of pseudo RH write(54) auxfield do nn=7,17 ! SM,SICE,SST,IVGTYP,ISLTYP,VEGFRA,SNOW,SMOIS,TSLB,TSK,SFCR if (nn==10 .or. nn==11) then read(53) iauxfield write(54) iauxfield else read(53) auxfield write(54) auxfield endif enddo endif call mpi_barrier(mpi_comm_world,ierror) aux(:,:)=0. do j=jsta,jend do i=ista,iend w0=bckgs(i,j,1)*bckgs(i,j,1) + & bckgs(i,j,2)*bckgs(i,j,2) if (w0 > 0.) w0=sqrt(w0) if ( bckgs(i,j,6) < w0) then aux(i,j)=w0 else aux(i,j)=bckgs(i,j,6) endif enddo enddo call mpi_allreduce(aux,auxfield,nx*ny, & !gust mpi_real4,mpi_sum,mpi_comm_world,ierror) if (mype==0) then write(54) auxfield !gust read(53) aux!advance sigges endif call mpi_barrier(mpi_comm_world,ierror) do nn=7,9 !vis,pblh,cldch aux(:,:)=0. do j=jsta,jend do i=ista,iend aux(i,j)=bckgs(i,j,nn) enddo enddo call mpi_allreduce(aux,auxfield,nx*ny, & mpi_real4,mpi_sum,mpi_comm_world,ierror) if (mype==0) then write(54) auxfield read(53) aux !advance sigges endif enddo call mpi_barrier(mpi_comm_world,ierror) !copy from sigges to sigcress if (mype==0) then do n=1,10 !wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas,uwnd10m,vwnd10m read(53) aux write(54) aux enddo close(53) close(54) endif call mpi_barrier(mpi_comm_world,ierror) ! PRINT*,'IN POST: CRESSMAN ANL DONE!' ! CALL MPI_FINALIZE(IERROR) ! STOP tdfname='sigcress' call td_flds(nlon,nlat,tdfname,mype,npe) !==> verification ! Determine if there is valid verf info: ! the assumption is that all valid models are listed before the first ! non-verification model (named chstar10) do n=1,20 if (trim(cvmodel(n)) .ne. chstar10) then if (mype==0) print*,' requested verf using model ',trim(cvmodel(n)) else exit endif enddo itotmodel=n-1 if (mype==0) print*,'itotmodel=',itotmodel if (itotmodel > 0) then do n=1,itotmodel call bias_rmse_cv_mpi(cgrid,itotmodel,n,cvmodel(n),rmusecv,rmuseb, & lpadjust,cblank1) enddo endif deallocate(aux) deallocate(auxfield) deallocate(iauxfield) end subroutine anl_quality !------------------------------------------------------ subroutine obj_anl4(rinovmax,lpadjust,usebckg,radii,spval, & nx,ny,kflds,igust,ivis,ipblh,idist, & npass,ista,iend,jsta,jend,rjbuffer_km,ds0,npe,mype) use mpi use kinds, only: i_kind,r_single use cressanl_common, only: nobsmax, xlocs, ylocs, & hobs, rmuses, oberrs, & obstypes, bckgs, xberrs, terrain implicit none include 'param.incl' !==>Declare passed varaibles integer(i_kind),intent(in) :: igust,ivis,ipblh,idist integer(i_kind),intent(in) :: npe,mype integer(i_kind),intent(in) :: ista,iend,jsta,jend integer(i_kind),intent(in) :: nx,ny integer(i_kind),intent(in) :: kflds,npass real(r_single),intent(in) :: rjbuffer_km real(r_single),intent(in) :: ds0 real(r_single),intent(in) :: rinovmax(kflds) logical,intent(in) :: lpadjust logical,intent(in) :: usebckg real(r_single),intent(in) :: radii(kflds) real(r_single),intent(in) :: spval !==>Declare local parameters integer(i_kind),parameter:: ndomains=10 !npe !ny/(20*nint(radii(nn))) !==>Declare local variables integer(i_kind) mype1 integer(i_kind) jbuffer integer(i_kind) ierror integer(i_kind) iter,nn,ij,npts,npts2 integer(i_kind) i,j,k,n,m,im0 integer(i_kind) nobssubmax integer(i_kind) nclose integer(i_kind),allocatable,dimension(:,:):: nobssub integer(i_kind),allocatable,dimension(:,:,:):: npointer real(r_single),allocatable,dimension(:,:):: g0,g real(r_single),allocatable,dimension(:):: fincrs,qcflgs integer(i_kind),allocatable,dimension(:):: j1,j2 integer(i_kind) nobsmax0,nobssubmax0,iproc integer(i_kind) ista0,iend0,jsta0,jend0 integer(i_kind),allocatable,dimension(:,:):: nobssub0 integer(i_kind),allocatable,dimension(:,:,:):: npointer0 real(r_single),allocatable,dimension(:):: xlocs0,ylocs0,fincrs0,qcflgs0 real(r_single) sumprod,wsum real(r_single) ylower,yupper real(r_single) wij real(r_single) xob,yob,elev real(r_single) xi,xj,weight4,dist2,radius real(r_single) xberr0,oberr0 real(r_single) u00,v00,t00,p00,q00,gust00,vis00,pblh00,dist00 !****************************************************************** !****************************************************************** allocate(fincrs(nobsmax)) allocate(qcflgs(nobsmax)) allocate(nobssub(ndomains,kflds)) allocate(j1(ndomains)) allocate(j2(ndomains)) mype1=mype+1 jbuffer=1000._r_single*rjbuffer_km/ds0 if (mype==0) print*,'in obj_anl4: rjbuffer_km,ds0,jbuffer=', & rjbuffer_km,ds0,jbuffer j1=-999 ; j2=-999 call para_range3(1,ny,ndomains,j1(1:ndomains),j2(1:ndomains)) nobssub(:,:)=0 do n=1,nobsmax if (trim(obstypes(n))=='ugrel-ob') then ; nn=1 elseif (trim(obstypes(n))=='vgrel-ob') then ; nn=2 elseif (trim(obstypes(n))=='t-ob') then ; nn=3 elseif (trim(obstypes(n))=='p-ob') then ; nn=4 elseif (trim(obstypes(n))=='q-ob') then ; nn=5 elseif (trim(obstypes(n))=='gust-ob') then ; nn=6 elseif (trim(obstypes(n))=='vis-ob') then ; nn=7 elseif (trim(obstypes(n))=='pblh-ob') then ; nn=8 elseif (trim(obstypes(n))=='dist-ob') then ; nn=9 else ; cycle ; endif do m=1,ndomains ylower = float(max(1,(j1(m)-jbuffer))) yupper = float(min(ny,(j2(m)+jbuffer))) if (ylocs(n) .ge. ylower .and. ylocs(n) .le. yupper) then nobssub(m,nn)=nobssub(m,nn)+1 endif enddo enddo do nn=1,kflds do m=1,ndomains print*,'in obj_anl4: mype,nn,m,nobssub(m,nn)=',mype,nn,m,nobssub(m,nn) enddo enddo nobssubmax=maxval(nobssub(:,:)) print*,'in obj_anl4: mype,nobssubmax=',nobssubmax allocate(npointer(nobssubmax,ndomains,kflds)) nobssub(:,:)=0 npointer(:,:,:)=0 do n=1,nobsmax if (trim(obstypes(n))=='ugrel-ob') then ; nn=1 elseif (trim(obstypes(n))=='vgrel-ob') then ; nn=2 elseif (trim(obstypes(n))=='t-ob') then ; nn=3 elseif (trim(obstypes(n))=='p-ob') then ; nn=4 elseif (trim(obstypes(n))=='q-ob') then ; nn=5 elseif (trim(obstypes(n))=='gust-ob') then ; nn=6 elseif (trim(obstypes(n))=='vis-ob') then ; nn=7 elseif (trim(obstypes(n))=='pblh-ob') then ; nn=8 elseif (trim(obstypes(n))=='dist-ob') then ; nn=9 else ; cycle ; endif do m=1,ndomains ylower = float(max(1,(j1(m)-jbuffer))) yupper = float(min(ny,(j2(m)+jbuffer))) if (ylocs(n) .ge. ylower .and. ylocs(n) .le. yupper) then nobssub(m,nn)=nobssub(m,nn)+1 im0=nobssub(m,nn) npointer(im0,m,nn)=n endif enddo enddo if (.not.usebckg) then ! bckgs=0. !original option call avghobs(mype,npe,u00,v00,t00,p00,q00, & gust00,vis00,pblh00,dist00) bckgs(:,:,1)=u00 bckgs(:,:,2)=v00 bckgs(:,:,3)=t00 bckgs(:,:,4)=p00 bckgs(:,:,5)=q00 bckgs(:,:,6)=gust00 bckgs(:,:,7)=vis00 bckgs(:,:,8)=pblh00 bckgs(:,:,9)=dist00 if (mype==0) print*,'in obj_anl4: u00,v00,t00,p00,q00=',& u00,v00,t00,p00,q00 if (mype==0) print*,'in obj_anl4: gust00,vis00,pblh00,dist00=',& gust00,vis00,pblh00,dist00 endif do iter=1,npass call computeincrs(rinovmax,lpadjust, & fincrs,qcflgs,spval,igust,ivis,ipblh,idist, & ista,iend,jsta,jend,nx,ny,kflds,nobsmax,mype,npe) do iproc=1,npe if (mype1==iproc) then nobsmax0=nobsmax nobssubmax0=nobssubmax ista0=ista iend0=iend jsta0=jsta jend0=jend endif call mpi_bcast (nobsmax0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) call mpi_bcast (nobssubmax0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) call mpi_bcast (ista0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) call mpi_bcast (iend0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) call mpi_bcast (jsta0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) call mpi_bcast (jend0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror) allocate(npointer0(nobssubmax0,ndomains,kflds)) allocate(nobssub0(ndomains,kflds)) allocate(xlocs0(nobsmax0)) allocate(ylocs0(nobsmax0)) allocate(fincrs0(nobsmax0)) allocate(qcflgs0(nobsmax0)) if (mype1==iproc) then npointer0(:,:,:)=npointer(:,:,:) nobssub0(:,:)=nobssub(:,:) xlocs0(:)=xlocs(:) ylocs0(:)=ylocs(:) fincrs0(:)=fincrs(:) qcflgs0(:)=qcflgs(:) endif call mpi_barrier(mpi_comm_world,ierror) npts=nobssubmax0*ndomains*kflds npts2=ndomains*kflds call mpi_bcast (npointer0, npts, mpi_integer, iproc-1, mpi_comm_world, ierror) call mpi_bcast (nobssub0, npts2, mpi_integer, iproc-1, mpi_comm_world, ierror) call mpi_bcast (xlocs0, nobsmax0, mpi_real, iproc-1, mpi_comm_world, ierror) call mpi_bcast (ylocs0, nobsmax0, mpi_real, iproc-1, mpi_comm_world, ierror) call mpi_bcast (fincrs0, nobsmax0, mpi_real, iproc-1, mpi_comm_world, ierror) call mpi_bcast (qcflgs0, nobsmax0, mpi_real, iproc-1, mpi_comm_world, ierror) allocate(g (ista0:iend0,jsta0:jend0)) allocate(g0(ista0:iend0,jsta0:jend0)) do nn=1,kflds radius=radii(nn) g0=0. ij=0 do m=1,ndomains do j=j1(m), j2(m) !these are the same on all processors xj=float(j) do i=ista0, iend0 !note these. they are the dims on the iproc processor xi=float(i) ij=ij+1 if (mype==mod(ij-1,npe)) then sumprod=0. wsum=0. nclose=0 do k=1,nobssub0(m,nn) n=npointer0(k,m,nn) if (qcflgs0(n) .lt. 1.) cycle xob=xlocs0(n) yob=ylocs0(n) dist2=(xi-xob)*(xi-xob)+(xj-yob)*(xj-yob) if (dist2 .le. (2.*radius)**2 ) then nclose=nclose+1 wij=weight4(xi,xj,xob,yob,radius) sumprod=sumprod+wij*fincrs0(n) wsum=wsum+wij endif enddo !k-loop !!!if (wsum.gt.1.e-5) then !!! g0(i,j)=sumprod/wsum !!!endif if (nclose .ge. 10) then !6 g0(i,j)=sumprod/wsum endif endif !mype-condition enddo !i-loop enddo !j-loop enddo !m-loop npts=(iend0-ista0+1)*(jend0-jsta0+1) call mpi_allreduce(g0(ista0:iend0,jsta0:jend0),g(ista0:iend0,jsta0:jend0), & npts,mpi_real4,mpi_sum,mpi_comm_world,ierror) if (mype1==iproc) then do j=jsta0,jend0 do i=ista0,iend0 bckgs(i,j,nn)=bckgs(i,j,nn)+g(i,j) enddo enddo endif call mpi_barrier(mpi_comm_world,ierror) if (mype==0) then print*,'in obj_anl4: iter,iproc,nn,g,min,max=',iter,iproc,nn, & minval(g(ista0:iend0,jsta0:jend0)), & maxval(g(ista0:iend0,jsta0:jend0)) endif call mpi_barrier(mpi_comm_world,ierror) enddo !nn-loop deallocate(npointer0) deallocate(nobssub0) deallocate(xlocs0) deallocate(ylocs0) deallocate(fincrs0) deallocate(qcflgs0) deallocate(g) deallocate(g0) enddo !iproc-loop enddo !iter-loop deallocate(fincrs) deallocate(qcflgs) deallocate(nobssub) deallocate(npointer) deallocate(j1) deallocate(j2) return end function weight4(xi,xj,xob,yob,r) use kinds, only: r_single implicit none real(r_single) xi,xj,xob,yob,r,radius2,dist2,weight4 radius2=r*r dist2=(xi-xob)*(xi-xob)+(xj-yob)*(xj-yob) ! weight4=max(0.,(radius2-dist2)/(radius2+dist2)) weight4=exp(-dist2/(2.*radius2)) ! weight4=exp(-4.*dist2/radius2) !Barnes type return end !------------------------------------------------------ subroutine avghobs(mype,npe,u00,v00,t00,p00,q00, & gust00,vis00,pblh00,dist00) use mpi use kinds, only: i_kind,r_single use cressanl_common, only: nobsmax, hobs, rmuses, & obstypes implicit none !==>Declare passed variables integer(i_kind),intent(in) :: mype,npe real(r_single),intent(out) :: u00,v00,t00,p00,q00, & gust00,vis00,pblh00,dist00 !Declare local parameters integer(i_kind),parameter::nt=9 !u,v,t,p,q,gust,vis,pblh,dist !Declare local variables integer(i_kind) n,nn,ierror integer(i_kind),allocatable::nokays(:),nokays2(:) real(r_single),allocatable::h00(:),h002(:) allocate(nokays(nt),nokays2(nt)) allocate(h00(nt),h002(nt)) u00=0. ; v00=0. ; t00=0. ; p00=0. ; q00=0. gust00=0. ; vis00=0. ; pblh00=0. ; dist00=0. h00(:)=0. nokays(:)=0 do n=1,nobsmax if (abs(rmuses(n)-1.).ge.1.e-03) cycle if (trim(obstypes(n))=='ugrel-ob') then ; nn=1 elseif (trim(obstypes(n))=='vgrel-ob') then ; nn=2 elseif (trim(obstypes(n))=='t-ob') then ; nn=3 elseif (trim(obstypes(n))=='p-ob') then ; nn=4 elseif (trim(obstypes(n))=='q-ob') then ; nn=5 elseif (trim(obstypes(n))=='gust-ob') then ; nn=6 elseif (trim(obstypes(n))=='vis-ob') then ; nn=7 elseif (trim(obstypes(n))=='pblh-ob') then ; nn=8 elseif (trim(obstypes(n))=='dist-ob') then ; nn=9 else ; cycle ; endif h00(nn)=h00(nn)+hobs(n) nokays(nn)=nokays(nn)+1 enddo h002(:)=h00 nokays2(:)=nokays(:) call mpi_allreduce(h002, h00, nt ,mpi_real4, mpi_sum, mpi_comm_world,ierror) call mpi_allreduce(nokays2, nokays, nt ,mpi_integer, mpi_sum, mpi_comm_world,ierror) if (mype==0) then do nn=1,nt print*,'nn,nokays(nn)=',nn,nokays(nn) enddo endif if (nokays(1) > 0 ) u00=h00(1)/real(nokays(1),r_single) if (nokays(2) > 0 ) v00=h00(2)/real(nokays(2),r_single) if (nokays(3) > 0 ) t00=h00(3)/real(nokays(3),r_single) if (nokays(4) > 0 ) p00=h00(4)/real(nokays(4),r_single) if (nokays(5) > 0 ) q00=h00(5)/real(nokays(5),r_single) if (nokays(6) > 0) gust00=h00(6)/real(nokays(6),r_single) if (nokays(7) > 0) vis00= h00(7)/real(nokays(7),r_single) if (nokays(8) > 0) pblh00=h00(8)/real(nokays(8),r_single) if (nokays(9) > 0) dist00=h00(9)/real(nokays(9),r_single) deallocate(nokays,nokays2) deallocate(h00,h002) end subroutine avghobs !------------------------------------------------------ subroutine computeincrs(rinovmax,lpadjust, & fincrs,qcflgs,spval,igust,ivis,ipblh,idist, & ista,iend,jsta,jend,nx,ny,kflds,nobsmax,mype,npe) use mpi use kinds, only: i_kind,r_single use cressanl_common, only: xlocs, ylocs, hgt0s, & hobs, rmuses, oberrs, & obstypes, insubdoms, dups, bckgs, & terrain, jpointer, & kts,kqs,kps,kugrds,kvgrds,kgusts, & kvis,kpblhs,kdists implicit none !==>Declare passed variables integer(i_kind),intent(in) :: igust,ivis,ipblh,idist integer(i_kind),intent(in) :: ista,iend,jsta,jend integer(i_kind),intent(in) :: nx,ny,kflds integer(i_kind),intent(in) :: mype,npe integer(i_kind),intent(in) :: nobsmax real(r_single),intent(in) :: rinovmax(kflds) logical,intent(in) :: lpadjust real(r_single),intent(in) :: spval real(r_single),intent(out) :: qcflgs(nobsmax) real(r_single),intent(out) :: fincrs(nobsmax) !Declare local variables integer(i_kind) n,nn,nn0,i,j,ii integer(i_kind) ierror,mype1 integer(i_kind) mts,mqs,mps,mus,mvs,mgusts,mvis,mpblhs,mdists integer(i_kind) mts2,mqs2,mps2,mus2,mvs2,mgusts2,mvis2,mpblhs2,mdists2 real(r_single),allocatable,dimension(:,:):: field,field2 real(r_single),allocatable,dimension(:,:):: bckgtv,tfield real(r_single) xx,yy,rmuse,fint real(r_single) tges,tob,qob,pges,zsges,hgt0 real(r_single) rdelz,half,half_tlapse,g_over_rd,rdp,rtvts half=0.5 half_tlapse=0.00325 g_over_rd=9.81/287.04 allocate(bckgtv(nx,ny)) allocate(tfield(nx,ny)) allocate(field(nx,ny)) allocate(field2(nx,ny)) field (:,:)=0. field2(:,:)=0. do j=jsta,jend do i=ista,iend field (i,j)=bckgs(i,j,3) ! temperature field2(i,j)=bckgs(i,j,3)*(1.+0.608*max(0.,bckgs(i,j,5))) !virtual temperature enddo enddo call mpi_allreduce(field,tfield,nx*ny, & mpi_real4,mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(field2,bckgtv,nx*ny, & mpi_real4,mpi_sum,mpi_comm_world,ierror) if (mype==0) then print*,'bckg-t, min,max=',minval(tfield),maxval(tfield) print*,'bckg-tv, min,max=',minval(bckgtv),maxval(bckgtv) endif mype1=mype+1 qcflgs(:)=0. fincrs(:)=0. mts=0 ; mqs=0 ; mps=0 ; mus=0 ; mvs=0 ; mgusts=0 ; mvis=0 mpblhs=0 ; mdists=0 do nn=1,kflds if (nn==1) nn0=kugrds if (nn==2) nn0=kvgrds if (nn==3) nn0=kts if (nn==4) nn0=kps if (nn==5) nn0=kqs if (nn==6) nn0=kgusts if (nn==7) nn0=kvis if (nn==8) nn0=kpblhs if (nn==9) nn0=kdists field2(:,:)=0. do j=jsta,jend do i=ista,iend field2(i,j)=bckgs(i,j,nn) enddo enddo call mpi_allreduce(field2,field,nx*ny, & mpi_real4,mpi_sum,mpi_comm_world,ierror) do ii=1,nn0 n=jpointer(ii,nn) xx=xlocs(n) yy=ylocs(n) rmuse=rmuses(n) if (abs(rmuse-1.).gt.1.e-03 .or. dups(n)) cycle call bilinear_2d0v2(field,1,nx,1,ny,fint,yy,xx)!Note the reverse order "yy,xx" if (nn==4 .and. lpadjust) then pges=fint call bilinear_2d0v2(terrain,1,nx,1,ny,zsges,yy,xx) rdp=0. rtvts=1. !hardwired to be sensible temperature for now if (rtvts==0.) then call bilinear_2d0v2(bckgtv,1,nx,1,ny,tges,yy,xx) else call bilinear_2d0v2(tfield,1,nx,1,ny,tges,yy,xx) endif hgt0=hgt0s(n) tob=spval !ideally would use observed temperature rdelz=hgt0-zsges if (tob /= spval) then tges = half*(tges+tob) else if(rdelz < 0.)then tges=tges-half_tlapse*rdelz endif endif rdp = g_over_rd*rdelz/tges !Adjust hydrostatically fint=1000.*exp(log(pges/1000.) - rdp) !/100. !hPa endif ! if (nn==3) then ! rtvts=1. !hardwired to be sensible temperature for now ! if (rtvts==0.) then ! call bilinear_2d0v2(bckgtv,1,nx,1,ny,fint,yy,xx) ! endif ! endif fincrs(n)=hobs(n)-fint if (abs(fincrs(n)).le.rinovmax(nn)) then qcflgs(n)=+1. if (insubdoms(n)) then if (nn==1) mus=mus+1 if (nn==2) mvs=mvs+1 if (nn==3) mts=mts+1 if (nn==4) mps=mps+1 if (nn==5) mqs=mqs+1 if (nn==6) mgusts=mgusts+1 if (nn==7) mvis=mvis+1 if (nn==8) mpblhs=mpblhs+1 if (nn==9) mdists=mdists+1 endif endif enddo enddo call mpi_allreduce(mus, mus2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mvs, mvs2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mts, mts2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mps, mps2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mqs, mqs2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mgusts, mgusts2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mvis, mvis2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mpblhs, mpblhs2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) call mpi_allreduce(mdists, mdists2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror) if (mype==0) then print*,'in computeincrs: #of usable u-obs=',mus2 print*,'in computeincrs: #of usable v-obs=',mvs2 print*,'in computeincrs: #of usable t-obs=',mts2 print*,'in computeincrs: #of usable p-obs=',mps2 print*,'in computeincrs: #of usable q-obs=',mqs2 print*,'in computeincrs: #of usable guts-obs=',mgusts2 print*,'in computeincrs: #of usable vis-obs=',mvis2 print*,'in computeincrs: #of usable pblh-obs=',mpblhs2 print*,'in computeincrs: #of usable dist-obs=',mdists2 endif deallocate(bckgtv) deallocate(tfield) deallocate(field) deallocate(field2) return end !------------------------------------------------------ !------------------------------------------------------ subroutine cvbasedrecal(cgrid,cmodel0) use mpi use kinds, only: i_kind,r_single,r_kind use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, & wspderr,wdirerr,wdirerr2,terr,tderr, & qerr,perr,gusterr,viserr,pblherr,disterr, & wspd10merr,td2merr,mxtmerr,mitmerr,pmslerr, & howverr,tcamterr,lcbaserr, & uwnd10merr,vwnd10merr 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,iuwnd10m,ivwnd10m implicit none !declare passed variables character(60),intent(in)::cgrid character(60),intent(in)::cmodel0 !declare local parameters real(r_single),parameter::epsilon=1.e-06_r_single character(1),parameter::ccomma=',' integer(i_kind),parameter::itotmodel=2 logical,parameter::lpadjust=.true. !declare local variables real(4) rmusecv,rmuseb character(60) cmodel integer(i_kind) nflds character(60) fnamesuffix character(80) rawstatsflname integer(i_kind),allocatable:: nobs(:) real(r_single),allocatable:: bias(:),abse(:),rmse(:) integer(i_kind),allocatable:: nobs2(:) real(r_single),allocatable:: bias2(:),abse2(:),rmse2(:) real(r_single),allocatable:: field(:,:) real(r_single) xx,yy,fint real(r_single) rmse0,rmse1,ratio logical lres1,lres2,lres3 integer(i_kind) nobsmax0 integer(i_kind) n,nn logical fexist if (mype==0) print*,'in cvbasedrecal: cmodel0=',trim(cmodel0) if (trim(cmodel0)=='sigfupdate02') then !applicable when miter=2. cv-obs put back at start of 2nd outer loop fnamesuffix='_fg02only' !this is the gsi output at the end of first outer loop rawstatsflname='rawstats.dat_for_rmuse_neg2'//trim(fnamesuffix) rmusecv=-2. rmuseb=+1. else if (trim(cmodel0)=='sigfupdate03') then !applicable when miter=3. cv-obs put back at start of 3rd outer loop fnamesuffix='_fg03only' !this is the gsi output at the end of second outer loop rawstatsflname='rawstats.dat_for_rmuse_neg3'//trim(fnamesuffix) rmusecv=-3. rmuseb=+1. else fnamesuffix='_anlonly' !applicable when miter=2, and cv-obs aren't put back into the analysis rawstatsflname='rawstats.dat_for_rmuse_neg3'//trim(fnamesuffix) rmusecv=-3. rmuseb=+1. endif if (mype==0) then print*,'in cvbasedrecal: fnamesuffix,rmusecv,rmuseb=',trim(fnamesuffix),rmusecv,rmuseb print*,'in cvbasedrecal: rawstatsflname=',trim(rawstatsflname) endif do n=1,itotmodel if (n==1) cmodel=cmodel0 if (n==2) cmodel='errfield.dat_precval' call bias_rmse_cv_mpi(cgrid,itotmodel,n,cmodel,rmusecv,rmuseb, & lpadjust,fnamesuffix) enddo inquire(file=trim(rawstatsflname),exist=fexist) if (.not.fexist) then if (mype==0) then print*,'in cvbasedrecal: no input file available' print*,'apparently no cross-validation was performed to create input file' endif return endif open (18,file=trim(rawstatsflname),form='unformatted') read(18) cmodel !this is either sigfupdate02, sigfupdate03, or siganl read(18) nflds,nobsmax0 allocate(nobs(nflds)) allocate(bias(nflds)) allocate(abse(nflds)) allocate(rmse(nflds)) read(18) nobs,bias,abse,rmse if (mype==0) then print*,'----------------------------' print*,'in cvbasedrecal: first read' print*,'cmodel=',trim(cmodel) print*,'nflds,nobsmax0=',nflds,nobsmax0 do nn=1,nflds !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist !wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas,uwnd10m,vwnd10m print*,'nn,nobs,bias,abse,rmse=',& nn,nobs(nn),bias(nn),abse(nn),rmse(nn) enddo print*,'----------------------------' endif lres1=any(bias(1:7) > 9998.) !any(bias(1:nflds) > 9998.) lres2=any(abse(1:7) > 9998.) !any(abse(1:nflds) > 9998.) lres3=any(rmse(1:7) > 9998.) !any(rmse(1:nflds) > 9998.) if (lres1 .or. lres2 .or. lres3) then if (mype==0) then print*,'in cvbasedrecal/ first read : apparently no cross-validation was done' print*,' do nothing and return' endif return endif open (19,file='cverrorupdate_input_current', form='formatted') write(19,'(a14)') '&cverrorupdate' write(19,'(a17,f14.6,a1)') 'pcverr = ' ,rmse(1),ccomma write(19,'(a17,f14.6,a1)') 'tcverr = ' ,rmse(2),ccomma write(19,'(a17,f14.6,a1)') 'qcverr = ' ,rmse(3),ccomma write(19,'(a17,f14.6,a1)') 'ucverr = ' ,rmse(4),ccomma write(19,'(a17,f14.6,a1)') 'u2cverr = ' ,rmse(4),ccomma write(19,'(a17,f14.6,a1)') 'vcverr = ' ,rmse(5),ccomma write(19,'(a17,f14.6,a1)') 'v2cverr = ' ,rmse(5),ccomma write(19,'(a17,f14.6,a1)') 'wspdcverr = ' ,rmse(6),ccomma write(19,'(a17,f14.6,a1)') 'wdircverr = ' ,rmse(8),ccomma write(19,'(a17,f14.6,a1)') 'tdcverr= ' ,rmse(9),ccomma do nn=10,nflds if (nobs(nn)<=0 .or. bias(nn)>9998. .or. abse(nn)>9998. .or. rmse(nn)>9998.) then cycle else if (nn==10) write(19,'(a17,f14.6,a1)') 'gustcverr = ' ,rmse(10),ccomma if (nn==11) write(19,'(a17,f14.6,a1)') 'viscverr = ' ,rmse(11),ccomma if (nn==12) write(19,'(a17,f14.6,a1)') 'pblhcverr = ' ,rmse(12),ccomma if (nn==13) write(19,'(a17,f14.6,a1)') 'distcverr = ' ,rmse(13),ccomma if (nn==14) write(19,'(a17,f14.6,a1)') 'wspd10mcverr = ' ,rmse(14),ccomma if (nn==15) write(19,'(a17,f14.6,a1)') 'td2mcverr = ' ,rmse(15),ccomma if (nn==16) write(19,'(a17,f14.6,a1)') 'mxtmcverr = ' ,rmse(16),ccomma if (nn==17) write(19,'(a17,f14.6,a1)') 'mitmcverr = ' ,rmse(17),ccomma if (nn==18) write(19,'(a17,f14.6,a1)') 'pmslcverr = ' ,rmse(18),ccomma if (nn==19) write(19,'(a17,f14.6,a1)') 'howvcverr = ' ,rmse(19),ccomma if (nn==20) write(19,'(a17,f14.6,a1)') 'tcamtcverr = ' ,rmse(20),ccomma if (nn==21) write(19,'(a17,f14.6,a1)') 'lcbascverr = ' ,rmse(21),ccomma if (nn==22) write(19,'(a17,f14.6,a1)') 'uwnd10mcverr = ' ,rmse(22),ccomma if (nn==23) write(19,'(a17,f14.6,a1)') 'vwnd10mcverr = ' ,rmse(23),ccomma endif enddo write(19,'(a1)') '/' close(19) read(18) cmodel !this is errfield.dat_precval read(18) nflds,nobsmax0 allocate(nobs2(nflds)) !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist allocate(bias2(nflds)) !wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas,uwnd10m,vwnd10m allocate(abse2(nflds)) allocate(rmse2(nflds)) read(18) nobs2,bias2,abse2,rmse2 if (mype==0) then print*,'----------------------------' print*,'in cvbasedrecal: second read' print*,'cmodel=',trim(cmodel) print*,'nflds,nobsmax0=',nflds,nobsmax0 do nn=1,nflds print*,'nn,nobs2,bias2,abse2,rmse2=',& nn,nobs2(nn),bias2(nn),abse2(nn),rmse2(nn) enddo print*,'----------------------------' endif lres1=any(bias2(1:7) > 9998.) !any(bias2(1:nflds) > 9998.) lres2=any(abse2(1:7) > 9998.) !any(abse2(1:nflds) > 9998.) lres3=any(rmse2(1:7) > 9998.) !any(rmse2(1:nflds) > 9998.) if (lres1 .or. lres2 .or. lres3) then if (mype==0) then print*,'in cvbasedrecal/ second read : apparently no cross-validation was done' print*,' do nothing and return' endif return endif allocate(field(nx,ny)) do nn=1,nflds !!!!!! if (nn >=12) cycle if (nobs(nn)<=0 .or. bias(nn)>9998. .or. abse(nn)>9998. .or. rmse(nn)>9998.) cycle rmse0=max(rmse2(nn),epsilon) rmse1=max(rmse(nn),epsilon) ratio=rmse1/rmse0 !!note: no sqrt needed if (nn==11) ratio=min(ratio,1.33) !avoid unacceptably large analysis errors for vis / 29Jun2012 if (mype==0) print*,'in cvbasedrecal:nn,nobs(nn),rmse1,rmse0,(rmse1/rmse0)=',nn,nobs(nn),rmse1,rmse0,ratio if (nn==1) perr=perr*ratio if (nn==2) terr=terr*ratio if (nn==3) qerr=qerr*ratio if (nn==4) then uerr2=uerr2*ratio uerr=uerr*ratio !reasonable, I think endif if (nn==5) then verr2=verr2*ratio verr=verr*ratio !reasonable, I think endif if (nn==6) wspderr=wspderr*ratio ! if (nn==8) then ! wdirerr2=wdirerr2*ratio ! wdirerr=wdirerr*ratio !reasonable, I think ! endif if (nn==9) tderr=tderr*ratio if (nn==10) gusterr=gusterr*ratio if (nn==11) viserr=viserr*ratio if (nn==12) pblherr=pblherr*ratio if (nn==13) disterr=disterr*ratio if (nn==14) wspd10merr=wspd10merr*ratio if (nn==15) td2merr=td2merr*ratio if (nn==16) mxtmerr=mxtmerr*ratio if (nn==17) mitmerr=mitmerr*ratio if (nn==18) pmslerr=pmslerr*ratio if (nn==19) howverr=howverr*ratio if (nn==20) tcamterr=tcamterr*ratio if (nn==21) lcbaserr=lcbaserr*ratio if (nn==22) uwnd10merr=uwnd10merr*ratio if (nn==23) vwnd10merr=vwnd10merr*ratio !no need for adjustment for nn=24,25,26 (altw, altw2, altwd) since we !aren't yet estimating the analysis error for these fields enddo !note: in future, consider disseminating uerr2,verr2,wdierr2 !instead of uerr,verr,wdierr close(18) deallocate(nobs,nobs2) deallocate(bias,bias2) deallocate(abse,abse2) deallocate(rmse,rmse2) deallocate(field) return end !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- subroutine td_flds(nx,ny,fname,mype,npe) implicit none integer(4),intent(in):: nx,ny,mype,npe character(60),intent(in):: fname integer(4) rtime(6),nlon,nlat,nsig integer(4) ioan1 integer(4) n real(4),allocatable,dimension(:,:)::field1,field2 real(4),allocatable,dimension(:,:)::psfc,t,q,td1,td2,td0 real(4),parameter::qmin=1.e-06 logical fexist allocate(field1(nx,ny)) allocate(field2(nx,ny)) allocate(psfc(nx,ny)) allocate(t(nx,ny)) allocate(q(nx,ny)) allocate(td1(nx,ny)) allocate(td2(nx,ny)) allocate(td0(nx,ny)) if (mype==0) print*,'in td_flds: fname=',trim(fname) if (trim(fname) /= 'sigges' .and. & trim(fname) /= 'siganl' .and. & trim(fname) /= 'sigfupdate02' .and. & trim(fname) /= 'sigfupdate03' .and. & trim(fname) /= 'sigcress' ) then if (mype==0) print*,'invalid file name ... returning' return endif ioan1=51 !==> original first guess td. this is the field created by smartinit open (ioan1,file='slabs2.dat',form='unformatted') do n=1,4 read(ioan1) td0 enddo close(ioan1) if (mype==0) print*,'in td_flds: td0,min,max=',minval(td0),maxval(td0) if (mype==0) & open (57,file='td.dat_'//trim(fname),form='unformatted') !output file if (trim(fname) == 'sigges') then if (mype==0) then write(57) td0 close(57) endif return endif !==> alternative first guess td as derived from first guess psfc, q, and t. ! it is used to compute the td increments, which are then added to td0 to ! yield the total analysis td field open (ioan1,file='sigges',form='unformatted') read(ioan1) rtime,nlon,nlat,nsig read(ioan1) field1,field2 !glat,dx read(ioan1) field1,field2 !glon,dy if (mype==0) print*,'in td_flds: for sigges: rtime=',rtime if (mype==0) print*,'in td_flds: for sigges: nlon,nlat,nsig=',nlon,nlat,nsig read(ioan1) psfc read(ioan1) field1 !fis read(ioan1) t read(ioan1) q close(ioan1) call get_dewpt(psfc,q,t,td1,1,nx,1,ny) if (mype==0) print*,'in td_flds: td1,min,max=',minval(td1),maxval(td1) !==> td for siganl, sigfupdate02, sigfupdate03, or sigcress inquire(file=trim(fname),exist=fexist) if (fexist) then open (ioan1,file=trim(fname),form='unformatted') read(ioan1) rtime,nlon,nlat,nsig read(ioan1) field1,field2 !glat,dx read(ioan1) field1,field2 !glon,dy if (mype==0) print*,'in td_flds: for ',trim(fname),': rtime=',rtime if (mype==0) print*,'in td_flds: for ',trim(fname),': nlon,nlat,nsig=',nlon,nlat,nsig read(ioan1) psfc read(ioan1) field1 !fis read(ioan1) t read(ioan1) q close(ioan1) q=max(qmin,q) call get_dewpt(psfc,q,t,td2,1,nx,1,ny) td2 = td2 - td1 td2(:,:)=max(-20.,td2(:,:)) ! call smther_one(td2,1,nx,1,ny,2) td2 = td2 + td0 !td1 if (mype==0) then print*,'in td_flds: for ',trim(fname),': td2,min,max=',minval(td2),maxval(td2) write(57) td2 close(57) endif endif deallocate(field1) deallocate(field2) deallocate(psfc) deallocate(t) deallocate(q) deallocate(td1) deallocate(td2) deallocate(td0) return end !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- subroutine stn_interpolator(cgrid) use mpi use mpitaskmod, only: mype,npe use mpitaskmod, only: nx => nlon use mpitaskmod, only: ny => nlat use controlvars, only: ictrwspd10m => iwspd10m implicit none ! Declare parameters integer(4), parameter:: nt=7 !order is: p,t,td,sqrt(u**2+v***2),gust,vis,cldch real(4), parameter:: spval=-9999. real(4), parameter :: gravity=9.81 ! Declare passed variables character(60),intent(in)::cgrid ! Declare local variables integer(4) ifield,i,j,k,m,n,nstns integer(4) lun integer(4) n0,nlen integer(4) ii,jj integer(4) ierror real(8) alat1,elon1,dx,elonv,alatan real(8) rlat8,rlon8,xx8,yy8 character(120) cstring,cstring2,bstring character(8) cstn0 character(60) tdfname logical fexist real(4) fint0,altsetting character(40) cfldname(nt) character(80)cinterp_p,cinterp_t,cinterp_td,cinterp_wind, & cinterp_gust,cinterp_vis,cinterp_cldch real(4),allocatable,dimension(:):: rlat,rlon real(4),allocatable,dimension(:):: xx,yy,elev character(8),allocatable,dimension(:):: cstation character(80),allocatable,dimension(:):: cinterp_method character(10),allocatable,dimension(:,:):: cfint logical,allocatable,dimension(:,:):: viable logical,allocatable,dimension(:):: inside character(4) cyyyy0 character(2) cmm0,cdd0,chh0,cmin0 character(3) chh0Z,cmonth0,cmin0Z character(4) cyyyy character(2) cmm,cdd,chh,cmin character(3) chhZ,cmonth,cminZ character(4) cyyyyp1 character(2) cmmp1,cddp1,chhp1,cminp1 character(3) chhp1Z,cmonthp1,cminp1Z integer(4) jtime(6),nlon,nlat,nsig real(4) rdelz,half,half_tlapse,g_over_rd,rdp,rtvts real(4) tges,tob,pges,zsges real(4),allocatable,dimension(:,:)::field1, field2 real(4),allocatable,dimension(:,:)::terrain,tfield,tvfield character(60) cgprocess integer(4) iyyyy0,imm0,idd0,ihh0,imin0 !report issued at this time integer(4) iyyyy,imm,idd,ihh,imin !time range beginning time integer(4) iyyyyp1,immp1,iddp1,ihhp1,iminp1 !time range ending time namelist/generatingprocess/cgprocess namelist/faa_timerange/iyyyy0,imm0,idd0,ihh0,imin0, & iyyyy,imm,idd,ihh,imin, & iyyyyp1,immp1,iddp1,ihhp1,iminp1 namelist/faa_stninterp_method/cinterp_p,cinterp_t,cinterp_td,cinterp_wind, & cinterp_gust,cinterp_vis,cinterp_cldch data cgprocess/'RTMA'/ data cinterp_p /'bilinear'/ data cinterp_t /'bilinear'/ data cinterp_td /'bilinear'/ data cinterp_wind /'bilinear'/ data cinterp_gust /'bilinear'/ data cinterp_vis /'bilinear'/ !nearest_neighbor data cinterp_cldch /'bilinear'/ data cfldname(1) /'SFCP'/ data cfldname(2) /'2mT'/ data cfldname(3) /'2mTD'/ data cfldname(4) /'10mW'/ data cfldname(5) /'10mGUST'/ data cfldname(6) /'VIS'/ data cfldname(7) /'CLDCH'/ inquire(file='select_stnlist.dat',exist=fexist) if (.not.fexist) then if (mype==0) print*,'No request made for analysis interpolation to select locations' return endif inquire(file='faa_timerange_input',exist=fexist) if(fexist) then open (55,file='faa_timerange_input',form='formatted') read (55,faa_timerange) close(55) if (mype==0) then print*,'in stn_interpolator : iyyyy0,imm0,idd0,ihh0,imin0=', & iyyyy0,imm0,idd0,ihh0,imin0 print*,'in stn_interpolator : iyyyy,imm,idd,ihh,imin=', & iyyyy,imm,idd,ihh,imin print*,'in stn_interpolator : iyyyyp1,immp1,iddp1,ihhp1,iminp1=', & iyyyyp1,immp1,iddp1,ihhp1,iminp1 endif else if (mype==0) print*,'Missing valid time range for analysis interpolation!!!' return endif inquire(file='generatingprocess_input',exist=fexist) if(fexist) then open (55,file='generatingprocess_input',form='formatted') read (55,generatingprocess) close(55) endif if (trim(cgprocess)=='rtma') cgprocess='RTMA' if (trim(cgprocess)=='urma') cgprocess='URMA' if (mype==0) then print*,'in stn_interpolator: cgprocess=',trim(cgprocess) endif !------------------------------------ write(cyyyy0,"(i4.4)") iyyyy0 write(cmm0,"(i2.2)") imm0 write(cdd0,"(i2.2)") idd0 write(chh0,"(i2.2)") ihh0 write(cmin0,"(i2.2)") imin0 ; cmin0Z=cmin0//'Z' call getmonth(imm0,cmonth0) !------------------------------------ write(cyyyy,"(i4.4)") iyyyy write(cmm,"(i2.2)") imm write(cdd,"(i2.2)") idd write(chh,"(i2.2)") ihh write(cmin,"(i2.2)") imin ; cminZ=cmin//'Z' call getmonth(imm,cmonth) !------------------------------------ write(cyyyyp1,"(i4.4)") iyyyyp1 write(cmmp1,"(i2.2)") immp1 write(cddp1,"(i2.2)") iddp1 write(chhp1,"(i2.2)") ihhp1 write(cminp1,"(i2.2)") iminp1 ; cminp1Z=cminp1//'Z' call getmonth(immp1,cmonthp1) !------------------------------------ call proj_info(cgrid,dx,alat1,elon1,elonv,alatan) if (mype==0) then print*,'in stn_interpolator: cgrid=',trim(cgrid) print*,'in stn_interpolator: nx,ny=',nx,ny print*,'in stn_interpolator: alat1=',alat1 print*,'in stn_interpolator: elon1=',elon1 print*,'in stn_interpolator: dx=',dx print*,'in stn_interpolator: elonv=',elonv print*,'in stn_interpolator: alatan=',alatan endif ! Determine number of stations on faa list open (55,file='select_stnlist.dat',form='formatted') do n=1,3 read(55,*) cstring enddo nstns=0 read_faastn: do read(55,'(a8,1x,f5.2,1x,f7.2)',end=101) nstns=nstns+1 enddo read_faastn 101 continue if (mype==0) print*,'in stn_interpolator: nstns=',nstns allocate(cstation(nstns)) allocate(rlat(nstns)) allocate(rlon(nstns)) allocate(xx(nstns)) allocate(yy(nstns)) allocate(elev(nstns)) allocate(inside(nstns)) allocate(cfint(nstns,nt)) allocate(viable(nstns,nt)) allocate(cinterp_method(nt)) rewind(55) do n=1,3 read(55,*) cstring enddo inside=.true. do n=1,nstns read(55,'(a8,1x,f5.2,1x,f7.2,1x,f7.2)') cstation(n),rlat(n),rlon(n),elev(n) rlat8=dble(rlat(n)) rlon8=dble(rlon(n)) ; if(rlon8 < 0._8) rlon8=rlon8+360._8 call latlon_to_grid(dx,alat1,elon1,elonv,alatan, & rlat8,rlon8,cgrid,xx8,yy8) xx(n)=xx8 yy(n)=yy8 if (xx(n)<1. .or. xx(n)>float(nx) .or. yy(n)<1. .or. yy(n)>float(ny)) inside(n)=.false. enddo close(55) viable=.true. inquire(file='non_viable_stnlocation_list.dat',exist=fexist) if(fexist) then open (55,file='non_viable_stnlocation_list.dat',form='formatted') do n=1,3 read(55,*) bstring enddo do read(55,'(a)',end=181) bstring nlen=len_trim(bstring) cstn0=' ' do k=1,min(8,nlen) if (bstring(k:k)==',') exit cstn0(k:k)=bstring(k:k) enddo n0=-1 do n=1,nstns if (cstn0==cstation(n)) then n0=n exit endif enddo if (n0 > 0) then do k=1,nlen if (k+5<=nlen) then if ( bstring(k:k+5)=='P=skip' ) viable(n0,1)=.false. if ( bstring(k:k+5)=='T=skip' ) viable(n0,2)=.false. endif if (k+6<=nlen) then if ( bstring(k:k+6)=='TD=skip' ) viable(n0,3)=.false. endif if (k+8<=nlen) then if ( bstring(k:k+8)=='WIND=skip' ) viable(n0,4)=.false. if ( bstring(k:k+8)=='GUST=skip' ) viable(n0,5)=.false. endif if (k+7<=nlen) then if ( bstring(k:k+7)=='VIS=skip' ) viable(n0,6)=.false. endif if (k+9<=nlen) then if ( bstring(k:k+9)=='CLDCH=skip' ) viable(n0,7)=.false. endif enddo endif enddo 181 continue close(55) endif inquire(file='faa_stninterp_method_input',exist=fexist) if(fexist) then open (55,file='faa_stninterp_method_input',form='formatted') read (55,faa_stninterp_method) close(55) endif if (mype==0) then print*,'in stn_interpolator :' print*,'cinterp_p=',trim(cinterp_p) print*,'cinterp_t=',trim(cinterp_t) print*,'cinterp_td=',trim(cinterp_td) print*,'cinterp_wind=',trim(cinterp_wind) print*,'cinterp_gust=',trim(cinterp_gust) print*,'cinterp_vis=',trim(cinterp_vis) print*,'cinterp_cldch=',trim(cinterp_cldch) endif cinterp_method(1)=cinterp_p cinterp_method(2)=cinterp_t cinterp_method(3)=cinterp_td cinterp_method(4)=cinterp_wind cinterp_method(5)=cinterp_gust cinterp_method(6)=cinterp_vis cinterp_method(7)=cinterp_cldch if (mype==0) then open (75,file='stn_analysis_values.dat',form='formatted') open (76,file='stn_outside_of_domain.dat',form='formatted') endif ! cstring2=trim(cgprocess)//' 2m-temperature (degrees Celsius) and surface pressure (inHg)' cstring2=trim(cgprocess)//' 2m-temperature (degrees Celsius) and altimeter setting (inHg)' if (mype==0) then write(75,'(a)') trim(cstring) write(75,'(a)') trim(cstring2) write(75,'(a10,a2,a3,1x,a2,1x,a3,1x,a4)') 'COMPUTED: ',chh0,cmin0Z,cdd0,cmonth0,cyyyy0 write(75,'(a7,a2,a3,1x,a2,1x,a3,1x,a4,a4,a2,a3,1x,a2,1x,a3,1x,a4)') 'VALID: ',chh,cminZ,cdd,cmonth,cyyyy,' to ',chhp1,cminp1Z,cddp1,cmonthp1,cyyyyp1 write(75,'(a)') trim(cstring) ! write(75,'(a)') 'station Lat Lon T' ! write(75,'(a)') 'station Lat Lon 2m-T PS' ! write(75,'(a)') 'station Lat Lon 2m-T SFC PRES' write(75,'(a)') 'station Lat Lon 2m-T ALT' write(75,*) endif allocate(field1(nx,ny)) allocate(field2(nx,ny)) allocate(terrain(nx,ny)) allocate(tfield(nx,ny)) allocate(tvfield(nx,ny)) half=0.5 half_tlapse=0.00325 g_over_rd=9.81/287.04 !==> prepare dew point field inquire(file='td.dat_siganl',exist=fexist) if (.not.fexist) then tdfname='siganl' ; call td_flds(nx,ny,tdfname,mype,npe) endif call mpi_barrier(mpi_comm_world,ierror) !==>read in terrain, temperature, and calculate virtual temp open (53,file='siganl',form='unformatted') read(53) jtime,nlon,nlat,nsig read(53) field1,field2! glat,dx read(53) field1,field2! glon,dy if (mype==0) then print*,'in stn_interpolator-1 / from analysis' print*,'in stn_interpolatorn-1 / jtime=',jtime print*,'in stn_interpolatorn-1 / nlon,nlat,nsig=',nlon,nlat,nsig print*,'**********************************************' endif read(53) field1 !pressure read(53) field1 ; terrain=field1/gravity !!!; terrain=max(0.,terrain) read(53) tfield read(53) field1 !sphum do j=1,ny do i=1,nx tvfield(i,j)=tfield(i,j)*(1.+0.608*max(0.,field1(i,j))) !virtual temperature enddo enddo if (mype==0) then print*,'in stn_interpolator: terrain,min,max=',minval(terrain),maxval(terrain) print*,'in stn_interpolator: tfield,min,max=',minval(tfield),maxval(tfield) print*,'in stn_interpolator: tvfield,min,max=',minval(tvfield),maxval(tvfield) endif !==>second read in from signal to do the actual interpolation rewind(53) read(53) jtime,nlon,nlat,nsig read(53) field1,field2! glat,dx read(53) field1,field2! glon,dy if (mype==0) then print*,'in stn_interpolator / from analysis' print*,'in stn_interpolator / jtime=',jtime print*,'in stn_interpolator / nlon,nlat,nsig=',nlon,nlat,nsig print*,'**********************************************' endif do ifield=1,nt if (mype==0) then print*,'======================================================================' print*,'======================================================================' print*,'in stn_interpolator: ifield,cfldname=',ifield,trim(cfldname(ifield)) endif if (ifield==5) then !must skip 12 records to get to gust do m=1,12 read(53) enddo endif read(53) field1 if (ifield==3) then !replace q with td open (54,file='td.dat_siganl',form='unformatted') read(54) field1 close(54) endif if (ifield==4) then !u-wind read(53) field2 !v-wind if (ictrwspd10m > 0) then inquire(file='uv_adjusted.dat_siganl',exist=fexist) !make sure cnv_to_grib2_mpi.f90 is writing out this if (fexist) then !file if you are sending wind values to the FAA. open (94,file='uv_adjusted.dat_siganl',form='unformatted') !relevant statment in cnv_to_grib2_mpi.f90 read(94) field1 !u-wind !is currently commented out read(94) field2 !v-wind close(94) endif endif do j=1,ny do i=1,nx field1(i,j)=field1(i,j)*field1(i,j)+field2(i,j)*field2(i,j) if (field1(i,j) > 0.) field1(i,j)=sqrt(field1(i,j)) !wind speed enddo enddo field2=field1 !store wind speed endif if (ifield==5) then !analyzed gust cannot be less than wind speed do j=1,ny do i=1,nx if (field1(i,j) < field2(i,j)) field1(i,j)=field2(i,j) enddo enddo endif if (mype==0) then print*,'in stn_interpolator:field,min,max=',minval(field1),maxval(field1) print* print*,'station rlat rlon xx yy fint0' endif do n=1,nstns if (inside(n)) then if (viable(n,ifield)) then if (ifield==1) call retrieve_tob(mype,cstation(n),tob,rtvts,spval) if (trim(cinterp_method(ifield))=='nearest_neighbor') then ii=max(1,min(nint(xx(n)),nx)) jj=max(1,min(nint(yy(n)),ny)) fint0=field1(ii,jj) if (ifield==1) then zsges=terrain(ii,jj) !if (abs(rtvts)<0.001) then tges=tvfield(ii,jj) !else ! tges=tfield(ii,jj) !endif endif else call bilinear_2d0v2(field1,1,nx,1,ny,fint0,yy(n),xx(n)) if (ifield==1) then call bilinear_2d0v2(terrain,1,nx,1,ny,zsges,yy(n),xx(n)) !if (abs(rtvts)<0.001) then call bilinear_2d0v2(tvfield,1,nx,1,ny,tges,yy(n),xx(n)) !else ! call bilinear_2d0v2(tfield,1,nx,1,ny,tges,yy(n),xx(n)) !endif endif endif if (ifield==1) then if (mype==0) then print*,'in stn_interpolator: pressure case / non-adjusted values' print*,'cstation(n),elev(n),zsges=',trim(cstation(n)),elev(n),zsges print*,'cstation(n),tob,tges,fint0=',trim(cstation(n)),tob,tges,fint0 endif rdp=0. rdelz=elev(n)-zsges if (abs(tob-spval)>0.001) then tges = half*(tges+tob) else if(rdelz < 0.)then tges=tges-half_tlapse*rdelz endif endif rdp = g_over_rd*rdelz/tges !Adjust hydrostatically pges=fint0 fint0=exp(log(pges) - rdp) !output in Pa if (mype==0) then print*,'in stn_interpolator: pressure case / adjusted values' print*,'cstation(n),tges,fint0=',tges,fint0 endif endif ! if (ifield==1) fint0=fint0*0.0002953 !sfc pressure in In Hg if (ifield==1) fint0=altsetting(fint0,elev(n)) !altimeter setting in In Hg if (ifield==2 .or. ifield==3) fint0=fint0-273.15 !temp and dew point in Celsius dg if (ifield==4 .or. ifield==5) fint0=fint0*1.94384 !wind and gust in kts if (ifield==6) fint0=fint0*0.000621371 !vis in mi if (ifield==7) fint0=fint0*3.28084 !ceil in ft write(cfint(n,ifield),'(f10.2)') fint0 else cfint(n,ifield)='N/A' cfint(n,ifield)=adjustr(cfint(n,ifield)) endif lun=75 else fint0=-9999.99 write(cfint(n,ifield),'(f10.2)') fint0 lun=76 endif if (mype==0) print*,trim(cstation(n)),rlat(n),rlon(n),xx(n),yy(n),fint0 enddo if (mype==0) then print*,'======================================================================' print*,'======================================================================' endif if (ifield==1 .or. ifield==6) read(53) !skip fis and pblh enddo close(53) do n=1,nstns if (inside(n)) then lun=75 else lun=76 endif ! if (mype==0) write(lun,'(a8,f6.2,f8.2, a10)') cstation(n),rlat(n),rlon(n),cfint(n,2) if (mype==0) write(lun,'(a8,f6.2,f8.2, a10,a10)') cstation(n),rlat(n),rlon(n),cfint(n,2),cfint(n,1) enddo if (mype==0) then close(75) close(76) endif deallocate(cstation) deallocate(rlat) deallocate(rlon) deallocate(xx) deallocate(yy) deallocate(elev) deallocate(cfint) deallocate(field1) deallocate(inside) deallocate(field2) deallocate(terrain) deallocate(tfield) deallocate(tvfield) deallocate(viable) deallocate(cinterp_method) return end !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- real(4) function altsetting(p,elev0) implicit none real(4),intent(in):: p,elev0 real(8) p8mb, elev08, result8, a8, b8, c8 p8mb=dble(p)/100._8 elev08=dble(elev0) a8=(1013.25_8**0.190284_8 * 0.0065_8)/288._8 b8=elev08/(p8mb-0.3_8)**0.190284_8 c8= (p8mb-0.3_8)*(1._8+a8*b8)**(1._8/0.190284_8) c8=c8*100._8 !convert to Pa altsetting=sngl(c8)*0.0002953 !ALT in In Hg end function altsetting !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- subroutine retrieve_tob(mype, this_cstn, this_tob, this_qtflg, spval) ! ! to compile: ! ifort -c -convert big_endian retrieve_tob.f90 ! ! abstract: ! read in tob and surface height z0 from the anxal diagnostic file implicit none !Declare passed variables character(8),intent(in):: this_cstn real(4),intent(in):: spval integer(4),intent(in):: mype real(4),intent(out):: this_tob,this_qtflg !Note: ! this_shgt0 :station height ! this_hgt0 :observation elevation ! this_z0 :model terrain at ob location !Declare local parameters character(8),allocatable,dimension(:):: cdiagbuf character(8),allocatable,dimension(:):: cprvstg character(8),allocatable,dimension(:):: csprvstg character(8) cstation character(3) otype integer(4) idate,nchar,nreal,i,ii,mypeout real(4),allocatable,dimension(:,:)::rdiagbuf real(4) rlat,rlon real(4) dtime,qtflg,rmuse integer(4) itype logical fexist, r15887_diagfile_fmt integer(4) nlen1, nlen2 real(4) rsign,dtime0 real(4) this_dtime_t, this_dtime_q, this_dtime_ps real(4) this_qob,this_psob,this_psmodel_orig,this_psmodel,this_shgt0,this_hgt0,this_z0 integer(4) this_itype character(8) this_cstation logical ltob,lqob,lpsob namelist/diagfile_fmt/r15887_diagfile_fmt data r15887_diagfile_fmt/.true./ inquire(file='diagfile_fmt_input',exist=fexist) if (fexist) then open (55,file='diagfile_fmt_input',form='formatted') read(55,*) r15887_diagfile_fmt close(55) endif if (mype==0) print*,'in retrieve_tob: r15887_diagfile_fmt=',r15887_diagfile_fmt itype=nint(spval) rlat=spval rlon=spval this_tob=spval this_qtflg=spval this_qob=spval this_psob=spval this_psmodel=spval this_psmodel_orig=spval this_dtime_t=huge(this_dtime_t) this_dtime_q=huge(this_dtime_q) this_dtime_ps=huge(this_dtime_ps) this_shgt0=spval this_hgt0=spval this_z0=spval this_itype=nint(spval) this_cstation="88888888" nlen1=len_trim(this_cstn) open (7,file='diag_conv.dat',form='unformatted') read(7,end=200) idate loop_read_obs: do read(7,end=200)otype,nchar,nreal,ii,mypeout ! if (mype==0) print*,'in retrieve_tob: idate=',idate ! if (mype==0) print*,'in retrieve_tob: otype,nchar,nreal,ii,mypeout=', & ! otype,nchar,nreal,ii,mypeout if (ii==0) then read(7,end=200) if (r15887_diagfile_fmt) read(7,end=200) cycle loop_read_obs endif allocate(cdiagbuf(ii)) allocate(cprvstg(ii)) allocate(csprvstg(ii)) allocate(rdiagbuf(nreal,ii)) if (r15887_diagfile_fmt) then read(7) cdiagbuf,rdiagbuf read(7) cprvstg,csprvstg else read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg endif ltob=otype(1:3)==' t' lqob=otype(1:3)==' q' lpsob=otype(2:3)=='ps' if (ltob.or.lqob.or.lpsob) then do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) dtime=rdiagbuf(8,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rsign=1. if (rdiagbuf(12,i) /= 0 ) rsign=rdiagbuf(12,i)/abs(rdiagbuf(12,i)) rmuse=rdiagbuf(12,i)+rsign*(rdiagbuf(11,i)-float(int(rdiagbuf(11,i)))) rmuse=rdiagbuf(12,i) endif if (ltob) then dtime0=this_dtime_t elseif (lqob) then dtime0=this_dtime_q elseif (lpsob) then dtime0=this_dtime_ps endif nlen2=len_trim(cstation) if ( (itype.eq.187.or.itype.eq.193) .and. (rmuse.gt.0.) .and. (nlen1.eq.nlen2) .and. & this_cstn(1:nlen1)==cstation(1:nlen1).and. (abs(dtime).lt.abs(dtime0)) ) then this_itype=itype !used for diagnostic purposes only this_cstation=cstation !used for diagnostic purposes only if (ltob) then rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) this_shgt0=rdiagbuf(5,i) this_hgt0=rdiagbuf(7,i) this_tob=rdiagbuf(17,i) this_qtflg=rdiagbuf(10,i) this_dtime_t=dtime this_z0=rdiagbuf(22,i) elseif (lqob) then this_qob=rdiagbuf(17,i) this_dtime_q=dtime elseif (lpsob) then this_psob=rdiagbuf(17,i) this_psmodel=rdiagbuf(17,i)-rdiagbuf(18,i) this_psmodel_orig=rdiagbuf(17,i)-rdiagbuf(19,i) this_dtime_ps=dtime endif endif enddo end if deallocate(cdiagbuf,rdiagbuf) deallocate(cprvstg,csprvstg) enddo loop_read_obs 200 continue close(7) if (mype==0) then print*,'in in retrieve_tob' print*,'this_cstn,this_cstation=',trim(this_cstn),trim(this_cstation) print*,'rlat,rlon=',rlat,rlon print*,'this_itype=',this_itype print*,'this_dtime_t,this_dtime_q,this_dtime_ps=',this_dtime_t,this_dtime_q,this_dtime_ps print*,'this_tob,this_qtflg=',this_tob,this_qtflg print*,'this_qob,this_psob,this_psmodel_orig,this_psmodel=',this_qob,this_psob,this_psmodel_orig,this_psmodel print*,'this_shgt0,this_hgt0,this_z0=',this_shgt0,this_hgt0,this_z0 endif return end !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- subroutine getmonth(n,month) implicit none integer(4),intent(in)::n character(3),intent(out)::month if (n==1) month='Jan' if (n==2) month='Feb' if (n==3) month='Mar' if (n==4) month='Apr' if (n==5) month='May' if (n==6) month='Jun' if (n==7) month='Jul' if (n==8) month='Aug' if (n==9) month='Sep' if (n==10) month='Oct' if (n==11) month='Nov' if (n==12) month='Dec' return end !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- subroutine mp_flush(i) implicit none integer(4),intent(in) :: i return end !--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------