Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
 Path: blob/devel/elmerice/examples/InverseMethods_OLD/src/USF_Init.F90
3206 views
1
        FUNCTION UIni( Model, nodenumber, dumy) RESULT(U)

2
        USE types

3


4
        implicit none

5
        TYPE(Model_t) :: Model

6
        Real(kind=dp) :: dumy,U

7
        INTEGER :: nodenumber

8


9
        Real(kind=dp),allocatable :: dem(:,:),xx(:),yy(:)

10
        Real(kind=dp) :: x,y

11
        Real(kind=dp) :: LinearInterp

12


13
        integer :: nx,ny

14
        integer :: i,j

15


16
        character(len=MAX_NAME_LEN) :: filin='Data/UDEM.xy'

17


18
        logical :: Firsttime=.true.

19


20
        SAVE dem,xx,yy,nx,ny

21
        SAVE Firsttime

22


23
        if (Firsttime) then

24
                Firsttime=.False.

25


26
        ! open file

27
                open(10,file=trim(filin))

28
                Read(10,*) nx

29
                Read(10,*) ny

30
                allocate(xx(nx),yy(ny))

31
                Allocate(dem(nx,ny))

32
                Do i=1,nx

33
                   Do j=1,ny

34
                      read(10,*) xx(i),yy(j),dem(i,j)

35
                   End Do

36
                End do

37
                close(10)

38
        End if

39


40
        ! position current point

41
        x=Model % Nodes % x (nodenumber)

42
        y=Model % Nodes % y (nodenumber)

43


44
        U=LinearInterp(dem,xx,yy,nx,ny,x,y)

45


46
        Return 

47
        End

48


49
        FUNCTION VIni( Model, nodenumber, dumy) RESULT(U)

50
        USE types

51


52
        implicit none

53
        TYPE(Model_t) :: Model

54
        Real(kind=dp) :: dumy,U

55
        INTEGER :: nodenumber

56


57
        Real(kind=dp),allocatable :: dem(:,:),xx(:),yy(:)

58
        Real(kind=dp) :: x,y

59
        Real(kind=dp) :: LinearInterp

60


61
        integer :: nx,ny

62
        integer :: i,j

63


64
        character(len=MAX_NAME_LEN) :: filin='Data/VDEM.xy'

65


66
        logical :: Firsttime=.true.

67


68
        SAVE dem,xx,yy,nx,ny

69
        SAVE Firsttime

70


71
        if (Firsttime) then

72
                Firsttime=.False.

73


74
        ! open file

75
                open(10,file=trim(filin))

76
                Read(10,*) nx

77
                Read(10,*) ny

78
                allocate(xx(nx),yy(ny))

79
                Allocate(dem(nx,ny))

80
                Do i=1,nx

81
                   Do j=1,ny

82
                      read(10,*) xx(i),yy(j),dem(i,j)

83
                   End Do

84
                End do

85
                close(10)

86
        End if

87


88
        ! position current point

89
        x=Model % Nodes % x (nodenumber)

90
        y=Model % Nodes % y (nodenumber)

91


92
        U=LinearInterp(dem,xx,yy,nx,ny,x,y)

93


94
        Return 

95
        End

96


97
        FUNCTION zsIni( Model, nodenumber, dumy) RESULT(U)

98
        USE types

99


100
        implicit none

101
        TYPE(Model_t) :: Model

102
        Real(kind=dp) :: dumy,U

103
        INTEGER :: nodenumber

104


105
        Real(kind=dp),allocatable :: dem(:,:),xx(:),yy(:)

106
        Real(kind=dp) :: x,y

107
        Real(kind=dp) :: LinearInterp

108


109
        integer :: nx,ny

110
        integer :: i,j

111


112
        character(len=MAX_NAME_LEN) :: filin='Data/zsDEM.xy'

113


114
        logical :: Firsttime=.true.

115


116
        SAVE dem,xx,yy,nx,ny

117
        SAVE Firsttime

118


119
        if (Firsttime) then

120
                Firsttime=.False.

121


122
        ! open file

123
                open(10,file=trim(filin))

124
                Read(10,*) nx

125
                Read(10,*) ny

126
                allocate(xx(nx),yy(ny))

127
                Allocate(dem(nx,ny))

128
                Do i=1,nx

129
                   Do j=1,ny

130
                      read(10,*) xx(i),yy(j),dem(i,j)

131
                   End Do

132
                End do

133
                close(10)

134
        End if

135


136
        ! position current point

137
        x=Model % Nodes % x (nodenumber)

138
        y=Model % Nodes % y (nodenumber)

139


140
        U=LinearInterp(dem,xx,yy,nx,ny,x,y)

141


142
        Return 

143
        End

144


145
        FUNCTION zbIni( Model, nodenumber, dumy) RESULT(U)

146
        USE types

147


148
        implicit none

149
        TYPE(Model_t) :: Model

150
        Real(kind=dp) :: dumy,U

151
        INTEGER :: nodenumber

152


153
        Real(kind=dp),allocatable :: dem(:,:),xx(:),yy(:)

154
        Real(kind=dp) :: x,y

155
        Real(kind=dp) :: LinearInterp

156


157
        integer :: nx,ny

158
        integer :: i,j

159


160
        character(len=MAX_NAME_LEN) :: filin='Data/zbDEM.xy'

161


162
        logical :: Firsttime=.true.

163


164
        SAVE dem,xx,yy,nx,ny

165
        SAVE Firsttime

166


167
        if (Firsttime) then

168
                Firsttime=.False.

169


170
        ! open file

171
                open(10,file=trim(filin))

172
                Read(10,*) nx

173
                Read(10,*) ny

174
                allocate(xx(nx),yy(ny))

175
                Allocate(dem(nx,ny))

176
                Do i=1,nx

177
                   Do j=1,ny

178
                      read(10,*) xx(i),yy(j),dem(i,j)

179
                   End Do

180
                End do

181
                close(10)

182
        End if

183


184
        ! position current point

185
        x=Model % Nodes % x (nodenumber)

186
        y=Model % Nodes % y (nodenumber)

187


188
        U=LinearInterp(dem,xx,yy,nx,ny,x,y)

189


190
        Return 

191
        End

192


193


194
        Function LinearInterp(dem,xx,yy,nx,ny,x,y) Result(InterP1)

195


196
        USE TYPES

197
        implicit none

198
        

199
        REAL(KIND=dp) :: dem(nx,ny),xx(nx),yy(ny)

200
        REAL(KIND=dp) :: Dx,Dy,DxDy

201
        Real(kind=dp) :: x,y,x_1,y_1,B(4)

202
        Real(kind=dp) :: InterP1

203
        integer :: nx,ny

204
        integer :: nx_1,ny_1

205


206
        Dx=(xx(nx)-xx(1))/(nx-1)

207
        Dy=(yy(ny)-yy(1))/(ny-1)

208
        DxDy=Dx*Dy

209


210
        ! lower left point in DEM

211
        nx_1=floor((x-xx(1))/Dx) + 1

212
        ny_1=floor((y-yy(1))/Dy) + 1

213
        nx_1=min(nx_1,nx-1)

214
        ny_1=min(ny_1,ny-1)

215


216
        x_1=xx(nx_1)

217
        y_1=yy(ny_1)

218


219


220
        ! DEM Value in surroundings points

221
        !       4 ----- 3

222
        !       |       |

223
        !       1 ----- 2

224
        B(1)=dem(nx_1,ny_1)

225
        B(2)=dem(nx_1+1,ny_1)

226
        B(3)=dem(nx_1+1,ny_1+1)

227
        B(4)=dem(nx_1,ny_1+1)

228


229
        

230
        ! Linear Interpolation at Point x,y

231
        InterP1=(x-x_1)*(y-y_1)*(B(3)+B(1)-B(2)-B(4))/DxDy

232
        InterP1=InterP1+(x-x_1)*(B(2)-B(1))/Dx+(y-y_1)*(B(4)-B(1))/Dy+B(1)

233


234
        Return

235
        End

236


237