C****************************************************************************** C * C Program to interpolate ADCIRC results from one FE grid to another FE grid. * C * C Input: 2 files in ADCIRC grid format (UNIT 14) * C - Only the standard first 2 lines (alpha label, ne,np) and the node * C and element tables are required for the grid being interpolated FROM* C - The standard first 2 lines (alpha label, ne,np), the node table, * C the element table and the boundary information are required for * C the grid being interpolated ONTO. * C - The depth field in the node table in both grids may be omitted. * C * C Output: Information necessary to interpolate results to the ONTO grid: * C - the node number of the ONTO grid * C - node #1 of the FROM grid element containing the ONTO node * C - node #2 of the FROM grid element containing the ONTO node * C - node #3 of the FROM grid element containing the ONTO node * C (these values are < 0 if the node lies outside the element) * C - the value the variable to be interpolated at node #1 should be * C multiplied by to perform the interpolation (STA1) * C - the value the variable to be interpolated at node #2 should be * C multiplied by to perform the interpolation (STA2) * C - the value the variable to be interpolated at node #3 should be * C multiplied by to perform the interpolation (STA3) * C * C * C Programmed by Rick Luettich, UNC Institute of Marine Sciences * C last modification 1/28/94 * C * C * C****************************************************************************** C PARAMETER STATEMENT VARIABLES * C * C MNPF = maximum number of nodal points in grid to be interpolated from * C MNEF = maximum number of element in grid to be interpolated from * C MNPO = maximum number of nodal points in grid to be interpolated onto * C * C****************************************************************************** C * C The following variables must also be set in the code (lines 73-76) * C * C SLAM0 = center longitude of CPP projection (deg) * C -79.0 for E.C. grid * C SFEA0 = center latitude of CPP projection (deg) * C 35.0 for E.C. grid * C * C NOTE: if SLAM0 = 0 and SFEA0 = 0, it is assumed the grid is NOT in * C lat,lon coordinates and no transformation is done. * C * C****************************************************************************** Program interp C C...Parameter Statements C PARAMETER(MNPF=195000) PARAMETER(MNEF=380000) PARAMETER(MNPO=55000) C C...Dimension arrays and define common blocks C implicit double precision (a-h,o-z) DIMENSION xf(mnpf),yf(mnpf),nmf1(mnef),nmf2(mnef),nmf3(mnef) DIMENSION xo(mnpo),yo(mnpo),sta1(mnpo),sta2(mnpo),sta3(mnpo), & badprox(mnpo),areaf(mnef) INTEGER badnode(mnpo),boundnode(mnpo),nnef(mnpo) CHARACTER*60 fgridf,ogridf,intoutf,intdiagf CHARACTER*24 agridf,agrido LOGICAL found C... C...DEFINE PI, SLAM0, SFEA0, COMPUTE SLAM0 AND SFEA0 IN RADIANS C... PI=3.141592653589793d0 SLAM0=-79.0d0 SFEA0=35.0d0 c SLAM0=0.0d0 c SFEA0=0.0d0 SLAM0R=SLAM0*PI/180.0d0 SFEA0R=SFEA0*PI/180.0d0 1055 FORMAT(a55) 1110 FORMAT(' File ',A55,' WAS NOT FOUND!',/,' **** Try again *****',/) 1111 FORMAT(' File ',A55,' WAS FOUND!',/) 31 write(*,*) ' Enter name of the grid file to be interpolated FROM' write(*,*) ' ' read(*,1055) fgridf inquire(file=fgridf,exist=found) if(found) goto 32 write(*,1110) fgridf goto 31 32 write(*,1111) fgridf 41 write(*,*) ' Enter name of the grid file to be interpolated ONTO' write(*,*) ' ' read(*,1055) ogridf inquire(file=ogridf,exist=found) if(found) goto 42 write(*,1110) ogridf goto 41 42 write(*,1111) ogridf write(*,*) ' Enter name of the interpolation output file' write(*,*) ' ' read(*,1055) intoutf write(*,*) ' Enter name of the interpolation diagnostics file' write(*,*) ' ' read(*,1055) intdiagf open(11,file=intdiagf) write(11,*) ' This file contains the nodes in the ONTO file that ' write(11,*) ' were not found within any of the elements in the ' write(11,*) ' FROM grid.' write(11,*) ' ' write(11,1999) 1999 format(' ONTO node',3x,'prox index',3x,'closest FROM ele',/) C... C...Open, read and store interpolated FROM grid information from UNIT 14 C... open(14,file=fgridf) write(*,*) ' Reading information in FROM grid file...............' write(*,*) ' ' read(14,'(A24)') agridf read(14,*) NEF,NPF if(nef.gt.mnef) then write(*,1000) write(*,1001) 1001 Format(' The # elements in the grid being interpolated FROM',/, & ' is > than the parameter MNEF. Respecify MNEF in ',/, & ' the parameter statement and recompile the code.') write(*,1010) stop endif if(npf.gt.mnpf) then write(*,1000) write(*,1002) 1002 Format(' The # nodes in the grid being interpolated FROM',/, & ' is > than the parameter MNPF. Respecify MNPF in ',/, & ' the parameter statement and recompile the code.') write(*,1010) stop endif 1000 format(' ***************** Dimensioning Error ******************') 1010 format(' ***************** Program Terminated ******************') C...FROM grid node table do i=1,npf read(14,*) jki,slam,sfea if((slam0.ne.0.).and.(sfea0.ne.0.)) then slam=pi*slam/180.d0 sfea=pi*sfea/180.d0 CALL CPP(xf(jki),yf(jki),slam,sfea,slam0r,sfea0r) else xf(jki)=slam yf(jki)=sfea endif enddo C...FROM grid element table do k=1,nef read(14,*) ijk,nhy,nmf1(ijk),nmf2(ijk),nmf3(ijk) x1=xf(nmf1(ijk)) x2=xf(nmf2(ijk)) x3=xf(nmf3(ijk)) y1=yf(nmf1(ijk)) y2=yf(nmf2(ijk)) y3=yf(nmf3(ijk)) areaf(ijk)=(x1-x3)*(y2-y3) - (x2-x3)*(y1-y3) end do close(14) C... C...Open, read and store interpolated ONTO grid information from UNIT 14 C... open(14,file=ogridf) write(*,*) ' Reading information in ONTO grid file...............' write(*,*) ' ' read(14,'(A24)') agrido read(14,*) neo,npo if(npo.gt.mnpo) then write(*,1000) write(*,1003) 1003 Format(' The # nodes in the grid being interpolated ONTO',/, & ' is > than the parameter MNPO. Respecify MNPO in ',/, & ' the parameter statement and recompile the code.') write(*,1010) endif C...ONTO grid node table do i=1,npo read(14,*) jki,slam,sfea if((slam0.ne.0.).and.(sfea0.ne.0.)) then slam=pi*slam/180.d0 sfea=pi*sfea/180.d0 CALL CPP(xo(jki),yo(jki),slam,sfea,slam0r,sfea0r) else xo(jki)=slam yo(jki)=sfea endif enddo C...skip over ONTO grid element table do k=1,neo read(14,*) end do C...ONTO grid boundary nodes nbn = 0 read(14,*) nob read(14,*) itnobn do i=1,nob read(14,*) nobn do j=1,nobn nbn = nbn + 1 read(14,*) boundnode(nbn) end do end do read(14,*) nlb read(14,*) itnlbn do i=1,nlb read(14,*) nlbn do j=1,nlbn nbn = nbn + 1 read(14,*) boundnode(nbn) end do end do close(14) C...Process ONTO grid locations write(*,*) ' Processing ONTO grid nodes................' write(*,*) ' ' ibn = 0 do i=1,npo aemin=1.0d+25 x4=xo(i) y4=yo(i) do k=1,nef x1=xf(nmf1(k)) x2=xf(nmf2(k)) x3=xf(nmf3(k)) y1=yf(nmf1(k)) y2=yf(nmf2(k)) y3=yf(nmf3(k)) a1=(x1-x4)*(y2-y4) - (x2-x4)*(y1-y4) a2=(x2-x4)*(y3-y4) - (x3-x4)*(y2-y4) a3=(x1-x3)*(y4-y3) - (x4-x3)*(y1-y3) aa=ABS(a1)+ABS(a2)+ABS(a3) ae=ABS(aa-areaf(k))/areaf(k) if(ae.lt.aemin) then aemin=ae nnef(i)=-k endif if(ae.lt.1.0E-10) then nnef(i)=k goto 100 ! LOCATION IS IN THIS ELEMENT endif end do ibn = ibn + 1 badnode(ibn) = i badprox(ibn) = aemin write(*,2000) i,aemin 2000 format(//,' WARNING: node ',I6,' in the ONTO grid is not', & ' inside any FROM grid element.',/, & ' Check the longitude and latitude for this', & ' location.',/,' Program will estimate nearest', & ' element',/,' The proximity index for this', & ' node equals ',E15.6,//) write(11,2001) badnode(ibn),badprox(ibn),abs(nnef(i)) 2001 format(1x,I10,E15.6,I10) 100 k=abs(nnef(i)) x1=xf(nmf1(k)) x2=xf(nmf2(k)) x3=xf(nmf3(k)) y1=yf(nmf1(k)) y2=yf(nmf2(k)) y3=yf(nmf3(k)) sta2(i)=((x4-x1)*(y3-y1)-(y4-y1)*(x3-x1))/areaf(k) sta3(i)=(-(x4-x1)*(y2-y1)+(y4-y1)*(x2-x1))/areaf(k) sta1(i)=1.0-sta2(i)-sta3(i) end do C...Open and write the interpolation output file (UNIT 2) open(2,file=intoutf) write(*,*) ' Writing interpolation output file................' write(*,*) ' ' write(2,'(2(a24,5x))') agridf,agrido write(2,*) npo do i=1,npo k=abs(nnef(i)) if(nnef(i).gt.0) then write(2,3000) i,nmf1(k),nmf2(k),nmf3(k),sta1(i),sta2(i), & sta3(i) else write(2,3000) i,-nmf1(k),-nmf2(k),-nmf3(k),sta1(i),sta2(i), & sta3(i) endif 3000 format(4(1x,I7),3(1x,e14.7)) end do C Determine if ONTO nodes that lie outside of the FROM grid are C boundary nodes in the ONTO grid if(ibn.gt.0) then write(*,*) ' ' write(*,*) 'Determining if bad nodes are boundary nodes ......' write(11,*) ' ' endif do 333 i=1,ibn do j=1,nbn if(badnode(i).eq.boundnode(j)) then write(11,*) ' Node ',badnode(i),' is a boundary node' goto 333 endif end do write(11,*) ' Node ',badnode(i),' is NOT a boundary node' 333 continue close(2) close(11) stop end SUBROUTINE CPP(X,Y,RLAMBDA,PHI,RLAMBDA0,PHI0) double precision x,y,rlambda,phi,rlambda0,phi0,r R=6378206.4d0 X=R*(RLAMBDA-RLAMBDA0)*COS(PHI0) Y=PHI*R RETURN END