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

Contents 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 - (show annotations) (download)
Thu Nov 18 04:00:05 2010 UTC (15 years, 1 month 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 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