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

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

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

Revision 1.9 - (show annotations) (download)
Fri Jun 21 18:36:05 2002 UTC (21 years, 10 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint46g_pre, checkpoint46f_post, checkpoint46d_pre, checkpoint46e_post, checkpoint46c_post, checkpoint46b_post, checkpoint46e_pre, checkpoint45d_post, checkpoint46c_pre, checkpoint46, checkpoint46a_post, checkpoint46a_pre, checkpoint46b_pre, checkpoint46d_post, checkpoint46g_post
Changes since 1.8: +6 -6 lines 
Added new parameter: deltaTfreesurf

Previously, the free-surface equation was intergrated forward
synchronously with the momentum equations. It is more consistent
to use the tracer time-step. This increases the number of
iterations required but strengthens the damping.

We *SHOULD* make the default time-step equal to the tracer time-step.
However, we don't for backward compatibility. At some point in the
future we need to change the default behaviour.

It turns out that the reason for the "reduced stability" encountered
in large-scale runs seems to be related to excess variability in
the free surface which in turn happens when the waves aren't damped.
Using a longer time-step fixes this.
1 C $Header: /u/gcmpack/MITgcm/model/src/calc_exact_eta.F,v 1.8 2002/02/10 00:44:37 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: CALC_EXACT_ETA
8 C !INTERFACE:
9 SUBROUTINE CALC_EXACT_ETA( UpdateEtaN_EtaH,
10 I bi,bj, uFld,vFld,
11 I myTime, myIter, myThid )
12 C !DESCRIPTION: \bv
13 C *==========================================================*
14 C | SUBROUTINE CALC_EXACT_ETA
15 C | o Compute again the surface "r-anomaly" (eta) to satisfy
16 C | exactly the convervation of the Total Volume
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 IMPLICIT NONE
22 C == Global variables
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "DYNVARS.h"
27 #include "GRID.h"
28 #include "SURFACE.h"
29 #include "FFIELDS.h"
30
31 C !INPUT/OUTPUT PARAMETERS:
32 C == Routine arguments ==
33 C UpdateEtaN_EtaH :: flag to distinguishe if this S/R is called
34 C at the end of a time step (TRUE) to update EtaN,
35 C or at the beginning of the time sptep (FALSE) to update EtaH
36 C uFld :: Zonal velocity ( m/s )
37 C vFld :: Meridional velocity ( m/s )
38 C bi,bj :: tile index
39 C myTime :: Current time in simulation
40 C myIter :: Current iteration number in simulation
41 C myThid :: Thread number for this instance of the routine.
42 _RL myTime
43 INTEGER myIter
44 INTEGER myThid
45 INTEGER bi,bj
46 LOGICAL UpdateEtaN_EtaH
47 _RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
48 _RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
49
50 C !LOCAL VARIABLES:
51 #ifdef EXACT_CONSERV
52 C Local variables in common block
53
54 C Local variables
55 C i,j,k :: Loop counters
56 C uTrans :: Volume transports ( uVel.xA )
57 C vTrans :: Volume transports ( vVel.yA )
58 INTEGER i,j,k
59 _RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 _RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61 CEOP
62
63 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
64
65 IF ( UpdateEtaN_EtaH .OR. myTime.EQ.startTime ) THEN
66
67 C-- Compute the Divergence of The Barotropic Flow :
68
69 C- Initialise
70 DO j=1-Oly,sNy+Oly
71 DO i=1-Olx,sNx+Olx
72 hDivFlow(i,j,bi,bj) = 0. _d 0
73 utrans(i,j) = 0. _d 0
74 vtrans(i,j) = 0. _d 0
75 ENDDO
76 ENDDO
77
78 DO k=1,Nr
79
80 C- Calculate velocity field "volume transports" through tracer cell faces
81 DO j=1,sNy+1
82 DO i=1,sNx+1
83 uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
84 & *drF(k)*_hFacW(i,j,k,bi,bj)
85 vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
86 & *drF(k)*_hFacS(i,j,k,bi,bj)
87 ENDDO
88 ENDDO
89
90 C- Integrate vertically the Horizontal Divergence
91 DO j=1,sNy
92 DO i=1,sNx
93 hDivFlow(i,j,bi,bj) = hDivFlow(i,j,bi,bj)
94 & +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
95 & +vTrans(i,j+1)-vTrans(i,j) )
96 ENDDO
97 ENDDO
98
99 C- End DO k=1,Nr
100 ENDDO
101
102 ENDIF
103
104 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105
106 IF ( UpdateEtaN_EtaH ) THEN
107
108 C-- Update etaN at the end of the time step :
109 C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
110 IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
111 DO j=1-Oly,sNy+Oly
112 DO i=1-Olx,sNx+Olx
113 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
114 ENDDO
115 ENDDO
116 ELSEIF (useRealFreshWaterFlux) THEN
117 DO j=1,sNy
118 DO i=1,sNx
119 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
120 & - implicDiv2Dflow*( EmPmR(i,j,bi,bj)
121 & +hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
122 & )*deltaTfreesurf
123 ENDDO
124 ENDDO
125 ELSE
126 DO j=1,sNy
127 DO i=1,sNx
128 etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
129 & - implicDiv2Dflow*hDivFlow(i,j,bi,bj)
130 & *recip_rA(i,j,bi,bj)*deltaTfreesurf
131 ENDDO
132 ENDDO
133 ENDIF
134
135 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
136
137 ELSE
138
139 #ifdef NONLIN_FRSURF
140 IF (useRealFreshWaterFlux .AND. nonlinFreeSurf.GT.0) THEN
141
142 C-- Called at the beginning of the time step :
143 C- keep present time EmPmR to compute later (S/R EXTERNAL_FORCING_SURF)
144 C tracers and momentum flux associated with fresh water input.
145
146 IF ( myTime.NE.startTime ) THEN
147 DO j=1-Oly,sNy+Oly
148 DO i=1-Olx,sNx+Olx
149 PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
150 ENDDO
151 ENDDO
152
153 ELSEIF( myTime .EQ. 0. _d 0 ) THEN
154 DO j=1-Oly,sNy+Oly
155 DO i=1-Olx,sNx+Olx
156 PmEpR(i,j,bi,bj) = 0. _d 0
157 ENDDO
158 ENDDO
159
160 ELSE
161 C needs previous time-step value of E-P-R, that has not been loaded
162 C and was not in pickup-file ; try to use etaN & etaH instead.
163 DO j=1,sNy
164 DO i=1,sNx
165 PmEpR(i,j,bi,bj) =
166 & hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
167 & + (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
168 & /(implicDiv2Dflow*deltaTfreesurf)
169 ENDDO
170 ENDDO
171 ENDIF
172
173 ENDIF
174 #endif /* NONLIN_FRSURF */
175
176 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
177
178 C-- Update etaH at the beginning of the time step :
179 C Incorporate the Explicit part of -Divergence(Barotropic_Flow)
180
181 IF ( useRealFreshWaterFlux .AND. myTime.EQ.startTime ) THEN
182 C needs previous time-step value of E-P-R, that has not been loaded
183 C and was not in pickup-file ; try to use etaN & etaH instead.
184 DO j=1-Oly,sNy+Oly
185 DO i=1-Olx,sNx+Olx
186 etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
187 & + (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
188 & *(1. - implicDiv2Dflow)/implicDiv2Dflow
189 ENDDO
190 ENDDO
191
192 ELSEIF (implicDiv2Dflow.EQ. 1. _d 0) THEN
193 DO j=1-Oly,sNy+Oly
194 DO i=1-Olx,sNx+Olx
195 etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
196 ENDDO
197 ENDDO
198
199 ELSEIF (useRealFreshWaterFlux) THEN
200 DO j=1,sNy
201 DO i=1,sNx
202 etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
203 & - (1. - implicDiv2Dflow)*( EmPmR(i,j,bi,bj)
204 & +hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
205 & )*deltaTfreesurf
206 ENDDO
207 ENDDO
208
209 ELSE
210 DO j=1,sNy
211 DO i=1,sNx
212 etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
213 & - (1. - implicDiv2Dflow)*hDivFlow(i,j,bi,bj)
214 & *recip_rA(i,j,bi,bj)*deltaTfreesurf
215 ENDDO
216 ENDDO
217 ENDIF
218
219 #ifdef ALLOW_OBCS
220 #ifdef NONLIN_FRSURF
221 C- note: 1) needs to apply OBC to etaH since viscous terms depend on hFacZ.
222 C that is not only function of boundaries hFac values.
223 C 2) has to be done before calc_surf_dr; but since obcs_calc is
224 C called later, hFacZ will lag 1 time step behind OBC update.
225 C 3) avoid also unrealistic value of etaH in OB regions that
226 C might produce many "WARNING" message from calc_surf_dr.
227 C-------
228 IF ( useOBCS .AND. nonlinFreeSurf.GT.0 )
229 & CALL OBCS_APPLY_ETA( bi, bj, etaH, myThid )
230 #endif /* NONLIN_FRSURF */
231 #endif /* ALLOW_OBCS */
232
233 ENDIF
234
235 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
236
237 #endif /* EXACT_CONSERV */
238
239 RETURN
240 END
  ViewVC Help
Powered by ViewVC 1.1.22