*IF DEF,C90_1A,OR,DEF,C90_2A,OR,DEF,C90_2B AAD2F404.298
C ******************************COPYRIGHT****************************** GTS2F400.10711
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.10712
C GTS2F400.10713
C Use, duplication or disclosure of this code is subject to the GTS2F400.10714
C restrictions as set forth in the contract. GTS2F400.10715
C GTS2F400.10716
C Meteorological Office GTS2F400.10717
C London Road GTS2F400.10718
C BRACKNELL GTS2F400.10719
C Berkshire UK GTS2F400.10720
C RG12 2SZ GTS2F400.10721
C GTS2F400.10722
C If no contract has been raised with this copy of the code, the use, GTS2F400.10723
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.10724
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.10725
C Modelling at the above address. GTS2F400.10726
C ******************************COPYRIGHT****************************** GTS2F400.10727
C GTS2F400.10728
CLL Subroutine TWBULB-------------------------------------------------- TWBULB1A.3
CLL TWBULB1A.4
CLL Purpose: To calculate from temperatures,humidities the wet bulb TWBULB1A.5
CLL potential temperature on model levels and output the wet bulb TWBULB1A.6
CCL freezing level height. TWBULB1A.7
CLL TWBULB1A.8
CLL TWBULB1A.9
CLL Model Modification history from model version 3.0: TWBULB1A.10
CLL version Date TWBULB1A.11
CLL 3.2 10/06/94 Written by S.A.Woltering TWBULB1A.12
CLL 4.0 08/09/95 Speed up by vectorising. RTHBarnes. ARB1F400.1
!LL 4.3 26/02/97 Add first & last points to arg.list. RTHBarnes. ARB2F403.137
CLL TWBULB1A.13
CLL Programming standard : UM Doc Paper no 3 TWBULB1A.14
CLL TWBULB1A.15
CLL External documentation : UMDP no TWBULB1A.16
CLL TWBULB1A.17
CLLEND ----------------------------------------------------------------- TWBULB1A.18
C TWBULB1A.19
C*L ARGUMENTS:---------------------------------------------------------- TWBULB1A.20
SUBROUTINE TWBULB 1,1TWBULB1A.21
1 (Q,PSTAR,T, ! IN Model field array. TWBULB1A.22
2 AK,BK, ! IN AK/BK array. TWBULB1A.23
3 P_FIELD,P_LEVELS,Q_LEVELS, ! IN Field scalars. TWBULB1A.24
4 TW ! OUT WBT array. TWBULB1A.25
5 ,FIRST_POINT,LAST_POINT) ! IN loop start & end ARB2F403.138
C TWBULB1A.27
IMPLICIT NONE TWBULB1A.28
C TWBULB1A.29
INTEGER TWBULB1A.30
* P_FIELD ! IN No of points on a field. TWBULB1A.31
*,P_LEVELS ! IN No of model levels. TWBULB1A.32
*,Q_LEVELS ! IN No of humidity levels. TWBULB1A.33
*,FIRST_POINT,LAST_POINT ! IN 1st & last pts for calc ARB2F403.139
C TWBULB1A.34
REAL TWBULB1A.35
* T(P_FIELD,P_LEVELS) ! IN Intial temperature at all points. TWBULB1A.36
*,PSTAR(P_FIELD) ! IN Surface pressure. TWBULB1A.37
*,Q(P_FIELD,Q_LEVELS) ! IN Specific humidity at full levels. TWBULB1A.38
*,AK(P_LEVELS) ! IN Value of "A" at model level. TWBULB1A.39
*,BK(P_LEVELS) ! IN Value of "B" at model level. TWBULB1A.40
C TWBULB1A.41
REAL TWBULB1A.42
* TW(P_FIELD,Q_LEVELS) ! OUT The WET BULB temperature at all TWBULB1A.43
C levels and points. TWBULB1A.44
C DEFINE LOCAL WORKSPACE ARRAYS----------------------------------------- TWBULB1A.45
REAL TWBULB1A.47
* P(P_FIELD) ! Pressure at each point. TWBULB1A.48
*,LATENT_HEAT(P_FIELD) ! Latent heat of evaporation (fn(T)). TWBULB1A.49
*,GG(P_FIELD) ! The 'G' used in equation (1.3). TWBULB1A.50
*,Q_SAT(P_FIELD) ! Saturation specific humidity ARB1F400.2
*,DIFF(P_FIELD) ! The difference between G(Tw) & G(Ti). ARB1F400.3
C*---------------------------------------------------------------------- TWBULB1A.59
C*L EXTERNAL SUBROUTINES CALLED----------------------------------------- TWBULB1A.60
EXTERNAL QSAT TWBULB1A.61
C*---------------------------------------------------------------------- TWBULB1A.62
C CALL COMDECKS. TWBULB1A.63
C----------------------------------------------------------------------- TWBULB1A.64
*CALL C_R_CP 
TWBULB1A.65
*CALL C_G TWBULB1A.66
*CALL C_0_DG_C TWBULB1A.67
*CALL C_LHEAT TWBULB1A.68
*CALL C_EPSLON TWBULB1A.69
C TWBULB1A.71
C----------------------------------------------------------------------- TWBULB1A.72
C DEFINE LOCAL SCALAR VARIABLES. TWBULB1A.73
C----------------------------------------------------------------------- TWBULB1A.74
INTEGER TWBULB1A.75
* I,K,L ! Loop counters. TWBULB1A.76
*,LOOP TWBULB1A.77
C TWBULB1A.78
REAL TWBULB1A.79
* COEFF ! Coeff used in latent heat calculation. TWBULB1A.80
*,CPV ! Specific heat for water vapour. TWBULB1A.81
*,MV ! Mol wt of water vapour. TWBULB1A.84
*,RSTAR ! Universal gas constant. TWBULB1A.85
*,GS ! The 'G' used in equation (1.3). TWBULB1A.86
*,DGBYDT ! The function DG/DT Eqn (1.6). TWBULB1A.87
*,TEMP1 ! TWBULB1A.88
*,TEMP2 ! TWBULB1A.89
C----------------------------------------------------------------------- TWBULB1A.90
C DEFINE PARAMETER STATEMENTS TWBULB1A.91
C COEFF - Used in the calculation of LATENT_H. TWBULB1A.92
C CPV - Specific heat of water vapour. TWBULB1A.93
C MV - Mol wt of water vapour KG/MOL. TWBULB1A.94
C RSTAR - Universal gas constant. TWBULB1A.95
C----------------------------------------------------------------------- TWBULB1A.96
PARAMETER (COEFF=2.34E3, CPV=1850.0, MV=0.01801, RSTAR=8.314) TWBULB1A.97
C----------------------------------------------------------------------- TWBULB1A.98
C Begin by looping over the wet-levels. TWBULB1A.99
C----------------------------------------------------------------------- TWBULB1A.100
DO L=1,Q_LEVELS ! Loop over all wet-levels. TWBULB1A.101
C----------------------------------------------------------------------- TWBULB1A.102
C Calculate pressure for all points on that wet-level. ARB1F400.4
C----------------------------------------------------------------------- TWBULB1A.104
DO I=1,P_FIELD ! Loop over points. TWBULB1A.105
P(I) = AK(L) + BK(L)*PSTAR(I) TWBULB1A.106
CL-------------------- SECTION 1 --------------------------------------- TWBULB1A.107
CL Calculate the function G for TA and QA (see eqn 3) doc no TWBULB1A.108
C G(Tw)=Qa(Ta)*L(Ta)+Ta(Cp+QaCpv) TWBULB1A.109
C Subscript a indicates initial values. TWBULB1A.110
C----------------------------------------------------------------------- TWBULB1A.111
LATENT_HEAT(I)=LC-COEFF*(T(I,L)-ZERODEGC) TWBULB1A.112
GG(I)=Q(I,L)*LATENT_HEAT(I)+T(I,L)*(CP+Q(I,L)*CPV) ! Eqn 1.2 TWBULB1A.113
C----------------------------------------------------------------------- TWBULB1A.114
TW(I,L)=T(I,L) ! Initialise TW. TWBULB1A.115
ENDDO ! End of points loop. TWBULB1A.116
CL-------------------- SECTION 2 --------------------------------------- TWBULB1A.117
CL Calculate the function DG/DT TWBULB1A.118
CL G'(T)=(Mv*L(Ta)**2*Qs(T)/R*T**2) +Cp+Qa*Cpv TWBULB1A.119
CL Iterate to find the Tw TWBULB1A.120
CL T(i+1)=T(i)+(G(Tw)-G(Ti))/DG/DT(i) TWBULB1A.121
C----------------------------------------------------------------------- TWBULB1A.122
CL---- Set loop counter. ARB1F400.5
LOOP=1 ARB1F400.6
1000 CONTINUE ARB1F400.7
C ARB1F400.8
CALL QSAT(Q_SAT,TW(1,L),P(1),P_FIELD) ARB1F400.9
C ARB1F400.10
DO I=FIRST_POINT,LAST_POINT ARB2F403.140
GS=Q_SAT(I)*LATENT_HEAT(I)+TW(I,L)*(CP+Q(I,L)*CPV) ! Eqn 1.3 ARB1F400.11
ARB1F400.12
CL--------------------- SECTION 2.1 ------------------------------------ TWBULB1A.132
CL Find the difference between G(Tw)-G(Ti) TWBULB1A.133
C----------------------------------------------------------------------- TWBULB1A.134
DIFF(I) = ABS(GG(I)-GS) ! Eqn 1.8 ARB1F400.13
IF (DIFF(I).GT.1.0) THEN ARB1F400.14
ARB1F400.15
CL--------------------- SECTION 2.2 ------------------------------------ TWBULB1A.138
CL Calculate the function DG/DT TWBULB1A.139
C----------------------------------------------------------------------- TWBULB1A.140
TEMP1 = RSTAR*TW(I,L)*TW(I,L) ARB1F400.16
TEMP2 = Q_SAT(I)*MV*LATENT_HEAT(I)*LATENT_HEAT(I) ARB1F400.17
DGBYDT = TEMP2/TEMP1 + CP + CPV*Q(I,L) ! Eqn 1.6 TWBULB1A.146
ARB1F400.18
CL--------------------- SECTION 2.3 ------------------------------------ TWBULB1A.151
CL Using the function DG/DT calculate an updated Temperature Tw TWBULB1A.152
C----------------------------------------------------------------------- TWBULB1A.153
TW(I,L) = TW(I,L) - (GS-GG(I)) / DGBYDT ! Eqn 1.7 TWBULB1A.154
ARB1F400.19
CL--------------------- SECTION 2.4 ------------------------------------ TWBULB1A.156
CL Using the new temperature Tw re-calculate GS First update LATENT_H TWBULB1A.157
C----------------------------------------------------------------------- TWBULB1A.158
LATENT_HEAT(I)=LC-COEFF*(TW(I,L)-ZERODEGC) TWBULB1A.159
ENDIF TWBULB1A.160
END DO ! I ARB1F400.20
ARB1F400.21
CL----Increment iteration loop counter. TWBULB1A.161
LOOP=LOOP+1 TWBULB1A.162
CL----Test for convergence. TWBULB1A.163
IF(LOOP.GT.10) THEN TWBULB1A.164
WRITE(6,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>' TWBULB1A.166
WRITE(6,*)'>>> TWBULB - Convergence failure, level ',L ARB1F400.22
WRITE(6,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>' TWBULB1A.169
GOTO 9999 TWBULB1A.170
ENDIF TWBULB1A.171
CL----Difference test. TWBULB1A.172
DO I=FIRST_POINT,LAST_POINT ARB2F403.141
IF (DIFF(I).GT.1.0) GOTO 1000 ARB1F400.24
ENDDO ! I ARB1F400.25
9999 CONTINUE ARB1F400.26
! Set points outside calculation range to sensible values. ARB2F403.142
DO I=1,FIRST_POINT-1 ARB2F403.143
TW(I,L) = TW(FIRST_POINT,L) ARB2F403.144
END DO ARB2F403.145
DO I=LAST_POINT+1,P_FIELD ARB2F403.146
TW(I,L) = TW(LAST_POINT,L) ARB2F403.147
END DO ARB2F403.148
ENDDO ! End of loop over levels. TWBULB1A.176
RETURN TWBULB1A.177
END TWBULB1A.178
*ENDIF TWBULB1A.179