/[MITgcm]/MITgcm_contrib/bbl/code/mypackage_calc_rhs.F
ViewVC logotype

Annotation of /MITgcm_contrib/bbl/code/mypackage_calc_rhs.F

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

Revision 1.1 - (hide annotations) (download)
Thu Nov 18 04:00:05 2010 UTC (14 years, 7 months ago) by dimitri
Branch: MAIN 
This is a first sketch of a bottom boundary layer parameterization
for MITgcm.  The hooks to main model currently reside with pkg/mypackage
and it is temporarily checked in MITgcm_contrib until it clears the
App Store vetting process.  Instructions on running a simple test
integration in a periodic channel are in MITgcm_contrib/bbl/readme.txt
and some output can be viewed using lookat_output.m in same directory.
1 dimitri 1.1 C $Header: $
2     C $Name: $
3    
4     #include "BBL_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: BBL_CALC_RHS
8    
9     C !INTERFACE:
10     SUBROUTINE MYPACKAGE_CALC_RHS(
11     I bi, bj, myTime, myIter, myThid )
12    
13     C !DESCRIPTION:
14     C Calculate tendency of tracers due to bottom boundary layer.
15     C Scheme is currently coded for dTtracerLev(k) .EQ. deltaT.
16    
17     C !USES:
18     IMPLICIT NONE
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "GRID.h"
23     #include "DYNVARS.h"
24     #include "BBL.h"
25    
26     C !INPUT PARAMETERS:
27     C bi,bj :: Tile indices
28     C myTime :: Current time in simulation
29     C myIter :: Current time-step number
30     C myThid :: my Thread Id number
31     INTEGER bi, bj
32     _RL myTime
33     INTEGER myIter, myThid
34    
35     C !OUTPUT PARAMETERS:
36    
37     C !LOCAL VARIABLES:
38     C i,j :: Loop indices
39     C kBot :: k index of bottommost wet grid box
40     C resThk :: thickness of bottommost wet grid box minus bbl_eta
41     C resTheta :: temperature of this residual volume
42     C resSalt :: salinity of this residual volume
43     C deltaRho :: density change
44     C deltaDpt :: depth change
45     C bbl_tend :: temporary variable for tendency terms
46     C sloc :: salinity of bottommost wet grid box
47     C tloc :: temperature of bottommost wet grid box
48     C rholoc :: Potential density of bottommost wet grid box
49     C bbl_rho :: Potential density of bottom boundary layer
50     C fZon :: Zonal flux
51     C fMer :: Meridional flux
52     INTEGER i, j, kBot
53     _RL resThk, resTheta, resSalt
54     _RL deltaRho, deltaDpt, bbl_tend
55     _RL sloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
56     _RL tloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
57     _RL rholoc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
58     _RL bbl_rho ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
59     _RL fZon ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
60     _RL fMer ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
61     CEOP
62    
63     C Initialize tendency terms and
64     C make local copy of bottomost temperature and salinity
65     DO j=1-OLy,sNy+OLy
66     DO i=1-OLx,sNx+OLx
67     bbl_TendTheta(i,j,bi,bj) = 0. _d 0
68     bbl_TendSalt (i,j,bi,bj) = 0. _d 0
69     kBot = max(1,kLowC(i,j,bi,bj))
70     tLoc(i,j) = theta(i,j,kBot,bi,bj)
71     sLoc(i,j) = salt (i,j,kBot,bi,bj)
72     ENDDO
73     ENDDO
74    
75     C Calculate potential density of bottommost wet grid box
76     CALL FIND_RHO_2D(
77     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
78     I tLoc, sLoc,
79     O rhoLoc,
80     I 1, bi, bj, myThid )
81    
82     C Calculate potential density of bbl
83     CALL FIND_RHO_2D(
84     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
85     I bbl_theta(1-OLx,1-OLy,bi,bj), bbl_salt(1-OLx,1-OLy,bi,bj),
86     O bbl_rho,
87     I 1, bi, bj, myThid )
88    
89     C==== Compute and apply vertical exchange between BBL and
90     C residual volume of botommost wet grid box.
91     C This operation does not change total tracer quantity
92     C in botommost wet grid box.
93    
94     DO j=1-Oly,sNy+Oly
95     DO i=1-Olx,sNx+Olx
96     kBot = kLowC(i,j,bi,bj)
97     IF ( kBot .GT. 0 ) THEN
98     resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj)
99    
100     C If bbl occupies the complete bottom model grid box or
101     C if model density is higher than BBL then mix instantly.
102     IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.bbl_rho(i,j)) ) THEN
103     bbl_theta(i,j,bi,bj) = tLoc(i,j)
104     bbl_salt (i,j,bi,bj) = sLoc(i,j)
105    
106     C If model density is lower than BBL, slowly diffuse upward.
107     ELSE
108     resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
109     & (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
110     resSalt = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
111     & (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
112     bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
113     & deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR
114     bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
115     & deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR
116     ENDIF
117     ENDIF
118     ENDDO
119     ENDDO
120    
121     C==== Compute zonal bbl exchange.
122     DO j=1-Oly,sNy+Oly
123     DO i=1-Olx,sNx+Olx-1
124     IF ((kLowC(i,j,bi,bj).GT.0).AND.(kLowC(i+1,j,bi,bj).GT.0)) THEN
125     deltaRho = bbl_rho(i+1,j) - bbl_rho(i,j)
126     deltaDpt = R_low(i ,j,bi,bj) + bbl_eta(i ,j,bi,bj) -
127     & R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
128    
129     C If heavy BBL water is higher than light BBL water,
130     C exchange properties laterally.
131     IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
132     bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
133     & ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) /
134     & bbl_RelaxH
135     bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) -
136     & ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) *
137     & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
138     & ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
139     bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
140     & ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) /
141     & bbl_RelaxH
142     bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) -
143     & ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) *
144     & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
145     & ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
146     ENDIF
147     ENDIF
148     ENDDO
149     ENDDO
150    
151     C==== Compute meridional bbl exchange.
152     DO j=1-Oly,sNy+Oly-1
153     DO i=1-Olx,sNx+Olx
154     IF ((kLowC(i,j,bi,bj).GT.0).AND.(kLowC(i,j+1,bi,bj).GT.0)) THEN
155     deltaRho = bbl_rho(i,j+1) - bbl_rho(i,j)
156     deltaDpt = R_low(i,j ,bi,bj) + bbl_eta(i,j ,bi,bj) -
157     & R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
158    
159     C If heavy BBL water is higher than light BBL water,
160     C exchange properties laterally.
161     IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
162     bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
163     & ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) /
164     & bbl_RelaxH
165     bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) -
166     & ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) *
167     & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
168     & ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) /
169     & bbl_RelaxH
170     bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
171     & ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) /
172     & bbl_RelaxH
173     bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) -
174     & ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) *
175     & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
176     & ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH )
177     ENDIF
178     ENDIF
179     ENDDO
180     ENDDO
181    
182     C==== Apply lateral BBL exchange then scale tendency term
183     C for botommost wet grid box.
184     DO j=1-Oly,sNy+Oly-1
185     DO i=1-Olx,sNx+Olx-1
186     kBot = kLowC(i,j,bi,bj)
187     IF ( kBot .GT. 0 ) THEN
188     bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
189     & deltaT * bbl_TendTheta(i,j,bi,bj)
190     bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
191     & deltaT * bbl_TendSalt (i,j,bi,bj)
192     bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) *
193     & bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
194     bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) *
195     & bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
196     ENDIF
197     ENDDO
198     ENDDO
199    
200     CALL EXCH_XY_RL( bbl_theta, myThid )
201     CALL EXCH_XY_RL( bbl_salt , myThid )
202    
203     RETURN
204     END
  ViewVC Help
Powered by ViewVC 1.1.22