PROGRAM GET_WEIGHTS !---------------------------------------------------------- ! CHANGE OF VARIABLE: ANALYSIS GRID => MODEL GRID ! GET HORIZONTAL INTERPOLATION WEIGHTS ! (interpolation from regular grid to irregular grid) ! ! The analysis grid is assumed larger than the model grid ! to facilitate the interpolation ! ! ETA MODEL : Arakawa E-grid ! ! OUTPUT: ! ! WGTU, IWGTU = u-point intpl weights ! WGTV, IWGTV = v-point intpl weights ! WGTT, IWGTT = T-point intpl weights ! WGTM, IWGTM = M-point intpl weights !------------------------------------------------------------ INCLUDE 'ntype.h' !! NTYPE= No. of cntrl vrbl types INTEGER iord,lhalf,lbig parameter(iord=1) ! order of intpl (1=lin, 2=quadratic) parameter(lhalf=1) ! Jim's improvement (lhalf=0) real,dimension(:,:,:),allocatable :: WGTU,WGTV,WGTT,WGTM integer,dimension(:,:,:),allocatable :: IWGTU,IWGTV integer,dimension(:,:,:),allocatable :: IWGTT,IWGTM real,dimension(:,:),allocatable :: x1in0,x2in0 !! output grid coord real,dimension(:,:),allocatable :: x1in0v real,dimension(:),allocatable :: x1grid0,x2grid0 !! input grid coord real,dimension(:),allocatable :: dxj real,dimension(:,:),allocatable :: PHI real dxeta,dyeta integer,dimension(:,:),allocatable :: iflagu,iflagv,iflagt,iflagm CHARACTER*20 FILENAME,FILEPHI REAL xstart,xstartv INTEGER nin,imeta,jmeta,n1grid,n2grid INCLUDE 'namelist.h' !!! define namelists (all) !------------------------------------------------------------- !--- read input files 'cntrl_index.name' (33) and 'cntrl_specs.name' (34) write(*,*) "START INTPL WEIGHTS " INCLUDE 'read_namelist' ! read all namelists ! in this version, model and analysis grids have same dimensions ! (it will be relaxed in later version) ! imeta,jmeta = model ! n1grid,n2grid = analysis nin=IM*JM imeta=IM jmeta=JM n1grid=IM n2grid=JM ! define lbig (number of interpolating points) if(lhalf.eq.0) then lbig = (iord*(iord+3)+2)/2 else lbig = (iord+1)**2 end if write(*,*) "order of interpolation is ",iord write(*,*) "number of interpolating points ",lbig allocate(DXJ(jmeta)) allocate(PHI(imeta,jmeta)) CALL INTLZ (IM,JM,dxeta0,dyeta0,DXJ,PHI) ! add dxeta to the left and right of the original grid ! add dyeta up and below the original grid xl_tot=(im-1)*dxeta0+2.*dxeta0 yl_tot=(jm-1)*dyeta0+2.*dyeta0 dxeta=xl_tot/(im-1) dyeta=yl_tot/(jm-1) dist=max(dxeta,dyeta) ! note: dxeta, dxj represent the A-grid distances allocate(x1in0(imeta,jmeta)) ; allocate(x2in0(imeta,jmeta)) allocate(x1in0v(imeta,jmeta)) allocate(x1grid0(n1grid)) ; allocate(x2grid0(n2grid)) ! define input (regular,analysis) grid coordinates ! NOTE: since A-grid, same for all variables do n=1,n1grid x1grid0(n)=(n-1)*dxeta/dist write(*,300) n,x1grid0(n) end do 300 format('original (regular) grid: i=',i3,' xgrid=',f7.2) do n=1,n2grid x2grid0(n)=(n-1)*dyeta/dist write(*,400) n,x2grid0(n) end do 400 format('original (regular) grid: j=',i3,' ygrid=',f7.2) ! define maximum x-axis length ! define output (irregular,model) grid coordinates ! There is a dyeta shift of the origin ! NOTE: for E-grid y-axis is same for all variables do j=1,jmeta do i=1,imeta x2in0(i,j)=dyeta0/dist + (j-1)*dyeta/dist end do end do do j=1,jmeta write(*,200) j,x2in0(1,j) end do 200 format('interpolated (irregular) grid: j=',i3,' yin=',f7.2) ! T-points (mass) ! There is a dxeta shift of the origin do j=1,jmeta nmod=MOD(J,2) XL=(im-1)*dxj(j) xstart=dxeta0/dist + 0.5*(XL_tot-XL)/dist + 0.5*(1-nmod)*dxj(j)/dist x1in0(1,j)=xstart do i=2,imeta x1in0(i,j)=xstart+(i-1)*dxj(j)/dist end do end do do i=1,imeta write(*,101) i,x1in0(i,1) end do 101 format('interpolated (irregular) T grid: i=',i3,' xin=',f7.2) ! V-points (wind) do j=1,jmeta nmod=MOD(J,2) XL=(im-1)*dxj(j) xstartv=dxeta0/dist + 0.5*(XL_max-XL)/dist + 0.5*nmod*dxj(j)/dist x1in0v(1,j)=xstartv do i=2,imeta x1in0v(i,j)=xstartv+(i-1)*dxj(j)/dist end do end do do i=1,imeta write(*,102) i,x1in0v(i,1) end do 102 format('interpolated (irregular) V grid: i=',i3,' xin=',f7.2) !-------------------------------------------------- !-------------------------------------------------- allocate(IWGTU(imeta,jmeta,lbig)) ; allocate(WGTU(imeta,jmeta,lbig)) allocate(IWGTV(imeta,jmeta,lbig)) ; allocate(WGTV(imeta,jmeta,lbig)) allocate(IWGTT(imeta,jmeta,lbig)) ; allocate(WGTT(imeta,jmeta,lbig)) allocate(IWGTM(imeta,jmeta,lbig)) ; allocate(WGTM(imeta,jmeta,lbig)) allocate(iflagu(imeta,jmeta)) allocate(iflagv(imeta,jmeta)) allocate(iflagt(imeta,jmeta)) allocate(iflagm(imeta,jmeta)) nf=1 ; ndx=0 ; ndy=0 ; ndxx=0 ndxy=0 ; ndyy=0 call simpin2f(wgtt,wgtt,wgtt,wgtt,wgtt,wgtt,iwgtt, & iflagt,x1in0,x2in0,nin,iord,lhalf,lbig,x1grid0,x2grid0, & n1grid,n2grid,nf,ndx,ndy,ndxx,ndxy,ndyy) ! for eta-model, just assign m=t (not used) wgtm (:,:,:)= wgtt (:,:,:) iwgtm (:,:,:)=iwgtt (:,:,:) iflagm(:,:)=iflagt(:,:) nf=1 ; ndx=0 ; ndy=0 ; ndxx=0 ndxy=0 ; ndyy=0 call simpin2f(wgtv,wgtv,wgtv,wgtv,wgtv,wgtv,iwgtv, & iflagv,x1in0v,x2in0,nin,iord,lhalf,lbig,x1grid0,x2grid0, & n1grid,n2grid,nf,ndx,ndy,ndxx,ndxy,ndyy) write(*,*) " " ! for eta-model, u and v points are the same wgtu (:,:,:)= wgtv (:,:,:) iwgtu (:,:,:)=iwgtv (:,:,:) iflagu(:,:)=iflagv(:,:) !========================================== !---------- write weights --------------- write(FILENAME,100) 100 format('intpl_wghts') nout=55 CLOSE(nout) OPEN(UNIT=nout,FILE=FILENAME,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' nout OPEN UNIT ERROR IER=',IER write(nout) lbig write(nout) wgtu,wgtv,wgtt,wgtm write(nout) iwgtu,iwgtv,iwgtt,iwgtm CLOSE(nout,STATUS='KEEP') !--- write dx,dy,phi ------------------- write(FILEPHI,301) 301 format('dxdy_phi') nphi=60 CLOSE(nphi) OPEN(UNIT=nphi,FILE=FILEPHI,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' nphi OPEN UNIT ERROR IER=',IER write(nphi) dxeta,dyeta,((PHI(i,j),i=1,IM),j=1,JM) CLOSE(nphi,STATUS='KEEP') !-------------------------------------- deallocate(WGTU) deallocate(WGTV) deallocate(WGTT) deallocate(WGTM) deallocate(IWGTU) deallocate(IWGTV) deallocate(IWGTT) deallocate(IWGTM) deallocate(iflagt) deallocate(iflagu) deallocate(iflagv) deallocate(iflagm) deallocate(x1in0) ; deallocate(x2in0) deallocate(x1in0v) deallocate(DXJ) deallocate(PHI) deallocate(x1grid0) ; deallocate(x2grid0) STOP END