*IF DEF,A03_3A,OR,DEF,A03_5A,OR,DEF,A03_5B,OR,DEF,A03_7A,OR,DEF,A03_6A ARN1F404.3
C ******************************COPYRIGHT****************************** GTS2F400.4447
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.4448
C GTS2F400.4449
C Use, duplication or disclosure of this code is subject to the GTS2F400.4450
C restrictions as set forth in the contract. GTS2F400.4451
C GTS2F400.4452
C Meteorological Office GTS2F400.4453
C London Road GTS2F400.4454
C BRACKNELL GTS2F400.4455
C Berkshire UK GTS2F400.4456
C RG12 2SZ GTS2F400.4457
C GTS2F400.4458
C If no contract has been raised with this copy of the code, the use, GTS2F400.4459
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.4460
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.4461
C Modelling at the above address. GTS2F400.4462
C ******************************COPYRIGHT****************************** GTS2F400.4463
C GTS2F400.4464
C*LL SUBROUTINE IMP_MIX ----------------------------------------------- IMPMIX2C.3
CLL IMPMIX2C.4
CLL Purpose: Calculate turbulent mixing increments for a passive tracer IMPMIX2C.5
CLL using an implicit numerical scheme. The tridiagonal IMPMIX2C.6
CLL matices are inverted using simple Gaussian elimination. IMPMIX2C.7
CLL IMPMIX2C.8
CLL IMPMIX2C.10
CLL Model Modification history : IMPMIX2C.11
CLL version Date IMPMIX2C.12
CLL 3.4 18/10/94 *DECK inserted into UM version 3.4. S Jackson IMPMIX2C.13
CLL 4.2 Oct. 96 T3E migration - *DEF CRAY removed GSS1F402.64
CLL S J Swarbrick GSS1F402.65
CLL 4.5 Jul. 98 Kill the IBM specific lines (JCThil) AJC1F405.143
CLL IMPMIX2C.14
CLL SDJ <- Programmers of some or all of previous code or changes IMPMIX2C.15
CLL IMPMIX2C.16
CLL Programming standard: UM Documentation Paper No 4, Version 2, IMPMIX2C.17
CLL dated 18/1/90 IMPMIX2C.18
CLL IMPMIX2C.19
CLL System component covered: P244 IMPMIX2C.20
CLL IMPMIX2C.21
CLL Project task: P24 IMPMIX2C.22
CLL IMPMIX2C.23
CLL Documentation: UM Documentation Paper No 24. IMPMIX2C.24
CLL IMPMIX2C.25
C*---------------------------------------------------------------------- IMPMIX2C.26
C*L Arguments :- IMPMIX2C.27
SUBROUTINE IMP_MIX ( 1,2IMPMIX2C.28
& P_FIELD,P1,P_POINTS,BL_LEVELS,DELTA_AK,DELTA_BK IMPMIX2C.29
&,GAMMA_RHOKH_RDZ,GAMMA_RHOK_DEP IMPMIX2C.30
&,PSTAR,TIMESTEP IMPMIX2C.31
&,F_FIELD,SURF_DEP_FLUX,FIELD IMPMIX2C.32
&,NRML IMPMIX2C.33
&,ERROR,LTIMER IMPMIX2C.34
& ) IMPMIX2C.35
IMPLICIT NONE IMPMIX2C.36
C IMPMIX2C.37
C Inputs :- IMPMIX2C.38
C IMPMIX2C.39
INTEGER IMPMIX2C.40
& P_FIELD ! IN No. of points in P-grid. IMPMIX2C.41
&,P1 ! IN First point to be processed in IMPMIX2C.42
C ! P-grid. IMPMIX2C.43
&,P_POINTS ! IN Number of P-grid points to be IMPMIX2C.44
C ! processed. IMPMIX2C.45
&,BL_LEVELS ! IN No. of atmospheric levels for IMPMIX2C.49
C ! which boundary layer fluxes are IMPMIX2C.50
C ! calculated. IMPMIX2C.51
REAL IMPMIX2C.52
& DELTA_AK(BL_LEVELS) ! IN Difference of hybrid 'A' across IMPMIX2C.53
C ! boundary layers (K-1/2 to K+1/2) IMPMIX2C.54
C ! (upper minus lower). IMPMIX2C.55
&,DELTA_BK(BL_LEVELS) ! IN Difference of hybrid 'B' across IMPMIX2C.56
C ! boundary layers (K-1/2 to K+1/2) IMPMIX2C.57
C ! (upper minus lower). IMPMIX2C.58
&,GAMMA_RHOKH_RDZ(P_FIELD,2:BL_LEVELS) IMPMIX2C.59
C ! IN Turbulent mixing coefs. above IMPMIX2C.60
C ! surface, =GAMMA(K)*RHOKH(,K) IMPMIX2C.61
C ! *RDZ(K) for K>=2 (from KMKH). IMPMIX2C.62
&,GAMMA_RHOK_DEP(P_FIELD) ! IN Surface exchange coefficient IMPMIX2C.63
C ! for surface deposition*GAMMA(1) IMPMIX2C.64
&,PSTAR(P_FIELD) ! IN Surface pressure (Pa). IMPMIX2C.65
&,TIMESTEP ! IN Timestep in seconds. IMPMIX2C.66
C IMPMIX2C.67
C Next 2 arrays are IN as "explicit" fluxes and OUT as "implicit" IMPMIX2C.68
C fluxes. IMPMIX2C.69
C IMPMIX2C.70
REAL IMPMIX2C.71
& F_FIELD(P_FIELD,BL_LEVELS) ! INOUT Flux of tracer IMPMIX2C.72
&,SURF_DEP_FLUX(P_FIELD) ! INOUT surface deposition flux IMPMIX2C.73
&,FIELD(P_FIELD,BL_LEVELS) ! INOUT Amount of tracer IMPMIX2C.74
IMPMIX2C.75
INTEGER IMPMIX2C.76
& NRML(P_FIELD) ! IN The number of model layers IMPMIX2C.77
C ! in the unstable rapidly mixing IMPMIX2C.78
C ! layer. Zero if surface layer IMPMIX2C.79
C ! is stable. IMPMIX2C.80
&,ERROR ! OUT 1 if bad arguments, else 0. IMPMIX2C.81
IMPMIX2C.82
LOGICAL IMPMIX2C.83
& LTIMER ! IN Logical switch for TIMER diagnos IMPMIX2C.84
IMPMIX2C.85
C* IMPMIX2C.86
C*L External references :- IMPMIX2C.87
IMPMIX2C.88
EXTERNAL TIMER IMPMIX2C.89
IMPMIX2C.90
C* IMPMIX2C.91
C*L Local and other symbolic constants :- IMPMIX2C.92
*CALL C_G 
IMPMIX2C.93
C* IMPMIX2C.94
C*L Workspace :- IMPMIX2C.95
C 4*BL_LEVELS + 4 blocks of real workspace are required. IMPMIX2C.96
REAL IMPMIX2C.98
& AF(P_FIELD,BL_LEVELS) ! Elements in rows in matrix IMPMIX2C.99
C ! equation (modified during IMPMIX2C.100
C ! Gaussian elimination calculations). IMPMIX2C.101
&,AF_RML(P_FIELD) ! Matrix element for field in IMPMIX2C.102
C ! rapidly mixing layer. IMPMIX2C.103
&,DELTAP(P_FIELD,BL_LEVELS) ! Vertical pressure difference across IMPMIX2C.104
C ! hybrid layers (upper minus lower) IMPMIX2C.105
C ! (Pa). IMPMIX2C.106
&,DELTAP_RML(P_FIELD) ! Vertical pressure difference across IMPMIX2C.107
C ! the rapidly mixing layer (Pa). IMPMIX2C.108
&,D_FIELD(P_FIELD,BL_LEVELS) ! Delta FIELD (tracer field) IMPMIX2C.109
C ! elements of vector on RHS, then IMPMIX2C.110
C ! LHS, of eqn P244.79. IMPMIX2C.111
&,DFLD_RML(P_FIELD) ! Delta FIELD for rapidly IMPMIX2C.112
C ! mixing layer. IMPMIX2C.113
&,DTRDZ(P_FIELD,BL_LEVELS) ! -g.dt/dp for bottom BL_LEVELS IMPMIX2C.114
C ! model layers (needed in P245). IMPMIX2C.115
&,DTRDZ_RML(P_FIELD) ! -g.dt/dp for the rapidly IMPMIX2C.116
C ! mixing layer (needed in P245). IMPMIX2C.117
C* IMPMIX2C.141
C Local scalars :- IMPMIX2C.142
REAL IMPMIX2C.143
& CF ! Matrix element for local increments to rml IMPMIX2C.144
&,CF_RML ! As above but for rapidly mixing layer increment. IMPMIX2C.145
&,RBF ! Reciprocal of B for local increments to rml IMPMIX2C.146
&,RBF_RML ! As above but for the rapidly mixing layer increment. IMPMIX2C.147
&,DTIMEG ! TIMESTEP * G (used in several places). IMPMIX2C.148
INTEGER IMPMIX2C.149
& BLM1 ! BL_LEVELS minus 1. IMPMIX2C.150
&,NRMLP1 ! NRML plus 1. IMPMIX2C.151
&,I ! Loop counter (horizontal field index). IMPMIX2C.152
&,J ! Offset version of I. IMPMIX2C.153
&,K ! Loop counter (vertical index). IMPMIX2C.154
&,KM1 ! K minus 1. IMPMIX2C.155
&,KP1 ! K plus 1. IMPMIX2C.156
C IMPMIX2C.157
IF (LTIMER) THEN IMPMIX2C.158
CALL TIMER('IMPMIX ',3) IMPMIX2C.159
ENDIF IMPMIX2C.160
IMPMIX2C.161
ERROR=0 IMPMIX2C.162
IMPMIX2C.163
DTIMEG = TIMESTEP * G IMPMIX2C.164
BLM1 = BL_LEVELS-1 IMPMIX2C.165
IMPMIX2C.166
C IMPMIX2C.167
CL---------------------------------------------------------------------- IMPMIX2C.168
CL (A) Calculations on P-grid. IMPMIX2C.169
CL---------------------------------------------------------------------- IMPMIX2C.170
C IMPMIX2C.171
C----------------------------------------------------------------------- IMPMIX2C.172
CL 1. Calculate pressure across layer (in hybrid coordinates), DELTAP, IMPMIX2C.173
CL and then -gdt/dP = dt/rho*dz for use throughout section (A) IMPMIX2C.174
C----------------------------------------------------------------------- IMPMIX2C.175
C IMPMIX2C.176
DO 1 K=1,BL_LEVELS IMPMIX2C.177
DO 11 I=P1,P1+P_POINTS-1 IMPMIX2C.178
DELTAP(I,K) = DELTA_AK(K) + PSTAR(I)*DELTA_BK(K) IMPMIX2C.179
DTRDZ(I,K) = -DTIMEG / DELTAP(I,K) IMPMIX2C.180
11 CONTINUE IMPMIX2C.181
1 CONTINUE IMPMIX2C.182
C IMPMIX2C.183
C----------------------------------------------------------------------- IMPMIX2C.184
CL 2. Calculate implicit FIELD increments due to local mixing within IMPMIX2C.185
CL the rapidly mixing layer (where it exists). IMPMIX2C.186
CL The surface fluxes F_FIELD(I,1) and SURF_DEP_FLUX(I) are used for IMPMIX2C.187
CL calculating the rapidly mixing layer increments so not here. IMPMIX2C.188
CL Therefore the matrix equation we must solve to find the implicit IMPMIX2C.189
CL FIELD increments due to local mixing within the rml does not IMPMIX2C.190
CL have a "surface" row. IMPMIX2C.191
C----------------------------------------------------------------------- IMPMIX2C.192
C IMPMIX2C.193
C----------------------------------------------------------------------- IMPMIX2C.194
CL 2.1 Start 'upward sweep' with lowest model layer, which will be the IMPMIX2C.195
CL bottom of the rapidly mixing layer (rml) if it exists. IMPMIX2C.196
C----------------------------------------------------------------------- IMPMIX2C.197
C IMPMIX2C.198
DO 21 I=P1,P1+P_POINTS-1 IMPMIX2C.199
IF (NRML(I) .GE. 2) THEN IMPMIX2C.200
C IMPMIX2C.201
C "Explicit" increments due to local mixing within the rml. IMPMIX2C.202
C IMPMIX2C.203
D_FIELD(I,1) = -DTRDZ(I,1) * F_FIELD(I,2) IMPMIX2C.204
C IMPMIX2C.205
C Define matrix elements (CF always zero for this case). IMPMIX2C.206
C IMPMIX2C.207
AF(I,1) = -DTRDZ(I,1) * GAMMA_RHOKH_RDZ(I,2) IMPMIX2C.208
RBF = 1.0 / ( 1.0 - AF(I,1) ) IMPMIX2C.209
C IMPMIX2C.210
C Now start Gaussian elimination IMPMIX2C.211
C IMPMIX2C.212
D_FIELD(I,1) = RBF * D_FIELD(I,1) IMPMIX2C.213
AF(I,1) = RBF * AF(I,1) IMPMIX2C.214
C IMPMIX2C.215
C Start calculating DELTAP_RML. Mid-level depths added in 2.2 below. IMPMIX2C.216
C IMPMIX2C.217
DELTAP_RML(I) = DELTAP(I,1) IMPMIX2C.218
ELSE ! No rapidly mixing layer calculations IMPMIX2C.219
DTRDZ_RML(I) = 1.E30 IMPMIX2C.220
DFLD_RML(I) = 1.E30 IMPMIX2C.221
AF_RML(I) = 1.E30 IMPMIX2C.222
DELTAP_RML(I) = 1.E30 IMPMIX2C.223
ENDIF IMPMIX2C.224
21 CONTINUE IMPMIX2C.225
C IMPMIX2C.226
C----------------------------------------------------------------------- IMPMIX2C.227
CL 2.2 Continue upward sweep through middle of the rapidly mixing layer IMPMIX2C.228
CL (if it exists) and to its top. NB NRML is always < or = BLM1. IMPMIX2C.229
C----------------------------------------------------------------------- IMPMIX2C.230
C IMPMIX2C.231
DO 22 K=2,BLM1 IMPMIX2C.232
KP1 = K+1 IMPMIX2C.233
KM1 = K-1 IMPMIX2C.234
DO 221 I=P1,P1+P_POINTS-1 IMPMIX2C.235
C IMPMIX2C.236
C If in the top rapidly mixing layer then do not include flux at its IMPMIX2C.237
C top in the calculation, i.e. F_FIELD(I,NRML+1) is not IMPMIX2C.238
C included here; it is taken account of in the non-local mixing IMPMIX2C.239
C through the "rapidly mixing layer". IMPMIX2C.240
C IMPMIX2C.241
IF ( K .EQ. NRML(I) ) THEN IMPMIX2C.242
C IMPMIX2C.243
C Add final DELTAP contribution to DELTAP_RML and then calculate IMPMIX2C.244
C DTRDZ_RML. Lower level contributions added in 2.1 and below. IMPMIX2C.245
C IMPMIX2C.246
DELTAP_RML(I) = DELTAP_RML(I) + DELTAP(I,K) IMPMIX2C.247
DTRDZ_RML(I) = -DTIMEG / DELTAP_RML(I) IMPMIX2C.248
C IMPMIX2C.249
C "Explicit" flux divergence across layer giving explicit IMPMIX2C.250
C increment due to the local mixing at the top of rml. IMPMIX2C.251
C IMPMIX2C.252
D_FIELD(I,K) = DTRDZ(I,K) * F_FIELD(I,K) IMPMIX2C.253
C IMPMIX2C.254
C Define matrix elements (AF always zero for this case). IMPMIX2C.255
C IMPMIX2C.256
CF = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K) IMPMIX2C.257
RBF = 1.0 / ( 1.0 - CF*( 1.0 + AF(I,KM1) ) ) IMPMIX2C.258
C IMPMIX2C.259
C Now start Gaussian elimination IMPMIX2C.260
C IMPMIX2C.261
D_FIELD(I,K) = RBF * ( D_FIELD(I,K) IMPMIX2C.262
& - CF*D_FIELD(I,KM1) ) IMPMIX2C.263
ELSEIF (K .LT. NRML(I)) THEN IMPMIX2C.264
C IMPMIX2C.265
C Add layer depths to form total rml depth. IMPMIX2C.266
C IMPMIX2C.267
DELTAP_RML(I) = DELTAP_RML(I) + DELTAP(I,K) IMPMIX2C.268
C IMPMIX2C.269
C "Explicit" flux divergence across layer giving explicit IMPMIX2C.270
C increment due to the local mixing. IMPMIX2C.271
C IMPMIX2C.272
D_FIELD(I,K) = -DTRDZ(I,K) * (F_FIELD(I,KP1) - F_FIELD(I,K)) IMPMIX2C.273
C IMPMIX2C.274
C Define matrix elements. IMPMIX2C.275
C IMPMIX2C.276
AF(I,K) = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,KP1) IMPMIX2C.277
CF = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K) IMPMIX2C.278
RBF = 1.0 / ( 1.0 - AF(I,K) - CF * ( 1.0 + AF(I,KM1) ) ) IMPMIX2C.279
C IMPMIX2C.280
C Now start Gaussian elimination IMPMIX2C.281
C IMPMIX2C.282
D_FIELD(I,K) = RBF * ( D_FIELD(I,K) - CF*D_FIELD(I,KM1) ) IMPMIX2C.283
AF(I,K) = RBF * AF(I,K) IMPMIX2C.284
ENDIF IMPMIX2C.285
221 CONTINUE IMPMIX2C.286
22 CONTINUE IMPMIX2C.287
C IMPMIX2C.288
C----------------------------------------------------------------------- IMPMIX2C.289
CL 2.3 Downward sweep through matrix. Add implicit increments due to IMPMIX2C.290
CL local mixing within the rapidly mixing layer. Update surface IMPMIX2C.291
CL deposition flux and the top-of-rml tracer flux using IMPMIX2C.292
CL local mixing increments for layers 1 and NRML respectively. IMPMIX2C.293
C----------------------------------------------------------------------- IMPMIX2C.294
C IMPMIX2C.295
DO 23 K=BLM1,1,-1 IMPMIX2C.296
KP1 = K + 1 IMPMIX2C.297
DO 231 I=P1,P1+P_POINTS-1 IMPMIX2C.298
IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRML(I))) THEN IMPMIX2C.299
FIELD(I,K) = FIELD(I,K) + D_FIELD(I,K) IMPMIX2C.300
F_FIELD(I,KP1) = F_FIELD(I,KP1) IMPMIX2C.301
& + GAMMA_RHOKH_RDZ(I,KP1)*D_FIELD(I,K) IMPMIX2C.302
ELSEIF ((NRML(I) .GE. 2) .AND. (K .LT. NRML(I))) THEN IMPMIX2C.303
D_FIELD(I,K) = D_FIELD(I,K) IMPMIX2C.304
& - AF(I,K)*D_FIELD(I,KP1) IMPMIX2C.305
FIELD(I,K) = FIELD(I,K) + D_FIELD(I,K) IMPMIX2C.306
ENDIF IMPMIX2C.307
IF ((NRML(I) .GE. 2) .AND. (K .EQ. 1)) THEN IMPMIX2C.308
SURF_DEP_FLUX(I) = SURF_DEP_FLUX(I) IMPMIX2C.309
& - GAMMA_RHOK_DEP(I) * D_FIELD(I,1) IMPMIX2C.310
ENDIF IMPMIX2C.311
231 CONTINUE IMPMIX2C.312
23 CONTINUE IMPMIX2C.313
C IMPMIX2C.314
C----------------------------------------------------------------------- IMPMIX2C.315
CL 4. Calculate those matrix and vector elements on the LHS of eqn IMPMIX2C.316
CL which are to do with implicit solution of the tracer IMPMIX2C.317
CL transport problem at the surface, above the rml (if it exists) IMPMIX2C.318
CL and between all levels if it does not. IMPMIX2C.319
CL Begin with "upward sweep" through lower half of matrix). IMPMIX2C.320
C----------------------------------------------------------------------- IMPMIX2C.321
C IMPMIX2C.322
DO 31 I=P1,P1+P_POINTS-1 IMPMIX2C.323
C----------------------------------------------------------------------- IMPMIX2C.324
CL 4.2 Lowest atmospheric layer FIELD row of matrix. IMPMIX2C.325
C----------------------------------------------------------------------- IMPMIX2C.326
IF (NRML(I) .GE. 2) THEN IMPMIX2C.327
NRMLP1 = NRML(I) + 1 IMPMIX2C.328
C IMPMIX2C.329
C "Explicit" rapidly mixing layer increment for FIELD. IMPMIX2C.330
C IMPMIX2C.331
DFLD_RML(I) = -DTRDZ_RML(I) * IMPMIX2C.332
& (F_FIELD(I,NRMLP1) - F_FIELD(I,1) IMPMIX2C.333
& - SURF_DEP_FLUX(I)) IMPMIX2C.334
AF_RML(I) = -DTRDZ_RML(I) * GAMMA_RHOKH_RDZ(I,NRMLP1) IMPMIX2C.335
CF_RML = -DTRDZ_RML(I) * GAMMA_RHOK_DEP(I) IMPMIX2C.336
RBF_RML = 1.0 / ( 1.0 - AF_RML(I) - CF_RML ) IMPMIX2C.337
DFLD_RML(I) = RBF_RML * DFLD_RML(I) IMPMIX2C.338
AF_RML(I) = RBF_RML * AF_RML(I) IMPMIX2C.339
ELSE IMPMIX2C.340
C IMPMIX2C.341
C "Explicit" increment to FIELD(1) when there is no rapidly mixing layer IMPMIX2C.342
C or it does not extend beyond the bottom model layer. IMPMIX2C.343
C IMPMIX2C.344
D_FIELD(I,1) = -DTRDZ(I,1) * IMPMIX2C.345
& ( F_FIELD(I,2) - F_FIELD(I,1) IMPMIX2C.346
& - SURF_DEP_FLUX(I) ) IMPMIX2C.347
IMPMIX2C.348
CF = -DTRDZ(I,1) * GAMMA_RHOK_DEP(I) IMPMIX2C.349
AF(I,1) = -DTRDZ(I,1) * GAMMA_RHOKH_RDZ(I,2) IMPMIX2C.350
RBF = 1.0 / ( 1.0 - AF(I,1) - CF ) IMPMIX2C.351
D_FIELD(I,1) = RBF * D_FIELD(I,1) IMPMIX2C.352
AF(I,1) = RBF * AF(I,1) IMPMIX2C.353
ENDIF IMPMIX2C.354
31 CONTINUE IMPMIX2C.355
C----------------------------------------------------------------------- IMPMIX2C.356
CL 4.3 Rows of matrix applying to FIELD transport into model layers in IMPMIX2C.357
CL the "middle" of the "boundary" layer, i.e. all but the bottom IMPMIX2C.358
CL layer and the top "boundary" layer. IMPMIX2C.359
C----------------------------------------------------------------------- IMPMIX2C.360
DO 43 K=2,BLM1 IMPMIX2C.361
KP1 = K+1 IMPMIX2C.362
KM1 = K-1 IMPMIX2C.363
DO 431 I=P1,P1+P_POINTS-1 IMPMIX2C.364
C IMPMIX2C.365
C "Explicit" flux divergence across layer giving explicit FIELD IMPMIX2C.366
C increment due to mixing above rml if it exists or for all levels if IMPMIX2C.367
C it does not. IMPMIX2C.368
C IMPMIX2C.369
NRMLP1 = NRML(I) + 1 IMPMIX2C.370
IF (K .GT. NRML(I)) THEN IMPMIX2C.371
D_FIELD(I,K) = -DTRDZ(I,K) * (F_FIELD(I,KP1) - F_FIELD(I,K)) IMPMIX2C.372
AF(I,K) = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,KP1) IMPMIX2C.373
CF = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K) IMPMIX2C.374
IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRMLP1)) THEN IMPMIX2C.375
RBF = 1.0 / ( 1.0 - AF(I,K) IMPMIX2C.376
& - CF * ( 1.0 + AF_RML(I) ) ) IMPMIX2C.377
D_FIELD(I,K) = RBF * ( D_FIELD(I,K) IMPMIX2C.378
& - CF*DFLD_RML(I) ) IMPMIX2C.379
ELSE IMPMIX2C.380
RBF = 1.0 / ( 1.0 - AF(I,K) IMPMIX2C.381
& - CF * ( 1.0 + AF(I,KM1) ) ) IMPMIX2C.382
D_FIELD(I,K) = RBF * ( D_FIELD(I,K) IMPMIX2C.383
& - CF*D_FIELD(I,KM1) ) IMPMIX2C.384
ENDIF IMPMIX2C.385
AF(I,K) = RBF * AF(I,K) IMPMIX2C.386
ENDIF IMPMIX2C.387
431 CONTINUE IMPMIX2C.388
43 CONTINUE IMPMIX2C.389
C----------------------------------------------------------------------- IMPMIX2C.390
CL 4.4 Top "boundary" layer FIELD row of matrix. FIELD for this layer IMPMIX2C.391
CL can then be, and is, updated. IMPMIX2C.392
C----------------------------------------------------------------------- IMPMIX2C.393
DO 44 I=P1,P1+P_POINTS-1 IMPMIX2C.394
D_FIELD(I,BL_LEVELS) = DTRDZ(I,BL_LEVELS) * F_FIELD(I,BL_LEVELS) IMPMIX2C.395
C IMPMIX2C.396
CF = -DTRDZ(I,BL_LEVELS) * GAMMA_RHOKH_RDZ(I,BL_LEVELS) IMPMIX2C.397
IF (NRML(I) .EQ. BLM1) THEN IMPMIX2C.398
RBF = 1.0 / ( 1.0 - CF*( 1.0 + AF_RML(I) ) ) IMPMIX2C.399
D_FIELD(I,BL_LEVELS) = RBF * ( D_FIELD(I,BL_LEVELS) IMPMIX2C.400
& - CF*DFLD_RML(I) ) IMPMIX2C.401
ELSE IMPMIX2C.402
RBF = 1.0 / ( 1.0 - CF*( 1.0 + AF(I,BLM1) ) ) IMPMIX2C.403
D_FIELD(I,BL_LEVELS) = RBF * ( D_FIELD(I,BL_LEVELS) IMPMIX2C.404
& - CF*D_FIELD(I,BLM1) ) IMPMIX2C.405
ENDIF IMPMIX2C.406
FIELD(I,BL_LEVELS) = FIELD(I,BL_LEVELS) + D_FIELD(I,BL_LEVELS) IMPMIX2C.407
44 CONTINUE IMPMIX2C.408
C IMPMIX2C.409
C----------------------------------------------------------------------- IMPMIX2C.410
CL 5. "Downward sweep" through whole matrix. FIELD is updated when IMPMIX2C.411
CL the final implicit increments have been calculated. IMPMIX2C.412
C----------------------------------------------------------------------- IMPMIX2C.413
CL 5.1 Remaining FIELD rows of matrix and add implicit increments above IMPMIX2C.414
CL the rml or at all levels if it is less than two layers thick. IMPMIX2C.415
C----------------------------------------------------------------------- IMPMIX2C.416
C IMPMIX2C.417
DO 51 K=BLM1,1,-1 IMPMIX2C.418
DO 511 I=P1,P1+P_POINTS-1 IMPMIX2C.419
IF ( (K .GT. NRML(I)) .OR. (NRML(I) .LT. 2) ) THEN IMPMIX2C.420
D_FIELD(I,K) = D_FIELD(I,K) - AF(I,K)*D_FIELD(I,K+1) IMPMIX2C.421
FIELD(I,K) = FIELD(I,K) + D_FIELD(I,K) IMPMIX2C.422
ENDIF IMPMIX2C.423
511 CONTINUE IMPMIX2C.424
51 CONTINUE IMPMIX2C.425
C IMPMIX2C.426
C----------------------------------------------------------------------- IMPMIX2C.427
CL 5.2 Rapidly mixing layer increments IMPMIX2C.428
C----------------------------------------------------------------------- IMPMIX2C.429
DO 52 I=P1,P1+P_POINTS-1 IMPMIX2C.430
IF ( NRML(I) .GE. 2 ) THEN IMPMIX2C.431
NRMLP1 = NRML(I) + 1 IMPMIX2C.432
DFLD_RML(I) = DFLD_RML(I) - AF_RML(I) * D_FIELD(I,NRMLP1) IMPMIX2C.433
FIELD(I,1) = FIELD(I,1) + DFLD_RML(I) IMPMIX2C.434
ENDIF IMPMIX2C.435
52 CONTINUE IMPMIX2C.436
C IMPMIX2C.437
C----------------------------------------------------------------------- IMPMIX2C.438
CL 5.4 Add implicit increments due to rapid mixing if in rapid mixing IMPMIX2C.439
CL layer. IMPMIX2C.440
C----------------------------------------------------------------------- IMPMIX2C.441
DO 54 K=2,BL_LEVELS IMPMIX2C.442
DO 541 I=P1,P1+P_POINTS-1 IMPMIX2C.443
IF (K .LE. NRML(I)) THEN IMPMIX2C.444
C IMPMIX2C.445
C Add the increments due to rapid mixing if in the rapidly mixing layer IMPMIX2C.446
C IMPMIX2C.447
FIELD(I,K) = FIELD(I,K) + DFLD_RML(I) IMPMIX2C.448
ENDIF IMPMIX2C.449
541 CONTINUE IMPMIX2C.450
54 CONTINUE IMPMIX2C.451
C IMPMIX2C.452
C----------------------------------------------------------------------- IMPMIX2C.453
CL 6. Calculate final implicit flux of tracer. IMPMIX2C.454
C----------------------------------------------------------------------- IMPMIX2C.455
CL 6.1 Surface fluxes. IMPMIX2C.456
C----------------------------------------------------------------------- IMPMIX2C.457
C IMPMIX2C.458
DO 61 I=P1,P1+P_POINTS-1 IMPMIX2C.459
IF ( NRML(I) .GE. 2 ) THEN IMPMIX2C.460
SURF_DEP_FLUX(I) = SURF_DEP_FLUX(I) - GAMMA_RHOK_DEP(I) * IMPMIX2C.461
& DFLD_RML(I) IMPMIX2C.462
ELSE IMPMIX2C.463
SURF_DEP_FLUX(I) = SURF_DEP_FLUX(I) - GAMMA_RHOK_DEP(I) * IMPMIX2C.464
& D_FIELD(I,1) IMPMIX2C.465
ENDIF IMPMIX2C.466
61 CONTINUE IMPMIX2C.467
C IMPMIX2C.468
C----------------------------------------------------------------------- IMPMIX2C.469
CL 6.2 Fluxes at layer interfaces above the surface. IMPMIX2C.470
C----------------------------------------------------------------------- IMPMIX2C.471
DO 62 K=2,BL_LEVELS IMPMIX2C.472
KM1 = K-1 IMPMIX2C.473
DO 621 I=P1,P1+P_POINTS-1 IMPMIX2C.474
C IMPMIX2C.475
C Calculate and store fluxes due to local mixing. IMPMIX2C.476
C F_FIELD(local mixing) stored in array AF. IMPMIX2C.477
C IMPMIX2C.478
NRMLP1 = NRML(I) + 1 IMPMIX2C.479
IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRMLP1)) THEN IMPMIX2C.480
AF(I,K) = F_FIELD(I,K) - GAMMA_RHOKH_RDZ(I,K) IMPMIX2C.481
& * ( D_FIELD(I,K) - DFLD_RML(I) ) IMPMIX2C.482
ELSE IMPMIX2C.483
AF(I,K) = F_FIELD(I,K) - GAMMA_RHOKH_RDZ(I,K) IMPMIX2C.484
& * ( D_FIELD(I,K) - D_FIELD(I,KM1) ) IMPMIX2C.485
ENDIF IMPMIX2C.486
C IMPMIX2C.487
C Now calculate the implicit fluxes including both local mixing and IMPMIX2C.488
C if appropriate also the fluxes due to rapid mixing through layers. IMPMIX2C.489
C IMPMIX2C.490
IF ( K .EQ. 2 ) THEN IMPMIX2C.491
IF ( NRML(I) .GE. 2 ) THEN IMPMIX2C.492
F_FIELD(I,K) = AF(I,K) IMPMIX2C.493
& + F_FIELD(I,KM1) + SURF_DEP_FLUX(I) IMPMIX2C.494
& - DFLD_RML(I) / DTRDZ(I,KM1) IMPMIX2C.495
ELSE IMPMIX2C.496
F_FIELD(I,K) = AF(I,K) IMPMIX2C.497
ENDIF IMPMIX2C.498
ELSEIF ( K .LE. NRML(I) ) THEN IMPMIX2C.499
F_FIELD(I,K) = AF(I,K) - AF(I,KM1) IMPMIX2C.500
& + F_FIELD(I,KM1) - DFLD_RML(I) / DTRDZ(I,KM1) IMPMIX2C.501
ELSE IMPMIX2C.502
F_FIELD(I,K) = AF(I,K) IMPMIX2C.503
ENDIF IMPMIX2C.504
621 CONTINUE IMPMIX2C.505
62 CONTINUE IMPMIX2C.506
C IMPMIX2C.507
999 CONTINUE ! Branch for error exit. IMPMIX2C.508
IMPMIX2C.509
IF (LTIMER) THEN IMPMIX2C.510
CALL TIMER('IMPMIX ',4) IMPMIX2C.511
ENDIF IMPMIX2C.512
IMPMIX2C.513
RETURN IMPMIX2C.514
END IMPMIX2C.515
*ENDIF IMPMIX2C.516