*IF DEF,A16_1A CONTRA1A.2
C ******************************COPYRIGHT****************************** GTS2F400.1153
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.1154
C GTS2F400.1155
C Use, duplication or disclosure of this code is subject to the GTS2F400.1156
C restrictions as set forth in the contract. GTS2F400.1157
C GTS2F400.1158
C Meteorological Office GTS2F400.1159
C London Road GTS2F400.1160
C BRACKNELL GTS2F400.1161
C Berkshire UK GTS2F400.1162
C RG12 2SZ GTS2F400.1163
C GTS2F400.1164
C If no contract has been raised with this copy of the code, the use, GTS2F400.1165
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.1166
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.1167
C Modelling at the above address. GTS2F400.1168
C ******************************COPYRIGHT****************************** GTS2F400.1169
C GTS2F400.1170
CLL SUBROUTINE CONTRAIL------------------------------------------------ CONTRA1A.3
CLL CONTRA1A.4
CLL PURPOSE: TO CALCULATE CONTRAIL UPPER AND LOWER HEIGHT CONTRA1A.5
CLL CONTRA1A.6
CLL WRITTEN BY C.M.ROBERTS CONTRA1A.7
CLL MODIFIED BY J.T.HEMING CONTRA1A.8
CLL CONTRA1A.9
CLL MODEL MODIFICATION HISTORY FROM MODEL VERSION 3.0: CONTRA1A.10
CLL VERSION DATE CONTRA1A.11
CLL 4.4 30/07/97 Change limit for TERM6 to prevent calculations GDR2F404.15
CLL going out of T3E real number range. D. Robinson GDR2F404.16
!LL 4.5 20/04/98 Start-end args added to enable dupicate halo GSM1F405.597
!LL calculations to be avoided. S.D.Mullerworth GSM1F405.598
CLL CONTRA1A.12
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4, CONTRA1A.13
CLL VERSION 2, DATED 18/01/90 CONTRA1A.14
CLL CONTRA1A.15
CLL LOGICAL COMPONENT NUMBER: D431 CONTRA1A.16
CLL CONTRA1A.17
CLL PROJECT TASK: CONTRA1A.18
CLL CONTRA1A.19
CLL DOCUMENTATION: CONTRA1A.20
CLL CONTRA1A.21
CLLEND------------------------------------------------------------------ CONTRA1A.22
C CONTRA1A.23
C*L ARGUMENTS:---------------------------------------------------------- CONTRA1A.24
SUBROUTINE CONTRAIL( 2,2CONTRA1A.25
C data in CONTRA1A.26
& P,T,PSTAR, CONTRA1A.27
& P_EXNER_HALF,THETA, CONTRA1A.28
C data out CONTRA1A.29
& UPPER_CONTRAIL,LOWER_CONTRAIL, CONTRA1A.30
C constants in CONTRA1A.31
& POINTS,P_LEVELS, GSM1F405.599
! Range of points to calculate GSM1F405.600
& START,END) GSM1F405.601
C----------------------------------------------------------------------- CONTRA1A.33
IMPLICIT NONE CONTRA1A.34
C----------------------------------------------------------------------- CONTRA1A.35
EXTERNAL ICAO_HT CONTRA1A.36
C CONTRA1A.37
INTEGER CONTRA1A.38
* POINTS ! IN NO OF POINTS CONTRA1A.39
*,P_LEVELS ! IN NO OF MODEL LEVELS CONTRA1A.40
&,START,END ! IN Range of points to calculate GSM1F405.602
C CONTRA1A.41
REAL CONTRA1A.42
* T(POINTS,P_LEVELS) ! IN TEMPERATURE AT FULL LEVELS CONTRA1A.43
*,P(POINTS,P_LEVELS) ! IN PRESSURE " " " CONTRA1A.44
*,PSTAR(POINTS) ! IN SURFACE PRESSURE CONTRA1A.45
*,P_EXNER_HALF(POINTS,P_LEVELS+1) ! IN EXNER PRESSURE AT MODEL CONTRA1A.46
* ! HALF LEVELS CONTRA1A.47
*,THETA(POINTS,P_LEVELS) ! IN POTENTIAL TEMPERATURE AT CONTRA1A.48
* ! MODEL FULL LEVELS CONTRA1A.49
*,UPPER_CONTRAIL(POINTS) ! OUT UPPER VALUE OF THE CONTRAIL CONTRA1A.50
*,LOWER_CONTRAIL(POINTS) ! OUT LOWER VALUE OF THE CONTRAIL CONTRA1A.51
C----------------------------------------------------------------------- CONTRA1A.52
C Local Variables CONTRA1A.53
C----------------------------------------------------------------------- CONTRA1A.54
INTEGER CONTRA1A.55
* I,J ! LOOP COUNTERS CONTRA1A.56
CONTRA1A.57
REAL CONTRA1A.58
* DEL_EXNER_JM1 ! EXNER PRESSURE DIFF BET LEVELS CONTRA1A.59
* ! J-3/2 AND J-1/2 CONTRA1A.60
*,DEL_EXNER_J ! " " " " "J-1/2 AND J+1/2 CONTRA1A.61
*,TERM1 ! \ CONTRA1A.62
*,TERM2 ! } TEMPORARY STORAGE VARIABLES CONTRA1A.63
*,TERM3 ! } CONTRA1A.64
*,TERM4 ! } USED IN SEVERAL DIFFERENT CONTRA1A.65
*,TERM5 ! } CALCULATIONS CONTRA1A.66
*,TERM6 ! / CONTRA1A.67
*,MINTRA_M_14(POINTS,P_LEVELS) !'MINTRA line minus 14 degrees' temp CONTRA1A.68
*,EXPONENT ! lapse rate between two levels *R/g CONTRA1A.69
*,P_M(POINTS) ! Pressure of point of intersection CONTRA1A.70
*,P_MAX(POINTS) ! Maximum pressure of intersection in a column CONTRA1A.71
*,P_MIN(POINTS) ! Minimum pressure of intersection in a column CONTRA1A.72
*,TEST1 ! Test variable CONTRA1A.73
*,TEST2 ! Test variable CONTRA1A.74
*,P_LIMIT ! Upper limit (lowest pressure) at which to CONTRA1A.75
* ! perform calculations CONTRA1A.76
C----------------------------------------------------------------------- CONTRA1A.77
C Logicals CONTRA1A.78
C----------------------------------------------------------------------- CONTRA1A.79
LOGICAL CONTRA1A.80
* INTERSECTION(POINTS) ! TRUE if intersection occurs. CONTRA1A.81
C----------------------------------------------------------------------- CONTRA1A.82
C Constants CONTRA1A.83
C----------------------------------------------------------------------- CONTRA1A.84
C* CONTRA1A.85
*CALL C_R_CP 
CONTRA1A.86
*CALL C_G CONTRA1A.87
C----------------------------------------------------------------------- CONTRA1A.88
REAL CP_OVER_TWO_G CONTRA1A.89
*,A0,A1 ! The 'MINTRA line minus 14 degrees' is CONTRA1A.90
C ! approximated by Tm=A0*(P**A1) CONTRA1A.91
PARAMETER(CP_OVER_TWO_G=CP/(2.*G) CONTRA1A.92
& ,A0=137.14816 ! This constant set for pressure in Pa CONTRA1A.93
& ,A1=0.046822 CONTRA1A.94
& ,P_LIMIT=5000.0) ! No calculations above this level in Pa CONTRA1A.95
C*---------------------------------------------------------------------- CONTRA1A.96
CL Initialise logical array and max and min pressure arrays CONTRA1A.97
C----------------------------------------------------------------------- CONTRA1A.98
DO I=START,END GSM1F405.603
INTERSECTION(I)=.FALSE. CONTRA1A.100
P_MAX(I)=0.0 CONTRA1A.101
P_MIN(I)=PSTAR(I) CONTRA1A.102
ENDDO CONTRA1A.103
C----------------------------------------------------------------------- CONTRA1A.104
CL Calculate the approximation to the 'MINTRA-line minus 14 degrees' CONTRA1A.105
CL at full levels for each point CONTRA1A.106
C----------------------------------------------------------------------- CONTRA1A.107
DO J=1,P_LEVELS CONTRA1A.108
DO I=START,END GSM1F405.604
MINTRA_M_14(I,J)=A0*P(I,J)**A1 CONTRA1A.110
ENDDO CONTRA1A.111
ENDDO CONTRA1A.112
C----------------------------------------------------------------------- CONTRA1A.113
CL Loop round all points and all levels CONTRA1A.114
C----------------------------------------------------------------------- CONTRA1A.115
DO 2 J=1,P_LEVELS CONTRA1A.116
DO 1 I=START,END GSM1F405.605
IF (J.EQ.1) THEN CONTRA1A.118
C----------------------------------------------------------------------- CONTRA1A.119
CL Check if level 1 temperature is less than MINTRA-14 CONTRA1A.120
C----------------------------------------------------------------------- CONTRA1A.121
IF (MINTRA_M_14(I,1).GT.T(I,1)) THEN CONTRA1A.122
P_M(I)=PSTAR(I) CONTRA1A.123
INTERSECTION(I)=.TRUE. CONTRA1A.124
ENDIF CONTRA1A.125
ELSEIF(P(I,J).GE.P_LIMIT)THEN CONTRA1A.126
C----------------------------------------------------------------------- CONTRA1A.127
CL Continue if below pressure level P_LIMIT CONTRA1A.128
CL Calculate 'MINTRA line minus 14 degrees' - 'environment curve' CONTRA1A.129
CL at J-1 and J CONTRA1A.130
C----------------------------------------------------------------------- CONTRA1A.131
TERM3=MINTRA_M_14(I,J)-T(I,J) CONTRA1A.132
TERM4=MINTRA_M_14(I,J-1)-T(I,J-1) CONTRA1A.133
C----------------------------------------------------------------------- CONTRA1A.134
CL If TERM3 and TERM4 have different signs or either is equal to zero CONTRA1A.135
CL there is a point of intersection between levels j-1 and j ; CONTRA1A.136
CL otherwise there is not. CONTRA1A.137
C----------------------------------------------------------------------- CONTRA1A.138
TEST1=TERM3*TERM4 CONTRA1A.139
IF (TEST1.LE.0.0) THEN CONTRA1A.140
INTERSECTION(I)=.TRUE. CONTRA1A.141
C----------------------------------------------------------------------- CONTRA1A.142
CL Exner pressure difference across layers CONTRA1A.143
C----------------------------------------------------------------------- CONTRA1A.144
DEL_EXNER_J=P_EXNER_HALF(I,J)-P_EXNER_HALF(I,J+1) CONTRA1A.145
DEL_EXNER_JM1=P_EXNER_HALF(I,J-1)-P_EXNER_HALF(I,J) CONTRA1A.146
C----------------------------------------------------------------------- CONTRA1A.147
CL Numerator CONTRA1A.148
C----------------------------------------------------------------------- CONTRA1A.149
TERM1=T(I,J-1)-T(I,J) CONTRA1A.150
C----------------------------------------------------------------------- CONTRA1A.151
CL Denominator CONTRA1A.152
C----------------------------------------------------------------------- CONTRA1A.153
TERM2=CP_OVER_TWO_G*(THETA(I,J-1)*DEL_EXNER_JM1 CONTRA1A.154
* +THETA(I,J)*DEL_EXNER_J) CONTRA1A.155
C----------------------------------------------------------------------- CONTRA1A.156
CL Lapse rate between level j-1 and j =TERM1/TERM2 CONTRA1A.157
CL Exponent=Lapse rate *(R/g) CONTRA1A.158
C----------------------------------------------------------------------- CONTRA1A.159
EXPONENT=TERM1/TERM2*R/G CONTRA1A.160
C----------------------------------------------------------------------- CONTRA1A.161
CL Find P_M, the pressure of intersection CONTRA1A.162
C----------------------------------------------------------------------- CONTRA1A.163
TERM5=T(I,J-1)/A0*P(I,J-1)**(-EXPONENT) CONTRA1A.164
TERM6=A1-EXPONENT CONTRA1A.165
IF (ABS(TERM6).LT.1.0E-5) TERM6=1.0E-5 GDR2F404.17
P_M(I)=TERM5**(1.0/TERM6) CONTRA1A.167
C----------------------------------------------------------------------- CONTRA1A.168
CL Check that P_M lies between P(I,J) and P(I,J-1) CONTRA1A.169
CL If not set P_M to the upper level P(I,J) CONTRA1A.170
C----------------------------------------------------------------------- CONTRA1A.171
TEST2=(P_M(I)-P(I,J))*(P_M(I)-P(I,J-1)) CONTRA1A.172
IF(TEST2.GT.0.0) P_M(I)=P(I,J) CONTRA1A.173
ENDIF CONTRA1A.174
ENDIF CONTRA1A.175
1 CONTINUE CONTRA1A.176
C----------------------------------------------------------------------- CONTRA1A.177
CL Is P_M a maximum or minimum pressure of intersection? CONTRA1A.178
C i.e. a solution to the equation at this point. CONTRA1A.179
C----------------------------------------------------------------------- CONTRA1A.180
DO I=START,END GSM1F405.606
IF (INTERSECTION(I)) THEN CONTRA1A.182
IF (P_M(I).GT.P_MAX(I)) P_MAX(I)=P_M(I) CONTRA1A.183
IF (P_M(I).LT.P_MIN(I)) P_MIN(I)=P_M(I) CONTRA1A.184
ENDIF CONTRA1A.185
ENDDO CONTRA1A.186
2 CONTINUE CONTRA1A.187
C----------------------------------------------------------------------- CONTRA1A.188
CL If only one one intersection set P_MIN to zero CONTRA1A.189
CL If no intersections set P_MAX and P_MIN to zero CONTRA1A.190
C----------------------------------------------------------------------- CONTRA1A.191
DO I=START,END GSM1F405.607
IF (INTERSECTION(I).AND.(P_MAX(I).EQ.P_MIN(I))) P_MIN(I)=0.0 CONTRA1A.193
IF (.NOT.INTERSECTION(I)) THEN CONTRA1A.194
P_MAX(I)=0.0 CONTRA1A.195
P_MIN(I)=0.0 CONTRA1A.196
ENDIF CONTRA1A.197
ENDDO CONTRA1A.198
C----------------------------------------------------------------------- CONTRA1A.199
CL Convert P_MAX and P_MIN into ICAO heights in thousands of feet CONTRA1A.200
C----------------------------------------------------------------------- CONTRA1A.201
CALL ICAO_HT(P_MAX(START),END-START+1,LOWER_CONTRAIL(START)) GSM1F405.608
CALL ICAO_HT(P_MIN(START),END-START+1,UPPER_CONTRAIL(START)) GSM1F405.609
C----------------------------------------------------------------------- CONTRA1A.204
RETURN CONTRA1A.205
END CONTRA1A.206
C----------------------------------------------------------------------- CONTRA1A.207
*ENDIF CONTRA1A.208