/[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.21 - (hide annotations) (download)
Fri May 20 01:09:37 2011 UTC (13 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint62y, checkpoint63g, checkpoint64, checkpoint63, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.20: +13 -5 lines
move call to OBCS_APPLY_R_STAR from update_r_star.F to calc_r_star.F
 and add argument etaFld:
this fixes missing EXCH + get consistent eta OB value (+ fix restart)

1 jmc 1.21 C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.20 2011/02/21 16:47:47 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.20 _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 jmc 1.14 INTEGER i,j,bi,bj
53 jmc 1.7 INTEGER numbWrite, numbWrMax
54 mlosch 1.11 INTEGER icntc1, icntc2, icntw, icnts
55     INTEGER iic, jjc, iic2, jjc2, iiw, jjw, iis, jjs
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 jmc 1.17 rStarAreaWeight = .TRUE.
66     C- Area-weighted average consistent with KE (& vert. advection):
67 jmc 1.18 IF ( vectorInvariantMomentum .AND.
68     & (selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3)
69     & ) rStarAreaWeight =.FALSE.
70 jmc 1.17
71 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
72    
73 jmc 1.4 #ifdef ALLOW_DEBUG
74     IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
75     #endif
76    
77 jmc 1.1 DO bj=myByLo(myThid), myByHi(myThid)
78     DO bi=myBxLo(myThid), myBxHi(myThid)
79 jmc 1.20
80 gforget 1.19 C-- before updating rStarFacC/S/W save current fields
81     DO j=1-Oly,sNy+Oly
82     DO i=1-Olx,sNx+Olx
83     rStarFacNm1C(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
84     rStarFacNm1S(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
85     rStarFacNm1W(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
86     ENDDO
87     ENDDO
88 jmc 1.20
89 jmc 1.1 C- 1rst bi,bj loop :
90    
91     C-- copy rStarFacX -> rStarExpX
92 jmc 1.13 DO j=1-Oly,sNy+Oly
93 jmc 1.12 DO i=1-Olx,sNx+Olx
94 jmc 1.1 rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
95     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
96     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
97     ENDDO
98 jmc 1.13 ENDDO
99 jmc 1.1
100     C-- Compute the new column thikness :
101     DO j=0,sNy+1
102     DO i=0,sNx+1
103 jmc 1.15 IF (kSurfC(i,j,bi,bj).LE.Nr ) THEN
104 jmc 1.12 rStarFacC(i,j,bi,bj) =
105 jmc 1.1 & (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
106     & *recip_Rcol(i,j,bi,bj)
107     ELSE
108     rStarFacC(i,j,bi,bj) = 1.
109     ENDIF
110     ENDDO
111     ENDDO
112 jmc 1.17 IF ( rStarAreaWeight ) THEN
113     C- Area weighted average
114 jmc 1.2 DO j=1,sNy
115 jmc 1.1 DO i=1,sNx+1
116 jmc 1.15 IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
117 jmc 1.16 tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
118 jmc 1.12 rStarFacW(i,j,bi,bj) =
119 jmc 1.1 & ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
120     & +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
121     & )*recip_rAw(i,j,bi,bj)
122     & +tmpfldW )/tmpfldW
123     ELSE
124     rStarFacW(i,j,bi,bj) = 1.
125     ENDIF
126 jmc 1.2 ENDDO
127     ENDDO
128     DO j=1,sNy+1
129     DO i=1,sNx
130 jmc 1.15 IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
131 jmc 1.16 tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
132 jmc 1.12 rStarFacS(i,j,bi,bj) =
133 jmc 1.1 & ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
134     & +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
135     & )*recip_rAs(i,j,bi,bj)
136     & +tmpfldS )/tmpfldS
137     ELSE
138     rStarFacS(i,j,bi,bj) = 1.
139 jmc 1.2 ENDIF
140     ENDDO
141     ENDDO
142 jmc 1.17 ELSE
143     C- Simple average
144     DO j=1,sNy
145     DO i=1,sNx+1
146     IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
147     tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
148     rStarFacW(i,j,bi,bj) =
149     & ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj) + etaFld(i,j,bi,bj) )
150     & +tmpfldW )/tmpfldW
151     ELSE
152     rStarFacW(i,j,bi,bj) = 1.
153     ENDIF
154     ENDDO
155     ENDDO
156     DO j=1,sNy+1
157     DO i=1,sNx
158     IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
159     tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
160     rStarFacS(i,j,bi,bj) =
161     & ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj) + etaFld(i,j,bi,bj) )
162     & +tmpfldS )/tmpfldS
163     ELSE
164     rStarFacS(i,j,bi,bj) = 1.
165     ENDIF
166     ENDDO
167     ENDDO
168     ENDIF
169 jmc 1.21 #ifdef ALLOW_OBCS
170     IF (useOBCS) THEN
171     CALL OBCS_APPLY_R_STAR(
172     I bi, bj, etaFld,
173     U rStarFacC, rStarFacW, rStarFacS,
174     I myTime, myIter, myThid )
175     ENDIF
176     #endif /* ALLOW_OBCS */
177 jmc 1.2
178     C- Needs to do something when r* ratio is too small ;
179     C for now, just stop
180 mlosch 1.11 icntc1 = 0
181     icntc2 = 0
182     icntw = 0
183     icnts = 0
184 jmc 1.2 DO j=1,sNy+1
185     DO i=1,sNx+1
186 jmc 1.12 IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
187 mlosch 1.11 icntc1 = icntc1 + 1
188     iic = i
189     jjc = j
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     iiw = i
194     jjw = j
195 jmc 1.2 ENDIF
196 jmc 1.12 IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
197 mlosch 1.11 icnts = icnts + 1
198     iis = i
199     jjs = j
200 jmc 1.2 ENDIF
201 mlosch 1.11 IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
202     icntc2 = icntc2 + 1
203     iic2 = i
204     jjc2 = j
205 jmc 1.1 ENDIF
206     ENDDO
207     ENDDO
208 jmc 1.12 IF ( icntc1+icnts+icntw .GT. 0 ) THEN
209     C- Print an error msg and then stop:
210     IF ( icntw .GT. 0 ) THEN
211 jmc 1.16 tmpfldW = rSurfW(iiw,jjw,bi,bj) - rLowW(iiw,jjw,bi,bj)
212 jmc 1.12 WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
213     & 'WARNING: r*FacW < hFacInf at',icntw,
214     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
215     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
216     & ' e.g. at i,j=',iiw,jjw,' ; rStarFac,H,eta =',
217     & rStarFacW(iiw,jjw,bi,bj), tmpfldW,
218     & etaFld(iiw-1,jjw,bi,bj), etaFld(iiw,jjw,bi,bj)
219     WRITE(errorMessageUnit,'(A)')
220     & 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
221     ENDIF
222     IF ( icnts .GT. 0 ) THEN
223 jmc 1.16 tmpfldS = rSurfS(iis,jjs,bi,bj) - rLowS(iis,jjs,bi,bj)
224 jmc 1.12 WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
225     & 'WARNING: r*FacS < hFacInf at',icnts,
226     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
227     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
228     & ' e.g. at i,j=',iis,jjs,' ; rStarFac,H,eta =',
229     & rStarFacS(iis,jjs,bi,bj), tmpfldS,
230     & etaFld(iis,jjs-1,bi,bj), etaFld(iis,jjs,bi,bj)
231     WRITE(errorMessageUnit,'(A)')
232     & 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
233     ENDIF
234     IF ( icntc1 .GT. 0 ) THEN
235     WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
236     & 'WARNING: r*FacC < hFacInf at',icntc1,
237     & ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
238     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
239     & ' e.g. at i,j=',iic,jjc,' ; rStarFac,H,eta =',
240     & rStarFacC(iic,jjc,bi,bj),
241     & Ro_surf(iic,jjc,bi,bj)-R_low(iic,jjc,bi,bj),
242     & etaFld(iic,jjc,bi,bj)
243     WRITE(errorMessageUnit,'(A)')
244     & 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
245     ENDIF
246 mlosch 1.11 STOP 'ABNORMAL END: S/R CALC_R_STAR'
247 jmc 1.12 ENDIF
248     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     WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
255     & ' e.g. at i,j=',iic2,jjc2,' ; rStarFac,H,eta =',
256     & rStarFacC(iic2,jjc2,bi,bj),
257     & Ro_surf(iic2,jjc2,bi,bj)-R_low(iic2,jjc2,bi,bj),
258 mlosch 1.11 & etaFld(iic2,jjc2,bi,bj)
259 jmc 1.12 ENDIF
260 jmc 1.1
261     C- end 1rst bi,bj loop.
262     ENDDO
263 jmc 1.3 ENDDO
264 jmc 1.1
265 jmc 1.12 _EXCH_XY_RL( rStarFacC, myThid )
266 jmc 1.1 CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
267    
268     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
269    
270     DO bj=myByLo(myThid), myByHi(myThid)
271     DO bi=myBxLo(myThid), myBxHi(myThid)
272     C- 2nd bi,bj loop :
273    
274 jmc 1.4 #ifdef ALLOW_EXCH2
275 jmc 1.7 #ifdef W2_FILL_NULL_REGIONS
276 jmc 1.4 C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
277     C in the corner regions of the tile (e.g.:[1-Olx:0,1-Oly:0])
278     C => need to add those lines (or to fix exch2):
279     DO j=1,Oly
280     DO i=1,Olx
281     ii = sNx+i
282     jj = sNy+j
283    
284 jmc 1.15 IF (kSurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
285     IF (kSurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
286     IF (kSurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
287     IF (kSurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.
288    
289     IF (kSurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
290     IF (kSurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
291     IF (kSurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
292     IF (kSurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.
293    
294     IF (kSurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
295     IF (kSurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
296     IF (kSurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
297     IF (kSurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.
298 jmc 1.7
299 jmc 1.4 ENDDO
300     ENDDO
301 jmc 1.7 #endif /* W2_FILL_NULL_REGIONS */
302 jmc 1.4 #endif /* ALLOW_EXCH2 */
303    
304 jmc 1.1 DO j=1-Oly,sNy+Oly
305 jmc 1.12 DO i=1-Olx,sNx+Olx
306     rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
307 jmc 1.1 & -rStarExpC(i,j,bi,bj))/deltaTfreesurf
308     rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
309     & -rStarExpW(i,j,bi,bj))/deltaTfreesurf
310     rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
311     & -rStarExpS(i,j,bi,bj))/deltaTfreesurf
312     rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
313     & / rStarExpC(i,j,bi,bj)
314     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
315     & / rStarExpW(i,j,bi,bj)
316     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
317     & / rStarExpS(i,j,bi,bj)
318     ENDDO
319     ENDDO
320    
321     C- end 2nd bi,bj loop.
322     ENDDO
323     ENDDO
324    
325 jmc 1.4 #ifdef ALLOW_DEBUG
326     IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
327     #endif
328    
329 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
330     #endif /* NONLIN_FRSURF */
331    
332     RETURN
333     END

  ViewVC Help
Powered by ViewVC 1.1.22