Subroutine STATSTEP( 1,14S_STATSP.38
! IN S_STATSP.39
& points, nlevs, nwet, ntrop, S_STATSP.40
! S_STATSP.41
& deltan, px, py, daysteps, stepcount, dayno, S_STATSP.42
& tr, vnr, vpr, qr, wr, tbar, tsd, tdash, dbar, dsd, ddash, S_STATSP.43
& vnbar, vpbar, vnsd, wbar, wsd, ctbar, ctsd, at, cdbar, S_STATSP.44
& cdsd, ad, cvnbar, cvnsd, avn, cwbar, cwsd, aw, S_STATSP.45
& press, rpress, u, v, t, q, prinstat, S_STATSP.46
& daycount, timestep, iv, ntab, iy, idum) S_STATSP.47
Implicit none S_STATSP.48
C S_STATSP.49
S_STATSP.50
S_STATSP.51
C S_STATSP.52
C--------------------------------------------------------------------- S_STATSP.53
C Arguments S_STATSP.54
C--------------------------------------------------------------------- S_STATSP.55
C S_STATSP.56
Integer S_STATSP.57
& points ! IN no of model columns S_STATSP.58
& ,nlevs ! IN no of levels. S_STATSP.59
& ,nwet ! IN no of model levels in which Q is S_STATSP.60
! set S_STATSP.61
& ,ntrop ! IN Max number of levels in the S_STATSP.62
! troposphere S_STATSP.63
& ,daycount ! IN Daynumber (1 represents S_STATSP.64
! 1st january) S_STATSP.65
& ,dayno ! IN Daynumber relative to winter S_STATSP.66
! solstice S_STATSP.67
& ,daysteps ! IN No. of timesteps in 1 day S_STATSP.68
& ,stepcount ! IN Timestep counter S_STATSP.69
& ,ntab ! IN Dimension of array used in S_STATSP.70
! random generator S_STATSP.71
& ,iv(ntab),iy,idum ! IN state of number generator saved S_STATSP.72
! on tape for continuation run S_STATSP.73
S_STATSP.74
Logical S_STATSP.75
& prinstat ! T if printout of increments S_STATSP.76
! due to statistical forcing S_STATSP.77
! required each timestep S_STATSP.78
Real S_STATSP.79
& ad(points,nwet-1) ! IN term a of equation 2.22 for dew S_STATSP.80
& ,at(points,nlevs-1) ! pt depression and temp. S_STATSP.81
& ,avn(points,nlevs-1) ! IN term a of equation 2.22 for S_STATSP.82
& ,aw(points,ntrop-1) ! horiz. and vertical velocities S_STATSP.83
S_STATSP.84
& ,cdbar(points,nwet) ! Mean and SD of random variable S_STATSP.85
& ,cdsd(points,nwet) ! for dew point depression S_STATSP.86
& ,ctbar(points,nlevs) ! IN Mean and SD of random variable S_STATSP.87
& ,ctsd(points,nlevs) ! for temp. S_STATSP.88
S_STATSP.89
& ,cvnbar(points,nlevs) ! IN Mean and SD of random variable S_STATSP.90
& ,cvnsd(points,nlevs) ! for velocity VN S_STATSP.91
& ,cwbar(points,ntrop) S_STATSP.92
& ,cwsd(points,ntrop) ! IN Mean and SD of random variable S_STATSP.93
! for vertical velocity S_STATSP.94
& ,dbar(points,nwet) ! IN Mean and SD dewpoint S_STATSP.95
& ,dsd(points,nwet) ! depression at daycount days S_STATSP.96
! from winter solstice (K) S_STATSP.97
& ,ddash(points,nwet) ! IN Dew pt. corrections S_STATSP.98
& ,deltan(points) ! IN Radius of area (m) S_STATSP.99
& ,dq1 ! IN Spec. humidity differences S_STATSP.100
& ,dq2 ! (Kg Kg^-1) S_STATSP.101
& ,dt1 ! IN Temp. differences (K) S_STATSP.102
& ,dt2 S_STATSP.103
& ,press(points,nlevs) ! IN Pressure coordinates (Pa) S_STATSP.104
& ,px(points,ntrop) ! IN Reciprocal log functions for S_STATSP.105
& ,py(points,ntrop-1) ! calc. of vert. advection S_STATSP.106
! used in eqns 2.12 and 2.13 S_STATSP.107
& ,q(points,nwet) ! OUT Specific humidity (Kg Kg^-1) S_STATSP.108
& ,qr(points,nwet,2) ! INOUT Randomly sampled humidity S_STATSP.109
! (Kg Kg^-1) S_STATSP.110
& ,rpress(points,nlevs) ! IN Reciprocal pressure for sigma S_STATSP.111
! levels (HPa or mb ^-1) S_STATSP.112
& ,t(points,nlevs) ! OUT Temp (K) S_STATSP.113
& ,tbar(points,nlevs) ! IN Mean and SD temperature at S_STATSP.114
! daycount days from S_STATSP.115
! winter solstice (K) S_STATSP.116
& ,tdash(points,nlevs) ! IN Temp. corrections (K) S_STATSP.117
& ,timestep ! IN model timestep (s) S_STATSP.118
& ,tsd(points,nlevs) ! IN SD of temp. at daycount days S_STATSP.119
! from winter solstice (K) S_STATSP.120
& ,tr(points,nlevs,2) ! INOUT Randomly sampled temp. (K) S_STATSP.121
& ,u(points,nlevs) S_STATSP.122
& ,v(points,nlevs) ! OUT Zonal and meridional winds S_STATSP.123
! (m s^-1) S_STATSP.124
& ,vnbar(points,nlevs) ! IN Mean and SD velocity VN at S_STATSP.125
! daycount days from S_STATSP.126
! winter solstice (m s^-1) S_STATSP.127
& ,vnr(points,nlevs,2) ! INOUT Randomly sampled horizontal S_STATSP.128
! velocity (m s^-1) S_STATSP.129
& ,vnsd(points,nlevs) ! IN Mean and SD velocity VN at S_STATSP.130
! daycount days from S_STATSP.131
! winter solstice (m s^-1) S_STATSP.132
& ,vpbar(points,nlevs) ! IN Mean velocity VP at S_STATSP.133
! daycount days from S_STATSP.134
! winter solstice (m s^-1) S_STATSP.135
& ,vpr(points,nlevs,2) ! INOUT Randomly sampled horizontal S_STATSP.136
! velocity (m s^-1) S_STATSP.137
& ,wbar(points,ntrop) S_STATSP.138
& ,wsd(points,ntrop) ! IN Mean and SD vertical S_STATSP.139
! velocity at daycount days S_STATSP.140
! from winter solstice (mb s^-1) S_STATSP.141
& ,wr(points,ntrop,2) ! INOUT Randomly sampled vertical S_STATSP.142
! velocity (mb s^-1) S_STATSP.143
C S_STATSP.144
C--------------------------------------------------------------------- S_STATSP.145
C Local variables S_STATSP.146
C--------------------------------------------------------------------- S_STATSP.147
C S_STATSP.148
Integer S_STATSP.149
& i, i1, j1, k, l ! Loop counters S_STATSP.150
Real S_STATSP.151
& cdd ! Randomly sampled variables for S_STATSP.152
& ,ct ! temp. and dew pt. depression S_STATSP.153
& ,cvn ! Randomly sampled variables for S_STATSP.154
& ,cw ! horizontal and vertical velocity S_STATSP.155
& ,d0 ! Randomly sampled dew pt. S_STATSP.156
! depression for 1st level S_STATSP.157
& ,dewpt(points,nwet,2) ! Dew-point S_STATSP.158
& ,dr ! Randomly sampled dew pt. depression S_STATSP.159
! (K) S_STATSP.160
& ,f1 ! Used in calc of advection term S_STATSP.161
! in equation 2.12 S_STATSP.162
& ,n1, n2 ! Constants for interpolation S_STATSP.163
& ,qinc(points,nwet) ! Specific humidity increments S_STATSP.164
! (Kg Kg^-1 s^-1) S_STATSP.165
& ,qk(points,nwet) ! Factor to prevent Q becoming S_STATSP.166
! negative S_STATSP.167
& ,qrint(points,nwet) ! Specific humidity (interpolated S_STATSP.168
! values) (Kg Kg^-1) S_STATSP.169
& ,rdaysteps, rstepcount ! Real values of timesteps in day S_STATSP.170
! and timestep counter S_STATSP.171
& ,t0 ! Randomly sampled temp. for 1st S_STATSP.172
! level (K) S_STATSP.173
& ,tinc(points,nlevs) ! Temp. increments (K s^-1) S_STATSP.174
& ,trint(points,nlevs) ! Temp. (interpolated values) (K) S_STATSP.175
& ,vnrint(points,nlevs) ! Horizontal velocities (linearly S_STATSP.176
& ,vprint(points,nlevs) ! interpolated values) (m s^-1) S_STATSP.177
& ,wrint(points,ntrop) ! Vertical velocity (linearly S_STATSP.178
! interpolated values) S_STATSP.179
! (mb s^-1) S_STATSP.180
& ,G05DDE ! Function that samples randomly S_STATSP.181
! from a Gaussian distribution S_STATSP.182
! with a given mean and SD S_STATSP.183
S_STATSP.184
C S_STATSP.185
C S_STATSP.186
If (stepcount .eq. 1) then S_STATSP.187
If (daycount .eq. 1) then S_STATSP.188
i1 = 1 S_STATSP.189
j1 = 2 S_STATSP.190
else S_STATSP.191
i1 = 2 S_STATSP.192
j1 = 2 S_STATSP.193
C S_STATSP.194
C Save state of Random Number Generator for continuation S_STATSP.195
C STATS run done from tape to allow for the first day of a S_STATSP.196
C STATS run, when G05DDE is used twice as many times (to set S_STATSP.197
C up 2 profiles) and so the variables after forcing on a S_STATSP.198
C continuation run would be different from an unbroken run. S_STATSP.199
C - daycount ne 1. S_STATSP.200
C S_STATSP.201
Call G05CFE
(idum,iv,iy) S_STATSP.202
Do l = 1, points S_STATSP.203
Do k = 1, nlevs S_STATSP.204
tr(l,k,j1-1) = tr(l,k,j1) S_STATSP.205
vnr(l,k,j1-1) = vnr(l,k,j1) S_STATSP.206
vpr(l,k,j1-1) = vpr(l,k,j1) S_STATSP.207
enddo S_STATSP.208
Do k = 1, nwet S_STATSP.209
qr(l,k,j1-1) = qr(l,k,j1) S_STATSP.210
enddo S_STATSP.211
Do k = 1, ntrop S_STATSP.212
wr(l,k,j1-1) = wr(l,k,j1) S_STATSP.213
enddo S_STATSP.214
enddo S_STATSP.215
endif S_STATSP.216
C S_STATSP.217
C--------------------------------------------------------------------- S_STATSP.218
C Create new profiles (2 for 1st day and 1 from then on) S_STATSP.219
C G05DDE(A,B) returns a pseudo-random real number taken from S_STATSP.220
C a normal (Gaussian) distribution with mean A and SD B. S_STATSP.221
C--------------------------------------------------------------------- S_STATSP.222
C S_STATSP.223
Do l = 1, points S_STATSP.224
Do k = i1, j1 S_STATSP.225
t0 = G05DDE(tbar(l,1),tsd(l,1)) S_STATSP.226
tr(l,1,k) = t0+tdash(l,1) S_STATSP.227
d0 = G05DDE(dbar(l,1),dsd(l,1)) S_STATSP.228
dr = d0+ddash(l,1) S_STATSP.229
dewpt(l,1,k) = tr(l,1,k)-dr S_STATSP.230
vnr(l,1,k) = G05DDE(vnbar(l,1),vnsd(l,1)) S_STATSP.231
wr(l,1,k) = G05DDE(wbar(l,1),wsd(l,1)) S_STATSP.232
Do i = 1, nlevs-1 S_STATSP.233
ct = G05DDE(ctbar(l,i),ctsd(l,i)) S_STATSP.234
t0 = at(l,i)*t0+ct S_STATSP.235
tr(l,i+1,k) = t0+tdash(l,i+1) S_STATSP.236
cvn = G05DDE(cvnbar(l,i),cvnsd(l,i)) S_STATSP.237
vnr(l,i+1,k) = avn(l,i)*vnr(l,i,k)+cvn S_STATSP.238
enddo S_STATSP.239
Do i = 1, nwet-1 S_STATSP.240
cdd = G05DDE(cdbar(l,i),cdsd(l,i)) S_STATSP.241
d0 = ad(l,i)*d0+cdd S_STATSP.242
dr = d0+ddash(l,i+1) S_STATSP.243
dewpt(l,i+1,k) = tr(l,i+1,k)-dr S_STATSP.244
enddo S_STATSP.245
Do i = 1, ntrop-1 S_STATSP.246
cw = G05DDE(cwbar(l,i),cwsd(l,i)) S_STATSP.247
wr(l,i+1,k) = aw(l,i)*wr(l,i,k)+cw S_STATSP.248
enddo S_STATSP.249
Do i = 1, nlevs S_STATSP.250
vpr(l,i,k) = G05DDE(vpbar(l,i),vnsd(l,i)) S_STATSP.251
enddo S_STATSP.252
C S_STATSP.253
C After the first profile is set up on the first day, S_STATSP.254
C save state of Random Number Generator for continuation S_STATSP.255
C STATS run done from tape to allow for the first day of a S_STATSP.256
C STATS run, when G05DDF is used twice as many times (to set S_STATSP.257
C up 2 profiles) and so the variables after forcing on a S_STATSP.258
C continuation run would be different from an unbroken run. S_STATSP.259
C - daycount EQ 1. S_STATSP.260
C S_STATSP.261
If (k .eq. 1 .and. l .eq. 1) S_STATSP.262
& Call G05CFE(idum,iv,iy) S_STATSP.263
enddo S_STATSP.264
enddo S_STATSP.265
Do k = i1, j1 S_STATSP.266
Call QSAT(qr(1,1,k), dewpt(1,1,k), press, (points*nwet)) S_STATSP.267
enddo S_STATSP.268
endif ! stepcount=1 S_STATSP.269
C S_STATSP.270
C Interpolate between 2 values S_STATSP.271
C S_STATSP.272
rdaysteps = real(daysteps) S_STATSP.273
rstepcount = real(stepcount) S_STATSP.274
n1 = (rdaysteps-rstepcount+1.) / rdaysteps S_STATSP.275
n2 = (rstepcount-1.) / rdaysteps S_STATSP.276
Do k = 1, nlevs S_STATSP.277
Do l = 1, points S_STATSP.278
trint(l,k) = n1 * tr(l,k,1) + n2 * tr(l,k,2) S_STATSP.279
vnrint(l,k) = n1 * vnr(l,k,1) + n2 * vnr(l,k,2) S_STATSP.280
vprint(l,k) = n1 * vpr(l,k,1) + n2 * vpr(l,k,2) S_STATSP.281
enddo S_STATSP.282
enddo S_STATSP.283
Do k = 1, nwet S_STATSP.284
Do l = 1, points S_STATSP.285
qrint(l,k) = n1 * qr(l,k,1) + n2 * qr(l,k,2) S_STATSP.286
enddo S_STATSP.287
enddo S_STATSP.288
Do k = 1, ntrop S_STATSP.289
Do l = 1, points S_STATSP.290
wrint(l,k) = n1 * wr(l,k,1) + n2 * wr(l,k,2) S_STATSP.291
enddo S_STATSP.292
enddo S_STATSP.293
C S_STATSP.294
C Set U and V S_STATSP.295
C S_STATSP.296
Do k = 1, nlevs S_STATSP.297
Do l = 1, points S_STATSP.298
u(l,k) = vnrint(l,k) S_STATSP.299
v(l,k) = vprint(l,k) S_STATSP.300
enddo S_STATSP.301
enddo S_STATSP.302
C S_STATSP.303
C--------------------------------------------------------------------- S_STATSP.304
C Add vertical advection increments to T and Q (eqns. 2 and 3) S_STATSP.305
C--------------------------------------------------------------------- S_STATSP.306
C S_STATSP.307
Do k = 1, nlevs S_STATSP.308
Do l = 1, points S_STATSP.309
tinc(l,k) = 0.0 S_STATSP.310
enddo S_STATSP.311
enddo S_STATSP.312
Do k = 1, nwet S_STATSP.313
Do l = 1, points S_STATSP.314
qinc(l,k) = 0.0 S_STATSP.315
enddo S_STATSP.316
enddo S_STATSP.317
Do l = 1, points S_STATSP.318
dt1 = t(l,2) - t(l,1) S_STATSP.319
dq1 = q(l,2) - q(l,1) S_STATSP.320
Do i = 2, ntrop S_STATSP.321
dt2 = t(l,i+1) - t(l,i) S_STATSP.322
dq2 = q(l,i+1) - q(l,i) S_STATSP.323
f1 = -wrint(l,i) * timestep * rpress(l,i) S_STATSP.324
tinc(l,i) = f1 * (dt2 * px(l,i) + dt1 * px(l,i-1) S_STATSP.325
& - (dt1 + dt2) * py(l,i-1) - .2856 * t(l,i)) S_STATSP.326
qinc(l,i) = f1 * (dq2 * px(l,i) + dq1 * px(l,i-1) S_STATSP.327
& - (dq1+dq2) * py(l,i-1)) S_STATSP.328
t(l,i) = t(l,i) + tinc(l,i) S_STATSP.329
dt1 = dt2 S_STATSP.330
dq1 = dq2 S_STATSP.331
enddo S_STATSP.332
enddo ! l S_STATSP.333
C S_STATSP.334
C This section prevents the vertical increment from allowing S_STATSP.335
C Q to become negative (lowest value Q can take is 1.E-6) S_STATSP.336
C S_STATSP.337
Do k = 2, ntrop S_STATSP.338
Do l = 1, points S_STATSP.339
qk(l,k) = -q(l,k) + 1.e-6 S_STATSP.340
If (qinc(l,k) .lt. qk(l,k)) then S_STATSP.341
qinc(l,k) = qk(l,k) S_STATSP.342
endif S_STATSP.343
q(l,k) = q(l,k) + qinc(l,k) S_STATSP.344
enddo S_STATSP.345
enddo S_STATSP.346
C S_STATSP.347
C Printout increments due to vert. advection if required S_STATSP.348
C S_STATSP.349
If (prinstat) S_STATSP.350
& Call PRINTSUB( S_STATSP.351
C ! IN S_STATSP.352
& points, nlevs, nwet, S_STATSP.353
C ! S_STATSP.354
& ' Temperatur/moisture profiles + incs due to vert advection ' S_STATSP.355
& ,stepcount , dayno, t, q, tinc, qinc) S_STATSP.356
C S_STATSP.357
C--------------------------------------------------------------------- S_STATSP.358
C Add horizontal increments to T and Q (eqns. 2 and 3) S_STATSP.359
C--------------------------------------------------------------------- S_STATSP.360
C S_STATSP.361
Do k = 1, nlevs S_STATSP.362
Do l = 1, points S_STATSP.363
tinc(l,k) = timestep * abs(vnrint(l,k)) S_STATSP.364
& * (trint(l,k)-t(l,k)) S_STATSP.365
& / deltan(l) S_STATSP.366
t(l,k) = t(l,k) + tinc(l,k) S_STATSP.367
enddo S_STATSP.368
enddo S_STATSP.369
C S_STATSP.370
C This section prevents the horizontal increment from allowing S_STATSP.371
C Q to become negative (lowest value Q can take is 1.E-6) S_STATSP.372
C S_STATSP.373
Do k = 1, nwet S_STATSP.374
Do l = 1, points S_STATSP.375
qinc(l,k) = timestep * abs(vnrint(l,k)) S_STATSP.376
& * (qrint(l,k)-q(l,k)) S_STATSP.377
& / deltan(l) S_STATSP.378
qk(l,k) = -1 * q(l,k) + 1.e-6 S_STATSP.379
If (qinc(l,k) .lt. qk(l,k)) then S_STATSP.380
qinc(l,k) = qk(l,k) S_STATSP.381
endif S_STATSP.382
q(l,k) = q(l,k) + qinc(l,k) S_STATSP.383
enddo S_STATSP.384
enddo S_STATSP.385
C S_STATSP.386
C Printout increments due to horiz. advection if required S_STATSP.387
C S_STATSP.388
If (prinstat) S_STATSP.389
& Call PRINTSUB( S_STATSP.390
C ! IN S_STATSP.391
& points, nlevs, nwet, S_STATSP.392
C ! S_STATSP.393
& ' Temperature/moisture profiles + incs due to horiz advection' S_STATSP.394
& ,stepcount, dayno, t, q, tinc, qinc) S_STATSP.395
Return S_STATSP.396
End ! Subroutine STATSTEP S_STATSP.397
*ENDIF S_STATSP.398