*IF DEF,A13_1B THQDIF1B.2
C ******************************COPYRIGHT****************************** THQDIF1B.3
C (c) CROWN COPYRIGHT 1996, METEOROLOGICAL OFFICE, All Rights Reserved. THQDIF1B.4
C THQDIF1B.5
C Use, duplication or disclosure of this code is subject to the THQDIF1B.6
C restrictions as set forth in the contract. THQDIF1B.7
C THQDIF1B.8
C Meteorological Office THQDIF1B.9
C London Road THQDIF1B.10
C BRACKNELL THQDIF1B.11
C Berkshire UK THQDIF1B.12
C RG12 2SZ THQDIF1B.13
C THQDIF1B.14
C If no contract has been raised with this copy of the code, the use, THQDIF1B.15
C duplication or disclosure of it is strictly prohibited. Permission THQDIF1B.16
C to do so must first be obtained in writing from the Head of Numerical THQDIF1B.17
C Modelling at the above address. THQDIF1B.18
C ******************************COPYRIGHT****************************** THQDIF1B.19
C THQDIF1B.20
CLL SUBROUTINE TH_Q_DIF ----------------------------------------- THQDIF1B.21
CLL THQDIF1B.22
CLL PURPOSE: CALCULATES DIFFUSION INCREMENTS FOR THETAL OR QT THQDIF1B.23
CLL IF STEEP SLOPE THEN EFFECTIVE DIFFUSION IS ZERO. THQDIF1B.24
CLL THQDIF1B.25
CLL NOT SUITABLE FOR SINGLE COLUMN USE. THQDIF1B.26
CLL THQDIF1B.27
CLL VERSION FOR CRAY Y-MP THQDIF1B.28
CLL THQDIF1B.29
CLL MM, TJ <- PROGRAMMER OF SOME OR ALL OF PREVIOUS CODE OR CHANGES THQDIF1B.30
CLL THQDIF1B.31
CLL MODEL MODIFICATION HISTORY FROM MODEL VERSION 4.1: THQDIF1B.32
CLL VERSION DATE THQDIF1B.33
!LL 4.2 28/11/96 New deck for HADCM2-specific section A13_1B, THQDIF1B.34
!LL as THQDIF1A but with reintroduced 'bug' in THQDIF1B.35
!LL value of SCALAR weight at poles. T.Johns THQDIF1B.36
CLL THQDIF1B.37
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4, THQDIF1B.38
CLL STANDARD B. VERSION 2, DATED 18/01/90 THQDIF1B.39
CLL THQDIF1B.40
CLL SYSTEM COMPONENTS COVERED: P131 THQDIF1B.41
CLL THQDIF1B.42
CLL SYSTEM TASK: P1 THQDIF1B.43
CLL THQDIF1B.44
CLL DOCUMENTATION: THE EQUATION USED IS (47) THQDIF1B.45
CLL IN UNIFIED MODEL DOCUMENTATION PAPER THQDIF1B.46
CLL NO. 10 M.J.P. CULLEN,T.DAVIES AND M.H.MAWSON THQDIF1B.47
CLL VERSION 16, DATED 09/01/91. THQDIF1B.48
CLLEND------------------------------------------------------------- THQDIF1B.49
THQDIF1B.50
C*L ARGUMENTS:--------------------------------------------------- THQDIF1B.51
SUBROUTINE TH_Q_DIF 4,5THQDIF1B.52
1 (FIELD,FIELD_INC, THQDIF1B.53
2 SEC_P_LATITUDE,ROW_LENGTH, THQDIF1B.54
*CALL ARGFLDPT 
THQDIF1B.55
4 P_FIELD,U_FIELD, THQDIF1B.56
5 DIFFUSION_EW,DIFFUSION_NS) THQDIF1B.57
THQDIF1B.58
THQDIF1B.59
IMPLICIT NONE THQDIF1B.60
THQDIF1B.61
INTEGER THQDIF1B.62
* U_FIELD !IN DIMENSION OF FIELDS ON VELOCITY GRID THQDIF1B.63
*, P_FIELD !IN DIMENSION OF FIELDS ON PRESSURE GRID THQDIF1B.64
*, ROW_LENGTH !IN NUMBER OF POINTS PER ROW THQDIF1B.65
THQDIF1B.66
! All TYPFLDPT arguments are intent IN THQDIF1B.67
*CALL TYPFLDPT THQDIF1B.68
THQDIF1B.69
REAL THQDIF1B.70
* FIELD(P_FIELD) !IN. THETAL OR QT FIELD. THQDIF1B.71
*, FIELD_INC(P_FIELD) !OUT DIFFUSION INCREMENT THQDIF1B.72
THQDIF1B.73
REAL THQDIF1B.74
* DIFFUSION_EW(P_FIELD) !IN HOLDS EAST-WEST EFFECTIVE DIFFUSION THQDIF1B.75
*,DIFFUSION_NS(P_FIELD) !IN HOLDS NORTH-SOUTH EFFECTIVE DIFFUSION THQDIF1B.76
*,SEC_P_LATITUDE(P_FIELD) !IN 1/COS(LAT) AT P POINTS THQDIF1B.77
THQDIF1B.78
C*--------------------------------------------------------------------- THQDIF1B.79
THQDIF1B.80
C*L DEFINE ARRAYS AND VARIABLES USED IN THIS ROUTINE----------------- THQDIF1B.81
C DEFINE LOCAL ARRAYS: 5 ARE REQUIRED THQDIF1B.82
THQDIF1B.83
REAL THQDIF1B.84
* FIELD1(P_FIELD) ! GENERAL WORKSPACE THQDIF1B.85
*,FIELD2(P_FIELD) ! GENERAL WORKSPACE THQDIF1B.86
*,FIELD3(P_FIELD) ! GENERAL WORKSPACE THQDIF1B.87
*,NP_FLUX(ROW_LENGTH) ! HOLDS NORTH POLAR FLUX THQDIF1B.88
*,SP_FLUX(ROW_LENGTH) ! HOLDS SOUTH POLAR FLUX THQDIF1B.89
C*--------------------------------------------------------------------- THQDIF1B.90
C DEFINE LOCAL VARIABLES THQDIF1B.91
THQDIF1B.92
C LOCAL REALS. THQDIF1B.93
REAL THQDIF1B.94
* SCALAR THQDIF1B.95
C COUNT VARIABLES FOR DO LOOPS ETC. THQDIF1B.96
INTEGER THQDIF1B.97
* I,J,IJ THQDIF1B.98
THQDIF1B.99
C*L EXTERNAL SUBROUTINE CALLS:--------------------------------------- THQDIF1B.100
EXTERNAL THQDIF1B.101
* POLAR THQDIF1B.102
C*--------------------------------------------------------------------- THQDIF1B.103
CL MAXIMUM VECTOR LENGTH ASSUMED IS END_P_UPDATE-START_P_UPDATE+1 THQDIF1B.104
CL--------------------------------------------------------------------- THQDIF1B.105
CL INTERNAL STRUCTURE. THQDIF1B.106
CL--------------------------------------------------------------------- THQDIF1B.107
CL THQDIF1B.108
CL--------------------------------------------------------------------- THQDIF1B.109
CL SECTION 1. DELTA LAMBDA TERMS THQDIF1B.110
CL--------------------------------------------------------------------- THQDIF1B.111
C---------------------------------------------------------------------- THQDIF1B.112
CL SECTION 1.1 CALCULATE DELTAPHILAMBDA*1/(DELTALAMBDA)SQUARED THQDIF1B.113
C---------------------------------------------------------------------- THQDIF1B.114
THQDIF1B.115
THQDIF1B.116
DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO-1 THQDIF1B.117
FIELD1(I)=FIELD(I+1)-FIELD(I) THQDIF1B.118
END DO THQDIF1B.119
THQDIF1B.120
THQDIF1B.121
*IF -DEF,MPP THQDIF1B.122
C RECALCULATE END-POINTS THQDIF1B.123
THQDIF1B.124
DO I=START_POINT_NO_HALO-1,END_P_POINT_NO_HALO,ROW_LENGTH THQDIF1B.125
IJ=I-ROW_LENGTH+1 THQDIF1B.126
FIELD1(I)=FIELD(IJ)-FIELD(I) THQDIF1B.127
END DO THQDIF1B.128
*ELSE THQDIF1B.129
! Set last point of field THQDIF1B.130
FIELD1(END_P_POINT_NO_HALO)=FIELD1(END_P_POINT_NO_HALO-1) THQDIF1B.131
*ENDIF THQDIF1B.132
THQDIF1B.133
C---------------------------------------------------------------------- THQDIF1B.134
CL SECTION 1.2 CALCULATE DELTA LAMBDA TERM THQDIF1B.135
C---------------------------------------------------------------------- THQDIF1B.136
THQDIF1B.137
DO I= START_POINT_NO_HALO+1,END_P_POINT_NO_HALO THQDIF1B.138
FIELD2(I)=(DIFFUSION_EW(I)*FIELD1(I)- THQDIF1B.139
& DIFFUSION_EW(I-1)*FIELD1(I-1))* THQDIF1B.140
& SEC_P_LATITUDE(I) THQDIF1B.141
END DO THQDIF1B.142
THQDIF1B.143
THQDIF1B.144
*IF -DEF,MPP THQDIF1B.145
C RECALCULATE END-POINTS THQDIF1B.146
DO I= START_POINT_NO_HALO,END_P_POINT_NO_HALO,ROW_LENGTH THQDIF1B.147
IJ=I+ROW_LENGTH-1 THQDIF1B.148
FIELD2(I)=(DIFFUSION_EW(I)*FIELD1(I)- THQDIF1B.149
& DIFFUSION_EW(IJ)*FIELD1(IJ))* THQDIF1B.150
& SEC_P_LATITUDE(I) THQDIF1B.151
END DO THQDIF1B.152
*ELSE THQDIF1B.153
FIELD2(START_POINT_NO_HALO)=FIELD2(START_POINT_NO_HALO+1) THQDIF1B.154
*ENDIF THQDIF1B.155
THQDIF1B.156
C---------------------------------------------------------------------- THQDIF1B.157
CL SECTION 2 CALCULATE PHI DIRECTION TERM AND ADD THQDIF1B.158
CL ONTO FIRST TERM TO GET TOTAL CORRECTION. THQDIF1B.159
C---------------------------------------------------------------------- THQDIF1B.160
THQDIF1B.161
C CALCULATE DELTA PHI TERMS THQDIF1B.162
THQDIF1B.163
DO I=START_POINT_NO_HALO-ROW_LENGTH,END_P_POINT_NO_HALO THQDIF1B.164
FIELD1(I)=FIELD(I)-FIELD(I+ROW_LENGTH) THQDIF1B.165
END DO THQDIF1B.166
THQDIF1B.167
C---------------------------------------------------------------------- THQDIF1B.168
CL SECTION 2.3 CALCULATE DELTAPHI TERM AND ADD ONTO DELTALAMBDA TERM THQDIF1B.169
C---------------------------------------------------------------------- THQDIF1B.170
THQDIF1B.171
DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO THQDIF1B.172
FIELD_INC(I)= (FIELD2(I)+ THQDIF1B.173
& DIFFUSION_NS(I-ROW_LENGTH)*FIELD1(I-ROW_LENGTH)- THQDIF1B.174
& DIFFUSION_NS(I)*FIELD1(I))*SEC_P_LATITUDE(I) THQDIF1B.175
END DO THQDIF1B.176
THQDIF1B.177
C---------------------------------------------------------------------- THQDIF1B.178
CL SECTION 3 CALCULATE DIFFUSION AT POLES THQDIF1B.179
C---------------------------------------------------------------------- THQDIF1B.180
THQDIF1B.181
*IF DEF,GLOBAL THQDIF1B.182
THQDIF1B.183
C GLOBAL MODEL CALCULATES POLAR DEL-SQUARED USING ACROSS-POLE DIFFERENCE THQDIF1B.184
C ASSUMING FLUX=0 AT HALF-GRID-LENGTH ON OTHER SIDE OF POLE THQDIF1B.185
THQDIF1B.186
THQDIF1B.187
! Note - 2.0* factor is incorrect but consistent with HADCM2 THQDIF1B.188
SCALAR=2.0*SEC_P_LATITUDE(TOP_ROW_START) THQDIF1B.189
THQDIF1B.190
! Do North Pole THQDIF1B.191
*IF DEF,MPP THQDIF1B.192
IF (at_top_of_LPG) THEN THQDIF1B.193
*ENDIF THQDIF1B.194
DO I=1,ROW_LENGTH THQDIF1B.195
J=TOP_ROW_START+I-1 THQDIF1B.196
FIELD3(J)=-DIFFUSION_NS(J)*FIELD1(J)*SCALAR THQDIF1B.197
NP_FLUX(I)=FIELD3(J) THQDIF1B.198
FIELD_INC(J)=0.0 THQDIF1B.199
ENDDO THQDIF1B.200
*IF DEF,MPP THQDIF1B.201
ENDIF THQDIF1B.202
*ENDIF THQDIF1B.203
THQDIF1B.204
! And South Pole THQDIF1B.205
*IF DEF,MPP THQDIF1B.206
! Note - 2.0* factor is incorrect but consistent with HADCM2 THQDIF1B.207
SCALAR=2.0*SEC_P_LATITUDE(P_BOT_ROW_START) THQDIF1B.208
THQDIF1B.209
IF (at_base_of_LPG) THEN THQDIF1B.210
*ENDIF THQDIF1B.211
DO I=1,ROW_LENGTH THQDIF1B.212
J=I+P_BOT_ROW_START-1 THQDIF1B.213
FIELD3(J)=DIFFUSION_NS(J-ROW_LENGTH)*FIELD1(J-ROW_LENGTH)* THQDIF1B.214
& SCALAR THQDIF1B.215
SP_FLUX(I)=FIELD3(J) THQDIF1B.216
FIELD_INC(J)=0.0 THQDIF1B.217
ENDDO THQDIF1B.218
*IF DEF,MPP THQDIF1B.219
ENDIF THQDIF1B.220
*ENDIF THQDIF1B.221
THQDIF1B.222
CL CALL POLAR TO UPDATE POLAR VALUES THQDIF1B.223
THQDIF1B.224
CALL POLAR(FIELD_INC,NP_FLUX,SP_FLUX, THQDIF1B.225
*CALL ARGFLDPT THQDIF1B.226
& P_FIELD,ROW_LENGTH,ROW_LENGTH, THQDIF1B.227
& 1,1,ROW_LENGTH,1) THQDIF1B.228
THQDIF1B.229
*ELSE THQDIF1B.230
CL LIMITED AREA MODEL ZEROES DEL-SQUARED ON BOUNDARIES. THQDIF1B.231
*IF DEF,MPP THQDIF1B.232
IF (at_left_of_LPG) THEN THQDIF1B.233
*ENDIF THQDIF1B.234
DO I=START_POINT_NO_HALO+FIRST_ROW_PT-1, THQDIF1B.235
& END_P_POINT_NO_HALO,ROW_LENGTH THQDIF1B.236
FIELD_INC(I)=0.0 THQDIF1B.237
ENDDO THQDIF1B.238
*IF DEF,MPP THQDIF1B.239
ENDIF THQDIF1B.240
THQDIF1B.241
IF (at_right_of_LPG) THEN THQDIF1B.242
*ENDIF THQDIF1B.243
DO I=START_POINT_NO_HALO+LAST_ROW_PT-1, THQDIF1B.244
& END_P_POINT_NO_HALO,ROW_LENGTH THQDIF1B.245
FIELD_INC(I)=0.0 THQDIF1B.246
ENDDO THQDIF1B.247
*IF DEF,MPP THQDIF1B.248
ENDIF THQDIF1B.249
*ENDIF THQDIF1B.250
*ENDIF THQDIF1B.251
CL END OF ROUTINE TH_Q_DIF THQDIF1B.252
THQDIF1B.253
RETURN THQDIF1B.254
END THQDIF1B.255
*ENDIF THQDIF1B.256