*IF DEF,A02_1B LWPLAN1B.2
C ******************************COPYRIGHT****************************** GTS2F400.5617
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.5618
C GTS2F400.5619
C Use, duplication or disclosure of this code is subject to the GTS2F400.5620
C restrictions as set forth in the contract. GTS2F400.5621
C GTS2F400.5622
C Meteorological Office GTS2F400.5623
C London Road GTS2F400.5624
C BRACKNELL GTS2F400.5625
C Berkshire UK GTS2F400.5626
C RG12 2SZ GTS2F400.5627
C GTS2F400.5628
C If no contract has been raised with this copy of the code, the use, GTS2F400.5629
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.5630
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.5631
C Modelling at the above address. GTS2F400.5632
C ******************************COPYRIGHT****************************** GTS2F400.5633
C GTS2F400.5634
CLL Subroutine LWPLAN ---------------------------------------------- LWPLAN1B.3
CLL LWPLAN1B.4
CLL Purpose : LWPLAN1B.5
CLL The version of routine LWPLAN used in LWPLAN1B.6
CLL version 1B (gaseous effects treated as Morcrette et al, 1986) of LWPLAN1B.7
CLL the UM LW code. LWPLAN1B.8
CLL Version 3 using the ECMWF quintic fits for the ECMWF LW bands, part LWPLAN1B.9
CLL of the alternative code giving ECMWF-like treatment of LW gaseous LWPLAN1B.10
CLL transmissivities. The changes are not great: temperatures are now LWPLAN1B.11
CLL normalized to a reference temperature before using the fits, the LWPLAN1B.12
CLL fits are quintic rather than cubic, the constants (in *COMDECK LWPLAN1B.13
CLL LWPFCO) are different, and the indentation is changed. LWPLAN1B.14
CLL Version 3 of LWPLAN was set up from version 2.2 to be part of LWPLAN1B.15
CLL version 1B (ECMWF-like gaseous transmissivities) of the LW from LWPLAN1B.16
CLL release 2.7 of the UM. William Ingram 22 June 1992 LWPLAN1B.17
CLL It calculates Planck ("black-body") fluxes in each longwave band as LWPLAN1B.18
CLL a quintic function of atmospheric (layer centre and boundary) and LWPLAN1B.19
CLL surface temperatures and returns their differences across LWPLAN1B.20
CLL half-layers (and optionally the surface black-body flux) for use by WI200893.32
CLL by LWMAST constructing longwave fluxes. WI200893.33
CLL LWPLAN1B.23
CLL Author: William Ingram LWPLAN1B.24
CLL LWPLAN1B.25
CLL Model Modification history from model version 3.0: LWPLAN1B.26
CLL version Date LWPLAN1B.27
CLL LWPLAN1B.28
CLL Programming standard : LWPLAN1B.29
CLL It conforms with standard A of version 3 (07/9/90) of UMDP 4, and LWPLAN1B.30
CLL contains no 8X-deprecated features. LWPLAN1B.31
CLL The code is standard FORTRAN 77 except for having ! comments. LWPLAN1B.32
CLL LWPLAN1B.33
CLL Logical components covered : P232, D23 LWPLAN1B.34
CLL P232 (longwave radiation) D23 (radiation diagnostics). LWPLAN1B.35
CLL LWPLAN1B.36
CLL Project task : P23 LWPLAN1B.37
CLL LWPLAN1B.38
CLL External documentation: UMDP 23. LWPLAN1B.39
CLL LWPLAN1B.40
CLLEND ----------------------------------------------------------------- LWPLAN1B.41
C*L LWPLAN1B.42
SUBROUTINE LWPLAN (TAC, PEXNER, PSTAR, AKH, BKH, TSTAR, AICE, 2LWPLAN1B.43
& SFUP, SFUPON, LWPLAN1B.44
& L2, NLEVS, L1, SEAFX, BHDB, THDB) LWPLAN1B.45
C* LWPLAN1B.46
*CALL LWNBANDS 
LWPLAN1B.47
C*L LWPLAN1B.48
INTEGER!, INTENT (IN) :: LWPLAN1B.49
& L1, ! First dimension of input arrays LWPLAN1B.50
& L2, ! Number of points to be treated LWPLAN1B.51
& NLEVS ! Number of levels LWPLAN1B.52
REAL!, INTENT(IN) :: LWPLAN1B.53
& TAC(L1,NLEVS), ! Atmospheric temperatures at layer LWPLAN1B.54
C ! centres LWPLAN1B.55
& PEXNER(L1,NLEVS+1), ! Exner function @ layer boundaries LWPLAN1B.56
& PSTAR(L1), ! Surface pressure LWPLAN1B.57
& AKH(NLEVS+1), ! As at layer boundaries LWPLAN1B.58
& BKH(NLEVS+1), ! Bs at layer boundaries LWPLAN1B.59
C ! The above four are used to get temperature at layer boundaries LWPLAN1B.60
& TSTAR(L1), ! Surface temperature LWPLAN1B.61
& AICE(L1) ! Sea-ice fraction LWPLAN1B.62
LOGICAL!, INTENT(IN) :: LWPLAN1B.63
& SFUPON LWPLAN1B.64
REAL!, INTENT(OUT) :: LWPLAN1B.65
& BHDB(L2,NLEVS,NBANDS), LWPLAN1B.66
& THDB(L2,NLEVS,NBANDS), LWPLAN1B.67
C ! BHDB holds the difference of the Planck functions across the LWPLAN1B.68
C ! bottom half of each layer and THDB across the top half. LWPLAN1B.69
& SEAFX(L1) ! Difference between the open-sea LWPLAN1B.70
C ! and sea-ice upward longwave fluxes at the surface. LWPLAN1B.71
& , SFUP(L1) ! Upward surface flux LWPLAN1B.72
C* LWPLAN1B.73
CL ! LWPLAN has no EXTERNAL calls and no dynamically allocated arrays LWPLAN1B.74
C LWPLAN1B.75
*CALL C_0_DG_C LWPLAN1B.76
*CALL C_R_CP LWPLAN1B.77
REAL TRPLAN ! Reference temperature to which LWPLAN1B.78
C ! temperatures are normalized before using the quintic fits. LWPLAN1B.79
PARAMETER ( TRPLAN = 250. ) LWPLAN1B.80
*CALL LWPFCO LWPLAN1B.81
C LWPLAN1B.82
INTEGER BAND, LEVEL, J ! Loopers over band, level & point LWPLAN1B.83
REAL TAB, ! Layer boundary temperature LWPLAN1B.84
& WTL, WTU, ! and weights for its calculation LWPLAN1B.85
& PL, ! Pressure at upper bdry of current LWPLAN1B.86
& PU, ! layer and next one up LWPLAN1B.87
& PLM1, ! Pressure at lower bdry of layer LWPLAN1B.88
& BBC, ! Black-body flux at a model layer LWPLAN1B.89
& BBB, ! centre or boundary LWPLAN1B.90
& BBTFS ! or at TFS LWPLAN1B.91
REAL TN, ! Simplify the calculation of LWPLAN1B.92
& PLANCK, ! Planck fluxes with statement LWPLAN1B.93
& T ! functions LWPLAN1B.94
*CALL P_EXNERC LWPLAN1B.95
TN(T) = T / TRPLAN - 1. LWPLAN1B.96
PLANCK(T,BAND) = LWPLAN1B.97
& ( ( ( ( PFCO(6,BAND) * TN(T) + PFCO(5,BAND) ) * TN(T) + LWPLAN1B.98
& PFCO(4,BAND) ) * TN(T) + PFCO(3,BAND) ) * TN(T) + LWPLAN1B.99
& PFCO(2,BAND) ) * TN(T) + PFCO(1,BAND) LWPLAN1B.100
C LWPLAN1B.101
CL ! Section 1 LWPLAN1B.102
CL ! ~~~~~~~~~ LWPLAN1B.103
CL ! First find Planck function for surface temperature: LWPLAN1B.104
C LWPLAN1B.105
DO 1 BAND=1, NBANDS LWPLAN1B.106
BBTFS = PLANCK(TFS,BAND) LWPLAN1B.107
Cfpp$ Select(CONCUR) LWPLAN1B.108
DO 11 J=1, L2 LWPLAN1B.109
IF ( AICE(J) .GT. 0. ) THEN LWPLAN1B.110
TAB = ( TSTAR(J) + (AICE(J)-1.) * TFS ) / AICE(J) LWPLAN1B.111
ELSE LWPLAN1B.112
TAB = TSTAR(J) LWPLAN1B.113
ENDIF LWPLAN1B.114
BBB = PLANCK(TAB,BAND) LWPLAN1B.115
SEAFX(J) = SEAFX(J) + BBTFS - BBB LWPLAN1B.116
IF ( AICE(J) .GT. 0. ) LWPLAN1B.117
& BBB = AICE(J) * BBB + (1.-AICE(J)) * BBTFS LWPLAN1B.118
BHDB(J,1,BAND) = BBB LWPLAN1B.119
IF ( SFUPON ) SFUP(J) = SFUP(J) + BBB LWPLAN1B.120
11 CONTINUE LWPLAN1B.121
1 CONTINUE LWPLAN1B.122
C LWPLAN1B.123
CL ! Section 2 LWPLAN1B.124
CL ! ~~~~~~~~~ LWPLAN1B.125
CL ! Loop over the model layers finding the dB across each half LWPLAN1B.126
CL ! of each one : LWPLAN1B.127
C LWPLAN1B.128
DO 2 LEVEL=1, NLEVS-1 LWPLAN1B.129
C ! LWPLAN1B.130
C ! First must find TAB for the top of the layer (same for each LWPLAN1B.131
C ! band) and store it in spare space in THDB: LWPLAN1B.132
Cfpp$ Select(CONCUR) LWPLAN1B.133
DO 20 J=1, L2 LWPLAN1B.134
PU = PSTAR(J)*BKH(LEVEL+2) + AKH(LEVEL+2) LWPLAN1B.135
PL = PSTAR(J)*BKH(LEVEL+1) + AKH(LEVEL+1) LWPLAN1B.136
PLM1 = PSTAR(J)*BKH(LEVEL) + AKH(LEVEL) LWPLAN1B.137
WTL = TAC(J,LEVEL+1)*( PEXNER(J,LEVEL+1) / LWPLAN1B.138
& P_EXNER_C(PEXNER(J,LEVEL+2),PEXNER(J,LEVEL+1),PU,PL,KAPPA) LWPLAN1B.139
& - 1.0 ) LWPLAN1B.140
WTU = TAC(J,LEVEL) * ( PEXNER(J,LEVEL) / LWPLAN1B.141
& P_EXNER_C(PEXNER(J,LEVEL+1),PEXNER(J,LEVEL),PL,PLM1,KAPPA) LWPLAN1B.142
& - 1.0 ) LWPLAN1B.143
TAB = LWPLAN1B.144
& ( WTL * TAC(J,LEVEL) + WTU * TAC(J,LEVEL+1) ) / ( WTL + WTU ) LWPLAN1B.145
THDB(J,NLEVS,1) = TAB LWPLAN1B.146
20 CONTINUE LWPLAN1B.147
C ! LWPLAN1B.148
C ! This loop finds the black-body fluxes at the middle & top of LWPLAN1B.149
C ! the layer and sets BHDB & THDB for the layer, using the LWPLAN1B.150
C ! black-body flux for the bottom of the layer which is already LWPLAN1B.151
C ! in BHDB and leaving the value at the top in the next level LWPLAN1B.152
C ! up of BHDB : LWPLAN1B.153
DO 21 BAND=1, NBANDS LWPLAN1B.154
Cfpp$ Select(CONCUR) LWPLAN1B.155
DO 22 J=1, L2 LWPLAN1B.156
BBC = PLANCK(TAC(J,LEVEL),BAND) LWPLAN1B.157
BHDB(J,LEVEL,BAND) = BBC - BHDB(J,LEVEL,BAND) LWPLAN1B.158
TAB = THDB(J,NLEVS,1) LWPLAN1B.159
BBB = PLANCK(TAB,BAND) LWPLAN1B.160
THDB(J,LEVEL,BAND) = BBB - BBC LWPLAN1B.161
BHDB(J,LEVEL+1,BAND) = BBB LWPLAN1B.162
22 CONTINUE LWPLAN1B.163
21 CONTINUE LWPLAN1B.164
2 CONTINUE LWPLAN1B.165
C LWPLAN1B.166
CL ! Section 3 LWPLAN1B.167
CL ! ~~~~~~~~~ LWPLAN1B.168
CL ! Finally, the top layer is treated slightly differently because LWPLAN1B.169
CL ! the black-body flux at the top is zero: LWPLAN1B.170
C LWPLAN1B.171
DO 3 BAND=1, NBANDS LWPLAN1B.172
Cfpp$ Select(CONCUR) LWPLAN1B.173
DO 30 J=1, L2 LWPLAN1B.174
BBC = PLANCK(TAC(J,NLEVS),BAND) LWPLAN1B.175
BHDB(J,NLEVS,BAND) = BBC - BHDB(J,NLEVS,BAND) LWPLAN1B.176
THDB(J,NLEVS,BAND) = - BBC LWPLAN1B.177
30 CONTINUE LWPLAN1B.178
3 CONTINUE LWPLAN1B.179
C LWPLAN1B.180
RETURN LWPLAN1B.181
END LWPLAN1B.182
*ENDIF A02_1B LWPLAN1B.183