*IF DEF,OCEAN @DYALLOC.4677
C ******************************COPYRIGHT****************************** GTS2F400.11521
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.11522
C GTS2F400.11523
C Use, duplication or disclosure of this code is subject to the GTS2F400.11524
C restrictions as set forth in the contract. GTS2F400.11525
C GTS2F400.11526
C Meteorological Office GTS2F400.11527
C London Road GTS2F400.11528
C BRACKNELL GTS2F400.11529
C Berkshire UK GTS2F400.11530
C RG12 2SZ GTS2F400.11531
C GTS2F400.11532
C If no contract has been raised with this copy of the code, the use, GTS2F400.11533
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.11534
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.11535
C Modelling at the above address. GTS2F400.11536
C ******************************COPYRIGHT****************************** GTS2F400.11537
C GTS2F400.11538
C*LL Subroutine VDIFCALT VDIFCALT.3
CLL VDIFCALT.4
CLL Can run on any compiler accepting lower case variables. VDIFCALT.5
CLL VDIFCALT.6
CLL The code must be precompiled by the UPDOC system. VDIFCALT.7
CLL Option L selects the code for Isopycnal diffusion. VDIFCALT.8
CLL Option V selects code necessary when I.P.D. is combined with VDIFCALT.9
CLL variable time-step with depth. VDIFCALT.10
CLL VDIFCALT.11
CLL Author: D.J.Carrington VDIFCALT.12
CLL Date: 12 December 1990 VDIFCALT.13
CLL Reviewer: R.A.Wood VDIFCALT.14
CLL Date: 19 December 1990 VDIFCALT.15
CLL Version 1.0 VDIFCALT.16
CLL VDIFCALT.17
! ORH1F305.5382
! Modification History: ORH1F305.5383
! Version Date Details ORH1F305.5384
! ------- ------- ------------------------------------------ ORH1F305.5385
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.5386
CLL 4.5 G.J.Rickard Introduce counter gradient flux terms for OOM1F405.864
CLL use with the full Large scheme OOM1F405.865
CLL OOM1F405.866
CLL External documentation: VDIFCALT.18
CLL Unified Model Documentation Paper No. 51. VDIFCALT.19
CLL VDIFCALT.20
CLL Naming convention of variables: Cox naming convention is used, VDIFCALT.21
CLL with the addition that lower-case variables are VDIFCALT.22
CLL introduced by the Isopycnal Diffusion scheme. VDIFCALT.23
CLL VDIFCALT.24
CLL Purpose of Subroutine: VDIFCALT.25
CLL Solves the vertical diffusion equation for tracers. VDIFCALT.26
CLL This separate treatment of the vertical diffusion allows a VDIFCALT.27
CLL greater time-step than would otherwise be possible. VDIFCALT.28
CLL Note that surface tracer fluxes are NOT introduced in this VDIFCALT.29
CLL routine; the routine is called after SFCADD has added them in. VDIFCALT.30
CLL VDIFCALT.31
CLL List of subroutines required for implementation of Isopycnal VDIFCALT.32
CLL Diffusion Scheme (in order of being called): VDIFCALT.33
CLL VERTCOFC * VDIFCALT.34
CLL VDIFCALC VDIFCALT.35
CLL VERTCOFT * VDIFCALT.36
CLL IPDCOFCL VDIFCALT.37
CLL IPDFLXCL VDIFCALT.38
CLL VDIFCALT * K-theory mixing scheme VDIFCALT.39
CLL VDIFCALT.40
CLLEND------------------------------------------------------------------ VDIFCALT.41
C* VDIFCALT.42
C*L---- Arguments ------------------------------------------------------ VDIFCALT.43
C VDIFCALT.44
SUBROUTINE VDIFCALT 2VDIFCALT.45
& ( J,IMT,IMTM1,KM,KMP1,KMM1,NT, VDIFCALT.46
& TA, VDIFCALT.47
& DZ2R,DZZ2RQ,DZ2RQ,C2DTTS, VDIFCALT.48
& DTTSA, ORH1F305.5389
& FM, VDIFCALT.54
& gnu,CGFLUX OOM1F405.867
& ) VDIFCALT.56
C VDIFCALT.57
C VDIFCALT.58
IMPLICIT NONE VDIFCALT.59
*CALL CNTLOCN 
ORH1F305.5387
*CALL OARRYSIZ ORH1F305.5388
C VDIFCALT.60
INTEGER VDIFCALT.61
& I,J,K,M,IMT,IMTM1,KM,KMP1,KMM1,NT VDIFCALT.62
C VDIFCALT.63
REAL VDIFCALT.64
& TA(IMT,KM,NT), VDIFCALT.65
& DZ2R(KM),DZ2RQ(IMT,KM),DZZ2RQ(IMT,KM), VDIFCALT.66
& FM(IMT,KM), VDIFCALT.67
& C2DTTS, VDIFCALT.68
& DTTSA(KM) , !Dimensioned when variable timestep used ORH1F305.5390
& gnu(IMT,KM) ! IN Vert diffusivity at TOP of box (cm2/s) VDIFCALT.74
&, CGFLUX(IMT,KM,NT) !IN counter gradient flux at TOP of box OOM1F405.868
C VDIFCALT.75
C VDIFCALT.76
C Declare local variables and arrays VDIFCALT.77
C VDIFCALT.78
REAL VDIFCALT.79
& aa(IMT,KM), ! c(1)A(k) } Eq.1.20 VDIFCALT.80
& cc(IMT,KM), ! c(2)A(k+1) } VDIFCALT.81
& bb(IMT,KM), ! aa+cc+1 VDIFCALT.82
& ee(IMT,KMP1), ! x(k-1) } Eq.1.18 VDIFCALT.83
& ff(IMT,KMP1,NT), ! y(k-1) } VDIFCALT.84
& efdr(IMT), ! 1/(bb-cc*ee) VDIFCALT.85
& tempa(IMT,KMP1), ! workspace VDIFCALT.86
& tempb(IMT,KMP1), ! workspace VDIFCALT.87
& CGTERM(IMT,KM,NT), !cgflux workspace on Tracer levels OOM1F405.869
& fxa, ! } VDIFCALT.88
& fxb, ! } local csts. VDIFCALT.89
& fxc, ! } VDIFCALT.90
& fxd ! } VDIFCALT.91
C* VDIFCALT.92
C VDIFCALT.93
C -------------------------------------------------------------- VDIFCALT.94
CL The arrays are set up here, VDIFCALT.95
CL in preparation for the recurrence method that follows. VDIFCALT.96
C -------------------------------------------------------------- VDIFCALT.97
C VDIFCALT.98
fxa=4. VDIFCALT.99
fxb=1. VDIFCALT.100
fxc=0. VDIFCALT.101
fxd=2.0 OOM1F405.870
C VDIFCALT.102
CL Certain quantities for the surface & bottom are set to zero VDIFCALT.103
CL to provide the boundary conditions VDIFCALT.104
C VDIFCALT.105
DO 210 I=1,IMT VDIFCALT.106
aa(I, 1)=fxc VDIFCALT.107
cc(I, KM)=fxc VDIFCALT.108
ee(I,KM+1)=fxc VDIFCALT.109
210 CONTINUE VDIFCALT.110
DO 220 M=1,NT VDIFCALT.111
DO 220 I=1,IMT VDIFCALT.112
ff(I,KM+1,M)=fxc VDIFCALT.113
220 CONTINUE VDIFCALT.114
C VDIFCALT.115
IF (L_OVARYT) THEN ORH1F305.5391
DO K = 1, KM ORH1F305.5392
DO I = 1, IMT ORH1F305.5393
! Adjust if using variable timestep with depth ORH1F305.5394
tempa(I,K)=DTTSA(K)*DZ2RQ(I,K)*fxa*FM(I,K) ORH1F305.5395
tempb(I,K)=gnu(I,K)*DZZ2RQ(I,K) ORH1F305.5396
ENDDO ! over I ORH1F305.5397
ENDDO ! over K ORH1F305.5398
ELSE ORH1F305.5399
DO K = 1, KM ORH1F305.5400
DO I = 1, IMT ORH1F305.5401
tempa(I,K)=C2DTTS*DZ2RQ(I,K)*fxa*FM(I,K) ORH1F305.5402
tempb(I,K)=gnu(I,K)*DZZ2RQ(I,K) ORH1F305.5403
ENDDO ! over I ORH1F305.5404
ENDDO ! over K ORH1F305.5405
ENDIF ORH1F305.5406
DO 240 K=2,KM VDIFCALT.129
DO 240 I=1,IMT VDIFCALT.130
aa(I,K)=tempa(I,K)*tempb(I,K) VDIFCALT.131
240 CONTINUE VDIFCALT.132
DO 250 K=1,KMM1 VDIFCALT.133
DO 250 I=1,IMT VDIFCALT.134
cc(I,K)=tempa(I,K)*tempb(I,K+1)*FM(I,K+1) VDIFCALT.135
250 CONTINUE VDIFCALT.136
DO 260 K=1,KM VDIFCALT.137
DO 260 I=1,IMT VDIFCALT.138
bb(I,K)=aa(I,K)+cc(I,K)+fxb VDIFCALT.139
260 CONTINUE VDIFCALT.140
C VDIFCALT.141
C OOM1F405.871
C COUNTER GRADIENT FLUX TERMS ARE FROM THE FULL LARGE SCHEME. OOM1F405.872
C SET UP EXTRA EXPLICIT TERMS FROM CGFLUX TO INCLUDE. OOM1F405.873
C BOUNDARY CONDITIONS IMPLY NEVER ANY CGFLUX ACROSS OCEAN BED. OOM1F405.874
C OOM1F405.875
DO M=1,NT OOM1F405.876
DO K=1,KMM1 OOM1F405.877
DO I=1,IMT OOM1F405.878
CGTERM(I,K,M)=(FM(I,K)*FXD)*(C2DTTS*DZ2RQ(I,K))* OOM1F405.879
& (CGFLUX(I,K,M)-CGFLUX(I,K+1,M)) OOM1F405.880
ENDDO OOM1F405.881
ENDDO OOM1F405.882
ENDDO OOM1F405.883
C OOM1F405.884
K=KM OOM1F405.885
DO M=1,NT OOM1F405.886
DO I=1,IMT OOM1F405.887
CGTERM(I,K,M)=0.0 OOM1F405.888
ENDDO OOM1F405.889
ENDDO OOM1F405.890
C OOM1F405.891
C -------------------------------------------------------------- VDIFCALT.142
CL The terms x(k-1) (contained in ee) and y(k-1) (contained in ff) VDIFCALT.143
CL in Equation (1.18) - and hence in Equation (1.20) - are calculated. VDIFCALT.144
C -------------------------------------------------------------- VDIFCALT.145
C VDIFCALT.146
DO 290 K=KM,1,-1 VDIFCALT.147
DO 270 I=2,IMTM1 VDIFCALT.148
efdr(I)=fxb/(bb(I,K)-cc(I,K)*ee(I,K+1)) VDIFCALT.149
270 CONTINUE VDIFCALT.150
DO 280 I=2,IMTM1 VDIFCALT.151
ee(I,K)=aa(I,K)*efdr(I) VDIFCALT.152
280 CONTINUE VDIFCALT.153
DO 290 M=1,NT VDIFCALT.154
DO 290 I=2,IMTM1 VDIFCALT.155
C ff(I,K,M)=efdr(I)*( TA(I,K,M)+(cc(I,K)*ff(I,K+1,M)) ) OOM1F405.892
ff(I,K,M)=(TA(I,K,M)+(cc(I,K)*ff(I,K+1,M)))*efdr(I) OOM1F405.893
& -(CGTERM(I,K,M)*efdr(I)) OOM1F405.894
C OOM1F405.895
290 CONTINUE VDIFCALT.157
C VDIFCALT.158
C -------------------------------------------------------------- VDIFCALT.159
CL The new tracer quantity for the top level model is VDIFCALT.160
CL calculated using Equation (1.22). VDIFCALT.161
CL Variable bb is re-used, this time with dimensions of the tracer. VDIFCALT.162
C -------------------------------------------------------------- VDIFCALT.163
C VDIFCALT.164
DO 350 M=1,NT VDIFCALT.165
DO 320 I=2,IMTM1 VDIFCALT.166
bb(I,1)= ff(I,1,M) VDIFCALT.167
320 CONTINUE VDIFCALT.168
C VDIFCALT.169
C -------------------------------------------------------------- VDIFCALT.170
CL The new tracer quantities terms are calculated for the rest of VDIFCALT.171
CL the model levels, using equation (1.20). VDIFCALT.172
C -------------------------------------------------------------- VDIFCALT.173
C VDIFCALT.174
DO 330 K=2,KM VDIFCALT.175
DO 330 I=2,IMTM1 VDIFCALT.176
bb(I,K)=ee(I,K)*bb(I,K-1)+ff(I,K,M) VDIFCALT.177
330 CONTINUE VDIFCALT.178
DO 340 K=1,KM VDIFCALT.179
DO 340 I=1,IMT VDIFCALT.180
TA(I,K,M)=bb(I,K)*FM(I,K) VDIFCALT.181
340 CONTINUE VDIFCALT.182
350 CONTINUE VDIFCALT.183
C VDIFCALT.184
RETURN VDIFCALT.185
END VDIFCALT.186
*ENDIF @DYALLOC.4678