*IF DEF,C90_1A,OR,DEF,C90_2A,OR,DEF,C90_2B AAD2F404.301
C ******************************COPYRIGHT****************************** GTS2F400.11737
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.11738
C GTS2F400.11739
C Use, duplication or disclosure of this code is subject to the GTS2F400.11740
C restrictions as set forth in the contract. GTS2F400.11741
C GTS2F400.11742
C Meteorological Office GTS2F400.11743
C London Road GTS2F400.11744
C BRACKNELL GTS2F400.11745
C Berkshire UK GTS2F400.11746
C RG12 2SZ GTS2F400.11747
C GTS2F400.11748
C If no contract has been raised with this copy of the code, the use, GTS2F400.11749
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.11750
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.11751
C Modelling at the above address. GTS2F400.11752
C ******************************COPYRIGHT****************************** GTS2F400.11753
C GTS2F400.11754
CLL SUBROUTINE VISBTY ----------------------------------------------- VISBTY1A.3
CLL VISBTY1A.4
CLL PURPOSE: VISBTY1A.5
CLL Process fields of temperature, specific humidity, cloud liquid APC0F405.117
CLL water or ThetaL and qt to give visibility in metres. APC0F405.118
CLL Calculated at model level (eg bottom eta level 25m) VISBTY1A.8
CLL or level within surface layer eg screen ht ( 1.5M ) VISBTY1A.9
CLL APC3F304.136
CLL APC3F304.138
CLL Model Modification history from model version 3.0: VISBTY1A.16
CLL version Date VISBTY1A.17
CLL 3.1 23/10/92 New deck author P.Smith VISBTY1A.18
CLL 3.1 20/01/93 New deck - used as mods at 2.7 & 2.8/3.0. VISBTY1A.19
CLL Interfacing done by R.T.H.Barnes. VISBTY1A.20
CLL 3.2 29/04/93 CCN and derived constants moved to MODECK C_VISBTY PC120793.175
CLL Programmer Pete Clark. PC120793.176
CLL 3.4 07/06/94 Aerosol field introduced. Programmer Pete Clark. APC3F304.146
CLL 4.0 05/09/95 Lower limit to aerosol introduced. Pete Clark. APC0F400.6
! 4.4 09/01/97 Only liquid water and not cloud ice is now used AYY1F404.112
! to calculate visibility. Damian Wilson. AYY1F404.113
!LL 4.5 29/04/98 Scheme replaced with NIMROD visibility diagnostic APC0F405.119
CLL VISBTY1A.21
CLL Programming standard: U M Doc. Paper No. 4 VISBTY1A.22
CLL VISBTY1A.23
CLL Logical components covered : VISBTY1A.24
CLL VISBTY1A.25
CLL Project task: VISBTY1A.26
CLL VISBTY1A.27
CLL External documentation VISBTY1A.28
CLL Forecasting Research Scientific Paper NO.4 VISBTY1A.29
CLL Diagnosis of visibility in the UK Met Office Mesoscale Model VISBTY1A.30
CLL and the use of a visibility analysis to constrain initial VISBTY1A.31
CLL conditions. SP Ballard, BJ Wright, BW Golding 1992 VISBTY1A.32
! NIMROD diagnostic: APC0F405.120
! Wright, B. J., 1997: Improvements to the Nimrod Visibility APC0F405.121
! Analysis/Forecast System. FR-Div. Tech. Rep., No. 217. APC0F405.122
! Wright, B. J., 1997: A New Visibility Analysis/Forecast System APC0F405.123
! for Nimrod. Met. Office FR Tech Rep., No. 222. APC0F405.124
CLL VISBTY1A.33
CLLEND---------------------------------------------------------------- VISBTY1A.34
C VISBTY1A.35
C*L Arguments:------------------------------------------------------- VISBTY1A.36
SUBROUTINE VISBTY 2,1VISBTY1A.37
+ (AK,BK,PSTAR,T,Q,QCL,QCF !INPUT APC0F405.125
& ,AEROSOL, PROB, RHCRIT, L_MURK !INPUT APC0F405.126
+ ,P_FIELD !INPUT VISBTY1A.39
+ ,VISIBILITY) !OUTPUT VISBTY1A.40
IMPLICIT NONE VISBTY1A.41
C--------------------------------------------------------------------- VISBTY1A.42
C Workspace usage:---------------------------------------------------- VISBTY1A.43
C 3 real arrays of size P_FIELD APC0F405.127
C*-------------------------------------------------------------------- VISBTY1A.45
C*L------------------------------------------------------------------- PC120793.177
C input variables----------------------------------------------------- VISBTY1A.47
C--------------------------------------------------------------------- VISBTY1A.48
INTEGER VISBTY1A.49
* P_FIELD ! IN NO. points in field. VISBTY1A.50
REAL VISBTY1A.51
& AK,BK ! IN Ak and Bk of level APC0F405.128
& ,PSTAR(P_FIELD) ! IN Surface pressure APC0F405.129
& ,T(P_FIELD) ! IN Temperature APC0F405.130
& ,Q(P_FIELD) ! IN Qt APC0F405.131
& ,QCL(P_FIELD) ! IN cloud water array. APC0F405.132
& ,QCF(P_FIELD) ! IN cloud ice array. APC0F405.133
& ,AEROSOL(P_FIELD) ! IN Aerosol mixing ratio(ug/kg) APC3F304.148
& ,PROB ! IN Probability level ( e.g 0.5 APC0F405.134
! corresponds to median). APC0F405.135
& ,RHCRIT ! IN Citical RH (determines APC0F405.136
! width of distribiution) APC0F405.137
LOGICAL APC3F304.149
& L_MURK ! IN : Aerosol present APC3F304.150
C--------------------------------------------------------------------- PC120793.178
C output variables---------------------------------------------------- VISBTY1A.56
C--------------------------------------------------------------------- VISBTY1A.57
REAL VISBTY1A.58
& VISIBILITY(P_FIELD) ! OUT visibility array. APC3F304.151
C*-------------------------------------------------------------------- VISBTY1A.60
C*L------------------------------------------------------------------- PC120793.179
C Local varables:----------------------------------------------------- VISBTY1A.61
C--------------------------------------------------------------------- VISBTY1A.62
REAL VISBTY1A.65
& QT(P_FIELD) ! total of cloud water and vapour APC0F405.138
& ,P(P_FIELD) ! pressure of level APC0F405.139
& ,Qs(P_FIELD) ! saturation vapour pressure APC0F405.140
C*L External subroutine called ---------------------------------------- APC0F405.141
EXTERNAL QSAT_WAT APC0F405.142
C*-------------------------------------------------------------------- PC120793.181
C constants for visibility calculation used to be set here but now PC120793.182
C set in MODECK. PC120793.183
C PI needed to set new constants. APC3F304.162
C--------------------------------------------------------------------- VISBTY1A.73
*CALL C_PI 
APC3F304.163
*CALL C_DENSTY APC0F405.143
*CALL C_VISBTY PC120793.184
! APC0F405.144
! Local parameter variables APC0F405.145
! APC0F405.146
REAL APC0F405.147
& OneThird ! 1/3 APC0F405.148
&, RHmax ! Maximum value of relative humidity APC0F405.149
& ! which is allowed to feed into the APC0F405.150
& ! calculation of the 'fog' droplet radius APC0F405.151
&, RHmin ! Minimum value of relative humidity which APC0F405.152
& ! is allowed to feed into the calculation APC0F405.153
& ! of the 'fog' droplet radius APC0F405.154
&, Weight ! Weighting on new value for iterative APC0F405.155
& ! solution of droplet radius APC0F405.156
&, Delta_radius_star! Convergence required for iterative APC0F405.157
& ! solution of droplet radius APC0F405.158
&, N ! Local number density APC0F405.159
&, qt_limit ! Smallest Qt value allowed APC0F405.160
&, radius_star_min ! APC0F405.161
&, radius_star_max ! APC0F405.162
&, radius_star_factor! APC0F405.163
APC0F405.164
INTEGER APC0F405.165
& Niterations ! Maximum number of iteration used to APC0F405.166
& ! estimate the water droplet radius APC0F405.167
PARAMETER ( OneThird = 1.0/3.0 APC0F405.168
&, RHmin = 0.001 APC0F405.169
&, RHmax = 0.99 APC0F405.170
&, Weight = 0.75 APC0F405.171
&, Delta_radius_star = 0.001 APC0F405.172
&, Niterations = 20 APC0F405.173
&, qt_limit = 0.0001 APC0F405.174
&, radius_star_min = 1.0 APC0F405.175
&, radius_star_max = 1000.0 APC0F405.176
&, radius_star_factor = 4.0 ) APC0F405.177
! APC0F405.178
! Local workspace variables APC0F405.179
! APC0F405.180
INTEGER APC0F405.181
& Point ! Loop variable for points APC0F405.182
&, Iteration ! Loop variable iterations used to estimate APC0F405.183
& ! the water droplet radius APC0F405.184
APC0F405.185
REAL APC0F405.186
& m_over_m0 ! Ratio of aerosol mass mixing ratio and APC0F405.187
& ! the standard aerosol mass mixing ratio APC0F405.188
&, RecipVis ! Recipirical of the visibility APC0F405.189
&, radius_dry ! Radius of dry aerosol particle (m) APC0F405.190
&, radius ! Radius of fog droplets (m) APC0F405.191
&, radius_star1 ! Previous estimate of water droplet radius APC0F405.192
& ! divided by the dry radius APC0F405.193
&, radius_star2 ! Current best estimate of water droplet APC0F405.194
& ! radius divided by the dry radius APC0F405.195
&, radius_act ! Activation droplet radius APC0F405.196
&, radius_star_act ! Activation droplet rad divided by dry rad APC0F405.197
&, A ! A0 divided by the dry radius APC0F405.198
&, RH_lim ! Limited RH value (fractional) APC0F405.199
&, Fn ! Value of droplet radius function APC0F405.200
&, Deriv ! Derivative of droplet radius function APC0F405.201
&, radius_star_diff ! Absolute value of radius_star1 minus APC0F405.202
& ! radius_star2 APC0F405.203
&, RHterm ! Relative humidity term in function to be APC0F405.204
& ! minimised to find the droplet radius APC0F405.205
&, qLterm ! Liquid water term in function to be APC0F405.206
& ! minimised to find the droplet radius APC0F405.207
&, RHderiv ! Derivative of relative humidity term APC0F405.208
&, qLderiv ! Derivative of liquid water term APC0F405.209
&, bs ! Width of distribution in total water APC0F405.210
& ! mixing ratio space (kg/kg) APC0F405.211
&, qt_mod ! Modified total water value based on the APC0F405.212
& ! probability of the value occurring APC0F405.213
& ! assuming a triangular distriubtion APC0F405.214
& ! of width bs. APC0F405.215
&, qt_mod_factor ! Factor to multiply bs to modify qt APC0F405.216
! Check Prob is legal APC0F405.217
IF ( Prob .LT. 0.0 .OR. Prob .GT. 1.0 ) THEN APC0F405.218
Write(6,*)"INVALID PROBABILITY VALUE in VISBTY",Prob APC0F405.219
Prob=MIN(MAX(Prob,0.0),1.0) APC0F405.220
ENDIF APC0F405.221
! Create factor to multiply bs by to modify qt APC0F405.222
IF ( Prob .EQ. 0.5 ) THEN APC0F405.223
qt_mod_factor = 0.0 APC0F405.224
ELSE IF ( Prob .GE. 0.0 .AND. Prob .LT. 0.5 ) THEN APC0F405.225
qt_mod_factor = ( 1.0 - SQRT( 2.0 * Prob ) ) APC0F405.226
ELSE IF ( Prob .GE. 0.5 .AND. Prob .LE. 1.0 ) THEN APC0F405.227
qt_mod_factor = - ( 1.0 - SQRT( 2.0 * (1.0-Prob) ) ) APC0F405.228
END IF APC0F405.229
! ---------------------------------------------------------------------- AYY1F404.116
! For the new cloud and precipitation scheme only use the liquid content AYY1F404.117
! 1. Calculate total of water vapour and liquid water contents, P and APC0F405.230
! limit aerosol APC0F405.231
! ---------------------------------------------------------------------- AYY1F404.119
DO Point=1,P_FIELD APC0F405.232
QT(Point) = Q(Point)+QCL(Point) APC0F405.233
END DO AYY1F404.125
APC0F405.234
! Make sure aerosol greater than lower bound APC0F405.235
DO Point=1,P_FIELD APC0F405.236
AEROSOL(Point)=MAX(AEROSOL(Point),AERO0) APC0F405.237
ENDDO VISBTY1A.87
APC0F405.238
! Calculate pressure APC0F405.239
DO Point=1,P_FIELD APC0F405.240
P(Point)=AK+BK*PSTAR(Point) APC0F405.241
ENDDO VISBTY1A.99
APC0F405.242
CALL QSAT_WAT (QS,T,P,P_FIELD) APC0F405.243
APC0F405.244
DO Point = 1 , P_FIELD APC0F405.245
APC0F405.246
!------------------------------------------------------------------- APC0F405.247
!* 2. Calculate the ratio of the aerosol mass mixing ratio to the APC0F405.248
!* standard mass mixing ratio, m_over_m0, and the aerosol number APC0F405.249
!* density, N, the dry radius, radius_dry: APC0F405.250
!* p APC0F405.251
!* (m ) APC0F405.252
!* r = r0 (--) APC0F405.253
!* d (m0) APC0F405.254
!* APC0F405.255
!* APC0F405.256
!* And the activation radius: APC0F405.257
!* APC0F405.258
!* 1/2 APC0F405.259
!* ( 3 ) APC0F405.260
!* ( 3 B0 r ) APC0F405.261
!* r = ( ------d-) APC0F405.262
!* act ( A0 ) APC0F405.263
!* APC0F405.264
!* and A (A0 divided by the dry radius). APC0F405.265
! N.B. AEROSOL is in ug/kg, m in kg/kg APC0F405.266
! If not available, use 10 ug/kg APC0F405.267
!------------------------------------------------------------------- APC0F405.268
APC0F405.269
if (L_MURK) then APC0F405.270
m_over_m0 = max(Aerosol(Point)/m0*1.0E-9, 0.0001) APC0F405.271
else APC0F405.272
m_over_m0 = max(10.0/m0*1.0E-9, 0.0001) APC0F405.273
endif APC0F405.274
APC0F405.275
N = N0 * m_over_m0**(1.0-3*power) APC0F405.276
APC0F405.277
radius_dry = radius0 * (m_over_m0)**power APC0F405.278
A = A0 / radius_dry APC0F405.279
APC0F405.280
radius_act = SQRT( (3 * B0 * radius_dry**3) / A0 ) APC0F405.281
radius_star_act = radius_act/radius_dry APC0F405.282
APC0F405.283
!------------------------------------------------------------- APC0F405.284
!* 3. Calculate the width of the total water APC0F405.285
!* distribution and a modified value of total water, based APC0F405.286
!* on a probability. APC0F405.287
!------------------------------------------------------------- APC0F405.288
APC0F405.289
bs = (1.0-RHcrit) * qs(Point) APC0F405.290
APC0F405.291
qt_mod = MAX( qt_limit, qt(Point)+ qt_mod_factor* bs) APC0F405.292
APC0F405.293
APC0F405.294
!==================================================================== APC0F405.295
!* 4. Use Newton-Raphson to iteratively improve on a first-guess APC0F405.296
!* droplet radius, using the droplet growth equation and the APC0F405.297
!* geometric relation between liquid water and droplet radius. APC0F405.298
!==================================================================== APC0F405.299
!* 4.1 Calculate a first guess relative humidity, qt/qs, but limit it APC0F405.300
!* to be in the range 0.001 -> 0.999. APC0F405.301
!* From this calculate a first-guess normalised radius using a APC0F405.302
!* simplified version of the droplet growth equation: APC0F405.303
!* APC0F405.304
!* 1/3 APC0F405.305
!* ( B0 ) APC0F405.306
!* r = ( 1 - ------ ) APC0F405.307
!* * ( ln(RH) ) APC0F405.308
!* APC0F405.309
!---------------------------------------------------------------------- APC0F405.310
APC0F405.311
RH_lim = MIN( MAX( qt_mod/qs(Point), RHmin ) , RHmax ) APC0F405.312
radius_star2 = (1.0-B0/LOG(RH_lim))**OneThird APC0F405.313
APC0F405.314
!---------------------------------------------------------------------- APC0F405.315
!* 4.2 Initialise the iteration counter, the normalised radius APC0F405.316
!* difference, and the updated normalised radius value. APC0F405.317
!---------------------------------------------------------------------- APC0F405.318
APC0F405.319
Iteration = 0 APC0F405.320
radius_star_diff = 1.0 APC0F405.321
radius_star1 = radius_star2 APC0F405.322
APC0F405.323
Do While ( Iteration .LT. Niterations .AND. APC0F405.324
& radius_star_diff .GT. Delta_radius_star ) APC0F405.325
APC0F405.326
!---------------------------------------------------------------------- APC0F405.327
!* 4.3 Update the iteration counter and the normalised radius value. APC0F405.328
!---------------------------------------------------------------------- APC0F405.329
APC0F405.330
Iteration = Iteration + 1 APC0F405.331
radius_star1 = Weight * radius_star2 APC0F405.332
& + ( 1.0 - Weight ) * radius_star1 APC0F405.333
APC0F405.334
!---------------------------------------------------------------------- APC0F405.335
!* 4.4 Calculate the relative humidity term: APC0F405.336
!* APC0F405.337
!* ( A B0 ) APC0F405.338
!* RHterm = exp( -- - ------ ) APC0F405.339
!* ( r 3 ) APC0F405.340
!* ( * r - 1 ) APC0F405.341
!* ( * ) APC0F405.342
!* APC0F405.343
!* and its derivative with respect to the normalised radius: APC0F405.344
!* APC0F405.345
!* ( 2 ) APC0F405.346
!* ( A 3 B0 r ) APC0F405.347
!* RHderiv = ( - -- + -------*- 2 ) * RHterm APC0F405.348
!* ( 2 ( 3 ) ) APC0F405.349
!* ( r ( r - 1 ) ) APC0F405.350
!* ( * ( * ) ) APC0F405.351
!* APC0F405.352
!---------------------------------------------------------------------- APC0F405.353
APC0F405.354
If ( radius_star1 .LT. radius_star_act ) then APC0F405.355
RHterm = EXP( A/radius_star1 APC0F405.356
& - B0/(radius_star1**3-1.0) )* qs(Point) APC0F405.357
RHderiv = - RHterm * ( -A/(radius_star1**2) APC0F405.358
& + (3.0*B0*radius_star1**2) APC0F405.359
& /(radius_star1**3-1.0)**2 ) APC0F405.360
Else APC0F405.361
RHterm = EXP( A/radius_star_act APC0F405.362
& - B0/(radius_star_act**3-1.0) ) * qs(Point) APC0F405.363
RHderiv = 0.0 APC0F405.364
Endif APC0F405.365
APC0F405.366
APC0F405.367
!---------------------------------------------------------------------- APC0F405.368
!* 4.5 Calculate the liquid water mixing ratio term: APC0F405.369
!* APC0F405.370
!* APC0F405.371
!* 4 3 ( 3 ) APC0F405.372
!* qLterm = - Pi rho_w N r ( r - 1 ) APC0F405.373
!* 3 d ( * ) APC0F405.374
!* APC0F405.375
!* and its derivative with respect to the normalised radius: APC0F405.376
!* APC0F405.377
!* 3 2 APC0F405.378
!* qLderiv = 4 Pi rho_w N r r APC0F405.379
!* d * APC0F405.380
!* APC0F405.381
!---------------------------------------------------------------------- APC0F405.382
APC0F405.383
qLterm = N * FourThirds * Pi * RHO_WATER * radius_dry**3 APC0F405.384
& * ( radius_star1**3 - 1.0 ) APC0F405.385
qLderiv = - N * 4.0 * Pi * RHO_WATER APC0F405.386
& * radius_dry**3 * radius_star1**2 APC0F405.387
APC0F405.388
!---------------------------------------------------------------------- APC0F405.389
!* 4.6 Calculate the function, Fn, and its derivative, Deriv, and APC0F405.390
!* an improved estimate of the normalised radius, APC0F405.391
!* using Newton Raphson: APC0F405.392
!* APC0F405.393
!* Fn = qt - RHterm - qLterm APC0F405.394
!* APC0F405.395
!* Deriv = RHderiv + qLderiv APC0F405.396
!* APC0F405.397
!* Fn APC0F405.398
!* r = r - ----- APC0F405.399
!* * new * Deriv APC0F405.400
!* APC0F405.401
!* The new estimate of the normalised radius is limited lie between APC0F405.402
!* prescribed maximum and minimum values and within a factor of the APC0F405.403
!* previous value to ensure that the soultion does not diverge. APC0F405.404
!---------------------------------------------------------------------- APC0F405.405
APC0F405.406
Fn = qt_mod - RHterm - qLterm APC0F405.407
Deriv = RHderiv + qLderiv APC0F405.408
APC0F405.409
radius_star2 = radius_star1 - Fn/Deriv APC0F405.410
APC0F405.411
IF ( radius_star2 .LT. radius_star_min ) APC0F405.412
& radius_star2 = radius_star_min APC0F405.413
IF ( radius_star2 .GT. radius_star_max ) APC0F405.414
& radius_star2 = radius_star_max APC0F405.415
IF ( radius_star2 .GT. radius_star_factor * radius_star1 ) APC0F405.416
& radius_star2 = radius_star_factor * radius_star1 APC0F405.417
IF ( radius_star2 .LT. radius_star1 / radius_star_factor ) APC0F405.418
& radius_star2 = radius_star1 / radius_star_factor APC0F405.419
APC0F405.420
!--------------------------------------------------------------------- APC0F405.421
!* 4.7 Calculate difference between the old and the new values of the APC0F405.422
!* normalised radius. APC0F405.423
!--------------------------------------------------------------------- APC0F405.424
APC0F405.425
radius_star_diff = ABS( radius_star1 - radius_star2 ) APC0F405.426
APC0F405.427
END DO APC0F405.428
APC0F405.429
!--------------------------------------------------------------------- APC0F405.430
!* 5. Calculate the radius from the final normalised radius. APC0F405.431
!--------------------------------------------------------------------- APC0F405.432
APC0F405.433
radius = radius_star2 * radius_dry APC0F405.434
APC0F405.435
!--------------------------------------------------------------------- APC0F405.436
!* 6. Calculate the visibility, Vis, using the equation: APC0F405.437
!* APC0F405.438
!* ln(liminal contrast) APC0F405.439
!* Vis = -------------2------ APC0F405.440
!* Beta0 N r APC0F405.441
!* APC0F405.442
!* (An extra term RecipVisAir is included in the recipical of APC0F405.443
!* visibility to limit visibilities to 100km in clean air). APC0F405.444
!--------------------------------------------------------------------- APC0F405.445
APC0F405.446
RecipVis = (N * radius**2) / VisFactor + RecipVisAir APC0F405.447
Visibility(Point) = 1/RecipVis APC0F405.448
APC0F405.449
END DO APC0F405.450
APC0F405.451
RETURN APC3F304.407
END APC3F304.408
C APC3F304.409
CLL SUBROUTINE VISTOQT ----------------------------------------------- APC0F405.452
CLL APC3F304.411
CLL PURPOSE: APC3F304.412
CLL Invert relationship between aerosol, visibility and water APC3F304.413
CLL content APC0F405.453
CLL This is needed for fog probability calculation. APC3F304.415
! Since 4.5, adopted NIMROD based code: APC0F405.454
CLL APC3F304.416
CLL APC3F304.418
CLL Model Modification history from model version 3.0: APC3F304.425
CLL version Date APC3F304.426
CLL 3.4 07/06/94 First written. Programmer Pete Clark. APC3F304.427
CLL 4.0 05/09/95 Diagnosed equivalent RH constrained to be >1%. APC0F400.19
CLL Lower limit to aerosol introduced. Pete Clark. APC0F400.20
! 4.5 30/04/98 NIMROD code adopted. APC0F405.455
! Pete Clark responsible for UM implementation of APC0F405.456
! Bruce Wright's NIMROD code. APC0F405.457
CLL APC3F304.428
CLL Programming standard: Unified Model Documentation Paper No 3, APC3F304.429
CLL Version 7, dated 11/3/93. APC3F304.430
CLL APC3F304.431
CLL Logical components covered : APC3F304.432
CLL APC3F304.433
CLL Project task: APC3F304.434
CLL APC3F304.435
CLL External documentation APC3F304.436
CLL Forecasting Research Scientific Paper NO.4 APC3F304.437
CLL Diagnosis of visibility in the UK Met Office Mesoscale Model APC3F304.438
CLL and the use of a visibility analysis to constrain initial APC3F304.439
CLL conditions. SP Ballard, BJ Wright, BW Golding 1992 APC3F304.440
! Wright, B. J., 1997: Improvements to the Nimrod Visibility APC0F405.458
! Analysis/Forecast System. FR-Div. Tech. Rep., No. 217. APC0F405.459
! Wright, B. J., 1997: A New Visibility Analysis/Forecast System APC0F405.460
! for Nimrod. Met. Office FR Tech Rep., No. 222. APC0F405.461
CLL APC3F304.441
CLLEND---------------------------------------------------------------- APC3F304.442
C APC3F304.443
C*L Arguments:------------------------------------------------------- APC3F304.444
SUBROUTINE VISTOQT 1APC0F405.462
& (VISIBILITY !INPUT APC3F304.446
& ,Qs !INPUT APC0F405.463
& ,AEROSOL !INPUT APC3F304.448
& ,L_MURK !INPUT APC3F304.449
& ,Npoints !INPUT APC0F405.464
& ,qt ) !OUTPUT APC0F405.465
IMPLICIT NONE APC3F304.451
C--------------------------------------------------------------------- APC3F304.452
C Workspace usage:---------------------------------------------------- APC3F304.453
C None APC3F304.454
C*-------------------------------------------------------------------- APC3F304.455
C*L------------------------------------------------------------------- APC3F304.456
C input variables----------------------------------------------------- APC3F304.457
C--------------------------------------------------------------------- APC3F304.458
INTEGER APC3F304.459
& Npoints ! IN NO. points in field. APC0F405.466
REAL APC3F304.461
& VISIBILITY ! IN visibility APC0F405.467
! NB Original code had Visibility APC0F405.468
! as an array - this feature was APC0F405.469
! not used so has been removed APC0F405.470
& ,Qs(Npoints) ! Saturated humidity mixing ratio APC0F405.471
& ,AEROSOL(Npoints) ! IN Aerosol mixing ratio(ug/kg) APC0F405.472
LOGICAL APC3F304.464
& L_MURK ! IN : Aerosol present APC3F304.465
C--------------------------------------------------------------------- APC3F304.466
C output variables---------------------------------------------------- APC3F304.467
C--------------------------------------------------------------------- APC3F304.468
REAL APC3F304.469
& qt(Npoints) ! OUT Total water mixing ratio (kg/kg) APC0F405.473
C*-------------------------------------------------------------------- APC3F304.471
C*L------------------------------------------------------------------- APC3F304.472
C--------------------------------------------------------------------- APC3F304.474
C*-------------------------------------------------------------------- APC3F304.486
C constants for visibility calculation used to be set here but now APC3F304.487
C set in COMDECK. APC3F304.488
C PI needed to set new constants. APC3F304.489
C--------------------------------------------------------------------- APC3F304.490
*CALL C_PI APC3F304.491
*CALL C_R_CP APC0F405.474
*CALL C_EPSLON APC0F405.475
*CALL C_DENSTY APC0F405.476
*CALL C_VISBTY APC3F304.492
! APC0F405.477
! Local parameter variables APC0F405.478
! APC0F405.479
REAL APC0F400.21
& qt_limit ! Smallest total water mixing ratio value allowed APC0F405.480
PARAMETER ( qt_limit = 0.0001 ) APC0F405.481
! APC0F405.482
! Local workspace variables APC0F405.483
! APC0F405.484
INTEGER APC0F405.485
& Point ! Loop variable for points APC0F405.486
APC0F405.487
REAL APC0F405.488
& qL ! Liquid water mixing ratio (Kg/Kg). APC0F405.489
&, radius_dry ! Dry particle radius for aerosol (m) APC0F405.490
&, radius ! Radius of fog droplets (m) APC0F405.491
&, radius_star ! Water droplet radius divided by dry radius APC0F405.492
&, radius_act ! Activation droplet radius APC0F405.493
&, radius_star_act ! Activation droplet radius divided by the dry APC0F405.494
&, radius_star_used ! Water droplet radius divided by the dry APC0F405.495
& ! radius actually used for the relative APC0F405.496
& ! humidity calculation APC0F405.497
&, RH ! Relative humidity derived from visibility APC0F405.498
&, A ! A0 divided by the dry radius APC0F405.499
&, m_over_m0 ! Ratio of the aerosol mass mixing ratio and APC0F405.500
& ! the standard aerosol mass mixing ratio APC0F405.501
&, N ! Number density of aerosol particles (/m3) APC0F405.502
APC0F405.503
APC0F405.504
!** APC0F405.505
APC0F405.506
Do Point = 1 , Npoints APC0F405.507
APC0F405.508
!--------------------------------------------------------------------- APC0F405.509
!* 1. Calculate the ratio of the aerosol mass mixing ratio to the APC0F405.510
!* standard mass mixing ratio, m_over_m0, and the aerosol number APC0F405.511
!* density, N: APC0F405.512
!* APC0F405.513
!* (1-3p) APC0F405.514
!* (m ) APC0F405.515
!* N = N0 (--) APC0F405.516
!* (m0) APC0F405.517
!* APC0F405.518
!* And the dry radius, radius_dry: APC0F405.519
!* p APC0F405.520
!* (m ) APC0F405.521
!* r = r0 (--) APC0F405.522
!* d (m0) APC0F405.523
!* APC0F405.524
!* And A (A0 divided by the dry radius). APC0F405.525
!* APC0F405.526
!* And the activation radius: APC0F405.527
!* APC0F405.528
!* 1/2 APC0F405.529
!* ( 3 ) APC0F405.530
!* ( 3 B0 r ) APC0F405.531
!* r = ( ------d-) APC0F405.532
!* act ( A0 ) APC0F405.533
!* APC0F405.534
! N.B. AEROSOL is in ug/kg, m in kg/kg APC0F405.535
! If not available, use 10 ug/kg APC0F405.536
!------------------------------------------------------------------- APC0F405.537
APC0F405.538
if (L_MURK) then APC0F405.539
m_over_m0 = max(Aerosol(Point)/m0*1.0E-9, 0.0001) APC0F405.540
else APC0F405.541
m_over_m0 = max(10.0/m0*1.0E-9, 0.0001) APC0F405.542
endif APC0F405.543
APC0F405.544
N = N0 * m_over_m0**(1.0-3*power) APC0F405.545
APC0F405.546
radius_dry = radius0 * (m_over_m0)**power APC0F405.547
A = A0 / radius_dry APC0F405.548
APC0F405.549
radius_act = SQRT( (3 * B0 * radius_dry**3) / A0 ) APC0F405.550
radius_star_act = radius_act/radius_dry APC0F405.551
APC0F405.552
!---------------------------------------------------------------------- APC0F405.553
!* 2. Calculate a water droplet radius, from the visibility: APC0F405.554
!* APC0F405.555
!* 1/2 APC0F405.556
!* ( ln( epsilon ) ) APC0F405.557
!* r = (---------------) APC0F405.558
!* ( Vis N Beta0 ) APC0F405.559
!* APC0F405.560
!* (An extra term RecipVisAir is included in the recipical of APC0F405.561
!* visibility to limit visibilities to 100km in clean air). APC0F405.562
!---------------------------------------------------------------------- APC0F405.563
APC0F405.564
radius = SQRT( (VisFactor/N) * APC0F405.565
& ((1.0/Visibility) - RecipVisAir) ) APC0F405.566
APC0F405.567
!---------------------------------------------------------------------- APC0F405.568
!* 3. Provided the diagnosed radius is greater than the dry radius, APC0F405.569
!* calculate the normalised droplet radius, and the saturated APC0F405.570
!* humidity mixing ratio. APC0F405.571
!---------------------------------------------------------------------- APC0F405.572
APC0F405.573
If ( radius .GT. radius_dry ) then APC0F405.574
APC0F405.575
radius_star = radius / radius_dry APC0F405.576
APC0F405.577
!---------------------------------------------------------------------- APC0F405.578
!* 5. Calculate the corresponding liquid water mixing ratio: APC0F405.579
!* APC0F405.580
!* APC0F405.581
!* 4 ( 3 3 ) APC0F405.582
!* qL = - Pi rho_w N ( r - r ) APC0F405.583
!* 3 ( d ) APC0F405.584
!* APC0F405.585
!---------------------------------------------------------------------- APC0F405.586
APC0F405.587
qL = FourThirds * Pi * rho_water * N * APC0F405.588
& ( radius**3 - radius_dry**3 ) APC0F405.589
APC0F405.590
!---------------------------------------------------------------------- APC0F405.591
!* 6. Calculate the relative humidity: APC0F405.592
!* APC0F405.593
!* ( A B0 ) APC0F405.594
!* RH = exp( -- - ------ ) APC0F405.595
!* ( r 3 ) APC0F405.596
!* ( * r - 1 ) APC0F405.597
!* ( * ) APC0F405.598
!* APC0F405.599
!---------------------------------------------------------------------- APC0F405.600
APC0F405.601
If ( radius_star .LT. radius_star_act ) then APC0F405.602
RH = EXP( A/radius_star APC0F405.603
& - B0 /( radius_star **3 - 1.0 ) ) APC0F405.604
Else APC0F405.605
RH = EXP( A/radius_star_act APC0F405.606
& - B0 /( radius_star_act **3 - 1.0 ) ) APC0F405.607
Endif APC0F405.608
APC0F405.609
!---------------------------------------------------------------------- APC0F405.610
!* 7. Calculate the total water mixing ratio: qt = RH * qs(T) + qL APC0F405.611
!---------------------------------------------------------------------- APC0F405.612
APC0F405.613
qt(Point) = MAX( RH * qs(Point) + qL, qt_limit ) APC0F405.614
APC0F405.615
!---------------------------------------------------------------------- APC0F405.616
!* 8. If the droplet radius is less than the dry radius, then set the APC0F405.617
!* total water mixing ratio to the minimum value. APC0F405.618
!---------------------------------------------------------------------- APC0F405.619
APC0F405.620
Else APC0F405.621
APC0F405.622
qt(Point) = qt_limit APC0F405.623
APC0F405.624
End if APC0F405.625
APC0F405.626
APC0F405.627
End Do APC0F405.628
APC0F405.629
RETURN PC120793.186
END PC120793.187
C PC120793.188
CLL SUBROUTINE FOG_FR------------------------------------------------ PC120793.189
CLL PC120793.190
CLL Purpose: Calculates fog fraction, using the large scale cloud PC120793.191
CLL scheme. The fog fraction is similar to the cloud fraction PC120793.192
CLL except it records the fraction of a grid box with RH PC120793.193
CLL greater than that required for the critical visibility PC120793.194
CLL (e.g 1 km). APC0F405.630
!LL Since 4.5, adopted NIMROD based code: APC0F405.631
!LL Calculates the fraction of a gridsquare with visibility APC0F405.632
!LL less than threshold, Vis_thresh, given the total water APC0F405.633
!LL mixing ratio, qt, temperature, T, pressure p, and the APC0F405.634
!LL aerosol mass mixing ratio, m, assuming a triangular APC0F405.635
!LL distribution of states about the median, characterised by APC0F405.636
!LL a critical relative humdity value, RHcrit. APC0F405.637
CLL NB: Throughout, levels are counted from the bottom up, PC120793.225
CLL i.e. the lowest level under consideration is level 1, the PC120793.226
CLL next lowest level 2, and so on. PC120793.227
CLL PC120793.228
CLL Suitable for single-column use. PC120793.229
CLL PC120793.230
CLL Model Modification history: PC120793.231
CLL version Date PC120793.232
CLL PC120793.233
CLL 3.2 04/05/93 Created by Pete Clark PC120793.234
!LL 3.4 04/08/95 LS_CLD replaced by GLUE_CLD. Andrew Bushell. AYY2F400.215
CLL 4.2 Oct. 96 T3E migration: *DEF CRAY removed (dynamic GSS9F402.35
CLL allocation now unconditional) GSS9F402.36
CLL S.J.Swarbrick GSS9F402.37
! 4.4 01/07/97 Calculation is now based on liquid water AYY1F404.126
! and not on ice. Calculates liquid fog fraction. AYY1F404.127
! This scheme is severely tied to the ideas behind AYY1F404.128
! the 1A cloud scheme. Damian Wilson. AYY1F404.129
!LL 4.5 30/04/98 NIMROD code adopted. Calculation is still based APC0F405.638
!LL on liquid water and not ice, but arguments changed APC0F405.639
!LL to pass T, q, qcl and qcf separately. This is to APC0F405.640
!LL make code independent of cloud scheme and capable APC0F405.641
!LL of future development to include ice properly. APC0F405.642
!LL Pete Clark responsible for UM implementation of APC0F405.643
!LL Bruce Wright's NIMROD code. APC0F405.644
CLL PC120793.235
CLL Programming standard: Unified Model Documentation Paper No 3, PC120793.236
CLL Version 5, dated 08/12/92. PC120793.237
CLL PC120793.238
CLL Documentation: APC0F405.645
!LL Wright, B. J., 1997: Improvements to the Nimrod Visibility APC0F405.646
!LL Analysis/Forecast System. FR-Div. Tech. Rep., No. 217. APC0F405.647
!LL Wright, B. J., 1997: A New Visibility Analysis/Forecast System APC0F405.648
!LL for Nimrod. Met. Office FR Tech Rep., No. 222. APC0F405.649
CLL PC120793.240
CLLEND---------------------------------------------------------------- PC120793.241
C PC120793.242
C*L PC120793.243
C*LArguments:--------------------------------------------------------- PC120793.244
SUBROUTINE FOG_FR( 2,2PC120793.245
+ AK,BK,PSTAR,RHCRIT,LEVELS,POINTS,PFIELD, PC120793.246
& T,AEROSOL,L_MURK,Q,QCL,QCF,VIS,FF,NVIS, APC0F405.650
& ERROR APC3F304.538
+) PC120793.249
IMPLICIT NONE PC120793.250
INTEGER PC120793.251
+ LEVELS ! IN No. of levels being processed. PC120793.252
+,POINTS ! IN No. of gridpoints being processed. PC120793.253
+,PFIELD ! IN No. of points in global field (at one PC120793.254
C ! vertical level). PC120793.255
&,NVIS ! IN No. of visibility thresholds APC0F405.651
REAL PC120793.256
+ PSTAR(PFIELD) ! IN Surface pressure (Pa). PC120793.257
+,RHCRIT(LEVELS) ! IN Critical relative humidity. See the PC120793.258
C ! the paragraph incorporating eqs P292.11 PC120793.259
C ! to P292.14; the values need to be tuned PC120793.260
C ! for the given set of levels. PC120793.261
+,AK(LEVELS) ! IN Hybrid "A" co-ordinate. PC120793.262
+,BK(LEVELS) ! IN Hybrid "B" co-ordinate. PC120793.263
+,Q(PFIELD,LEVELS) ! IN Specific Humidity APC0F405.652
C ! (kg per kg air). PC120793.265
&,QCL(PFIELD,LEVELS) ! Cloud liquid water content at APC0F405.653
C ! processed levels (kg per kg air). APC0F405.654
&,QCF(PFIELD,LEVELS) ! Cloud ice content at processed levels APC0F405.655
C ! (kg per kg air). APC0F405.656
&,T(PFIELD,LEVELS) ! IN Temperature (K). APC0F405.657
&,AEROSOL(PFIELD,LEVELS) ! IN Aerosol mixing ratio(ug/kg) APC3F304.539
&,VIS(NVIS) ! Visibility thresholds APC0F405.658
LOGICAL APC3F304.541
& L_MURK ! IN : Aerosol present APC3F304.542
APC3F304.543
REAL PC120793.268
+ FF(PFIELD,LEVELS,NVIS) ! OUT Vis prob at processed levels APC0F405.659
C ! (decimal fraction). PC120793.270
INTEGER ERROR ! OUT 0 if OK; 1 if bad arguments. PC120793.271
C PC120793.272
C*-------------------------------------------------------------------- PC120793.273
C*L Workspace usage---------------------------------------------------- PC120793.274
REAL ! "Automatic" arrays on Cray. PC120793.277
& P(POINTS) APC0F405.660
&,QT(POINTS) ! total of cloud water and vapour APC0F405.661
&,QS(POINTS) ! Saturated spec humidity for temp T APC0F405.662
&,qt_thresh(POINTS) ! modified qt APC0F405.663
&,bs APC0F405.664
C*L External subroutine called ---------------------------------------- PC120793.299
EXTERNAL QSAT_WAT,VISTOQT APC0F405.665
C* Local, including SAVE'd, storage------------------------------------ PC120793.301
C PC120793.302
INTEGER K,I,J ! Loop counters: K - vertical level index. APC0F405.666
C ! I - horizontal field index. PC120793.308
! J - Vis threshold index. APC0F405.667
C*-------------------------------------------------------------------- PC120793.309
C* Local and other physical constants---------------------------------- PC120793.310
*CALL C_PI APC3F304.550
*CALL C_VISBTY PC120793.315
! AYY1F404.135
! (a) Scalars effectively expanded to workspace by the Cray (using AYY1F404.136
! vector registers). AYY1F404.137
C----------------------------------------------------------------------- PC120793.326
C Check input arguments for potential over-writing problems. PC120793.327
C----------------------------------------------------------------------- PC120793.328
ERROR=0 PC120793.329
IF(POINTS.GT.PFIELD)THEN PC120793.330
ERROR=1 PC120793.331
GOTO1000 PC120793.332
ENDIF PC120793.333
C PC120793.334
C PC120793.345
C----------------------------------------------------------------------- PC120793.346
CL Subroutine structure : PC120793.347
CL Loop round levels to be processed. PC120793.348
C----------------------------------------------------------------------- PC120793.349
C PC120793.350
DO K=1,LEVELS PC120793.351
C PC120793.352
C----------------------------------------------------------------------- PC120793.353
CL 1. Calculate Pressure and initialise temporary arrays PC120793.354
C----------------------------------------------------------------------- PC120793.355
C PC120793.356
DO I=1,POINTS PC120793.357
P(I)=AK(K)+PSTAR(I)*BK(K) PC120793.358
QT(I)=Q(I,K)+QCL(I,K) APC0F405.668
ENDDO ! Loop over points PC120793.362
APC0F405.669
!----------------------------------------------------------------------- APC0F405.670
!* 2. Calculate total water threshold corresponding to visibility APC0F405.671
! Since Qs is needed more than once, pre-calculate and pass it APC0F405.672
!----------------------------------------------------------------------- APC0F405.673
APC0F405.674
CALL QSAT_WAT (Qs,T(1,K),P,POINTS) APC0F405.675
APC0F405.676
DO J=1,NVIS APC0F405.677
APC0F405.678
Call VISTOQT( VIS(J), Qs, AEROSOL(1,K), L_MURK, APC0F405.679
& points, qt_thresh ) APC0F405.680
APC0F405.681
APC0F405.682
!----------------------------------------------------------------------- APC0F405.683
!* 3. Calculate the width of the distribution in total water space, bs: APC0F405.684
!* APC0F405.685
!* bs = ( 1 - RHcrit ) * qs(T) APC0F405.686
!* APC0F405.687
!----------------------------------------------------------------------- APC0F405.688
APC0F405.689
Do I = 1 , points APC0F405.690
APC0F405.691
bs = (1.0-RHcrit(K)) * qs(I) APC0F405.692
APC0F405.693
!======================================================================= APC0F405.694
!* 4. Calculate the fraction of states in a triangular APC0F405.695
!* distribution which exceed the total water threshold. APC0F405.696
!======================================================================= APC0F405.697
APC0F405.698
!----------------------------------------------------------------------- APC0F405.699
!* 4.1 If total water threshold value is less than the total water value APC0F405.700
!* minus the width of the distribution, then all of the states have APC0F405.701
!* a total water value exceeding the threshold, so set the APC0F405.702
!* visibility fraction to 1.0 APC0F405.703
!----------------------------------------------------------------------- APC0F405.704
APC0F405.705
if ( qt_thresh(I) .LE. qt(I)-bs ) then APC0F405.706
APC0F405.707
FF(I,K,J) = 1.0 APC0F405.708
APC0F405.709
!----------------------------------------------------------------------- APC0F405.710
!* 4.2 If total water threshold value is greater than the total water APC0F405.711
!* value minus the width of the distribution, but less than the APC0F405.712
!* total water value then the visibility fraction, VF, is given by: APC0F405.713
!* APC0F405.714
!* 2 APC0F405.715
!* ( qt - qt + bs ) APC0F405.716
!* VF = 1.0 - 0.5 * ( thresh ) APC0F405.717
!* ( ------------------- ) APC0F405.718
!* ( bs ) APC0F405.719
!* APC0F405.720
!----------------------------------------------------------------------- APC0F405.721
APC0F405.722
Else if ( qt_thresh(I) .GT. qt(I)-bs .AND. APC0F405.723
& qt_thresh(I) .LE. qt(I) ) then APC0F405.724
APC0F405.725
FF(I,K,J) = 1.0 - 0.5 * APC0F405.726
& (( qt_thresh(I) - qt(I) + bs )/ bs)**2 APC0F405.727
APC0F405.728
!----------------------------------------------------------------------- APC0F405.729
!* 4.3 If total water threshold value is greater than the total water APC0F405.730
!* value, but less than the total water value plus the width of the APC0F405.731
!* distribution, then the visibility fraction, VF, is given by: APC0F405.732
!* APC0F405.733
!* 2 APC0F405.734
!* ( qt + bs - qt ) APC0F405.735
!* VF = 0.5 * ( thresh ) APC0F405.736
!* ( ------------------- ) APC0F405.737
!* ( bs ) APC0F405.738
!* APC0F405.739
!----------------------------------------------------------------------- APC0F405.740
APC0F405.741
Else if ( qt_thresh(I) .GT. qt(I) .AND. APC0F405.742
& qt_thresh(I) .LE. qt(I)+bs ) then APC0F405.743
APC0F405.744
FF(I,K,J)= 0.5 * (( qt(I) + bs - qt_thresh(I))/bs)**2 APC0F405.745
APC0F405.746
!----------------------------------------------------------------------- APC0F405.747
!* 4.4 If total water threshold value is greater than the total water APC0F405.748
!* value plus the width of the distribution, then non of the states APC0F405.749
!* have a total water value exceeding the threshold, so set the APC0F405.750
!* visibility fraction to 0.0 APC0F405.751
!----------------------------------------------------------------------- APC0F405.752
APC0F405.753
Else APC0F405.754
APC0F405.755
FF(I,K,J) = 0.0 APC0F405.756
APC0F405.757
End if APC0F405.758
APC0F405.759
End Do ! Loop over Points I APC0F405.760
APC0F405.761
End Do ! Loop over VIS J APC0F405.762
APC0F405.763
ENDDO ! Loop over levels PC120793.388
C PC120793.389
1000 CONTINUE ! Error exit PC120793.394
RETURN VISBTY1A.106
END VISBTY1A.107
*ENDIF VISBTY1A.108