*IF DEF,A16_1A THETAW1A.2
C ******************************COPYRIGHT****************************** GTS2F400.10207
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.10208
C GTS2F400.10209
C Use, duplication or disclosure of this code is subject to the GTS2F400.10210
C restrictions as set forth in the contract. GTS2F400.10211
C GTS2F400.10212
C Meteorological Office GTS2F400.10213
C London Road GTS2F400.10214
C BRACKNELL GTS2F400.10215
C Berkshire UK GTS2F400.10216
C RG12 2SZ GTS2F400.10217
C GTS2F400.10218
C If no contract has been raised with this copy of the code, the use, GTS2F400.10219
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.10220
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.10221
C Modelling at the above address. GTS2F400.10222
C ******************************COPYRIGHT****************************** GTS2F400.10223
C GTS2F400.10224
CLL Subroutine THETAW-------------------------------------------------- THETAW1A.3
CLL THETAW1A.4
CLL Purpose: To calculate from temperatures,humidities the wet bulb THETAW1A.5
CLL potential temperature. THETAW1A.6
CLL THETAW1A.7
CLL D.Robinson <- programmer of some or all of previous code or changes THETAW1A.8
CLL THETAW1A.9
CLL Model Modification history from model version 3.0: THETAW1A.10
CLL version Date THETAW1A.11
CLL 3.2 13/07/93 Changed CHARACTER*(*) to CHARACTER*(80) for TS150793.204
CLL portability. Author Tracey Smith. TS150793.205
!LL 4.5 15/04/98 Added START,END arguments to enable duplicate GSM1F405.799
!LL halo calculations to be avoided. S.D.Mullerworth GSM1F405.800
CLL THETAW1A.12
CLL Programming standard : UM Doc Paper no 3 THETAW1A.13
CLL THETAW1A.14
CLL Logical components covered : D482 THETAW1A.15
CLL THETAW1A.16
CLL Project task : D482 THETAW1A.17
CLL THETAW1A.18
CLL External documentation : UMDP no THETAW1A.19
CLL THETAW1A.20
CLLEND ----------------------------------------------------------------- THETAW1A.21
C THETAW1A.22
C*L ARGUMENTS:----------------------------------------------------- THETAW1A.23
SUBROUTINE THETAW 1,4THETAW1A.24
1 (PRESS_REQD,THETA,Q,P,PSTAR,P_EXNER_HALF,TW, THETAW1A.25
2 AK,BK,AKH,BKH,P_FIELD,P_LEVELS,Q_LEVELS,T_REF, THETAW1A.26
& START,END, GSM1F405.801
3 CMESSAGE,ICODE) THETAW1A.27
C THETAW1A.28
IMPLICIT NONE THETAW1A.29
C THETAW1A.30
INTEGER THETAW1A.31
* P_FIELD ! IN No of points on a field THETAW1A.32
*,P_LEVELS ! IN NO OF MODEL LEVELS THETAW1A.33
*,Q_LEVELS ! IN NO OF HUMIDITY LEVELS THETAW1A.34
*,T_REF ! IN LEVEL OF MODEL USED TO CALC HEIGHT THETAW1A.35
&,START,END ! IN Range of points to calculate GSM1F405.802
*,ICODE ! OUT RETURN CODE THETAW1A.36
C THETAW1A.37
CHARACTER*(80) CMESSAGE TS150793.206
C THETAW1A.39
REAL THETAW1A.40
* P(P_FIELD,P_LEVELS) ! IN Pressure at each level and point THETAW1A.41
*,THETA(P_FIELD,P_LEVELS)! IN Potential temperature on all pts THETAW1A.42
*,P_EXNER_HALF(P_FIELD,P_LEVELS+1) ! IN EXNER Pressure at model THETAW1A.43
* ! half levels THETAW1A.44
*,PSTAR(P_FIELD) ! IN Surface pressure THETAW1A.45
*,Q(P_FIELD,Q_LEVELS) ! IN Specific humidity at full levels THETAW1A.46
*,MODEL_HALF_HEIGHT(P_FIELD,P_LEVELS+1)! IN Height of 1/2 levels THETAW1A.47
*,AKH(P_LEVELS+1) ! IN Value of "A" at mid layer THETAW1A.48
*,BKH(P_LEVELS+1) ! IN Value of "B" at mid layer THETAW1A.49
*,AK (P_LEVELS) ! IN Value of "A" at model level THETAW1A.50
*,BK (P_LEVELS) ! IN Value of "B" at model level THETAW1A.51
*,PRESS_REQD ! IN Pressure level required for THETAW THETAW1A.52
C THETAW1A.53
REAL THETAW1A.54
* TW(P_FIELD) ! OUT The WET BULB temperature. On THETAW1A.55
C ! output the TW at 1000mb ie thetaw THETAW1A.56
C THETAW1A.57
C*--------------------------------------------------------------- THETAW1A.58
C THETAW1A.59
C*L WORKSPACE USAGE---------------------------------------------- THETAW1A.60
C 10*P_FIELD THETAW1A.61
C DEFINE LOCAL WORKSPACE ARRAYS THETAW1A.62
REAL THETAW1A.63
* T(P_FIELD) ! Intial temperature at Press_reqd (Ta) THETAW1A.64
*,TW_TEMP(P_FIELD) ! A temporary store for temperature. THETAW1A.65
*,QAD(P_FIELD) ! SPECIFIC HUMIDITY AT PRESS_REQD THETAW1A.66
*,Q_SAT(P_FIELD) ! Saturation specific humidityEQD THETAW1A.67
*,PRESS_REQD_HORIZ(P_FIELD) ! PRESS_REQD BROADCAST OVER DOMAIN THETAW1A.68
*,LATENT_HEAT(P_FIELD) ! Latent heat of evaporation (fn(T)) THETAW1A.69
*,GG(P_FIELD) ! The 'G' used in equation (1.3) THETAW1A.70
*,GS ! The 'G' used in equation (1.3) THETAW1A.71
*,DIFF(P_FIELD) ! The difference between G(Tw) & G(Ti) THETAW1A.72
*,DGBYDT ! The function DG/DT Eqn (1.6) THETAW1A.73
*,FF(P_FIELD) ! The function DT/DP Eqn (1.5) THETAW1A.74
*,FF_NEXT(P_FIELD) ! The function DT/DP THETAW1A.75
C*--------------------------------------------------------------- THETAW1A.76
C THETAW1A.77
C*L EXTERNAL SUBROUTINES CALLED---------------------------------- THETAW1A.78
EXTERNAL V_INT,V_INT_T,QSAT THETAW1A.79
C*--------------------------------------------------------------- THETAW1A.80
C THETAW1A.81
*CALL C_R_CP 
THETAW1A.82
*CALL C_G THETAW1A.83
*CALL C_0_DG_C THETAW1A.84
*CALL C_LHEAT THETAW1A.85
*CALL C_EPSLON THETAW1A.86
C THETAW1A.87
C---------------------------------------------------------------- THETAW1A.88
C DEFINE LOCAL VARIABLES THETAW1A.89
C---------------------------------------------------------------- THETAW1A.90
INTEGER THETAW1A.91
* I,K ! LOOP COUNTERS THETAW1A.92
*,J ! LOWEST ETA LEVEL WITH TEMP<T0 THETAW1A.93
*,LOOP THETAW1A.94
*,NLOOP THETAW1A.95
*,LOOP_FIELD(P_FIELD) THETAW1A.96
C THETAW1A.97
REAL THETAW1A.98
* DUMMY1 THETAW1A.99
*,DUMMY2 THETAW1A.100
*,COEFF ! COEFF USED IN THE CALCULATION OF LATENT HEAT THETAW1A.101
*,CPV ! SPECIFIC HEAT FOR WATER VAPOUR THETAW1A.102
*,MV ! Mol wt of water vapour THETAW1A.103
*,MD ! Mol wt of dry air THETAW1A.104
*,RSTAR ! Universal gas constant THETAW1A.105
*,TEMP1 ! THETAW1A.106
*,TEMP2 ! THETAW1A.107
*,PP ! THETAW1A.108
*,DELTAP ! Pressure increment THETAW1A.109
*,P_TARGET ! The target pressure (1000mb) THETAW1A.110
*,REAL_NLOOP ! Real NLOOP THETAW1A.111
*,REM ! Temp local variable THETAW1A.112
THETAW1A.113
DATA COEFF/2.34E3/ ! Used in the calculation of LATENT_H THETAW1A.114
DATA CPV/1850.0/ ! Specific heat of water vapour THETAW1A.115
DATA MV/0.01801/ ! Mol wt of water vapour KG/MOL THETAW1A.116
DATA MD/0.02896/ ! Mol wt of dry air KG/MOL THETAW1A.117
DATA RSTAR/8.314/ ! Universal gas constant THETAW1A.118
C THETAW1A.119
C THETAW1A.120
C------------------------------------------------------------------ THETAW1A.121
CL 1. Calculate the temperature and specific hum at the PRESS_REQ THETAW1A.122
C------------------------------------------------------------------ THETAW1A.123
THETAW1A.124
DO 1 I=START,END GSM1F405.803
PRESS_REQD_HORIZ(I)=PRESS_REQD THETAW1A.126
1 CONTINUE THETAW1A.127
THETAW1A.128
CALL V_INT_T(T,PRESS_REQD_HORIZ,P(1,T_REF),PSTAR, THETAW1A.129
& P_EXNER_HALF,THETA,P_FIELD,P_LEVELS,T_REF,AKH,BKH GSM1F405.804
& ,START,END) GSM1F405.805
THETAW1A.131
THETAW1A.132
CALL V_INT(P,PRESS_REQD_HORIZ,Q,QAD,P_FIELD,Q_LEVELS, THETAW1A.133
& DUMMY1,DUMMY2,.FALSE.,START,END) GSM1F405.806
THETAW1A.135
THETAW1A.136
C--------------------------------------------------------------------- THETAW1A.137
CL---------------------SECTION 2 ------------------------------------- THETAW1A.138
CL Calculate the function G for TA and QA (see eqn 3) doc no THETAW1A.139
C G(Tw)=Qa(Ta)*L(Ta)+Ta(Cp+QaCpv) THETAW1A.140
C Subscript a indicates initial values THETAW1A.141
C--------------------------------------------------------------------- THETAW1A.142
THETAW1A.143
THETAW1A.144
DO 20 I=START,END GSM1F405.807
LATENT_HEAT(I)=LC-COEFF*(T(I)-ZERODEGC) THETAW1A.146
GG(I)=QAD(I)*LATENT_HEAT(I)+T(I)*(CP+QAD(I)*CPV) ! Eqn 1.2 THETAW1A.147
20 CONTINUE THETAW1A.148
THETAW1A.149
C--------------------------------------------------------------------- THETAW1A.150
CL---------------------SECTION 3 ------------------------------------- THETAW1A.151
CL Calculate the function DG/DT (see eqn 5) doc no THETAW1A.152
CL G'(T)=(Mv*L(Ta)**2*Qs(T)/R*T**2) +Cp+Qa*Cpv THETAW1A.153
CL Iterate to find the Tw using 6 THETAW1A.154
CL T(i+1)=T(i)+(G(Tw)-G(Ti))/DG/DT(i) THETAW1A.155
C--------------------------------------------------------------------- THETAW1A.156
THETAW1A.157
DO 30 I=START,END GSM1F405.808
TW(I)=T(I) THETAW1A.159
30 CONTINUE THETAW1A.160
THETAW1A.161
LOOP=1 THETAW1A.162
10000 CONTINUE THETAW1A.163
THETAW1A.164
CALL QSAT(Q_SAT(START),TW(START),PRESS_REQD_HORIZ(START) GSM1F405.809
& ,END-START+1) !Qs @ required Ta & P GSM1F405.810
C THETAW1A.166
DO 31 I=START,END GSM1F405.811
GS=Q_SAT(I)*LATENT_HEAT(I)+TW(I)*(CP+QAD(I)*CPV) ! Eqn 1.3 THETAW1A.168
THETAW1A.169
CL---------------------SECTION 3.1------------------------------------ THETAW1A.170
CL Find the difference between G(Tw)-G(Ti) THETAW1A.171
C--------------------------------------------------------------------- THETAW1A.172
THETAW1A.173
DIFF(I)=ABS(GG(I)-GS) ! Eqn 1.8 THETAW1A.174
IF(DIFF(I).GT.1.0) THEN THETAW1A.175
THETAW1A.176
CL---------------------SECTION 3.1------------------------------------ THETAW1A.177
CL Calculate the function DG/DT THETAW1A.178
C--------------------------------------------------------------------- THETAW1A.179
THETAW1A.180
TEMP1=RSTAR*TW(I)**2 THETAW1A.181
TEMP2=Q_SAT(I)*MV*LATENT_HEAT(I)**2 THETAW1A.182
DGBYDT=TEMP2/TEMP1 + CP + CPV*QAD(I) ! Eqn 1.6 THETAW1A.183
THETAW1A.184
CL---------------------SECTION 3.2------------------------------------ THETAW1A.185
CL Using the function DG/DT calculate an updated Temperature Tw THETAW1A.186
C--------------------------------------------------------------------- THETAW1A.187
THETAW1A.188
TW(I) = TW(I) - (GS-GG(I)) / DGBYDT ! Eqn 1.7 THETAW1A.189
THETAW1A.190
CL---------------------SECTION 3.3------------------------------------ THETAW1A.191
CL Using the new temperature Tw re-calculate GS First update LATENT_H THETAW1A.192
C--------------------------------------------------------------------- THETAW1A.193
THETAW1A.194
LATENT_HEAT(I)=LC-COEFF*(TW(I)-ZERODEGC) THETAW1A.195
THETAW1A.196
ENDIF THETAW1A.197
31 CONTINUE THETAW1A.198
THETAW1A.199
LOOP=LOOP+1 THETAW1A.200
THETAW1A.201
IF(LOOP.GT.10) THEN THETAW1A.202
ICODE=1 THETAW1A.203
CMESSAGE='THETAW Convergence failure in THETAW' THETAW1A.204
GOTO 9999 THETAW1A.205
ENDIF THETAW1A.206
THETAW1A.207
DO 32 I=START,END GSM1F405.812
IF(DIFF(I).GT.1.0) GOTO 10000 THETAW1A.209
32 CONTINUE THETAW1A.210
THETAW1A.211
C--------------------------------------------------------------------- THETAW1A.212
CL---------------------SECTION 4 ------------------------------------- THETAW1A.213
CL Having calculated now have to descend the sat adiabat down THETAW1A.214
CL to 1000mb solving:- THETAW1A.215
CL DT=(L*EPS/P*Es +R*T/Md)/(Cp*P+(EPS*Mv*L**2*Es/(R*T**2)) *DP THETAW1A.216
C--------------------------------------------------------------------- THETAW1A.217
THETAW1A.218
PP=PRESS_REQD THETAW1A.219
DELTAP=5000.0 ! Pressure increment in PASCALS THETAW1A.220
P_TARGET=100000.0 ! 1000 mb THETAW1A.221
REAL_NLOOP=(P_TARGET-PP)/DELTAP + 0.000001 THETAW1A.222
NLOOP=REAL_NLOOP THETAW1A.223
REM=REAL_NLOOP-NLOOP THETAW1A.224
IF(REM.GT.0.00001) THEN THETAW1A.225
NLOOP=NLOOP+1 THETAW1A.226
DELTAP=(P_TARGET-PP)/NLOOP THETAW1A.227
ENDIF THETAW1A.228
THETAW1A.229
CL Calculate DT/DP = F(T,P) THETAW1A.230
GSM1F405.813
DO 41 K=1,NLOOP THETAW1A.232
DO 40 I=START,END GSM1F405.814
TEMP1=(Q_SAT(I)*PP*MV*LATENT_HEAT(I)**2/(RSTAR*TW(I)**2)) + CP*PP THETAW1A.234
TEMP2=Q_SAT(I)*LATENT_HEAT(I)+RSTAR*TW(I)/MD THETAW1A.235
FF(I)=TEMP2/TEMP1 ! Eqn 1.5 THETAW1A.236
TW_TEMP(I)=FF(I)*DELTAP + TW(I) ! Eqn 1.10 THETAW1A.237
40 CONTINUE THETAW1A.238
THETAW1A.239
PP=PP+DELTAP THETAW1A.240
DO 42 I=START,END GSM1F405.815
PRESS_REQD_HORIZ(I)=PP THETAW1A.242
LATENT_HEAT(I)=LC-COEFF*(TW_TEMP(I)-ZERODEGC) THETAW1A.243
42 CONTINUE THETAW1A.244
THETAW1A.245
CALL QSAT(Q_SAT(START),TW_TEMP(START),PRESS_REQD_HORIZ(START) GSM1F405.816
& ,END-START+1) !Qs @ Tw & P GSM1F405.817
THETAW1A.247
DO 43 I=START,END GSM1F405.818
TEMP1=(Q_SAT(I)*PP*MV*LATENT_HEAT(I)**2/(RSTAR*TW_TEMP(I)**2)) THETAW1A.249
* + CP*PP THETAW1A.250
TEMP2=Q_SAT(I)*LATENT_HEAT(I)+RSTAR*TW_TEMP(I)/MD THETAW1A.251
FF_NEXT(I)=TEMP2/TEMP1 THETAW1A.252
TW(I)=(FF(I)+FF_NEXT(I))*0.5*DELTAP + TW(I) ! Eqn 1.11 THETAW1A.253
LATENT_HEAT(I)=LC-COEFF*(TW(I)-ZERODEGC) THETAW1A.254
43 CONTINUE THETAW1A.255
41 CONTINUE THETAW1A.256
THETAW1A.257
9999 CONTINUE THETAW1A.258
RETURN THETAW1A.259
END THETAW1A.260
C THETAW1A.261
C THETAW1A.262
*ENDIF THETAW1A.263