-
Notifications
You must be signed in to change notification settings - Fork 0
/
module_init_utilities.f90
95 lines (79 loc) · 2.68 KB
/
module_init_utilities.f90
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
MODULE module_init_utilities
CONTAINS
real function interp_0( v_in, &
z_in, z_out, nz_in )
implicit none
integer nz_in, nz_out
real v_in(nz_in), z_in(nz_in)
real z_out
integer kp, k, im, ip
logical interp, increasing_z
real height, w1, w2
logical debug
parameter ( debug = .false. )
! does vertical coordinate increase or decrease with increasing k?
! set offset appropriately
height = z_out
if(debug) write(6,*) ' height in interp_0 ',height
if (z_in(nz_in) .gt. z_in(1)) then
if(debug) write(6,*) ' monotonic increase in z in interp_0 '
IF (height > z_in(nz_in)) then
if(debug) write(6,*) ' point 1 in interp_0 '
w2 = (z_in(nz_in)-height)/(z_in(nz_in)-z_in(nz_in-1))
w1 = 1.-w2
interp_0 = w1*v_in(nz_in) + w2*v_in(nz_in-1)
ELSE IF (height < z_in(1)) then
if(debug) write(6,*) ' point 2 in interp_0 '
w2 = (z_in(2)-height)/(z_in(2)-z_in(1))
w1 = 1.-w2
interp_0 = w1*v_in(2) + w2*v_in(1)
ELSE
if(debug) write(6,*) ' point 3 in interp_0 '
interp = .false.
kp = nz_in
DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) )
IF( ((z_in(kp) .ge. height) .and. &
(z_in(kp-1) .le. height)) ) THEN
w2 = (height-z_in(kp))/(z_in(kp-1)-z_in(kp))
w1 = 1.-w2
interp_0 = w1*v_in(kp) + w2*v_in(kp-1)
if(debug) write(6,*) ' interp data, kp, w1, w2 ',kp, w1, w2
if(debug) write(6,*) ' interp data, v_in(kp), v_in(kp-1), interp_0 ', &
v_in(kp), v_in(kp-1), interp_0
interp = .true.
END IF
kp = kp-1
ENDDO
ENDIF
else
if(debug) write(6,*) ' monotonic decrease in z in interp_0 '
IF (height < z_in(nz_in)) then
if(debug) write(6,*) ' point 1 in interp_0 '
w2 = (z_in(nz_in)-height)/(z_in(nz_in)-z_in(nz_in-1))
w1 = 1.-w2
interp_0 = w1*v_in(nz_in) + w2*v_in(nz_in-1)
ELSE IF (height > z_in(1)) then
if(debug) write(6,*) ' point 2 in interp_0 '
w2 = (z_in(2)-height)/(z_in(2)-z_in(1))
w1 = 1.-w2
interp_0 = w1*v_in(2) + w2*v_in(1)
ELSE
if(debug) write(6,*) ' point 3 in interp_0 '
interp = .false.
kp = nz_in
height = z_out
DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) )
IF( ((z_in(kp) .le. height) .and. &
(z_in(kp-1) .ge. height)) ) THEN
w2 = (height-z_in(kp))/(z_in(kp-1)-z_in(kp))
w1 = 1.-w2
interp_0 = w1*v_in(kp) + w2*v_in(kp-1)
interp = .true.
END IF
kp = kp-1
ENDDO
ENDIF
end if
return
END FUNCTION interp_0
END MODULE module_init_utilities