subroutine cnv_to_grib2_mpi(cgrid,skipanlerr) !*********************************************************************** ! abstract: convert rtma background, analysis, and analysis error * ! from binary format to grib2 format * ! * ! program history log: * ! 2005-10-08 pondeca * !*********************************************************************** use mpi use mpitaskmod, only: mype,npe use mpitaskmod, only: is => ista use mpitaskmod, only: ie => iend use mpitaskmod, only: js => jsta use mpitaskmod, only: je => jend use mpitaskmod, only: nx => nlon use mpitaskmod, only: ny => nlat use controlvars, only: ictrwspd10m => iwspd10m !to distinguish from original wspd10m use controlvars, only: ictrtd2m => itd2m !and td2m fields use controlvars, only: igust,ivis,ipblh,idist, & imxtm,imitm, & ipmsl,ihowv,itcamt,ilcbas, & iuwnd10m,ivwnd10m implicit none !Declare passed variables character(60),intent(in):: cgrid logical, intent(in) :: skipanlerr !Option determining whether or not to ! skip computing/outputting the analysis error !Declare local parameters integer(4),parameter:: nvarmax=30 !maximum number of rtma flds that can be processed integer(4),parameter:: n0=9 !# of non-optional bckg, anl and anlerr flds integer(4),parameter:: newvars=14 !# of optional ctr vars. see above "use controlvars" statements real(4),parameter::cldch_bitmapval=1.e+21 !Declare local variables integer(4) nvar2d character(60) varname0(n0) character(60),allocatable::varname(:) character(60),allocatable::typeproc(:) character(30) fname character(60) newvarname(newvars) logical ctrvarpresent(newvars) character(60) filename,filename2,filename3 integer(4) nngust,nnctrwspd10m_0,nnctrwspd10m integer(4) nndist_0,nndist integer(4) nnuwnd10m_0,nnuwnd10m integer(4) nnvwnd10m_0,nnvwnd10m integer(4) nnhowv_0,nnhowv integer(4) n,ncase,i,j,nn,k,i0,j0,maxncase integer(4) itime01 (6) !fcst cycle for guess (ie. reference time) integer(4) itimem1 (6) !fgat integer(4) itimep1 (6) !fgat integer(4) itime0 (6) !analysis integer(4) itime (6) integer(4) ifcsthr,ifcsthrm1,ifcsthrp1 !fcst hours for guess, and fgat fields logical guessfcst,guessfcstm1,guessfcstp1 !.true. for forecast fields logical fexist logical lbiascor logical usefgat logical usectrwspd10m logical use_uwnd10m_vwnd10m integer(4) isign !significance of reference time !!!(0 for analysis and anlerr, and 1 for forecast) integer(4) ityped !type of processed data !!!(0 for analysis and anlerr, and 1 for forecast) integer(4) mx,my integer(4) iopt,iopthgt,npass,npasshgt real(4) smcf1,smcf2 character(60) cgridhalfres logical lhalfres logical lmapgrid logical one2one integer(4) ix,jy character(60) cgridnew integer(4),parameter::nx_cohres=2145 !grid 184 (true ndfd grid) integer(4),parameter::ny_cohres=1377 !grid 184 (true ndfd grid) integer(4),parameter::nx_nwrfc=709 integer(4),parameter::ny_nwrfc=795 integer(4),parameter::i0wexp=200 integer(4) ierror integer(4) is2,ie2,js2,je2 integer(4) is3,ie3,js3,je3 character(60) cgprocess real(4),allocatable,dimension(:,:):: field,field2 real(4),allocatable,dimension(:,:,:):: gndfd,gndfd9 real(4),allocatable,dimension(:,:):: cldchmap namelist/ptimeinfo/itime01,ifcsthr,guessfcst, & usefgat,lbiascor,usectrwspd10m, & use_uwnd10m_vwnd10m namelist/ptimeinfom1/itimem1,ifcsthrm1,guessfcstm1 namelist/ptimeinfop1/itimep1,ifcsthrp1,guessfcstp1 namelist/hgridoptions/lhalfres,lmapgrid,iopt,iopthgt, & smcf1,smcf2,one2one,npass,npasshgt namelist/generatingprocess/cgprocess ! !=================================================================================================== !==> use data statements to set default values for various variables !=================================================================================================== ! data cgprocess/'rtma'/ data lhalfres/.false./ data lmapgrid/.false./ data iopt/1/ data iopthgt/0/ data smcf1/0.50/ data smcf2/0.52/ data one2one/.false./ data npass/1/ data npasshgt/0/ data itime01/2005, 6, 17, 12, 0, 0/ data ifcsthr/1/ data guessfcst/.true./ data usefgat/.false./ data lbiascor/.false./ data usectrwspd10m/.false./ data use_uwnd10m_vwnd10m/.false./ data itimem1/2005, 6, 17, 12, 0, 0/ data ifcsthrm1/0/ data guessfcstm1/.false./ data itimep1/2005, 6, 17, 12, 0, 0/ data ifcsthrp1/2/ data guessfcstp1/.true./ data varname0/ & 'hgt', & 'psfc', & 't2m', & 'td2m', & 'u10m', & 'v10m', & 'q2m', & 'wdir10m', & 'wspd10m'/ ! !=================================================================================================== !==> load arrays newvarname & ctrvarpresent newvarname(1) = 'gust' ; ctrvarpresent(1) = igust>0 ; nngust=1 newvarname(2) = 'vis' ; ctrvarpresent(2) = ivis>0 newvarname(3) = 'pblh' ; ctrvarpresent(3) = ipblh>0 newvarname(4) = 'dist' ; ctrvarpresent(4) = idist>0 ; nndist_0=4 newvarname(5) = 'ctrwspd10m' ; ctrvarpresent(5) = ictrwspd10m>0 ; nnctrwspd10m_0=5 newvarname(6) = 'ctrtd2m' ; ctrvarpresent(6) = ictrtd2m>0 newvarname(7) = 'mxtm' ; ctrvarpresent(7) = imxtm>0 newvarname(8) = 'mitm' ; ctrvarpresent(8) = imitm>0 newvarname(9) = 'pmsl' ; ctrvarpresent(9) = ipmsl>0 newvarname(10) = 'howv' ; ctrvarpresent(10) = ihowv>0 ; nnhowv_0=10 newvarname(11) = 'tcamt' ; ctrvarpresent(11) = itcamt>0 newvarname(12) = 'lcbas' ; ctrvarpresent(12) = ilcbas>0 newvarname(13) = 'uwnd10m' ; ctrvarpresent(13) = iuwnd10m>0 ;nnuwnd10m_0=13 newvarname(14) = 'vwnd10m' ; ctrvarpresent(14) = ivwnd10m>0 ;nnvwnd10m_0=14 newvarname(4) = 'cldch' !rename to cloud ceiling height newvarname(13) = 'u20m' !rename for convenience newvarname(14) = 'v20m' !rename for convenience ! !=================================================================================================== !==> print some useful information and read in namelist ptimeinfo !=================================================================================================== if (mype==0) print*,'in cnv_to_grib2: igust,idist,ivis,ipblh=',igust,idist,ivis,ipblh if (mype==0) print*,'in cnv_to_grib2: ictrwspd10m,ictrtd2m,imxtm,imitm=',ictrwspd10m,ictrtd2m,imxtm,imitm if (mype==0) print*,'in cnv_to_grib2: ihowv,itcamt,ilcbas=',ihowv,itcamt,ilcbas if (mype==0) print*,'in cnv_to_grib2: iuwnd10m,ivwnd10m=',iuwnd10m,ivwnd10m ; if (mype==0) print*,'in cnv_to_grib2: kipanlerr=',skipanlerr ; print* read(9,ptimeinfo) !fetch lbiascor, usefgat, usectrwspd10m and use_uwnd10m_vwnd10m close(9) usefgat=.false. !force this for now /MPondeca, 3Feb2016 lbiascor=.false. !force this for now /MPondeca, 3Feb2016 if (mype==0) then print*,'in cnv_to_grib2: cgrid=',trim(cgrid) print*,'in cnv_to_grib2: nx,ny=',nx,ny print*,'in cnv_to_grib2: lbiascor=',lbiascor print*,'in cnv_to_grib2: usefgat=',usefgat print*,'in cnv_to_grib2: usectrwspd10m=',usectrwspd10m print*,'in cnv_to_grib2: use_uwnd10m_vwnd10m=',use_uwnd10m_vwnd10m endif !=================================================================================================== !==> load varname array !=================================================================================================== ! nndist=-1 nnctrwspd10m=-1 nnuwnd10m=-1 nnvwnd10m=-1 nnhowv=-1 n=0 do nn=1,newvars if (ctrvarpresent(nn)) then n=n+1 if (nn==nndist_0) nndist=n if (nn==nnctrwspd10m_0) nnctrwspd10m=n if (nn==nnuwnd10m_0) nnuwnd10m=n if (nn==nnvwnd10m_0) nnvwnd10m=n if (nn==nnhowv_0) nnhowv=n endif enddo if (.not.ctrvarpresent(nnctrwspd10m_0)) then if (mype==0) then print*,'in cnv_to_grib2: ctrwspd10m not analyzed. set usectrwspd10m to .false.' endif usectrwspd10m=.false. endif if (.not.ctrvarpresent(nnuwnd10m_0) .or. .not.ctrvarpresent(nnvwnd10m_0)) then if (mype==0) then print*,'in cnv_to_grib2: either uwnd10m or vwnd10m or both were not analyzed. set use_uwnd10m_vwnd10m to .false.' endif use_uwnd10m_vwnd10m=.false. endif nvar2d=n0+n if (usectrwspd10m) nvar2d=nvar2d+3 !slots for (psi,chi)-based u10m, v10m, wdir. keep as backup or for diagnostic purposes if (mype==0) print*,'in cnv_to_grib2: nvar2d=',nvar2d if (mype==0) print*,'in cnv_to_grib2: nngust=',nngust if (mype==0) print*,'in cnv_to_grib2: nndist=',nndist if (mype==0) print*,'in cnv_to_grib2: usectrwspd10m=',usectrwspd10m if (mype==0) print*,'in cnv_to_grib2: use_uwnd10m_vwnd10m=',use_uwnd10m_vwnd10m if (mype==0) print*,'in cnv_to_grib2: nnctrwspd10m=',nnctrwspd10m if (mype==0) print*,'in cnv_to_grib2: nnuwnd10m,nnvwnd10m=',nnuwnd10m,nnvwnd10m if (mype==0) print*,'in cnv_to_grib2: nnhowv=',nnhowv allocate(varname(nvar2d)) varname (1:n0) = varname0 (1:n0) n=0 do nn=1,newvars if (ctrvarpresent(nn)) then n=n+1 varname (n0+n) = newvarname(nn) endif enddo if (mype==0) print*,'in cnv_to_grib2: n0+n=',n0+n if (usectrwspd10m) then varname(nvar2d-2)='usfc' varname(nvar2d-1)='vsfc' varname(nvar2d )='wdirsfc' endif ! !=================================================================================================== !==> read in namelist for generating process (rtma or urma) !=================================================================================================== ! inquire(file='generatingprocess_input',exist=fexist) if(fexist) then open (55,file='generatingprocess_input',form='formatted') read (55,generatingprocess) close(55) endif if (mype==0) then print*,'in cnv_to_grib2: cgprocess=',trim(cgprocess) endif ! !=================================================================================================== !==> read in namelist of smoothing options for half-resolution grids !=================================================================================================== ! inquire(file='hgridoptions_input',exist=fexist) if(fexist) then open (55,file='hgridoptions_input',form='formatted') read (55,hgridoptions) close(55) endif if (mype==0) print*,'in cnv_to_grib2: lhalfres=',lhalfres if (mype==0) print*,'in cnv_to_grib2: lmapgrid=',lmapgrid if (lhalfres) then if (mype==0) then print*,'in cnv_to_grib2: iopt=',iopt print*,'in cnv_to_grib2: iopthgt=',iopthgt print*,'in cnv_to_grib2: smcf1,smcf2=',smcf1,smcf2 print*,'in cnv_to_grib2: one2one=',one2one print*,'in cnv_to_grib2: npass=',npass print*,'in cnv_to_grib2: npasshgt=',npasshgt endif endif ! !=================================================================================================== !==> retrieve the analysis time open (55,file='siganl',form='unformatted') read(55) itime0 close(55) if (mype==0) print*,'in cnv_to_grib2: itime0=',itime0 ! !=================================================================================================== !==> populate gndfd with background fields ! (the wind components are relative to the grid). !=================================================================================================== ! allocate(gndfd(is:ie,js:je,nvar2d)) allocate(typeproc(nvar2d)) allocate(field(nx,ny)) allocate(cldchmap(is:ie,js:je)) cldchmap(is:ie,js:je)=+1. if (ctrvarpresent(nndist_0)) then inquire(file='cldch_map_input.dat',exist=fexist) if(fexist) then open (55,file='cldch_map_input.dat',form='unformatted') read (55)field close(55) if (mype==0) print*,'in cnv_to_grib2: cldchmap,min,max=', & minval(field),maxval(field) do j=js,je do i=is,ie cldchmap(i,j)=field(i,j) enddo enddo endif endif maxncase=3 if (skipanlerr) maxncase=2 do 100 ncase=1,maxncase ! ncase= bckg , analysis , anlerr if (mype==0) print*,'in cnv_to_grib2: ncase=',ncase if (ncase==1) fname='sigges' if (ncase==2) fname='siganl' if (ncase==3) fname='anlerr' call retrieve_gesanl(cgrid,nvar2d,n0,gndfd, & fname,newvarname,ctrvarpresent,newvars,usectrwspd10m,use_uwnd10m_vwnd10m, & nngust,nnctrwspd10m,nnuwnd10m,nnvwnd10m,nnhowv_0,nnhowv) ! if (nndist > 0) then ! PRINT*,'FIRST: CLDCH,min,max=',minval(gndfd(:,:,n0+nndist)),maxval(gndfd(:,:,n0+nndist)) ! do j=js,je ! do i=is,ie !! if (cldchmap(i,j)<0. .and. gndfd(i,j,n0+nndist)==0.1) then ! if (gndfd(i,j,n0+nndist)>=20000.) then ! gndfd(i,j,n0+nndist)=cldch_bitmapval ! endif ! enddo ! enddo ! PRINT*,'SECOND: CLDCH,min,max=',minval(gndfd(:,:,n0+nndist)),maxval(gndfd(:,:,n0+nndist)) ! endif !=================================================================================================== !==> prepare to write grib2 file !=================================================================================================== if (ncase.eq.1) then read(9,ptimeinfo) if (mype==0) & print*,'in cnv_to_grib2,itime01,ifcsthr,guessfcst=', & itime01,ifcsthr,guessfcst close(9) if (guessfcst) then itime(1:6)=itime01(1:6) isign=1 ityped=1 typeproc(1:nvar2d)='forecast' else itime(1:6)=itime0(1:6) isign=0 ityped=0 typeproc(1:nvar2d)='analysis' endif if (.not.ctrvarpresent(nnctrwspd10m_0)) usectrwspd10m=.false. if (.not.ctrvarpresent(nnuwnd10m_0) .or. .not.ctrvarpresent(nnvwnd10m_0)) use_uwnd10m_vwnd10m=.false. endif if (ncase.eq.2) then itime(1:6)=itime0(1:6) isign=0 ityped=0 typeproc(1:nvar2d)='analysis' endif if (ncase.eq.3) then itime(1:6)=itime0(1:6) isign=0 ityped=0 typeproc(1:nvar2d)='analysis error' varname(2:nvar2d)='anlerr-'//varname(2:nvar2d) endif if (ncase.eq.1) filename='bckg.grib2' if (ncase.eq.2) filename='anl.grib2' if (ncase.eq.3) filename='anlerr.grib2' !=================================================================================================== !==> write grib2 file !=================================================================================================== ! first create binary output for further processing if desired if ((lmapgrid .or. lhalfres) .and. mype==0) then open (41,file=trim(filename)//'.bin',form='unformatted') write(41) nvar2d write(41) varname,typeproc write(41) isign,ityped,ifcsthr write(41) itime endif allocate(field2(nx,ny)) do n=1,nvar2d field(:,:)=0. do j=js,je do i=is,ie field(i,j)=gndfd(i,j,n) enddo enddo field2=0. call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) if (mype==0) write(41) field2 enddo close(41) deallocate(field2) call write_ndfd_grib2(gndfd,is,ie,js,je,nx,ny,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthr,cgrid,cgprocess,filename,mype) if (ncase.eq.1) then if (lbiascor) then allocate(gndfd9(is:ie,js:je,nvar2d)) fname='slabs2_nobiasc.dat' call retrieve_gesanl(cgrid,nvar2d,n0,gndfd9, & fname,newvarname,ctrvarpresent,newvars,usectrwspd10m,use_uwnd10m_vwnd10m, & nngust,nnctrwspd10m,nnuwnd10m,nnvwnd10m,nnhowv_0,nnhowv) filename='bckg.grib2_nobiasc' call write_ndfd_grib2(gndfd9,is,ie,js,je,nx,ny,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthr,cgrid,cgprocess,filename,mype) deallocate(gndfd9) endif if (usefgat) then read(19,ptimeinfom1) close(19) if (mype==0) & print*,'in cnv_to_grib2,itimem1,ifcsthrm1,guessfcstm1=',& itimem1,ifcsthrm1,guessfcstm1 itime(1:6)=itimem1(1:6) if (guessfcstm1) then typeproc(1:nvar2d)='forecast' isign=1 ityped=1 else typeproc(1:nvar2d)='analysis' isign=0 ityped=0 endif allocate(gndfd9(is:ie,js:je,nvar2d)) fname='slabsm1.dat' call retrieve_gesanl(cgrid,nvar2d,n0,gndfd9, & fname,newvarname,ctrvarpresent,newvars,usectrwspd10m,use_uwnd10m_vwnd10m, & nngust,nnctrwspd10m,nnuwnd10m,nnvwnd10m,nnhowv_0,nnhowv) filename='bckg.grib2_m1' call write_ndfd_grib2(gndfd9,is,ie,js,je,nx,ny,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthrm1,cgrid,cgprocess,filename,mype) read(29,ptimeinfop1) close(29) if (mype==0) & print*,'in cnv_to_grib2,itimep1,ifcsthrp1,guessfcstp1=',& itimep1,ifcsthrp1,guessfcstp1 itime(1:6)=itimep1(1:6) typeproc(1:nvar2d)='forecast' isign=1 ityped=1 fname='slabsp1.dat' call retrieve_gesanl(cgrid,nvar2d,n0,gndfd9, & fname,newvarname,ctrvarpresent,newvars,usectrwspd10m,use_uwnd10m_vwnd10m, & nngust,nnctrwspd10m,nnuwnd10m,nnvwnd10m,nnhowv_0,nnhowv) filename='bckg.grib2_p1' call write_ndfd_grib2(gndfd9,is,ie,js,je,nx,ny,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthrp1,cgrid,cgprocess,filename,mype) deallocate(gndfd9) endif endif 100 continue deallocate(gndfd) deallocate(typeproc) deallocate(varname) !=================================================================================================== !==> deallocate fields !=================================================================================================== deallocate(field) call mpi_barrier(mpi_comm_world,ierror) ! !=================================================================================================== !==> if cgrid is 'cohresext' or 'cohreswexp', map output to true conus2p5 ndfd grid and nwrfc grid !=================================================================================================== if ((trim(cgrid)=='cohresext' .or. trim(cgrid)=='cohreswexp') .and. lmapgrid ) then if (trim(cgrid)=='cohresext') i0=0 if (trim(cgrid)=='cohreswexp') i0=i0wexp do ncase=1,3 ! bckg , analysis , anlerr if (ncase.eq.1) filename='bckg.grib2' if (ncase.eq.2) filename='anl.grib2' if (ncase.eq.3) filename='anlerr.grib2' open (41,file=trim(filename)//'.bin',form='unformatted') do nn=1,2 !map to cohres and nwrfc if (nn==1) then ix=nx_cohres jy=ny_cohres cgridnew='cohres' filename2=trim(filename)//'_map2cohres' cgridhalfres='conus' endif if (nn==2) then ix=nx_nwrfc jy=ny_nwrfc cgridnew='nwrfc' filename2=trim(filename)//'_map2nwrfc' cgridhalfres='nwrfc5p0' !never requested by NWRFC endif call horiz_domain_partition(ix,jy,mype,npe,is2,ie2,js2,je2) ! print*,'in cnv_to_grib2:filename2,is2,ie2,js2,je2=', & ! trim(filename2),is2,ie2,js2,je2 allocate(gndfd9(is2:ie2,js2:je2,nvar2d)) allocate(field(nx,ny)) rewind(41) read(41) nvar2d allocate(typeproc(nvar2d)) allocate(varname(nvar2d)) read(41) varname,typeproc read(41) isign,ityped,ifcsthr read(41) itime do k=1,nvar2d read(41) field if (nn==1) then do j=js2,je2 do i=is2,ie2 gndfd9(i,j,k)=field(i0+i,j) enddo enddo endif if (nn==2) then j0=ny-jy do j=js2,je2 do i=is2,ie2 gndfd9(i,j,k)=field(i0+i,j0+j) !note that js2=1 and je2=ny_nwrfc enddo enddo endif enddo if (nn==2) then if (mype==0) print*,'in cnv_to_grib2: for nwrfc grid, j0=',j0 endif call write_ndfd_grib2(gndfd9,is2,ie2,js2,je2,ix,jy,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthr,cgridnew,cgprocess,filename2,mype) if (nn==1 .and. lhalfres) then !note: no need for half grid for nwrfc mx=(ix-1+0.0001)/2+1 my=(jy-1+00.001)/2+1 if (mype==0) print*,'in cnv_to_grib2: & &dimensions for half-resolution grid ', & trim(cgridhalfres),' are mx,my=',mx,my call horiz_domain_partition(mx,my,mype,npe,is3,ie3,js3,je3) ! print*,'in cnv_to_grib2: cgridhalfres,is3,ie3,js3,je3=', & ! trim(cgridhalfres),is3,ie3,js3,je3 allocate(gndfd(is3:ie3,js3:je3,nvar2d)) call resolution_halfing(gndfd9,ix,jy,nvar2d,is2,ie2,js2,je2, & gndfd,mx,my,is3,ie3,js3,je3,varname,iopt,iopthgt, & smcf1,smcf2,one2one,npass,npasshgt,mype) filename3=trim(filename2)//'_halfres' call write_ndfd_grib2(gndfd,is3,ie3,js3,je3,mx,my,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthr,cgridhalfres,cgprocess,filename3,mype) deallocate(gndfd) endif deallocate(typeproc) deallocate(varname) deallocate(gndfd9) deallocate(field) enddo close(41) enddo endif ! !=================================================================================================== !==> generate half resolution analysis when cgrid is cohres or akhres. ! note that for cgrid=cohresext and cohreswexp it's taken care of in the previous block !=================================================================================================== if ((trim(cgrid)=='cohres' .or. trim(cgrid)=='akhres') .and. lhalfres) then if (trim(cgrid)=='cohres') cgridhalfres='conus' if (trim(cgrid)=='akhres') cgridhalfres='alaska' do ncase=1,3 ! bckg or analysis , anlerr if (ncase.eq.1) filename='bckg.grib2' if (ncase.eq.2) filename='anl.grib2' if (ncase.eq.3) filename='anlerr.grib2' open (41,file=trim(filename)//'.bin',form='unformatted') read(41) nvar2d allocate(typeproc(nvar2d)) allocate(varname(nvar2d)) read(41) varname,typeproc read(41) isign,ityped,ifcsthr read(41) itime allocate(gndfd9(is:ie,js:je,nvar2d)) allocate(field(nx,ny)) do k=1,nvar2d read(41) field do j=js,je do i=is,ie gndfd9(i,j,k)=field(i,j) enddo enddo enddo close(41) mx=(nx-1+0.0001)/2+1 my=(ny-1+00.001)/2+1 if (mype==0) print*,'in cnv_to_grib2: & &dimensions for half-resolution grid ', & trim(cgridhalfres),' are mx,my=',mx,my call horiz_domain_partition(mx,my,mype,npe,is3,ie3,js3,je3) ! print*,'in cnv_to_grib2: cgridhalfres,is3,ie3,js3,je3=', & ! trim(cgridhalfres),is3,ie3,js3,je3 allocate(gndfd(is3:ie3,js3:je3,nvar2d)) call resolution_halfing(gndfd9,nx,ny,nvar2d,is,ie,js,je, & gndfd,mx,my,is3,ie3,js3,je3,varname,iopt,iopthgt, & smcf1,smcf2,one2one,npass,npasshgt,mype) filename=trim(filename)//'_halfres' call write_ndfd_grib2(gndfd,is3,ie3,js3,je3,mx,my,nvar2d,& varname,typeproc,isign,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6),& ityped,ifcsthr,cgridhalfres,cgprocess,filename,mype) deallocate(field) deallocate(gndfd9) deallocate(typeproc) deallocate(varname) deallocate(gndfd) enddo endif deallocate(cldchmap) end subroutine cnv_to_grib2_mpi !=================================================================================================== !=================================================================================================== subroutine retrieve_gesanl(cgrid,nvar2d,n0,g, & fname,newvarname,ctrvarpresent,newvars,usectrwspd10m,use_uwnd10m_vwnd10m, & nngust,nnctrwspd10m,nnuwnd10m,nnvwnd10m,nnhowv_0,nnhowv) use mpi use mpitaskmod, only: mype,npe use mpitaskmod, only: is => ista use mpitaskmod, only: ie => iend use mpitaskmod, only: js => jsta use mpitaskmod, only: je => jend use mpitaskmod, only: nx => nlon use mpitaskmod, only: ny => nlat implicit none !Declare passed variables character(60),intent(in):: cgrid integer(4),intent(in):: nvar2d,n0,newvars integer(4),intent(in):: nngust,nnctrwspd10m integer(4),intent(in):: nnuwnd10m,nnvwnd10m integer(4),intent(in):: nnhowv_0,nnhowv character(30),intent(in):: fname character(60),intent(in):: newvarname(newvars) logical,intent(in)::ctrvarpresent(newvars) logical,intent(in):: usectrwspd10m logical,intent(in):: use_uwnd10m_vwnd10m real(4),intent(out)::g(is:ie,js:je,nvar2d) !Declare local parameters integer(4),parameter:: mndelmax=5 !Declare local variables real(4), allocatable,dimension(:,:),save::td0,td1,fis0,glon real(4), allocatable,dimension(:,:):: field,field2 real(4), allocatable,dimension(:,:):: td2 real(4), allocatable,dimension(:,:):: ue,ve real(4), allocatable,dimension(:,:):: ufield,vfield,wfield integer(4), parameter:: ioan1=10 real(4),parameter::qmin=1.e-06 integer(4) n,nn,ierror integer(4) m,m1,m2,n1,n2,mndel integer(4) itime0(6),itime(6),nlon,nlat,nsig integer(4) istat real(4) u0,v0,w0,wdir0,w2 real(4) u00,v00,w00 real(4) amin1,amin2 real(4) amax1,amax2 integer(4) i,j logical fexist !=================================================================================================== print*,'in retrieve_gesanl: fname=',fname if (.not. allocated(td0)) allocate(td0(is:ie,js:je)) if (.not. allocated(td1)) allocate(td1(is:ie,js:je)) if (.not. allocated(fis0)) allocate(fis0(is:ie,js:je)) if (.not. allocated(glon)) allocate(glon(is:ie,js:je)) allocate( field (nx,ny) ) allocate( field2(nx,ny) ) allocate( td2 (is:ie,js:je) ) allocate( ue (is:ie,js:je) ) allocate( ve (is:ie,js:je) ) ! !=================================================================================================== !==> readin td0 !=================================================================================================== if (trim(fname)=='sigges') then open (11,file='slabs2.dat',form='unformatted') do n=1,4 read(11) field if (n==2) fis0(is:ie,js:je)=field(is:ie,js:je) !saved field if (n==4) td0(is:ie,js:je)=field(is:ie,js:je) !saved field enddo close(11) endif !=================================================================================================== !==> read in guess and analysis fields !=================================================================================================== open (ioan1,file=trim(fname),form='unformatted') if (trim(fname)=='sigges'.or.trim(fname)=='siganl') then read(ioan1) itime0,nlon,nlat,nsig read(ioan1) field,field2 !glat,dx read(ioan1) field,field2 !glon,dy if (trim(fname)=='sigges') glon (is:ie,js:je) = field(is:ie,js:je) if (mype==0) then print*,'in retrieve_gesanl, fname, itime0=',trim(fname),itime0 print*,'in retrieve_gesanl, fanme, nlon,nlat,nsig=', & trim(fname),nlon,nlat,nsig endif read(ioan1) field ; g(is:ie,js:je,2) = field(is:ie,js:je) !psfc read(ioan1) field ; g(is:ie,js:je,1) = fis0(is:ie,js:je) !field is geopght !use fis0, which is terrain read(ioan1) field ; g(is:ie,js:je,3) = field(is:ie,js:je) !t read(ioan1) field ; g(is:ie,js:je,7) = field(is:ie,js:je) !q if (trim(fname)=='siganl') then do j=js,je do i=is,ie g(i,j,7)=max(qmin,g(i,j,7)) enddo enddo endif !--------------------------------------------------------------------- !psfc, q, t, td call get_dewpt(g(:,:,2),g(:,:,7),g(:,:,3),g(:,:,4),is,ie,js,je) amin1=minval(g(:,:,4)) amax1=maxval(g(:,:,4)) call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'in retrieve_gesanl: fname, td,min,max=',trim(fname),amin2,amax2 !--------------------------------------------------------------------- if (trim(fname)=='sigges') then td1(:,:)=g(:,:,4) !used to refine td analysis g(:,:,4)=td0(:,:) !this is the actual field used as background for td elseif (trim(fname)=='siganl') then call refine_tdanl(td0,td1,g(:,:,3),g(:,:,4),is,ie,js,je,nx,ny,mype) endif read(ioan1) field ; g(is:ie,js:je,5) = field(is:ie,js:je) !u read(ioan1) field ; g(is:ie,js:je,6) = field(is:ie,js:je) !v !--------------------------------------------------------------------- !u, v, wspd, wdir call wspdwdir(cgrid,glon,g(:,:,5),g(:,:,6),ue,ve,g(:,:,9),g(:,:,8), & is,ie,js,je,mype) !--------------------------------------------------------------------- !must jump records to get to gust do n=1,11 read(ioan1) enddo if (trim(fname)=='siganl') read(ioan1) !jump one more record n=n0 do nn=1,newvars read(ioan1) field !field always present even if ctr var isn't being used if (ctrvarpresent(nn)) then if (mype==0) print*,'in retrieve_gesanl: fname,nn, ctrvarpresent(nn),field,min,max=',& &trim(fname),nn, ctrvarpresent(nn),minval(field),maxval(field) n=n+1 g(is:ie,js:je,n)=field(is:ie,js:je) endif enddo if (mype==0) print*,'in retrieve_gesanl: n0,nnctrwspd10m,nnuwnd10m,nnvwnd10m=', & n0,nnctrwspd10m,nnuwnd10m,nnvwnd10m !========================================================================================= !==> if requested, scale u and v and swap (psi,chi)-based wind speed with ctrwspd10m !========================================================================================= if (usectrwspd10m) then allocate( ufield (nx,ny) ) allocate( vfield (nx,ny) ) allocate( wfield (nx,ny) ) !--------------------------------------- field=0. if (use_uwnd10m_vwnd10m) then field(is:ie,js:je)=g(is:ie,js:je,n0+nnuwnd10m) else field(is:ie,js:je)=g(is:ie,js:je,5) endif call mpi_allreduce(field,ufield,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) !--------------------------------------- !--------------------------------------- field=0. if (use_uwnd10m_vwnd10m) then field(is:ie,js:je)=g(is:ie,js:je,n0+nnvwnd10m) else field(is:ie,js:je)=g(is:ie,js:je,6) endif call mpi_allreduce(field,vfield,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) !--------------------------------------- !--------------------------------------- field=0. if (use_uwnd10m_vwnd10m) then do j=js,je do i=is,ie field(i,j)=g(i,j,n0+nnuwnd10m)*g(i,j,n0+nnuwnd10m)+ & g(i,j,n0+nnvwnd10m)*g(i,j,n0+nnvwnd10m) if (field(i,j).ne.0.0) field(i,j)=sqrt(field(i,j)) enddo enddo else field(is:ie,js:je)=g(is:ie,js:je,9) endif call mpi_allreduce(field,wfield,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) !--------------------------------------- !for diagnostic purposes only, load (psi,chi)-based wind fields !and the independently analyzed u and v fields g(is:ie,js:je,nvar2d-2)=g(is:ie,js:je,5) !usfc g(is:ie,js:je,nvar2d-1)=g(is:ie,js:je,6) !vsfc g(is:ie,js:je,nvar2d )=g(is:ie,js:je,8) !wdirsfc !now determine the final u and v fields do j=js,je do i=is,ie !swap (psi,chi) based wind speed with ctrwspd10m w2=g(i,j,n0+nnctrwspd10m) g(i,j,n0+nnctrwspd10m)=g(i,j,9) g(i,j,9)=w2 !prepare to rescale u,v if (use_uwnd10m_vwnd10m) then u0=g(i,j,n0+nnuwnd10m) v0=g(i,j,n0+nnvwnd10m) else u0=g(i,j,5) v0=g(i,j,6) endif w0=u0*u0+v0*v0 ; if (w0.ne.0.0) w0=sqrt(w0) !rescale u,v g(i,j,5)=u0 g(i,j,6)=v0 if (w0>0.) then g(i,j,5)=(u0/w0)*w2 g(i,j,6)=(v0/w0)*w2 elseif(w0==0 .and. w2 >0.) then w00=1000. do mndel=1,mndelmax m1=max(1,i-mndel) ; m2=min(nx,i+mndel) n1=max(1,j-mndel) ; n2=min(ny,j+mndel) do n=n1,n2 do m=m1,m2 if (wfield(m,n)>0 .and. wfield(m,n) analyzed gust cannot be less than wind speed !========================================================================================= if (trim(fname)=='siganl' .and. ctrvarpresent(nngust)) then do j=js,je do i=is,ie if (g(i,j,n0+nngust) < g(i,j,9)) g(i,j,n0+nngust)=g(i,j,9) enddo enddo endif !========================================================================================= !==> remove residual amplitude from howv analysis !========================================================================================= if (mype==0) print*,'in retrieve_gesanl: nnhowv_0,ctrvarpresent(nnhowv_0)=', & nnhowv_0,ctrvarpresent(nnhowv_0) if ((trim(fname)=='sigges' .or. trim(fname)=='siganl') .and. ctrvarpresent(nnhowv_0)) then open (94,file='rtma_slmask_howv.dat',form='unformatted') read(94) field close(94) do j=js,je do i=is,ie if (field(i,j) >= 0.5) g(i,j,n0+nnhowv)=9.999e+20 enddo enddo endif elseif (trim(fname)=='anlerr') then inquire(file='bckg_z.dat',exist=fexist) if (fexist) then deallocate(field2) allocate(field2(ny,nx)) !transposed open(94,file='bckg_z.dat',form='unformatted') read(94) field2 field=transpose(field2) g(is:ie,js:je,1)=field(is:ie,js:je) close(94) deallocate(field2) else g(:,:,1)=fis0(:,:) endif open (53,file='errfield.dat',form='unformatted') do n=2,9 !perr,terr,tderr,uerr,verr,qerr,wdirerr2,wspderr read(53) field g(is:ie,js:je,n)=field(is:ie,js:je) enddo read(53) field !skip uerr2 read(53) field !skip verr2 read(53) field !skip wdirerr do n=n0+1,nvar2d !gusterr,viserr,pblherr,disterr, ... lcbas,uwnd10m,vwnd10m read(53,iostat=istat) field !field only present if ctr var is being used if (istat/=0) exit g(is:ie,js:je,n)=field(is:ie,js:je) enddo close(53) if (usectrwspd10m) then do j=js,je do i=is,ie g(i,j,nvar2d-2)=g(i,j,5) !anlerr-usfc g(i,j,nvar2d-1)=g(i,j,6) !anlerr-vsfc g(i,j,nvar2d )=g(i,j,8) !anlerr-wdirsfc field(i,j)=g(i,j,9) g(i,j,9)=g(i,j,n0+nnctrwspd10m) g(i,j,n0+nnctrwspd10m)=field(i,j) if (use_uwnd10m_vwnd10m) then g(i,j,5)=g(i,j,n0+nnuwnd10m) g(i,j,6)=g(i,j,n0+nnvwnd10m) endif enddo enddo endif if (ctrvarpresent(nnhowv_0)) then open (94,file='rtma_slmask_howv.dat',form='unformatted') read(94) field close(94) do j=js,je do i=is,ie if (field(i,j) >= 0.5) g(i,j,n0+nnhowv)=9.999e20 enddo enddo endif elseif (trim(fname)=='slabs2_nobiasc.dat' .or. & trim(fname)=='slabsm1.dat' .or. & trim(fname)=='slabsp1.dat') then read(ioan1) field ; g(is:ie,js:je,2) = field(is:ie,js:je) !psfc read(ioan1) field ; g(is:ie,js:je,1) = field(is:ie,js:je) !fis read(ioan1) field ; g(is:ie,js:je,3) = field(is:ie,js:je) !t read(ioan1) field ; g(is:ie,js:je,4) = field(is:ie,js:je) !td read(ioan1) field ; g(is:ie,js:je,5) = field(is:ie,js:je) !u read(ioan1) field ; g(is:ie,js:je,6) = field(is:ie,js:je) !v read(ioan1) field ; g(is:ie,js:je,7) = field(is:ie,js:je) !q call wspdwdir(cgrid,glon,g(:,:,5),g(:,:,6), & ue,ve,g(:,:,9),g(:,:,8),is,ie,js,je,mype) n=n0 do nn=1,newvars read(ioan1) field !field always present even if ctr var isn't being used if (ctrvarpresent(nn)) then print*,'in retrieve_gesanl: fname,nn, ctrvarpresent(nn),field,min,max=',& &trim(fname),nn, ctrvarpresent(nn),minval(field),maxval(field) n=n+1 g(is:ie,js:je,n)=field(is:ie,js:je) endif enddo if (usectrwspd10m) then g(is:ie,js:je,nvar2d-2)=g(is:ie,js:je,5) !usfc g(is:ie,js:je,nvar2d-1)=g(is:ie,js:je,6) !vsfc g(is:ie,js:je,nvar2d )=g(is:ie,js:je,8) !wdirsfc endif else print*,'unknown fname ... aborting in subroutine retrieve_gesanl' call mpi_finalize(ierror) stop endif close(ioan1) deallocate(field) if (allocated(field2)) deallocate(field2) deallocate(td2) deallocate(ue) deallocate(ve) end subroutine retrieve_gesanl !=================================================================================================== !=================================================================================================== subroutine refine_tdanl(td0,td1,t2,td2,is,ie,js,je,nx,ny,mype) use mpi implicit none !Declare passed variables integer(4),intent(in):: is,ie,js,je,nx,ny,mype real(4),dimension(is:ie,js:je),intent(in)::td0,td1,t2 real(4),dimension(is:ie,js:je),intent(inout)::td2 !Declare local variables integer(4) i,j,ierror real(4) amin1,amin2 real(4) amax1,amax2 real(4),allocatable,dimension(:,:)::field,field2 ! !=================================================================================================== !==> recompute the analysis dew point to be the guess dew point td0 plus the dew point increment, ! whereby the increment is the difference between the calculated analysis dew point td2 (via ! subroutine get_dewpt) and the calculated guess dew point td1. A two-pass smoothing is applied ! to that difference field !=================================================================================================== ! allocate(field(nx,ny)) allocate(field2(nx,ny)) !dew point increment do j=js,je do i=is,ie td2(i,j) = max ( -20., td2(i,j)-td1(i,j) ) enddo enddo amin1=minval(td2(:,:)) amax1=maxval(td2(:,:)) call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'in refine_tdanl: td2 inc/min,max=',amin2,amax2 !----------------------------------------------------- field(:,:)=0. do j=js,je do i=is,ie field(i,j)=td2(i,j) enddo enddo field2(:,:)=0. call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) ! call smther_one(field2,1,nx,1,ny,2) do j=js,je do i=is,ie td2(i,j)=field2(i,j) enddo enddo do j=js,je do i=is,ie td2(i,j) = min ( t2(i,j), td2(i,j)+td0(i,j) ) !td1 enddo enddo amin1=minval(td2(is:ie,js:je)) amax1=maxval(td2(is:ie,js:je)) call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror) call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror) if (mype==0) print*,'in refine_tdanl: td2 anl/min,max=',amin2,amax2 deallocate(field) deallocate(field2) end subroutine refine_tdanl !=================================================================================================== !=================================================================================================== !*************************************************************************************************** subroutine resolution_halfing(gin,nx,ny,nk,is,ie,js,je, & gout,mx,my,m1,m2,n1,n2,varname,iopt,iopthgt, & smcf1,smcf2,one2one,npass,npasshgt,mype) !*************************************************************************************************** !abstract: half the resolution of the input hres grid ! iopt=0 ==> select every other grid point of hres grid ! iopt=1 ==> 0.5x+0.125*(values at four nearest grid points of hres grid) ! iopt=2 ==> 0.5x+0.125*(values at four vertice points of ! enclosing square of hres grid) ! note: x is the hres field value at the targetted grid point ! note: for wdir10m, code always used iopt=0 ! ! iopthgt: for terrain height, override iopt with iopthgt ! ! smcf1,smcf2 : coefficients used with smoother-desmoother ! one2one : if .true. perform 1-2-1 smoothing rather than smoothing-desmoothing ! npass : number of soothing passes ! npasshgt : for terrain height, override npass with npasshgt !*************************************************************************************************** use mpi implicit none !Declare passed variables integer(4),intent(in)::nx,ny,nk,is,ie,js,je integer(4),intent(in)::mx,my,m1,m2,n1,n2 integer(4),intent(in)::iopt,iopthgt,npass,npasshgt integer(4),intent(in)::mype real(4),intent(in)::gin(is:ie,js:je,nk) real(4),intent(in)::smcf1,smcf2 character(60),intent(in)::varname(nk) logical,intent(in)::one2one real(4),intent(out)::gout(m1:m2,n1:n2,nk) !Declare local variables integer(4) iopt0,npass0 integer(4) i,j,i2,j2,k integer(4) im,ip,jm,jp integer(4) ierror real(4),allocatable,dimension(:,:):: hfine,hfine2 real(4),allocatable,dimension(:,:):: hcoarse,hcoarse2 !--------------------------------------------------------------------------------------------------- allocate(hfine(nx,ny)) allocate(hfine2(nx,ny)) allocate(hcoarse(mx,my)) allocate(hcoarse2(mx,my)) do k=1,nk if (mype==0) print*,'in resolution_halfing: k, varname=',k,varname(k) hfine2(:,:)=0. do j=js,je do i=is,ie hfine2(i,j)=gin(i,j,k) enddo enddo call mpi_allreduce(hfine2,hfine,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) hcoarse(:,:)=0. if (trim(varname(k))=='hgt') then iopt0=iopthgt else if(trim(varname(k))=='wdir10m') then iopt0=0 !won't try to average wdir else if(trim(varname(k))=='cldch') then iopt0=0 !won't try to average cldch due to bitmap else iopt0=iopt endif if (iopt0.ne.0 .and. iopt0.ne.1 .and. iopt0.ne.2) then if (mype==0) then print*,'in resolution_halfing: warning ... invalid iopt0=',iopt0,' value' print*,'in resolution_halfing: reset to iopt0=1' endif iopt0=1 endif do j2=n1,n2 j=2*j2-1 jm=max(j-1,1) jp=min(j+1,ny) do i2=m1,m2 i=2*i2-1 im=max(i-1,1) ip=min(i+1,nx) if (iopt0==0) then hcoarse(i2,j2)=hfine(i,j) elseif (iopt0==1) then hcoarse(i2,j2)=0.5*hfine(i,j) + & 0.125*(hfine(i,jm)+hfine(i,jp)+hfine(im,j)+hfine(ip,j)) elseif (iopt0==2) then hcoarse(i2,j2)=0.5*hfine(i,j) + & 0.125*(hfine(im,jm)+hfine(ip,jm)+hfine(im,jp)+hfine(ip,jp)) endif enddo enddo if (trim(varname(k))=='hgt') then npass0=npasshgt else if(trim(varname(k))=='wdir10m') then npass0=0 else npass0=npass endif call mpi_allreduce(hcoarse,hcoarse2,mx*my,mpi_real4,mpi_sum,mpi_comm_world,ierror) call smoother_desmoother(hcoarse2,1,mx,1,my,smcf1,smcf2,one2one,npass0) do j=n1,n2 do i=m1,m2 gout(i,j,k)=hcoarse2(i,j) enddo enddo enddo deallocate(hfine) deallocate(hfine2) deallocate(hcoarse) deallocate(hcoarse2) return end !*************************************************************** subroutine get_earth_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, & olon,dg2rad) ! ! Abstract: Rotates wind relative to grid on Lambert Conformal ! or polar streo projection to get winds with respect ! to the true north. ! Adapted from Stan Benjamin's code snippet implicit none real(4) un,vn,ut,vt,ROTCON_P,LON_XX_P,olon, & dg2rad,angle2,sinx2,cosx2 !** ROTCON_P R WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO !** AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL !** LAT_TAN_P R LATITUDE AT LAMBERT CONFORMAL PROJECTION !** IS TRUE (DEG) !** LON_XX_P R MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG) angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad sinx2 = sin(angle2) cosx2 = cos(angle2) un = cosx2*ut+sinx2*vt vn =-sinx2*ut+cosx2*vt return end !*************************************************************** subroutine get_grid_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, & olon,dg2rad) ! ! Abstract: Rotates wind relative to true north to ! get winds with respect to grid on Lambert Conformal ! or polar streo projection. ! Adapted from Stan Benjamin's code snippet which performs the ! inverse operation implicit none real(4) un,vn,ut,vt,ROTCON_P,LON_XX_P,olon, & dg2rad,angle2,sinx2,cosx2 !** ROTCON_P R WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO !** AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL !** LAT_TAN_P R LATITUDE AT LAMBERT CONFORMAL PROJECTION !** IS TRUE (DEG) !** LON_XX_P R MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG) angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad sinx2 = sin(angle2) cosx2 = cos(angle2) ut = cosx2*un-sinx2*vn vt = sinx2*un+cosx2*vn return end !******************************************************************************* !******************************************************************************* subroutine wspdwdir(cgrid,glon,u,v,ue,ve,w,wd,i1,i2,j1,j2,mype) implicit none ! Declare passed variables character(60),intent(in):: cgrid integer(4),intent(in)::i1,i2,j1,j2 integer(4),intent(in)::mype real(4),dimension(i1:i2,j1:j2),intent(in) :: glon,u,v real(4),dimension(i1:i2,j1:j2),intent(out):: ue,ve,w,wd ! Declare local variables real(8) alat18,elon18,dx8,elonv8,alatan8 real(4) dg2rad,elonv4,alatan4,xn real(4) ylon integer(4) i,j logical lambconform logical polarstereo call proj_info(cgrid,dx8,alat18,elon18,elonv8,alatan8) if (mype==0) then print*,'in wspdwdir: cgrid=',trim(cgrid) print*,'in wspdwdir: alat18=',alat18 print*,'in wspdwdir: elon18=',elon18 print*,'in wspdwdir: dx8=',dx8 print*,'in wspdwdir: elonv8=',elonv8 print*,'in wspdwdir: alatan8=',alatan8 endif lambconform=trim(cgrid)=='conus'.or.trim(cgrid)=='cohres'.or. & trim(cgrid)=='hrrr'.or.trim(cgrid)=='cohresext'.or.trim(cgrid)=='cohreswexp' polarstereo=trim(cgrid)=='alaska'.or.trim(cgrid)=='akhres'.or. & trim(cgrid)=='juneau' dg2rad=atan(1.)/45. elonv4=elonv8 alatan4=alatan8 xn=1. if (lambconform) xn=sin(alatan4*dg2rad) if (mype==0) print*,'in wspdwdir: lambconform,polarstereo,xn=', & lambconform,polarstereo,xn if (lambconform .or. polarstereo) then do j=j1,j2 do i=i1,i2 ylon=glon(i,j)/dg2rad !CHECK SIGN OF GLON FOR CONUS + ALASKA call get_earth_winds(ue(i,j),ve(i,j),u(i,j),v(i,j),xn,elonv4-360.,ylon,dg2rad) enddo enddo else ue(i1:i2,j1:j2)=u(i1:i2,j1:j2) ve(i1:i2,j1:j2)=v(i1:i2,j1:j2) endif do j=j1,j2 do i=i1,i2 call speeddir(ue(i,j),ve(i,j),wd(i,j),w(i,j)) enddo enddo return end !******************************************************************************* !******************************************************************************* subroutine speeddir(ue,ve,wdir,wspd) ! implicit none !Declare passed variables real(4),intent(in):: ue,ve real(4),intent(out):: wdir,wspd !Declare local variables real(4) dg2rad,aearth dg2rad=atan(1.)/45. wspd=ue*ue+ve*ve if (wspd.ne.0.0) wspd=sqrt(wspd) if (wspd .eq. 0.0) then wdir=0.! calm !99999.0 else if (ve.eq.0.0) then if (ue.gt.0.0) wdir = 270. if (ue.lt.0.0) wdir = 90. else aearth = atan(ue/ve)/dg2rad if (ue.le.0.0 .and. ve.le.0.0 ) wdir = aearth if (ue.le.0.0 .and. ve.ge.0.0 ) wdir = aearth + 180.0 if (ue.ge.0.0 .and. ve.ge.0.0 ) wdir = aearth + 180.0 if (ue.ge.0.0 .and. ve.le.0.0 ) wdir = aearth + 360.0 endif endif return end !******************************************************************************* !******************************************************************************* subroutine get_dewpt(p,q,t,td,is,ie,js,je) ! !==> input variables ! q-specific humidity ; p-pressure in Pa ! t-temperature in K !==> output variable ! td-dewpoint in K implicit none integer(4),intent(in):: is,ie,js,je real(4),dimension(is:ie,js:je),intent(in)::q,p,t real(4),dimension(is:ie,js:je),intent(out)::td real(4),parameter::eps=0.62197 !=Rd/Rv integer(4) i,j real(4) e, qv, eln do j=js,je do i=is,ie qv=q(i,j)/(1.-q(i,j)) e=p(i,j)/100.*qv/(eps+qv) eln=alog(e) td(i,j) = (243.5*eln-440.8)/(19.48-eln)+273.15 td(i,j) = min(t(i,j),td(i,j)) enddo enddo return end !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine write_ndfd_grib2(g,i1,i2,j1,j2,nx,ny,nflds,varname,typeproc, & isign,y0,m0,d0,h0,min0,s0,ityped,ifcsthr,clocation,cgprocess,filename,mype) !****************************************************************** ! prgmmr: pondeca org: np20 date: 2006-03-03 * ! * ! abstract: * ! use steve gilbert's gribcreate, addgrid, addfield, and gribend * ! subroutines to write NDFD files in Grib2 format * ! * ! program history log: * ! 2006-03-03 pondeca * ! * ! input argument list: * ! * ! 1. g(i1:i2,j1:j2,nflds): floating point array that contains all * ! the nflds fields on the mpi subdomains of dimensions * ! (i1:i2,j1:j2). The global dimensions are (nx,ny) * ! * ! 2. varname(nflds): character array that conatains the names of * ! all the nflds fields from g(nx,ny,nflds) * ! * ! 3. typeproc(nflds): character array that contains process types * ! (analysis, forecast, or analysis error) * ! * ! 4. y0,m0,d0,h0,min0,s0: reference time (year, month, day, hour, * ! minutes and seconds) * ! * ! 5. isign: significance of reference time. Values of interest to * ! NDFD are 0 for analysis and 1 for start of forecast. * ! * ! 6. ityped: type of processed data. Values of interest to * ! NDFD are 0 for analysis products and 1 for forecast products * ! * ! 7. ifcsthr: forecast hour from reference time. use negative * ! value for forecast time in minutes * ! * ! 8. clocation: string variable. denotes grid name. supported * ! grids are "conus", "alaska" , "hawaii", "guam", "prico". * ! "cohres", "akhres" , "hrrr" , "juneau" , "cohresext" , * ! "conusext" , and "nwrfc" * ! * ! 9. cgprocess: generating process: rtma or urma * ! * ! output argument list: * ! 1. filename: Name of output Grib2 file * ! * ! attributes: * ! language: f90 * ! machine: ibm RS/6000 SP * !****************************************************************** use mpi implicit none integer(4),intent(in):: mype integer(4),intent(in):: i1,i2,j1,j2 integer(4),intent(in):: nx,ny,nflds,isign, & y0,m0,d0,h0,min0,s0,ityped,ifcsthr integer(4), parameter::max_bytes=200*130000 integer(4), parameter::idefnum=1 ! integer(4), parameter::ipdstmplen=15 integer(4), parameter::idrstmplen=5 real(4),parameter::cldch_bitmapval=1.e+21 integer(4) y1_minmaxt,m1_minmaxt,d1_minmaxt, & h1_minmaxt,min1_minmaxt,s1_minmaxt, & ifcsthr_minmaxt integer(4) y2_minmaxt,m2_minmaxt,d2_minmaxt, & h2_minmaxt,min2_minmaxt,s2_minmaxt integer(4) listsec0(2) integer(4) listsec1(13) integer(4) ierr,n,i,j,ij integer(4) lengrib,dscal integer(4) lunout integer(4) igds(5),igdstmplen integer(4) ideflist(idefnum) ! integer(4) ipdstmpl(ipdstmplen) integer(4) ipdstmplen integer(4), allocatable,dimension(:):: ipdstmpl integer(4) idrstmpl(idrstmplen) integer(4) idrsnum,ibmap,numcoord,ipdsnum integer(4),allocatable,dimension(:):: igdstmpl character*60 filename character*60 varname(nflds) character*60 typeproc(nflds) character*1 cgrib(max_bytes) character*60 clocation character*60 cgprocess logical*1 bmap(nx,ny) logical fexist logical lminmaxt real(4) g(i1:i2,j1:j2,nflds) real(4) coordlist(1) real(4),allocatable,dimension(:):: fld real(4),allocatable,dimension(:,:):: field,field2 integer(4) ierror namelist/minmaxt_range/y1_minmaxt,m1_minmaxt,d1_minmaxt, & h1_minmaxt,min1_minmaxt,s1_minmaxt, & ifcsthr_minmaxt, & y2_minmaxt,m2_minmaxt,d2_minmaxt, & h2_minmaxt,min2_minmaxt,s2_minmaxt allocate(fld(nx*ny)) allocate(field(nx,ny)) allocate(field2(nx,ny)) if (mype==0) print*,'---- starting write_ndfd_grib2 -------' !==>determine ndfd grid location and set igdstmplen dimension if (mype==0) print*,'grid name is=',trim(clocation) if (mype==0) print*,'generating process is=',trim(cgprocess) if (trim(clocation)=="puerto_rico") clocation="prico" if (trim(clocation) /= "conus" .and. & trim(clocation) /= "alaska" .and. & trim(clocation) /= "hawaii" .and. & trim(clocation) /= "guam" .and. & trim(clocation) /= "prico" .and. & trim(clocation) /= "cohres" .and. & trim(clocation) /= "akhres" .and. & trim(clocation) /= "hrrr" .and. & trim(clocation) /= "juneau" .and. & trim(clocation) /= "cohresext" .and. & trim(clocation) /= "conusext" .and. & trim(clocation) /= "cohreswexp" .and. & trim(clocation) /= "nwrfc" ) then if (mype==0) print*,'grid not supported. aborting ...' call abort endif y1_minmaxt=y0 ; y2_minmaxt=y0 m1_minmaxt=m0 ; m2_minmaxt=m0 d1_minmaxt=d0 ; d2_minmaxt=d0 h1_minmaxt=h0 ; h2_minmaxt=h0 min1_minmaxt=min0 ; min2_minmaxt=min0 s1_minmaxt=s0 ; s2_minmaxt=s0 ifcsthr_minmaxt=ifcsthr inquire(file='minmaxt_range_input',exist=fexist) if(fexist) then open (55,file='minmaxt_range_input',form='formatted') read (55,minmaxt_range) close(55) endif if (mype.eq.0) then print*,'y1_minmaxt,m1_minmaxt,d1_minmaxt,h1_minmaxt,min1_minmaxt,s1_minmaxt=',& y1_minmaxt,m1_minmaxt,d1_minmaxt,h1_minmaxt,min1_minmaxt,s1_minmaxt ; print* print*,'ifcsthr_minmaxt=',ifcsthr_minmaxt ; print* print*,'y2_minmaxt,m2_minmaxt,d2_minmaxt,h2_minmaxt,min2_minmaxt,s2_minmaxt=',& y2_minmaxt,m2_minmaxt,d2_minmaxt,h2_minmaxt,min2_minmaxt,s2_minmaxt ; print* endif if (trim(clocation)=="conus") igdstmplen=22 if (trim(clocation)=="alaska") igdstmplen=18 if (trim(clocation)=="hawaii") igdstmplen=19 if (trim(clocation)=="guam") igdstmplen=19 if (trim(clocation)=="prico" ) igdstmplen=19 if (trim(clocation)=="cohres") igdstmplen=22 if (trim(clocation)=="akhres") igdstmplen=18 if (trim(clocation)=="hrrr" ) igdstmplen=22 if (trim(clocation)=="juneau") igdstmplen=18 if (trim(clocation)=="cohresext") igdstmplen=22 if (trim(clocation)=="conusext") igdstmplen=22 if (trim(clocation)=="cohreswexp") igdstmplen=22 if (trim(clocation)=="nwrfc") igdstmplen=22 allocate(igdstmpl(igdstmplen)) !==>open file that will store grib2 messages if (mype==0) then call getlun90(lunout,1) print*,'about to open',lunout,trim(filename) call baopenw(lunout,trim(filename),ierr) print*,'write_ndfd_grib2: opened ',lunout, & 'for grib2 data ',trim(filename), & 'return code is ',ierr endif do n=1,nflds lminmaxt=trim(varname(n)).eq.'mxtm'.or.trim(varname(n)).eq.'anlerr-mxtm'.or. & trim(varname(n)).eq.'mitm'.or.trim(varname(n)).eq.'anlerr-mitm' !==>initialize new GRIB2 message and pack ! GRIB2 sections 0 (Indicator Section) and 1 (Identification ! Section) listsec0(1)=0 ! Discipline-GRIB Master Table Number/ Meteorological products (see Code Table 0.0) if (varname(n)=='howv' .or. varname(n)=='anlerr-howv') listsec0(1)=10 ! Oceanographic products listsec0(2)=2 ! GRIB Edition Number (currently 2) listsec1(1)=7 ! Id of orginating centre (Common Code Table C-1) listsec1(2)=4 !"EMC"! Id of orginating sub-centre (local table)/Table C of ON388 listsec1(3)=1 ! GRIB Master Tables Version Number (Code Table 1.0) listsec1(4)=1 !per Brent! GRIB Local Tables Version Number (Code Table 1.1) listsec1(5)=isign ! Significance of Reference Time (Code Table 1.2) if (lminmaxt) then listsec1(6)=y1_minmaxt ! Reference Time - Year (4 digits) listsec1(7)=m1_minmaxt ! Reference Time - Month listsec1(8)=d1_minmaxt ! Reference Time - Day listsec1(9)=h1_minmaxt ! Reference Time - Hour listsec1(10)=min1_minmaxt ! Reference Time - Minute listsec1(11)=s1_minmaxt ! Reference Time - Second else listsec1(6)=y0 ! Reference Time - Year (4 digits) listsec1(7)=m0 ! Reference Time - Month listsec1(8)=d0 ! Reference Time - Day listsec1(9)=h0 ! Reference Time - Hour listsec1(10)=min0 ! Reference Time - Minute listsec1(11)=s0 ! Reference Time - Second endif listsec1(12)=0 ! Production status of data (Code Table 1.3) listsec1(13)=ityped ! Type of processed data (Code Table 1.4) call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr) if (mype==0) print*,'gribcreate status=',ierr !==> Pack up Grid Definition Section (Section 3) add to GRIB2 message. igds(1)=0 !Source of grid definition (see Code Table 3.0) igds(2)=nx*ny !Number of grid points in the defined grid. igds(3)=0 !Number of octets needed for each additional grid points definition igds(4)=0 !Interpretation of list for optional points definition (Code Table 3.11) if (trim(clocation)=="conus") igds(5)=30 !Grid Definition Template Number (Code Table 3.1) if (trim(clocation)=="alaska") igds(5)=20 if (trim(clocation)=="hawaii") igds(5)=10 if (trim(clocation)=="guam" ) igds(5)=10 if (trim(clocation)=="prico") igds(5)=10 if (trim(clocation)=="cohres") igds(5)=30 if (trim(clocation)=="akhres") igds(5)=20 if (trim(clocation)=="hrrr") igds(5)=30 if (trim(clocation)=="juneau") igds(5)=20 if (trim(clocation)=="cohresext") igds(5)=30 if (trim(clocation)=="conusext") igds(5)=30 if (trim(clocation)=="cohreswexp") igds(5)=30 if (trim(clocation)=="nwrfc") igds(5)=30 if (trim(clocation)=="conus") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="alaska") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="hawaii") call apply_template_310_ndfd(igdstmpl,igdstmplen,"h") if (trim(clocation)=="guam") call apply_template_310_ndfd(igdstmpl,igdstmplen,"g") if (trim(clocation)=="prico") call apply_template_310_ndfd(igdstmpl,igdstmplen,"p") if (trim(clocation)=="cohres") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="akhres") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="hrrr") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="juneau") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="cohresext") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="conusext") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="cohreswexp") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) if (trim(clocation)=="nwrfc") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) ideflist=0 !Used if igds(3) .ne. 0. Dummy array otherwise call addgrid(cgrib,max_bytes,igds,igdstmpl,igdstmplen,ideflist,idefnum,ierr) if (mype==0) print*,'addgrid status=',ierr !==> pack up sections 4 through 7 for a given field and add them to a GRIB2 message. ! They are Product Definition Section, Data Representation Section, Bit-Map Section ! and Data Section, respectively. if (lminmaxt) then ipdsnum=8 else ipdsnum=0 !Product Definition Template Number ( see Code Table 4.0) endif if (mype==0) print*,'write_ndfd_grib2:, n,varname,typeproc,lminmaxt=', & n,varname(n),typeproc(n),lminmaxt if (lminmaxt) then if (mype.eq.0) print*,'start use of template_4p8' ipdstmplen=29 allocate (ipdstmpl(ipdstmplen)) call apply_template_4p8_ndfd(ipdstmpl,ipdstmplen,ifcsthr_minmaxt,& y2_minmaxt,m2_minmaxt,d2_minmaxt, & h2_minmaxt,min2_minmaxt,s2_minmaxt, & varname(n),typeproc(n),cgprocess) if (mype.eq.0) print*,'end use of template_4p8' else ipdstmplen=15 allocate (ipdstmpl(ipdstmplen)) call apply_template_40_ndfd(ipdstmpl,ipdstmplen,ifcsthr,& varname(n),typeproc(n),cgprocess) endif numcoord=0 coordlist=0. !needed for hybrid vertical coordinate idrsnum=0 !Data Representation Template Number ( see Code Table 5.0 ) call apply_template_50_ndfd(idrstmpl,idrstmplen,varname(n)) dscal=idrstmpl(3) !NO NEED I BELIEVE /MPondeca, May 15, 2007 field=0. do j=j1,j2 do i=i1,i2 field(i,j)=g(i,j,n) enddo enddo call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror) ij=0 do j=1,ny do i=1,nx ij=ij+1 fld(ij)=field2(i,j) enddo enddo ibmap=255 ! Bitmap indicator ( see Code Table 6.0 ) bmap(:,:)=.false. if (varname(n)=='howv' .or. varname(n)=='anlerr-howv') then ibmap=0 ij=0 do j=1,ny do i=1,nx ij=ij+1 bmap(i,j)=.true. if (fld(ij) .gt. 1.e20) bmap(i,j)=.false. enddo enddo endif ! if (varname(n).eq.'cldch') then ! ibmap=0 ! ! ij=0 ! do j=1,ny ! do i=1,nx ! ij=ij+1 ! bmap(i,j)=.true. !! if (abs(fld(ij))>0.1*cldch_bitmapval) bmap(i,j)=.false. ! if (abs(fld(ij))>=20000.) bmap(i,j)=.false. ! enddo ! enddo ! endif call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmplen, & coordlist,numcoord,idrsnum,idrstmpl, & idrstmplen,fld,nx*ny,ibmap,bmap,ierr) if (mype==0) print*,'addfield status=',ierr deallocate(ipdstmpl) ! enddo !CORRECT WAY. UNCOMMENT! !==> finalize GRIB message after all grids ! and fields have been added. It adds the End Section ( "7777" ) call gribend(cgrib,max_bytes,lengrib,ierr) if (mype==0) print*,'gribend status=',ierr if (mype==0) print*,'length of the final GRIB2 message in octets =',lengrib if (mype==0) call wryte(lunout, lengrib, cgrib) enddo !INCORRET WAY. DELETE THIS LINE if (mype==0) print*,'---- exiting write_ndfd_grib2 -------' deallocate(igdstmpl) deallocate(fld) deallocate(field) deallocate(field2) call mpi_barrier(mpi_comm_world,ierror) return end !=========================================================================== !=========================================================================== subroutine apply_template_310_ndfd(ifield3,len3,c1) implicit none integer(4) nx,ny,lat1,lon1,lat2,lon2,lad,ds1 integer(4) len3 integer(4) ifield3(len3) character(1) c1 if (c1=="h") then !hawaii nx=321 ny=225 ! lat1=18066780 !lat of 1st grd pt in micro-deg ! lon1=198374755 !east-long of 1st grd pt in micro-deg ! lat2=23082000 !lat of last grd pt in micro-deg ! lon2=206031000 !east-long of last grd pt in micro-deg lat1=18072699 !lat of 1st grd pt in micro-deg lon1=198474999 !east-long of 1st grd pt in micro-deg lat2=23087799 !lat of last grd pt in micro-deg lon2=206130999 !east-long of last grd pt in micro-deg lad=20000000 !lat at which projection intersects earth ds1=2500000 !grid spacing in x and y elseif (c1=="g") then !guam nx=193 ny=193 lat1=12349884 lon1=143686538 lat2=16794399 lon2=148280000 lad=20000000 ds1=2500000 elseif (c1=="p") then !puerto rico ! nx=177 ! ny=129 ! lat1=16828685 ! lon1=291804687 ! lat2=19747399 ! lon2=296027600 ! lad=20000000 ! ds1=2500000 nx=353 ny=257 lat1=16828700 lon1=291804700 lat2=19736200 lon2=296015500 lad=20000000 ds1=1250000 endif ifield3(1) = 1 !Earth assumed spherical with radius specified by data producer ifield3(2) = 0 ifield3(3) = 6371200 ifield3(4) = 0 ifield3(5) = 0 ifield3(6) = 0 ifield3(7) = 0 ifield3(8) = nx ifield3(9) = ny ifield3(10) = lat1 ifield3(11) = lon1 ifield3(12) = 56! 8! 56 !CHECK/MPondeca !file from NDFD site using "0" ifield3(13) = lad ifield3(14) = lat2 ifield3(15) = lon2 ifield3(16) = 64! Why is NDFD using 80? ifield3(17) = 0 ifield3(18) = ds1 ifield3(19) = ds1 return end !=========================================================================== !=========================================================================== subroutine apply_template_320_ndfd(ifield3,len3,cgrid) implicit none integer(4) len3 integer(4) ifield3(len3) character*60 cgrid integer(4) nx,ny !number of grid points in x and y-direction integer(4) ds1 !grid spacing in x and y integer(4) lat1 !lat of 1st grd pt in micro-deg integer(4) lon1 !east-long of 1st grd pt in micro-deg integer(4) latan1 !true lat in micro-deg integer(4) lonv !y-axis || to long-circle at this long if (trim(cgrid)=='alaska') then nx=825 ny=553 lat1=40530101 lon1=181429000 latan1=60000000 lonv=210000000 ds1=5953125 elseif (trim(cgrid)=='akhres') then nx=1649 ny=1105 lat1=40530101 lon1=181429000 latan1=60000000 lonv=210000000 ds1=2976563 elseif (trim(cgrid)=='juneau') then nx=655 ny=855 lat1=51500000 lon1=217500000 latan1=60000000 lonv=225000000 ds1=1448281 endif ifield3(1) = 1 !Earth assumed spherical with radius specified by data producer ifield3(2) = 0 ifield3(3) = 6371200 ifield3(4) = 0 ifield3(5) = 0 ifield3(6) = 0 ifield3(7) = 0 ifield3(8) = nx ifield3(9) = ny ifield3(10) = lat1 ifield3(11) = lon1 ifield3(12) = 8! 56 !CHECK/MPondeca !NDFD file from Web uses 0 ifield3(13) = latan1 ifield3(14) = lonv ifield3(15) = ds1 ifield3(16) = ds1 ifield3(17) = 0 ifield3(18) = 64! Why is NDFD using 80? return end !=========================================================================== !=========================================================================== !=========================================================================== subroutine apply_template_330_ndfd(ifield3,len3,cgrid) implicit none integer(4) len3 integer(4) ifield3(len3) character*60 cgrid integer(4) nx,ny !number of grid points in x and y-direction integer(4) ds1 !grid spacing in mm integer(4) lat1 !lat of 1st grd pt in micro-deg integer(4) lon1 !east-long of 1st grd pt in micro-deg integer(4) latan1 !true lat in micro-deg integer(4) lonv !y-axis || to long-circle at this long if (trim(cgrid)=='conus') then nx=1073 ny=689 ds1=5079406 lat1=20191999 lon1=238445999 latan1=25000000 lonv=265000000 else if (trim(cgrid)=='cohres') then nx=2145 ny=1377 ds1=2539703 lat1=20191999 lon1=238445999 latan1=25000000 lonv=265000000 else if (trim(cgrid)=='hrrr') then ! hrrr nx=1799 ny=1059 lat1=21138000 lon1=237280000 latan1=38500000 lonv=262500000 ds1=3000000 else if (trim(cgrid)=='cohresext') then nx=2145 ny=1597 ds1=2539703 lat1=20191999 lon1=238445999 latan1=25000000 lonv=265000000 else if (trim(cgrid)=='conusext') then nx=1073 ny=799 ds1=5079406 lat1=20191999 lon1=238445999 latan1=25000000 lonv=265000000 else if (trim(cgrid)=='cohreswexp') then nx=2345 ny=1597 ds1=2539703 lat1=19228976 lon1=233723448 latan1=25000000 lonv=265000000 else if (trim(cgrid)=='nwrfc') then nx=709 ny=795 ds1=2539703 lat1=37979684 lon1=234042704 latan1=25000000 lonv=265000000 end if ifield3(1) = 1 !Earth assumed spherical with radius specified by data producer ifield3(2) = 0 ifield3(3) = 6371200 ifield3(4) = 0 ifield3(5) = 0 ifield3(6) = 0 ifield3(7) = 0 ifield3(8) = nx ifield3(9) = ny ifield3(10) = lat1 ifield3(11) = lon1 ifield3(12) = 8! 56 ifield3(13) = latan1 ifield3(14) = lonv ifield3(15) = ds1 ifield3(16) = ds1 ifield3(17) = 0 ifield3(18) = 64! Why is NDFD using 80? ifield3(19) = latan1 ifield3(20) = latan1 ifield3(21) = -90000000 ifield3(22) = 0 return end !=========================================================================== subroutine apply_template_40_ndfd(ifield4,len4,ifcsthr,var,typep,cgprocess) implicit none integer(4) len4,ifcsthr integer(4) ifield4(len4) character*60 var character*60 typep character*60 cgprocess !==> ifield4(1):parameter category (see Code Table 4.1) !==> ifield4(2):parameter number (see Code Table 4.2) if (trim(var) .eq. 'psfc' .or. & trim(var) .eq. 'anlerr-psfc') then ifield4(1) = 3 ifield4(2) = 0 elseif (trim(var) .eq. 't2m' .or. & trim(var) .eq. 'anlerr-t2m' ) then ifield4(1) = 0 ifield4(2) = 0 elseif (trim(var) .eq. 'td2m' .or. & trim(var) .eq. 'anlerr-td2m' ) then ifield4(1) = 0 ifield4(2) = 6 elseif (trim(var) .eq. 'u10m' .or. & trim(var) .eq. 'anlerr-u10m') then ifield4(1) = 2 ifield4(2) = 2 elseif (trim(var) .eq. 'v10m' .or. & trim(var) .eq. 'anlerr-v10m') then ifield4(1) = 2 ifield4(2) = 3 elseif (trim(var) .eq. 'wdir10m' .or. & trim(var) .eq. 'anlerr-wdir10m') then ifield4(1) = 2 ifield4(2) = 0 elseif (trim(var) .eq. 'wspd10m' .or. & trim(var) .eq. 'anlerr-wspd10m') then ifield4(1) = 2 ifield4(2) = 1 elseif (trim(var) .eq. 'q2m' .or. & trim(var) .eq. 'anlerr-q2m') then ifield4(1) = 1 ifield4(2) = 0 elseif (trim(var) .eq. 'hgt')then ifield4(1) = 3 ifield4(2) = 5 elseif (trim(var) .eq. 'gust' .or. & trim(var) .eq. 'anlerr-gust')then ifield4(1) = 2 ifield4(2) = 22 elseif (trim(var) .eq. 'vis' .or. & trim(var) .eq. 'anlerr-vis')then ifield4(1) = 19 ifield4(2) = 0 elseif (trim(var) .eq. 'pblh' .or. & trim(var) .eq. 'anlerr-pblh')then ifield4(1) = 19 ifield4(2) = 12 elseif (trim(var) .eq. 'ctrwspd10m' .or. & trim(var) .eq. 'anlerr-ctrwspd10m')then ifield4(1) = 2 ifield4(2) = 1 elseif (trim(var) .eq. 'usfc' .or. & trim(var) .eq. 'anlerr-usfc') then ifield4(1) = 2 ifield4(2) = 2 elseif (trim(var) .eq. 'vsfc' .or. & trim(var) .eq. 'anlerr-vsfc') then ifield4(1) = 2 ifield4(2) = 3 elseif (trim(var) .eq. 'wdirsfc' .or. & trim(var) .eq. 'anlerr-wdirsfc') then ifield4(1) = 2 ifield4(2) = 0 elseif (trim(var) .eq. 'wspdsfc' .or. & trim(var) .eq. 'anlerr-wspdsfc') then ifield4(1) = 2 ifield4(2) = 1 elseif (trim(var) .eq. 'u20m' .or. & trim(var) .eq. 'anlerr-u20m') then ifield4(1) = 2 ifield4(2) = 2 elseif (trim(var) .eq. 'v20m' .or. & trim(var) .eq. 'anlerr-v20m') then ifield4(1) = 2 ifield4(2) = 3 elseif (trim(var) .eq. 'ctrtd2m' .or. & trim(var) .eq. 'anlerr-ctrtd2m' ) then ifield4(1) = 0 ifield4(2) = 6 elseif (trim(var) .eq. 'pmsl' .or. & trim(var) .eq. 'anlerr-pmsl') then ifield4(1) = 3 ifield4(2) = 1 elseif (trim(var) .eq. 'howv' .or. & trim(var) .eq. 'anlerr-howv') then ifield4(1) = 0 ifield4(2) = 3 elseif (trim(var) .eq. 'tcamt' .or. & trim(var) .eq. 'anlerr-tcamt') then ifield4(1) = 6 ifield4(2) = 1 elseif (trim(var) .eq. 'lcbas' .or. & trim(var) .eq. 'anlerr-lcbas') then ifield4(1) = 6 ifield4(2) = 11 elseif (trim(var) .eq. 'cldch' .or. & trim(var) .eq. 'anlerr-cldch') then ifield4(1) = 6 !3 ifield4(2) = 13 !5 else print*,'define param category and param number for field ',var print*,'sorry, ... aborting in apply_template_40_ndfd' call abort endif !==> ifield4(3):type of generating process (see Code Table 4.3) if (trim(typep) .eq. 'analysis') then ifield4(3) = 0 elseif (trim(typep) .eq. 'forecast') then ifield4(3) = 2 elseif (trim(typep) .eq. 'analysis error') then ifield4(3) = 7 else print*,'for RTMA,typep must be either analysis, forecast, & & or analysis error' print*,'sorry, ... aborting in apply_template_40_ndfd' call abort endif !==> ifield4(4):background generating process identifier ! (defined by originating Center) ifield4(4) = 0 !hasn't been defined yet !==> ifield4(5):analysis or forecast generating process identifier ! (defined by originating Center) ifield4(5) = 109 !default is rtma if (trim(cgprocess)=='urma' .or. trim(cgprocess)=='URMA') ifield4(5) = 118 !==> ifield4(6):hours of observational data cutoff after reference time !==> ifield4(7):minutes of observational data cutoff after reference time ifield4(6) = 0 !per steve ifield4(7) = 0 !==> ifield4(8):indicator of unit of time range (see Code Table 4.4) ifield4(8) = 1 if (ifcsthr < 0) ifield4(8) = 0 !==> ifield4(9):forecast time in units defined by ifield4(8) ifield4(9) = 0 if (trim(typep) .eq. 'forecast') ifield4(9) = abs(ifcsthr) !==> ifield4(10):type of first fixed surface (see Code Table 4.5) !==> ifield4(11):scale factor of first fixed surface !==> ifield4(12):scaled value of first fixed surface ifield4(11) = 0 !Because not saving any precision if (trim(var) .eq. 'psfc' .or. trim(var) .eq. 'anlerr-psfc' .or. & trim(var) .eq. 'hgt' .or. & trim(var) .eq. 'vis' .or. trim(var) .eq. 'anlerr-vis' ) then ifield4(10) = 1 ifield4(12) = 0 else ifield4(10) = 103 if (trim(var) .eq. 't2m' .or. & trim(var) .eq. 'anlerr-t2m') ifield4(12) = 2 if (trim(var) .eq. 'td2m' .or. & trim(var) .eq. 'anlerr-td2m') ifield4(12) = 2 if (trim(var) .eq. 'u10m' .or. & trim(var) .eq. 'anlerr-u10m') ifield4(12) = 10 if (trim(var) .eq. 'v10m' .or. & trim(var) .eq. 'anlerr-v10m') ifield4(12) = 10 if (trim(var) .eq. 'wdir10m' .or. & trim(var) .eq. 'anlerr-wdir10m') ifield4(12) = 10 if (trim(var) .eq. 'wspd10m' .or. & trim(var) .eq. 'anlerr-wspd10m') ifield4(12) = 10 if (trim(var) .eq. 'q2m' .or. & trim(var) .eq. 'anlerr-q2m') ifield4(12) = 2 if (trim(var) .eq. 'gust' .or. & trim(var) .eq. 'anlerr-gust') ifield4(12) = 10 if (trim(var) .eq. 'pblh' .or. & trim(var) .eq. 'anlerr-pblh') ifield4(12) = 220 endif !add new vars / MPondeca, 29Jul2014 if (trim(var) .eq. 'ctrwspd10m' .or. trim(var) .eq. 'anlerr-ctrwspd10m') then ifield4(10) = 1 !assume ground level ifield4(12) = 0 endif if (trim(var) .eq. 'u20m' .or. trim(var) .eq. 'anlerr-u20m' .or. & trim(var) .eq. 'v20m' .or. trim(var) .eq. 'anlerr-v20m') then ifield4(10) = 103 ifield4(12) = 20 endif if (trim(var) .eq. 'ctrtd2m' .or. trim(var) .eq. 'anlerr-ctrtd2m') then ifield4(10) = 1 !assume ground level ifield4(12) = 0 endif if (trim(var) .eq. 'pmsl' .or. trim(var) .eq. 'anlerr-pmsl') then ifield4(10) = 101 ifield4(12) = 0 !is this value correct? / MPondeca, 22Jul2014 endif if (trim(var) .eq. 'howv' .or. trim(var) .eq. 'anlerr-howv') then ifield4(10) = 1 ifield4(12) = 1 endif if (trim(var) .eq. 'tcamt' .or. trim(var) .eq. 'anlerr-tcamt') then ifield4(10) = 200 !is this value correct? / MPondeca, 22Jul2014 ifield4(12) = 0 endif if (trim(var) .eq. 'lcbas' .or. trim(var) .eq. 'anlerr-lcbas') then ifield4(10) = 2 ifield4(12) = 0 endif if (trim(var) .eq. 'cldch' .or. trim(var) .eq. 'anlerr-cldch') then ifield4(10) = 215 ifield4(12) = 0 !find out what to write here /MPondeca, 22Jul2014 endif if (trim(var) .eq. 'usfc' .or. trim(var) .eq. 'anlerr-usfc') then ifield4(10) = 1 ifield4(12) = 0 endif if (trim(var) .eq. 'vsfc' .or. trim(var) .eq. 'anlerr-vsfc') then ifield4(10) = 1 ifield4(12) = 0 endif if (trim(var) .eq. 'wdirsfc' .or. trim(var) .eq. 'anlerr-wdirsfc') then ifield4(10) = 1 ifield4(12) = 0 endif if (trim(var) .eq. 'wspdsfc' .or. trim(var) .eq. 'anlerr-wspdsfc') then ifield4(10) = 1 ifield4(12) = 0 endif !==> ifield4(13):type of second fixed surface(See Code Table 4.5) !==> ifield4(14):scale factor of second fixed surface !==> ifield4(15):scaled value of second fixed surface ifield4(13) = 255 ifield4(14) = 0 ifield4(15) = 0 return end !=========================================================================== !=========================================================================== subroutine apply_template_4p8_ndfd(ifield4,len4,ifcsthr, & y2,m2,d2,h2,min2,s2,var,typep,cgprocess) implicit none integer(4) len4,ifcsthr integer(4) ifield4(len4) integer(4) y2,m2,d2,h2,min2,s2 character*60 var character*60 typep character*60 cgprocess !==> ifield4(1):parameter category (see Code Table 4.1) !==> ifield4(2):parameter number (see Code Table 4.2) if (trim(var) .eq. 'mxtm' .or. & trim(var) .eq. 'anlerr-mxtm' ) then ifield4(1) = 0 ifield4(2) = 4 elseif (trim(var) .eq. 'mitm' .or. & trim(var) .eq. 'anlerr-mitm' ) then ifield4(1) = 0 ifield4(2) = 5 else print*,'define param category and param number for field ',var print*,'sorry, ... aborting in apply_template_40_ndfd' call abort endif !==> ifield4(3):type of generating process (see Code Table 4.3) if (trim(typep) .eq. 'analysis') then ifield4(3) = 0 elseif (trim(typep) .eq. 'forecast') then ifield4(3) = 2 elseif (trim(typep) .eq. 'analysis error') then ifield4(3) = 7 else print*,'for RTMA,typep must be either analysis, forecast, & & or analysis error' print*,'sorry, ... aborting in apply_template_40_ndfd' call abort endif !==> ifield4(4):background generating process identifier ! (defined by originating Center) ifield4(4) = 0 !hasn't been defined yet !==> ifield4(5):analysis or forecast generating process identifier ! (defined by originating Center) ifield4(5) = 109 !default is rtma if (trim(cgprocess)=='urma' .or. trim(cgprocess)=='URMA') ifield4(5) = 118 !==> ifield4(6):hours of observational data cutoff after reference time !==> ifield4(7):minutes of observational data cutoff after reference time ifield4(6) = 0 !per steve ifield4(7) = 0 !==> ifield4(8):indicator of unit of time range (see Code Table 4.4) ifield4(8) = 1 if (ifcsthr < 0) ifield4(8) = 0 !==> ifield4(9):forecast time in units defined by ifield4(8) ifield4(9) = 0 if (trim(typep) .eq. 'forecast') ifield4(9) = abs(ifcsthr) !==> ifield4(10):type of first fixed surface (see Code Table 4.5) !==> ifield4(11):scale factor of first fixed surface !==> ifield4(12):scaled value of first fixed surface ifield4(11) = 0 !Because not saving any precision if (trim(var) .eq. 'mxtm' .or. trim(var) .eq. 'anlerr-mxtm' .or. & trim(var) .eq. 'mitm' .or. trim(var) .eq. 'anlerr-mitm' ) then ifield4(10) = 103 ifield4(12) = 2 endif !==> ifield4(13):type of second fixed surface(See Code Table 4.5) !==> ifield4(14):scale factor of second fixed surface !==> ifield4(15):scaled value of second fixed surface ifield4(13) = 255 ifield4(14) = 0 ifield4(15) = 0 ifield4(16) = y2 ifield4(17) = m2 ifield4(18) = d2 ifield4(19) = h2 ifield4(20) = min2 ifield4(21) = s2 ifield4(22) = 1 !Nmber of range specifications describing the !time intervals used to calculate statistically processed field ifield4(23) = 0 ! Total number of data values missing in statistical process if (trim(var) .eq. 'mxtm' .or. & trim(var) .eq. 'anlerr-mxtm' ) then ifield4(24) = 2 !statistical process used to calculate the processed field endif if (trim(var) .eq. 'mitm' .or. & trim(var) .eq. 'anlerr-mitm' ) then ifield4(24) = 3 endif ifield4(25) = 255 ifield4(26) = 1 !255 ifield4(27) = 12 !255 ifield4(28) = 255 ifield4(29) = 255 return end !=========================================================================== !=========================================================================== subroutine apply_template_50_ndfd(ifield5,len5,var) !******************************************************************************** ! ifield5(1): reference value(R) (IEEE 32-bit floating-point value) * ! ifield5(2): binary scale factor (E) * ! ifield5(3): decimal scale factor (D) * ! ifield5(4): number of bits used for each packed value for simple packing * ! or for each group reference value for complex packing or * ! spatial differencing * ! ifield5(5): type of original field values (See Code Table 5.1) * !******************************************************************************** implicit none integer(4) len5 integer(4) ifield5(len5) character*60 var ifield5(1)=0 !Any value. Will be overwritten ifield5(2)=0 if (trim(var) .eq. 'psfc') then ifield5(3) = 0! -1 elseif (trim(var) .eq. 't2m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'td2m') then ifield5(3) = 2! 2 elseif (trim(var) .eq. 'u10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'v10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'wdir10m') then ifield5(3) = 1! 0 elseif (trim(var) .eq. 'wspd10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'q2m') then ifield5(3) = 6 elseif (trim(var) .eq. 'gust') then ifield5(3) = 2 elseif (trim(var) .eq. 'vis') then ifield5(3) = 0 elseif (trim(var) .eq. 'pblh') then ifield5(3) = 0 elseif (trim(var) .eq. 'ctrwspd10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'usfc') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'vsfc') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'wdirsfc') then ifield5(3) = 1! 0 elseif (trim(var) .eq. 'wspdsfc') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'u20m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'v20m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'ctrtd2m') then ifield5(3) = 2! 2 elseif (trim(var) .eq. 'mxtm') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'mitm') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'pmsl') then ifield5(3) = 0! -1 elseif (trim(var) .eq. 'howv') then ifield5(3) = 2 elseif (trim(var) .eq. 'tcamt') then ifield5(3) = 2 elseif (trim(var) .eq. 'lcbas') then ifield5(3) = 1 elseif (trim(var) .eq. 'cldch') then ifield5(3) = 1 elseif (trim(var) .eq. 'anlerr-psfc') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-t2m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-td2m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-u10m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-v10m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-q2m') then ifield5(3) = 6 elseif (trim(var) .eq. 'anlerr-wdir10m') then ifield5(3) = 1! 0 elseif (trim(var) .eq. 'anlerr-wspd10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'hgt') then ifield5(3) = 1 elseif (trim(var) .eq. 'anlerr-gust') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-vis') then ifield5(3) = 0 elseif (trim(var) .eq. 'anlerr-pblh') then ifield5(3) = 0 elseif (trim(var) .eq. 'anlerr-ctrwspd10m') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'anlerr-usfc') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-vsfc') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-wdirsfc') then ifield5(3) = 1! 0 elseif (trim(var) .eq. 'anlerr-wspdsfc') then ifield5(3) = 2! 1 elseif (trim(var) .eq. 'anlerr-u20m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-v20m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-ctrtd2m') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-mxtm') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-mitm') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-pmsl') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-howv') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-tcamt') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-lcbas') then ifield5(3) = 2 elseif (trim(var) .eq. 'anlerr-cldch') then ifield5(3) = 2 else print*,'please define decimal scale factor ',var print*,'sorry, ... aborting in apply_template_50_ndfd' call abort endif ifield5(4) = 0 !Must reset to 0 ifield5(5) = 0 return end !=========================================================================== !=========================================================================== !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: GETLUN GET UNIQUE LOGICAL UNIT NUMBERS ! PRGMMR: SMITH, TRACY ORG: FSL/PROFS DATE: 90-06-15 ! ! ABSTRACT: THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE ! OR RETURNS THEM TO THE POOL FOR CLFILE. ! ! PROGRAM HISTORY LOG: ! FORTRAN 90 VERSION IS GETLUN90: PONDECA, DATE: 2006-03-08 ! ! USAGE: CALL GETLUN(LUN,OPTN) ! INPUT ARGUMENT LIST: ! LUN - INTEGER LOGICAL UNIT NUMBER ! OPTN - INTEGER CNCT=1, DSCT=2. ! IF CONNECTING A FILE(CNCT) SET THE NUMBER TO ! NEGATIVE SO IT WON'T BE USED UNTIL AFTER ! DSCT SETS IT POSITIVE. ! ! OUTPUT ARGUMENT LIST: ! LUN - INTEGER LOGICAL UNIT NUMBER ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN-90 ! MACHINE: NAS-9000 !$$$ SUBROUTINE GETLUN90(LUN,OPTN) !* THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE !* OR RETURNS THEM TO THE POOL FOR CLFILE IMPLICIT NONE INTEGER LUN,NUM(80),OPTN,CNCT,DSCT,I PARAMETER (CNCT=1,DSCT=2) SAVE NUM ! DATA NUM/99,98,97,96,95,94,93,92,91,90, & 89,88,87,86,85,84,83,82,81,80, & 79,78,77,76,75,74,73,72,71,70, & 69,68,67,66,65,64,63,62,61,60, & 59,58,57,56,55,54,53,52,51,50, & 49,48,47,46,45,44,43,42,41,40, & 39,38,37,36,35,34,33,32,31,30, & 29,28,27,26,25,24,23,22,21,20/ !* START IF(OPTN.EQ.CNCT) THEN DO 10 I=1,80 IF(NUM(I).GT.0) THEN LUN=NUM(I) NUM(I)=-NUM(I) GOTO 20 ENDIF 10 CONTINUE PRINT*, 'NEED MORE THAN 80 UNIT NUMBERS' 20 CONTINUE ELSE IF(OPTN.EQ.DSCT) THEN !* MAKE THE NUMBER AVAILABLE BY SETTING POSITIVE DO 30 I=1,80 IF(LUN.EQ.-NUM(I)) NUM(I)=ABS(NUM(I)) 30 CONTINUE END IF RETURN END !=========================================================================== !======================================================================= !======================================================================= subroutine smther_one_same(g1,is,ie,js,je,ns) !$$$ subprogram documentation block ! . . . . ! subprogram: smther_one ! prgmmr: pondeca org: np22 date: 2006-08-01 ! ! abstract: apply 1-2-1 smoother in each direction of data slab ! ! program history log: ! 2006-08-01 pondeca ! ! input argument list: ! g1 - 2d array of field to be smoothed ! is,ie - first and last i values of data slab ! js,je - first and last j values of data slab ! ns - number of passes ! ! output argument list: ! g1 - smoothed 2d field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! implicit none integer(4) is, ie, js, je integer(4) i,j,l,ip,im,jp,jm integer(4), intent(in) :: ns real(4), dimension(is:ie, js:je), intent(inout) :: g1 ! on input: original data slab ! on ouput: filtered data slab real(4), allocatable:: g2(:,:) allocate(g2(is:ie,js:je)) do l=1,ns do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=.25*(g1(ip,j)+g1(im,j))+.5*g1(i,j) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=.25*(g2(i,jp)+g2(i,jm))+.5*g2(i,j) end do end do end do deallocate(g2) return end subroutine smther_one_same !======================================================================= !======================================================================= subroutine smoother_desmoother(g1,is,ie,js,je,smcf1,smcf2,one2one,ns) !$$$ subprogram documentation block ! . . . . ! subprogram: smoother_desmoother ! prgmmr: pondeca org: np22 date: 2006-08-01 ! ! abstract: apply smoother-desmoother to 2-d data slab ! ! program history log: ! 2011-09-15 pondeca ! ! input argument list: ! g1 - 2d array of field to be smoothed ! is,ie - first and last i values of data slab ! js,je - first and last j values of data slab ! ns - number of passes ! ! output argument list: ! g1 - smoothed 2d field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block implicit none integer, parameter :: r_single = selected_real_kind(6) ! single precision integer, parameter :: i_kind = selected_int_kind(8) ! long integer integer(i_kind) , intent(in ) :: is, ie, js, je integer(i_kind) , intent(in ) :: ns real(r_single), dimension(is:ie, js:je), intent(inout) :: g1 ! on input: original data slab ! on ouput: filtered data slab real(r_single), intent(in) :: smcf1 ! smoothing coefficient real(r_single), intent(in) :: smcf2 ! desmoothing coefficient logical,intent(in) :: one2one ! if .true., perform 1-2-1 smoothing instead integer(i_kind) i,j,l,ip,im,jp,jm,istep real(r_single), allocatable:: g2(:,:) real(r_single) smcf0 real(r_single),parameter::half=0.5_r_single real(r_single),parameter::quarter=0.25_r_single allocate(g2(is:ie,js:je)) do l=1,ns if (.not.one2one) then do istep=1,2 if (istep==1) smcf0=+abs(smcf1) if (istep==2) smcf0=-abs(smcf2) do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=g1(i,j)+smcf0*(0.5*(g1(ip,j)+g1(im,j))-g1(i,j)) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=g2(i,j)+smcf0*(0.5*(g2(i,jp)+g2(i,jm))-g2(i,j)) end do end do end do elseif (one2one) then do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=quarter*(g1(ip,j)+g1(im,j))+half*g1(i,j) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=quarter*(g2(i,jp)+g2(i,jm))+half*g2(i,j) end do end do endif end do deallocate(g2) return end subroutine smoother_desmoother !======================================================================= !=======================================================================