*IF DEF,C92_1A,OR,DEF,MAKEBC UIE3F404.20
C ******************************COPYRIGHT****************************** GTS2F400.3853
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.3854
C GTS2F400.3855
C Use, duplication or disclosure of this code is subject to the GTS2F400.3856
C restrictions as set forth in the contract. GTS2F400.3857
C GTS2F400.3858
C Meteorological Office GTS2F400.3859
C London Road GTS2F400.3860
C BRACKNELL GTS2F400.3861
C Berkshire UK GTS2F400.3862
C RG12 2SZ GTS2F400.3863
C GTS2F400.3864
C If no contract has been raised with this copy of the code, the use, GTS2F400.3865
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.3866
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.3867
C Modelling at the above address. GTS2F400.3868
C ******************************COPYRIGHT****************************** GTS2F400.3869
C GTS2F400.3870
CLL SUBROUTINE H_INT_CO---------------------------------------------- HINTCO1A.3
CLL HINTCO1A.4
CLL Purpose: Calculates bi-linear horizontal interpolation HINTCO1A.5
CLL coefficients and gather indices for interpolating HINTCO1A.6
CLL between generalised latitude-longitude grids (eg HINTCO1A.7
CLL global, regional or rotated lat-lon grid) in which the HINTCO1A.8
CLL gridlength may vary with latitude and/or longitude. The HINTCO1A.9
CLL interpolation is carried out by subroutine HINTCO1A.10
CLL H_INT. Gather indices point to bottom left hand HINTCO1A.11
CLL corner and bottom right hand corner of each grid box on HINTCO1A.12
CLL source grid enclosing a target point. Two indices are HINTCO1A.13
CLL needed to cater for east-west (lambda direction) cyclic HINTCO1A.14
CLL boundaries when the source data is global. If a target po HINTCO1A.15
CLL falls outside the domain of the source data, one sided HINTCO1A.16
CLL differencing is used. The source latitude coordinates HINTCO1A.17
CLL must be supplied in decreasing order. The source long- HINTCO1A.18
CLL itude coordinates must be supplied in increasing order, HINTCO1A.19
CLL starting at any value, but not wrapping round. The HINTCO1A.20
CLL target points may be specified in any order. HINTCO1A.21
CLL HINTCO1A.22
CLL A.Dickinson <- programmer of some or all of previous code or changes HINTCO1A.23
CLL J.Gregory <- programmer of some or all of previous code or changes HINTCO1A.24
CLL HINTCO1A.25
CLL Model Modification history from model version 3.0: HINTCO1A.26
CLL version Date HINTCO1A.27
CLL 3.4 24/06/94 New checks to keep within array bounds for ADR1F304.96
CLL non_cyclic cases. D. Robinson ADR1F304.97
! 4.5 29/07/98 Optimisation changes for T3E UDG5F405.56
! Author D.M. Goddard UDG5F405.57
CLL HINTCO1A.28
CLL Programming standard: HINTCO1A.29
CLL Unified Model Documentation Paper No 3 HINTCO1A.30
CLL Version No 1 15/1/90 HINTCO1A.31
CLL HINTCO1A.32
CLL System component: S121,S122 HINTCO1A.33
CLL HINTCO1A.34
CLL System task: S1 HINTCO1A.35
CLL HINTCO1A.36
CLL Documentation: HINTCO1A.37
CLL The interpolation formulae are described in HINTCO1A.38
CLL unified model on-line documentation paper S1. HINTCO1A.39
CLL HINTCO1A.40
CLL ----------------------------------------------------------------- HINTCO1A.41
C HINTCO1A.42
C*L Arguments:------------------------------------------------------- HINTCO1A.43
SUBROUTINE H_INT_CO 16HINTCO1A.44
*(INDEX_B_L,INDEX_B_R,WEIGHT_T_R,WEIGHT_B_R,WEIGHT_T_L,WEIGHT_B_L HINTCO1A.45
*,LAMBDA_SRCE,PHI_SRCE,LAMBDA_TARG,PHI_TARG HINTCO1A.46
*,POINTS_LAMBDA_SRCE,POINTS_PHI_SRCE,POINTS,CYCLIC) HINTCO1A.47
HINTCO1A.48
IMPLICIT NONE HINTCO1A.49
HINTCO1A.50
INTEGER HINTCO1A.51
* POINTS_LAMBDA_SRCE !IN Number of lambda points on source grid HINTCO1A.52
*,POINTS_PHI_SRCE !IN Number of phi points on source grid HINTCO1A.53
*,POINTS !IN Total number of points on target grid HINTCO1A.54
*,INDEX_B_L(POINTS) !OUT Index of bottom lefthand corner HINTCO1A.55
* ! of source gridbox HINTCO1A.56
*,INDEX_B_R(POINTS) !OUT Index of bottom righthand corner HINTCO1A.57
* ! of source gridbox HINTCO1A.58
HINTCO1A.59
REAL HINTCO1A.60
* LAMBDA_TARG(POINTS) !IN Lambda coords of target grid in degrees HINTCO1A.61
* ! using same rotation as source grid HINTCO1A.62
*,PHI_TARG(POINTS) !IN Phi coords of target grid in degrees HINTCO1A.63
* ! using same rotation as source grid HINTCO1A.64
*,WEIGHT_T_R(POINTS) !OUT Weight applied to value at top right HINTCO1A.65
* ! hand corner of source gridbox HINTCO1A.66
*,WEIGHT_B_L(POINTS) !OUT Weight applied to value at bottom left HINTCO1A.67
* ! hand corner of source gridbox HINTCO1A.68
*,WEIGHT_B_R(POINTS) !OUT Weight applied to value at bottom right HINTCO1A.69
* ! hand corner of source gridbox HINTCO1A.70
*,WEIGHT_T_L(POINTS) !OUT Weight applied to value at top left HINTCO1A.71
* ! hand corner of source gridbox HINTCO1A.72
*,LAMBDA_SRCE(POINTS_LAMBDA_SRCE) !IN Lambda coords of source grid HINTCO1A.73
* ! in degrees HINTCO1A.74
*,PHI_SRCE(POINTS_PHI_SRCE) !IN Phi coords of target grid in degree HINTCO1A.75
HINTCO1A.76
LOGICAL HINTCO1A.77
* CYCLIC !IN =T, then source data is cyclic HINTCO1A.78
* ! =F, then source data is non-cyclic HINTCO1A.79
HINTCO1A.80
C Local arrays:--------------------------------------------------------- HINTCO1A.81
REAL HINTCO1A.82
* T_LAMBDA(POINTS) !Local value of target longitude HINTCO1A.83
HINTCO1A.84
INTEGER HINTCO1A.85
* IXP1(POINTS) !Longitudinal index plus 1 HINTCO1A.86
*,IX(POINTS) !Longitudinal index HINTCO1A.87
*,IY(POINTS) !Latitudinal index HINTCO1A.88
C External subroutines called:------------------------------------------ HINTCO1A.89
C None HINTCO1A.90
C*---------------------------------------------------------------------- HINTCO1A.91
C*L Local variables:--------------------------------------------------- HINTCO1A.92
REAL HINTCO1A.93
* A !Longitudinal weight HINTCO1A.94
*,B !Latitudinal weight HINTCO1A.95
HINTCO1A.96
INTEGER HINTCO1A.97
* I,J !Loop indices HINTCO1A.98
C ---------------------------------------------------------------------- HINTCO1A.99
HINTCO1A.100
CL 1. Initialise arrays HINTCO1A.101
HINTCO1A.102
IF(CYCLIC)THEN HINTCO1A.103
HINTCO1A.104
C 1.1 Cyclic case HINTCO1A.105
HINTCO1A.106
DO I=1,POINTS HINTCO1A.107
HINTCO1A.108
C Scale target longitude so that it falls between LAMBDA_SRCE(1) HINTCO1A.109
C and LAMBDA_SRCE(1)+360 HINTCO1A.110
T_LAMBDA(I)=MOD(LAMBDA_TARG(I)-LAMBDA_SRCE(1)+720.,360.) HINTCO1A.111
& +LAMBDA_SRCE(1) HINTCO1A.112
HINTCO1A.113
C Initialise longitudinal & latitudinal indices HINTCO1A.114
IX(I)=0 HINTCO1A.115
IY(I)=1 HINTCO1A.116
ENDDO HINTCO1A.117
HINTCO1A.118
ELSE HINTCO1A.119
HINTCO1A.120
C 1.2 Non cyclic case HINTCO1A.121
HINTCO1A.122
DO I=1,POINTS HINTCO1A.123
HINTCO1A.124
C Assign target longitude to local array for use below HINTCO1A.125
T_LAMBDA(I)=LAMBDA_TARG(I) HINTCO1A.126
HINTCO1A.127
C Initialise longitudinal & latitudinal indices HINTCO1A.128
IX(I)=1 HINTCO1A.129
IY(I)=1 HINTCO1A.130
HINTCO1A.131
ENDDO HINTCO1A.132
HINTCO1A.133
ENDIF HINTCO1A.134
HINTCO1A.135
CL 2. Calculate lat and lon index of bottom left hand corner of HINTCO1A.136
CL source grid box enclosing each target point. HINTCO1A.137
HINTCO1A.138
C Longitude HINTCO1A.139
DO I=1,POINTS UDG5F405.58
DO J=1,POINTS_LAMBDA_SRCE UDG5F405.59
IF(LAMBDA_SRCE(J).LE.T_LAMBDA(I))THEN UDG5F405.60
IX(I)=J UDG5F405.61
ELSE UDG5F405.62
GOTO 120 UDG5F405.63
ENDIF UDG5F405.64
ENDDO UDG5F405.65
120 CONTINUE UDG5F405.66
ENDDO UDG5F405.67
HINTCO1A.145
C Latitude HINTCO1A.146
DO I=1,POINTS UDG5F405.68
DO J=1,POINTS_PHI_SRCE UDG5F405.69
IF(PHI_SRCE(J).GE.PHI_TARG(I))THEN UDG5F405.70
IY(I)=J UDG5F405.71
ELSE UDG5F405.72
GOTO 140 UDG5F405.73
ENDIF UDG5F405.74
ENDDO UDG5F405.75
140 CONTINUE UDG5F405.76
ENDDO UDG5F405.77
UDG5F405.78
UDG5F405.79
HINTCO1A.152
CL 3. Correct 1-D indices for wrap around etc and then calculate HINTCO1A.153
CL 2-D indices of bottom left and bottom right hand corner HINTCO1A.154
CL of each grid box. HINTCO1A.155
HINTCO1A.156
IF(CYCLIC)THEN HINTCO1A.157
C 3.1 Cyclic case HINTCO1A.158
HINTCO1A.159
DO I=1,POINTS HINTCO1A.160
HINTCO1A.161
C Set index for cyclic wrap around HINTCO1A.162
IF(IX(I).LT.1)THEN HINTCO1A.163
IX(I)=POINTS_LAMBDA_SRCE HINTCO1A.164
T_LAMBDA(I)=T_LAMBDA(I)+360. HINTCO1A.165
ENDIF HINTCO1A.166
HINTCO1A.167
C Set index for one sided difference if target point to north or HINTCO1A.168
C south of source area. HINTCO1A.169
IY(I)=MAX(IY(I),1) HINTCO1A.170
IY(I)=MIN(IY(I),POINTS_PHI_SRCE-1) HINTCO1A.171
HINTCO1A.172
C 2-D indices HINTCO1A.173
INDEX_B_L(I)=IX(I)+IY(I)*POINTS_LAMBDA_SRCE HINTCO1A.174
INDEX_B_R(I)=INDEX_B_L(I)+1 HINTCO1A.175
HINTCO1A.176
C Correct for cyclic boundaries if target point outside source grid. HINTCO1A.177
HINTCO1A.178
IXP1(I)=IX(I)+1 HINTCO1A.179
IF(IX(I).EQ.POINTS_LAMBDA_SRCE)THEN HINTCO1A.180
INDEX_B_R(I)=INDEX_B_R(I)-POINTS_LAMBDA_SRCE HINTCO1A.181
IXP1(I)=IXP1(I)-POINTS_LAMBDA_SRCE HINTCO1A.182
ENDIF HINTCO1A.183
HINTCO1A.184
ENDDO HINTCO1A.185
HINTCO1A.186
ELSE HINTCO1A.187
HINTCO1A.188
C 3.2 Non cyclic case HINTCO1A.189
HINTCO1A.190
DO I=1,POINTS HINTCO1A.191
HINTCO1A.192
C Set index for one sided difference if outside source area HINTCO1A.193
IX(I)=MAX(IX(I),1) HINTCO1A.194
IX(I)=MIN(IX(I),POINTS_LAMBDA_SRCE-1) HINTCO1A.195
IF (IX(I).LT.1) THEN ! IX(I) < 1 if POINTS_LAMBDA_SRCE = 1 ADR1F304.98
IX(I)=1 ADR1F304.99
ENDIF ADR1F304.100
ADR1F304.101
IXP1(I)=IX(I)+1 HINTCO1A.196
IXP1(I)=MIN(IXP1(I),POINTS_LAMBDA_SRCE) ADR1F304.102
HINTCO1A.197
C Set index for one sided difference if outside source area HINTCO1A.198
IY(I)=MAX(IY(I),1) HINTCO1A.199
IY(I)=MIN(IY(I),POINTS_PHI_SRCE-1) HINTCO1A.200
IF (IY(I).LT.1) THEN ! IY(I) < 1 if POINTS_PHI_SRCE = 1 ADR1F304.103
IY(I)=1 ADR1F304.104
ENDIF ADR1F304.105
ADR1F304.106
HINTCO1A.201
C 2-D indices HINTCO1A.202
INDEX_B_L(I)=IX(I)+IY(I)*POINTS_LAMBDA_SRCE HINTCO1A.203
INDEX_B_R(I)=INDEX_B_L(I)+1 HINTCO1A.204
HINTCO1A.205
ENDDO HINTCO1A.206
HINTCO1A.207
ENDIF HINTCO1A.208
HINTCO1A.209
CL 4. Compute interpolation weights HINTCO1A.210
HINTCO1A.211
DO I=1,POINTS HINTCO1A.212
HINTCO1A.213
C Calculate basic weights (equation 2.2) HINTCO1A.214
A=(AMOD(360.+LAMBDA_SRCE(IXP1(I))-LAMBDA_SRCE(IX(I)),360.)) HINTCO1A.215
IF(A.NE.0.)THEN HINTCO1A.216
A=(T_LAMBDA(I)-LAMBDA_SRCE(IX(I)))/A HINTCO1A.217
ELSE HINTCO1A.218
A=0. HINTCO1A.219
ENDIF HINTCO1A.220
HINTCO1A.221
B=(PHI_TARG(I)-PHI_SRCE(IY(I)+1))/ HINTCO1A.222
* (PHI_SRCE(IY(I))-PHI_SRCE(IY(I)+1)) HINTCO1A.223
HINTCO1A.224
C Calculate bi-linear interpolation weights as per equation (2.1) HINTCO1A.225
HINTCO1A.226
WEIGHT_T_R(I)=A*B HINTCO1A.227
WEIGHT_B_L(I)=(1.-A)*(1.-B) HINTCO1A.228
WEIGHT_T_L(I)=(1.-A)*B HINTCO1A.229
WEIGHT_B_R(I)=A*(1.-B) HINTCO1A.230
HINTCO1A.231
ENDDO HINTCO1A.232
HINTCO1A.233
RETURN HINTCO1A.234
END HINTCO1A.235
*ENDIF HINTCO1A.236