SUBROUTINE HorizontalInterp( ,5HORINT1A.3
& Scheme, LGlobal, LVector, VSB2F405.1
& Len1In, Len2In, HORINT1A.5
& Len1Out, Len2Out, HORINT1A.6
& LambdaOrigin, HORINT1A.7
& PhiOrigin, HORINT1A.8
& DLambdaIn, HORINT1A.9
& DPhiIn, HORINT1A.10
& LambdaOut, HORINT1A.11
& PhiOut, HORINT1A.12
& DataIn, HORINT1A.13
& DataOut, HORINT1A.14
& ErrorStatus,ErrorMessage) HORINT1A.15
HORINT1A.16
! Description: HORINT1A.17
! Interpolate by linear, cubic or quintic lagrange interpolation, a 2d HORINT1A.18
! field to a 2d set of points. Input can be on a sphere or within a HORINT1A.19
! rectangular box. Requested output points must lie within the region HORINT1A.20
! defined for the input data. HORINT1A.21
! An option for Monotonicity can be enforced (currently disabled). HORINT1A.22
! HORINT1A.23
! i=x=lambda, j=y=phi HORINT1A.24
! assumption that nth point defined by: HORINT1A.25
! (LambdaOrigin+(n-1)*DLambdaIn , PhiOrigin+(n-1)*DPhiIn) HORINT1A.26
! assumption that input and output in radians in range HORINT1A.27
! (lambda 0 to 2pi ; phi -pi/2 to +pi/2) VSB1F401.47
! assumption that check for output positions outside input grid HORINT1A.29
! domain is done externally HORINT1A.30
! HORINT1A.31
! Method: This is a modified version of the routine interpolation HORINT1A.32
! written by Mark Mawson and described in: HORINT1A.33
! The proposed semi-Lagrangian advection scheme for the HORINT1A.34
! semi-Implicit Unified Model integration scheme. HORINT1A.35
! F.R. Division working paper No 162. HORINT1A.36
! Mark H. Mawson HORINT1A.37
! HORINT1A.38
! HORINT1A.39
! Owner: Stuart Bell HORINT1A.40
! HORINT1A.41
! History: HORINT1A.42
! Version Date Comment HORINT1A.43
! ------- ---- ------- HORINT1A.44
! 4.0 6/6/95 Equiv. to VAR code as at time of build:Stuart Bell HORINT1A.45
! 4.4 23/4/97 Relax checks so C90 dump OK on T3E:Stuart Bell VSB0F404.7
! 4.5 23/04/98 Add LVector arg. to enable Polar wind correction:SB VSB2F405.2
! HORINT1A.46
! Code Description: HORINT1A.47
! Language: Fortran 77 plus HORINT1A.48
! Software Standards: "UM and Met O standards". HORINT1A.49
! HORINT1A.50
! HORINT1A.51
! Declarations: HORINT1A.52
HORINT1A.53
IMPLICIT NONE HORINT1A.54
HORINT1A.55
*CALL C_PI 
HORINT1A.56
HORINT1A.57
!* Subroutine arguments HORINT1A.58
! Scalar arguments with :intent(in) HORINT1A.59
INTEGER Scheme ! a code saying which order HORINT1A.60
! ! of scheme to use: HORINT1A.61
! ! 1 = linear;3 = cubic;5 = quintic HORINT1A.62
LOGICAL LGlobal ! True, if global HORINT1A.63
LOGICAL LVector ! True, if wind component VSB2F405.3
HORINT1A.64
INTEGER Len1In ! Dimension of DataIn in i direction. HORINT1A.65
INTEGER Len2In ! Dimension of DataIn in j direction. HORINT1A.66
INTEGER Len1Out ! Dimension of DataOut in i direction. HORINT1A.67
INTEGER Len2Out ! Dimension of DataOut in j direction. HORINT1A.68
! ErrorStatus HORINT1A.69
INTEGER ErrorStatus HORINT1A.70
CHARACTER*256 ErrorMessage HORINT1A.71
HORINT1A.72
REAL DLambdaIn ! Holds spacing between points in HORINT1A.73
! ! the i direction for the input data field. HORINT1A.74
REAL DPhiIn ! Holds spacing between points in HORINT1A.75
! ! the j direction for the input field. HORINT1A.76
REAL LambdaOrigin ! Position of 1st point in HORINT1A.77
! ! the i direction for the input data field. HORINT1A.78
REAL PhiOrigin ! Position of 1st point in HORINT1A.79
! ! the j direction for the input field. HORINT1A.80
HORINT1A.81
! Array arguments with :intent(in) HORINT1A.82
REAL DataIn (Len1In,Len2In) ! Data to be interpolated. HORINT1A.83
HORINT1A.84
REAL LambdaOut (Len1Out,Len2Out) ! Lambda coordinate of HORINT1A.85
! ! output data on input grid HORINT1A.86
REAL PhiOut (Len1Out,Len2Out) ! Phi coordinate of HORINT1A.87
! ! output data on input grid HORINT1A.88
HORINT1A.89
! Array arguments with :intent(out) HORINT1A.90
REAL DataOut (Len1Out,Len2Out) ! Data interpolated HORINT1A.91
! ! to desired locations. HORINT1A.92
HORINT1A.93
!* End of Subroutine arguments HORINT1A.94
HORINT1A.95
! Local named constants: HORINT1A.96
INTEGER Extend ! extra points for DataExt HORINT1A.97
LOGICAL L_mono ! True, if interpolation required to be monotone HORINT1A.98
PARAMETER (Extend=3) HORINT1A.99
PARAMETER (L_mono=.false.) HORINT1A.100
HORINT1A.101
! Local scalars: HORINT1A.102
INTEGER LowerBound ! Lower Bound of DataExt in both directions HORINT1A.103
INTEGER i !} Loop HORINT1A.104
INTEGER j !} indices. HORINT1A.105
HORINT1A.106
INTEGER HalfLen1In ! Len1In/2 HORINT1A.107
INTEGER VIorder ! order of interpolation HORINT1A.108
INTEGER SwapPole ! 1 if no sign swap, -1 otherwise VSB2F405.4
HORINT1A.109
REAL HalfPi ! pi/2 HORINT1A.110
REAL TwoPi ! pi*2 HORINT1A.111
REAL RecipDLambdaIn ! 1/DLambdaIn HORINT1A.112
REAL RecipDPhiIn ! 1/DPhiIn HORINT1A.113
REAL HalfDPhiIn ! DPhiIn/2 HORINT1A.114
HORINT1A.115
REAL Test1,Test2 !used in origin checking HORINT1A.116
LOGICAL FirstOnPole ! " " " HORINT1A.117
HORINT1A.118
! Local arrays: HORINT1A.119
INTEGER IOut (Len1Out,Len2Out) HORINT1A.120
INTEGER JOut (Len1Out,Len2Out) HORINT1A.121
HORINT1A.122
REAL DataExt(1-Extend:Len1In+Extend, HORINT1A.123
& 1-Extend:Len2In+Extend) ! Data in extended HORINT1A.124
REAL DataHigh (Len1Out,Len2Out) HORINT1A.125
REAL DataMono (Len1Out,Len2Out) HORINT1A.126
REAL WtLambda (Len1Out,Len2Out) HORINT1A.127
REAL WtPhi (Len1Out,Len2Out) HORINT1A.128
REAL P1,P2 ! used for relaxing real number checks VSB0F404.8
LOGICAL LEQR ! 'Logical EQual Real' VSB0F404.9
HORINT1A.129
HORINT1A.130
!- End of header ------------------------------------------------------- HORINT1A.131
HORINT1A.132
! ---------------------------------------------------------------------- HORINT1A.133
! Section 1. Initialize HORINT1A.134
! ---------------------------------------------------------------------- HORINT1A.135
LEQR(P1,P2) = ((ABS(P1-P2)) .LT. (1.E-6*ABS(P1+P2))) VSB0F404.10
ErrorStatus = 0 HORINT1A.136
ErrorMessage= " " HORINT1A.137
HORINT1A.138
! 0.0 Set derived variables HORINT1A.139
HalfPi = Pi * 0.5 HORINT1A.140
TwoPi = Pi * 2.0 HORINT1A.141
RecipDLambdaIn = 1.0 / DLambdaIn HORINT1A.142
RecipDPhiIn = 1.0 / DPhiIn HORINT1A.143
HalfDPhiIn = DPhiIn * 0.5 HORINT1A.144
HORINT1A.145
! 1.0 Error trap unsupported options. HORINT1A.146
HORINT1A.147
! 1.1 Check that a valid interpolation scheme has been specified. HORINT1A.148
If (( Scheme .ne. 1 .and. HORINT1A.149
& Scheme .ne. 3 .and. HORINT1A.150
& Scheme .ne. 5 ) ) then HORINT1A.151
HORINT1A.152
ErrorStatus=-1 HORINT1A.153
ErrorMessage="Invalid value of 'Scheme':linear interp. done " HORINT1A.154
HORINT1A.155
VIorder = 1 HORINT1A.156
HORINT1A.157
Else HORINT1A.158
HORINT1A.159
VIorder = Scheme HORINT1A.160
HORINT1A.161
End If HORINT1A.162
HORINT1A.163
! 1.2 Check code matches Extend parameter HORINT1A.164
IF (Extend .ne. 3) THEN HORINT1A.165
HORINT1A.166
ErrorStatus=1 HORINT1A.167
ErrorMessage="Don't change 'Extend' without amending code " HORINT1A.168
GOTO 999 HORINT1A.169
HORINT1A.170
END IF HORINT1A.171
HORINT1A.172
! 1.3 For global case check grid offset from south pole (PhiOrigin) HORINT1A.173
! using the local function LEQR to test near equal reals VSB0F404.11
IF (LGlobal) THEN HORINT1A.174
Test1=-HalfPi VSB0F404.12
Test2=-HalfPi+HalfDPhiIn VSB0F404.13
IF(LEQR(PhiOrigin,Test1)) THEN VSB0F404.14
! origin at pole HORINT1A.180
FirstOnPole = .True. HORINT1A.181
ELSEIF(LEQR(PhiOrigin,Test2)) THEN VSB0F404.15
! origin half a grid length from pole HORINT1A.183
FirstOnPole = .False. HORINT1A.184
ELSE VSB0F404.16
! unexpected origin HORINT1A.186
ErrorStatus=1 HORINT1A.187
ErrorMessage="HORINT1A: Incorrect grid " VSB0F404.17
GOTO 999 HORINT1A.189
END IF HORINT1A.190
END IF HORINT1A.191
HORINT1A.192
! 1.4 Get lower bound for DataExt HORINT1A.193
LowerBound=1-Extend HORINT1A.194
! Len1In = SIZE (DataIn, 1) HORINT1A.195
! Len2In = SIZE (DataIn, 2) HORINT1A.196
! Len1Out = SIZE (DataOut, 1) HORINT1A.197
! Len2Out = SIZE (DataOut, 2) HORINT1A.198
HORINT1A.199
!----------------------------------------------------------------------- HORINT1A.200
! 2.0 Extend input data array to bigger area to allow interpolation to HORINT1A.201
! be done without having to redo any end points. HORINT1A.202
!----------------------------------------------------------------------- HORINT1A.203
HORINT1A.204
HalfLen1In = Len1In/2 HORINT1A.205
HORINT1A.206
! 2.1 Extend in i direction. HORINT1A.207
IF (LGlobal) THEN ! 2.1.1 Global HORINT1A.208
DO j = 1, Len2In HORINT1A.209
DataExt (-2,j) = DataIn (Len1In-2,j) HORINT1A.210
DataExt (-1,j) = DataIn (Len1In-1,j) HORINT1A.211
DataExt (0,j) = DataIn (Len1In,j) HORINT1A.212
DataExt (Len1In+1,j) = DataIn (1,j) HORINT1A.213
DataExt (Len1In+2,j) = DataIn (2,j) HORINT1A.214
DataExt (Len1In+3,j) = DataIn (3,j) HORINT1A.215
END DO HORINT1A.216
HORINT1A.217
ELSE ! 2.1.2 Limited Area. HORINT1A.218
DO j = 1, Len2In HORINT1A.219
DataExt (-2,j) = DataIn (1,j) HORINT1A.220
DataExt (-1,j) = DataIn (1,j) HORINT1A.221
DataExt (0,j) = DataIn (1,j) HORINT1A.222
DataExt (Len1In+1,j) = DataIn (Len1In,j) HORINT1A.223
DataExt (Len1In+2,j) = DataIn (Len1In,j) HORINT1A.224
DataExt (Len1In+3,j) = DataIn (Len1In,j) HORINT1A.225
END DO HORINT1A.226
HORINT1A.227
END IF HORINT1A.228
HORINT1A.229
DO j = 1, Len2In HORINT1A.230
DO i = 1, Len1In HORINT1A.231
DataExt (i,j) = DataIn (i,j) HORINT1A.232
END DO HORINT1A.233
END DO HORINT1A.234
HORINT1A.235
! 2.2 Extend in j direction. HORINT1A.236
HORINT1A.237
IF (LGlobal) THEN ! 2.2.1 Global. HORINT1A.238
SwapPole = 1 VSB2F405.5
IF (LVector)SwapPole = -1 VSB2F405.6
HORINT1A.239
IF (.not.FirstOnPole) THEN HORINT1A.240
! 2.2.1.1.3 v points. HORINT1A.241
DO i = -2, HalfLen1In+2 HORINT1A.242
DataExt (i,-2) =SwapPole*DataExt(i+HalfLen1In,3) VSB2F405.7
DataExt (i,-1) =SwapPole*DataExt(i+HalfLen1In,2) VSB2F405.8
DataExt (i,0) =SwapPole*DataExt(i+HalfLen1In,1) VSB2F405.9
DataExt (i,Len2In+1) =SwapPole*DataExt(i+HalfLen1In,Len2In) VSB2F405.10
DataExt (i,Len2In+2) =SwapPole*DataExt(i+HalfLen1In,Len2In-1) VSB2F405.11
DataExt (i,Len2In+3) =SwapPole*DataExt(i+HalfLen1In,Len2In-2) VSB2F405.12
END DO HORINT1A.249
HORINT1A.250
DO i = HalfLen1In+3, Len1In+3 HORINT1A.251
DataExt (i,-2) =SwapPole*DataExt(i-HalfLen1In,3) VSB2F405.13
DataExt (i,-1) =SwapPole*DataExt(i-HalfLen1In,2) VSB2F405.14
DataExt (i,0) =SwapPole*DataExt(i-HalfLen1In,1) VSB2F405.15
DataExt (i,Len2In+1) =SwapPole*DataExt(i-HalfLen1In,Len2In) VSB2F405.16
DataExt (i,Len2In+2) =SwapPole*DataExt(i-HalfLen1In,Len2In-1) VSB2F405.17
DataExt (i,Len2In+3) =SwapPole*DataExt(i-HalfLen1In,Len2In-2) VSB2F405.18
END DO HORINT1A.258
HORINT1A.259
ELSE HORINT1A.260
! 2.2.1.1.4 Other points. HORINT1A.261
HORINT1A.262
DO i = -2, HalfLen1In+2 HORINT1A.263
DataExt (i,-2) =SwapPole*DataExt(i+HalfLen1In,4) VSB2F405.19
DataExt (i,-1) =SwapPole*DataExt(i+HalfLen1In,3) VSB2F405.20
DataExt (i,0) =SwapPole*DataExt(i+HalfLen1In,2) VSB2F405.21
DataExt (i,Len2In+1) =SwapPole*DataExt(i+HalfLen1In,Len2In-1) VSB2F405.22
DataExt (i,Len2In+2) =SwapPole*DataExt(i+HalfLen1In,Len2In-2) VSB2F405.23
DataExt (i,Len2In+3) =SwapPole*DataExt(i+HalfLen1In,Len2In-3) VSB2F405.24
END DO HORINT1A.270
HORINT1A.271
DO i = HalfLen1In+3, Len1In+3 HORINT1A.272
DataExt (i,-2) =SwapPole*DataExt(i-HalfLen1In,4) VSB2F405.25
DataExt (i,-1) =SwapPole*DataExt(i-HalfLen1In,3) VSB2F405.26
DataExt (i,0) =SwapPole*DataExt(i-HalfLen1In,2) VSB2F405.27
DataExt (i,Len2In+1) =SwapPole*DataExt(i-HalfLen1In,Len2In-1) VSB2F405.28
DataExt (i,Len2In+2) =SwapPole*DataExt(i-HalfLen1In,Len2In-2) VSB2F405.29
DataExt (i,Len2In+3) =SwapPole*DataExt(i-HalfLen1In,Len2In-3) VSB2F405.30
END DO HORINT1A.279
HORINT1A.280
END IF ! END IF .not.FirstOnPole HORINT1A.281
HORINT1A.282
ELSE HORINT1A.283
! 2.2.2 Limited Area code. HORINT1A.284
HORINT1A.285
DO i = -2, Len1In+3 HORINT1A.286
DataExt (i,-2) = DataExt (i,1) HORINT1A.287
DataExt (i,-1) = DataExt (i,1) HORINT1A.288
DataExt (i,0) = DataExt (i,1) HORINT1A.289
DataExt (i,Len2In+1) = DataExt (i,Len2In) HORINT1A.290
DataExt (i,Len2In+2) = DataExt (i,Len2In) HORINT1A.291
DataExt (i,Len2In+3) = DataExt (i,Len2In) HORINT1A.292
END DO HORINT1A.293
HORINT1A.294
END IF ! END IF LGlobal HORINT1A.295
HORINT1A.296
!----------------------------------------------------------------------- HORINT1A.297
! 3.0 For each output point find i and j so that the point on the HORINT1A.298
! output grid lies between i and i+1 and j and j+1. HORINT1A.299
! and compute the weights HORINT1A.300
!----------------------------------------------------------------------- HORINT1A.301
HORINT1A.302
! 3.1 Get output coordinates and weights where LAM straddling meridian VSB1F401.48
VSB1F401.49
IF (.NOT.LGlobal. AND. VSB1F401.50
& LambdaOrigin+len1In*DLambdaIn.GT.TwoPi) THEN VSB1F401.51
VSB1F401.52
DO j = 1, Len2Out HORINT1A.303
DO i = 1, Len1Out HORINT1A.304
WtLambda(i,j) = (LambdaOut(i,j) - LambdaOrigin) VSB1F401.53
& * RecipDLambdaIn VSB1F401.54
VSB1F401.55
IF (LambdaOut(i,j).LT.LambdaOrigin) VSB1F401.56
& WtLambda(i,j) = WtLambda(i,j) + (TwoPi * RecipDLambdaIn) VSB1F401.57
VSB1F401.58
IOut(i,j) = INT(WtLambda(i,j) + 1) HORINT1A.306
WtLambda(i,j) = WtLambda(i,j) - (IOut(i,j) - 1.0) HORINT1A.307
HORINT1A.308
WtPhi(i,j) = (PhiOut(i,j) - PhiOrigin) * RecipDPhiIn HORINT1A.309
JOut(i,j) = INT(WtPhi(i,j) + 1) HORINT1A.310
WtPhi(i,j) = WtPhi(i,j) - (JOut(i,j) - 1.0) HORINT1A.311
HORINT1A.312
END DO HORINT1A.313
END DO HORINT1A.314
VSB1F401.59
! 3.2 Get output coordinates and weights for normal case VSB1F401.60
VSB1F401.61
ELSE VSB1F401.62
DO j = 1, Len2Out VSB1F401.63
DO i = 1, Len1Out VSB1F401.64
WtLambda(i,j) = (LambdaOut(i,j) - LambdaOrigin) VSB1F401.65
& * RecipDLambdaIn VSB1F401.66
IOut(i,j) = INT(WtLambda(i,j) + 1) VSB1F401.67
WtLambda(i,j) = WtLambda(i,j) - (IOut(i,j) - 1.0) VSB1F401.68
VSB1F401.69
WtPhi(i,j) = (PhiOut(i,j) - PhiOrigin) * RecipDPhiIn VSB1F401.70
JOut(i,j) = INT(WtPhi(i,j) + 1) VSB1F401.71
WtPhi(i,j) = WtPhi(i,j) - (JOut(i,j) - 1.0) VSB1F401.72
VSB1F401.73
END DO VSB1F401.74
END DO VSB1F401.75
END IF VSB1F401.76
HORINT1A.315
!----------------------------------------------------------------------- HORINT1A.316
! 4.0 Perform required interpolations. HORINT1A.317
!----------------------------------------------------------------------- HORINT1A.318
HORINT1A.319
! 4.1 Call the specified interpolation scheme: HORINT1A.320
IF(VIorder.eq.1)Then HORINT1A.321
CALL HorizontalInterpLinear (LowerBound, HORINT1A.322
& Len1In, Len2In, HORINT1A.323
& Len1Out, Len2Out, HORINT1A.324
& DataExt, HORINT1A.325
& WtLambda, HORINT1A.326
& WtPhi, HORINT1A.327
& IOut, HORINT1A.328
& JOut, HORINT1A.329
& DataOut) HORINT1A.330
HORINT1A.331
ELSEIF(VIorder.eq.3)Then VSB2F405.31
CALL HorizontalInterpCubic (LowerBound, HORINT1A.333
& Len1In, Len2In, HORINT1A.334
& Len1Out, Len2Out, HORINT1A.335
& DataExt, HORINT1A.336
& WtLambda, HORINT1A.337
& WtPhi, HORINT1A.338
& IOut, HORINT1A.339
& JOut, HORINT1A.340
& DataOut) HORINT1A.341
HORINT1A.342
HORINT1A.343
ELSEIF(VIorder.eq.5)Then VSB2F405.32
CALL HorizontalInterpQuintic (LowerBound, HORINT1A.345
& Len1In, Len2In, HORINT1A.346
& Len1Out, Len2Out, HORINT1A.347
& DataExt, HORINT1A.348
& WtLambda, HORINT1A.349
& WtPhi, HORINT1A.350
& IOut, HORINT1A.351
& JOut, HORINT1A.352
& DataOut) HORINT1A.353
EndIF HORINT1A.354
HORINT1A.355
!----------------------------------------------------------------------- HORINT1A.356
! 5.0 Perform montonicity enforcement if high order scheme used and HORINT1A.357
! monotonicity required. HORINT1A.358
! section commented out, montonicity code not in vertical interpolation HORINT1A.359
! and as such is currently disabled here (so adjoint not calculated) HORINT1A.360
!----------------------------------------------------------------------- HORINT1A.361
HORINT1A.362
IF(L_mono.AND.VIorder.gt.1)THEN HORINT1A.363
CALL HorizontalInterpLinear (LowerBound, HORINT1A.364
& Len1In, Len2In, HORINT1A.365
& Len1Out, Len2Out, HORINT1A.366
& DataExt, HORINT1A.367
& WtLambda, HORINT1A.368
& WtPhi, HORINT1A.369
& IOut, HORINT1A.370
& JOut, HORINT1A.371
& DataOut) HORINT1A.372
HORINT1A.373
HORINT1A.374
DO j = 1, Len2Out HORINT1A.375
DO i = 1, Len1Out HORINT1A.376
DataHigh (i,j) = DataOut (i,j) HORINT1A.377
END DO HORINT1A.378
END DO HORINT1A.379
HORINT1A.380
CALL HorizontalInterpMonotone HORINT1A.381
& (LowerBound, HORINT1A.382
& Len1In, Len2In, HORINT1A.383
& Len1Out, Len2Out, HORINT1A.384
& DataExt, HORINT1A.385
& DataHigh, HORINT1A.386
& DataMono, HORINT1A.387
& IOut, HORINT1A.388
& JOut, HORINT1A.389
& DataOut) HORINT1A.390
HORINT1A.391
END IF HORINT1A.392
HORINT1A.393
HORINT1A.394
! End of routine. HORINT1A.395
999 CONTINUE HORINT1A.396
RETURN HORINT1A.397
END HORINT1A.398
*ENDIF HORINT1A.399