DIFFERENCES FOR anberror.f90 30a31,32 > ! def ims - > ! def ime - 34a37,38 > ! def jms - > ! def jme - 38a43,44 > ! def jms - > ! def kme - 83,85d88 < ! def rtma_subdomain_option - if true, then apply recursive filters in subdomain space < ! (currently constructed to work only when running 2d analysis < ! and analysis and filter grid are the same) 103a107 > integer(i_kind) ims,ime,jms,jme,kms,kme ! subdomain + halo lat, lon, vert grid indices (full units) 105c109 < integer(i_kind),allocatable::idvar(:),jdvar(:),kvar_start(:),kvar_end(:),levs_jdvar(:) --- > integer(i_kind),allocatable::idvar(:),jdvar(:),kvar_start(:),kvar_end(:) 108c112 < type(filter_cons),save:: filter_all(7) --- > type(filter_cons) filter_all(7) 110c114 < integer(i_long) ifilt_ord,npass,ngauss,normal,nsmooth,nsmooth_shapiro --- > integer(i_long) ifilt_ord,npass,ngauss,normal 116,119d119 < logical rtma_subdomain_option < integer(i_kind),allocatable:: nrecv_halo(:),ndrecv_halo(:),nsend_halo(:),ndsend_halo(:) < integer(i_kind) nrecv_halo_loc,nsend_halo_loc < integer(i_kind),allocatable:: info_send_halo(:,:),info_recv_halo(:,:) 154,156c154,156 < ids=0 ; ide=0 ; ips=0 ; ipe=0 < jds=0 ; jde=0 ; jps=0 ; jpe=0 < kds=0 ; kde=0 ; kps=0 ; kpe=0 --- > ids=0 ; ide=0 ; ims=0 ; ime=0 ; ips=0 ; ipe=0 > jds=0 ; jde=0 ; jms=0 ; jme=0 ; jps=0 ; jpe=0 > kds=0 ; kde=0 ; kms=0 ; kme=0 ; kps=0 ; kpe=0 177,178d176 < nsmooth=0 < nsmooth_shapiro=0 190d187 < rtma_subdomain_option=.false. 280d276 < use gridmod, only: istart,jstart,nlon,nlat 289,292c285 < if(rtma_subdomain_option) then < call halo_update_reg0(mype) < ! call test_halo_update(mype) < call anberror_vert_partition_subdomain_option(mype) --- > call anberror_vert_partition(mype) 294,302c287 < nlatf=nlat ; nlonf=nlon < ids=1 ; ide=nlat < jds=1 ; jde=nlon < ! following without halo < ips=max(ids,min(istart(mype+1),ide)) ; ipe=max(ids,min(lat2+istart(mype+1)-3,ide)) < jps=max(jds,min(jstart(mype+1),jde)) ; jpe=max(jds,min(lon2+jstart(mype+1)-3,jde)) < < else < call anberror_vert_partition(mype) --- > ! initialize fgrid2agrid interpolation constants 304c289 < ! initialize fgrid2agrid interpolation constants --- > call create_fgrid2agrid(grid_ratio) 306c291,294 < call create_fgrid2agrid(grid_ratio) --- > ids=1 ; ide=nlatf > jds=1 ; jde=nlonf > ims=ids ; ime=ide ; ips=ids ; ipe=ide > jms=jds ; jme=jde ; jps=jds ; jpe=jde 308,314d295 < ids=1 ; ide=nlatf < jds=1 ; jde=nlonf < ips=ids ; ipe=ide < jps=jds ; jpe=jde < < end if < 423c404,405 < write(6,*)' in anberror_vert_partition, kps,kpe=',kps,kpe --- > kms=kps ; kme=kpe > write(6,*)' in anberror_vert_partition, kps,kms,kpe,kme=',kps,kms,kpe,kme 453,704d434 < subroutine anberror_vert_partition_subdomain_option(mype) < < ! to use anisotropic filter code directly from subdomains, need different < ! vertical definition < < use kinds, only: i_kind < use gridmod, only: nsig < implicit none < < integer(i_kind),intent(in):: mype < < integer(i_kind) idvar_last,k,kk,vlevs < < ! vlevs=6*nsig+4 ! all variables < vlevs=4*nsig+1 ! for rtma, currently do only u,v,psfc,t,q < kds=1 ; kde=vlevs < kps=kds ; kpe=kde < < ! initialize nvars,idvar,kvar_start,kvar_end < ! Determine how many vertical levels each mpi task will < ! handle in the horizontal smoothing < ! nvars=10 < nvars=5 < allocate(idvar(kds:kde),jdvar(kds:kde),kvar_start(nvars),kvar_end(nvars)) < allocate(var_names(nvars),levs_jdvar(kds:kde)) < var_names( 1)="st" < var_names( 2)="vp" < var_names( 3)="ps" < var_names( 4)="tv" < var_names( 5)="q" < ! var_names( 6)="oz" < ! var_names( 7)="sst" < ! var_names( 8)="stl" < ! var_names( 9)="sti" < ! var_names(10)="cw" < < ! idvar jdvar < ! 1 1 stream function < ! 2 2 velocity potential < ! 3 3 surface pressure < ! 4 4 virtual temperature < ! 5 5 specific humidity < !! 6 6 ozone < !! 7 7 sst < !! 10 8 cloud water < !! 8 9 surface temp (land) < !! 9 10 surface temp (ice) < idvar(1 :nsig ) =1 ! st < idvar(nsig+1 :2*nsig) =2 ! vp < idvar(2*nsig+1) =3 ! ps < idvar(2*nsig+2:3*nsig+1)=4 ! tv < idvar(3*nsig+2:4*nsig+1)=5 ! q < jdvar=idvar < < kk=0 < do k=1,nsig < kk=kk+1 < levs_jdvar(kk)=k ! st < end do < do k=1,nsig < kk=kk+1 < levs_jdvar(kk)=k ! vp < end do < kk=kk+1 < levs_jdvar(kk)=1 ! ps < do k=1,nsig < kk=kk+1 < levs_jdvar(kk)=k ! tv < end do < do k=1,nsig < kk=kk+1 < levs_jdvar(kk)=k ! q < end do < < kvar_start(1) =1 < kvar_end(1) =nsig < kvar_start(2) =kvar_end(1)+1 < kvar_end(2) =kvar_end(1)+nsig < kvar_start(3) =kvar_end(2)+1 < kvar_end(3) =kvar_end(2)+1 < kvar_start(4) =kvar_end(3)+1 < kvar_end(4) =kvar_end(3)+nsig < kvar_start(5) =kvar_end(4)+1 < kvar_end(5) =kvar_end(4)+nsig < < if(mype.eq.0) then < do k=kds,kde < write(6,*)' in anberror_vert_partition_subdomain_option, k,idvar,jdvar,levs_jdvar=', & < k,idvar(k),jdvar(k),levs_jdvar(k) < end do < do k=1,nvars < write(6,*)' k,kvar_start,end(k)=',k,kvar_start(k),kvar_end(k) < end do < end if < write(6,*)' in anberror_vert_partition_subdomain_option, kps,kpe=',kps,kpe < < end subroutine anberror_vert_partition_subdomain_option < < subroutine halo_update_reg0(mype) < < use kinds, only: r_kind,i_kind < use gridmod, only: lat2,lon2,istart,jstart,nlat,nlon < use mpimod, only: npe,mpi_integer4,mpi_sum,mpi_comm_world,ierror < use raflib, only: indexxi4 < implicit none < < integer(i_kind),intent(in):: mype < < integer(i_kind) i,ii,j,k,mm1,mpe,iglob,jglob,mpi_string1 < integer(i_kind) ijglob_pe(nlat,nlon),ijglob_pe0(nlat,nlon) < integer(i_kind) iorigin(3*(lat2+lon2)),indx(3*(lat2+lon2)),iwork(3*(lat2+lon2)) < < allocate(nrecv_halo(0:npe-1),ndrecv_halo(0:npe),nsend_halo(0:npe-1),ndsend_halo(0:npe)) < allocate(info_send_halo(2,3*(lat2+lon2)),info_recv_halo(2,3*(lat2+lon2))) < < < if(npe.eq.1) return < mm1=mype+1 < < ijglob_pe0=0 < do j=2,lon2-1 < jglob=j+jstart(mm1)-2 < do i=2,lat2-1 < iglob=i+istart(mm1)-2 < ijglob_pe0(iglob,jglob)=mype < end do < end do < call mpi_allreduce(ijglob_pe0,ijglob_pe,nlat*nlon,mpi_integer4,mpi_sum,mpi_comm_world,ierror) < < ! create list of all points to be received with global i,j coordinates < ii=0 < nrecv_halo=0 < !ierror=0 < do j=1,lon2,lon2-1 < jglob=j+jstart(mm1)-2 < if(jglob.lt.1.or.jglob.gt.nlon) cycle < do i=1,lat2 < iglob=i+istart(mm1)-2 < if(iglob.lt.1.or.iglob.gt.nlat) cycle < ii=ii+1 < info_recv_halo(1,ii)=iglob ; info_recv_halo(2,ii)=jglob < iorigin(ii)=ijglob_pe(iglob,jglob) < nrecv_halo(ijglob_pe(iglob,jglob))=nrecv_halo(ijglob_pe(iglob,jglob))+1 < ! if(iorigin(ii).eq.mype) ierror=ierror+1 < end do < end do < do i=1,lat2,lat2-1 < iglob=i+istart(mm1)-2 < if(iglob.lt.1.or.iglob.gt.nlat) cycle < do j=2,lon2-1 ! already have corner points < jglob=j+jstart(mm1)-2 < if(jglob.lt.1.or.jglob.gt.nlon) cycle < ii=ii+1 < info_recv_halo(1,ii)=iglob ; info_recv_halo(2,ii)=jglob < iorigin(ii)=ijglob_pe(iglob,jglob) < nrecv_halo(ijglob_pe(iglob,jglob))=nrecv_halo(ijglob_pe(iglob,jglob))+1 < ! if(iorigin(ii).eq.mype) ierror=ierror+1 < end do < end do < < ndrecv_halo(0)=0 < do mpe=1,npe < ndrecv_halo(mpe)=ndrecv_halo(mpe-1)+nrecv_halo(mpe-1) < end do < < call mpi_alltoall(nrecv_halo,1,mpi_integer4,nsend_halo,1,mpi_integer4,mpi_comm_world,ierror) < ndsend_halo(0)=0 < do mpe=1,npe < ndsend_halo(mpe)=ndsend_halo(mpe-1)+nsend_halo(mpe-1) < end do < nsend_halo_loc=ndsend_halo(npe) < nrecv_halo_loc=ndrecv_halo(npe) < < ! sort origin pe numbers from smallest to largest < if(ii.gt.0) then < call indexxi4(ii,iorigin,indx) < < ! use sort index to reorder < do j=1,2 < do i=1,ii < iwork(i)=info_recv_halo(j,indx(i)) < end do < do i=1,ii < info_recv_halo(j,i)=iwork(i) < end do < end do < end if < < call mpi_type_contiguous(2,mpi_integer4,mpi_string1,ierror) < call mpi_type_commit(mpi_string1,ierror) < call mpi_alltoallv(info_recv_halo,nrecv_halo,ndrecv_halo,mpi_string1, & < info_send_halo,nsend_halo,ndsend_halo,mpi_string1,mpi_comm_world,ierror) < call mpi_type_free(mpi_string1,ierror) < < < ! convert info arrays back to local coordinate units < ! ierror=0 < do i=1,nsend_halo_loc < iglob=info_send_halo(1,i) < info_send_halo(1,i)=iglob-istart(mm1)+2 < ! if(info_send_halo(1,i).gt.lat2-1.or.info_send_halo(1,i).lt.2) ierror=ierror+1 < jglob=info_send_halo(2,i) < info_send_halo(2,i)=jglob-jstart(mm1)+2 < ! if(info_send_halo(2,i).gt.lon2-1.or.info_send_halo(2,i).lt.2) ierror=ierror+1 < end do < ! ierror=0 < do i=1,nrecv_halo_loc < iglob=info_recv_halo(1,i) < info_recv_halo(1,i)=iglob-istart(mm1)+2 < ! if(info_recv_halo(1,i).gt.lat2.or.info_recv_halo(1,i).lt.1) ierror=ierror+1 < jglob=info_recv_halo(2,i) < info_recv_halo(2,i)=jglob-jstart(mm1)+2 < ! if(info_recv_halo(2,i).gt.lon2.or.info_recv_halo(2,i).lt.1) ierror=ierror+1 < end do < < end subroutine halo_update_reg0 < < subroutine halo_update_reg(f,nvert) < < use kinds, only: r_kind,i_kind < use gridmod, only: lat2,lon2 < use mpimod, only: npe,mpi_rtype,mpi_comm_world,ierror < implicit none < < integer(i_kind),intent(in):: nvert < real(r_kind),intent(inout):: f(lat2,lon2,nvert) < < integer(i_kind) i,k,mpi_string2 < real(r_kind) bufsend(nvert,nsend_halo_loc),bufrecv(nvert,nrecv_halo_loc) < < if(npe.eq.1) return < < ! now gather up points to send < do i=1,nsend_halo_loc < do k=1,nvert < bufsend(k,i)=f(info_send_halo(1,i),info_send_halo(2,i),k) < end do < end do < call mpi_type_contiguous(nvert,mpi_rtype,mpi_string2,ierror) < call mpi_type_commit(mpi_string2,ierror) < call mpi_alltoallv(bufsend,nsend_halo,ndsend_halo,mpi_string2, & < bufrecv,nrecv_halo,ndrecv_halo,mpi_string2,mpi_comm_world,ierror) < call mpi_type_free(mpi_string2,ierror) < ! finally distribute points back < do i=1,nrecv_halo_loc < do k=1,nvert < f(info_recv_halo(1,i),info_recv_halo(2,i),k)=bufrecv(k,i) < end do < end do < < end subroutine halo_update_reg < 706,754d435 < < subroutine test_halo_update(mype) < < use kinds,only: r_kind,i_kind < use anberror, only: halo_update_reg < use gridmod, only: lat2,lon2,istart,jstart,nlon,nlat < implicit none < < integer(i_kind),intent(in):: mype < < real(r_kind) f(lat2,lon2,2) < integer(i_kind) i,iglob,j,jglob,ierror < < f=0._r_kind < do j=2,lon2-1 < jglob=j+jstart(mype+1)-2 < do i=2,lat2-1 < iglob=i+istart(mype+1)-2 < f(i,j,1)=iglob < f(i,j,2)=jglob < end do < end do < ierror=0 < do j=1,lon2 < jglob=j+jstart(mype+1)-2 < if(jglob.lt.1.or.jglob.gt.nlon) cycle < do i=1,lat2 < iglob=i+istart(mype+1)-2 < if(iglob.lt.1.or.iglob.gt.nlat) cycle < if(nint(f(i,j,1)).ne.iglob.or.nint(f(i,j,2)).ne.jglob) ierror=ierror+1 < end do < end do < write(0,*)' in test_halo_update, before halo update ierror=',ierror < call halo_update_reg(f,2) < !call halo_update_reg(f(1,1,1),1) < !call halo_update_reg(f(1,1,2),1) < ierror=0 < do j=1,lon2 < jglob=j+jstart(mype+1)-2 < if(jglob.lt.1.or.jglob.gt.nlon) cycle < do i=1,lat2 < iglob=i+istart(mype+1)-2 < if(iglob.lt.1.or.iglob.gt.nlat) cycle < if(nint(f(i,j,1)).ne.iglob.or.nint(f(i,j,2)).ne.jglob) ierror=ierror+1 < end do < end do < write(0,*)' in test_halo_update, after halo update ierror=',ierror < < end subroutine test_halo_update DIFFERENCES FOR anberror.f90_double_precision DIFFERENCES FOR anberror.f90_no_changes DIFFERENCES FOR anbkerror_reg.f90 118,119d117 < use anberror, only: rtma_subdomain_option < use constants, only: zero 151c149,151 < if(rtma_subdomain_option) then --- > ! Convert from subdomain to full horizontal field distributed among processors > iflg=1 > call sub2grid(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp,iflg) 153,158c153,154 < oz=zero < cwmr=zero < sst=zero < slndt=zero < sicet=zero < call ansmoothrf_reg_subdomain_option(t,p,q,st,vp,mype) --- > ! Apply horizontal smoother for number of horizontal scales > call ansmoothrf_reg(hwork,mype) 160c156,157 < else --- > ! Put back onto subdomains > call grid2sub(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp) 162,173d158 < ! Convert from subdomain to full horizontal field distributed among processors < iflg=1 < call sub2grid(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp,iflg) < < ! Apply horizontal smoother for number of horizontal scales < call ansmoothrf_reg(hwork,mype) < < ! Put back onto subdomains < call grid2sub(hwork,t,p,q,oz,sst,slndt,sicet,cwmr,st,vp) < < end if < 287c272 < filter_all,ngauss --- > ims,ime,jms,jme,kms,kme,filter_all,ngauss 297c282 < real(r_kind),dimension(nlat,nlon,kps:kpe),intent(inout):: work --- > real(r_kind),dimension(nlat,nlon,kms:kme),intent(inout):: work 301,302c286,287 < real(r_kind),dimension(ips:ipe,jps:jpe,kps:kpe):: worka < real(r_single),dimension(ngauss,ips:ipe,jps:jpe,kps:kpe):: workb --- > real(r_kind),dimension(ims:ime,jms:jme,kms:kme):: worka > real(r_single),dimension(ngauss,ims:ime,jms:jme,kms:kme):: workb 306,307c291,292 < do k=kps,kpe < call tfgrid2agrid(work(1,1,k),worka(ips,jps,k)) --- > do k=kms,kme > call tfgrid2agrid(work(1,1,k),worka(ims,jms,k)) 312,314c297,299 < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe --- > do k=kms,kme > do j=jms,jme > do i=ims,ime 324,325c309,312 < call raf4(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe) < call raf4_ad(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe) --- > call raf4(workb,filter_all,ngauss, & > ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe,ims,ime,jms,jme,kms,kme,mype,npe) > call raf4_ad(workb,filter_all,ngauss, & > ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe,ims,ime,jms,jme,kms,kme,mype,npe) 328,330c315,317 < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe --- > do k=kms,kme > do j=jms,jme > do i=ims,ime 340,341c327,328 < do k=kps,kpe < call fgrid2agrid(worka(ips,jps,k),work(1,1,k)) --- > do k=kms,kme > call fgrid2agrid(worka(ims,jms,k),work(1,1,k)) 467,649d453 < < subroutine ansmoothrf_reg_subdomain_option(t,p,q,st,vp,mype) < !$$$ subprogram documentation block < ! . . . . < ! subprogram: ansmoothrf_reg_subdomain_option anisotropic rf for regional mode < ! prgmmr: parrish org: np22 date: 2005-02-14 < ! < ! abstract: apply anisotropic rf for regional mode (using subdomains instead of slabs) < ! NOTE: only works if filter grid is same as analysis grid < ! < ! program history log: < ! 2005-02-14 parrish < ! < ! input argument list: < ! t,p,q,oz,st,slndt,sicet,cwmr,st,vp - fields to be smoothed < ! < ! output argument list: < ! t,p,q,oz,st,slndt,sicet,cwmr,st,vp - smoothed fields < ! < ! attributes: < ! language: f90 < ! machine: ibm RS/6000 SP < !$$$ < use kinds, only: r_kind,i_kind,r_single < use anberror, only: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe, & < filter_all,ngauss,halo_update_reg < use mpimod, only: npe < use constants, only: zero < use gridmod, only: lat2,lon2,istart,jstart,nsig < use raflib, only: raf4_ad,raf4 < implicit none < < ! Declare passed variables < integer(i_kind), intent(in):: mype < real(r_kind),dimension(lat2,lon2),intent(inout):: p < real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: t,q,vp,st < < ! Declare local variables < integer(i_kind) i,igauss,iloc,j,jloc,k,kk,mm1 < real(r_single),dimension(ngauss,ips:ipe,jps:jpe,kps:kpe):: workb < < mm1=mype+1 < < ! transfer variables to ngauss copies < kk=0 < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < do igauss=1,ngauss < workb(igauss,i,j,kk)=st(iloc,jloc,k) < end do < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < do igauss=1,ngauss < workb(igauss,i,j,kk)=vp(iloc,jloc,k) < end do < end do < end do < end do < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < do igauss=1,ngauss < workb(igauss,i,j,kk)=p(iloc,jloc) < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < do igauss=1,ngauss < workb(igauss,i,j,kk)=t(iloc,jloc,k) < end do < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < do igauss=1,ngauss < workb(igauss,i,j,kk)=q(iloc,jloc,k) < end do < end do < end do < end do < < ! apply recursive filter < < call raf4(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe) < call raf4_ad(workb,filter_all,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npe) < < ! add together ngauss copies < kk=0 < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < st(iloc,jloc,k)=zero < do igauss=1,ngauss < st(iloc,jloc,k)=st(iloc,jloc,k)+workb(igauss,i,j,kk) < end do < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < vp(iloc,jloc,k)=zero < do igauss=1,ngauss < vp(iloc,jloc,k)=vp(iloc,jloc,k)+workb(igauss,i,j,kk) < end do < end do < end do < end do < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < p(iloc,jloc)=zero < do igauss=1,ngauss < p(iloc,jloc)=p(iloc,jloc)+workb(igauss,i,j,kk) < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < t(iloc,jloc,k)=zero < do igauss=1,ngauss < t(iloc,jloc,k)=t(iloc,jloc,k)+workb(igauss,i,j,kk) < end do < end do < end do < end do < do k=1,nsig < kk=kk+1 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < q(iloc,jloc,k)=zero < do igauss=1,ngauss < q(iloc,jloc,k)=q(iloc,jloc,k)+workb(igauss,i,j,kk) < end do < end do < end do < end do < < ! update halos: < call halo_update_reg(p,1) < call halo_update_reg(t,nsig) < call halo_update_reg(q,nsig) < call halo_update_reg(vp,nsig) < call halo_update_reg(st,nsig) < < end subroutine ansmoothrf_reg_subdomain_option DIFFERENCES FOR anbkerror_reg.f90_double_precision DIFFERENCES FOR anbkerror_reg.f90_no_changes DIFFERENCES FOR anbkerror_reg.f90_with_shapiro DIFFERENCES FOR anisofilter.f90 47a48 > ims,ime,jms,jme,kms,kme, & 50c51 < nvars,idvar,jdvar,kvar_start,kvar_end,levs_jdvar, & --- > nvars,idvar,jdvar,kvar_start,kvar_end, & 53,54c54 < ngauss,rgauss,an_amp,an_vs,an_flen_u,an_flen_t,an_flen_z,& < rtma_subdomain_option,nsmooth,nsmooth_shapiro --- > ngauss,rgauss,an_amp,an_vs,an_flen_u,an_flen_t,an_flen_z 76d75 < use mpimod, only: mpi_real4 84c83 < logical,parameter:: lreadnorm = .false. --- > logical,parameter:: lreadnorm = .true. 175,179c174 < if(rtma_subdomain_option) then < call get2berr_reg_subdomain_option(mype) < else < call get2berr_reg(mype) < end if --- > call get2berr_reg(mype) 231d225 < integer(i_kind) kk 233d226 < real(4) aspect_out(nlon,nlat) 249,250d241 < character(80) fname < 341,361d331 < !?????????????????????????? < ! output map of aspect(.,.,.,levs_id(1) < !do k=1,nsig1o < ! if(nvar_id(k).eq.1) then < ! do kk=1,7 < ! write(0,*)' output map of aspect(',kk,') for st' < ! do j=1,nlonf < ! do i=1,nlatf < ! aspect_out(j,i)=aspect(kk,i,j,k) < ! end do < ! end do < ! write(fname,'("aspect_st0_",i1)')kk < ! call outgrads1(aspect_out,nlon,nlat,fname) < ! end do < ! end if < !end do < ! if(mype.gt.-100) then < ! call mpi_finalize(i) < ! stop < ! end if < 405d374 < nsmooth,nsmooth_shapiro, & 408a378 > ims, ime, jms, jme, kms, kme, & ! memory indices 412a383 > 421,422d391 < if(i.eq.253.and.j.eq.341) write(0,*)' slab run i253,j341,igauss=', & < igauss,trim(chvarname(ivar)),' amp=',filter_all(1)%amp(igauss,i,j,k) 431,434d399 < ! if(mype.gt.-100) then < ! call mpi_finalize(i) < ! stop < ! end if 481,483c446,448 < ! if (j.eq.(jps+2) .and. (i.eq.(ips+2) .or. i.eq.(ips+8))) then < ! print*,'in anisofilter,mype,k,ivar,i,j,factk=',mype,k,ivar,i,j,factk < ! endif --- > if (j.eq.(jps+2) .and. (i.eq.(ips+2) .or. i.eq.(ips+8))) then > print*,'in anisofilter,mype,k,ivar,i,j,factk=',mype,k,ivar,i,j,factk > endif 485,486d449 < if(i.ge.641.and.i.le.651.and.j.eq.75.and.ivar.eq.4) write(0,*)' slab run i,j,igauss=', & < i,j,igauss,trim(chvarname(ivar)),' factor,amp=',factor,filter_all(1)%amp(igauss,i,j,k) 755d717 < nsmooth,nsmooth_shapiro, & 758a721 > ims, ime, jms, jme, kms, kme, & ! memory indices 1077d1039 < integer(i_kind) nsmooth_smooth,nsmooth_shapiro_smooth 1305,1306d1266 < nsmooth_smooth=0 < nsmooth_shapiro_smooth=0 1310d1269 < nsmooth_smooth,nsmooth_shapiro_smooth, & 1313a1273 > ims, ime, jms, jme, kms, kme, & ! memory indices 1331,1334c1291,1310 < call raf_sm4(field0f,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) < call raf_sm4_ad(field0f,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) < call raf_sm4(field0zf,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) < call raf_sm4_ad(field0zf,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) --- > call raf_sm4(field0f,filter_all,ngauss_smooth, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npe) > call raf_sm4_ad(field0f,filter_all,ngauss_smooth, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npe) > call raf_sm4(field0zf,filter_all,ngauss_smooth, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npe) > call raf_sm4_ad(field0zf,filter_all,ngauss_smooth, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npe) 1350,1351c1326,1335 < ! call raf_sm4(z0f,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) < ! call raf_sm4_ad(z0f,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps,kpe,mype,npe) --- > ! call raf_sm4(z0f,filter_all,ngauss_smooth, & > ! ids, ide, jds, jde, kds, kde, & ! domain indices > ! ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ! ims, ime, jms, jme, kms, kme, & ! memory indices > ! mype, npe) > ! call raf_sm4_ad(z0f,filter_all,ngauss_smooth, & > ! ids, ide, jds, jde, kds, kde, & ! domain indices > ! ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ! ims, ime, jms, jme, kms, kme, & ! memory indices > ! mype, npe) 1847,2763d1830 < < subroutine get2berr_reg_subdomain_option(mype) < !$$$ subprogram documentation block < ! . . . . < ! subprogram: get2berr_reg < ! prgmmr: parrish org: np22 date: 2005-02-08 < ! < ! abstract: setup everything for anisotropic regional background < ! error (2dvar case). < ! < ! program history log: < ! 2005-02-08 parrish < ! 2005-03-28 wu - replace mlath with mlat and setup qoption=2 < ! 2005-05-24 pondeca - add reconstruction of aspect tensor via < ! eigenvector decomposition < ! 2005-11-29 derber - unify ozone variance calculations < ! 2006-02-02 treadon - rename prsl as ges_prsl < ! 2006-04-17 treadon - use rlsig from call rdgstat_reg; replace sigl < ! with ges_prslavg/ges_psfcavg < ! 2006-07-01 pondeca - add terrain following covariances. < ! remove eigenvector decomposition < ! 2007-05-30 h.liu - remove ozmz, use factoz < ! < ! input argument list: < ! mype - mpi task id < ! < ! output argument list: < ! < ! < ! attributes: < ! language: f90 < ! machine: ibm RS/6000 SP < ! < !$$$ < implicit none < < ! Declare passed variables < integer(i_kind),intent(in):: mype < < ! Declare local parameters < real(r_kind),parameter:: r25 = 1.0_r_kind/25.0_r_kind !this was put in 200708 code < !no need. r25 already defined < < ! Declare local variables < integer(i_kind) i,j,k,l,lp,k1,kvar,ivar,im,ip,jm,jp,mm1,iloc,iploc,imloc,jloc,jploc,jmloc < integer(i_kind) iglob,jglob < integer(i_kind) kk < < real(4) aspect_out(nlon,nlat),aspect_out0(nlon,nlat) < real(r_kind) d,dl1,dl2,factk,factor,hwll_loc < real(r_kind) a1,a2,a3,a4,a5,a6,detai < real(r_kind) biga1,biga2,biga3,biga4,biga5,biga6 < real(r_kind) fx2,fx1,fx3,dxi,dyi,dzi < real(r_kind) asp1,asp2,asp3,factoz < real(r_kind) detratio,deta0,deta1 < < real(r_kind),allocatable,dimension(:,:)::bckgvar,zaux,zsmooth < real(r_kind),allocatable,dimension(:,:,:):: bckgvar0f < real(r_single),allocatable,dimension(:,:)::bckgvar4,bckgvar4a,zsmooth4,zsmooth4a < real(r_single),allocatable,dimension(:,:)::region_dx4,region_dy4,psg4,psg4a < character*10 chvarname(5) < character(80) fname < real(r_single) dummy_amp < < data chvarname/'psi', 'chi', 'ps', 't', 'pseudorh'/ < < < call read_bckgstats_subdomain_option(mype) < call get_background_subdomain_option(mype) < < !-----optionally, remove latitudinal dependence from statistics < l=int(rllatf(nlatf/2,nlonf/2)) < ! do j=0,mlat+1 < ! if (j>=1 .and. j<= mlat) corz(j,:,:)=corz(l,:,:) < ! if (j>=1 .and. j<= mlat) corp(j)=corp(l) < ! hwll(j,:,:)=hwll(l,:,:) < ! hwllp(j)=hwllp(l) < ! end do < < !-----define the anisotropic aspect tensor----------------------------- < !------------------------------------------------------------------------------------- < ! Set up scales < < ! This first loop for nsig1o will be if we aren't dealing with < ! surface pressure, skin temperature, or ozone < < rltop_wind=900._r_kind*6._r_kind! 5._r_kind! 2.5_r_kind! *5._r_kind !2. < rltop_temp=100._r_kind! *1.75_r_kind! *2._r_kind < rltop_q=100._r_kind! *1.75_r_kind! *2._r_kind < rltop_psfc=150._r_kind! *2._r_kind < < allocate(aspect(7,ips:ipe,jps:jpe,kps:kpe)) < aspect(:,:,:,:)=zero < mm1=mype+1 < < do k=kds,kde < < k1=levs_jdvar(k) < call isotropic_scales_subdomain_option(asp10f,asp20f,asp30f,k,mype) < < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < < asp1=asp10f(iloc,jloc) < asp2=asp20f(iloc,jloc) < asp3=asp30f(iloc,jloc) < < if(i.eq.nlatf/2.and.j.eq.nlonf/2) then < write(6,*)'("at domain center, var,k1,asp1,asp2,asp3 =",2i4,3f11.3)', & < jdvar(k),k1,asp1,asp2,asp3 < write(6,*)'("at domain center, var,k1,dxf,dyf =",2i4,3f11.3)', & < jdvar(k),k1,dxf(iloc,jloc),dyf(iloc,jloc) < end if < jp=min(nlonf,j+1) ; jm=max(1,j-1) < jploc=jp-jstart(mm1)+2 < jmloc=jm-jstart(mm1)+2 < dxi=one/(jp-jm) < < ip=min(nlatf,i+1) ; im=max(1,i-1) < iploc=ip-istart(mm1)+2 < imloc=im-istart(mm1)+2 < dyi=one/(ip-im) < < fx1= dyi*(z0f(iploc,jloc,k1)-z0f(imloc,jloc,k1)) < fx2= dxi*(z0f(iloc,jploc,k1)-z0f(iloc,jmloc,k1)) < fx3= zero < < rltop=rltop_wind < afact=zero < < if (jdvar(k)==1) rltop=rltop_wind < if (jdvar(k)==2) rltop=rltop_wind < if (jdvar(k)==3) rltop=rltop_psfc < if (jdvar(k)==4) rltop=rltop_temp < if (jdvar(k)==5) rltop=rltop_q < if (jdvar(k) <= 5) afact=one !(use "zero" for isotropic computations) < < if (afact==one) then < asp1=scalex1*asp1 < asp2=scalex2*asp2 < endif < < < aspect(1,i,j,k) = 1./asp1**2 + afact*fx1*fx1/rltop**2 ! 1st (y) direction x1*x1 < aspect(2,i,j,k) = 1./asp2**2 + afact*fx2*fx2/rltop**2 ! 2nd (x) direction x2*x2 < aspect(3,i,j,k) = 1./asp3**2 ! 3rd (z) direction x3*x3 < aspect(4,i,j,k) = afact*fx3*fx2/rltop**2 ! x3*x2 < aspect(5,i,j,k) = afact*fx3*fx1/rltop**2 ! x3*x1 < aspect(6,i,j,k) = afact*fx2*fx1/rltop**2 ! x2*x1 < aspect(7,i,j,k)= zero < < deta0=one/(asp1*asp2)**2 < deta1=aspect(1,i,j,k)*aspect(2,i,j,k)-aspect(6,i,j,k)*aspect(6,i,j,k) < aspect(1,i,j,k) = aspect(1,i,j,k)/sqrt(deta1)*sqrt(deta0) < aspect(2,i,j,k) = aspect(2,i,j,k)/sqrt(deta1)*sqrt(deta0) < aspect(6,i,j,k) = aspect(6,i,j,k)/sqrt(deta1)*sqrt(deta0) < < end do < end do < end do < < !??????????????????????////// < ! output map of aspect(.,.,.,1) for examination < ! aspect_out0=zero < ! do k=1,7 < ! do j=jps,jpe < ! do i=ips,ipe < ! aspect_out0(j,i)=aspect(k,i,j,1) < ! end do < ! end do < ! call mpi_reduce(aspect_out0,aspect_out,nlat*nlon,mpi_real4,mpi_sum,0,mpi_comm_world,i) < ! if(mype.eq.0) then < ! write(fname,'("aspect_st_",i1)')k < ! call outgrads1(aspect_out,nlon,nlat,fname) < ! end if < ! end do < ! if(mype.gt.-100) then < ! call mpi_finalize(i) < ! stop < ! end if < < < < ! Invert to get true aspect tensor < < do k=kps,kpe < if(levs_jdvar(k).eq.0 ) cycle < do j=jps,jpe < do i=ips,ipe < a1=aspect(1,i,j,k) < a2=aspect(2,i,j,k) < a3=aspect(3,i,j,k) < a4=aspect(4,i,j,k) < a5=aspect(5,i,j,k) < a6=aspect(6,i,j,k) < biga1=a2*a3-a4*a4 < biga2=a1*a3-a5*a5 < biga3=a1*a2-a6*a6 < biga4=a5*a6-a1*a4 < biga5=a4*a6-a2*a5 < biga6=a4*a5-a3*a6 < detai=one/(a1*biga1+a6*biga6+a5*biga5) < aspect(1,i,j,k)=biga1*detai < aspect(2,i,j,k)=biga2*detai < aspect(3,i,j,k)=biga3*detai < aspect(4,i,j,k)=biga4*detai < aspect(5,i,j,k)=biga5*detai < aspect(6,i,j,k)=biga6*detai < end do < end do < end do < < < if(mype.eq.0) write(6,*)'rltop_wind,rltop_temp,rltop_q,rltop_psfc=',& < rltop_wind,rltop_temp,rltop_q,rltop_psfc < < if(mype.eq.0) write(6,*)' in get2berr_reg, nlat,nlon,nlatf,nlonf=',nlat,nlon,nlatf,nlonf < < if(mype.eq.0) write(6,*)' in get2berr_reg, ids,ide=',ids,ide < if(mype.eq.0) write(6,*)' in get2berr_reg, jds,jde=',jds,jde < write(6,*)'in get2berr_reg, mype,ips,ipe,jps,jpe,kps,kpe=',mype,ips,ipe,jps,jpe,kps,kpe < < if(lreadnorm) normal=0 < !?????????????????????special test-- turn off normalization, then read in from slab run < ! normal=0 < < call init_raf4(aspect,triad4,ngauss,rgauss,npass,normal,binom,ifilt_ord,filter_all, & < nsmooth,nsmooth_shapiro, & < nvars,idvar,kvar_start,kvar_end,var_names, & < ids, ide, jds, jde, kds, kde, & ! domain indices < ips, ipe, jps, jpe, kps, kpe, & ! patch indices < mype, npe) < < call antest_maps0_subdomain_option(mype,theta0f,z0f) < < ! write(fname,'("fltnorm.dat_",i4.4)')mype < ! open (94,file=trim(fname),form='unformatted') < ! if(lreadnorm) then < ! read(94) filter_all(1)%amp < ! else < ! write(94) filter_all(1)%amp < ! end if < ! close(94) < < ! do k=1,5 < ! if(k.eq.1) fname='fltnorm.dat_psi' < ! if(k.eq.2) fname='fltnorm.dat_chi' < ! if(k.eq.3) fname='fltnorm.dat_ps' < ! if(k.eq.4) fname='fltnorm.dat_t' < ! if(k.eq.5) fname='fltnorm.dat_pseudorh' < ! open (94,file=trim(fname),form='unformatted') < ! do j=jds,jde < ! do i=ids,ide < ! read(94) dummy_amp < ! if(i.ge.ips.and.i.le.ipe.and.j.ge.jps.and.j.le.jpe) then < ! filter_all(1)%amp(1,i,j,k)=dummy_amp < ! if(i.eq.253.and.j.eq.341) write(0,*)' subdomain run i253,j341,',trim(fname),' amp=',dummy_amp < ! end if < ! end do < ! end do < ! close(94) < ! end do < ! if(mype.gt.-100) then < ! call mpi_finalize(i) < ! stop < ! end if < < allocate(bckgvar0f(ips:ipe,jps:jpe,kps:kpe)) < bckgvar0f=zero < < ! filter normalized to unit amplitude. now adjust amplitude to correspond < ! to desired error amplitude. < ! < ! (corz, corp contain sqrt(error variance) from background error file) < ! (an_amp are input tuneable error amplitude parameters) < < factoz = 0.0002_r_kind*r25 < do k=kps,kpe < if(k.eq.1) fname='fltnorm.dat_psi' < if(k.eq.2) fname='fltnorm.dat_chi' < if(k.eq.3) fname='fltnorm.dat_ps' < if(k.eq.4) fname='fltnorm.dat_t' < if(k.eq.5) fname='fltnorm.dat_pseudorh' < ivar=jdvar(k) < kvar=levs_jdvar(k) < do j=jps,jpe < do i=ips,ipe < l=max(min(int(rllatf(i,j)),mlat),1) < lp=min((l+1),mlat) < dl2=rllatf(i,j)-float(l) < dl1=one-dl2 < if(ivar.eq. 1) factk=dl1*corz(l,kvar,1)+dl2*corz(lp,kvar,1) ! streamfunction < if(ivar.eq. 2) factk=dl1*corz(l,kvar,2)+dl2*corz(lp,kvar,2) ! velocity potential < if(ivar.eq. 3) factk=dl1*corp(l)+dl2*corp(lp) ! ps < if(ivar.eq. 4) factk=dl1*corz(l,kvar,3)+dl2*corz(lp,kvar,3) ! temperature < if(ivar.eq. 5) factk=dl1*corz(l,kvar,4)+dl2*corz(lp,kvar,4) ! specific humidity < < ! Note: ideally, correct an_amp should come from < ! namelist anbkgerr < < if (ivar == 1) an_amp(:,ivar)=0.35_r_double < if (ivar == 2) an_amp(:,ivar)=0.35_r_double*2.063_r_double < if (ivar == 3) an_amp(:,ivar)=0.70_r_double*2._r_double < if (ivar == 4) an_amp(:,ivar)=1.00_r_double*2._r_double < if (ivar == 5) an_amp(:,ivar)=0.5_r_double*1.5_r_double*2._r_double < < do igauss=1,ngauss < if (j.eq.(jps+2) .and. (i.eq.(ips+2) .or. i.eq.(ips+8))) then < print*,'in anisofilter,mype,k,ivar,i,j,factk=',mype,k,ivar,i,j,factk < endif < factor=factk*an_amp(igauss,ivar) < if(i.ge.641.and.i.le.651.and.j.eq.75.and.ivar.eq.4) write(0,*)' sub run i,j,igauss=', & < i,j,igauss,trim(chvarname(ivar)),' factor,amp=',factor,filter_all(1)%amp(igauss,i,j,k) < filter_all(1)%amp(igauss,i,j,k)=factor*filter_all(1)%amp(igauss,i,j,k) < bckgvar0f(i,j,k)=factor < end do < end do < end do < end do < < < ! write out a few bckg fields and the error variances. For layer constant < ! statistics, these are good for post-processing purposes on the real < ! model grid. Interface between filter grid and model domain needed < ! for non-constant statistics! < < allocate(bckgvar4a(nlat,nlon)) < allocate(bckgvar4(nlat,nlon)) < < < do k=kps,kpe < bckgvar4a=zero < do j=jps,jpe < do i=ips,ipe < bckgvar4a(i,j)=bckgvar0f(i,j,k) < end do < end do < call mpi_reduce(bckgvar4a,bckgvar4,nlat*nlon,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < if(mype.eq.0) then < ivar=jdvar(k) < open (94,file='bckgvar.dat_'//trim(chvarname(ivar)),form='unformatted') < write(94) bckgvar4 < close(94) < end if < enddo < < allocate(zsmooth4a(nlat,nlon)) < allocate(zsmooth4(nlat,nlon)) < zsmooth4a=zero < do j=2,lon2-1 < jglob=j+jstart(mm1)-2 < do i=2,lat2-1 < iglob=i+istart(mm1)-2 < zsmooth4a(iglob,jglob)=z0f(i,j,1) < end do < end do < call mpi_reduce(zsmooth4a,zsmooth4,nlat*nlon,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < < allocate(psg4a(nlat,nlon)) < allocate(psg4(nlat,nlon)) < psg4a=zero < do j=2,lon2-1 < jglob=j+jstart(mm1)-2 < do i=2,lat2-1 < iglob=i+istart(mm1)-2 < psg4a(iglob,jglob)=psg(i,j,1) < end do < end do < call mpi_reduce(psg4a,psg4,nlat*nlon,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < < allocate(region_dy4(nlat,nlon),region_dx4(nlat,nlon)) < < if (mype.eq.0) then < region_dx4=region_dx < region_dy4=region_dy < open (94,file='bckg_dxdy.dat',form='unformatted') < write(94) region_dx4,region_dy4 < close(94) < < open (94,file='bckg_psfc.dat',form='unformatted') < write(94) psg4 < close(94) < < open (94,file='bckg_z.dat',form='unformatted') < write(94) zsmooth4 < close(94) < end if < < < deallocate(bckgvar0f) < deallocate(bckgvar4a) < deallocate(bckgvar4) < deallocate(region_dy4,region_dx4) < deallocate(corz,corp,hwll,hwllp,vz,aspect) < deallocate(dxf,dyf,rllatf,theta0f) < deallocate(u0f,v0f,z0f) < deallocate(zsmooth4a,psg4a) < deallocate(zsmooth4,psg4) < < end subroutine get2berr_reg_subdomain_option < !======================================================================= < !======================================================================= < subroutine read_bckgstats_subdomain_option(mype) < !$$$ subprogram documentation block < ! . . . . < ! subprogram: read_bckgstats < ! prgmmr: pondeca org: np22 date: 2006-08-01 < ! < ! abstract: read in background error statistics. built from parrish's < ! old anprewgt_reg < ! < ! < ! program history log: < ! 2006-08-01 pondeca < ! < ! input argument list: < ! mype - mpi task id < ! < ! output argument list: < ! < ! < ! attributes: < ! language: f90 < ! machine: ibm RS/6000 SP < ! < implicit none < < ! Declare passed variables < integer(i_kind),intent(in):: mype < < ! Declare local variables < integer(i_kind) j,k,l,n,inerr < < real(r_kind),allocatable::rlsig(:),dlsig(:) < < ! Read dimension of stats file < inerr=22 < open(inerr,file='berror_stats',form='unformatted') < rewind inerr < read(inerr)msig,mlat < < ! Allocate arrays in stats file < allocate ( corz(1:mlat,1:nsig,1:4) ) < allocate ( corp(1:mlat) ) < allocate ( hwll(0:mlat+1,1:nsig,1:4),hwllp(0:mlat+1) ) < allocate ( vz(1:nsig,0:mlat+1,1:6) ) < allocate ( agvi(0:mlat+1,1:nsig,1:nsig) ) < allocate ( bvi(0:mlat+1,1:nsig),wgvi(0:mlat+1,1:nsig) ) < < ! Read in background error stats and interpolate in vertical < ! to that specified in namelist < < allocate(rlsig(nsig)) < allocate(dlsig(nsig)) < < call rdgstat_reg(msig,mlat,inerr,& < hwll,hwllp,vz,agvi,bvi,wgvi,corz,corp,rlsig) < deallocate(agvi,bvi,wgvi) < close(inerr) < < write(6,*)'in read_bckgstats,mlat=',mlat < < allocate(vzimax(1:nsig,1:6)) < allocate(vzimin(1:nsig,1:6)) < allocate(vziavg(1:nsig,1:6)) < allocate(corzavg(1:nsig,1:4)) < allocate(hwllavg(1:nsig,1:4)) < < < if(qoption==2)then < ! varq(1:mlat,1:nsig)=corz(1:mlat,1:nsig,4) < ! corz(1:mlat,1:nsig,4)=1. < varq(1:mlat,1:nsig)=min(max(corz(1:mlat,1:nsig,4),0.015_r_kind),one) < corz(llmin:llmax,1:nsig,4)=one < endif < < ! Normalize vz with del sigma and convert to vertical grid units! < if (.not.twodvar_regional) then < dlsig(1)=rlsig(1)-rlsig(2) < do k=2,nsig-1 < dlsig(k)=half*(rlsig(k-1)-rlsig(k+1)) < enddo < dlsig(nsig)=rlsig(nsig-1)-rlsig(nsig) < else < dlsig=one ! Really no meaning for 2dvar. Set to 1.0 to avoid < ! division by zero below. Array vz is reset for 2dvar < ! case, so vz calculation below is not truly need. < end if < < do k=1,nsig < vz(k,0:mlat+1,1:6)=vz(k,0:mlat+1,1:6)*dlsig(k) < end do < < call create_anberror_vars_reg(mype) < < !----- apply scaling to vertical length scales. < ! note: parameter vs needs to be inverted < < write(6,*)'in read_bckgstats,an_vs=',an_vs < an_vs=one/an_vs < vz=vz/an_vs < if (twodvar_regional) vz(:,:,:)=sqrt(one) < < < !-----compute and print out diagnostics for < ! vertical length scales < < do n=1,6 < do k=1,nsig < vzimax(k,n)=maxval(one/vz(k,0:mlat+1,n)) < vzimin(k,n)=minval(one/vz(k,0:mlat+1,n)) < vziavg(k,n)=sum((one/vz(k,0:mlat+1,n)))/float(mlat+2) < end do < < if(mype.eq.0) then < do k=1,nsig < write(6,*)'(" var,k,max,min,avg vert corlen =",2i4,3f11.3)', & < n,k,vzimax(k,n),vzimin(k,n),vziavg(k,n) < end do < end if < end do < < !-----optionally, remove latitudinal dependence from statistics < if (.not.latdepend) then < < do n=1,4 < do k=1,nsig < corzavg(k,n)=sum(corz(1:mlat,k,n))/float(mlat) < hwllavg(k,n)=sum(hwll(0:mlat+1,k,n))/float(mlat+2) < end do < end do < corpavg=sum(corp(1:mlat))/float(mlat) < hwllpavg=sum(hwllp(0:mlat+1))/float(mlat+2) < < do j=1,mlat < corz(j,1:nsig,1:4)=corzavg(1:nsig,1:4) < corp(j)=corpavg < end do < do j=0,mlat+1 < hwll(j,1:nsig,1:4)=hwllavg(1:nsig,1:4) < hwllp(j)=hwllpavg < vz(1:nsig,j,1:6)=one/vziavg(1:nsig,1:6) < end do < < end if < < ! hybrid sigma level structure calculated in rdgstat_reg < ! ks used to load constant horizontal scales for SF/VP < ! above sigma level 0.15 < ! loop l for diff variable of each PE. < < psfc015=r015*ges_psfcavg < allocate (ks(kds:kde)) < do l=kds,kde < ks(l)=nsig+1 < if(jdvar(l)<3)then < k_loop: do k=1,nsig < if (ges_prslavg(k)< psfc015) then < ks(l)=k < exit k_loop < end if < enddo k_loop < endif < end do < < end subroutine read_bckgstats_subdomain_option < !======================================================================= < !======================================================================= < subroutine get_background_subdomain_option(mype) < !$$$ subprogram documentation block < ! . . . . < ! subprogram: get_background < ! prgmmr: pondeca org: np22 date: 2006-08-01 < ! < ! abstract: compute background fields and their spatial derivatives < ! on filter grid for use in anisotropic covariance model. < ! built from parrish's old anprewgt_reg < ! < ! < ! program history log: < ! 2006-08-01 pondeca < ! < ! input argument list: < ! mype - mpi task id < ! < ! output argument list: < ! < ! < ! attributes: < ! language: f90 < ! machine: ibm RS/6000 SP < ! < use anberror, only: halo_update_reg < implicit none < < ! Declare passed variables < integer(i_kind),intent(in):: mype < < ! Declare local variables < integer(i_kind) i,j,k,l,ip,im,jp,jm,kp,km,mm1,iflg,k1,ivar,lp < integer(i_kind) n,iloc,jloc < integer(i_kind) kds0,kde0,kps0,kpe0 < integer(i_kind) nvars0 < integer(i_kind),allocatable::idvar0(:),kvar_start0(:),kvar_end0(:) < character(80),allocatable::var_names0(:) < integer(i_kind) nsmooth_smooth,nsmooth_shapiro_smooth < < real(r_kind) d,dl1,dl2,factk < real(r_kind) asp1,asp2,asp3 < < real(r_kind),allocatable,dimension(:,:,:)::field2 < real(r_single),allocatable,dimension(:,:,:)::field < < ! get dxf,dyf,rllatf < ! note: if filter grid coarser than analysis grid, then normalized < ! adjoint of filter to analysis interpolation used to transfer fields < ! from analysis grid to filter grid. otherwise, normal interpolation < ! is done. this is transparent at this level. it appears in the < ! definition of the interpolation and adjoint of interpolation < ! weights. check for accuracy.(done). < < if(nlat.ne.nlatf.or.nlon.ne.nlonf) then < if(mype.eq.0) & < write(6,*)' rtma_subdomain_option true, nlat ne nlatf and/or nlon ne nlonf, nlat,nlatf,nlon,nlonf=', & < nlat,nlatf,nlon,nlonf < call mpi_finalize(i) < stop < end if < < allocate(dxf(nlatf,nlonf),dyf(nlatf,nlonf),rllatf(nlatf,nlonf)) < < dxf=region_dx < dyf=region_dy < rllatf=rllat < if(mype.eq.0) then < write(6,*)' at 11.28 in anprewgt_reg, nlat,nlon,nlatf,nlonf=',nlat,nlon,nlatf,nlonf < write(6,*)' at 11.29 in anprewgt_reg, min,max(rllat)=',minval(rllat),maxval(rllat) < write(6,*)' at 11.3 in anprewgt_reg, min,max(rllatf)=',minval(rllatf),maxval(rllatf) < write(6,*)' at 11.31 in anprewgt_reg, min,max(grid_ratio_lon*dx)=', & < minval(grid_ratio_lon*region_dx),maxval(grid_ratio_lon*region_dx) < write(6,*)' at 11.32 in anprewgt_reg, min,max(dxf)=',minval(dxf),maxval(dxf) < write(6,*)' at 11.33 in anprewgt_reg, min,max(grid_ratio_lat*dy)=', & < minval(grid_ratio_lat*region_dy),maxval(grid_ratio_lat*region_dy) < write(6,*)' at 11.34 in anprewgt_reg, min,max(dyf)=',minval(dyf),maxval(dyf) < end if < < mm1=mype+1 < < ! define local vertical index for smoothing single variables < kds0=1 < kde0=nsig < kps0=kds0 < kpe0=kde0 < nvars0=1 < allocate(idvar0(kds0:kde0),kvar_start0(nvars0),kvar_end0(nvars0)) < allocate(var_names0(nvars0)) < idvar0=1 < kvar_start0=1 < kvar_end0=nsig < var_names0(1)='background' < ! ------------------------------------------------------------ < ! ------------ in this section, set up isotropic filter for < ! ------------ generating smoothed guess < ! ------------------------------------------------------------ < < allocate(aspect(7,ips:ipe,jps:jpe,kps0:kpe0)) < < do k=kps0,kpe0 < do j=jps,jpe < do i=ips,ipe < < aspect(1,i,j,k)=1.5_r_kind**2 !3.! 10. < aspect(2,i,j,k)=1.5_r_kind**2 !3.! 10. < aspect(3,i,j,k)=1._r_kind**2 < aspect(4:7,i,j,k)=zero < < end do < end do < end do < < ngauss_smooth=1 < rgauss_smooth=one < npass_smooth=1 < normal_smooth=0 < ifilt_ord_smooth=4 < nsmooth_smooth=0 < nsmooth_shapiro_smooth=0 < call init_raf4(aspect,triad4,ngauss_smooth,rgauss_smooth, & < npass_smooth,normal_smooth,binom, & < ifilt_ord_smooth,filter_all, & < nsmooth_smooth,nsmooth_shapiro_smooth, & < nvars0,idvar0,kvar_start0,kvar_end0,var_names0, & < ids, ide, jds, jde, kds0, kde0, & ! domain indices < ips, ipe, jps, jpe, kps0, kpe0, & ! patch indices < mype, npe) < < it=ntguessig < < allocate(field(ips:ipe,jps:jpe,kps0:kpe0),field2(lat2,lon2,nsig)) < allocate(theta0f(lat2,lon2,nsig)) < allocate( u0f(lat2,lon2,nsig)) < allocate( v0f(lat2,lon2,nsig)) < allocate( z0f(lat2,lon2,nsig)) < do n=1,4 < < do k=kps0,kpe0 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < if (n==1) field(i,j,k)=ges_tv(iloc,jloc,k,it)/(ges_prsl(iloc,jloc,k,it)/r100)**rd_over_cp < if (n==2) field(i,j,k)=ges_u(iloc,jloc,k,it) < if (n==3) field(i,j,k)=ges_v(iloc,jloc,k,it) < if (n==4) field(i,j,k)=ges_z(iloc,jloc,it) < end do < end do < end do < call raf_sm4(field,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps0,kpe0,mype,npe) < call raf_sm4_ad(field,filter_all,ngauss_smooth,ips,ipe,jps,jpe,kps0,kpe0,mype,npe) < do k=kps0,kpe0 < do j=jps,jpe < jloc=j-jstart(mm1)+2 < do i=ips,ipe < iloc=i-istart(mm1)+2 < field2(iloc,jloc,k)=field(i,j,k) < end do < end do < end do < call halo_update_reg(field2,nsig) < do k=1,nsig < do j=1,lon2 < do i=1,lat2 < if (n==1) theta0f(i,j,k)=field2(i,j,k) < if (n==2) u0f(i,j,k)=field2(i,j,k) < if (n==3) v0f(i,j,k)=field2(i,j,k) < if (n==4) z0f(i,j,k)=field2(i,j,k) < end do < end do < end do < end do < allocate(psg(lat2,lon2,nsig)) < do k=1,nsig < do j=1,lon2 < do i=1,lat2 < psg(i,j,k)=1000._r_single*ges_ps(i,j,it) < end do < end do < end do < < allocate(asp10f(lat2,lon2)) < allocate(asp20f(lat2,lon2)) < allocate(asp30f(lat2,lon2)) < < deallocate(field) < deallocate(field2) < deallocate(aspect) < deallocate(idvar0,kvar_start0,kvar_end0,var_names0) < < end subroutine get_background_subdomain_option < !======================================================================= < !======================================================================= < subroutine isotropic_scales_subdomain_option(scale1,scale2,scale3,k,mype) < !$$$ subprogram documentation block < ! . . . . < ! subprogram: isotropic_scales_subdomain_option < ! prgmmr: pondeca org: np22 date: 2008-01-27 < ! < ! abstract: compute isotropic length scales of auto-correlation model. < ! built from parrish's old anprewgt_reg < ! < ! program history log: < ! 2006-08-01 pondeca < ! < ! input argument list: < ! mype - mpi task id < ! k - level number of field in subdomain mode < ! < ! output argument list: < ! scale1 - 2d field of correlations lengths in the x-direction < ! scale2 - 2d field of correlations lengths in the y-direction < ! scale3 - 2d field of correlations lengths in the vertical direction < ! < ! < ! attributes: < ! language: f90 < ! machine: ibm RS/6000 SP < ! < implicit none < < !Declare passed variables < integer(i_kind) k, mype < < !Declare local variables < integer(i_kind) i,j,k1,ivar,l,lp,iglob,jglob,mm1 < < real(r_kind) scale1(lat2,lon2) < real(r_kind) scale2(lat2,lon2) < real(r_kind) scale3(lat2,lon2) < real(r_kind) dl1,dl2,factk,hwll_loc < < < mm1=mype+1 < if (jdvar(k)==1) ivar=1 ! streamfunction < if (jdvar(k)==2) ivar=2 ! velocity potential < if (jdvar(k)==4) ivar=3 ! temperature < if (jdvar(k)==5) ivar=4 ! specific humidity < if (jdvar(k)==6) ivar=5 ! Ozone < if (jdvar(k)==7) ivar=1 ! SST < if (jdvar(k)==8) ivar=4 ! cloud water < if (jdvar(k)==9) ivar=1 ! surface temp (land) < if (jdvar(k)==10) ivar=1 ! surface temp (ice) < < k1=levs_jdvar(k) < < do j=1,lon2 < jglob=j+jstart(mm1)-2 < if(jglob.lt.1.or.jglob.gt.nlonf) cycle < do i=1,lat2 < iglob=i+istart(mm1)-2 < if(iglob.lt.1.or.iglob.gt.nlatf) cycle < < if( jdvar(k)==1 .or. jdvar(k)==2 .or. jdvar(k)==4 .or. & < jdvar(k)==5 .or. jdvar(k)==8 ) then < if(k1 >= ks(k))then < l=int(rllatf(nlatf/2,nlonf/2)) < hwll_loc=hwll(l,k1,ivar) < else < l=int(rllatf(iglob,jglob)) < lp=l+1 < dl2=rllatf(iglob,jglob)-float(l) < dl1=one-dl2 < hwll_loc=dl1*hwll(l,k1,ivar)+dl2*hwll(lp,k1,ivar) < end if < scale3(i,j)=one/vz(k1,l,ivar) < < else if (jdvar(k)==3) then !surface pressure < l=int(rllatf(iglob,jglob)) < lp=l+1 < dl2=rllatf(iglob,jglob)-float(l) < dl1=one-dl2 < hwll_loc=dl1*hwllp(l)+dl2*hwllp(lp) < scale3(i,j)=one < < else if (jdvar(k)==6) then < if(k1 <= nsig*3/4)then < hwll_loc=r400000 < else < hwll_loc=(r800000-r400000*(nsig-k1)/(nsig-nsig*3/4)) < end if < l=int(rllatf(nlatf/2,nlonf/2)) < scale3(i,j)=one/vz(k1,l,ivar) < < else if (jdvar(k)==7) then ! SST < l=int(rllatf(iglob,jglob)) < lp=l+1 < dl2=rllatf(iglob,jglob)-float(l) < dl1=one-dl2 < hwll_loc=half*(dl1*hwll(l,1,ivar)+dl2*hwll(lp,1,ivar)) < scale3(i,j)=one < < else if (jdvar(k)==9) then ! surface temp (land) < l=int(rllatf(iglob,jglob)) < lp=l+1 < dl2=rllatf(iglob,jglob)-float(l) < dl1=one-dl2 < hwll_loc=zero_25*(dl1*hwll(l,1,ivar)+dl2*hwll(lp,1,ivar)) < scale3(i,j)=one < < else if (jdvar(k)==10) then ! surface temp (ice) < l=int(rllatf(iglob,jglob)) < lp=l+1 < dl2=rllatf(iglob,jglob)-float(l) < dl1=one-dl2 < hwll_loc=zero_25*(dl1*hwll(l,1,ivar)+dl2*hwll(lp,1,ivar)) < scale3(i,j)=one < end if < < scale1(i,j)=hwll_loc/dyf(iglob,jglob) < scale2(i,j)=hwll_loc/dxf(iglob,jglob) < < < if (twodvar_regional) then < < !Rescaling < < < if (jdvar(k)==4 .or. jdvar(k)==5) then < scale1(i,j)=0.3_r_kind*scale1(i,j) !0.6! !0.4! 0.75 < scale2(i,j)=0.3_r_kind*scale2(i,j) !0.6! !0.4! 0.75 < else < scale1(i,j)=0.3_r_kind*scale1(i,j)*75._r_kind/100._r_kind !0.4! !0.5 < scale2(i,j)=0.3_r_kind*scale2(i,j)*75._r_kind/100._r_kind !0.4! !0.5 < endif < < else < < if (.not.latdepend) then < scale1(i,j)=max(scale1(i,j),scale2(i,j)) < scale2(i,j)=scale1(i,j) < endif < < !rescaling to roughly match original analysis from purely isotropic < !option, ie.. when anisotropic=.false. in namelist "anbkgerr". < < scale1(i,j)=rfact0h(jdvar(k))*scale1(i,j) < scale2(i,j)=rfact0h(jdvar(k))*scale2(i,j) < if (jdvar(k) /=3 .and. jdvar(k) /=7 .and. & < jdvar(k) /=9 .and. jdvar(k) /=10 ) & < scale3(i,j)=rfact0v(jdvar(k))*scale3(i,j) < endif < enddo < enddo < return < end subroutine isotropic_scales_subdomain_option < DIFFERENCES FOR anisofilter.f90_15Jan2008 DIFFERENCES FOR anisofilter.f90_no_changes DIFFERENCES FOR antest_maps0.f90 23d22 < character(80) plotname 47c46 < i_plotcor=460 !440! 430 --- > i_plotcor=440! 430 66,67c65,66 < ! open(lunin,file='cormaps_'//trim(var_plotcor),form='unformatted') < ! rewind lunin --- > open(lunin,file='cormaps_'//trim(var_plotcor),form='unformatted') > rewind lunin 85,86c84,85 < !if(mype.eq.0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, & < ! nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1 --- > if(mype.eq.0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, & > nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1 112,113c111 < !if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'theta') --- > if(mype.eq.0) write(lunin) outwork0 134,135c132 < ! if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'sm_theta') --- > if(mype.eq.0) write(lunin) outwork0 151,155c148 < ! if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) then < write(plotname,'("slab_",a)')trim(var_plotcor) < call outgrads1(outwork0,nlon,nlat,plotname) < end if --- > if(mype.eq.0) write(lunin) outwork0 170,171c163 < ! if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'z') --- > if(mype.eq.0) write(lunin) outwork0 192,193c184 < ! if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'sm_z') --- > if(mype.eq.0) write(lunin) outwork0 196c187 < !close(lunin) --- > close(lunin) 197a189,192 > ! if(mype.gt.-1000) then > ! call mpi_finalize(i) > ! stop > ! end if 202,479d196 < subroutine antest_maps0_subdomain_option(mype,theta0f,z0f) < < ! this routine creates output maps of background error correlations < < use kinds, only: r_kind,i_kind,r_single < use anberror, only: nvars,kvar_start,kvar_end,var_names,kps,kpe < use gridmod, only: nsig,nlon,nlat,istart,jstart,lat2,lon2 < use constants, only: zero,one,rd_over_cp < use mpimod, only: npe,ierror,mpi_integer4,mpi_real4,mpi_real8,mpi_max,mpi_min,mpi_sum,mpi_comm_world < use guess_grids, only: ges_tv,ges_z,ntguessig,ges_prsl < use fgrid2agrid_mod, only: nlatf,nlonf,fgrid2agrid < implicit none < < integer(i_kind) mype < real(r_single) theta0f(lat2,lon2,nsig) < real(r_single) z0f(lat2,lon2,nsig) < < real(r_kind),dimension(lat2,lon2,nsig):: twork,qwork,stwork,vpwork < real(r_kind),dimension(lat2,lon2):: pwork < real(r_kind) tempf(nlat,nlon),tempc(nlatf,nlonf) < real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat) < character(80) ref_plotcor < character(80) var_plotcor < character(80) plotname < integer(i_kind) i_plotcor,j_plotcor,k_plotcor < integer(i_kind) iloc,jloc,kloc < < real(r_kind)h00,h000 < integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot,k_plot < integer(i_kind) it,mm1 < real(r_kind),parameter:: r100=100.0_r_kind < integer(i_kind) lvar < < !********************************************************************* < ! variable names expected for var_plotcor are < ! < ! st -- stream function < ! vp -- velocity potential < ! ps -- surface pressure < ! tv -- virtual temperature < ! q -- specific humidity < ! oz -- ozone < ! sst -- sea surface temperature < ! stl -- skin temp over land < ! sti -- skin temp over ice < ! cw -- cloud water < !********************************************************************* < ! Make choices here! < i_plotcor=460 !440! 430 < j_plotcor=265! 250 < k_plotcor=1 < iloc=i_plotcor-istart(mype+1)+2 < jloc=j_plotcor-jstart(mype+1)+2 < ! var_plotcor='st' < !Note: Must call this subroutine from anprewgt_reg.f90 < !Make sure statement has been uncommented! < ! End of choice section < !********************************************************************* < < ref_plotcor='theta' < it=ntguessig < do 200 lvar=1,5 < if (lvar.eq.1) var_plotcor='st' < if (lvar.eq.2) var_plotcor='vp' < if (lvar.eq.3) var_plotcor='ps' < if (lvar.eq.4) var_plotcor='tv' < if (lvar.eq.5) var_plotcor='q' < lunin=1 < if(mype.eq.0) then < ! open(lunin,file='cormaps_'//trim(var_plotcor),form='unformatted') < ! rewind lunin < end if < ivar_plot=0 < do ivar=1,nvars < if(trim(var_names(ivar)).eq.trim(var_plotcor)) then < ivar_plot=ivar < exit < end if < end do < if(ivar_plot.eq.0) then < write(0,*)' in antest_maps0, variable ',trim(var_plotcor),' not found. program stops' < call mpi_finalize(ierror) < stop < end if < kloc=k_plotcor < twork=zero < pwork=zero < qwork=zero < stwork=zero < vpwork=zero < if(i_plotcor.ge.ips.and.i_plotcor.le.ipe.and. & < j_plotcor.ge.jps.and.j_plotcor.le.jpe) then < if(ivar_plot.eq.1) stwork(iloc,jloc,kloc)=one < if(ivar_plot.eq.2) vpwork(iloc,jloc,kloc)=one < if(ivar_plot.eq.3) pwork(iloc,jloc)=one < if(ivar_plot.eq.4) twork(iloc,jloc,kloc)=one < if(ivar_plot.eq.5) qwork(iloc,jloc,kloc)=one < end if < < < call ansmoothrf_reg_subdomain_option(twork,pwork,qwork,stwork,vpwork,mype) < !if(mype.eq.0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, & < ! nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1 < if(mype.eq.0) write(0,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), & < ' i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', & < nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1 < < !---------------------in case we haven't normalized, divide by value of correlation point < h00=zero < if(i_plotcor.ge.ips.and.i_plotcor.le.ipe.and. & < j_plotcor.ge.jps.and.j_plotcor.le.jpe) then < if(ivar_plot.eq.1) h00=stwork(iloc,jloc,kloc) < if(ivar_plot.eq.2) h00=vpwork(iloc,jloc,kloc) < if(ivar_plot.eq.3) h00=pwork(iloc,jloc) < if(ivar_plot.eq.4) h00=twork(iloc,jloc,kloc) < if(ivar_plot.eq.5) h00=qwork(iloc,jloc,kloc) < end if < call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror) < stwork=stwork/h000 < vpwork=vpwork/h000 < pwork=pwork/h000 < twork=twork/h000 < qwork=qwork/h000 < < < ! output original pot temp (slow way to get full 2d field) -- this is reference field < < mm1=mype+1 < do k=1,nsig < outwork=0._r_single < do j=2,lon2-1 < jglob=jstart(mm1)-2+j < do i=2,lat2-1 < iglob=istart(mm1)-2+i < outwork(jglob,iglob)=ges_tv(i,j,k,it)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp < end do < end do < call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < !if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'theta') < end do < < ! output "smoothed pot temp" < < do k=1,nsig < outwork=0._r_single < do j=2,lon2-1 < jglob=jstart(mm1)-2+j < do i=2,lat2-1 < iglob=istart(mm1)-2+i < outwork(jglob,iglob)=theta0f(i,j,k) < end do < end do < call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < !if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'sm_theta') < end do < < do k=kvar_start(ivar_plot),kvar_end(ivar_plot) < kloc=levs_jdvar(k) < outwork=0._r_single < do j=2,lon2-1 < jglob=jstart(mm1)-2+j < do i=2,lat2-1 < iglob=istart(mm1)-2+i < if(ivar_plot.eq.1) outwork(jglob,iglob)=stwork(i,j,kloc) < if(ivar_plot.eq.2) outwork(jglob,iglob)=vpwork(i,j,kloc) < if(ivar_plot.eq.3) outwork(jglob,iglob)=pwork(i,j) < if(ivar_plot.eq.4) outwork(jglob,iglob)=twork(i,j,kloc) < if(ivar_plot.eq.5) outwork(jglob,iglob)=qwork(i,j,kloc) < end do < end do < ! very slow way to move field from local processor to processor 0 < < call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < !if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) then < write(plotname,'("sub_",a)')trim(var_plotcor) < call outgrads1(outwork0,nlon,nlat,plotname) < end if < end do < < ! output original terrain (slow way to get full 2d field) -- this is reference field < < mm1=mype+1 < outwork=0._r_single < do j=2,lon2-1 < jglob=jstart(mm1)-2+j < do i=2,lat2-1 < iglob=istart(mm1)-2+i < outwork(jglob,iglob)=ges_z(i,j,it) < end do < end do < call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < !if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'z') < < ! output "smoothed terrain" < < PRINT*,'IN ANPREWGT_REG_subdomain_option,KPS,KPE=',KPS,KPE < < do k=1,1 < outwork=0._r_single < do j=2,lon2-1 < jglob=jstart(mm1)-2+j < do i=2,lat2-1 < iglob=istart(mm1)-2+i < outwork(jglob,iglob)=z0f(i,j,k) < end do < end do < call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror) < ! if(mype.eq.0) write(lunin) outwork0 < if(mype.eq.0) call outgrads1(outwork0,nlon,nlat,'sm_z') < end do < < !close(lunin) < 200 continue < < end subroutine antest_maps0_subdomain_option < < !------------------------------------------------------------------------------------- < < subroutine outgrads1(f,nx,ny,label) < < character(*) label < integer(4) nx,ny < real(4) f(nx,ny) < < integer(4) i,l,next,last,np,ntime,ioutdat,ioutcor,koutmax < real(4) rlonmap0,undef,dlonmap,pinc,startp,rlatmap0,dlatmap < character(80) dsdes,dsdat < character(80) datdes(1000) < character(1) blank < data blank/' '/ < data undef/-9.99e33/ < < ioutcor=10 < ioutdat=11 < < write(dsdes,'(a,".des")')trim(label) < write(dsdat,'(a,".dat")')trim(label) < open(unit=ioutcor,file=dsdes,form='formatted') < open(unit=ioutdat,file=dsdat,form='unformatted') < ntime=1 < rlonmap0=1. < dlonmap=1. < rlatmap0=1. < dlatmap=1. < startp=1. < pinc=1. < koutmax=1 < do i=1,1000 < write(datdes(i),'(80a1)')(blank,l=1,80) < end do < write(datdes(1),'("DSET ",a)')trim(dsdat) < write(datdes(2),'("options big_endian sequential")') < write(datdes(3),'("TITLE ",a)')trim(label) < write(datdes(4),'("UNDEF ",e11.2)')undef < write(datdes(5),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nx,rlonmap0,dlonmap < write(datdes(6),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')ny,rlatmap0,dlatmap < next=7 < write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc < next=next+1 < write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax < next=next+1 < write(datdes(next),'("VARS 1")') < next=next+1 < write(datdes(next),'("f ",i5," 99 f ")')np < next=next+1 < write(datdes(next),'("ENDVARS")') < last=next < write(ioutcor,'(a80)')(datdes(i),i=1,last) < close(ioutcor) < < write(ioutdat) f < close(ioutdat) < < end subroutine outgrads1 DIFFERENCES FOR antest_maps0.f90_single_variable DIFFERENCES FOR b_l.f DIFFERENCES FOR balmod.f90 DIFFERENCES FOR balmod.f90_no_changes DIFFERENCES FOR berror.f90 DIFFERENCES FOR bkerror.f90 DIFFERENCES FOR bkgcov.f90 DIFFERENCES FOR bkgvar.f90 DIFFERENCES FOR bkgvar_rewgt.f90 DIFFERENCES FOR calctends.f90 DIFFERENCES FOR calctends_ad.f90 DIFFERENCES FOR calctends_tl.f90 DIFFERENCES FOR combine_radobs.f90 DIFFERENCES FOR compact_diffs.f90 DIFFERENCES FOR compute_derived.f90 DIFFERENCES FOR compute_fact10.f90 DIFFERENCES FOR constants.f90 DIFFERENCES FOR converr.f90 DIFFERENCES FOR convinfo.f90 DIFFERENCES FOR convthin.f90 DIFFERENCES FOR deter_subdomain.f90 DIFFERENCES FOR dprodx.F90 DIFFERENCES FOR dtast.f90 DIFFERENCES FOR fgrid2agrid_mod.f90 DIFFERENCES FOR fill_mass_grid2.f90 DIFFERENCES FOR fill_nmm_grid2.f90 DIFFERENCES FOR fpvsx_ad.f90 DIFFERENCES FOR gengrid_vars.f90 DIFFERENCES FOR genqsat.f90 DIFFERENCES FOR genstats_gps.f90 DIFFERENCES FOR gesinfo.f90 DIFFERENCES FOR get_derivatives.f90 DIFFERENCES FOR get_semimp_mats.f90 DIFFERENCES FOR get_streamlined_stats_2009.f DIFFERENCES FOR get_streamlined_stats_2009_new.f DIFFERENCES FOR get_tend_derivs.f90 DIFFERENCES FOR getprs.f90 DIFFERENCES FOR getstvp.f90 DIFFERENCES FOR getuv.f90 DIFFERENCES FOR getvvel.f90 DIFFERENCES FOR glbsoi.f90 DIFFERENCES FOR grdcrd.f90 DIFFERENCES FOR grid2sub.f90 DIFFERENCES FOR gridmod.f90 DIFFERENCES FOR gscond_ad.f90 DIFFERENCES FOR gsi_io.f90 DIFFERENCES FOR gsimain.F90 346,347c346 < grid_ratio,an_flen_u,an_flen_t,an_flen_z,& < rtma_subdomain_option,nsmooth,nsmooth_shapiro --- > grid_ratio,an_flen_u,an_flen_t,an_flen_z 538,548d536 < ! rtma_subdomain_option - if true, then call alternative code which calls recursive filter < ! directly from subdomain mode, bypassing transition to/from < ! horizontal slabs. This is mainly to improve efficiency for < ! 2d rtma analysis. at the moment, this only works for < ! twodvar_regional=.true. rtma_subdomain_option will be forced < ! to false when twodvar_regional=.false. < ! nsmooth - number of 1-2-1 smoothing passes before and after background error application < ! nsmooth_shapiro - number of 2nd moment preserving (shapiro) smoothing passes before and after < ! background error application. < ! NOTE: default for nsmooth and nsmooth_shapiro is 0. < ! if both are > 0, then nsmooth will be forced to zero. 552,553c540 < grid_ratio,nord_f2a,an_flen_u,an_flen_t,an_flen_z,rtma_subdomain_option, & < nsmooth,nsmooth_shapiro --- > grid_ratio,nord_f2a,an_flen_u,an_flen_t,an_flen_z 731,733d717 < < ! make sure rtma_subdomain_option turned on only for 2d analysis < rtma_subdomain_option=rtma_subdomain_option.and.twodvar_regional DIFFERENCES FOR gsisub.f90 DIFFERENCES FOR guess_grids.f90 DIFFERENCES FOR half_nmm_grid2.f90 DIFFERENCES FOR hilbert_curve.f90 DIFFERENCES FOR init_commvars.f90 DIFFERENCES FOR init_var_tl.f90 DIFFERENCES FOR intall.f90 DIFFERENCES FOR intdw.f90 DIFFERENCES FOR intgps.f90 DIFFERENCES FOR intjc.f90 DIFFERENCES FOR intlimq.f90 DIFFERENCES FOR intoz.f90 DIFFERENCES FOR intpcp.f90 DIFFERENCES FOR intps.f90 DIFFERENCES FOR intpw.f90 DIFFERENCES FOR intq.f90 DIFFERENCES FOR intrad.f90 DIFFERENCES FOR intrp2a.f90 DIFFERENCES FOR intrp3oz.f90 DIFFERENCES FOR intrppx.f90 DIFFERENCES FOR intrw.f90 DIFFERENCES FOR intspd.f90 DIFFERENCES FOR intsrw.f90 DIFFERENCES FOR intsst.f90 DIFFERENCES FOR intt.f90 DIFFERENCES FOR intw.f90 DIFFERENCES FOR jcdivtends.f90 DIFFERENCES FOR jcmod.f90 DIFFERENCES FOR jfunc.f90 DIFFERENCES FOR jfunc_tl.f90 DIFFERENCES FOR kinds.f90 DIFFERENCES FOR lagmod.f90 DIFFERENCES FOR linbal.f90 DIFFERENCES FOR m_fvAnaGrid.F90 DIFFERENCES FOR m_gsiCheck.F90 DIFFERENCES FOR m_gsiGuess.F90 DIFFERENCES FOR m_utests.F90 DIFFERENCES FOR mod_inmi.f90 DIFFERENCES FOR mod_strong.f90 DIFFERENCES FOR mod_vtrans.f90 DIFFERENCES FOR mp_compact_diffs_mod1.f90 DIFFERENCES FOR mp_compact_diffs_support.f90 DIFFERENCES FOR mpi_bufr_mod.F90 DIFFERENCES FOR mpimod.F90 DIFFERENCES FOR ncepgfs_io.f90 DIFFERENCES FOR nlmsas_ad.f90 DIFFERENCES FOR normal_rh_to_q.f90 DIFFERENCES FOR obs_para.f90 DIFFERENCES FOR obsmod.F90 DIFFERENCES FOR obsmod_tl.f90 DIFFERENCES FOR omegas_ad.f90 DIFFERENCES FOR oneobmod.F90 DIFFERENCES FOR ozinfo.f90 DIFFERENCES FOR pcgsoi.f90 DIFFERENCES FOR pcgsoi.f90_no_changes DIFFERENCES FOR pcp_k.f90 DIFFERENCES FOR pcpinfo.f90 DIFFERENCES FOR penal.f90 DIFFERENCES FOR phil.f90 DIFFERENCES FOR phil1.f90 DIFFERENCES FOR plib8.f90 DIFFERENCES FOR polcarf.f90 DIFFERENCES FOR precpd_ad.f90 DIFFERENCES FOR prepbykx.f DIFFERENCES FOR prewgt.f90 DIFFERENCES FOR prewgt_reg.f90 DIFFERENCES FOR projmethod_support.f90 102d101 < real(8) dplev_mask 112c111 < prd0=dplev_mask(gy(1,iter),gx(1,i),mype) --- > prd0=dplev(gy(1,iter),gx(1,i),mype) 117c116 < prd0=dplev_mask(gx(1,iter),gy(1,iter),mype) --- > prd0=dplev(gx(1,iter),gy(1,iter),mype) 125c124 < prd0=dplev_mask(gradx,grady,mype) --- > prd0=dplev(gradx,grady,mype) 279,441d277 < < real(8) function dplev_mask(dx,dy,mype) < < use kinds, only: r_kind,i_kind < use jfunc, only: nval_levs < use gridmod, only: lat2,lon2,twodvar_regional < implicit none < < ! Declar passed variables < real(r_kind),dimension(lat2,lon2,nval_levs),intent(in)::dx,dy < integer(i_kind),intent(in)::mype < < ! Declare local variables < logical mask(nval_levs) < logical fast < real(8) fast_dplev,dplev5 < < fast=.false. < mask=.true. < ! set fast to .true. for twodvar_regional, < ! substantially faster, but no roundoff error reduction and < ! results differ for different number of processors. < if(twodvar_regional) then < fast=.true. < mask(5)=.false. < mask(6)=.false. < mask(8)=.false. < end if < < if(fast) then < dplev_mask=fast_dplev(dx,dy,mype,mask) < else < dplev_mask=dplev5(dx,dy,mype,mask) < end if < < end function dplev_mask < < real(8) function fast_dplev(dx,dy,mype,mask) < < use kinds, only: r_kind,i_kind < use jfunc, only: nval_levs < use gridmod, only: lat2,lon2 < use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype < use constants, only: zero < implicit none < < ! Declar passed variables < real(r_kind),dimension(lat2,lon2,nval_levs),intent(in)::dx,dy < integer(i_kind),intent(in)::mype < logical,intent(in):: mask(nval_levs) < < ! Declare local variables < real(r_kind),dimension(npe):: sumall < < integer(i_kind) i,j,k < real(r_kind) sum < < sum=zero < do k=1,nval_levs < if(.not.mask(k)) cycle < do j=2,lon2-1 < do i=2,lat2-1 < sum=sum+dx(i,j,k)*dy(i,j,k) < end do < end do < end do < < call mpi_allgather(sum,1,mpi_rtype,sumall,1,mpi_rtype,mpi_comm_world,ierror) < fast_dplev=zero < do i=1,npe < fast_dplev=fast_dplev+sumall(i) < end do < < end function fast_dplev < < real(8) function dplev5(dx,dy,mype,mask) < < !$$$ subprogram documentation block < ! . . . . < ! subprogram: dplev 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 jfunc, only: nval_levs < use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,& < ltosi,ltosj,iglobal,itotsub,ijn,displs_g < use mpimod, only: mpi_comm_world,ierror,mpi_rtype,strip < use constants, only: zero < implicit none < < ! Declar passed variables < real(r_kind),dimension(lat2,lon2,nval_levs),intent(in)::dx,dy < integer(i_kind),intent(in)::mype < logical,intent(in):: mask(nval_levs) < < ! Declare local variables < real(r_kind),dimension(lat1*lon1):: zsm < real(r_kind),dimension(iglobal):: work1 < real(r_kind),dimension(lat2,lon2):: sum < real(r_kind),dimension(nlat,nlon):: sumall < < integer(i_kind) i,j,k,mm1,kk < real(r_kind) e,y,temp < < mm1=mype+1 < sum=zero < do k=1,nval_levs < if(.not.mask(k)) cycle < do j=1,lon2 < do i=1,lat2 < sum(i,j)=sum(i,j)+dx(i,j,k)*dy(i,j,k) < end do < end do < end do < do j=1,lon1*lat1 < zsm(j)=zero < end do < < call strip(sum,zsm,1) < < call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& < work1,ijn,displs_g,mpi_rtype,& < mpi_comm_world,ierror) < < do k=1,iglobal < i=ltosi(k) ; j=ltosj(k) < sumall(i,j)=work1(k) < end do < dplev5=zero < e=zero < do j=1,nlon < do i=1,nlat < ! Compensated summation version of sum < temp=dplev5 < y=sumall(i,j)+e < dplev5=temp+y < e=(temp-dplev5)+y < ! dplev=dplev+sumall(i,j) < end do < end do < < end function dplev5 DIFFERENCES FOR psichi2uv_reg.f90 DIFFERENCES FOR psichi2uvt_reg.f90 DIFFERENCES FOR q_diag.f90 DIFFERENCES FOR qcmod.f90 DIFFERENCES FOR qcssmi.f90 DIFFERENCES FOR radinfo.f90 DIFFERENCES FOR raflib.f90 50c50 < use mpimod,only: mpi_comm_world,mpi_integer,mpi_integer1,mpi_integer2,mpi_integer4, & --- > use mpimod,only: mpi_comm_world,mpi_integer,mpi_integer2,mpi_integer4, & 59d58 < logical new_factorization 68,102c67,78 < integer(i_long) nsmooth < integer(i_long) nsmooth_shapiro < integer(i_long) nrows < integer(i_long) nsend_halox_loc < integer(i_long) nrecv_halox_loc < integer(i_long) nsend_haloy_loc < integer(i_long) nrecv_haloy_loc < integer(i_long),pointer:: nrecv_halox(:) => NULL() < integer(i_long),pointer:: ndrecv_halox(:) => NULL() < integer(i_long),pointer:: nsend_halox(:) => NULL() < integer(i_long),pointer:: ndsend_halox(:) => NULL() < integer(i_long),pointer:: info_send_halox(:,:) => NULL() < integer(i_long),pointer:: info_recv_halox(:,:) => NULL() < integer(i_long),pointer:: nrecv_haloy(:) => NULL() < integer(i_long),pointer:: ndrecv_haloy(:) => NULL() < integer(i_long),pointer:: nsend_haloy(:) => NULL() < integer(i_long),pointer:: ndsend_haloy(:) => NULL() < integer(i_long),pointer:: info_send_haloy(:,:) => NULL() < integer(i_long),pointer:: info_recv_haloy(:,:) => NULL() < integer(i_long),pointer::istart(:) => NULL() < integer(i_long),pointer::ib(:) => NULL() < integer(i_long),pointer::nrecv(:) => NULL() < integer(i_long),pointer::ndrecv(:) => NULL() < integer(i_long),pointer::nsend(:) => NULL() < integer(i_long),pointer::ndsend(:) => NULL() < real(r_single),pointer:: gnorm_halox(:) => NULL() < real(r_single),pointer:: gnorm_haloy(:) => NULL() < real(r_single),pointer::lnf(:,:,:,:) => NULL() < real(r_single),pointer::bnf(:,:,:) => NULL() < real(r_single),pointer::fmat(:,:,:,:,:) => NULL() < real(r_single),pointer::fmat0(:,:,:,:) => NULL() < real(r_single),pointer::amp(:,:,:,:) => NULL() < integer(i_short),pointer::ia(:) => NULL() < integer(i_short),pointer::ja(:) => NULL() < integer(i_short),pointer::ka(:) => NULL() --- > integer(i_long),pointer::istart(:) > integer(i_long),pointer::ib(:) > integer(i_long),pointer::nrecv(:) > integer(i_long),pointer::ndrecv(:) > integer(i_long),pointer::nsend(:) > integer(i_long),pointer::ndsend(:) > real(r_single),pointer::lnf(:,:,:,:) > real(r_single),pointer::bnf(:,:,:) > real(r_single),pointer::amp(:,:,:,:) > integer(i_short),pointer::ia(:) > integer(i_short),pointer::ja(:) > integer(i_short),pointer::ka(:) 106d81 < 109c84,88 < subroutine adjoint_check4(filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > subroutine adjoint_check4(filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 114,115c93,95 < INTEGER(i_long), INTENT(IN) :: ips,ipe,jps,jpe,kps,kpe < INTEGER(i_long), INTENT(IN) :: mype,npes,ngauss --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 116a97,99 > INTEGER(i_long), INTENT(IN) :: & > mype, npes,ngauss > 119,121c102,104 < real(r_single) xvec( ngauss,ips:ipe, jps:jpe, kps:kpe ) < real(r_single) yvec( ngauss,ips:ipe, jps:jpe, kps:kpe ) < real(r_single) zvec( ngauss,ips:ipe, jps:jpe, kps:kpe ) --- > real(r_single) xvec( ngauss,ims:ime, jms:jme, kms:kme ) > real(r_single) yvec( ngauss,ims:ime, jms:jme, kms:kme ) > real(r_single) zvec( ngauss,ims:ime, jms:jme, kms:kme ) 134c117,121 < call raf_sm4(yvec,filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > call raf_sm4(yvec,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 136c123,127 < call raf_sm4_ad(zvec,filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > call raf_sm4_ad(zvec,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 164c155,159 < SUBROUTINE raf4_ad(g,filter,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > SUBROUTINE raf4_ad(g,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 170c165,167 < INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, ips, ipe, jps, jpe, kps, kpe --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 172c169,170 < INTEGER(i_long), INTENT(IN) :: ngauss,mype, npes --- > INTEGER(i_long), INTENT(IN) :: & > ngauss,mype, npes 174c172 < real(r_single), DIMENSION( ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION( ngauss,ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 179c177 < integer(i_long) i,icolor,igauss,ipass,j,k,iadvance,iback --- > integer(i_long) i,icolor,igauss,ipass,j,k 180a179 > real(r_single), DIMENSION( ims:ime, jms:jme ) :: gwrk 181a181,184 > nsmooth=0 ! hardwired values for > nsmooth_shapiro=1 !testing purposes. Will change to namelist variables > > 187,194c190 < if(filter(icolor)%npointsmaxall.gt.0) then < if(filter(1)%new_factorization) then < iadvance=2 < iback=1 < call one_color4_new_factorization(g,filter(icolor),ngauss,ipass,filter(1)%ifilt_ord, & < filter(icolor)%nstrings,filter(icolor)%istart, & < ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < else --- > if(filter(icolor)%npointsmaxall.gt.0) & 196,198c192,196 < filter(icolor)%nstrings,filter(icolor)%istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) < end if < end if --- > filter(icolor)%nstrings,filter(icolor)%istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 205,212c203,211 < ! apply smoothing if nsmooth or nsmooth_shapiro > 0 < if(kps.le.kpe) then < do i=1,max(filter(1)%nsmooth,filter(1)%nsmooth_shapiro) < call smther(filter(1),g,filter(1)%nrows,ngauss,filter(1)%nsmooth,filter(1)%nsmooth_shapiro, & < ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe, & < mype,npes,filter(1)%gnorm_halox,filter(1)%gnorm_haloy,.true.) < end do < end if --- > if (nsmooth > 0 ) then > do igauss=1,ngauss > do k=kms,kme > gwrk=g(igauss,ims:ime,jms:jme,k) > call smther_one(gwrk,ims,ime,jms,jme,nsmooth)! self_adjoint > g(igauss,ims:ime,jms:jme,k)=gwrk > enddo > enddo > endif 213a213,222 > if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then > do igauss=1,ngauss > do k=kms,kme > gwrk=g(igauss,ims:ime,jms:jme,k) > call tsmther_two(gwrk,ims,ime,jms,jme,nsmooth_shapiro) > g(igauss,ims:ime,jms:jme,k)=gwrk > enddo > enddo > endif > 228c237,241 < SUBROUTINE raf_sm4_ad(g,filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > SUBROUTINE raf_sm4_ad(g,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 235c248,250 < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 237c252,253 < INTEGER(i_long), INTENT(IN) :: mype, npes,ngauss --- > INTEGER(i_long), INTENT(IN) :: & > mype, npes,ngauss 239c255 < real(r_single), DIMENSION( ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION( ngauss,ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 244c260 < integer(i_long) icolor,ipass,iadvance,iback --- > integer(i_long) icolor,ipass 251,258c267 < if(filter(icolor)%npointsmaxall.gt.0) then < if(filter(1)%new_factorization) then < iadvance=2 < iback=1 < call one_color4_new_factorization(g,filter(icolor),ngauss,ipass,filter(1)%ifilt_ord, & < filter(icolor)%nstrings,filter(icolor)%istart, & < ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < else --- > if(filter(icolor)%npointsmaxall.gt.0) & 260,262c269,273 < filter(icolor)%nstrings,filter(icolor)%istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) < end if < end if --- > filter(icolor)%nstrings,filter(icolor)%istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 271c282,286 < SUBROUTINE rad_sm24_ad(g,filter,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > SUBROUTINE rad_sm24_ad(g,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 279c294,296 < INTEGER(i_long), INTENT(IN) :: ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 281c298,299 < INTEGER(i_long), INTENT(IN) :: mype,npes,ngauss --- > INTEGER(i_long), INTENT(IN) :: & > mype, npes,ngauss 283c301 < real(r_single), DIMENSION(2,ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION(2,ngauss,ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 288,289c306 < integer(i_long) icolor,ipass,i,ii,j,k,kk,n,iadvance,iback < real(r_single) gwork(ngauss,ips:ipe,jps:jpe,kps:kpe) --- > integer(i_long) icolor,ipass 296,303c313 < if(filter(icolor)%npointsmaxall.gt.0) then < if(filter(1)%new_factorization) then < iadvance=2 < iback=1 < call one_color24_new_factorization(g,filter(icolor),ngauss,ipass,filter(1)%ifilt_ord, & < filter(icolor)%nstrings,filter(icolor)%istart, & < ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < else --- > if(filter(icolor)%npointsmaxall.gt.0) & 305,308c315,319 < filter(icolor)%nstrings,filter(icolor)%istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) < end if < end if < --- > filter(icolor)%nstrings,filter(icolor)%istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 315,343d325 < ! apply smoothing if nsmooth or nsmooth_shapiro > 0 < if(max(filter(1)%nsmooth,filter(1)%nsmooth_shapiro) > 0.and.kps.le.kpe) then < do kk=1,2 < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < gwork(n,i,j,k)=g(kk,n,i,j,k) < end do < end do < end do < end do < do ii=1,max(filter(1)%nsmooth,filter(1)%nsmooth_shapiro) < call smther(filter(1),gwork,filter(1)%nrows,ngauss,filter(1)%nsmooth,filter(1)%nsmooth_shapiro, & < ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe, & < mype,npes,filter(1)%gnorm_halox,filter(1)%gnorm_haloy,.true.) < end do < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(kk,n,i,j,k)=gwork(n,i,j,k) < end do < end do < end do < end do < end do < end if < 348c330,334 < lenbar,lenmax,lenmin,npoints1,mype,npes,nvars) --- > lenbar,lenmax,lenmin,npoints1, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes,nvars) 355c341,345 < INTEGER(i_long), INTENT(IN) :: mype,npes,nvars --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices > INTEGER(i_long), INTENT(IN) :: & > mype, npes,nvars 378d367 < ! if(ifilt_ord.lt.0) write(0,*)' in alpha_beta4, ifilt_ord,npass=',ifilt_ord,npass 724c713 < write(6,*)'GETHEX: ALL ITERATIONS USED UP. This should never happen' --- > write(6,*)'INDEXXI4: ALL ITERATIONS USED UP. This should never happen' 858d846 < nsmooth,nsmooth_shapiro, & 861a850 > ims, ime, jms, jme, kms, kme, & ! memory indices 872c861,862 < ips, ipe, jps, jpe, kps, kpe ! patch indices --- > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 890,894d879 < ! NOTE: if ifilt_ord = -4, this is the special case < ! of a new approximate 4th order factorization that < ! only requires 2nd order amount of work. < integer(i_long), intent(inout) :: nsmooth ! number of 1-2-1 smoothing passes to apply < integer(i_long), intent(in) :: nsmooth_shapiro ! number of 2nd moment preserving shapiro smoothers 951d935 < integer(i_long) istat 953,973d936 < ! nsmooth_shapiro has priority over nsmooth < if(nsmooth_shapiro.gt.0) nsmooth=0 < ! set up halo update (only if nsmooth and/or nsmooth_shapiro > 0) < filter(1)%nrows=0 < if(nsmooth.gt.0) filter(1)%nrows=1 < if(nsmooth_shapiro.gt.0) filter(1)%nrows=3 < filter(1)%nsmooth=nsmooth < filter(1)%nsmooth_shapiro=nsmooth_shapiro < call add_halox0(filter(1),filter(1)%nrows,ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes) < call add_haloy0(filter(1),filter(1)%nrows,ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes) < deallocate(filter(1)%gnorm_halox,stat=istat) < allocate(filter(1)%gnorm_halox(ips:ipe)) < filter(1)%gnorm_halox=1. < deallocate(filter(1)%gnorm_haloy,stat=istat) < allocate(filter(1)%gnorm_haloy(jps:jpe)) < filter(1)%gnorm_haloy=1. < if(nsmooth_shapiro.gt.0) then < call smther_two_gnorm(filter(1)%gnorm_halox,ids,ide,ips,ipe) < call smther_two_gnorm(filter(1)%gnorm_haloy,jds,jde,jps,jpe) < end if < 976,983c939 < write(0,*)' at 1 in init_raf4, ifilt_ord,npass=',ifilt_ord,npass < if(ifilt_ord.eq.-4) then < filter(1)%ifilt_ord=2 < filter(1)%new_factorization=.true. < else < filter(1)%ifilt_ord=ifilt_ord < filter(1)%new_factorization=.false. < end if --- > filter(1)%ifilt_ord=ifilt_ord 1243c1199 < mype,npes,icolor) --- > ims, ime, jms, jme, kms, kme,mype,npes,icolor) ! memory indices 1259,1265c1215,1217 < deallocate(filter(icolor)%nsend,stat=istat) < deallocate(filter(icolor)%ndsend,stat=istat) < deallocate(filter(icolor)%nrecv,stat=istat) < deallocate(filter(icolor)%ndrecv,stat=istat) < deallocate(filter(icolor)%ia,stat=istat) < deallocate(filter(icolor)%ja,stat=istat) < deallocate(filter(icolor)%ka,stat=istat) --- > deallocate(filter(icolor)%nsend,filter(icolor)%ndsend) > deallocate(filter(icolor)%nrecv,filter(icolor)%ndrecv) > deallocate(filter(icolor)%ia,filter(icolor)%ja,filter(icolor)%ka) 1279c1231,1234 < ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1283c1238 < deallocate(filter(icolor)%ib,stat=istat) --- > deallocate(filter(icolor)%ib) 1286c1241,1246 < call sort_strings4(info_string,aspect_full,npoints_recv(mype),filter(icolor)%ib,mype,npes,nvars) --- > call sort_strings4(info_string,aspect_full, & > npoints_recv(mype),filter(icolor)%ib, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1290c1250 < call count_strings(info_string,filter(icolor)%nstrings,nstrings_var,nvars,npoints_recv(mype)) --- > call count_strings(info_string,filter(icolor)%nstrings,nstrings_var,nvars,npoints_recv(mype)) 1294,1295c1254 < deallocate(filter(icolor)%istart,stat=istat) < allocate(filter(icolor)%istart(filter(icolor)%nstrings+1)) --- > deallocate(filter(icolor)%istart,filter(icolor)%lnf,filter(icolor)%bnf) 1297,1306c1256,1269 < if(filter(1)%new_factorization) then < deallocate(filter(icolor)%fmat,stat=istat) < deallocate(filter(icolor)%fmat0,stat=istat) < allocate(filter(icolor)%fmat(2,npoints_recv(mype),npass,ngauss,2)) < allocate(filter(icolor)%fmat0(npoints_recv(mype),npass,ngauss,2)) < do igauss=1,ngauss < call new_alpha_beta4(info_string,aspect_full,rgauss(igauss), & < filter(icolor)%fmat,filter(icolor)%fmat0,igauss,ngauss, & < filter(icolor)%istart,npoints_recv(mype),binomial,npass, & < lenbar,lenmax,lenmin,npoints1,mype,npes,nvars) --- > allocate(filter(icolor)%istart(filter(icolor)%nstrings+1)) > allocate(filter(icolor)%lnf(ifilt_ord,max(1,npoints_recv(mype)),npass,ngauss)) > allocate(filter(icolor)%bnf(max(1,npoints_recv(mype)),npass,ngauss)) > do igauss=1,ngauss > if(npoints_recv(mype).gt.0) & > call alpha_beta4(info_string,aspect_full,rgauss(igauss), & > filter(icolor)%lnf,filter(icolor)%bnf,igauss,ngauss, & > filter(icolor)%istart,npoints_recv(mype),binomial,npass, & > ifilt_ord, & > lenbar,lenmax,lenmin,npoints1, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes,nvars) 1308,1321c1271 < end do < else < deallocate(filter(icolor)%lnf,stat=istat) < deallocate(filter(icolor)%bnf,stat=istat) < allocate(filter(icolor)%lnf(ifilt_ord,npoints_recv(mype),npass,ngauss)) < allocate(filter(icolor)%bnf(npoints_recv(mype),npass,ngauss)) < do igauss=1,ngauss < call alpha_beta4(info_string,aspect_full,rgauss(igauss), & < filter(icolor)%lnf,filter(icolor)%bnf,igauss,ngauss, & < filter(icolor)%istart,npoints_recv(mype),binomial,npass, & < ifilt_ord,lenbar,lenmax,lenmin,npoints1,mype,npes,nvars) < < end do < end if --- > end do 1322a1273 > 1393,1394c1344,1349 < ! call adjoint_check4(filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) < call normalize2_raf4(filter,filter(1)%ngauss,normal, & --- > ! call adjoint_check4(filter, & > ! ids, ide, jds, jde, kds, kde, & ! domain indices > ! ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ! ims, ime, jms, jme, kms, kme, & ! memory indices > ! mype, npes) > call normalize_raf4(filter,filter(1)%ngauss,normal, & 1396a1352 > ims, ime, jms, jme, kms, kme, & ! memory indices 1405a1362 > ims, ime, jms, jme, kms, kme, & ! memory indices 1410c1367,1368 < ips, ipe, jps, jpe, kps, kpe ! patch indices --- > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 1417,1418c1375,1376 < real(r_single) ranvec(2, ngauss,ips:ipe, jps:jpe, kps:kpe ) < real(r_single) bigg( ngauss,ips:ipe, jps:jpe, kps:kpe ) --- > real(r_single) ranvec(2, ngauss,ims:ime, jms:jme, kms:kme ) > real(r_single) bigg( ngauss,ims:ime, jms:jme, kms:kme ) 1420a1379,1380 > ! real(r_double) timef > ! real(r_double) t0,t1,timerand,timetot,timerand0,timetot0 1425c1385 < real(4) work(2,ids:ide,jds:jde) --- > real(4) work(2,ips:ipe,jps:jpe) 1430,1431d1389 < integer(i_long) istat < real(8) timef,time0 1433c1391 < deallocate(filter(1)%amp,stat=istat) --- > deallocate(filter(1)%amp) 1436,1437d1393 < if(mype.eq.0) write(0,*)' entering normalize_raf4, normal=',normal < time0=timef() 1459a1416,1417 > ! timerand=0._r_double > ! t0=timef() 1461a1420 > ! t1=timef() 1477,1478d1435 < if(ips.eq.1.and.jps.eq.1.and.k.eq.1) write(0,'(" sample,ranvec=",i5,2f9.1)') & < loop,ranvec(1,1,1,1,1),ranvec(2,1,1,1,1) 1481c1438,1444 < call rad_sm24_ad(ranvec,filter,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > ! timerand=timerand+timef()-t1 > > call rad_sm24_ad(ranvec,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1491a1455,1459 > ! timetot=timef()-t0 > ! call mpi_reduce(timetot,timetot0,1,mpi_real8,mpi_max,0,mpi_comm_world,ierror) > ! call mpi_reduce(timerand,timerand0,1,mpi_real8,mpi_max,0,mpi_comm_world,ierror) > ! if(mype.eq.0) write(6,'(" total time for normalization = ",f15.3)').001_8*timetot0 > ! if(mype.eq.0) write(6,'(" total time for random_number = ",f15.3)').001_8*timerand0 1502,1506d1469 < if(mype.eq.0) write(0,*)' leaving normalize_raf4, elapsed time = ',.001*(timef()-time0) < if(mype.gt.-100) then < call mpi_finalize(i) < stop < end if 1510,1601d1472 < subroutine normalize2_raf4(filter,ngauss,normal, & < ids, ide, jds, jde, kds, kde, & ! domain indices < ips, ipe, jps, jpe, kps, kpe, & ! patch indices < mype, npes) < < < INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices < ips, ipe, jps, jpe, kps, kpe ! patch indices < < INTEGER(i_long), INTENT(IN) :: & < ngauss,normal,mype, npes < < TYPE(filter_cons) filter(7) < < real(r_single) ranvec(2, ngauss,ips:ipe, jps:jpe, kps:kpe ) < real(r_single) bigg( ngauss,ips:ipe, jps:jpe, kps:kpe ) < < integer(i_long) i,igauss,j,k,loop,nsamples,ierror < real(4) this_one1,this_one2 < < integer(i_long) istat,ii(kps:kpe),jj(kps:kpe) < integer(1) flips(2**27-8) < real(8) timef,time0 < < deallocate(filter(1)%amp,stat=istat) < allocate(filter(1)%amp(ngauss,ips:ipe,jps:jpe,kps:kpe)) < < if(normal.eq.0) then < filter(1)%amp=1. < return < end if < < ! read in fixed set of random flips < if(mype.eq.0) then < open(997433,file='random_flips',form='unformatted') < read(997433) flips < close(997433) < end if < call mpi_bcast(flips,2**27-8,mpi_integer1,0,mpi_comm_world,ierror) < < nsamples=abs(normal)/2 < < bigg=0._r_double < ii=-1 < jj=(2**27-8)*mype/npes < do loop=1,nsamples < ranvec=0._r_single < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < ii(k)=mod(ii(k)+1,8) < if(ii(k).eq.0) jj(k)=jj(k)+1 < if(jj(k).gt.2**27-8) jj(k)=1 < this_one1=-1. < if(btest(flips(jj(k)),ii(k))) this_one1=1. < ii(k)=mod(ii(k)+1,8) < if(ii(k).eq.0) jj(k)=jj(k)+1 < if(jj(k).gt.2**27-8) jj(k)=1 < this_one2=-1. < if(btest(flips(jj(k)),ii(k))) this_one2=1. < do igauss=1,ngauss < ranvec(1,igauss,i,j,k)=this_one1 < ranvec(2,igauss,i,j,k)=this_one2 < end do < end do < end do < end do < < call rad_sm24_ad(ranvec,filter,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do igauss=1,ngauss < bigg(igauss,i,j,k)=bigg(igauss,i,j,k)+ranvec(1,igauss,i,j,k)**2+ranvec(2,igauss,i,j,k)**2 < end do < end do < end do < end do < end do < < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do igauss=1,ngauss < filter(1)%amp(igauss,i,j,k)=1._r_double/sqrt(bigg(igauss,i,j,k)/(2._r_double*nsamples)) < end do < end do < end do < end do < < end subroutine normalize2_raf4 < 1603c1474,1478 < nstrings,istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > nstrings,istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1608c1483,1485 < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 1617c1494 < real(r_single), DIMENSION( ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION( ngauss,ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 1709c1586,1590 < nstrings,istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > nstrings,istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1714c1595,1597 < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 1723c1606 < real(r_single), DIMENSION(2,ngauss, ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION(2,ngauss, ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 1821c1704,1708 < SUBROUTINE raf4(g,filter,ngauss,ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > SUBROUTINE raf4(g,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1827c1714,1716 < INTEGER(i_long), INTENT(IN) :: ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 1829c1718,1719 < INTEGER(i_long), INTENT(IN) :: mype, npes,ngauss --- > INTEGER(i_long), INTENT(IN) :: & > mype, npes,ngauss 1831c1721 < real(r_single), DIMENSION(ngauss, ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION(ngauss, ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 1836c1726,1728 < integer(i_long) i,icolor,ipass,igauss,j,k,iadvance,iback --- > real(r_single), DIMENSION( ims:ime, jms:jme ) :: gwrk > integer(i_long) i,icolor,ipass,igauss,j,k > INTEGER(i_long) nsmooth,nsmooth_shapiro 1837a1730,1752 > nsmooth=0 ! hardwired values for > nsmooth_shapiro=1 !testing purposes. Will change to namelist variables > > if (nsmooth > 0 ) then > do igauss=1,ngauss > do k=kms,kme > gwrk=g(igauss,ims:ime,jms:jme,k) > call smther_one(gwrk,ims,ime,jms,jme,nsmooth) > g(igauss,ims:ime,jms:jme,k)=gwrk > enddo > enddo > endif > > if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then > do igauss=1,ngauss > do k=kms,kme > gwrk=g(igauss,ims:ime,jms:jme,k) > call smther_two(gwrk,ims,ime,jms,jme,nsmooth_shapiro) > g(igauss,ims:ime,jms:jme,k)=gwrk > enddo > enddo > endif > 1839a1755,1756 > ! write(6,*)' at 1 in raf4, ngauss,npass,i,j,kps,pe=',ngauss,filter(1)%npass,& > ! ips,ipe,jps,jpe,kps,kpe 1840a1758,1759 > ! write(6,*)' k,min,max g=',k,minval(g(1:ngauss,ips:ipe,jps:jpe,k)), & > ! maxval(g(1:ngauss,ips:ipe,jps:jpe,k)) 1850,1858d1768 < ! apply smoothing if nsmooth or nsmooth_shapiro > 0 < if(kps.le.kpe) then < do i=1,max(filter(1)%nsmooth,filter(1)%nsmooth_shapiro) < call smther(filter(1),g,filter(1)%nrows,ngauss,filter(1)%nsmooth,filter(1)%nsmooth_shapiro, & < ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe, & < mype,npes,filter(1)%gnorm_halox,filter(1)%gnorm_haloy,.false.) < end do < end if < 1864,1871c1774,1776 < if(filter(icolor)%npointsmaxall.gt.0) then < if(filter(1)%new_factorization) then < iadvance=1 < iback=2 < call one_color4_new_factorization(g,filter(icolor),ngauss,ipass,filter(1)%ifilt_ord, & < filter(icolor)%nstrings,filter(icolor)%istart, & < ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < else --- > ! if(mype.eq.0) write(6,*)' at 2 in raf4, icolor,npointsmaxall=', & > ! icolor,filter(icolor)%npointsmaxall > if(filter(icolor)%npointsmaxall.gt.0) & 1873,1875c1778,1782 < filter(icolor)%nstrings,filter(icolor)%istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) < end if < end if --- > filter(icolor)%nstrings,filter(icolor)%istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1884c1791,1795 < SUBROUTINE raf_sm4(g,filter,ngauss,ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > SUBROUTINE raf_sm4(g,filter,ngauss, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1890c1801,1803 < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 1895c1808 < real(r_single), DIMENSION( ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & --- > real(r_single), DIMENSION( ngauss,ims:ime, jms:jme, kms:kme ), INTENT(INOUT) :: & 1900c1813 < integer(i_long) icolor,ipass,iadvance,iback --- > integer(i_long) icolor,ipass 1907,1914c1820 < if(filter(icolor)%npointsmaxall.gt.0) then < if(filter(1)%new_factorization) then < iadvance=1 < iback=2 < call one_color4_new_factorization(g,filter(icolor),ngauss,ipass,filter(1)%ifilt_ord, & < filter(icolor)%nstrings,filter(icolor)%istart, & < ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < else --- > if(filter(icolor)%npointsmaxall.gt.0) & 1916,1918c1822,1826 < filter(icolor)%nstrings,filter(icolor)%istart,ips,ipe,jps,jpe,kps,kpe,mype,npes) < end if < end if --- > filter(icolor)%nstrings,filter(icolor)%istart, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1927c1835,1840 < subroutine sort_strings4(info_string,aspect_full,npoints_recv,ib,mype,npes,nvars) --- > subroutine sort_strings4(info_string,aspect_full, & > npoints_recv,ib, & > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 1933,1934c1846,1850 < INTEGER(i_long), INTENT(IN) :: mype, npes < integer(i_long),intent(in):: nvars --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices > INTEGER(i_long), INTENT(IN) :: & > mype, npes 1946c1862 < integer(i_long),intent(out):: ib(npoints_recv) --- > integer(i_long) ib(npoints_recv) 1948,1949c1864 < integer(i_long) ib0(npoints_recv),ib1(npoints_recv) < integer(i_llong) ij_origin(npoints_recv) --- > integer(i_long) ij_origin(npoints_recv) 1951,1955c1866 < integer(i_short) info2(8,npoints_recv) < integer(i_long) istart(npoints_recv+1) < integer(i_long) numvar < integer(i_long) icountvar(nvars),istartvar(nvars) < real(r_single) aspect2(npoints_recv) --- > real(r_single) work(npoints_recv) 1957,1963c1868,1873 < integer(i_long) i,j < integer(i_llong) idist,idistlen,idistmax,idistmin,idjxlen,idjxylen,idjxyzlen,idjxyzx0len,idjxyzxy0len < integer(i_llong) idjxyzxyz0len < integer(i_long) ii,k < integer(i_llong) ix0,ix0len,ix0max,ix0min,iy0,iy0len,iy0max,iy0min,iz0,iz0len,iz0max < integer(i_llong) iz0min,jumpx,jumpxlen,jumpxmax,jumpxmin,jumpy,jumpylen,jumpymax,jumpymin < integer(i_llong) jumpz,jumpzlen,jumpzmax,jumpzmin --- > integer(i_long) i,idist,idistlen,idistmax,idistmin,idjxlen,idjxylen,idjxyzlen,idjxyzx0len,idjxyzxy0len > integer(i_long) idjxyzxyz0len > integer(i_long) ivar,ivarmax,ivarmin > integer(i_long) ix0,ix0len,ix0max,ix0min,iy0,iy0len,iy0max,iy0min,iz0,iz0len,iz0max > integer(i_long) iz0min,j,jumpx,jumpxlen,jumpxmax,jumpxmin,jumpy,jumpylen,jumpymax,jumpymin > integer(i_long) jumpz,jumpzlen,jumpzmax,jumpzmin 1965,1991d1874 < ! sort by var first (but not with indexxi8, because destroys order of points within a subgroup) < < ii=0 < icountvar=0 < istartvar=0 < ib0=0 < do k=1,nvars < do i=1,npoints_recv < if(info_string(8,i).eq.k) then < icountvar(k)=icountvar(k)+1 < ii=ii+1 < if(icountvar(k).eq.1) istartvar(k)=ii < ib0(ii)=i < do j=1,8 < info2(j,ii)=info_string(j,i) < end do < aspect2(ii)=aspect_full(i) < end if < end do < end do < if(minval(ib0).eq.0) then < write(6,*)' error in sort_strings4, problem with ib0' < call stop2(68) < end if < < ! now sort each section in usual way < 1994,1996d1876 < ij_origin=0 < do k=1,nvars < if(icountvar(k).eq.0) cycle 2004,2007c1884,1888 < do i=istartvar(k),istartvar(k)+icountvar(k)-1 < jumpx=info2(5,i) ; jumpy=info2(6,i) ; jumpz=info2(7,i) < ix0=info2(2,i) ; iy0=info2(3,i) ; iz0=info2(4,i) < idist=info2(1,i) --- > ivarmin=jumpxmin ; ivarmax=jumpxmax > do i=1,npoints_recv > jumpx=info_string(5,i) ; jumpy=info_string(6,i) ; jumpz=info_string(7,i) > ix0=info_string(2,i) ; iy0=info_string(3,i) ; iz0=info_string(4,i) > idist=info_string(1,i) ; ivar=info_string(8,i) 2010a1892 > ivarmin=min(ivar,ivarmin) 2023a1906 > idjxyzxyz0len=idjxyzxy0len*iz0len 2025,2028c1908,1911 < do i=istartvar(k),istartvar(k)+icountvar(k)-1 < jumpx=info2(5,i) ; jumpy=info2(6,i) ; jumpz=info2(7,i) < ix0=info2(2,i) ; iy0=info2(3,i) ; iz0=info2(4,i) < idist=info2(1,i) --- > do i=1,npoints_recv > jumpx=info_string(5,i) ; jumpy=info_string(6,i) ; jumpz=info_string(7,i) > ix0=info_string(2,i) ; iy0=info_string(3,i) ; iz0=info_string(4,i) > idist=info_string(1,i) ; ivar=info_string(8,i) 2035c1918,1919 < +idjxyzxy0len*(iz0-iz0min) --- > +idjxyzxy0len*(iz0-iz0min) & > +idjxyzxyz0len*(ivar-ivarmin) 2037,2039c1921,1929 < call indexxi8(icountvar(k),ij_origin(istartvar(k)),ib1(istartvar(k))) < do i=istartvar(k),istartvar(k)+icountvar(k)-1 < ib1(i)=ib1(i)+istartvar(k)-1 --- > call indexxi4(npoints_recv,ij_origin,ib) > > do j=1,8 > do i=1,npoints_recv > iwork(i)=info_string(j,ib(i)) > end do > do i=1,npoints_recv > info_string(j,i)=iwork(i) > end do 2041d1930 < end do 2043c1932 < ib(i)=ib0(ib1(i)) --- > work(i)=aspect_full(ib(i)) 2046,2049c1935 < do j=1,8 < info_string(j,i)=info2(j,ib1(i)) < end do < aspect_full(i)=aspect2(ib1(i)) --- > aspect_full(i)=work(i) 2057c1943,1946 < ips,ipe,jps,jpe,kps,kpe,mype,npes) --- > ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme, & ! memory indices > mype, npes) 2062,2063c1951,1955 < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices < INTEGER(i_long), INTENT(IN) :: mype, npes --- > INTEGER(i_long), INTENT(IN) :: ids, ide, jds, jde, kds, kde, & ! domain indices > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices > INTEGER(i_long), INTENT(IN) :: & > mype, npes 2093a1986,1988 > integer(i_long) idest(npoints_send) > integer(i_long) indx(npoints_send) > integer(i_short) iwork(npoints_send) 2098d1992 < integer(i_long) idestlast 2106d1999 < idestlast=-1 2114,2118d2006 < if(idestpe.lt.idestlast) then < write(6,*)'STRING_ASSEMBLE4: ***PROBLEM*** destination pes out of order for label_string' < call stop2(68) < end if < idestlast=idestpe 2125a2014 > idest(mbuf)=idestpe 2138a2028,2069 > ! sort destination pe numbers from smallest to largest > > if(mbuf.gt.0) then > call indexxi4(mbuf,idest,indx) > > ! use sort index to reorder everything > > do j=1,8 > do i=1,mbuf > iwork(i)=string_info(j,indx(i)) > end do > do i=1,mbuf > string_info(j,i)=iwork(i) > end do > end do > do i=1,mbuf > iwork(i)=ia(indx(i)) > end do > do i=1,mbuf > ia(i)=iwork(i) > end do > do i=1,mbuf > iwork(i)=ja(indx(i)) > end do > do i=1,mbuf > ja(i)=iwork(i) > end do > do i=1,mbuf > iwork(i)=ka(indx(i)) > end do > do i=1,mbuf > ka(i)=iwork(i) > end do > > do i=1,mbuf > work(i)=full_aspect(indx(i)) > end do > do i=1,mbuf > full_aspect(i)=work(i) > end do > end if > 2179c2110 < mype, npes,icolor ) --- > ims, ime, jms, jme, kms, kme, mype, npes,icolor ) ! memory indices 2187c2118,2119 < ips, ipe, jps, jpe, kps, kpe ! patch indices --- > ips, ipe, jps, jpe, kps, kpe, & ! patch indices > ims, ime, jms, jme, kms, kme ! memory indices 2191c2123 < INTEGER(i_short), DIMENSION( 3, * ), INTENT(INout) :: & --- > INTEGER(i_short), DIMENSION( 3, * ), INTENT(IN) :: & 2194c2126 < INTEGER(i_short), DIMENSION( 5, * ), INTENT(INout) :: & --- > INTEGER(i_short), DIMENSION( 5, * ), INTENT(IN) :: & 2210,2212d2141 < INTEGER(i_short), DIMENSION( 6, (ipe-ips+1)*(jpe-jps+1)*(kpe-kps+1) ) :: label_string2 < INTEGER(i_short), DIMENSION( 3, (ipe-ips+1)*(jpe-jps+1)*(kpe-kps+1) ) :: i1filter2 < INTEGER(i_short), DIMENSION( 5, (ipe-ips+1)*(jpe-jps+1)*(kpe-kps+1) ) :: i2filter2 2218a2148 > ! if(mype.eq.0.and.icolor.eq.3) write(6,*)' at 1 in string_label, nstrings=',nstrings 2231,2232c2161,2162 < label_string2(1,n)=i ; label_string2(2,n)=j ; label_string2(3,n)=k < label_string2(4,n)=idist ; label_string2(6,n)=ivar --- > label_string(1,n)=i ; label_string(2,n)=j ; label_string(3,n)=k > label_string(4,n)=idist ; label_string(6,n)=ivar 2237,2238c2167,2168 < label_string2(1,n)=i ; label_string2(2,n)=j ; label_string2(3,n)=k < label_string2(4,n)=idist ; label_string2(6,n)=ivar --- > label_string(1,n)=i ; label_string(2,n)=j ; label_string(3,n)=k > label_string(4,n)=idist ; label_string(6,n)=ivar 2243,2244c2173,2174 < label_string2(1,n)=i ; label_string2(2,n)=j ; label_string2(3,n)=k < label_string2(4,n)=idist ; label_string2(6,n)=ivar --- > label_string(1,n)=i ; label_string(2,n)=j ; label_string(3,n)=k > label_string(4,n)=idist ; label_string(6,n)=ivar 2250,2251c2180,2181 < ivar=label_string2(6,n) < labelijk(n)=label_string2(1,n)+(ide-ids+1)*(label_string2(2,n)-1+(jde-jds+1)*(label_string2(3,n)-1 & --- > ivar=label_string(6,n) > labelijk(n)=label_string(1,n)+(ide-ids+1)*(label_string(2,n)-1+(jde-jds+1)*(label_string(3,n)-1 & 2319c2249 < label_string2(5,i)=mpe --- > label_string(5,i)=mpe 2321,2344d2250 < < !---- reorder label_string in order of increasing pe number < allocate(index(nstrings)) < call indexxi8(nstrings,labelijk,index) < do i=1,nstrings < do j=1,6 < label_string(j,i)=label_string2(j,index(i)) < end do < do j=1,3 < i1filter2(j,i)=i1filter(j,index(i)) < end do < do j=1,5 < i2filter2(j,i)=i2filter(j,index(i)) < end do < end do < do i=1,nstrings < do j=1,3 < i1filter(j,i)=i1filter2(j,i) < end do < do j=1,5 < i2filter(j,i)=i2filter2(j,i) < end do < end do < deallocate(index) 2454,3217d2359 < subroutine add_halox0(filter,nrows,ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes) < < ! setup everything needed for adding halo rows when required < < implicit none < < TYPE(filter_cons), INTENT(inOUT) :: filter < integer(i_long), intent(in):: nrows < integer(i_long), intent(in):: ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes < < integer(i_long) ijglob_pe(ids:ide,jds:jde),ijglob_pe0(ids:ide,jds:jde) < integer(i_long) nrecv_halo(0:npes-1),ndrecv_halo(0:npes) < integer(i_long) nsend_halo(0:npes-1),ndsend_halo(0:npes) < integer(i_long) i,i0,ii,ir,istat,j,mpe < integer(i_long) info_recv_halo(2,npes*2*nrows*(ide-ids+1)) < integer(i_long) info_send_halo(2,npes*2*nrows*(ide-ids+1)) < integer(i_long) iorigin(npes*2*nrows*(ide-ids+1)) < integer(i_long) indx(npes*2*nrows*(ide-ids+1)) < integer(i_long) iwork(npes*2*nrows*(ide-ids+1)) < integer(i_long) nrecv_halo_loc,nsend_halo_loc,mpi_string1,ierror < < if(npes.eq.1.or.nrows.eq.0) return < if(ips.eq.ids.and.ipe.eq.ide) return < < ijglob_pe0=0 < do j=jps,jpe < do i=ips,ipe < ijglob_pe0(i,j)=mype < end do < end do < call mpi_allreduce(ijglob_pe0,ijglob_pe,(ide-ids+1)*(jde-jds+1),mpi_integer4,mpi_sum,mpi_comm_world,ierror) < < ! create list of all points to be recieved < ii=0 < nrecv_halo=0 < do i0=ips-nrows,ipe+1,ipe+1-ips+nrows < do ir=0,nrows-1 < i=i0+ir < if(i.lt.ids.or.i.gt.ide) cycle < do j=jps,jpe < ii=ii+1 < info_recv_halo(1,ii)=i ; info_recv_halo(2,ii)=j < iorigin(ii)=ijglob_pe(i,j) < nrecv_halo(ijglob_pe(i,j))=nrecv_halo(ijglob_pe(i,j))+1 < end do < end do < end do < < ndrecv_halo(0)=0 < do mpe=1,npes < ndrecv_halo(mpe)=ndrecv_halo(mpe-1)+nrecv_halo(mpe-1) < end do < < call mpi_alltoall(nrecv_halo,1,mpi_integer4,nsend_halo,1,mpi_integer4,mpi_comm_world,ierror) < ndsend_halo(0)=0 < do mpe=1,npes < ndsend_halo(mpe)=ndsend_halo(mpe-1)+nsend_halo(mpe-1) < end do < nsend_halo_loc=ndsend_halo(npes) < nrecv_halo_loc=ndrecv_halo(npes) < < ! sort origin pe numbers from smallest to largest < if(ii.gt.0) then < call indexxi4(ii,iorigin,indx) < < ! use sort index to reorder < do j=1,2 < do i=1,ii < iwork(i)=info_recv_halo(j,indx(i)) < end do < do i=1,ii < info_recv_halo(j,i)=iwork(i) < end do < end do < end if < < call mpi_type_contiguous(2,mpi_integer4,mpi_string1,ierror) < call mpi_type_commit(mpi_string1,ierror) < call mpi_alltoallv(info_recv_halo,nrecv_halo,ndrecv_halo,mpi_string1, & < info_send_halo,nsend_halo,ndsend_halo,mpi_string1,mpi_comm_world,ierror) < call mpi_type_free(mpi_string1,ierror) < < filter%nsend_halox_loc=nsend_halo_loc < filter%nrecv_halox_loc=nrecv_halo_loc < deallocate(filter%nrecv_halox,stat=istat) < allocate(filter%nrecv_halox(0:npes-1)) < do i=0,npes-1 < filter%nrecv_halox(i)=nrecv_halo(i) < end do < deallocate(filter%ndrecv_halox,stat=istat) < allocate(filter%ndrecv_halox(0:npes)) < do i=0,npes < filter%ndrecv_halox(i)=ndrecv_halo(i) < end do < deallocate(filter%nsend_halox,stat=istat) < allocate(filter%nsend_halox(0:npes-1)) < do i=0,npes-1 < filter%nsend_halox(i)=nsend_halo(i) < end do < deallocate(filter%ndsend_halox,stat=istat) < allocate(filter%ndsend_halox(0:npes)) < do i=0,npes < filter%ndsend_halox(i)=ndsend_halo(i) < end do < deallocate(filter%info_send_halox,stat=istat) < allocate(filter%info_send_halox(2,nsend_halo_loc)) < do i=1,nsend_halo_loc < do j=1,2 < filter%info_send_halox(j,i)=info_send_halo(j,i) < end do < end do < deallocate(filter%info_recv_halox,stat=istat) < allocate(filter%info_recv_halox(2,nrecv_halo_loc)) < do i=1,nrecv_halo_loc < do j=1,2 < filter%info_recv_halox(j,i)=info_recv_halo(j,i) < end do < end do < < end subroutine add_halox0 < < subroutine add_haloy0(filter,nrows,ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes) < < ! setup everything needed for adding halo rows when required < < implicit none < < TYPE(filter_cons), INTENT(inOUT) :: filter < integer(i_long), intent(in):: nrows < integer(i_long), intent(in):: ids,ide,jds,jde,ips,ipe,jps,jpe,mype,npes < < integer(i_long) ijglob_pe(ids:ide,jds:jde),ijglob_pe0(ids:ide,jds:jde) < integer(i_long) nrecv_halo(0:npes-1),ndrecv_halo(0:npes) < integer(i_long) nsend_halo(0:npes-1),ndsend_halo(0:npes) < integer(i_long) i,ii,istat,j,j0,jr,mpe < integer(i_long) info_recv_halo(2,npes*2*nrows*(jde-jds+1)) < integer(i_long) info_send_halo(2,npes*2*nrows*(jde-jds+1)) < integer(i_long) iorigin(npes*2*nrows*(jde-jds+1)) < integer(i_long) indx(npes*2*nrows*(jde-jds+1)) < integer(i_long) iwork(npes*2*nrows*(jde-jds+1)) < integer(i_long) nrecv_halo_loc,nsend_halo_loc,mpi_string1,ierror < < if(npes.eq.1.or.nrows.eq.0) return < if(jps.eq.jds.and.jpe.eq.jde) return < < ijglob_pe0=0 < do j=jps,jpe < do i=ips,ipe < ijglob_pe0(i,j)=mype < end do < end do < call mpi_allreduce(ijglob_pe0,ijglob_pe,(ide-ids+1)*(jde-jds+1),mpi_integer4,mpi_sum,mpi_comm_world,ierror) < < ! create list of all points to be recieved < ii=0 < nrecv_halo=0 < do j0=jps-nrows,jpe+1,jpe+1-jps+nrows < do jr=0,nrows-1 < j=j0+jr < if(j.lt.jds.or.j.gt.jde) cycle < do i=ips,ipe < ii=ii+1 < info_recv_halo(1,ii)=i ; info_recv_halo(2,ii)=j < iorigin(ii)=ijglob_pe(i,j) < nrecv_halo(ijglob_pe(i,j))=nrecv_halo(ijglob_pe(i,j))+1 < end do < end do < end do < < ndrecv_halo(0)=0 < do mpe=1,npes < ndrecv_halo(mpe)=ndrecv_halo(mpe-1)+nrecv_halo(mpe-1) < end do < < call mpi_alltoall(nrecv_halo,1,mpi_integer4,nsend_halo,1,mpi_integer4,mpi_comm_world,ierror) < ndsend_halo(0)=0 < do mpe=1,npes < ndsend_halo(mpe)=ndsend_halo(mpe-1)+nsend_halo(mpe-1) < end do < nsend_halo_loc=ndsend_halo(npes) < nrecv_halo_loc=ndrecv_halo(npes) < < ! sort origin pe numbers from smallest to largest < if(ii.gt.0) then < call indexxi4(ii,iorigin,indx) < < ! use sort index to reorder < do j=1,2 < do i=1,ii < iwork(i)=info_recv_halo(j,indx(i)) < end do < do i=1,ii < info_recv_halo(j,i)=iwork(i) < end do < end do < end if < < call mpi_type_contiguous(2,mpi_integer4,mpi_string1,ierror) < call mpi_type_commit(mpi_string1,ierror) < call mpi_alltoallv(info_recv_halo,nrecv_halo,ndrecv_halo,mpi_string1, & < info_send_halo,nsend_halo,ndsend_halo,mpi_string1,mpi_comm_world,ierror) < call mpi_type_free(mpi_string1,ierror) < < filter%nsend_haloy_loc=nsend_halo_loc < filter%nrecv_haloy_loc=nrecv_halo_loc < deallocate(filter%nrecv_haloy,stat=istat) < allocate(filter%nrecv_haloy(0:npes-1)) < do i=0,npes-1 < filter%nrecv_haloy(i)=nrecv_halo(i) < end do < deallocate(filter%ndrecv_haloy,stat=istat) < allocate(filter%ndrecv_haloy(0:npes)) < do i=0,npes < filter%ndrecv_haloy(i)=ndrecv_halo(i) < end do < deallocate(filter%nsend_haloy,stat=istat) < allocate(filter%nsend_haloy(0:npes-1)) < do i=0,npes-1 < filter%nsend_haloy(i)=nsend_halo(i) < end do < deallocate(filter%ndsend_haloy,stat=istat) < allocate(filter%ndsend_haloy(0:npes)) < do i=0,npes < filter%ndsend_haloy(i)=ndsend_halo(i) < end do < deallocate(filter%info_send_haloy,stat=istat) < allocate(filter%info_send_haloy(2,nsend_halo_loc)) < do i=1,nsend_halo_loc < do j=1,2 < filter%info_send_haloy(j,i)=info_send_halo(j,i) < end do < end do < deallocate(filter%info_recv_haloy,stat=istat) < allocate(filter%info_recv_haloy(2,nrecv_halo_loc)) < do i=1,nrecv_halo_loc < do j=1,2 < filter%info_recv_haloy(j,i)=info_recv_halo(j,i) < end do < end do < < end subroutine add_haloy0 < < subroutine add_halo_x(f,g,filter,ngauss,nrows,ids,ide,ips,ipe,jps,jpe,kps,kpe,mype,npes) < < implicit none < < TYPE(filter_cons), INTENT(in) :: filter < integer(i_long), intent(in):: ngauss,nrows < integer(i_long), intent(in):: ids,ide,ips,ipe,jps,jpe,kps,kpe,mype,npes < real(r_single),intent(in):: f(ngauss,ips:ipe,jps:jpe,kps:kpe) < real(r_single),intent(out):: g(ngauss,max(ids,ips-nrows):min(ide,ipe+nrows),jps:jpe,kps:kpe) < < integer(i_long) i,j,k,mpi_string2,n,ierror < real(r_single),allocatable:: bufsend(:,:,:),bufrecv(:,:,:) < < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=f(n,i,j,k) < end do < end do < end do < end do < < if(npes.eq.1.or.nrows.eq.0.or.(ips.eq.ids.and.ipe.eq.ide)) return < < ! now gather up points to send < < allocate(bufsend(ngauss,kps:kpe,filter%nsend_halox_loc),bufrecv(ngauss,kps:kpe,filter%nrecv_halox_loc)) < do i=1,filter%nsend_halox_loc < do k=kps,kpe < do n=1,ngauss < bufsend(n,k,i)=f(n,filter%info_send_halox(1,i),filter%info_send_halox(2,i),k) < end do < end do < end do < call mpi_type_contiguous(ngauss*(kpe-kps+1),mpi_real4,mpi_string2,ierror) < call mpi_type_commit(mpi_string2,ierror) < call mpi_alltoallv(bufsend,filter%nsend_halox,filter%ndsend_halox,mpi_string2, & < bufrecv,filter%nrecv_halox,filter%ndrecv_halox,mpi_string2,mpi_comm_world,ierror) < call mpi_type_free(mpi_string2,ierror) < ! finally distribute points back < do i=1,filter%nrecv_halox_loc < do k=kps,kpe < do n=1,ngauss < g(n,filter%info_recv_halox(1,i),filter%info_recv_halox(2,i),k)=bufrecv(n,k,i) < end do < end do < end do < deallocate(bufsend,bufrecv) < < end subroutine add_halo_x < < subroutine add_halo_y(f,g,filter,ngauss,nrows,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) < < implicit none < < TYPE(filter_cons), INTENT(in) :: filter < integer(i_long), intent(in):: ngauss,nrows < integer(i_long), intent(in):: jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes < real(r_single),intent(in):: f(ngauss,ips:ipe,jps:jpe,kps:kpe) < real(r_single),intent(out):: g(ngauss,ips:ipe,max(jds,jps-nrows):min(jde,jpe+nrows),kps:kpe) < < integer(i_long) i,j,k,mpi_string2,n,ierror < real(r_single),allocatable:: bufsend(:,:,:),bufrecv(:,:,:) < < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=f(n,i,j,k) < end do < end do < end do < end do < < if(npes.eq.1.or.nrows.eq.0.or.(jps.eq.jds.and.jpe.eq.jde)) return < < ! now gather up points to send < < allocate(bufsend(ngauss,kps:kpe,filter%nsend_haloy_loc),bufrecv(ngauss,kps:kpe,filter%nrecv_haloy_loc)) < do i=1,filter%nsend_haloy_loc < do k=kps,kpe < do n=1,ngauss < bufsend(n,k,i)=f(n,filter%info_send_haloy(1,i),filter%info_send_haloy(2,i),k) < end do < end do < end do < call mpi_type_contiguous(ngauss*(kpe-kps+1),mpi_real4,mpi_string2,ierror) < call mpi_type_commit(mpi_string2,ierror) < call mpi_alltoallv(bufsend,filter%nsend_haloy,filter%ndsend_haloy,mpi_string2, & < bufrecv,filter%nrecv_haloy,filter%ndrecv_haloy,mpi_string2,mpi_comm_world,ierror) < call mpi_type_free(mpi_string2,ierror) < ! finally distribute points back < do i=1,filter%nrecv_haloy_loc < do k=kps,kpe < do n=1,ngauss < g(n,filter%info_recv_haloy(1,i),filter%info_recv_haloy(2,i),k)=bufrecv(n,k,i) < end do < end do < end do < deallocate(bufsend,bufrecv) < < end subroutine add_halo_y < < !><><>>><><><><><><<><><<<>< following is new code for implementation of approximate 4th order factorization < !><><><<><><><><><><><>><><>< which yields almost a 4th order result for the cost of a 2nd order filter. < !><><><><><>><><><<> this is because the approximate factorization breaks the 4th order operator < !><><><><>><><><><>>< into 4 factors, which are then recombined in such a way that the last < !><><<><><<><>< < subroutine one_color4_new_factorization(g,filter,ngauss,ipass,ifilt_ord, & < nstrings,istart,ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < < ! apply one forward-backward recursive filter for one color. < ! use new approximate 4th order factorization < < < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices < < INTEGER(i_long), INTENT(IN) :: & < mype, npes,ngauss,iadvance,iback < < INTEGER(i_long), INTENT(IN) :: & < ipass ! total number of contiguous string points < integer(i_long),intent(in):: ifilt_ord < < real(r_single), DIMENSION( ngauss,ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & < g ! input--field on grid, output--filtered field on grid < < integer(i_long),intent(in):: nstrings < integer(i_long),intent(in):: istart(nstrings+1) < type(filter_cons) filter < < real(r_single) work(ngauss,max(1,filter%npointsmax),2) < real(r_single) work2(max(1,filter%npointsmax)) < < integer(i_long) i,ierr,igauss,ishort_end,j,l,mpi_string < < !-- gather up strings < < call mpi_type_contiguous(ngauss,mpi_real4,mpi_string,ierr) < call mpi_type_commit(mpi_string,ierr) < < if(filter%npoints_send.gt.0) then < do i=1,filter%npoints_send < do igauss=1,ngauss < work(igauss,i,1)=g(igauss,filter%ia(i),filter%ja(i),filter%ka(i)) < end do < end do < end if < if(npes.eq.1.and.filter%npoints_recv.gt.0) then < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(igauss,i,2)=work(igauss,i,1) < end do < end do < else < call mpi_alltoallv(work(1,1,1),filter%nsend,filter%ndsend,mpi_string, & < work(1,1,2),filter%nrecv,filter%ndrecv,mpi_string,mpi_comm_world,ierr) < end if < if(filter%npoints_recv.gt.0) then < do igauss=1,ngauss < do i=1,filter%npoints_recv < work2(i)=work(igauss,filter%ib(i),2) < end do < < do j=1,nstrings < do i=istart(j),istart(j+1)-1 < do l=1,min(ifilt_ord,i-istart(j)) < work2(i)=work2(i)-filter%fmat(l,i,ipass,igauss,iadvance)*work2(i-l) < end do < work2(i)=filter%fmat0(i,ipass,igauss,iadvance)*work2(i) < end do < do i=istart(j+1)-1,istart(j),-1 < do l=1,min(ifilt_ord,istart(j+1)-i-1) < work2(i)=work2(i)-filter%fmat(l,i+l,ipass,igauss,iback)*work2(i+l) < end do < work2(i)=filter%fmat0(i,ipass,igauss,iback)*work2(i) < end do < end do < < do i=1,filter%npoints_recv < work(igauss,i,1)=work2(i) < end do < end do < < !-- send strings back < < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(igauss,filter%ib(i),2)=work(igauss,i,1) < end do < end do < end if < if(npes.eq.1.and.filter%npoints_recv.gt.0) then < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(igauss,i,1)=work(igauss,i,2) < end do < end do < else < call mpi_alltoallv(work(1,1,2),filter%nrecv,filter%ndrecv,mpi_string, & < work(1,1,1),filter%nsend,filter%ndsend,mpi_string,mpi_comm_world,ierr) < end if < if(filter%npoints_send.gt.0) then < do i=1,filter%npoints_send < do igauss=1,ngauss < g(igauss,filter%ia(i),filter%ja(i),filter%ka(i))=work(igauss,i,1) < end do < end do < end if < < call mpi_type_free(mpi_string,ierr) < < end subroutine one_color4_new_factorization < < subroutine one_color24_new_factorization(g,filter,ngauss,ipass,ifilt_ord, & < nstrings,istart,ips,ipe,jps,jpe,kps,kpe,mype,npes,iadvance,iback) < < ! apply one forward-backward recursive filter for one color < < < INTEGER(i_long), INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe ! patch indices < < INTEGER(i_long), INTENT(IN) :: & < mype, npes,ngauss,iadvance,iback < < INTEGER(i_long), INTENT(IN) :: & < ipass ! total number of contiguous string points < integer(i_long),intent(in):: ifilt_ord < < real(r_single), DIMENSION(2,ngauss, ips:ipe, jps:jpe, kps:kpe ), INTENT(INOUT) :: & < g ! input--field on grid, output--filtered field on grid < < integer(i_long),intent(in):: nstrings < integer(i_long),intent(in):: istart(nstrings+1) < type(filter_cons) filter < < real(r_single) work(2,ngauss,max(1,filter%npointsmax),2) < real(r_single) work2(max(1,filter%npointsmax)) < < integer(i_long) i,ierr,igauss,ii,ishort_end,j,l,mpi_string < < !-- gather up strings < < call mpi_type_contiguous(2*ngauss,mpi_real4,mpi_string,ierr) < call mpi_type_commit(mpi_string,ierr) < < if(filter%npoints_send.gt.0) then < do i=1,filter%npoints_send < do igauss=1,ngauss < work(1,igauss,i,1)=g(1,igauss,filter%ia(i),filter%ja(i),filter%ka(i)) < work(2,igauss,i,1)=g(2,igauss,filter%ia(i),filter%ja(i),filter%ka(i)) < end do < end do < end if < if(npes.eq.1.and.filter%npoints_recv.gt.0) then < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(1,igauss,i,2)=work(1,igauss,i,1) < work(2,igauss,i,2)=work(2,igauss,i,1) < end do < end do < else < call mpi_alltoallv(work(1,1,1,1),filter%nsend,filter%ndsend,mpi_string, & < work(1,1,1,2),filter%nrecv,filter%ndrecv,mpi_string,mpi_comm_world,ierr) < end if < if(filter%npoints_recv.gt.0) then < do igauss=1,ngauss < do ii=1,2 < do i=1,filter%npoints_recv < work2(i)=work(ii,igauss,filter%ib(i),2) < end do < < do j=1,nstrings < do i=istart(j),istart(j+1)-1 < do l=1,min(ifilt_ord,i-istart(j)) < work2(i)=work2(i)-filter%fmat(l,i,ipass,igauss,iadvance)*work2(i-l) < end do < work2(i)=filter%fmat0(i,ipass,igauss,iadvance)*work2(i) < end do < do i=istart(j+1)-1,istart(j),-1 < do l=1,min(ifilt_ord,istart(j+1)-i-1) < work2(i)=work2(i)-filter%fmat(l,i+l,ipass,igauss,iback)*work2(i+l) < end do < work2(i)=filter%fmat0(i,ipass,igauss,iback)*work2(i) < end do < end do < < do i=1,filter%npoints_recv < work(ii,igauss,i,1)=work2(i) < end do < end do < end do < < !-- send strings back < < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(1,igauss,filter%ib(i),2)=work(1,igauss,i,1) < work(2,igauss,filter%ib(i),2)=work(2,igauss,i,1) < end do < end do < end if < if(npes.eq.1.and.filter%npoints_recv.gt.0) then < do i=1,filter%npoints_recv < do igauss=1,ngauss < work(1,igauss,i,1)=work(1,igauss,i,2) < work(2,igauss,i,1)=work(2,igauss,i,2) < end do < end do < else < call mpi_alltoallv(work(1,1,1,2),filter%nrecv,filter%ndrecv,mpi_string, & < work(1,1,1,1),filter%nsend,filter%ndsend,mpi_string,mpi_comm_world,ierr) < end if < if(filter%npoints_send.gt.0) then < do i=1,filter%npoints_send < do igauss=1,ngauss < g(1,igauss,filter%ia(i),filter%ja(i),filter%ka(i))=work(1,igauss,i,1) < g(2,igauss,filter%ia(i),filter%ja(i),filter%ka(i))=work(2,igauss,i,1) < end do < end do < end if < < call mpi_type_free(mpi_string,ierr) < < end subroutine one_color24_new_factorization < < subroutine new_alpha_beta4(info_string,aspect_full,rgauss,fmat,fmat0,igauss,ngauss, & < istart_out,npoints_mype,binomial,npass, & < lenbar,lenmax,lenmin,npoints1,mype,npes,nvars) < < < ! compute recursion constants alpha and beta along unbroken strings < < IMPLICIT NONE < < INTEGER(i_long), INTENT(IN) :: mype, npes,nvars < < INTEGER(i_long), INTENT(IN) :: npoints_mype,npass < < REAL(r_double), DIMENSION( 10, 10 ), INTENT(IN) :: binomial < < INTEGER(i_short), DIMENSION( 8, npoints_mype ), INTENT(IN) :: & < info_string ! 1---- distance from origin to current point < ! 2,3,4-- origin coordinates < ! 5,6,7,8-- jumpx,jumpy,jumpz,ivar for this string < real(r_single), DIMENSION( npoints_mype ) , INTENT(IN) :: & < aspect_full < integer(i_long) igauss,ngauss < real(r_double) rgauss < real(r_single) fmat(2,npoints_mype,npass,ngauss,2),fmat0(npoints_mype,npass,ngauss,2) < integer(i_long) istart_out(*) < integer(i_long) lenmax(nvars),lenmin(nvars),npoints1(nvars) ! diagnostic output--to look at string < real(r_double) lenbar(nvars) < < integer(i_long) i,iend,ipass,istart,ivar < integer(i_long) nstrings < < nstrings=0 < < istart=1 < if(npoints_mype.gt.1) then < do i=2,npoints_mype < if(info_string(1,i).ne.info_string(1,i-1)+1.or. & < info_string(2,i).ne.info_string(2,i-1).or. & < info_string(3,i).ne.info_string(3,i-1).or. & < info_string(4,i).ne.info_string(4,i-1).or. & < info_string(5,i).ne.info_string(5,i-1).or. & < info_string(6,i).ne.info_string(6,i-1).or. & < info_string(7,i).ne.info_string(7,i-1).or. & < info_string(8,i).ne.info_string(8,i-1)) then < iend=i-1 < if(igauss.eq.1) then < ivar=info_string(8,iend) < lenbar(ivar)=lenbar(ivar)+(iend-istart+1) < lenmax(ivar)=max(iend-istart+1,lenmax(ivar)) < lenmin(ivar)=min(iend-istart+1,lenmin(ivar)) < if(iend.eq.istart) npoints1(ivar)=npoints1(ivar)+1 < end if < nstrings=nstrings+1 < istart_out(nstrings)=istart < istart_out(nstrings+1)=iend+1 < < do ipass=1,npass < call new_alpha_betaa4(aspect_full(istart),rgauss,iend-istart+1,binomial(ipass,npass), & < fmat(1,istart,ipass,igauss,1),fmat0(istart,ipass,igauss,1), & < fmat(1,istart,ipass,igauss,2),fmat0(istart,ipass,igauss,2)) < end do < < istart=iend+1 < end if < end do < end if < iend=npoints_mype < if(igauss.eq.1) then < ivar=info_string(8,iend) < lenbar(ivar)=lenbar(ivar)+(iend-istart+1) < lenmax(ivar)=max(iend-istart+1,lenmax(ivar)) < lenmin(ivar)=min(iend-istart+1,lenmin(ivar)) < if(iend.eq.istart) npoints1(ivar)=npoints1(ivar)+1 < end if < nstrings=nstrings+1 < istart_out(nstrings)=istart < istart_out(nstrings+1)=iend+1 < do ipass=1,npass < < call new_alpha_betaa4(aspect_full(istart),rgauss,iend-istart+1,binomial(ipass,npass), & < fmat(1,istart,ipass,igauss,1),fmat0(istart,ipass,igauss,1), & < fmat(1,istart,ipass,igauss,2),fmat0(istart,ipass,igauss,2)) < end do < < end subroutine new_alpha_beta4 < < subroutine new_alpha_betaa4(aspect,rgauss,ng,binomial,fmat1,fmat01,fmat2,fmat02) < < ! compute constants for special Purser 4th order factorization < < ! --> aspect: correlation scale (squared, i think), grid units < ! --> rgauss: multiplier of aspect, for multiple gaussians--used for fat-tail filter < ! --> ng: length of string < ! --> binomial: weighting factors (perhaps not needed with high-order filter) < ! <-- fmat1,fmat01,fmat2,fmat02 : filter parameters for special factorization < ! < < IMPLICIT NONE < < INTEGER(i_long), INTENT(IN) :: ng < < REAL(r_double), INTENT(IN) :: binomial < < real(r_single), DIMENSION( ng ), INTENT(IN) :: aspect < real(r_double) rgauss < real(r_single) fmat1(2,ng),fmat01(ng) < real(r_single) fmat2(2,ng),fmat02(ng) < < real(r_double) sig(0:ng-1) < real(r_double) fmat(0:ng-1,-2:0,2) < integer(i_long) i,j < < do i=1,ng < sig(i-1)=sqrt(rgauss*aspect(i)*binomial) < end do < call stringop(ng-1,sig,fmat) < < do i=1,ng < fmat1(2,i)=fmat(i-1,-2,1) < fmat1(1,i)=fmat(i-1,-1,1) < fmat01(i)=1._r_double/fmat(i-1,0,1) < fmat2(2,i)=fmat(i-1,-2,2) < fmat2(1,i)=fmat(i-1,-1,2) < fmat02(i)=1._r_double/fmat(i-1,0,2) < end do < < end subroutine new_alpha_betaa4 < < !============================================================================= < subroutine stringop(n,sig,fmat) < !============================================================================= < implicit none < integer(4), intent(IN ):: n < real(8),dimension(0:n), intent(IN ):: sig < real(8),dimension(0:n,-2:0,2),intent(OUT):: fmat < !---------------------------------------------------------------------------- < complex(8),dimension(2) :: rts < real(8), dimension(0:n,-1:1) :: dmat < real(8) :: dmatt < integer(4) :: i,im < data rts/ & < (-3.4588884621354108_8, 1.7779487522437312_8), & < (-0.54111153786458903_8, 5.0095518087248694_8)/ < !============================================================================= < dmat=0 < do i=1,n < im=i-1 < dmatt=sig(im)*sig(i) < dmat(im,1)=-dmatt < dmat(i,-1)=-dmatt < dmat(im,0)= dmat(im,0)+dmatt < dmat(i,0 )= dmat(i,0 )+dmatt < enddo < do i=1,2 < call quadfil(n,dmat,rts(i),fmat(:,:,i)) < enddo < end subroutine stringop < < !============================================================================= < subroutine quadfil(n,dmat,rts,fmat) < !============================================================================= < use module_pmat2,only: mulbb,l1lb < implicit none < integer(4), intent(IN ):: n < real(8),dimension(0:n,-1:1),intent(IN ):: dmat < complex(8), intent(IN ):: rts < real(8),dimension(0:n,-2:0),intent(OUT):: fmat < !----------------------------------------------------------------------------- < real(8),dimension(0:n,-1:1) :: tmat < real(8),dimension(0:n,-2:2) :: umat < real(8) :: a1,a2,rrtsi,qrtsi < complex(8) :: rtsi < integer :: np < !============================================================================= < np=n+1 < rtsi=1/rts < rrtsi=real(rtsi) < qrtsi=imag(rtsi) < a1=-2*rrtsi < a2=rrtsi**2+qrtsi**2 < tmat(0:n,-1:1)=a2*dmat < tmat(0:n,0)=tmat(0:n,0)+a1 < call mulbb(dmat,tmat(0:n,-1:1),umat(0:n,-2:2), np,np, 1,1, 1,1, 2,2) < umat(0:n,0)=umat(0:n,0)+1 < if(n.eq.0) then < fmat=0 < fmat(0,0)=umat(0,0) < else < call l1lb(umat,fmat,np,2) < end if < < end subroutine quadfil < !><><<><><<><>< < !><><>><><><< end of new code < 3554c2696 < DO it=1,100 ! 4000 ! this should be ample --- > DO it=1,4000 ! this should be ample 3673,3805c2815 < subroutine smther(filter,g,nrows,ngauss,nsmooth,nsmooth_shapiro, & < ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes,gnormx,gnormy,adjoint) < < use kinds,only: r_single,i_long < use raflib,only: add_halo_x,add_halo_y,filter_cons < implicit none < < TYPE(filter_cons) filter ! structure defining recursive filter < integer(i_long) nrows,ngauss,nsmooth,nsmooth_shapiro < logical adjoint < integer(i_long) ids,ide,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes < real(r_single),intent(inout):: g(ngauss,ips:ipe,jps:jpe,kps:kpe) < real(r_single) gnormx(ips:ipe),gnormy(jps:jpe) < < integer(i_long) i,j,k,n < real(r_single),allocatable:: gt(:,:,:,:) < < ! smoothing in x direction < allocate(gt(ngauss,max(ids,ips-nrows):min(ide,ipe+nrows),jps:jpe,kps:kpe)) < if(nsmooth_shapiro.gt.0.and..not.adjoint) then < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=gnormx(i)*g(n,i,j,k) < end do < end do < end do < end do < end if < call add_halo_x(g,gt,filter,ngauss,nrows,ids,ide,ips,ipe,jps,jpe,kps,kpe,mype,npes) < < if(nsmooth.gt.0) & < call smther_one_x(gt,g,ngauss,ids,ide,ips,ipe,jps,jpe,kps,kpe) < if(nsmooth_shapiro.gt.0) & < call smther_two_x(gt,g,ngauss,ids,ide,ips,ipe,jps,jpe,kps,kpe) < deallocate(gt) < if(nsmooth_shapiro.gt.0.and.adjoint) then < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=gnormx(i)*g(n,i,j,k) < end do < end do < end do < end do < end if < < ! smoothing in y direction < allocate(gt(ngauss,ips:ipe,max(jds,jps-nrows):min(jde,jpe+nrows),kps:kpe)) < if(nsmooth_shapiro.gt.0.and..not.adjoint) then < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=gnormy(j)*g(n,i,j,k) < end do < end do < end do < end do < end if < call add_halo_y(g,gt,filter,ngauss,nrows,jds,jde,ips,ipe,jps,jpe,kps,kpe,mype,npes) < if(nsmooth.gt.0) & < call smther_one_y(gt,g,ngauss,jds,jde,ips,ipe,jps,jpe,kps,kpe) < if(nsmooth_shapiro.gt.0) & < call smther_two_y(gt,g,ngauss,jds,jde,ips,ipe,jps,jpe,kps,kpe) < deallocate(gt) < if(nsmooth_shapiro.gt.0.and.adjoint) then < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=gnormy(j)*g(n,i,j,k) < end do < end do < end do < end do < end if < < end subroutine smther < < subroutine smther_one_x(gx,g,ngauss,ids,ide,ips,ipe,jps,jpe,kps,kpe) < < use kinds,only: r_single,i_long < implicit none < < integer(i_long) ngauss < integer(i_long) ids,ide,ips,ipe,jps,jpe,kps,kpe < real(r_single),intent(in):: gx(ngauss,max(ids,ips-1):min(ide,ipe+1),jps:jpe,kps:kpe) < real(r_single),intent(out):: g(ngauss,ips:ipe,jps:jpe,kps:kpe) < < integer(i_long) i,im,ip,j,k,n < < do k=kps,kpe < do j=jps,jpe < do i=ips,ipe < ip=min(i+1,ide) ; im=max(ids,i-1) < do n=1,ngauss < g(n,i,j,k)=.25*(gx(n,ip,j,k)+gx(n,im,j,k)) +.5*gx(n,i,j,k) < end do < end do < end do < end do < < end subroutine smther_one_x < < subroutine smther_one_y(gy,g,ngauss,jds,jde,ips,ipe,jps,jpe,kps,kpe) < < use kinds,only: r_single,i_long < implicit none < < integer(i_long) ngauss < integer(i_long) jds,jde,ips,ipe,jps,jpe,kps,kpe < real(r_single),intent(in):: gy(ngauss,ips:ipe,max(jds,jps-1):min(jde,jpe+1),kps:kpe) < real(r_single),intent(out):: g(ngauss,ips:ipe,jps:jpe,kps:kpe) < < integer(i_long) i,j,jm,jp,k,n < < do k=kps,kpe < do j=jps,jpe < jp=min(j+1,jde) ; jm=max(jds,j-1) < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.25*(gy(n,i,jp,k)+gy(n,i,jm,k)) +.5*gy(n,i,j,k) < end do < end do < end do < end do < < end subroutine smther_one_y < < subroutine smther_two_x(gx,g,ngauss,ids,ide,ips,ipe,jps,jpe,kps,kpe) --- > subroutine smther_one(g1,is,ie,js,je,ns) 3808c2818 < ! subprogram: smther_two --- > ! subprogram: smther_one 3811,3812c2821,2822 < ! abstract: apply second moment preserving "shapiro smoother" < ! in each direction of data slab --- > ! abstract: apply 1-2-1 smoother in each direction of data slab > ! note: this subroutine is self-adjoint 3834,3835c2844,2846 < integer(i_long) ngauss < integer(i_long) ids,ide,ips,ipe,jps,jpe,kps,kpe --- > integer(i_long) is, ie, js, je > integer(i_long) i,j,l,ip,im,jp,jm > integer(i_long), intent(in) :: ns 3837,3838c2848 < real(r_single) gx(ngauss,max(ids,ips-3):min(ide,ipe+3),jps:jpe,kps:kpe) < real(r_single) g(ngauss,ips:ipe,jps:jpe,kps:kpe) --- > real(r_single), dimension(is:ie, js:je), intent(inout) :: g1 3840a2851 > 3842c2853 < integer(i_long) i,im,im3,ip,ip3,istart,iend,j,k,n --- > real(r_single), allocatable:: g2(:,:) 3844,3851c2855,2861 < istart=ips ; iend=ipe < if(ips.eq.ids) then < i=ids ; ip=i+1 ; ip3=i+3 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*gx(n,ip,j,k)-.03125*gx(n,ip3,j,k) < end do --- > 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) 3853,3859c2863,2868 < end do < i=ids+1 ; ip=i+1 ; ip3=i+3 ; im=i-1 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*(gx(n,ip,j,k)+gx(n,im,j,k))-.03125*gx(n,ip3,j,k) < 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) 3861,3908c2870 < end do < i=ids+2 ; ip=i+1 ; ip3=i+3 ; im=i-1 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*(gx(n,ip,j,k)+gx(n,im,j,k))-.03125*gx(n,ip3,j,k) < end do < end do < end do < istart=ids+3 < end if < if(ipe.eq.ide) then < i=ide-2 ; ip=i+1 ; im3=i-3 ; im=i-1 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*(gx(n,ip,j,k)+gx(n,im,j,k))-.03125*gx(n,im3,j,k) < end do < end do < end do < i=ide-1 ; ip=i+1 ; im3=i-3 ; im=i-1 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*(gx(n,ip,j,k)+gx(n,im,j,k))-.03125*gx(n,im3,j,k) < end do < end do < end do < i=ide ; im3=i-3 ; im=i-1 < do k=kps,kpe < do j=jps,jpe < do n=1,ngauss < g(n,i,j,k)=.5*gx(n,i,j,k)+.28125*gx(n,im,j,k)-.03125*gx(n,im3,j,k) < end do < end do < end do < iend=ide-3 < end if < do k=kps,kpe < do j=jps,jpe < do i=istart,iend < ip=i+1 ; im=i-1 ; ip3=i+3 ; im3=i-3 < do n=1,ngauss < g(n,i,j,k)=.28125*(gx(n,ip,j,k)+gx(n,im,j,k))+.5*gx(n,i,j,k)-.03125*(gx(n,ip3,j,k)+gx(n,im3,j,k)) < end do < end do < end do < end do --- > end do 3910c2872,2873 < end subroutine smther_two_x --- > end do > deallocate(g2) 3912c2875,2878 < subroutine smther_two_y(gy,g,ngauss,jds,jde,ips,ipe,jps,jpe,kps,kpe) --- > return > end subroutine smther_one > > subroutine smther_two(g1,is,ie,js,je,ns) 3941,3942c2907,2909 < integer(i_long) ngauss < integer(i_long) jds,jde,ips,ipe,jps,jpe,kps,kpe --- > integer(i_long) is, ie, js, je > integer(i_long) i,j,l,ip,im,jp,jm,ip3,im3,jp3,jm3 > integer(i_long), intent(in) :: ns 3944,3945c2911 < real(r_single) gy(ngauss,ips:ipe,max(jds,jps-3):min(jde,jpe+3),kps:kpe) < real(r_single) g(ngauss,ips:ipe,jps:jpe,kps:kpe) --- > real(r_single), dimension(is:ie, js:je), intent(inout) :: g1 3949c2915 < integer(i_long) i,jm,jm3,jp,jp3,jstart,jend,j,k,n --- > real(r_single), allocatable:: g2(:,:) 3951,3958c2917,2924 < jstart=jps ; jend=jpe < if(jps.eq.jds) then < j=jds ; jp=j+1 ; jp3=j+3 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*gy(n,i,jp,k)-.03125*gy(n,i,jp3,k) < end do --- > 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) > ip3=min(i+3,ie) ; im3=max(is,i-3) > g2(i,j)=.28125*(g1(ip,j)+g1(im,j))+.5*g1(i,j)-.03125*(g1(ip3,j)+g1(im3,j)) 3960,3966c2926,2932 < end do < j=jds+1 ; jp=j+1 ; jp3=j+3 ; jm=j-1 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*(gy(n,i,jp,k)+gy(n,i,jm,k))-.03125*gy(n,i,jp3,k) < end do --- > end do > > do i=is,ie > do j=js,je > jp=min(j+1,je) ; jm=max(js,j-1) > jp3=min(j+3,je) ; jm3=max(js,j-3) > g1(i,j)=.28125*(g2(i,jp)+g2(i,jm))+.5*g2(i,j)-.03125*(g2(i,jp3)+g2(i,jm3)) 3968,4015c2934 < end do < j=jds+2 ; jp=j+1 ; jp3=j+3 ; jm=j-1 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*(gy(n,i,jp,k)+gy(n,i,jm,k))-.03125*gy(n,i,jp3,k) < end do < end do < end do < jstart=jds+3 < end if < if(jpe.eq.jde) then < j=jde-2 ; jp=j+1 ; jm3=j-3 ; jm=j-1 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*(gy(n,i,jp,k)+gy(n,i,jm,k))-.03125*gy(n,i,jm3,k) < end do < end do < end do < j=jde-1 ; jp=j+1 ; jm3=j-3 ; jm=j-1 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*(gy(n,i,jp,k)+gy(n,i,jm,k))-.03125*gy(n,i,jm3,k) < end do < end do < end do < j=jde ; jm3=j-3 ; jm=j-1 < do k=kps,kpe < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.5*gy(n,i,j,k)+.28125*gy(n,i,jm,k)-.03125*gy(n,i,jm3,k) < end do < end do < end do < jend=jde-3 < end if < do k=kps,kpe < do j=jstart,jend < jp=j+1 ; jm=j-1 ; jp3=j+3 ; jm3=j-3 < do i=ips,ipe < do n=1,ngauss < g(n,i,j,k)=.28125*(gy(n,i,jp,k)+gy(n,i,jm,k))+.5*gy(n,i,j,k)-.03125*(gy(n,i,jp3,k)+gy(n,i,jm3,k)) < end do < end do < end do < end do --- > end do 4017c2936,2937 < end subroutine smther_two_y --- > end do > deallocate(g2) 4019c2939,2940 < subroutine smther_two_gnorm(gnorm,ids,ide,ips,ipe) --- > return > end subroutine smther_two 4021c2942,2946 < use kinds,only: r_single,r_double,i_long --- > subroutine tsmther_two(g1,is,ie,js,je,ns) > > ! adjoint of tsmther_two > > use kinds,only: r_single,i_long 4024c2949,2951 < integer(i_long) ids,ide,ips,ipe --- > integer(i_long) is, ie, js, je > integer(i_long) i,j,l,ip,im,jp,jm,ip3,im3,jp3,jm3 > integer(i_long), intent(in) :: ns 4026c2953,2955 < real(r_single) gnorm(ips:ipe) --- > real(r_single), dimension(is:ie, js:je), intent(inout) :: g1 > ! on input: original data slab > ! on ouput: filtered data slab 4028,4029c2957 < integer(i_long) i,ip,ip3,istart,iend < real(r_double) const81,const82 --- > real(r_single), allocatable:: g2(:,:) 4031,4054c2959,2960 < const81=4._8/3._8 < const82=32._8/33._8 < istart=ips ; iend=ipe < if(ips.eq.ids) then < i=ids < gnorm(i)=const81 ! 1./(.5+.28125-.03125) < i=ids+1 < gnorm(i)=const82 ! 1./(.5+.28125*2.-.03125) < i=ids+2 < gnorm(i)=const82 ! 1./(.5+.28125*2.-.03125) < istart=ids+3 < end if < if(ipe.eq.ide) then < i=ide-2 < gnorm(i)=const82 ! 1./(.5+.28125*2.-.03125) < i=ide-1 < gnorm(i)=const82 ! 1./(.5+.28125*2.-.03125) < i=ide < gnorm(i)=const81 ! 1./(.5+.28125-.03125) < iend=ide-3 < end if < do i=istart,iend < gnorm(i)=1. ! 1./(.28125*2.+.5-.03125*2.) < end do --- > allocate(g2(is:ie,js:je)) > do l=1,ns 4056c2962,2992 < end subroutine smther_two_gnorm --- > g2(is:ie,js:je)=0._r_single > do i=is,ie > do j=js,je > jp=min(j+1,je) ; jm=max(js,j-1) > jp3=min(j+3,je) ; jm3=max(js,j-3) > g2(i,jm3)=g2(i,jm3)-.03125*g1(i,j) > g2(i,jp3)=g2(i,jp3)-.03125*g1(i,j) > g2(i,j)=g2(i,j)+.5*g1(i,j) > g2(i,jm)=g2(i,jm)+.28125*g1(i,j) > g2(i,jp)=g2(i,jp)+.28125*g1(i,j) > end do > end do > > g1(is:ie,js:je)=0._r_single > do j=js,je > do i=is,ie > ip=min(i+1,ie) ; im=max(is,i-1) > ip3=min(i+3,ie) ; im3=max(is,i-3) > g1(im3,j)=g1(im3,j)-.03125*g2(i,j) > g1(ip3,j)=g1(ip3,j)-.03125*g2(i,j) > g1(i,j)=g1(i,j)+.5*g2(i,j) > g1(im,j)=g1(im,j)+.28125*g2(i,j) > g1(ip,j)=g1(ip,j)+.28125*g2(i,j) > end do > end do > > end do > deallocate(g2) > > return > end subroutine tSmther_two DIFFERENCES FOR raflib.f90_double_precision DIFFERENCES FOR raflib.f90_no_changes DIFFERENCES FOR raflib.f90_with_shapiro DIFFERENCES FOR rdgrbsst.f90 DIFFERENCES FOR rdgstat_reg.f90 DIFFERENCES FOR read_airs.f90 DIFFERENCES FOR read_amsre.f90 DIFFERENCES FOR read_avhrr.f90 DIFFERENCES FOR read_avhrr_navy.f90 DIFFERENCES FOR read_bufrtovs.f90 DIFFERENCES FOR read_files.f90 DIFFERENCES FOR read_goesimg.f90 DIFFERENCES FOR read_goesndr.f90 DIFFERENCES FOR read_gps.f90 DIFFERENCES FOR read_guess.f90 DIFFERENCES FOR read_l2bufr_mod.f90 DIFFERENCES FOR read_lidar.f90 DIFFERENCES FOR read_modsbufr.f90 DIFFERENCES FOR read_obs.f90 DIFFERENCES FOR read_ozone.F90 DIFFERENCES FOR read_pcp.f90 DIFFERENCES FOR read_prepbufr.f90 1346c1346 < elseif (trim(cgrid) == 'puerto_rico') then --- > elseif (trim(cgrid) == 'prico') then 1394c1394 < elseif (trim(cgrid) == 'puerto_rico') then --- > elseif (trim(cgrid) == 'prico') then 1437c1437 < if (trim(cgrid)=='hawaii' .or. trim(cgrid)=='guam' .or. trim(cgrid)=='puerto_rico') then --- > if (trim(cgrid)=='hawaii' .or. trim(cgrid)=='guam' .or. trim(cgrid)=='prico') then DIFFERENCES FOR read_prepbufr.f90_2Nov2007 DIFFERENCES FOR read_prepbufr.f90_4Jan2008 DIFFERENCES FOR read_prepbufr.f90_no_changes DIFFERENCES FOR read_radar.f90 DIFFERENCES FOR read_ssmi.f90 DIFFERENCES FOR read_ssmis.f90 DIFFERENCES FOR read_superwinds.f90 DIFFERENCES FOR read_wrf_mass_files.f90 DIFFERENCES FOR read_wrf_mass_guess.F90 DIFFERENCES FOR read_wrf_nmm_files.f90 DIFFERENCES FOR read_wrf_nmm_guess.F90 DIFFERENCES FOR regional_io.f90 DIFFERENCES FOR ret_ssmis.f90 DIFFERENCES FOR retrieval_amsre.f90 DIFFERENCES FOR retrieval_mi.f90 DIFFERENCES FOR rfdpar.f90 DIFFERENCES FOR rsearch.F90 DIFFERENCES FOR satthin.F90 DIFFERENCES FOR setupbend.f90 DIFFERENCES FOR setupdw.f90 DIFFERENCES FOR setupoz.f90 DIFFERENCES FOR setuppcp.f90 DIFFERENCES FOR setupps.f90 DIFFERENCES FOR setupps.f90_no_changes DIFFERENCES FOR setuppw.f90 DIFFERENCES FOR setupq.f90 DIFFERENCES FOR setupq.f90_no_changes DIFFERENCES FOR setuprad.f90 DIFFERENCES FOR setupref.f90 DIFFERENCES FOR setuprhsall.f90 DIFFERENCES FOR setuprw.f90 DIFFERENCES FOR setupspd.f90 DIFFERENCES FOR setupsrw.f90 DIFFERENCES FOR setupsst.f90 DIFFERENCES FOR setupt.f90 DIFFERENCES FOR setupt.f90_21Dec2007 DIFFERENCES FOR setupt.f90_3Jan2008 DIFFERENCES FOR setupt.f90_no_changes DIFFERENCES FOR setupw.f90 DIFFERENCES FOR setupw.f90_no_changes DIFFERENCES FOR sfc_model.f90 DIFFERENCES FOR simpin1.f90 DIFFERENCES FOR simpin1_init.f90 DIFFERENCES FOR smooth_polcarf.f90 DIFFERENCES FOR smoothrf.f90 DIFFERENCES FOR smoothwwrf.f90 DIFFERENCES FOR smoothzrf.f90 DIFFERENCES FOR specmod.f90 DIFFERENCES FOR spectral_transforms.f90 DIFFERENCES FOR sst_retrieval.f90 DIFFERENCES FOR statsconv.f90 DIFFERENCES FOR statsoz.f90 DIFFERENCES FOR statspcp.f90 DIFFERENCES FOR statsrad.f90 DIFFERENCES FOR stop1.f90 DIFFERENCES FOR stpcalc.f90 DIFFERENCES FOR stpcld.f90 DIFFERENCES FOR stpdw.f90 DIFFERENCES FOR stpgps.f90 DIFFERENCES FOR stpjc.f90 DIFFERENCES FOR stplimq.f90 DIFFERENCES FOR stpoz.f90 DIFFERENCES FOR stppcp.f90 DIFFERENCES FOR stpps.f90 DIFFERENCES FOR stppw.f90 DIFFERENCES FOR stpq.f90 DIFFERENCES FOR stprad.f90 DIFFERENCES FOR stprw.f90 DIFFERENCES FOR stpspd.f90 DIFFERENCES FOR stpsrw.f90 DIFFERENCES FOR stpsst.f90 DIFFERENCES FOR stpt.f90 DIFFERENCES FOR stpw.f90 DIFFERENCES FOR strong_bal_correction.f90 DIFFERENCES FOR strong_baldiag_inc.f90 DIFFERENCES FOR strong_fast_global_mod.f90 DIFFERENCES FOR strong_slow_global_mod.f90 DIFFERENCES FOR stvp2uv_reg.f90 DIFFERENCES FOR sub2grid.f90 DIFFERENCES FOR support_2dvar.f90 DIFFERENCES FOR support_2dvar.f90_15Nov2007 DIFFERENCES FOR support_2dvar.f90_19Nov2007 DIFFERENCES FOR support_2dvar.f90_20Nov2007 DIFFERENCES FOR support_2dvar.f90_standby DIFFERENCES FOR support_2dvar.f90_wait_20Nov2007 DIFFERENCES FOR tendsmod.f90 DIFFERENCES FOR tintrp2a.f90 DIFFERENCES FOR tintrp3.f90 DIFFERENCES FOR tpause.f90 DIFFERENCES FOR tpause_t.F90 DIFFERENCES FOR transform.f90 DIFFERENCES FOR tstvp2uv_reg.f90 DIFFERENCES FOR turb.f90 DIFFERENCES FOR turbmod.f90 DIFFERENCES FOR turbvars.f90 DIFFERENCES FOR tv_to_tsen.f90 DIFFERENCES FOR unfill_mass_grid2.f90 DIFFERENCES FOR unfill_nmm_grid2.f90 DIFFERENCES FOR unhalf_nmm_grid2.f90 DIFFERENCES FOR update_guess.f90 DIFFERENCES FOR ut_fvAnaGrid.F90 DIFFERENCES FOR wrf_binary_interface.F90 DIFFERENCES FOR wrf_netcdf_interface.F90 DIFFERENCES FOR write_all.f90 DIFFERENCES FOR write_bkgvars_grid.f90 DIFFERENCES FOR wrwrfmassa.F90 DIFFERENCES FOR wrwrfnmma.F90 DIFFERENCES FOR zrnmi_mod.f90