Subroutine vert_interp 12VINTPF1A.22
& (data_in, data_points, VINTPF1A.23
& data_levels, desired_r, VINTPF1A.24
& r_at_data, interp_order, VINTPF1A.25
& data_out ) VINTPF1A.26
VINTPF1A.27
! Purpose: VINTPF1A.28
! Performs vertical interpolation of a field to a VINTPF1A.29
! desired r surface given the r value at each of VINTPF1A.30
! the data points. Where the desired surface is below/above VINTPF1A.31
! the bottom/top data point a missing data indicator is VINTPF1A.32
! returned. VINTPF1A.33
! r can be any quantity that is monotonic increasing with VINTPF1A.34
! model level, such as height, or potential temperature VINTPF1A.35
! (assuming that the model is statically stable). VINTPF1A.36
! VINTPF1A.37
! Method: VINTPF1A.38
! Cubic Lagrange interpolation. VINTPF1A.39
! VINTPF1A.40
! Original Progammer: Mark H. Mawson VINTPF1A.41
! Current code owner: I Edmond VINTPF1A.42
! VINTPF1A.43
! History: VINTPF1A.44
! Date Version Comment VINTPF1A.45
! ---- ------- ------- VINTPF1A.46
! VINTPF1A.47
! 22/5/96 4.1 Linear extrapolation added to routine to VINTPF1A.48
! obtain data at points above/below VINTPF1A.49
! top/bottom levels. Ian Edmond VINTPF1A.50
! 29/07/98 4.5 Optimisation changes for T3E UDG5F405.594
! Author D.M. Goddard UDG5F405.595
! VINTPF1A.51
! Code Description: VINTPF1A.52
! Language: FORTRAN 77 + CRAY extensions VINTPF1A.53
! This code is written to UMDP3 programming standards. VINTPF1A.54
! VINTPF1A.55
! System component covered: ?? VINTPF1A.56
! System Task: ?? VINTPF1A.57
! VINTPF1A.58
Implicit None VINTPF1A.59
VINTPF1A.60
! Arguments with Intent IN. ie: Input variables. VINTPF1A.61
VINTPF1A.62
Integer VINTPF1A.63
& data_points ! number of rows of data VINTPF1A.64
&, data_levels ! number of levels of data VINTPF1A.65
&, interp_order ! 1 = linear, 3 = cubic, 5=quintic VINTPF1A.66
VINTPF1A.67
Real VINTPF1A.68
& data_in (data_points, data_levels) VINTPF1A.69
&, r_at_data (data_points, data_levels) VINTPF1A.70
&, desired_r(data_points) ! desired value to which VINTPF1A.71
& ! data should be interpolated t VINTPF1A.72
VINTPF1A.73
! Arguments with Intent OUT. ie: Output variables. VINTPF1A.74
Real VINTPF1A.75
& data_out (data_points) VINTPF1A.76
VINTPF1A.77
! Local variables VINTPF1A.78
VINTPF1A.79
Integer VINTPF1A.80
& j,k VINTPF1A.81
VINTPF1A.82
Integer last UDG5F405.596
Integer VINTPF1A.83
& level_below(data_points) VINTPF1A.84
VINTPF1A.85
Real VINTPF1A.86
& r_here VINTPF1A.87
&, r_here_plus VINTPF1A.88
&, r_here_plus2 VINTPF1A.89
&, r_here_minus VINTPF1A.90
&, r_here_plus3 VINTPF1A.91
&, r_here_minus2 VINTPF1A.92
VINTPF1A.93
! ---------------------------------------------------------------------- VINTPF1A.94
! Section 1. Find level below which desired surface is. VINTPF1A.95
! ---------------------------------------------------------------------- VINTPF1A.96
VINTPF1A.97
last=1 UDG5F405.597
Do j=1, data_points UDG5F405.598
If(r_at_data(j,last) .le .desired_r(j))then UDG5F405.599
Do k=last+1, data_levels UDG5F405.600
If(r_at_data(j,k) .gt .desired_r(j))then UDG5F405.601
level_below(j)=k-1 UDG5F405.602
go to 9876 UDG5F405.603
Endif UDG5F405.604
Enddo UDG5F405.605
level_below(j)=data_levels UDG5F405.606
Else UDG5F405.607
Do k=last-1, 1, -1 UDG5F405.608
If(r_at_data(j,k) .le .desired_r(j))then UDG5F405.609
level_below(j)=k UDG5F405.610
go to 9876 UDG5F405.611
Endif UDG5F405.612
Enddo UDG5F405.613
level_below(j)=0 UDG5F405.614
Endif UDG5F405.615
9876 continue UDG5F405.616
LAST=Max(level_below(j),1) UDG5F405.617
End do UDG5F405.618
VINTPF1A.109
VINTPF1A.110
Do j = 1, data_points VINTPF1A.111
VINTPF1A.112
! If requested level is above top of model, do linear VINTPF1A.113
! extrapolation using data on top and second top levels. VINTPF1A.114
If ( desired_r(j) .gt. r_at_data(j,data_levels) ) Then VINTPF1A.115
data_out(j) = data_in(j,data_levels) + (desired_r(j) VINTPF1A.116
& - r_at_data(j,data_levels)) * (data_in(j,data_levels) VINTPF1A.117
& - data_in(j,data_levels-1))/(r_at_data(j,data_levels) VINTPF1A.118
& - r_at_data(j,data_levels-1)) VINTPF1A.119
End If VINTPF1A.120
VINTPF1A.121
! If requested level is below bottom of model, do linear VINTPF1A.122
! extrapolation using data on first and second levels. VINTPF1A.123
If ( desired_r(j) .lt. r_at_data(j,1) ) Then VINTPF1A.124
data_out(j) = data_in(j,1) + (desired_r(j) VINTPF1A.125
& - r_at_data(j,1)) * (data_in(j,1) VINTPF1A.126
& - data_in(j,2))/(r_at_data(j,1) - r_at_data(j,2)) VINTPF1A.127
End If VINTPF1A.128
VINTPF1A.129
End Do VINTPF1A.130
VINTPF1A.131
! ---------------------------------------------------------------------- VINTPF1A.132
! Section 2. Vertical interpolation. VINTPF1A.133
! ---------------------------------------------------------------------- VINTPF1A.134
VINTPF1A.135
Do j = 1, data_points VINTPF1A.136
VINTPF1A.137
If (level_below(j) .eq. 0.or. UDG5F405.619
& level_below(j) .eq. data_levels) Then UDG5F405.620
UDG5F405.621
UDG5F405.622
data_out (j) = data_out (j) VINTPF1A.139
VINTPF1A.140
Else if (level_below(j) .eq. 1 .or. VINTPF1A.141
& level_below(j) .eq. data_levels - 1 VINTPF1A.142
& .or. interp_order .eq. 1 ) Then VINTPF1A.143
! linearly interpolate VINTPF1A.144
data_out (j) = ( (desired_r(j) - VINTPF1A.145
& r_at_data(j,level_below(j)) ) VINTPF1A.146
& * data_in (j,level_below(j)+1) VINTPF1A.147
& -(desired_r(j) - VINTPF1A.148
& r_at_data(j,level_below(j)+1)) * VINTPF1A.149
& data_in (j,level_below(j)) ) / VINTPF1A.150
& ( r_at_data(j,level_below(j)+1) - VINTPF1A.151
& r_at_data(j,level_below(j)) ) VINTPF1A.152
VINTPF1A.153
Else if (level_below(j) .eq. 2 .or. VINTPF1A.154
& level_below(j) .eq. data_levels - 2 VINTPF1A.155
& .or. interp_order .eq. 3 ) Then VINTPF1A.156
VINTPF1A.157
! cubicly interpolate VINTPF1A.158
VINTPF1A.159
r_here_minus = r_at_data(j,level_below(j)-1) VINTPF1A.160
r_here = r_at_data(j,level_below(j)) VINTPF1A.161
r_here_plus = r_at_data(j,level_below(j)+1) VINTPF1A.162
r_here_plus2 = r_at_data(j,level_below(j)+2) VINTPF1A.163
VINTPF1A.164
data_out (j) = ( (desired_r(j) - r_here) * VINTPF1A.165
& (desired_r(j) - r_here_plus )* VINTPF1A.166
& (desired_r(j) - r_here_plus2 ) ) / VINTPF1A.167
& ( (r_here_minus - r_here) * VINTPF1A.168
& (r_here_minus - r_here_plus )* VINTPF1A.169
& (r_here_minus - r_here_plus2 ) ) * VINTPF1A.170
& data_in (j,level_below(j)-1) + VINTPF1A.171
& ( (desired_r(j) - r_here_minus) * VINTPF1A.172
& (desired_r(j) - r_here_plus )* VINTPF1A.173
& (desired_r(j) - r_here_plus2 ) ) / VINTPF1A.174
& ( (r_here - r_here_minus) * VINTPF1A.175
& (r_here - r_here_plus )* VINTPF1A.176
& (r_here - r_here_plus2 ) ) * VINTPF1A.177
& data_in (j,level_below(j)) + VINTPF1A.178
& ( (desired_r(j) - r_here_minus) * VINTPF1A.179
& (desired_r(j) - r_here )* VINTPF1A.180
& (desired_r(j) - r_here_plus2 ) ) / VINTPF1A.181
& ( (r_here_plus - r_here_minus) * VINTPF1A.182
& (r_here_plus - r_here )* VINTPF1A.183
& (r_here_plus - r_here_plus2 ) ) * VINTPF1A.184
& data_in (j,level_below(j)+1) + VINTPF1A.185
& ( (desired_r(j) - r_here_minus) * VINTPF1A.186
& (desired_r(j) - r_here )* VINTPF1A.187
& (desired_r(j) - r_here_plus ) ) / VINTPF1A.188
& ( (r_here_plus2 - r_here_minus) * VINTPF1A.189
& (r_here_plus2 - r_here )* VINTPF1A.190
& (r_here_plus2 - r_here_plus ) ) * VINTPF1A.191
& data_in (j,level_below(j)+2) VINTPF1A.192
VINTPF1A.193
Else VINTPF1A.194
! interpolate quinticly VINTPF1A.195
VINTPF1A.196
r_here_minus2 = r_at_data(j,level_below(j)-2) VINTPF1A.197
r_here_minus = r_at_data(j,level_below(j)-1) VINTPF1A.198
r_here = r_at_data(j,level_below(j)) VINTPF1A.199
r_here_plus = r_at_data(j,level_below(j)+1) VINTPF1A.200
r_here_plus2 = r_at_data(j,level_below(j)+2) VINTPF1A.201
r_here_plus3 = r_at_data(j,level_below(j)+3) VINTPF1A.202
VINTPF1A.203
Data_out (j) = ((desired_r(j) - r_here_minus) * VINTPF1A.204
& (desired_r(j) - r_here )* VINTPF1A.205
& (desired_r(j) - r_here_plus )* VINTPF1A.206
& (desired_r(j) - r_here_plus2 )* VINTPF1A.207
& (desired_r(j) - r_here_plus3 ))/ VINTPF1A.208
& ( (r_here_minus2 - r_here_minus) * VINTPF1A.209
& (r_here_minus2 - r_here )* VINTPF1A.210
& (r_here_minus2 - r_here_plus )* VINTPF1A.211
& (r_here_minus2 - r_here_plus2 )* VINTPF1A.212
& (r_here_minus2 - r_here_plus3 ) ) * VINTPF1A.213
& data_in (j,level_below(j)-2) + VINTPF1A.214
& ((desired_r(j) - r_here_minus2) * VINTPF1A.215
& (desired_r(j) - r_here )* VINTPF1A.216
& (desired_r(j) - r_here_plus )* VINTPF1A.217
& (desired_r(j) - r_here_plus2 )* VINTPF1A.218
& (desired_r(j) - r_here_plus3 ))/ VINTPF1A.219
& ( (r_here_minus - r_here_minus2) * VINTPF1A.220
& (r_here_minus - r_here )* VINTPF1A.221
& (r_here_minus - r_here_plus )* VINTPF1A.222
& (r_here_minus - r_here_plus2 )* VINTPF1A.223
& (r_here_minus - r_here_plus3 ) ) * VINTPF1A.224
& data_in (j,level_below(j)-1) + VINTPF1A.225
& ((desired_r(j) - r_here_minus2) * VINTPF1A.226
& (desired_r(j) - r_here_minus )* VINTPF1A.227
& (desired_r(j) - r_here_plus )* VINTPF1A.228
& (desired_r(j) - r_here_plus2 )* VINTPF1A.229
& (desired_r(j) - r_here_plus3 ))/ VINTPF1A.230
& ( (r_here - r_here_minus2) * VINTPF1A.231
& (r_here - r_here_minus )* VINTPF1A.232
& (r_here - r_here_plus )* VINTPF1A.233
& (r_here - r_here_plus2 )* VINTPF1A.234
& (r_here - r_here_plus3 ) ) * VINTPF1A.235
& data_in (j,level_below(j)) + VINTPF1A.236
& ((desired_r(j) - r_here_minus2) * VINTPF1A.237
& (desired_r(j) - r_here_minus )* VINTPF1A.238
& (desired_r(j) - r_here )* VINTPF1A.239
& (desired_r(j) - r_here_plus2 )* VINTPF1A.240
& (desired_r(j) - r_here_plus3 ))/ VINTPF1A.241
& ( (r_here_plus - r_here_minus2) * VINTPF1A.242
& (r_here_plus - r_here_minus )* VINTPF1A.243
& (r_here_plus - r_here )* VINTPF1A.244
& (r_here_plus - r_here_plus2 )* VINTPF1A.245
& (r_here_plus - r_here_plus3 ) ) * VINTPF1A.246
& data_in (j,level_below(j)+1) + VINTPF1A.247
& ((desired_r(j) - r_here_minus2) * VINTPF1A.248
& (desired_r(j) - r_here_minus )* VINTPF1A.249
& (desired_r(j) - r_here )* VINTPF1A.250
& (desired_r(j) - r_here_plus )* VINTPF1A.251
& (desired_r(j) - r_here_plus3 ))/ VINTPF1A.252
& ( (r_here_plus2 - r_here_minus2) * VINTPF1A.253
& (r_here_plus2 - r_here_minus )* VINTPF1A.254
& (r_here_plus2 - r_here )* VINTPF1A.255
& (r_here_plus2 - r_here_plus )* VINTPF1A.256
& (r_here_plus2 - r_here_plus3 ) ) * VINTPF1A.257
& data_in (j,level_below(j)+2) + VINTPF1A.258
& ((desired_r(j) - r_here_minus2) * VINTPF1A.259
& (desired_r(j) - r_here_minus )* VINTPF1A.260
& (desired_r(j) - r_here )* VINTPF1A.261
& (desired_r(j) - r_here_plus )* VINTPF1A.262
& (desired_r(j) - r_here_plus2 ))/ VINTPF1A.263
& ( (r_here_plus3 - r_here_minus2) * VINTPF1A.264
& (r_here_plus3 - r_here_minus )* VINTPF1A.265
& (r_here_plus3 - r_here )* VINTPF1A.266
& (r_here_plus3 - r_here_plus )* VINTPF1A.267
& (r_here_plus3 - r_here_plus2 ) ) * VINTPF1A.268
& data_in (j,level_below(j)+3) VINTPF1A.269
VINTPF1A.270
End If VINTPF1A.271
VINTPF1A.272
End Do VINTPF1A.273
VINTPF1A.274
! end of routine VINTPF1A.275
VINTPF1A.276
Return VINTPF1A.277
End VINTPF1A.278
*ENDIF VINTPF1A.279