Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
! Created by Andrey Debolskiy on 15.08.2024.
module external_forcing
implicit none
!type declarations
type forcing_1d
real, allocatable :: x(:) ! 0 - closed, 1 - opened
real, allocatable :: val(:)
integer :: xsize
end type forcing_1d
type forcing_2d
real, allocatable :: x(:) ! 0 - closed, 1 - opened
real, allocatable :: y(:) ! 0 - closed, 1 - opened
real, allocatable :: val(:,:)
integer :: xsize, ysize
end type forcing_2d
public :: ext_forc_dealloc_1d, ext_forc_dealloc_2d
public :: get_1d_interpolated, get_1d_to_1d_interpolated, &
read_2d_focing
contains
real function get_1d_interpolated(xq, forcing1d) result(val)
real, intent(in) :: xq
type (forcing_1d), intent(in) :: forcing1d
!local
integer i
real x1,x2, val1, val2
if (xq < forcing1d%x(1)) then
val = forcing1d%val(1)
return
end if
do i = 1, forcing1d%xsize
if (forcing1d%x(i) >= xq) exit
end do
if (i == forcing1d%xsize) then
val = forcing1d%val(i)
return
else
x1 = forcing1d%x(i-1)
x2 = forcing1d%x(i)
val1 = forcing1d%val(i-1)
val2 = forcing1d%val(i)
val = val1 + (val2 - val1) * (xq - x1) / (x2 -x1)
return
end if
end function get_1d_interpolated
subroutine get_1d_to_1d_interpolated(valq, xq, nq, forcing1d)
integer, intent(in):: nq
real, dimension(nq), intent(in) :: xq
real, dimension(nq), intent(out) :: valq
type (forcing_1d), intent(in) :: forcing1d
!local
integer k
do k = 1, nq
valq(k) = get_1d_interpolated(xq(k), forcing1d)
enddo
end subroutine get_1d_to_1d_interpolated
subroutine read_2d_focing(out_forc, filename, filename_x, filename_y)
implicit none
type (forcing_2d), intent(out):: out_forc
character(*), intent(in):: filename, filename_x, filename_y
!local vars
integer nx,ny, i, j
real x,y
open(10, FILE=trim(filename))
read(10,*) nx
read(10,*) ny
out_forc%xsize=nx
out_forc%ysize=ny
allocate(out_forc%val(nx,ny))
do i = 1, nx
READ(10,*) (out_forc%val(i,j),j=1,ny)
end do
close(10)
! read coords
allocate(out_forc%x(nx))
open(10, FILE=trim(filename_x))
do i=1,nx
read(10,*) out_forc%x(i)
end do
close(10)
allocate(out_forc%y(ny))
open(10, FILE=trim(filename_y))
do i=1,ny
read(10,*) out_forc%y(i)
end do
close(10)
end subroutine read_2d_focing
subroutine ext_forc_dealloc_1d(forc_1d)
implicit none
type (forcing_1d), intent(inout) :: forc_1d
deallocate(forc_1d%val)
deallocate(forc_1d%x)
end subroutine ext_forc_dealloc_1d
subroutine ext_forc_dealloc_2d(forc_2d)
implicit none
type (forcing_2d), intent(inout) :: forc_2d
deallocate(forc_2d%val)
deallocate(forc_2d%x)
deallocate(forc_2d%y)
end subroutine ext_forc_dealloc_2d
end module external_forcing