/[MITgcm]/MITgcm/model/src/calc_isoslopes.F
ViewVC logotype

Contents of /MITgcm/model/src/calc_isoslopes.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.16 - (show annotations) (download)
Wed Jun 21 18:59:27 2000 UTC (23 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: HEAD
Changes since 1.15: +1 -1 lines
FILE REMOVED
Removed calc_isoslopes as part of packaging of GM/Redi scheme.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_isoslopes.F,v 1.15 2000/06/09 02:45:04 heimbach Exp $
2 #include "CPP_OPTIONS.h"
3
4 CStartOfInterface
5 SUBROUTINE CALC_ISOSLOPES(
6 I bi, bj, iMin, iMax, jMin, jMax, K,
7 I rhoKm1, rhoKp1, rhotmp,
8 O K13, K23, K33, KapGM,
9 I myThid )
10 C /==========================================================\
11 C | SUBROUTINE CALC_ISOSLOPES |
12 C | o Calculate isentropic surface slopes for GM scheme. |
13 C |==========================================================|
14 C | The small-angle approixmation of the full-slope tensor |
15 C | is used here. A maximum slope is also imposed to prevent |
16 C | unphysical values in homogenised regions. |
17 C | This routine works on a given XY slab of fluid. |
18 C \==========================================================/
19 IMPLICIT NONE
20
21 C == Global variables ==
22 #include "SIZE.h"
23 #include "GRID.h"
24 #include "DYNVARS.h"
25 #include "EEPARAMS.h"
26 #include "PARAMS.h"
27
28 C == Routine arguments ==
29 C
30 _RL rhoKm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
31 _RL rhoKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
32 _RL rhotmp(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
33 _RL K13(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
34 _RL K23(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
35 _RL K33(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
36 _RL KapGM(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37 INTEGER bi,bj,iMin,iMax,jMin,jMax,K
38 INTEGER myThid
39 CEndOfInterface
40
41 C == Local variables ==
42 INTEGER i,j
43 _RL dSigmaDx,dSigmaDy,dSigmaDr,SlopeX,SlopeY,dRdSigma,N2
44 _RL gradSmod,stratLimit,rMaxSlope
45 _RL Small_Number
46 PARAMETER(Small_Number=1. _d -12)
47 _RS deltaH,zero_rs
48 PARAMETER(zero_rs=0.)
49
50 #ifdef ALLOW_AUTODIFF_TAMC
51 DO i=1-OLx,sNx+OLx
52 K13 (i, 1-OLy,k) = 0.0
53 K23 (i, 1-OLy,k) = 0.0
54 K33 (i, 1-OLy,k) = 0.0
55 K13 (i,sNy+OLy,k) = 0.0
56 K23 (i,sNy+OLy,k) = 0.0
57 K33 (i,sNy+OLy,k) = 0.0
58 ENDDO
59 DO j=1-OLy,sNy+OLy
60 K13 ( 1-OLx,j,k) = 0.0
61 K23 ( 1-OLx,j,k) = 0.0
62 K33 ( 1-OLx,j,k) = 0.0
63 K13 (sNx+OLx,j,k) = 0.0
64 K23 (sNx+OLx,j,k) = 0.0
65 K33 (sNx+OLx,j,k) = 0.0
66 ENDDO
67 #endif
68
69 C Calculate terms for mean Richardson number
70 C which is used in the "variable K" parameterisaton.
71 C Distance between interface above layer and the integration depth (GMdepth)
72 deltaH=abs(GMdepth)-abs(rF(k))
73 C If positive we limit this to the layer thickness
74 deltaH=min(deltaH,drF(k))
75 C If negative then we are below the integration level
76 deltaH=max(deltaH,zero_rs)
77 C Now we convert deltaH to a non-dimensional fraction
78 deltaH=deltaH/GMdepth
79
80 rMaxSlope=0.
81 if (GMmaxSlope.ne.0.) rMaxSlope=1. _d 0 / GMmaxSlope
82
83 #ifdef ALLOW_AUTODIFF_TAMC
84 !HPF$ INDEPENDENT
85 #endif
86 DO j=1-Oly+1,sNy+Oly-1
87 #ifdef ALLOW_AUTODIFF_TAMC
88 !HPF$ INDEPENDENT
89 #endif
90 DO i=1-Olx+1,sNx+Olx-1
91
92 C Horizontal gradient of Sigma at rVel points
93 dSigmaDx=0.25*
94 & (_maskW(i+1,j, k ,bi,bj)
95 & *_recip_dxC(i+1,j,bi,bj)
96 & *(rhoKp1(i+1,j)-rhoKp1( i ,j))
97 & +_maskW( i ,j, k ,bi,bj)
98 & *_recip_dxC(i,j,bi,bj)
99 & *(rhoKp1( i ,j)-rhoKp1(i-1,j))
100 & +_maskW(i+1,j,k-1,bi,bj)
101 & *_recip_dxC(i+1,j,bi,bj)
102 & *(rhoKm1(i+1,j)-rhoKm1( i ,j))
103 & +_maskW( i ,j,k-1,bi,bj)
104 & *_recip_dxC(i,j,bi,bj)
105 & *(rhoKm1( i ,j)-rhoKm1(i-1,j)) )
106 dSigmaDy=0.25*
107 & (_maskS(i,j+1, k ,bi,bj)
108 & *_recip_dyC(i,j+1,bi,bj)
109 & *(rhoKp1(i,j+1)-rhoKp1(i, j ))
110 & +_maskS(i, j, k ,bi,bj)
111 & *_recip_dyC(i,j,bi,bj)
112 & *(rhoKp1(i, j )-rhoKp1(i,j-1))
113 & +_maskS(i,j+1,k-1,bi,bj)
114 & *_recip_dyC(i,j+1,bi,bj)
115 & *(rhoKm1(i,j+1)-rhoKm1(i, j ))
116 & +_maskS(i, j ,k-1,bi,bj)
117 & *_recip_dyC(i,j,bi,bj)
118 & *(rhoKm1(i, j )-rhoKm1(i,j-1)) )
119
120 C Calculate magnitude of horizontal gradient
121 #ifdef ALLOW_AUTODIFF_TAMC
122 if (dSigmaDx*dSigmaDx+dSigmaDy*dSigmaDy .eq. 0) then
123 gradSmod=0.0
124 else
125 gradSmod=sqrt(dSigmaDx*dSigmaDx+dSigmaDy*dSigmaDy)
126 endif
127 #else
128 gradSmod=sqrt(dSigmaDx*dSigmaDx+dSigmaDy*dSigmaDy)
129 #endif
130
131 C Limit Stratification
132 stratLimit=-Small_Number-gradSmod*rMaxSlope
133 dSigmaDr=recip_drC(k)*rkFac*(rhotmp(i,j)-rhoKp1(i,j))
134 C Altern -> if (dSigmaDz.gt.stratLimit) dSigmaDz=stratLimit
135 dSigmaDr=min(dSigmaDr,stratLimit)
136 N2=(-Gravity*recip_Rhonil)*dSigmaDr
137 dRdSigma=1. _d 0 / dSigmaDr
138
139 C Iso-neutral slopes
140 SlopeX=-dSigmaDx*dRdSigma
141 SlopeY=-dSigmaDy*dRdSigma
142 if (hFacC(i,j,k,bi,bj).eq.0.) then
143 SlopeX=0.
144 SlopeY=0.
145 endif
146
147 C Components of Redi/GM tensor
148 K13(i,j,k)=2. _d 0 * SlopeX
149 K23(i,j,k)=2. _d 0 * SlopeY
150 K33(i,j,k)=SlopeX*SlopeX+SlopeY*SlopeY
151
152 C Depth average of M^2/N^2 * N
153 KapGM(i,j)=KapGM(i,j)+deltaH*
154 & GMalpha*GMlength*GMlength* sqrt(K33(i,j,k)*N2)
155
156 C Limit range that KapGM can take
157 KapGM(i,j)=min(KapGM(i,j),GMmaxval)
158
159 ENDDO
160 ENDDO
161
162 RETURN
163 END

  ViewVC Help
Powered by ViewVC 1.1.22