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

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

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

Revision 1.2 - (show annotations) (download)
Tue Dec 10 03:00:59 2002 UTC (21 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint47e_post, checkpoint47c_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47i_post, checkpoint47d_post, checkpoint47g_post, checkpoint48a_post, checkpoint47j_post, branch-exfmods-tag, checkpoint48c_post, checkpoint47f_post, checkpoint48, checkpoint47h_post
Branch point for: branch-exfmods-curt
Changes since 1.1: +31 -1 lines 
 * OCEANICP & realFreshWater: include P-E direct effect on wVel ;
   NOTES: requires option NONLIN_FRSURF to be "#define".
1 C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.1 2002/10/07 16:17:09 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: INTEGR_CONTINUITY
8 C !INTERFACE:
9 SUBROUTINE INTEGR_CONTINUITY(
10 I bi, bj, uFld, vFld,
11 I myTime, myIter, myThid )
12 C !DESCRIPTION: \bv
13 C *==========================================================*
14 C | SUBROUTINE INTEGR_CONTINUITY
15 C | o Integrate the continuity Eq : compute vertical velocity
16 C | and free surface "r-anomaly" (etaN) to satisfy
17 C | exactly the convervation of the Total Volume
18 C *==========================================================*
19 C \ev
20
21 C !USES:
22 IMPLICIT NONE
23 C == Global variables
24 #include "SIZE.h"
25 #include "EEPARAMS.h"
26 #include "PARAMS.h"
27 #include "DYNVARS.h"
28 #include "GRID.h"
29 #include "SURFACE.h"
30 #include "FFIELDS.h"
31
32 C !INPUT/OUTPUT PARAMETERS:
33 C == Routine arguments ==
34 C uFld :: Zonal velocity ( m/s )
35 C vFld :: Meridional velocity ( m/s )
36 C bi,bj :: tile index
37 C myTime :: Current time in simulation
38 C myIter :: Current iteration number in simulation
39 C myThid :: Thread number for this instance of the routine.
40 _RL myTime
41 INTEGER myIter
42 INTEGER myThid
43 INTEGER bi,bj
44 LOGICAL UpdateEtaN_EtaH
45 _RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46 _RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
47
48 C !LOCAL VARIABLES:
49 C Local variables in common block
50
51 C Local variables
52 C i,j,k :: Loop counters
53 C uTrans :: Volume transports ( uVel.xA )
54 C vTrans :: Volume transports ( vVel.yA )
55 INTEGER i,j,k
56 _RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 _RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58 _RL wSurf
59 CEOP
60
61 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
62
63 #ifdef EXACT_CONSERV
64 IF (exactConserv) THEN
65
66 C-- Compute the Divergence of The Barotropic Flow :
67
68 C- Initialise
69 DO j=1-Oly,sNy+Oly
70 DO i=1-Olx,sNx+Olx
71 hDivFlow(i,j,bi,bj) = 0. _d 0
72 utrans(i,j) = 0. _d 0
73 vtrans(i,j) = 0. _d 0
74 ENDDO
75 ENDDO
76
77 DO k=1,Nr
78
79 C- Calculate velocity field "volume transports" through tracer cell faces
80 DO j=1,sNy+1
81 DO i=1,sNx+1
82 uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
83 & *drF(k)*_hFacW(i,j,k,bi,bj)
84 vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
85 & *drF(k)*_hFacS(i,j,k,bi,bj)
86 ENDDO
87 ENDDO
88
89 C- Integrate vertically the Horizontal Divergence
90 DO j=1,sNy
91 DO i=1,sNx
92 hDivFlow(i,j,bi,bj) = hDivFlow(i,j,bi,bj)
93 & +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
94 & +vTrans(i,j+1)-vTrans(i,j) )
95 ENDDO
96 ENDDO
97
98 C- End DO k=1,Nr
99 ENDDO
100
101 ENDIF
102
103 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104
105 IF (exactConserv .AND. myTime.NE.startTime) THEN
106 C-- Update etaN at the end of the time step :
107 C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
108
109 IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
110 DO j=1-Oly,sNy+Oly
111 DO i=1-Olx,sNx+Olx
112 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
113 ENDDO
114 ENDDO
115 ELSEIF (useRealFreshWaterFlux) THEN
116 DO j=1,sNy
117 DO i=1,sNx
118 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
119 & - implicDiv2Dflow*( convertEmP2rUnit*EmPmR(i,j,bi,bj)
120 & +hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
121 & )*deltaTfreesurf
122 ENDDO
123 ENDDO
124 ELSE
125 DO j=1,sNy
126 DO i=1,sNx
127 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
128 & - implicDiv2Dflow*hDivFlow(i,j,bi,bj)
129 & *recip_rA(i,j,bi,bj)*deltaTfreesurf
130 ENDDO
131 ENDDO
132 ENDIF
133
134 ENDIF
135 #endif /* EXACT_CONSERV */
136
137 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
138
139 DO k=Nr,1,-1
140 C-- Integrate continuity vertically for vertical velocity
141
142 CALL INTEGRATE_FOR_W(
143 I bi, bj, k, uFld, vFld,
144 O wVel,
145 I myThid )
146
147 #ifdef NONLIN_FRSURF
148 IF ( k.EQ.Nr .AND. myTime.NE. 0. _d 0 .AND.
149 & useRealFreshWaterFlux .AND.
150 & buoyancyRelation .EQ. 'OCEANICP' ) THEN
151
152 IF ( myTime.NE.startTime ) THEN
153 DO j=1,sNy
154 DO i=1,sNx
155 wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
156 & -convertEmP2rUnit*EmPmR(i,j,bi,bj)*maskC(i,j,k,bi,bj)
157 ENDDO
158 ENDDO
159 ELSEIF (nonlinFreeSurf.GE.0) THEN
160 C needs previous time-step value of E-P-R, that has not been loaded
161 C and was not in pickup-file ; try to use etaN & etaH instead.
162 DO j=1,sNy
163 DO i=1,sNx
164 wSurf = hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
165 & + (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
166 & /(implicDiv2Dflow*deltaTfreesurf)
167 wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
168 & + wSurf*maskC(i,j,k,bi,bj)
169 ENDDO
170 ENDDO
171 ENDIF
172
173 ENDIF
174 #endif /* NONLIN_FRSURF */
175
176 #ifdef ALLOW_OBCS
177 #ifdef ALLOW_NONHYDROSTATIC
178 C-- Apply OBC to W if in N-H mode
179 IF (useOBCS.AND.nonHydrostatic) THEN
180 CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
181 ENDIF
182 #endif /* ALLOW_NONHYDROSTATIC */
183 #endif /* ALLOW_OBCS */
184
185 C- End DO k=Nr,1,-1
186 ENDDO
187
188 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
189
190 RETURN
191 END
  ViewVC Help
Powered by ViewVC 1.1.22