*IF DEF,OCEAN ISOPYC_M.2
C ******************************COPYRIGHT****************************** ISOPYC_M.3
C (c) CROWN COPYRIGHT 1998, METEOROLOGICAL OFFICE, All Rights Reserved. ISOPYC_M.4
C ISOPYC_M.5
C Use, duplication or disclosure of this code is subject to the ISOPYC_M.6
C restrictions as set forth in the contract. ISOPYC_M.7
C ISOPYC_M.8
C Meteorological Office ISOPYC_M.9
C London Road ISOPYC_M.10
C BRACKNELL ISOPYC_M.11
C Berkshire UK ISOPYC_M.12
C RG12 2SZ ISOPYC_M.13
C ISOPYC_M.14
C If no contract has been raised with this copy of the code, the use, ISOPYC_M.15
C duplication or disclosure of it is strictly prohibited. Permission ISOPYC_M.16
C to do so must first be obtained in writing from the Head of Numerical ISOPYC_M.17
C Modelling at the above address. ISOPYC_M.18
C ******************************COPYRIGHT****************************** ISOPYC_M.19
C ISOPYC_M.20
CLL Subroutine ISOPYC_M ISOPYC_M.21
CLL ISOPYC_M.22
CLL Author: M J Roberts ISOPYC_M.23
CLL ISOPYC_M.24
CLL Description: Mirrors subroutine IPDCOFCL for the Griffies ISOPYC_M.25
CLL isopycnal diffusion. Sets up the diffusion matrices ISOPYC_M.26
CLL ready for isopycnal diffusion, done in ISOFLUX ISOPYC_M.27
CLL ISOPYC_M.28
CLL External documentation: ISOPYC_M.29
CLL Modular Ocean Model 2 manual and code ISOPYC_M.30
CLL ISOPYC_M.31
CLL Implemented at UM vn 4.5 26 June 1998 ISOPYC_M.32
CLL ISOPYC_M.33
CLL Modification History : ISOPYC_M.34
CLL Version Date Comment & Name ISOPYC_M.35
CLL ------- -------- -------------------------------------------- ISOPYC_M.36
CLL ISOPYC_M.37
CLL Subroutine dependencies. ISOPYC_M.38
CLL Called by TRACER ISOPYC_M.39
CLL Calls subroutines ELEMENTS, AI_CALC, ISOPYC_A ISOPYC_M.40
CLL ISOPYC_M.41
CLLEND------------------------------------------------------------------ ISOPYC_M.42
C ISOPYC_M.43
SUBROUTINE ISOPYC_M( 2,9ISOPYC_M.44
*CALL ARGSIZE 
ISOPYC_M.45
*CALL COCAWRKA ISOPYC_M.46
& ,j,cstr,dyur,dxur,dz2r,dzz,dzz2r,athkdftu,athkdftv,ahi,athkdf, ISOPYC_M.47
& slope_max,dz,dyu,dxu,csu,dxt4r,dyt4r,slopec,dslope,dxtr,dytr, ISOPYC_M.48
& cst,csr,Ai_ez,Ai_nz,Ai_bx,Ai_by,K11,K22,K33,gnu,adv_vetiso, ISOPYC_M.49
& adv_vbtiso,adv_fbiso,J_OFFSET,athkdf_bi, ISOPYC_M.50
& csjm,dyurjm,j_1,srho) ISOPYC_M.51
C Input variables ISOPYC_M.52
ISOPYC_M.53
*CALL OARRYSIZ ISOPYC_M.54
*CALL TYPSIZE ISOPYC_M.55
*CALL CNTLOCN ISOPYC_M.56
*CALL COCTWRKA ISOPYC_M.57
*CALL OTIMER ISOPYC_M.58
ISOPYC_M.59
INTEGER i,k,ip,kr,j,jq,J_OFFSET,j_1 ISOPYC_M.60
ISOPYC_M.61
C Input variables ISOPYC_M.62
REAL slmxr,dtxsqr(km),dzt2r(km),dzz2r(km+1), ISOPYC_M.63
& dzw(0:km),csu(jmt),dyu(jmt),dxu(imt), ISOPYC_M.64
& cst(jmt),dytr(jmt),cstr(jmt),dxt4r(imt),dyt4r(jmt), ISOPYC_M.65
& slopec,dslope,dzwr(0:km),dz2r(km),dzz(km+1), ISOPYC_M.66
& dxtr(imt),dztr(km),dxur(imt),dyur(jmt),dz(km), ISOPYC_M.67
& gnu(imt,km),csr(jmt),ahi(km),athkdf(km) ISOPYC_M.68
REAL athkdftu(imt_vis,jmt_vis),athkdftv(imt_vis,jmt_vis) ISOPYC_M.69
REAL csjm,dyurjm,athkdf_bi ISOPYC_M.70
ISOPYC_M.71
C full tensor variables ISOPYC_M.72
REAL delta_iso,s_minus,s_plus ISOPYC_M.73
ISOPYC_M.74
C isopycnal diffusion variables defined below ISOPYC_M.75
REAL Ai_ez(imt_iso,km_iso,0:1,0:1),Ai_bx(imt_iso,km_iso,0:1,0:1), ISOPYC_M.76
& Ai_by(imt_iso,km_iso,0:1,0:1),Ai_nz(imt_iso,km_iso,0:1,0:1), ISOPYC_M.77
& K11(imt_iso,km_iso),K22(imt_iso,km_iso),K33(imt_iso,km_iso) ISOPYC_M.78
ISOPYC_M.79
C GM variables ISOPYC_M.80
REAL adv_vetiso(imt_gmm,km_gmm), ISOPYC_M.81
& adv_vbtiso(imt_gmm,0:km_gmm),adv_fbiso(imt_gmm,0:km_gmm) ISOPYC_M.82
REAL srho(imt,km) ISOPYC_M.83
ISOPYC_M.84
C ISOPYC_M.85
C calls subroutines: ISOPYC_M.86
C elements - initialises the gradients and expansion ISOPYC_M.87
C coefficients and updates the row ISOPYC_M.88
C ISOPYC_M.89
C ai_calc - calculates the diffusion coefficients ISOPYC_M.90
C Ai ISOPYC_M.91
C ISOPYC_M.92
C isopyc_a - calculates the advection velocities due ISOPYC_M.93
C to the Gent and McWilliams scheme ISOPYC_M.94
c ISOPYC_M.95
c isopycnal diffusion variables: ISOPYC_M.96
c ISOPYC_M.97
c ahi(km) = isopycnal tracer mixing coefficient (cm**2/sec) ISOPYC_M.98
c drodx = d(rho)/dx local to east face of T cell ISOPYC_M.99
c drody = d(rho)/dy local to north face of T cell ISOPYC_M.100
c drodz = d(rho)/dz local to bottom face of T cell ISOPYC_M.101
c Ai_e = diffusion coefficient on eastern face of T cell ISOPYC_M.102
c Ai_n = diffusion coefficient on northern face of T cell ISOPYC_M.103
c Ai_bx = diffusion coefficient on bottom face of T cell ISOPYC_M.104
c Ai_by = diffusion coefficient on bottom face of T cell ISOPYC_M.105
c ISOPYC_M.106
c K11 = (xx) component of isopycnal diffusion tensor ISOPYC_M.107
c K22 = (yy) component of isopycnal diffusion tensor ISOPYC_M.108
c K33 = (zz) component of isopycnal diffusion tensor ISOPYC_M.109
c ISOPYC_M.110
c slmxr = reciprocal of maximum allowable slope of isopycnals for ISOPYC_M.111
c small angle approximation ISOPYC_M.112
c ISOPYC_M.113
c slopec = } Values used to control amount of tapering ISOPYC_M.114
c dslope = } of isopycnal slopes with option L_OISOTAPER ISOPYC_M.115
ISOPYC_M.116
c if OISOGM ISOPYC_M.117
c adv_vetiso = zonal isopycnal mixing velocity computed at the ISOPYC_M.118
c center of the eastern face of the "t" cells ISOPYC_M.119
c adv_vntiso = meridional isopycnal mixing velocity computed at ISOPYC_M.120
c the center of the northern face of the "t" cells ISOPYC_M.121
c (Note: this includes the cosine as in "adv_vnt") ISOPYC_M.122
c adv_vbtiso = vertical isopycnal mixing velocity computed at the ISOPYC_M.123
c center of the top face of the "t" cells ISOPYC_M.124
c adv_fbiso = "adv_vbtiso" * (tracer) evaluated at the center of ISOPYC_M.125
c the bottom face of the "t" cells ISOPYC_M.126
c athkdf = isopycnal thickness diffusivity (cm**2/sec) ISOPYC_M.127
c endif ISOPYC_M.128
ISOPYC_M.129
c ISOPYC_M.130
c----------------------------------------------------------------------- ISOPYC_M.131
c set local constants ISOPYC_M.132
c----------------------------------------------------------------------- ISOPYC_M.133
C MOM variables set relative to Cox variables ISOPYC_M.134
do k=0,km ISOPYC_M.135
dzwr(k)=1./dzz(k+1) ISOPYC_M.136
dzw(k)=dzz(k+1) ISOPYC_M.137
enddo ISOPYC_M.138
ISOPYC_M.139
do k=1,km ISOPYC_M.140
dzt2r(k)=dz2r(k) ISOPYC_M.141
dztr(k)=2.*dz2r(k) ISOPYC_M.142
enddo ISOPYC_M.143
ISOPYC_M.144
c ISOPYC_M.145
c----------------------------------------------------------------------- ISOPYC_M.146
c estimate alpha, beta, and gradients on sides of T cells for ISOPYC_M.147
c the next row and updates the current row ISOPYC_M.148
c----------------------------------------------------------------------- ISOPYC_M.149
c ISOPYC_M.150
IF (L_OTIMER) CALL TIMER('ELEMENTS',103) ISOPYC_M.151
ISOPYC_M.152
call elements ( ISOPYC_M.153
*CALL ARGSIZE ISOPYC_M.154
*CALL COCAWRKA ISOPYC_M.155
& ,j,cstr,dyur,dxur,dz2r,dzz2r,dzw,dtxsqr,delta_iso, ISOPYC_M.156
& s_minus,s_plus,dzt2r,Ai_ez,Ai_nz,Ai_bx,Ai_by,K11,K22,K33, ISOPYC_M.157
& slope_max,slmxr,adv_vetiso,adv_vbtiso,adv_fbiso,J_OFFSET) ISOPYC_M.158
ISOPYC_M.159
IF (L_OTIMER) CALL TIMER('ELEMENTS',104) ISOPYC_M.160
c ISOPYC_M.161
c----------------------------------------------------------------------- ISOPYC_M.162
c compute Ai_ez centered on eastern face of T cells ISOPYC_M.163
c compute Ai_nz centered on the northern face of T cells ISOPYC_M.164
c compute Ai_bx & Ai_by centered on bottom face of T cells ISOPYC_M.165
c----------------------------------------------------------------------- ISOPYC_M.166
c ISOPYC_M.167
IF (L_OTIMER) CALL TIMER('AI_CALC',103) ISOPYC_M.168
ISOPYC_M.169
call ai_calc ( ISOPYC_M.170
*CALL ARGSIZE ISOPYC_M.171
*CALL COCAWRKA ISOPYC_M.172
& ,j,cstr,dyur,dxur,dxtr,dz2r,ahi,csr,csjm,dyurjm,j_1, ISOPYC_M.173
& slmxr,dtxsqr,dyu,dxu,dzt2r,dzw,csu,dxt4r,dyt4r,dytr,cst, ISOPYC_M.174
& Ai_ez,Ai_nz,Ai_bx,Ai_by,K11,K22,K33,gnu(1,1), ISOPYC_M.175
& delta_iso,s_plus,s_minus,athkdftu,athkdftv,athkdf) ISOPYC_M.176
ISOPYC_M.177
IF (L_OTIMER) CALL TIMER('AI_CALC',104) ISOPYC_M.178
c ISOPYC_M.179
c ISOPYC_M.180
IF (L_OISOGM) THEN ISOPYC_M.181
c ISOPYC_M.182
c----------------------------------------------------------------------- ISOPYC_M.183
c compute isopycnal advective velocities for tracers ISOPYC_M.184
c----------------------------------------------------------------------- ISOPYC_M.185
c ISOPYC_M.186
ISOPYC_M.187
IF (L_OTIMER) CALL TIMER('ISOPYC_A',103) ISOPYC_M.188
ISOPYC_M.189
call isopyc_a ( ISOPYC_M.190
*CALL ARGSIZE ISOPYC_M.191
*CALL COCAWRKA ISOPYC_M.192
& ,j,dtxsqr,slmxr,dz2r,csu,dz,athkdftu,athkdftv,athkdf, ISOPYC_M.193
& cstr,dxtr,dytr,adv_vetiso,adv_vbtiso,slopec,dslope,dxur,dyur, ISOPYC_M.194
& cst,csr,J_OFFSET,athkdf_bi,srho) ISOPYC_M.195
c ISOPYC_M.196
IF (L_OTIMER) CALL TIMER('ISOPYC_A',104) ISOPYC_M.197
ISOPYC_M.198
ENDIF ISOPYC_M.199
c ISOPYC_M.200
ISOPYC_M.201
RETURN ISOPYC_M.202
END ISOPYC_M.203
ISOPYC_M.204
*ENDIF ISOPYC_M.205
ISOPYC_M.206
ISOPYC_M.207