*IF DEF,C92_1A VINTTH1A.2
C ******************************COPYRIGHT****************************** GTS2F400.11665
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.11666
C GTS2F400.11667
C Use, duplication or disclosure of this code is subject to the GTS2F400.11668
C restrictions as set forth in the contract. GTS2F400.11669
C GTS2F400.11670
C Meteorological Office GTS2F400.11671
C London Road GTS2F400.11672
C BRACKNELL GTS2F400.11673
C Berkshire UK GTS2F400.11674
C RG12 2SZ GTS2F400.11675
C GTS2F400.11676
C If no contract has been raised with this copy of the code, the use, GTS2F400.11677
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.11678
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.11679
C Modelling at the above address. GTS2F400.11680
C ******************************COPYRIGHT****************************** GTS2F400.11681
C GTS2F400.11682
CLL SUBROUTINE V_INT_TH---------------------------------------------- VINTTH1A.3
CLL VINTTH1A.4
CLL Purpose: Calculates potential temperature of model layers using VINTTH1A.5
CLL temperatures input on an arbitray set of pressure levels VINTTH1A.6
CLL VINTTH1A.7
CLL R.Swinbank <- programmer of some or all of previous code or changes VINTTH1A.8
CLL D.Robinson <- programmer of some or all of previous code or changes VINTTH1A.9
CLL VINTTH1A.10
CLL Model Modification history from model version 3.0: VINTTH1A.11
CLL version Date VINTTH1A.12
CLL VINTTH1A.13
CLL 4.2 01/08/96 *DEF CRAY removed. Use of 32-bit log and exp GSS5F402.200
CLL for CRAY shared memory computers no longer GSS5F402.201
CLL supported. GSS5F402.202
CLL GSS5F402.203
CLL Author: A. Dickinson Reviewer: R. Rawlins GSS5F402.204
CLL GSS5F402.205
CLL Programming standard : VINTTH1A.14
CLL VINTTH1A.15
CLL Logical components covered : S113 VINTTH1A.16
CLL VINTTH1A.17
CLL Project task : VINTTH1A.18
CLL VINTTH1A.19
CLL Documentation: The interpolation formulae are described in VINTTH1A.20
CLL unified model on-line documentation paper S1. VINTTH1A.21
CLL VINTTH1A.22
CLLEND ----------------------------------------------------------------- VINTTH1A.23
C VINTTH1A.24
C*L Arguments:------------------------------------------------------- VINTTH1A.25
SUBROUTINE V_INT_TH VINTTH1A.26
*(P_HALF,P_EXNER_HALF,THETA,T_SRCE,POINTS,P_SRCE,P_LEVELS, VINTTH1A.27
* SRCE_LEVELS,PSTAR,AKH,BKH) VINTTH1A.28
VINTTH1A.29
IMPLICIT NONE VINTTH1A.30
VINTTH1A.31
INTEGER VINTTH1A.32
* POINTS !IN Number of points to be processed VINTTH1A.33
*,P_LEVELS !IN Number of model levels VINTTH1A.34
*,SRCE_LEVELS !IN Number of input levels VINTTH1A.35
VINTTH1A.36
REAL VINTTH1A.37
* P_HALF(POINTS,P_LEVELS+1) !IN Pressure of model half levels VINTTH1A.38
*,P_EXNER_HALF(POINTS,P_LEVELS+1) !IN Exner pressure at model VINTTH1A.39
* ! half levels VINTTH1A.40
*,PSTAR(POINTS) !IN surface pressure VINTTH1A.41
*,AKH(P_LEVELS+1) !IN Hybrid coord. A at half levels VINTTH1A.42
*,BKH(P_LEVELS+1) !IN Hybrid coord. B at half levels VINTTH1A.43
*,THETA(POINTS,P_LEVELS) !OUT Potential temperature at full levels VINTTH1A.44
*,T_SRCE(POINTS,SRCE_LEVELS) !IN Input temperature fields VINTTH1A.45
*,P_SRCE(POINTS,SRCE_LEVELS) !IN Pressure of input temperature VINTTH1A.46
* ! fields VINTTH1A.47
VINTTH1A.48
C Workspace usage:----------------------------------------------------- VINTTH1A.49
REAL TL(POINTS),PL(POINTS) VINTTH1A.50
INTEGER KI(POINTS) VINTTH1A.51
C External subroutines called:----------------------------------------- VINTTH1A.52
C None GSS5F402.206
C*--------------------------------------------------------------------- VINTTH1A.54
C Define local variables:---------------------------------------------- VINTTH1A.55
INTEGER I,J,K,JK,JI VINTTH1A.56
REAL PUB ! Pressure at top of layer VINTTH1A.57
REAL PLB ! Pressure at bottom of layer VINTTH1A.58
REAL P_EXNER_FULL ! Exner pressure at full model levels VINTTH1A.59
REAL PKPH,TKPH,THICK,PJ,PJM1,TJ,TJM1 VINTTH1A.60
C---------------------------------------------------------------------- VINTTH1A.64
C Constants from comdecks:--------------------------------------------- VINTTH1A.65
*CALL C_G 
VINTTH1A.66
*CALL C_R_CP VINTTH1A.67
*CALL C_LAPSE VINTTH1A.68
REAL LAPSE_R_OVER_G VINTTH1A.69
PARAMETER(LAPSE_R_OVER_G=LAPSE*R/G) VINTTH1A.70
C---------------------------------------------------------------------- VINTTH1A.71
VINTTH1A.72
*CALL P_EXNERC VINTTH1A.73
VINTTH1A.74
CL 1. Initialise arrays PL,TL,KI and THETA VINTTH1A.75
VINTTH1A.76
DO I=1,POINTS VINTTH1A.77
KI(I)=1 VINTTH1A.78
TL(I)=0.0 VINTTH1A.79
PL(I)=P_HALF(I,1) VINTTH1A.80
ENDDO VINTTH1A.81
VINTTH1A.82
C THETA is used to accumulate model layer thicknesses (x2) in section 2 VINTTH1A.83
C The values are converted to potential temperatures in section 3. VINTTH1A.84
DO J=1,P_LEVELS VINTTH1A.85
DO I=1,POINTS VINTTH1A.86
THETA(I,J)=0.0 VINTTH1A.87
ENDDO VINTTH1A.88
ENDDO VINTTH1A.89
VINTTH1A.90
CL 2. Loop over JK. Each time around we either increment one input VINTTH1A.91
CL level (J) or one model level (K), as appropriate; JK=J+K. VINTTH1A.92
CL For each combination j,k we calculate the thickness of the VINTTH1A.93
CL overlap between the interval P_SRCE(j-1 to j) and model layer VINTTH1A.94
CL k (i.e. P_HALF(k to k+1)). VINTTH1A.95
VINTTH1A.96
DO 100 JK=2,P_LEVELS+SRCE_LEVELS+1 VINTTH1A.97
VINTTH1A.98
CDIR$ IVDEP VINTTH1A.99
! Fujitsu vectorization directive GRB0F405.557
!OCL NOVREC GRB0F405.558
DO I=1,POINTS VINTTH1A.100
VINTTH1A.101
C Get J (input level) VINTTH1A.102
C Max. value of J is SRCE_LEVELS+1 (above top input level). VINTTH1A.103
C JI is J value for interpolation - limited to SRCE_LEVELS, since VINTTH1A.104
C we extrapolate above top level. VINTTH1A.105
VINTTH1A.106
K=KI(I) VINTTH1A.107
J=MIN(JK-K,SRCE_LEVELS+1) VINTTH1A.108
JI=MIN(J,SRCE_LEVELS) VINTTH1A.109
VINTTH1A.110
C Gather P,T values for interpolation VINTTH1A.111
PJ=P_SRCE(I,JI) VINTTH1A.112
TJ=T_SRCE(I,JI) VINTTH1A.113
PJM1=P_SRCE(I,JI-1) VINTTH1A.114
TJM1=T_SRCE(I,J-1) VINTTH1A.115
VINTTH1A.116
C If K>P_LEVELS integration is complete; VINTTH1A.117
C if input pressure is greater than range of model levels integration is VINTTH1A.118
C not yet started. VINTTH1A.119
C In both cases no action is required. VINTTH1A.120
IF (K.LE.P_LEVELS .AND. PJ.LT.P_HALF(I,1)) THEN VINTTH1A.121
VINTTH1A.122
C If TL is not set, integration has just started, so we require T at VINTTH1A.123
C model level 1/2 (P_HALF(1)), which must be in current input pressure VINTTH1A.124
C interval. VINTTH1A.125
VINTTH1A.126
IF(TL(I).EQ.0.0) THEN VINTTH1A.127
IF(J.EQ.1) THEN VINTTH1A.128
C If J=1, extrapolate below input level 1 VINTTH1A.129
TL(I)=TJ*(P_HALF(I,1)/PJ)**LAPSE_R_OVER_G GSS5F402.207
ELSE VINTTH1A.135
C Otherwise interpolate between levels J-1 and J VINTTH1A.136
TL(I)=TJM1+ALOG(P_HALF(I,1)/PJM1)*(TJ-TJM1) VINTTH1A.141
* /ALOG(PJ/PJM1) VINTTH1A.142
END IF VINTTH1A.144
END IF VINTTH1A.145
VINTTH1A.146
PKPH=P_HALF(I,K+1) VINTTH1A.147
VINTTH1A.148
C Check whether next level up (from PL) is a model layer boundary, or VINTTH1A.149
C an input level; then integrate from PL to this level. VINTTH1A.150
VINTTH1A.151
IF(PJ.GE.PKPH.AND.J.LE.SRCE_LEVELS) THEN VINTTH1A.152
C Integrate from PL to input level J VINTTH1A.153
GSS5F402.208
THICK=(TL(I)+TJ)*ALOG(PL(I)/PJ) VINTTH1A.157
GSS5F402.209
C Save P & T values at current level VINTTH1A.159
PL(I)=PJ VINTTH1A.160
TL(I)=TJ VINTTH1A.161
VINTTH1A.162
ELSE VINTTH1A.163
C Get T at model level K+1/2 VINTTH1A.164
IF(J.EQ.1) THEN VINTTH1A.165
C If J=1, extrapolate below input level 1 VINTTH1A.166
GSS5F402.210
TKPH=TJ*(PKPH/PJ)**LAPSE_R_OVER_G VINTTH1A.170
GSS5F402.211
ELSE VINTTH1A.172
C Otherwise interpolate between levels J-1 and J VINTTH1A.173
GSS5F402.212
TKPH=TJM1+ALOG(PKPH/PJM1)*(TJ-TJM1)/ALOG(PJ/PJM1) VINTTH1A.177
GSS5F402.213
END IF VINTTH1A.179
VINTTH1A.180
C Integrate from PL to model level k+1/2 VINTTH1A.181
GSS5F402.214
THICK=(TL(I)+TKPH)*ALOG(PL(I)/PKPH) GSS5F402.215
C Save P & T values and increment K, so next model level is used VINTTH1A.187
C on next iteration of JK VINTTH1A.188
PL(I)=PKPH VINTTH1A.189
TL(I)=TKPH VINTTH1A.190
KI(I)=K+1 VINTTH1A.191
END IF VINTTH1A.192
VINTTH1A.193
C Integration results are accumulated in THETA; they are converted VINTTH1A.194
C to theta values in section 3. VINTTH1A.195
THETA(I,K)=THETA(I,K)+THICK VINTTH1A.196
VINTTH1A.197
END IF VINTTH1A.198
VINTTH1A.199
ENDDO VINTTH1A.200
100 CONTINUE VINTTH1A.201
VINTTH1A.202
VINTTH1A.203
CL 3. Convert thickness values to potential temperature values. VINTTH1A.204
VINTTH1A.205
DO J=1,P_LEVELS VINTTH1A.206
DO I=1,POINTS VINTTH1A.207
VINTTH1A.208
C Model layer above source data coverage: equation (3.28) VINTTH1A.209
IF(P_HALF(I,J).LE.P_SRCE(I,SRCE_LEVELS))THEN VINTTH1A.210
PUB=PSTAR(I)*BKH(J+1) + AKH(J+1) VINTTH1A.211
PLB=PSTAR(I)*BKH(J) + AKH(J) VINTTH1A.212
THETA(I,J)=T_SRCE(I,SRCE_LEVELS)/ VINTTH1A.213
& P_EXNER_C( P_EXNER_HALF(I,J+1),P_EXNER_HALF(I,J), VINTTH1A.214
& PUB,PLB,KAPPA ) VINTTH1A.215
VINTTH1A.216
ELSE VINTTH1A.217
C Convert layer thickness to theta: finalise equation (3.27) VINTTH1A.218
THETA(I,J)=KAPPA*THETA(I,J) VINTTH1A.219
* /(2.*(P_EXNER_HALF(I,J)-P_EXNER_HALF(I,J+1))) VINTTH1A.220
ENDIF VINTTH1A.221
VINTTH1A.222
C If model layer straddles top srce level and interpolated temperature VINTTH1A.223
C falls outside the range of the top two source level temperatures, VINTTH1A.224
C then use equ 3.28 VINTTH1A.225
VINTTH1A.226
IF((P_HALF(I,J).GT.P_SRCE(I,SRCE_LEVELS)).AND. VINTTH1A.227
* (P_HALF(I,J+1).LE.P_SRCE(I,SRCE_LEVELS)))THEN VINTTH1A.228
VINTTH1A.229
PUB = PSTAR(I)*BKH(J+1) + AKH(J+1) VINTTH1A.230
PLB = PSTAR(I)*BKH(J) + AKH(J) VINTTH1A.231
P_EXNER_FULL = P_EXNER_C VINTTH1A.232
* ( P_EXNER_HALF(I,J+1),P_EXNER_HALF(I,J),PUB,PLB,KAPPA ) VINTTH1A.233
TJ=THETA(I,J)*P_EXNER_FULL VINTTH1A.234
IF((TJ.GT.T_SRCE(I,SRCE_LEVELS).AND. VINTTH1A.235
* TJ.GT.T_SRCE(I,SRCE_LEVELS-1)).OR. VINTTH1A.236
* (TJ.LT.T_SRCE(I,SRCE_LEVELS).AND. VINTTH1A.237
* TJ.LT.T_SRCE(I,SRCE_LEVELS-1)))THEN VINTTH1A.238
THETA(I,J)=T_SRCE(I,SRCE_LEVELS)/P_EXNER_FULL VINTTH1A.239
ENDIF VINTTH1A.240
VINTTH1A.241
ENDIF VINTTH1A.242
ENDDO VINTTH1A.243
ENDDO VINTTH1A.244
VINTTH1A.245
RETURN VINTTH1A.246
END VINTTH1A.247
*ENDIF VINTTH1A.248