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

Annotation of /MITgcm/model/src/calc_r_star.F

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


Revision 1.23 - (hide annotations) (download)
Fri Aug 15 08:21:57 2014 UTC (9 years, 10 months ago) by atn
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65c, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.22: +44 -47 lines
print out all points pertaining to r_star crash

1 atn 1.23 C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.22 2014/04/29 21:07:39 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4 jmc 1.4 #include "PACKAGES_CONFIG.h"
5 jmc 1.1 #include "CPP_OPTIONS.h"
6 jmc 1.7 #ifdef ALLOW_EXCH2
7     # include "W2_OPTIONS.h"
8     #endif
9 jmc 1.1
10     CBOP
11     C !ROUTINE: CALC_R_STAR
12     C !INTERFACE:
13     SUBROUTINE CALC_R_STAR( etaFld,
14     I myTime, myIter, myThid )
15     C !DESCRIPTION: \bv
16     C *==========================================================*
17 jmc 1.12 C | SUBROUTINE CALC_R_STAR
18 jmc 1.1 C | o Calculate new column thickness & scaling factor for r*
19 jmc 1.20 C | according to the surface r-position (Non-Lin Free-Surf)
20 jmc 1.1 C *==========================================================*
21     C \ev
22    
23     C !USES:
24     IMPLICIT NONE
25     C == Global variables
26     #include "SIZE.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29     #include "GRID.h"
30     #include "SURFACE.h"
31    
32     C !INPUT/OUTPUT PARAMETERS:
33     C == Routine arguments ==
34 jmc 1.21 C etaFld :: current eta field used to update the hFactor
35     C myTime :: current time in simulation
36     C myIter :: current iteration number in simulation
37     C myThid :: thread number for this instance of the routine.
38 jmc 1.22 _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
39 jmc 1.1 _RL myTime
40     INTEGER myIter
41     INTEGER myThid
42    
43     #ifdef NONLIN_FRSURF
44    
45     C !LOCAL VARIABLES:
46     C Local variables
47 jmc 1.17 C rStarAreaWeight :: use area weighted average for rStarFac at U & V point
48     C i,j, bi,bj :: loop counter
49     C numbWrite :: count the Number of warning written on STD-ERR file
50     C numbWrMax :: maximum Number of warning written on STD-ERR file
51     LOGICAL rStarAreaWeight
52 atn 1.23 _RL maxhFacC
53 jmc 1.14 INTEGER i,j,bi,bj
54 jmc 1.7 INTEGER numbWrite, numbWrMax
55 mlosch 1.11 INTEGER icntc1, icntc2, icntw, icnts
56 jmc 1.1 _RL tmpfldW, tmpfldS
57     CEOP
58 jmc 1.17
59 jmc 1.20 #ifdef W2_FILL_NULL_REGIONS
60     INTEGER ii,jj
61     #endif
62     DATA numbWrite / 0 /
63     numbWrMax = Nx*Ny
64    
65 atn 1.23 maxhFacC = 0. _d 0
66    
67 jmc 1.17 rStarAreaWeight = .TRUE.
68     C- Area-weighted average consistent with KE (& vert. advection):
69 jmc 1.18 IF ( vectorInvariantMomentum .AND.
70     & (selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3)
71     & ) rStarAreaWeight =.FALSE.
72 jmc 1.17
73 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
74    
75 jmc 1.4 #ifdef ALLOW_DEBUG
76     IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
77     #endif
78    
79 jmc 1.1 DO bj=myByLo(myThid), myByHi(myThid)
80     DO bi=myBxLo(myThid), myBxHi(myThid)
81 jmc 1.20
82 gforget 1.19 C-- before updating rStarFacC/S/W save current fields
83 jmc 1.22 DO j=1-OLy,sNy+OLy
84     DO i=1-OLx,sNx+OLx
85 gforget 1.19 rStarFacNm1C(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
86     rStarFacNm1S(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
87     rStarFacNm1W(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
88     ENDDO
89     ENDDO
90 jmc 1.20
91 jmc 1.1 C- 1rst bi,bj loop :
92    
93     C-- copy rStarFacX -> rStarExpX
94 jmc 1.22 DO j=1-OLy,sNy+OLy
95     DO i=1-OLx,sNx+OLx
96 jmc 1.1 rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
97     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
98     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
99     ENDDO
100 jmc 1.13 ENDDO
101 jmc 1.1
102     C-- Compute the new column thikness :
103     DO j=0,sNy+1
104     DO i=0,sNx+1
105 jmc 1.15 IF (kSurfC(i,j,bi,bj).LE.Nr ) THEN
106 jmc 1.12 rStarFacC(i,j,bi,bj) =
107 jmc 1.1 & (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
108     & *recip_Rcol(i,j,bi,bj)
109     ELSE
110     rStarFacC(i,j,bi,bj) = 1.
111     ENDIF
112     ENDDO
113     ENDDO
114 jmc 1.17 IF ( rStarAreaWeight ) THEN
115     C- Area weighted average
116 jmc 1.2 DO j=1,sNy
117 jmc 1.1 DO i=1,sNx+1
118 jmc 1.15 IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
119 jmc 1.16 tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
120 jmc 1.12 rStarFacW(i,j,bi,bj) =
121 jmc 1.1 & ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
122     & +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
123     & )*recip_rAw(i,j,bi,bj)
124     & +tmpfldW )/tmpfldW
125     ELSE
126     rStarFacW(i,j,bi,bj) = 1.
127     ENDIF
128 jmc 1.2 ENDDO
129     ENDDO
130     DO j=1,sNy+1
131     DO i=1,sNx
132 jmc 1.15 IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
133 jmc 1.16 tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
134 jmc 1.12 rStarFacS(i,j,bi,bj) =
135 jmc 1.1 & ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
136     & +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
137     & )*recip_rAs(i,j,bi,bj)
138     & +tmpfldS )/tmpfldS
139     ELSE
140     rStarFacS(i,j,bi,bj) = 1.
141 jmc 1.2 ENDIF
142     ENDDO
143     ENDDO
144 jmc 1.17 ELSE
145     C- Simple average
146     DO j=1,sNy
147     DO i=1,sNx+1
148     IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
149     tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
150     rStarFacW(i,j,bi,bj) =
151     & ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj) + etaFld(i,j,bi,bj) )
152     & +tmpfldW )/tmpfldW
153     ELSE
154     rStarFacW(i,j,bi,bj) = 1.
155     ENDIF
156     ENDDO
157     ENDDO
158     DO j=1,sNy+1
159     DO i=1,sNx
160     IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
161     tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
162     rStarFacS(i,j,bi,bj) =
163     & ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj) + etaFld(i,j,bi,bj) )
164     & +tmpfldS )/tmpfldS
165     ELSE
166     rStarFacS(i,j,bi,bj) = 1.
167     ENDIF
168     ENDDO
169     ENDDO
170     ENDIF
171 jmc 1.21 #ifdef ALLOW_OBCS
172     IF (useOBCS) THEN
173     CALL OBCS_APPLY_R_STAR(
174     I bi, bj, etaFld,
175     U rStarFacC, rStarFacW, rStarFacS,
176     I myTime, myIter, myThid )
177     ENDIF
178     #endif /* ALLOW_OBCS */
179 jmc 1.2
180     C- Needs to do something when r* ratio is too small ;
181     C for now, just stop
182 mlosch 1.11 icntc1 = 0
183     icntc2 = 0
184     icntw = 0
185     icnts = 0
186 jmc 1.2 DO j=1,sNy+1
187     DO i=1,sNx+1
188 jmc 1.12 IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
189 mlosch 1.11 icntc1 = icntc1 + 1
190 jmc 1.2 ENDIF
191 jmc 1.12 IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
192 mlosch 1.11 icntw = icntw + 1
193 jmc 1.2 ENDIF
194 jmc 1.12 IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
195 mlosch 1.11 icnts = icnts + 1
196 jmc 1.2 ENDIF
197 mlosch 1.11 IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
198     icntc2 = icntc2 + 1
199 atn 1.23 maxhFacC = max(rStarFacC(i,j,bi,bj),maxhFacC)
200 jmc 1.1 ENDIF
201     ENDDO
202     ENDDO
203 atn 1.23
204 jmc 1.12 IF ( icntc1+icnts+icntw .GT. 0 ) THEN
205     C- Print an error msg and then stop:
206 atn 1.23 DO j=1,sNy+1
207     DO i=1,sNx+1
208     IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
209     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
210     & ' fail at i,j=',i,j,' ; rStarFacC,H,eta =',
211     & rStarFacC(i,j,bi,bj),
212     & Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj),
213     & etaFld(i,j,bi,bj)
214     ENDIF
215     IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
216     tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
217     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
218     & ' fail at i,j=',i,j,' ; rStarFacW,H,eta =',
219     & rStarFacW(i,j,bi,bj), tmpfldW,
220     & etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
221     ENDIF
222     IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
223     tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
224     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
225     & ' fail at i,j=',i,j,' ; rStarFacS,H,eta =',
226     & rStarFacS(i,j,bi,bj), tmpfldS,
227     & etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
228     ENDIF
229     ENDDO
230     ENDDO
231     IF ( icntc1 .GT. 0 )
232     & WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
233     & 'WARNING: r*FacC < hFacInf at',icntc1,
234     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
235     IF ( icntw .GT. 0 )
236     & WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
237 jmc 1.12 & 'WARNING: r*FacW < hFacInf at',icntw,
238     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
239 atn 1.23 IF ( icnts .GT. 0 )
240     & WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
241 jmc 1.12 & 'WARNING: r*FacS < hFacInf at',icnts,
242     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
243 atn 1.23 WRITE(errorMessageUnit,'(A)')
244     & 'STOP in CALC_R_STAR : too SMALL rStarFac[C,W,S] !'
245 mlosch 1.11 STOP 'ABNORMAL END: S/R CALC_R_STAR'
246 jmc 1.12 ENDIF
247 atn 1.23
248 jmc 1.12 C-- Usefull warning when r*Fac becomes very large:
249     IF ( icntc2.GT.0 .AND. numbWrite.LE.numbWrMax ) THEN
250     numbWrite = numbWrite + 1
251     WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
252     & 'WARNING: r*FacC > hFacSup at',icntc2,
253     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
254 atn 1.23 WRITE(errorMessageUnit,'(A,E14.6)')
255     & 'WARNING: max(hFacC) is ',maxhFacC
256 jmc 1.12 ENDIF
257 jmc 1.1
258     C- end 1rst bi,bj loop.
259     ENDDO
260 jmc 1.3 ENDDO
261 jmc 1.1
262 jmc 1.12 _EXCH_XY_RL( rStarFacC, myThid )
263 jmc 1.1 CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
264    
265     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
266    
267     DO bj=myByLo(myThid), myByHi(myThid)
268     DO bi=myBxLo(myThid), myBxHi(myThid)
269     C- 2nd bi,bj loop :
270    
271 jmc 1.4 #ifdef ALLOW_EXCH2
272 jmc 1.7 #ifdef W2_FILL_NULL_REGIONS
273 jmc 1.4 C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
274 jmc 1.22 C in the corner regions of the tile (e.g.:[1-OLx:0,1-OLy:0])
275 jmc 1.4 C => need to add those lines (or to fix exch2):
276 jmc 1.22 DO j=1,OLy
277     DO i=1,OLx
278 jmc 1.4 ii = sNx+i
279     jj = sNy+j
280    
281 jmc 1.15 IF (kSurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
282     IF (kSurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
283     IF (kSurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
284     IF (kSurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.
285    
286     IF (kSurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
287     IF (kSurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
288     IF (kSurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
289     IF (kSurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.
290    
291     IF (kSurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
292     IF (kSurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
293     IF (kSurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
294     IF (kSurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.
295 jmc 1.7
296 jmc 1.4 ENDDO
297     ENDDO
298 jmc 1.7 #endif /* W2_FILL_NULL_REGIONS */
299 jmc 1.4 #endif /* ALLOW_EXCH2 */
300    
301 jmc 1.22 DO j=1-OLy,sNy+OLy
302     DO i=1-OLx,sNx+OLx
303 jmc 1.12 rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
304 jmc 1.22 & -rStarExpC(i,j,bi,bj))/deltaTFreeSurf
305 jmc 1.1 rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
306 jmc 1.22 & -rStarExpW(i,j,bi,bj))/deltaTFreeSurf
307 jmc 1.1 rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
308 jmc 1.22 & -rStarExpS(i,j,bi,bj))/deltaTFreeSurf
309 jmc 1.1 rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
310     & / rStarExpC(i,j,bi,bj)
311     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
312     & / rStarExpW(i,j,bi,bj)
313     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
314     & / rStarExpS(i,j,bi,bj)
315     ENDDO
316     ENDDO
317    
318 jmc 1.22 IF ( fluidIsAir ) THEN
319     DO j=1-OLy,sNy+OLy
320     DO i=1-OLx,sNx+OLx
321     pStarFacK(i,j,bi,bj) = rStarFacC(i,j,bi,bj)**atm_kappa
322     ENDDO
323     ENDDO
324     #ifdef ALLOW_AUTODIFF
325     ELSE
326     DO j=1-OLy,sNy+OLy
327     DO i=1-OLx,sNx+OLx
328     pStarFacK(i,j,bi,bj) = 1. _d 0
329     ENDDO
330     ENDDO
331     #endif
332     ENDIF
333    
334 jmc 1.1 C- end 2nd bi,bj loop.
335     ENDDO
336     ENDDO
337    
338 jmc 1.4 #ifdef ALLOW_DEBUG
339     IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
340     #endif
341    
342 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
343     #endif /* NONLIN_FRSURF */
344    
345     RETURN
346     END

  ViewVC Help
Powered by ViewVC 1.1.22