/[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.3 - (hide annotations) (download)
Tue Jul 6 01:05:53 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57g_pre, checkpoint57b_post, checkpoint57g_post, checkpoint56b_post, checkpoint54d_post, checkpoint54e_post, checkpoint57d_post, checkpoint57i_post, checkpoint55, checkpoint57, checkpoint56, checkpoint54f_post, checkpoint55i_post, checkpoint55c_post, checkpoint57f_post, checkpoint57a_post, checkpoint57h_pre, checkpoint54b_post, checkpoint57h_post, checkpoint55g_post, checkpoint57c_post, checkpoint55d_post, checkpoint55d_pre, checkpoint57c_pre, checkpoint55j_post, checkpoint54a_post, checkpoint55h_post, checkpoint57e_post, checkpoint55b_post, checkpoint55f_post, eckpoint57e_pre, checkpoint56a_post, checkpoint57h_done, checkpoint57f_pre, checkpoint56c_post, checkpoint57a_pre, checkpoint55a_post, checkpoint55e_post, checkpoint54c_post
Changes since 1.2: +2 -5 lines 
re-write staggerTimeStep: step forward momentum 1rst and then T,S
1 jmc 1.3 C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.2 2003/04/11 13:02:37 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: CALC_R_STAR
8     C !INTERFACE:
9     SUBROUTINE CALC_R_STAR( etaFld,
10     I myTime, myIter, myThid )
11     C !DESCRIPTION: \bv
12     C *==========================================================*
13     C | SUBROUTINE CALC_R_STAR
14     C | o Calculate new column thickness & scaling factor for r*
15     C | according to the surface r-position (Non-Lin Free-Surf)
16     C *==========================================================*
17     C \ev
18    
19     C !USES:
20     IMPLICIT NONE
21     C == Global variables
22     #include "SIZE.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "GRID.h"
26     #include "SURFACE.h"
27    
28     C !INPUT/OUTPUT PARAMETERS:
29     C == Routine arguments ==
30     C myTime :: Current time in simulation
31     C myIter :: Current iteration number in simulation
32     C myThid :: Thread number for this instance of the routine.
33     C etaFld :: current eta field used to update the hFactor
34     _RL myTime
35     INTEGER myIter
36     INTEGER myThid
37     _RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
38    
39     #ifdef NONLIN_FRSURF
40    
41     C !LOCAL VARIABLES:
42     C Local variables
43     C i,j,k,bi,bj :: loop counter
44 jmc 1.2 C numbWrite :: count the Number of warning written on STD-ERR file
45     C numbWrMax :: maximum Number of warning written on STD-ERR file
46 jmc 1.1 INTEGER i,j,k,bi,bj
47 jmc 1.2 INTEGER km, numbWrite, numbWrMax
48 jmc 1.1 _RL tmpfldW, tmpfldS
49     c CHARACTER*(MAX_LEN_MBUF) suff
50     CEOP
51 jmc 1.2 DATA numbWrite / 0 /
52     numbWrMax = Nx*Ny
53 jmc 1.1
54     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
55    
56     IF (groundAtK1) THEN
57     km = 1
58     ELSE
59     km = Nr
60     ENDIF
61    
62     DO bj=myByLo(myThid), myByHi(myThid)
63     DO bi=myBxLo(myThid), myBxHi(myThid)
64     C- 1rst bi,bj loop :
65    
66     IF (myIter.EQ.-1) THEN
67     C-- Initialise arrays :
68     DO j=1-Oly,sNy+Oly
69     DO i=1-Olx,sNx+Olx
70     rStarFacC(i,j,bi,bj) = 1.
71     rStarFacW(i,j,bi,bj) = 1.
72     rStarFacS(i,j,bi,bj) = 1.
73     rStarExpC(i,j,bi,bj) = 1.
74     rStarExpW(i,j,bi,bj) = 1.
75     rStarExpS(i,j,bi,bj) = 1.
76     rStarDhCDt(i,j,bi,bj) = 0.
77     rStarDhWDt(i,j,bi,bj) = 0.
78     rStarDhSDt(i,j,bi,bj) = 0.
79     PmEpR(i,j,bi,bj) = 0.
80     ENDDO
81     ENDDO
82     DO k=1,Nr
83     DO j=1-Oly,sNy+Oly
84     DO i=1-Olx,sNx+Olx
85     h0FacC(i,j,k,bi,bj) = hFacC(i,j,k,bi,bj)
86     h0FacW(i,j,k,bi,bj) = hFacW(i,j,k,bi,bj)
87     h0FacS(i,j,k,bi,bj) = hFacS(i,j,k,bi,bj)
88     ENDDO
89     ENDDO
90     ENDDO
91     ELSE
92     C-- copy rStarFacX -> rStarExpX
93     DO j=1-Oly,sNy+Oly
94     DO i=1-Olx,sNx+Olx
95     rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
96     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
97     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
98     ENDDO
99     ENDDO
100     ENDIF
101    
102     C-- Compute the new column thikness :
103     DO j=0,sNy+1
104     DO i=0,sNx+1
105     IF (maskH(i,j,bi,bj).EQ.1. ) THEN
106     rStarFacC(i,j,bi,bj) =
107     & (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.2 DO j=1,sNy
115 jmc 1.1 DO i=1,sNx+1
116     IF (maskW(i,j,km,bi,bj).EQ.1. ) THEN
117     tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
118     & - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
119     rStarFacW(i,j,bi,bj) =
120     & ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
121     & +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
122     & )*recip_rAw(i,j,bi,bj)
123     & +tmpfldW )/tmpfldW
124     ELSE
125     rStarFacW(i,j,bi,bj) = 1.
126     ENDIF
127 jmc 1.2 ENDDO
128     ENDDO
129     DO j=1,sNy+1
130     DO i=1,sNx
131 jmc 1.1 IF (maskS(i,j,km,bi,bj).EQ.1. ) THEN
132     tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
133     & - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
134     rStarFacS(i,j,bi,bj) =
135     & ( 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    
145     C- Needs to do something when r* ratio is too small ;
146     C for now, just stop
147     DO j=1,sNy+1
148     DO i=1,sNx+1
149     IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
150     numbWrite = numbWrite + 1
151     WRITE(errorMessageUnit,'(2A,5I4,I10)')
152     & 'WARNING: r*FacC < hFacInf at:',
153     & ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
154     WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
155     & 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
156     & Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
157     WRITE(errorMessageUnit,'(A)')
158     & 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
159     STOP 'ABNORMAL END: S/R CALC_SURF_DR'
160     ENDIF
161     IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
162     numbWrite = numbWrite + 1
163     tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
164     & - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
165     WRITE(errorMessageUnit,'(2A,5I4,I10)')
166     & 'WARNING: r*FacW < hFacInf at:',
167     & ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
168     WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
169     & 'rStarFac,H,eta =', rStarFacW(i,j,bi,bj), tmpfldW,
170     & etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
171     WRITE(errorMessageUnit,'(A)')
172     & 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
173     STOP 'ABNORMAL END: S/R CALC_SURF_DR'
174     ENDIF
175     IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
176     numbWrite = numbWrite + 1
177     tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
178     & - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
179     WRITE(errorMessageUnit,'(2A,5I4,I10)')
180     & 'WARNING: r*FacS < hFacInf at:',
181     & ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
182     WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
183     & 'rStarFac,H,eta =', rStarFacS(i,j,bi,bj), tmpfldS,
184     & etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
185     WRITE(errorMessageUnit,'(A)')
186     & 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
187     STOP 'ABNORMAL END: S/R CALC_R_STAR'
188     ENDIF
189     C-- Usefull warning when r*Fac becomes very large:
190     IF ( numbWrite.LE.numbWrMax .AND.
191     & rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
192     numbWrite = numbWrite + 1
193     WRITE(errorMessageUnit,'(2A,5I4,I10)')
194     & 'WARNING: hFacC > hFacSup at:',
195     & ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
196     WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
197     & 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
198     & Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
199 jmc 1.1 ENDIF
200     ENDDO
201     ENDDO
202    
203     C- end 1rst bi,bj loop.
204     ENDDO
205 jmc 1.3 ENDDO
206 jmc 1.1
207     _EXCH_XY_RL( rStarFacC, myThid )
208     CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
209    
210     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211    
212     DO bj=myByLo(myThid), myByHi(myThid)
213     DO bi=myBxLo(myThid), myBxHi(myThid)
214     C- 2nd bi,bj loop :
215    
216     DO j=1-Oly,sNy+Oly
217     DO i=1-Olx,sNx+Olx
218     rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
219     & -rStarExpC(i,j,bi,bj))/deltaTfreesurf
220     rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
221     & -rStarExpW(i,j,bi,bj))/deltaTfreesurf
222     rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
223     & -rStarExpS(i,j,bi,bj))/deltaTfreesurf
224     rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
225     & / rStarExpC(i,j,bi,bj)
226     rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
227     & / rStarExpW(i,j,bi,bj)
228     rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
229     & / rStarExpS(i,j,bi,bj)
230     ENDDO
231     ENDDO
232    
233     C- end 2nd bi,bj loop.
234     ENDDO
235     ENDDO
236    
237     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
238     #endif /* NONLIN_FRSURF */
239    
240     RETURN
241     END
  ViewVC Help
Powered by ViewVC 1.1.22