*IF DEF,A16_1A PHYDIA1A.2
C ******************************COPYRIGHT****************************** GTS2F400.7219
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.7220
C GTS2F400.7221
C Use, duplication or disclosure of this code is subject to the GTS2F400.7222
C restrictions as set forth in the contract. GTS2F400.7223
C GTS2F400.7224
C Meteorological Office GTS2F400.7225
C London Road GTS2F400.7226
C BRACKNELL GTS2F400.7227
C Berkshire UK GTS2F400.7228
C RG12 2SZ GTS2F400.7229
C GTS2F400.7230
C If no contract has been raised with this copy of the code, the use, GTS2F400.7231
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.7232
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.7233
C Modelling at the above address. GTS2F400.7234
C ******************************COPYRIGHT****************************** GTS2F400.7235
C GTS2F400.7236
CLL SUBROUTINE PHY_DIAG------------------------------------------------ PHYDIA1A.3
CLL PHYDIA1A.4
CLL PURPOSE: Calculation of physical diagnostics. These diagnostics PHYDIA1A.5
CLL are those based on temperatures humidities and pressure. PHYDIA1A.6
CLL Those diagnostics involving the winds are in the other PHYDIA1A.7
CLL diagnostic routine DYN_DIAG PHYDIA1A.8
CLL Date: 08/06/92 PHYDIA1A.9
CLL Exner at model level made consistent with the geopotential eqn. PHYDIA1A.10
CLL PHYDIA1A.11
CLL Model Modification history from model version 3.0: PHYDIA1A.12
CLL version Date PHYDIA1A.13
CLL 3.1 10/01/93 Get model half Hts for all dependant routines. PS100193.1
CLL 3.1 09/02/93 CW090293.1
CLL Modified by C.Wilson CW090293.2
CLL Add alternate call to V_INT_T for temperatures on pressure levels CW090293.3
CLL to call V_INT_TP if QV_INT_TP is true CW090293.4
CLL so as to reduce bias in stratospheric temperature output CW090293.5
CLL 3.2 13/07/93 Changed CHARACTER*(*) to CHARACTER*(80) for TS150793.108
CLL portability. Author Tracey Smith. TS150793.109
CLL 3.2 07/06/93 Correct indexing of WORK1 for use in DIA_THADV RB070693.1
CLL 3.4 21/9/94 Calculate model level geopotential heights ASW3F304.24
CLL (in m). Author S.A.Woltering. ASW3F304.25
CLL 3.4 15/07/94 Redimension WORK1 from P_TO_UV for input to GRR0F304.15
CLL DIA_THADV using new utility routine. Rick Rawlins. GRR0F304.16
CLL 4.0 22/11/95 Correction to Model level geopotential heights. ARB1F400.27
CLL Author Steve Woltering ARB1F400.28
CLL 4.1 15/05/96 Put tracer fields on pressure levels. D.Podd ADP0F401.7
CLL 4.4 16/10/97 Extra swapbounds for P_EXNER_HALF. D Robinson. GDR2F404.24
!LL 4.5 15/04/98 Added start-end arguments to V_INT, V_INT_T and GSM1F405.748
!LL V_INT_Z routines. S.D.Mullerworth GSM1F405.749
CLL 4.5 3/09/98 Corrected relative humidity diagnostic. ADM0F405.297
CLL D Wilson. ADM0F405.298
CLL 4.5 05/06/98 New argument L_VINT_TP. Use to call V_INT_TP or GDR4F405.15
CLL V_INT_T. Also L_LSPICE for correcting relative GDR4F405.16
CLL humidity diagnostic. D. Robinson. GDR4F405.17
CLL PHYDIA1A.14
CLL CW090293.6
CLL Programming standard: U M DOC Paper NO. 4, PHYDIA1A.15
CLL PHYDIA1A.16
CLL SYSTEM COMPONENTS D471 D421 D481 D482 D432 D423 D422 D431 D472 PHYDIA1A.17
CLL D435 D484 D441 PHYDIA1A.18
CLL PHYDIA1A.19
CLL SYSTEM TASK: D4 PHYDIA1A.20
CLL PHYDIA1A.21
CLL External documentation PHYDIA1A.22
CLL PHYDIA1A.23
CLLEND------------------------------------------------------------------ PHYDIA1A.24
C PHYDIA1A.25
C*L ARGUMENTS:--------------------------------------------------------- PHYDIA1A.26
SUBROUTINE PHY_DIAG( 2,23PHYDIA1A.27
*CALL ARGFLDPT 
GSM1F405.750
C primary data in PHYDIA1A.28
& PSTAR,U,V,Q,THETA,OROG,P_EXNER_HALF,LAND,TSTAR,TRACERS, ADP0F401.8
C primary data constants PHYDIA1A.30
& U_ROWS,P_ROWS,ROW_LENGTH,P_LEVELS,Q_LEVELS,P_FIELD, PHYDIA1A.31
& U_FIELD,AK,BK,AKH,BKH,EW_SPACE,NS_SPACE,SEC_U_LATITUDE, PHYDIA1A.32
& TR_LEVELS,TR_VARS,TR_P_FIELD_DA, ADP0F401.9
C Input pressure for required interpolation PHYDIA1A.33
& T_P_PRESS,HTS_PRESS,REL_HUMID_PRESS,WBPT_PRESS,TH_ADV_PRESS, PHYDIA1A.34
& TR_PRESS, ADP0F401.10
C Input pressure for indices for product field PHYDIA1A.35
& H2_IND, PHYDIA1A.36
C DIAGNOSTICS OUT PHYDIA1A.37
& HEIGHTS, T_P, REL_HUMID_P, WBPT, SNPROB, MINUS20_ICAO, PHYDIA1A.38
& MINUS20_P, FREEZE_ICAO, FREEZE_Z, FREEZE_P, CONTRAIL_LOWERHT, PHYDIA1A.39
& CONTRAIL_UPPERHT, TROP_P, TROP_T, TROP_Z, TROP_ICAO, PHYDIA1A.40
& MINUS20_Z, TH_ADV_SINGLE, DUCT_HEIGHT, DUCT_INT, PHYDIA1A.41
& P_MSL, TH_ADV_MEAN, HEIGHTS2, MODHT, ASW3F304.26
& TRACERWORK,INT16,PT_TRACER, ADP0F401.11
C diagnostic lengths PHYDIA1A.43
& T_P_LEVS,HTS_LEVS,REL_HUMID_P_LEVS,WBPT_LEVS,TH_ADV_P_LEVS, PHYDIA1A.44
& H2_P_LEVS, PHYDIA1A.45
& TR_PRESS_LEVS,NUM_STASH_LEVELSDA, ADP0F401.12
C diagnostic logical indicators PHYDIA1A.46
& QMODEL_HALF_HEIGHT,QHEIGHTS, QT_P, QREL_HUMID_P, PHYDIA1A.47
& QWBPT, QSNPROB, QMINUS20_ICAO, QMINUS20_P, PHYDIA1A.48
& QFREEZE_ICAO, QFREEZE_Z, QFREEZE_P, QCONTRA_LOWERHT, PHYDIA1A.49
& QCONTRA_UPPERHT, QTROP_P, QTROP_T, QTROP_Z, PHYDIA1A.50
& QTROP_ICAO, QMINUS20_Z, QTH_ADV_SINGLE, QDUCT_HEIGHT, PHYDIA1A.51
& QDUCT_INT, QPMSL, QTH_ADV_MEAN, QHTS2, QMODHT, ASW3F304.27
& SF_TRACER, ADP0F401.13
& Z_REF, L_VINT_TP, L_LSPICE, GDR4F405.18
& ICODE,CMESSAGE) GDR4F405.19
C PHYDIA1A.55
IMPLICIT NONE PHYDIA1A.56
*CALL C_EPSLON ASW3F304.28
*CALL C_R_CP PHYDIA1A.57
*CALL C_G PHYDIA1A.58
*CALL TYPFLDPT GSM1F405.751
PHYDIA1A.59
INTEGER PHYDIA1A.60
* P_FIELD !IN 1ST DIMENSION OF FIELD OF PSTAR PHYDIA1A.61
*, U_FIELD !IN 1ST DIMENSION OF FIELD OF U,V PHYDIA1A.62
*, U_ROWS !IN NUMBER OF ROWS FOR U,V FIELD PHYDIA1A.63
*, P_ROWS !IN NUMBER OF ROWS FOR P,T FIELD PHYDIA1A.64
*, ROW_LENGTH !IN NUMBER OF POINTS PER ROW PHYDIA1A.65
*, LEVELS !IN NUMBER OF MODEL LEVELS PHYDIA1A.66
*, P_LEVELS !IN NUMBER OF PRESSURE LEVELS PHYDIA1A.67
*, Q_LEVELS !IN NUMBER OF WET LEVELS PHYDIA1A.68
*, ICODE !IN RETURN CODE : IRET=0 NORMAL EXIT PHYDIA1A.69
*, NUM_STASH_LEVELSDA !IN NUMBER OF STASH LEVELS ADP0F401.14
*, TR_LEVELS !IN NUMBER OF TRACER LEVELS ADP0F401.15
*, TR_VARS !IN NUMBER OF TRACERS ADP0F401.16
*, TR_P_FIELD_DA !IN P_FIELD for DA of tracer arrays ADP0F401.17
*, INT16 !IN Dummy variable for STASH_MAXLEN(16) ADP0F401.18
*, PT_TRACER(TR_VARS+1) ADP0F401.19
* !IN POINTERS TO TRACERS IN STASHWORK ADP0F401.20
*, TR_PRESS_LEVS(TR_VARS+1) ADP0F401.21
* !IN NO OF LEVS ON WHICH TO INTERP TRACERS ADP0F401.22
*, T_P_LEVS !IN NO OF LEVS ON WHICH TO INTERP T_P PHYDIA1A.70
*, REL_HUMID_P_LEVS !IN NO OF LEVS ON WHICH TO INTERP REL HUM PHYDIA1A.71
*, HTS_LEVS !IN NO OF LEVS ON WHICH TO INTERP HEIGHTS PHYDIA1A.72
*, H2_P_LEVS !IN NO OF LEVS ON WHICH TO INTERP H**2 PHYDIA1A.73
*, WBPT_LEVS !IN NO OF LEVS ON WHICH TO CALCULATE WBPT PHYDIA1A.74
*, TH_ADV_P_LEVS !IN NO OF LEVS ON WHICH TO CALCULATE THADV PHYDIA1A.75
*, Z_REF !IN LEVEL OF MODEL USED TO CALCULATE PMSL PHYDIA1A.76
*, BOTTOM_QC_LEVEL ! Tropopause quality control level PHYDIA1A.77
*, TOP_QC_LEVEL ! " " " " PHYDIA1A.78
INTEGER PHYDIA1A.79
& H2_IND(H2_P_LEVS) !IN indices for pressure levels for h**2 PHYDIA1A.80
C PHYDIA1A.81
C PHYDIA1A.82
LOGICAL PHYDIA1A.83
* QT_P !IN LOGICAL FLAG FOR PRESS INTER TEMPERTATURE PHYDIA1A.84
*,QHEIGHTS !IN " " " HEIGHTS ON ANY P SURFACE PHYDIA1A.85
*,QREL_HUMID_P !IN " " " HUMIDITIES ANY P SURFACE PHYDIA1A.86
*,QWBPT !IN " " " WET BULB POTENTIAL TEMP PHYDIA1A.87
*,QSNPROB !IN " " " SNOW PROBABILITY PHYDIA1A.88
*,QMINUS20_ICAO !IN " " " -20 deg C LEVEL ICAO HT PHYDIA1A.89
*,QMINUS20_P !IN " " " -20 deg C LEVEL PRESSURE PHYDIA1A.90
*,QFREEZE_ICAO !IN " " " FREEZING LEVEL ICAO HT PHYDIA1A.91
*,QFREEZE_Z !IN " " " FREEZING LEVEL HT PHYDIA1A.92
*,QFREEZE_P !IN " " " FREEZING LEVEL PRESSURE PHYDIA1A.93
*,QCONTRA_LOWERHT!IN " " " HEIGHT OF LOWER CONTRAIL PHYDIA1A.94
*,QCONTRA_UPPERHT!IN " " " HEIGHT OF UPPER CONTRAIL PHYDIA1A.95
*,QMODEL_HALF_HEIGHT!IN " " " HEIGHTS ON MODEL LAYER PHYDIA1A.96
*,LAND(P_FIELD) !IN LAND SEA MASK SET TRUE FOR LAND PHYDIA1A.97
LOGICAL PHYDIA1A.98
* QTROP_P !IN " " " TROPOPAUSE PRESSURE PHYDIA1A.99
*,QTROP_T !IN " " " TROPOPAUSE TEMPERATURE PHYDIA1A.100
*,QTROP_Z !IN " " " TROPOPAUSE HEIGHT PHYDIA1A.101
*,QTROP_ICAO !IN " " " TROPOPAUSE ICAO HEIGHT PHYDIA1A.102
*,QMINUS20_Z !IN " " " -20 deg C LEVEL HEIGHT PHYDIA1A.103
*,QTH_ADV_SINGLE !IN " " " THERMAL ADV SINGLE LEV PHYDIA1A.104
*,QDUCT_HEIGHT !IN " " " RADIO DUCT HEIGHT PHYDIA1A.105
*,QDUCT_INT !IN " " " RADIO DUCT INTENSITY PHYDIA1A.106
*,QPMSL !IN " " " PRESSURE AT MEAN SEA LVL PHYDIA1A.107
*,QTH_ADV_MEAN !IN " " " THERMAL ADVECTION MEAN PHYDIA1A.108
*,QHTS2 !IN " " " heights**2 on pressure lev PHYDIA1A.109
*,QMODHT !IN " " " Model level geopot heights ASW3F304.29
*,SF_TRACER(TR_VARS+1) !IN LOGICAL FLAGS FOR TRACERS ADP0F401.23
&, L_VINT_TP !IN Switch to control Vert Int for Output of GDR4F405.20
! Temperature on model levels. D Robinson. GDR4F405.21
&, L_LSPICE !IN Switch for New cloud/precip microphysics. GDR4F405.22
C PHYDIA1A.110
CHARACTER*80 CMESSAGE TS150793.110
C PHYDIA1A.113
REAL PHYDIA1A.114
* PSTAR(P_FIELD) !IN PRIMARY MODEL ARRAY FOR PSTAR FIELD PHYDIA1A.115
*,TSTAR(P_FIELD) !IN PRIMARY MODEL ARRAY FOR TSTAR FIELD PHYDIA1A.116
*,OROG(P_FIELD) !IN PRIMARY MODEL OROGRAPHY PHYDIA1A.117
*,P_EXNER_HALF(P_FIELD,P_LEVELS+1) !IN EXNER PRESS ON 1/2 LVLS PHYDIA1A.118
*,THETA(P_FIELD,P_LEVELS) !IN PRIMARY MODEL ARRAY FOR THETA FIELD PHYDIA1A.119
*,U(U_FIELD,P_LEVELS) !INT PRIMARY MODEL ARRAY FOR U FIELD PHYDIA1A.120
*,V(U_FIELD,P_LEVELS) !IN PRIMARY MODEL ARRAY FOR V FIELD PHYDIA1A.121
*,Q(P_FIELD,Q_LEVELS) !IN PRIMARY MODEL ARRAY FOR HUMIDITY PHYDIA1A.122
*,TRACERS(TR_P_FIELD_DA,TR_LEVELS,TR_VARS+1) ADP0F401.24
* !IN PRIMARY MODEL ARRAY FOR TRACERS ADP0F401.25
*,HEIGHTS(P_FIELD,HTS_LEVS) ! OUTPUT HEIGHTS ON ANY P SURFACE PHYDIA1A.123
*,T_P(P_FIELD,T_P_LEVS) ! OUTPUT TEMPS ON ANY P SURFACE PHYDIA1A.124
*,P_MSL(P_FIELD) ! OUTPUT PMSL PHYDIA1A.125
*,TROP_T(P_FIELD) ! OUTPUT TEMPS OF TROPOPAUSE PHYDIA1A.126
*,TROP_P(P_FIELD) ! OUTPUT PRESSURE OF TROPOPAUSE PHYDIA1A.127
*,TROP_Z(P_FIELD) ! OUTPUT HEIGHT OF TROPOPAUSE PRESSURE SURFACE PHYDIA1A.128
*,TROP_ICAO(P_FIELD) ! OUTPUT ICAO HT OF TROPOPAUSE PRESSURE PHYDIA1A.129
*,HEIGHTS2(P_FIELD,H2_P_LEVS) ! OUTPUT height**2 on any P surface PHYDIA1A.130
*,MODHT(P_FIELD,P_LEVELS) ! OUTPUT GEOPOT. MODEL LEV HEIGHTS ASW3F304.30
C PHYDIA1A.131
REAL PHYDIA1A.132
* REL_HUMID_P(P_FIELD,REL_HUMID_P_LEVS) PHYDIA1A.133
* ! OUTPUT HUMIDITIES ON ANY P SURFACE PHYDIA1A.134
*,WBPT(P_FIELD,WBPT_LEVS) ! WET BULB POTTEMP ON ANY P SURFACE PHYDIA1A.135
*,SNPROB(P_FIELD) ! SNOW PROBABILITY PHYDIA1A.136
*,MINUS20_ICAO(P_FIELD) ! MINUS 20 LEVEL ICAO HEIGHT PHYDIA1A.137
*,MINUS20_P(P_FIELD) ! MINUS 20 LEVEL PRESSURE PHYDIA1A.138
*,FREEZE_Z(P_FIELD) ! HEIGHT OF THE FREEZING LEVEL PHYDIA1A.139
*,FREEZE_ICAO(P_FIELD) ! ICAO HEIGHT OF THE FREEZING LEVEL PHYDIA1A.140
*,FREEZE_P(P_FIELD) ! PRESSURE OF THE FREEZING LEVEL PHYDIA1A.141
*,CONTRAIL_UPPERHT(P_FIELD) ! UPPER VALUE OF THE CONTRAIL PHYDIA1A.142
*,CONTRAIL_LOWERHT(P_FIELD) ! LOWER VALUE OF THE CONTRAIL PHYDIA1A.143
*,MINUS20_Z(P_FIELD) ! MINUS 20 LEVEL TRUE HEIGHT PHYDIA1A.144
*,TH_ADV_SINGLE(P_FIELD,TH_ADV_P_LEVS)! THERMAL ADV SINGLE LEVELS PHYDIA1A.145
*,DUCT_HEIGHT(P_FIELD) ! RADIO DUCT HEIGHT PHYDIA1A.146
*,DUCT_INT(P_FIELD) ! RADIO DUCT INTENSITY PHYDIA1A.147
*,TH_ADV_MEAN(P_FIELD) ! THERMAL ADVECTION MEAN 850/700/500 PHYDIA1A.148
*,TRACER_P(TR_P_FIELD_DA) ! TRACER FIELD ON ANY P SURFACE ADP0F401.26
*,TRACERWORK(INT16) ! STASHWORK FOR SUBSTITUTION OF TRACERS ADP0F401.27
C---------------------------------------------------------------------- PHYDIA1A.149
C AK,BK DEFINE HYBRID VERTICAL COORDINATES P=A+BP*, PHYDIA1A.150
C DELTA_AK,DELTA_BK DEFINE LAYER PRESSURE THICKNESS PD=AD+BDP*, PHYDIA1A.151
C---------------------------------------------------------------------- PHYDIA1A.152
REAL PHYDIA1A.153
* AKH(P_LEVELS+1) !IN value at layer boundary PHYDIA1A.154
*,BKH(P_LEVELS+1) !IN value at layer boundary PHYDIA1A.155
*,AK (P_LEVELS) !IN VALUE AT LAYER CENTRE PHYDIA1A.156
*,BK (P_LEVELS) !IN VALUE AT LAYER CENTRE PHYDIA1A.157
*,EW_SPACE !IN LONGITUDE GRID SPACING PHYDIA1A.158
*,NS_SPACE !IN LATITUDE GRID SPACING PHYDIA1A.159
*,SEC_U_LATITUDE(U_FIELD) !IN 1/COS(LAT) AT U POINTS PHYDIA1A.160
*,TIMESTEP !IN TIMESTEP PHYDIA1A.161
*,T_P_PRESS(T_P_LEVS) !IN Pressures reqd for Temp inter PHYDIA1A.162
*,HTS_PRESS(HTS_LEVS) !IN " " " HTS PHYDIA1A.163
*,REL_HUMID_PRESS(REL_HUMID_P_LEVS)!IN " " " REL Humidity PHYDIA1A.164
*,WBPT_PRESS(WBPT_LEVS) !IN " " " W.B.P.T PHYDIA1A.165
*,TH_ADV_PRESS(TH_ADV_P_LEVS) !IN " " " Thermal adv PHYDIA1A.166
*,TR_PRESS(TR_VARS+1,NUM_STASH_LEVELSDA) ADP0F401.28
* !IN Pressures reqd for Tracer interp ADP0F401.29
C*---------------------------------------------------------------------- PHYDIA1A.167
PHYDIA1A.168
C*L WORKSPACE USAGE:--------------------------------------------------- PHYDIA1A.169
C REAL ARRAYS REQUIRED AT FULL FIELD LENGTH 4*(P_LEVELS)+7+Q_LEVELS PHYDIA1A.170
C----------------------------------------------------------------------- PHYDIA1A.171
REAL PHYDIA1A.172
* MODEL_HALF_HEIGHT(P_FIELD,P_LEVELS+1) !OUT HEIGHTS OF MODEL HALF PHYDIA1A.173
*, WORK1(P_FIELD,P_LEVELS+1) ! Will hold Rel Humid and P_HALF PHYDIA1A.174
* ! and temperatures for FREEZE PHYDIA1A.175
* ! and pressures for THADV PHYDIA1A.176
*,WORK2(P_FIELD) ! Workspace for Mean thermal advection PHYDIA1A.177
*, PHI_STAR(P_FIELD) ! Geopotential PHYDIA1A.178
*, P(P_FIELD,P_LEVELS) ! Pressure ARRAY PHYDIA1A.179
*, PP ! Intermediate PRESSURE PHYDIA1A.180
*, P_EXNER_FULL_L ! Full Exner pressure for a point/level PHYDIA1A.181
*, BOTTOM_QC_PRESS ! pressure of the bottom QC level PHYDIA1A.182
*, TOP_QC_PRESS ! pressure of the top QC level PHYDIA1A.183
*, PZ(P_FIELD) ! Pressure surface in which results reQD PHYDIA1A.184
*, T(P_FIELD) ! Temperature array PHYDIA1A.185
*, Q_SAT(P_FIELD) ! PHYDIA1A.186
*, DUMMY1 PHYDIA1A.187
*, DUMMY2 PHYDIA1A.188
C PHYDIA1A.189
C*--------------------------------------------------------------------- PHYDIA1A.190
C PHYDIA1A.191
C*L EXTERNAL SUBROUTINES CALLED------------%-------------------------- PHYDIA1A.192
EXTERNAL V_INT_Z,V_INT_ZH,PMSL,V_INT_T,P_TO_UV,TROP,V_INT,QSAT PHYDIA1A.193
EXTERNAL QCTROP,THETAW,ICAO_HT,SNOWPR,FREEZE,DUCT,CONTRAIL PHYDIA1A.194
EXTERNAL DIA_THADV PHYDIA1A.195
EXTERNAL V_INT_TP CW090293.7
C*------------------------------------------------------------------ PHYDIA1A.196
CL MAXIMUM VECTOR LENGTH ASSUMED IS (ROWS+1) * ROWLENGTH PHYDIA1A.197
C---------------------------------------------------------------------- PHYDIA1A.198
C*L DEFINE LOCAL VARIABLES PHYDIA1A.199
INTEGER PHYDIA1A.200
* P_POINTS ! NUMBER OF P POINTS NEEDED PHYDIA1A.201
*, ROWS_P ! NUMBER OF P ROWS NEEDED PHYDIA1A.202
*, U_POINTS ! NUMBER OF U POINTS UPDATED PHYDIA1A.203
*, START_U ! START POSITION FOR U POINTS UPDATED PHYDIA1A.204
*, VERT_POINTS ! NUMBER OF POINTS NON-ZERO DIFFUSION COEFFS PHYDIA1A.205
*, ISL ! LOCAL COUNTER PHYDIA1A.206
*, BL ! BOTTOM LEVEL REQD PHYDIA1A.207
*, TL ! TOP LEVEL REQD PHYDIA1A.208
*, T_REF ! REF level for below sfce Temp extrapola PHYDIA1A.209
*, IPOINT ! Reference point for checking output PHYDIA1A.210
&, P_FLD_VALID ! No of P points excluding halos & unused rows GSM1F405.752
PHYDIA1A.211
REAL PHYDIA1A.212
* PRESS_REQD ! The pressure required for THETAW and THADV PHYDIA1A.213
*, T0 ! Temperature required for FREEZE PHYDIA1A.214
*, CP_OVER_G ! Used in model level heights calculation ASW3F304.31
C----------------------------------------------------------------------- PHYDIA1A.215
C R IS GAS CONSTANT FOR DRY AIR PHYDIA1A.216
C CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE PHYDIA1A.217
C PREF IS REFERENCE SURFACE PRESSURE PHYDIA1A.218
C----------------------------------------------------------------------- PHYDIA1A.219
INTEGER K,I,II,IK,ITR ! LOOP COUNTERS IN ROUTINE ADP0F401.30
PHYDIA1A.221
REAL PHYDIA1A.222
& PU,PL PHYDIA1A.223
PARAMETER (CP_OVER_G = CP / G) ASW3F304.33
*IF DEF,MPP GDR2F404.25
*CALL PARVARS GDR2F404.26
*ENDIF GDR2F404.27
*CALL P_EXNERC PHYDIA1A.224
C*---------------------------------------------------------------------- PHYDIA1A.226
PHYDIA1A.227
T_REF=2 ! This value is used in the vertical interpolation and is PHYDIA1A.228
C set to the 2nd model level. This is not dependant on vert res PHYDIA1A.229
ICODE=0 PHYDIA1A.230
PHYDIA1A.231
! Size of p-fields excluding halos and unused rows GSM1F405.753
P_FLD_VALID=LAST_P_FLD_PT-FIRST_FLD_PT+1 GSM1F405.754
GSM1F405.755
*IF DEF,MPP GDR2F404.28
CALL SWAPBOUNDS(P_EXNER_HALF,ROW_LENGTH,P_ROWS, GDR2F404.29
& Offx,Offy,P_LEVELS+1) GDR2F404.30
*ENDIF GDR2F404.31
C----------------------------------------------------------------------- PHYDIA1A.232
C Calculation of pressure at ETA levels PHYDIA1A.233
C----------------------------------------------------------------------- PHYDIA1A.234
DO K=1,P_LEVELS PHYDIA1A.235
DO I=1,P_FIELD PHYDIA1A.236
P(I,K)=AK(K)+BK(K)*PSTAR(I) PHYDIA1A.237
ENDDO PHYDIA1A.238
ENDDO PHYDIA1A.239
C----------------------------------------------------------------------- PHYDIA1A.240
DO I=1,P_FIELD PHYDIA1A.241
PHI_STAR(I)=OROG(I)*G PHYDIA1A.242
ENDDO PHYDIA1A.243
C set logical if model_half_heights needed for later PS100193.2
IF(QHEIGHTS .OR. PS100193.3
* QMODHT .OR. ASW3F304.34
* QTROP_Z .OR. PS100193.4
* QSNPROB .OR. PS100193.5
* QFREEZE_Z .OR. PS100193.6
* QMINUS20_Z ) QMODEL_HALF_HEIGHT=.TRUE. PS100193.7
C PHYDIA1A.245
IF(QMODEL_HALF_HEIGHT) THEN PHYDIA1A.246
CALL V_INT_ZH(P_EXNER_HALF,THETA,Q,PHI_STAR, PHYDIA1A.247
* MODEL_HALF_HEIGHT,P_FIELD,P_LEVELS,Q_LEVELS) PHYDIA1A.248
ELSEIF(QHEIGHTS) THEN PHYDIA1A.249
WRITE(6,100) PHYDIA1A.250
WRITE(6,100) PHYDIA1A.251
100 FORMAT(' WARNING MODEL_HALF_HEIGHTS NOT SWITCHED ON') PHYDIA1A.252
ICODE=1 PHYDIA1A.253
CMESSAGE='PHY_DIAG:to calculate HTs,model half HTs needed first' PHYDIA1A.254
GOTO 1000 PHYDIA1A.255
ENDIF ASW3F304.35
C----------------------------------------------------------------------- ASW3F304.36
CL--------------------------Section 0----------------------------------- ASW3F304.37
CL------------------GEOPOTENTIAL MODEL LEVEL HEIGHTS ITEM 281----------- ASW3F304.38
CL ASW3F304.39
CL section 0 calculation of geopotential height of a model level ASW3F304.40
CL surface. Heights in M ASW3F304.41
C----------------------------------------------------------------------- ASW3F304.42
IF (QMODHT) THEN ASW3F304.43
DO K=1,Q_LEVELS ASW3F304.44
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.756
PU = AKH(K+1)+BKH(K+1)*PSTAR(I) ASW3F304.46
PL = AKH(K)+BKH(K)*PSTAR(I) ASW3F304.47
P_EXNER_FULL_L = P_EXNER_C( P_EXNER_HALF(I,K+1), ASW3F304.48
& P_EXNER_HALF(I,K),PU,PL,KAPPA ) ASW3F304.49
MODHT(I,K) = MODEL_HALF_HEIGHT(I,K) + CP_OVER_G* ARB1F400.29
& (1.0+C_VIRTUAL*Q(I,K))*THETA(I,K) ASW3F304.51
& *(P_EXNER_HALF(I,K) - P_EXNER_FULL_L) ASW3F304.52
ENDDO ASW3F304.54
ENDDO ASW3F304.55
IF(P_LEVELS .GT. Q_LEVELS) THEN ASW3F304.56
DO K=Q_LEVELS+1,P_LEVELS ASW3F304.57
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.757
PU = AKH(K+1)+BKH(K+1)*PSTAR(I) ASW3F304.59
PL = AKH(K)+BKH(K)*PSTAR(I) ASW3F304.60
P_EXNER_FULL_L = P_EXNER_C( P_EXNER_HALF(I,K+1), ASW3F304.61
& P_EXNER_HALF(I,K),PU,PL,KAPPA ) ASW3F304.62
MODHT(I,K) = MODEL_HALF_HEIGHT(I,K) + CP_OVER_G* ARB1F400.30
& THETA(I,K) ASW3F304.64
& *(P_EXNER_HALF(I,K) - P_EXNER_FULL_L) ASW3F304.65
ENDDO ASW3F304.67
ENDDO ASW3F304.68
END IF ASW3F304.69
ENDIF PHYDIA1A.256
C----------------------------------------------------------------------- PHYDIA1A.257
CL--------------------------Section 1----------------------------------- PHYDIA1A.258
CL------------------GEOPOTENTIAL HEIGHTS ITEM 202----------------------- PHYDIA1A.259
CL PHYDIA1A.260
CL section 2 calculation of geopotential height of an arbitary PHYDIA1A.261
CL pressure surface. Pressure input is in Pascals and heights in M PHYDIA1A.262
C----------------------------------------------------------------------- PHYDIA1A.263
C First check the height codes are valid and point to STASH_PRESSURE PHYDIA1A.264
C----------------------------------------------------------------------- PHYDIA1A.265
IF(QHEIGHTS) THEN PHYDIA1A.266
DO K=1,HTS_LEVS PHYDIA1A.267
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.758
PZ(I)=HTS_PRESS(K)*100.0 ! convert to pascals PHYDIA1A.269
ENDDO PHYDIA1A.270
CALL V_INT_Z(PZ,P(1,Z_REF),PSTAR,P_EXNER_HALF,THETA,Q, PHYDIA1A.271
& MODEL_HALF_HEIGHT,HEIGHTS(1,K),P_FIELD,P_LEVELS,Q_LEVELS, PHYDIA1A.272
& Z_REF,AKH,BKH,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.759
ENDDO PHYDIA1A.274
ENDIF PHYDIA1A.275
C----------------------------------------------------------------------- PHYDIA1A.276
CL--------------------------Section 2----------------------------------- PHYDIA1A.277
CL------------------MEAN SEA LEVEL PRESSURE ITEM 222-------------------- PHYDIA1A.278
C----------------------------------------------------------------------- PHYDIA1A.279
IF(QPMSL) THEN PHYDIA1A.280
CALL PMSL(P_MSL,P(1,Z_REF),PSTAR,P_EXNER_HALF,THETA,Q,PHI_STAR GSM1F405.760
& ,P_FIELD,P_LEVELS,Q_LEVELS,Z_REF,AKH,BKH,FIRST_FLD_PT GSM1F405.761
& ,LAST_P_FLD_PT) GSM1F405.762
ENDIF PHYDIA1A.283
C----------------------------------------------------------------------- PHYDIA1A.284
CL--------------------------Section 3----------------------------------- PHYDIA1A.285
CL------------------TEMPERATURES ON PRESSURE ITEM 203------------------- PHYDIA1A.286
CL PHYDIA1A.287
CL Pressure input is in Pascals and potential temperature is input PHYDIA1A.288
CL in Kelvin PHYDIA1A.289
C----------------------------------------------------------------------- PHYDIA1A.290
IF(QT_P) THEN PHYDIA1A.291
DO K=1,T_P_LEVS PHYDIA1A.292
DO I=1,P_FIELD PHYDIA1A.293
PZ(I)=T_P_PRESS(K)*100.0 ! convert to pascals PHYDIA1A.294
ENDDO PHYDIA1A.295
! From Vn 4.5, L_VINT_TP is set through UMUI GDR4F405.23
IF (L_VINT_TP) THEN GDR4F405.24
CALL V_INT_TP(T_P(1,K),PZ,P(1,T_REF),PSTAR,P_EXNER_HALF, CW090293.12
& THETA,P_FIELD,P_LEVELS,T_REF,AKH,BKH) CW090293.13
ELSE CW090293.14
CALL V_INT_T(T_P(1,K),PZ,P(1,T_REF),PSTAR,P_EXNER_HALF, CW090293.15
& THETA,P_FIELD,P_LEVELS,T_REF,AKH,BKH,FIRST_FLD_PT GSM1F405.763
& ,LAST_P_FLD_PT) GSM1F405.764
ENDIF CW090293.17
ENDDO PHYDIA1A.298
ENDIF PHYDIA1A.299
C----------------------------------------------------------------------- PHYDIA1A.300
CL--------------------------Section 4----------------------------------- PHYDIA1A.301
CL------------------RELATIVE HUM ON PRESSURE ITEM 204------------------- PHYDIA1A.302
C----------------------------------------------------------------------- PHYDIA1A.303
IF (QREL_HUMID_P) THEN PHYDIA1A.304
C----------------------------------------------------------------------- PHYDIA1A.305
C Calculation of temperature at theta points PHYDIA1A.306
C Calculation of relative humidity at pressure points PHYDIA1A.307
C----------------------------------------------------------------------- PHYDIA1A.308
DO K=1,Q_LEVELS PHYDIA1A.309
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.765
PU=PSTAR(I)*BKH(K+1) + AKH(K+1) PHYDIA1A.311
PL=PSTAR(I)*BKH(K) + AKH(K) PHYDIA1A.312
P_EXNER_FULL_L= PHYDIA1A.313
& P_EXNER_C(P_EXNER_HALF(I,K+1),P_EXNER_HALF(I,K), PHYDIA1A.314
& PU,PL,KAPPA) PHYDIA1A.315
T(I)=THETA(I,K)*P_EXNER_FULL_L PHYDIA1A.316
ENDDO PHYDIA1A.317
CALL QSAT(Q_SAT(FIRST_FLD_PT),T(FIRST_FLD_PT) GSM1F405.766
& ,P(FIRST_FLD_PT,K),P_FLD_VALID) GSM1F405.767
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.768
WORK1(I,K)=Q(I,K)/Q_SAT(I)*100.0 ! Rel humidity in WORK1 PHYDIA1A.320
! Relative humidity should NOT be limited to 100% but do so in old ADM0F405.299
! versions of the cloud scheme to maintain bit comparison. ADM0F405.300
IF(WORK1(I,K).GT.100.0 .AND. .NOT. L_LSPICE) THEN ADM0F405.301
WORK1(I,K)=100.0 PHYDIA1A.323
ELSEIF(WORK1(I,K).LT.0.) THEN PHYDIA1A.324
WORK1(I,K)=0.0 PHYDIA1A.325
ENDIF PHYDIA1A.326
ENDDO PHYDIA1A.327
ENDDO PHYDIA1A.328
C----------------------------------------------------------------------- PHYDIA1A.329
CL Interpolate relative humidity on a pressure surface PHYDIA1A.330
C----------------------------------------------------------------------- PHYDIA1A.331
DO K=1,REL_HUMID_P_LEVS PHYDIA1A.332
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.769
PZ(I)=REL_HUMID_PRESS(K)*100.0 ! convert to Pascals PHYDIA1A.334
ENDDO PHYDIA1A.335
CALL V_INT(P,PZ,WORK1,REL_HUMID_P(1,K), PHYDIA1A.336
& P_FIELD,Q_LEVELS,DUMMY1,DUMMY2,.FALSE. GSM1F405.770
& ,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.771
ENDDO PHYDIA1A.338
ENDIF PHYDIA1A.339
C----------------------------------------------------------------------- PHYDIA1A.340
C------------------TROPOPAUSE ITEMS 214 215 216 217-------------------- PHYDIA1A.341
CL----------Section 5 calculation of tropopause temp/ht and pressure---- PHYDIA1A.342
C----------------------------------------------------------------------- PHYDIA1A.343
IF (QTROP_T.AND.QTROP_P.AND.QTROP_Z) THEN PHYDIA1A.344
IF(QMODEL_HALF_HEIGHT) THEN PHYDIA1A.345
C----------------------------------------------------------------------- PHYDIA1A.346
C Call tropopause program with quality control PHYDIA1A.347
C----------------------------------------------------------------------- PHYDIA1A.348
BOTTOM_QC_PRESS=60000.0 ! pressure of the bottom QC level PHYDIA1A.349
TOP_QC_PRESS=10000.0 ! pressure of the top QC level PHYDIA1A.350
BOTTOM_QC_LEVEL=1 PHYDIA1A.351
TOP_QC_LEVEL=P_LEVELS PHYDIA1A.352
DO K=1,P_LEVELS PHYDIA1A.353
PP= AK(K) + BK(K) * 100000 ! Pstar assumed to be 1000MB PHYDIA1A.354
IF(BOTTOM_QC_PRESS.LT.PP) THEN PHYDIA1A.355
BOTTOM_QC_LEVEL=K PHYDIA1A.356
ENDIF PHYDIA1A.357
IF(TOP_QC_PRESS.LT.PP) THEN PHYDIA1A.358
TOP_QC_LEVEL=K PHYDIA1A.359
ENDIF PHYDIA1A.360
ENDDO PHYDIA1A.361
DO K=1,P_LEVELS+1 PHYDIA1A.362
DO I=1,P_FIELD PHYDIA1A.363
WORK1(I,K)=AKH(K)+BKH(K)*PSTAR(I) ! Calculate P_HALF PHYDIA1A.364
ENDDO PHYDIA1A.365
ENDDO PHYDIA1A.366
CALL QCTROP(THETA,WORK1,P_EXNER_HALF PHYDIA1A.367
& ,MODEL_HALF_HEIGHT,TROP_T,TROP_P,TROP_Z,P_FIELD,P_LEVELS PHYDIA1A.368
& ,P(1,Z_REF),P(1,T_REF),PSTAR,Q,Q_LEVELS,Z_REF,T_REF,Z_REF PHYDIA1A.369
& ,BOTTOM_QC_LEVEL,TOP_QC_LEVEL,AKH,BKH GSM1F405.772
& ,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.773
C Last Z_REF used as Min TROP L PHYDIA1A.371
PHYDIA1A.372
C----------------------------------------------------------------------- PHYDIA1A.373
C This piece of code is to ensure sensible values of TROP ht and TEMP. PHYDIA1A.374
C Although,the TROP height is also produced it can be rubbish so for PHYDIA1A.375
C the moment it is best to call V_INTZ PHYDIA1A.376
C----------------------------------------------------------------------- PHYDIA1A.377
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.774
IF(TROP_T(I).LT.0) THEN PHYDIA1A.379
TROP_T(I)=199.0 PHYDIA1A.380
ENDIF PHYDIA1A.381
IF(TROP_P(I).LT.0) THEN PHYDIA1A.382
TROP_P(I)=111.0 PHYDIA1A.383
ENDIF PHYDIA1A.384
ENDDO PHYDIA1A.385
C CALL V_INT_Z(TROP_P,P(1,Z_REF),PSTAR,P_EXNER_HALF,THETA,Q, PHYDIA1A.386
C & MODEL_HALF_HEIGHT,TROP_Z,P_FIELD,P_LEVELS,Q_LEVELS,Z_REF, PHYDIA1A.387
C & AKH,BKH,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.775
C----------------------------------------------------------------------- PHYDIA1A.389
C Calculates ICAO heights from Trop pressure field PHYDIA1A.390
C----------------------------------------------------------------------- PHYDIA1A.391
IF(QTROP_ICAO) THEN PHYDIA1A.392
CALL ICAO_HT(TROP_P(FIRST_FLD_PT),P_FLD_VALID GSM1F405.776
& ,TROP_ICAO(FIRST_FLD_PT)) GSM1F405.777
ENDIF PHYDIA1A.394
C PHYDIA1A.395
ELSE PHYDIA1A.396
WRITE(6,444) PHYDIA1A.397
WRITE(6,444) PHYDIA1A.398
444 FORMAT(' Subroutine TROP not called No MODEL_HALF_HEIGHTS') PHYDIA1A.399
ENDIF PHYDIA1A.400
ELSEIF(QTROP_T.NEQV.QTROP_P.OR.QTROP_T.NEQV.QTROP_Z.OR.QTROP_P. PHYDIA1A.401
& NEQV.QTROP_Z) THEN PHYDIA1A.402
WRITE(6,*)'Subroutine TROP not called-tropopause temperature,' GIE0F403.464
WRITE(6,*)' pressure and height all need to be selected' GIE0F403.465
ENDIF PHYDIA1A.405
C----------------------------------------------------------------------- PHYDIA1A.406
CL-------------------SECTION 6 ITEM 205--------------------------------- PHYDIA1A.407
CL Calculation of Wet Bulb Potential Temperature PHYDIA1A.408
C----------------------------------------------------------------------- PHYDIA1A.409
IF (QWBPT) THEN PHYDIA1A.410
DO K=1,WBPT_LEVS PHYDIA1A.411
PRESS_REQD=WBPT_PRESS(K)*100.0 ! convert to Pascals PHYDIA1A.412
ICODE=0 PHYDIA1A.413
CALL THETAW PHYDIA1A.414
1 (PRESS_REQD,THETA,Q,P,PSTAR,P_EXNER_HALF,WBPT(1,K), PHYDIA1A.415
2 AK,BK,AKH,BKH,P_FIELD,P_LEVELS,Q_LEVELS,T_REF, PHYDIA1A.416
& FIRST_FLD_PT,LAST_P_FLD_PT, GSM1F405.778
3 ICODE,CMESSAGE) PHYDIA1A.417
IF(ICODE.NE.0) THEN PHYDIA1A.418
ICODE=1 PHYDIA1A.419
CMESSAGE='PHY_DIAG error in calculation of WBPT' PHYDIA1A.420
GOTO 1000 PHYDIA1A.421
ENDIF PHYDIA1A.422
ENDDO PHYDIA1A.423
ENDIF PHYDIA1A.424
C----------------------------------------------------------------------- PHYDIA1A.425
CL---------------------SECTION 7 ITEM 206------------------------------- PHYDIA1A.426
CL Calculation of Snow probability PHYDIA1A.427
C----------------------------------------------------------------------- PHYDIA1A.428
IF (QSNPROB)THEN PHYDIA1A.429
CALL SNOWPR PHYDIA1A.430
1 (P,PSTAR,P_EXNER_HALF,THETA,Q,MODEL_HALF_HEIGHT, PHYDIA1A.431
2 P_FIELD,P_LEVELS,Q_LEVELS,Z_REF,AKH,BKH, PHYDIA1A.432
& SNPROB,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.779
ENDIF PHYDIA1A.434
C----------------------------------------------------------------------- PHYDIA1A.435
CL--------------------SECTION 8 ITEMS 207/8/18-------------------------- PHYDIA1A.436
CL Calculation of heights and pressure of -20 isotherm PHYDIA1A.437
C----------------------------------------------------------------------- PHYDIA1A.438
IF (QMINUS20_P.AND.QMINUS20_Z) THEN PHYDIA1A.439
IF(QMODEL_HALF_HEIGHT)THEN PHYDIA1A.440
C----------------------------------------------------------------------- PHYDIA1A.441
C Set T0 to 253.16K (-20 deg C) PHYDIA1A.442
C----------------------------------------------------------------------- PHYDIA1A.443
T0=253.16 PHYDIA1A.444
C----------------------------------------------------------------------- PHYDIA1A.445
C Calculate temperatures at full levels (stored in WORK1) PHYDIA1A.446
C----------------------------------------------------------------------- PHYDIA1A.447
DO K=1,P_LEVELS PHYDIA1A.448
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.780
PU=PSTAR(I)*BKH(K+1) + AKH(K+1) PHYDIA1A.450
PL=PSTAR(I)*BKH(K) + AKH(K) PHYDIA1A.451
P_EXNER_FULL_L= PHYDIA1A.452
& P_EXNER_C(P_EXNER_HALF(I,K+1),P_EXNER_HALF(I,K), PHYDIA1A.453
& PU,PL,KAPPA) PHYDIA1A.454
WORK1(I,K)=THETA(I,K)*P_EXNER_FULL_L PHYDIA1A.455
ENDDO PHYDIA1A.456
ENDDO PHYDIA1A.457
C----------------------------------------------------------------------- PHYDIA1A.458
CALL FREEZE PHYDIA1A.459
1 (T0,P,THETA,WORK1,P_EXNER_HALF,PSTAR,Q,MODEL_HALF_HEIGHT, PHYDIA1A.460
2 OROG,MINUS20_Z,MINUS20_P, PHYDIA1A.461
& P_FIELD,P_LEVELS,Q_LEVELS,Z_REF,AKH,BKH GSM1F405.781
& ,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.782
C----------------------------------------------------------------------- PHYDIA1A.463
IF(QMINUS20_ICAO)THEN PHYDIA1A.464
CALL ICAO_HT(MINUS20_P(FIRST_FLD_PT),P_FLD_VALID GSM1F405.783
& ,MINUS20_ICAO(FIRST_FLD_PT)) GSM1F405.784
ENDIF PHYDIA1A.466
C----------------------------------------------------------------------- PHYDIA1A.467
ELSE PHYDIA1A.468
WRITE(6,556) PHYDIA1A.469
WRITE(6,556) PHYDIA1A.470
556 FORMAT(' Subroutine FREEZE not called-no MODEL_HALF_HEIGHTs') PHYDIA1A.471
ENDIF PHYDIA1A.472
ELSEIF(QMINUS20_P.NEQV.QMINUS20_Z) THEN PHYDIA1A.473
WRITE(6,*)'Subroutine FREEZE not called-both -20 level pressure' GIE0F403.466
WRITE(6,*)'and height need to be selected' GIE0F403.467
ENDIF PHYDIA1A.476
C----------------------------------------------------------------------- PHYDIA1A.477
CL-------------------SECTION 9 ITEMS 209-211---------------------------- PHYDIA1A.478
CL Calculation of Freezing level heights and pressure PHYDIA1A.479
C----------------------------------------------------------------------- PHYDIA1A.480
IF (QFREEZE_Z.AND.QFREEZE_P)THEN PHYDIA1A.481
IF(QMODEL_HALF_HEIGHT)THEN PHYDIA1A.482
C----------------------------------------------------------------------- PHYDIA1A.483
C Set T0 to 273.16K (0 deg C) PHYDIA1A.484
C----------------------------------------------------------------------- PHYDIA1A.485
T0=273.16 PHYDIA1A.486
C----------------------------------------------------------------------- PHYDIA1A.487
C Calculate temperatures at full levels (stored in WORK1) PHYDIA1A.488
C----------------------------------------------------------------------- PHYDIA1A.489
DO K=1,P_LEVELS PHYDIA1A.490
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.785
PU=PSTAR(I)*BKH(K+1) + AKH(K+1) PHYDIA1A.492
PL=PSTAR(I)*BKH(K) + AKH(K) PHYDIA1A.493
P_EXNER_FULL_L= PHYDIA1A.494
& P_EXNER_C(P_EXNER_HALF(I,K+1),P_EXNER_HALF(I,K), PHYDIA1A.495
& PU,PL,KAPPA) PHYDIA1A.496
WORK1(I,K)=THETA(I,K)*P_EXNER_FULL_L PHYDIA1A.497
ENDDO PHYDIA1A.498
ENDDO PHYDIA1A.499
C----------------------------------------------------------------------- PHYDIA1A.500
CALL FREEZE PHYDIA1A.501
1 (T0,P,THETA,WORK1,P_EXNER_HALF,PSTAR,Q,MODEL_HALF_HEIGHT, PHYDIA1A.502
2 OROG,FREEZE_Z,FREEZE_P, PHYDIA1A.503
& P_FIELD,P_LEVELS,Q_LEVELS,Z_REF,AKH,BKH GSM1F405.786
& ,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.787
IF(QFREEZE_ICAO)THEN PHYDIA1A.505
CALL ICAO_HT(FREEZE_P(FIRST_FLD_PT),P_FLD_VALID GSM1F405.788
& ,FREEZE_ICAO(FIRST_FLD_PT)) GSM1F405.789
ENDIF PHYDIA1A.507
C----------------------------------------------------------------------- PHYDIA1A.508
ELSE PHYDIA1A.509
WRITE(6,555) PHYDIA1A.510
WRITE(6,555) PHYDIA1A.511
555 FORMAT(' Subroutine FREEZE not called-no MODEL_HALF_HEIGHTs') PHYDIA1A.512
ENDIF PHYDIA1A.513
ELSEIF(QFREEZE_P.NEQV.QFREEZE_Z) THEN PHYDIA1A.514
WRITE(6,*)'Subroutine FREEZE not called - both freezing level' GIE0F403.468
WRITE(6,*)'pressure and height need to be selected' GIE0F403.469
ENDIF PHYDIA1A.517
C----------------------------------------------------------------------- PHYDIA1A.518
CL-----------------SECTION 10 ITEMS 220/1------------------------------- PHYDIA1A.519
CL Calculation of duct height and intensity PHYDIA1A.520
C----------------------------------------------------------------------- PHYDIA1A.521
IF(QDUCT_INT.AND.QDUCT_HEIGHT)THEN PHYDIA1A.522
CALL DUCT PHYDIA1A.523
1 (PSTAR,TSTAR,THETA,Q,U,V, PHYDIA1A.524
2 AK,BK,AKH,BKH,LAND,DUCT_HEIGHT,DUCT_INT,P_EXNER_HALF, PHYDIA1A.525
3 P_LEVELS,Q_LEVELS,ROW_LENGTH,P_ROWS,U_ROWS,P_FIELD,U_FIELD) PHYDIA1A.526
ELSEIF(QDUCT_INT.NEQV.QDUCT_HEIGHT)THEN PHYDIA1A.527
WRITE(6,*)'Subroutine DUCT not called - both duct height and' GIE0F403.470
WRITE(6,*)'intensity need to be selected' GIE0F403.471
ENDIF PHYDIA1A.530
C----------------------------------------------------------------------- PHYDIA1A.531
CL-----------------SECTION 11 ITEMS 212/3------------------------------- PHYDIA1A.532
CL Calculation of contrail lower and upper height PHYDIA1A.533
C----------------------------------------------------------------------- PHYDIA1A.534
IF(QCONTRA_LOWERHT.AND.QCONTRA_UPPERHT)THEN PHYDIA1A.535
C----------------------------------------------------------------------- PHYDIA1A.536
C Calculate temperatures at full levels (stored in WORK1) PHYDIA1A.537
C----------------------------------------------------------------------- PHYDIA1A.538
DO K=1,P_LEVELS PHYDIA1A.539
DO I=1,P_FIELD PHYDIA1A.540
PU=PSTAR(I)*BKH(K+1) + AKH(K+1) PHYDIA1A.541
PL=PSTAR(I)*BKH(K) + AKH(K) PHYDIA1A.542
P_EXNER_FULL_L= PHYDIA1A.543
& P_EXNER_C(P_EXNER_HALF(I,K+1),P_EXNER_HALF(I,K), PHYDIA1A.544
& PU,PL,KAPPA) PHYDIA1A.545
WORK1(I,K)=THETA(I,K)*P_EXNER_FULL_L PHYDIA1A.546
ENDDO PHYDIA1A.547
ENDDO PHYDIA1A.548
C----------------------------------------------------------------------- PHYDIA1A.549
CALL CONTRAIL PHYDIA1A.550
1 (P,WORK1,PSTAR,P_EXNER_HALF,THETA,CONTRAIL_UPPERHT, PHYDIA1A.551
2 CONTRAIL_LOWERHT,P_FIELD,P_LEVELS,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.790
ELSEIF(QCONTRA_UPPERHT.NEQV.QCONTRA_LOWERHT)THEN PHYDIA1A.553
WRITE(6,*)'Subroutine CONTRAIL not called - both lower height and' GIE0F403.472
WRITE(6,*)'upper height need to be selected' GIE0F403.473
ENDIF PHYDIA1A.556
C----------------------------------------------------------------------- PHYDIA1A.557
C----------------------------------------------------------------------- PHYDIA1A.558
CL-----------------SECTION 12 ITEMS 219/223----------------------------- PHYDIA1A.559
CL Calculation of thermal advection PHYDIA1A.560
C----------------------------------------------------------------------- PHYDIA1A.561
IF(QTH_ADV_SINGLE.OR.QTH_ADV_MEAN)THEN PHYDIA1A.562
C----------------------------------------------------------------------- PHYDIA1A.563
CL Interpolate P field from P/T points to U/V points for all levels. RB070693.2
C----------------------------------------------------------------------- PHYDIA1A.566
DO K=1,P_LEVELS PHYDIA1A.567
CALL P_TO_UV PHYDIA1A.568
1 (P(1,K),WORK1(1,K),P_FIELD,U_FIELD,ROW_LENGTH, GRR0F304.17
2 P_ROWS) RB070693.5
ENDDO PHYDIA1A.570
ENDIF PHYDIA1A.571
C Change effective dimension of WORK1 for input to DIA_THADV GRR0F304.18
C (only first U_FIELD points defined for each level from P_TO_UV) GRR0F304.19
CALL CHANGE_DIMENS(WORK1,P_FIELD,U_FIELD,P_LEVELS,ICODE) GRR0F304.20
C----------------------------------------------------------------------- PHYDIA1A.572
CL Item 219 single level thermal advection--------------------------- PHYDIA1A.573
C----------------------------------------------------------------------- PHYDIA1A.574
IF(QTH_ADV_SINGLE)THEN PHYDIA1A.575
C----------------------------------------------------------------------- PHYDIA1A.576
CL Loop over required levels PHYDIA1A.577
C----------------------------------------------------------------------- PHYDIA1A.578
DO K=1,TH_ADV_P_LEVS PHYDIA1A.579
PRESS_REQD=TH_ADV_PRESS(K)*100.0 ! Convert to pascals PHYDIA1A.580
CALL DIA_THADV PHYDIA1A.581
1 (U,V,THETA,P,WORK1,PSTAR,PRESS_REQD,P_EXNER_HALF, PHYDIA1A.582
2 TH_ADV_SINGLE(1,K),P_FIELD,U_FIELD,P_LEVELS,ROW_LENGTH, PHYDIA1A.583
3 P_ROWS,SEC_U_LATITUDE,EW_SPACE,NS_SPACE,AKH,BKH) PHYDIA1A.584
PHYDIA1A.585
ENDDO PHYDIA1A.586
ENDIF PHYDIA1A.587
C----------------------------------------------------------------------- PHYDIA1A.588
CL Item 223 mean thermal advection over levels 850/700/500mb- PHYDIA1A.589
C----------------------------------------------------------------------- PHYDIA1A.590
IF(QTH_ADV_MEAN)THEN PHYDIA1A.591
C----------------------------------------------------------------------- PHYDIA1A.592
C Initialise array PHYDIA1A.593
C----------------------------------------------------------------------- PHYDIA1A.594
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.791
TH_ADV_MEAN(I)=0.0 PHYDIA1A.596
ENDDO PHYDIA1A.597
C----------------------------------------------------------------------- PHYDIA1A.598
CL Loop over levels 850/700/500mb PHYDIA1A.599
C----------------------------------------------------------------------- PHYDIA1A.600
DO K=1,3 PHYDIA1A.601
IF(K.EQ.1)PRESS_REQD=85000.0 PHYDIA1A.602
IF(K.EQ.2)PRESS_REQD=70000.0 PHYDIA1A.603
IF(K.EQ.3)PRESS_REQD=50000.0 PHYDIA1A.604
CALL DIA_THADV PHYDIA1A.605
1 (U,V,THETA,P,WORK1,PSTAR,PRESS_REQD,P_EXNER_HALF,WORK2, PHYDIA1A.606
2 P_FIELD,U_FIELD,P_LEVELS,ROW_LENGTH,P_ROWS, PHYDIA1A.607
3 SEC_U_LATITUDE,EW_SPACE,NS_SPACE,AKH,BKH) PHYDIA1A.608
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.792
TH_ADV_MEAN(I)=TH_ADV_MEAN(I)+WORK2(I) PHYDIA1A.610
ENDDO PHYDIA1A.611
ENDDO PHYDIA1A.612
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.793
TH_ADV_MEAN(I)=TH_ADV_MEAN(I)/3.0 PHYDIA1A.614
ENDDO PHYDIA1A.615
ENDIF PHYDIA1A.616
C----------------------------------------------------------------------- PHYDIA1A.617
CL Section 13 ITEMS 224 products of other fields PHYDIA1A.618
C----------------------------------------------------------------------- PHYDIA1A.619
CL Item 224 - height**2 on pressure levels PHYDIA1A.620
CL Must be requested on identical levels to height on pressure PHYDIA1A.621
CL levels. PHYDIA1A.622
C PHYDIA1A.623
IF (QHTS2) THEN PHYDIA1A.624
DO K=1,H2_P_LEVS PHYDIA1A.625
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.794
HEIGHTS2(I,K) = HEIGHTS(I,H2_IND(K))*HEIGHTS(I,H2_IND(K)) PHYDIA1A.627
ENDDO PHYDIA1A.628
ENDDO PHYDIA1A.629
ENDIF ADP0F401.31
C----------------------------------------------------------------------- ADP0F401.32
CL-----------------SECTION 14 ITEMS 226-254----------------------------- ADP0F401.33
CL Interpolate tracers onto pressure levels ADP0F401.34
C----------------------------------------------------------------------- ADP0F401.35
IF (TR_VARS.GT.0) THEN ADP0F401.36
DO ITR=1,TR_VARS ADP0F401.37
IF (SF_TRACER(ITR)) THEN ADP0F401.38
DO K=1,TR_PRESS_LEVS(ITR) ADP0F401.39
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.795
PZ(I)=TR_PRESS(ITR,K)*100.0 ! convert to Pascals ADP0F401.41
C----------------------------------------------------------------------- ADP0F401.42
C Extract current tracer field from STASHWORK/TRACERWORK ADP0F401.43
C----------------------------------------------------------------------- ADP0F401.44
TRACER_P(I)=TRACERWORK( PT_TRACER(ITR) - 1 + I + ADP0F401.45
& ( P_FIELD * (K-1) ) ) ADP0F401.46
ENDDO ADP0F401.47
CALL V_INT(P,PZ,TRACERS(1,1,ITR),TRACER_P, ADP0F401.48
& P_FIELD,TR_LEVELS,DUMMY1,DUMMY2,.FALSE. GSM1F405.796
& ,FIRST_FLD_PT,LAST_P_FLD_PT) GSM1F405.797
C----------------------------------------------------------------------- ADP0F401.50
C Copy interpolated tracer field back to STASHWORK/TRACERWORK ADP0F401.51
C----------------------------------------------------------------------- ADP0F401.52
DO I=FIRST_FLD_PT,LAST_P_FLD_PT GSM1F405.798
TRACERWORK( PT_TRACER(ITR) - 1 + I + ADP0F401.54
& ( P_FIELD * (K-1) ) ) = TRACER_P(I) ADP0F401.55
ENDDO ADP0F401.56
ENDDO ADP0F401.57
ENDIF ADP0F401.58
ENDDO ADP0F401.59
ENDIF PHYDIA1A.630
C----------------------------------------------------------------------- PHYDIA1A.631
1000 CONTINUE PHYDIA1A.632
RETURN PHYDIA1A.633
END PHYDIA1A.634
*ENDIF PHYDIA1A.635