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

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

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

Revision 1.8 - (hide annotations) (download)
Fri Apr 11 13:00:21 2003 UTC (21 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint50c_post, checkpoint52d_pre, checkpoint50c_pre, checkpoint52j_pre, checkpoint51o_pre, checkpoint51l_post, checkpoint52l_post, checkpoint52k_post, checkpoint54, checkpoint51, checkpoint53, checkpoint52, checkpoint50d_post, checkpoint52f_post, checkpoint50b_pre, checkpoint51f_post, checkpoint51d_post, checkpoint51t_post, checkpoint51n_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint51j_post, checkpoint52e_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint53d_post, checkpoint52b_pre, checkpoint51l_pre, checkpoint52m_post, checkpoint51q_post, checkpoint51b_pre, checkpoint52b_post, checkpoint52c_post, checkpoint51h_pre, checkpoint50f_post, checkpoint50f_pre, checkpoint52f_pre, checkpoint54a_pre, checkpoint53c_post, branchpoint-genmake2, checkpoint51r_post, checkpoint51i_post, checkpoint51b_post, checkpoint51c_post, checkpoint53a_post, checkpoint52d_post, checkpoint53g_post, checkpoint50g_post, checkpoint52a_pre, checkpoint50h_post, checkpoint52i_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre, checkpoint52h_pre, checkpoint53f_post, checkpoint52j_post, checkpoint50e_post, branch-netcdf, checkpoint50d_pre, checkpoint52n_post, checkpoint53b_pre, checkpoint51e_post, checkpoint51o_post, checkpoint51f_pre, checkpoint53b_post, checkpoint52a_post, checkpoint51g_post, ecco_c52_e35, checkpoint50b_post, checkpoint51m_post, checkpoint53d_pre, checkpoint51a_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-genmake2, branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.7: +5 -7 lines 
limit number of WARNING when hFac becomes very large (> hFacSup);
1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/model/src/calc_surf_dr.F,v 1.7 2002/08/09 22:55:35 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6 cnh 1.3 CBOP
7     C !ROUTINE: CALC_SURF_DR
8     C !INTERFACE:
9 jmc 1.1 SUBROUTINE CALC_SURF_DR(etaFld,
10     I myTime, myIter, myThid )
11 cnh 1.3 C !DESCRIPTION: \bv
12     C *==========================================================*
13     C | SUBROUTINE CALC_SURF_DR
14     C | o Calculate the new surface level thickness according to
15     C | the surface r-position (Non-Linear Free-Surf)
16     C | o take decision if grid box becomes too thin or too thick
17     C *==========================================================*
18     C \ev
19    
20     C !USES:
21 jmc 1.1 IMPLICIT NONE
22     C == Global variables
23     #include "SIZE.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "GRID.h"
27     #include "SURFACE.h"
28    
29 cnh 1.3 C !INPUT/OUTPUT PARAMETERS:
30 jmc 1.1 C == Routine arguments ==
31 cnh 1.3 C myTime :: Current time in simulation
32     C myIter :: Current iteration number in simulation
33     C myThid :: Thread number for this instance of the routine.
34     C etaFld :: current eta field used to update the hFactor
35 jmc 1.1 _RL myTime
36     INTEGER myIter
37     INTEGER myThid
38     _RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
39    
40     #ifdef NONLIN_FRSURF
41    
42 cnh 1.3 C !LOCAL VARIABLES:
43 jmc 1.2 C Local variables in common block
44 cnh 1.3 C Rmin_surf :: minimum r_value of the free surface position
45 jmc 1.2 C that satisfy the hFacInf criteria
46     COMMON /LOCAL_CALC_SURF_DR/ Rmin_surf
47     _RL Rmin_surf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48 jmc 1.1 C Local variables
49 cnh 1.3 C i,j,k,bi,bj :: loop counter
50     C rSurftmp :: free surface r-position that is used to compute hFac_surf
51     C adjust_nb_pt :: Nb of grid points where rSurf is adjusted (hFactInf)
52     C adjust_volum :: adjustment effect on the volume (domain size)
53 jmc 1.7 C numbWrite :: count the Number of warning written on STD-ERR file
54     C numbWrMax :: maximum Number of warning written on STD-ERR file
55 jmc 1.6 INTEGER i,j,bi,bj
56 jmc 1.7 INTEGER ks, numbWrite, numbWrMax
57 jmc 1.2 _RL hFacInfMOM, Rmin_tmp, hFactmp, adjust_nb_pt, adjust_volum
58     _RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59 jmc 1.1 _RS hhm, hhp
60     CHARACTER*(MAX_LEN_MBUF) suff
61 cnh 1.3 CEOP
62 jmc 1.7 DATA numbWrite / 0 /
63     numbWrMax = Nx*Ny
64 jmc 1.1
65     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
66    
67 jmc 1.4 IF (myIter.EQ.-1) THEN
68 jmc 1.2
69     hFacInfMOM = hFacInf
70    
71 jmc 1.1 DO bj=myByLo(myThid), myByHi(myThid)
72     DO bi=myBxLo(myThid), myBxHi(myThid)
73 jmc 1.2
74     C-- Initialise arrays :
75 jmc 1.1 DO j=1-Oly,sNy+Oly
76     DO i=1-Olx,sNx+Olx
77     hFac_surfC(i,j,bi,bj) = 0.
78     hFac_surfW(i,j,bi,bj) = 0.
79     hFac_surfS(i,j,bi,bj) = 0.
80 jmc 1.6 PmEpR(i,j,bi,bj) = 0.
81 jmc 1.2 Rmin_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
82 jmc 1.1 ENDDO
83     ENDDO
84 jmc 1.2
85     C-- Compute the mimimum value of r_surf (used for computing hFac_surfC)
86     DO j=1,sNy
87     DO i=1,sNx
88     ks = ksurfC(i,j,bi,bj)
89     IF (ks.LE.Nr) THEN
90     Rmin_tmp = rF(ks+1)
91     IF ( ks.EQ.ksurfW(i,j,bi,bj))
92     & Rmin_tmp = MAX(Rmin_tmp, R_low(i-1,j,bi,bj))
93     IF ( ks.EQ.ksurfW(i+1,j,bi,bj))
94     & Rmin_tmp = MAX(Rmin_tmp, R_low(i+1,j,bi,bj))
95     IF ( ks.EQ.ksurfS(i,j,bi,bj))
96     & Rmin_tmp = MAX(Rmin_tmp, R_low(i,j-1,bi,bj))
97     IF ( ks.EQ.ksurfS(i,j+1,bi,bj))
98     & Rmin_tmp = MAX(Rmin_tmp, R_low(i,j+1,bi,bj))
99    
100     Rmin_surf(i,j,bi,bj) =
101     & MAX( MAX(rF(ks+1),R_low(i,j,bi,bj)) + hFacInf*drF(ks),
102     & Rmin_tmp + hFacInfMOM*drF(ks)
103     & )
104     ENDIF
105     ENDDO
106     ENDDO
107    
108     C- end bi,bj loop.
109 jmc 1.1 ENDDO
110     ENDDO
111 jmc 1.2
112     _EXCH_XY_R8( Rmin_surf, myThid )
113    
114     C- end of initialization block
115 jmc 1.1 ENDIF
116    
117     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
118    
119 jmc 1.2 adjust_nb_pt = 0.
120     adjust_volum = 0.
121    
122 jmc 1.1 DO bj=myByLo(myThid), myByHi(myThid)
123     DO bi=myBxLo(myThid), myBxHi(myThid)
124 jmc 1.2
125     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
126     C-- Compute the new fractional thickness of surface level (ksurfC):
127    
128     DO j=0,sNy+1
129     DO i=0,sNx+1
130     rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
131 jmc 1.1 ks = ksurfC(i,j,bi,bj)
132     IF (ks.LE.Nr) THEN
133 jmc 1.2 IF (rSurftmp(i,j) .LT. Rmin_surf(i,j,bi,bj)) THEN
134     C-- Needs to do something :
135 jmc 1.7 IF (numbWrite.LE.numbWrMax) THEN
136     numbWrite = numbWrite + 1
137 jmc 1.2 hFactmp = ( rSurftmp(i,j)-MAX(rF(ks+1),R_low(i,j,bi,bj))
138     & )*recip_drF(ks)
139 jmc 1.1 IF (hFactmp.LT.hFacInf) THEN
140 jmc 1.7 WRITE(errorMessageUnit,'(2A,6I4,I10)')
141     & 'WARNING: hFacC < hFacInf at:',
142 jmc 1.1 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
143 jmc 1.2 ELSE
144 jmc 1.7 WRITE(errorMessageUnit,'(2A,6I4,I10)')
145     & 'WARNING: hFac < hFacInf near:',
146 jmc 1.2 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
147 jmc 1.1 ENDIF
148 jmc 1.7 WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
149     & 'hFac_n-1,hFac_n,eta =',
150 jmc 1.2 & hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
151 jmc 1.7 ENDIF
152 jmc 1.2 C-- Decide to STOP :
153 jmc 1.7 c WRITE(errorMessageUnit,'(A)')
154     c & 'STOP in CALC_SURF_DR : too SMALL hFac !'
155 jmc 1.2 c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
156     C----------
157    
158     C-- Continue with Rmin_surf:
159     IF ( i.GE.1.AND.i.LE.sNx .AND.
160     & j.GE.1.AND.j.LE.sNy ) THEN
161     adjust_nb_pt = adjust_nb_pt + 1.
162     adjust_volum = adjust_volum
163     & + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
164     ENDIF
165     rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
166     C----------
167     ENDIF
168    
169     C-- Set hFac_surfC :
170     hFac_surfC(i,j,bi,bj) =
171     & ( rSurftmp(i,j) - MAX(rF(ks+1), R_low(i,j,bi,bj))
172     & )*recip_drF(ks)*maskC(i,j,ks,bi,bj)
173    
174     C-- Usefull warning when hFac becomes very large:
175 jmc 1.8 IF ( numbWrite.LE.numbWrMax .AND.
176     & hFac_surfC(i,j,bi,bj).GT.hFacSup ) THEN
177     numbWrite = numbWrite + 1
178 jmc 1.7 WRITE(errorMessageUnit,'(2A,6I4,I10)')
179     & 'WARNING: hFacC > hFacSup at:',
180 jmc 1.1 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
181 jmc 1.7 WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
182     & 'hFac_n-1,hFac_n,eta =', hfacC(i,j,ks,bi,bj),
183     & hFac_surfC(i,j,bi,bj), etaFld(i,j,bi,bj)
184 jmc 1.8 ENDIF
185 jmc 1.1 C----------
186     ENDIF
187 jmc 1.2
188 jmc 1.1 ENDDO
189     ENDDO
190    
191     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192     C-- Compute fractional thickness of surface level, for U & V point:
193    
194     DO j=1,sNy
195     DO i=1,sNx+1
196     ks = ksurfW(i,j,bi,bj)
197     IF (ks.LE.Nr) THEN
198 jmc 1.2 hhm = rF(ks)
199     IF(ks.EQ.ksurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
200     hhp = rF(ks)
201     IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
202     hFac_surfW(i,j,bi,bj) =
203     & ( MIN(hhm,hhp)
204     & - MAX(rF(ks+1),R_low(i-1,j,bi,bj),R_low(i,j,bi,bj))
205     & )*recip_drF(ks)*maskW(i,j,ks,bi,bj)
206 jmc 1.1 ENDIF
207     ENDDO
208     ENDDO
209 jmc 1.2
210 jmc 1.1 DO j=1,sNy+1
211     DO i=1,sNx
212     ks = ksurfS(i,j,bi,bj)
213     IF (ks.LE.Nr) THEN
214 jmc 1.2 hhm = rF(ks)
215     IF(ks.EQ.ksurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
216     hhp = rF(ks)
217     IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
218     hFac_surfS(i,j,bi,bj) =
219     & ( MIN(hhm,hhp)
220     & - MAX(rF(ks+1),R_low(i,j-1,bi,bj),R_low(i,j,bi,bj))
221     & )*recip_drF(ks)*maskS(i,j,ks,bi,bj)
222 jmc 1.1 ENDIF
223     ENDDO
224     ENDDO
225 jmc 1.2
226     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
227    
228     C- end bi,bj loop.
229 jmc 1.1 ENDDO
230     ENDDO
231    
232 jmc 1.2 C-- Global diagnostic :
233     _GLOBAL_SUM_R8( adjust_nb_pt , myThid )
234     _GLOBAL_SUM_R8( adjust_volum , myThid )
235     IF (adjust_nb_pt .GE.1.) THEN
236     _BEGIN_MASTER( myThid )
237 jmc 1.7 WRITE(standardMessageUnit,'(2(A,I10),1PE16.8)')
238     & ' SURF_ADJUSTMENT: Iter=', myIter,
239     & ' Nb_pts,Vol=', nint(adjust_nb_pt), adjust_volum
240     _END_MASTER( myThid )
241 jmc 1.2 ENDIF
242    
243     _EXCH_XY_R4(hFac_surfC, myThid )
244 jmc 1.1 CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)
245 jmc 1.5
246     IF (useRealFreshWaterFlux .AND. myTime.EQ.startTime)
247     & _EXCH_XY_R4( PmEpR, myThid )
248 jmc 1.1
249 jmc 1.2 C-----
250     C Note: testing ksurfW,S is equivalent to a full height mask
251     C ==> no need for applying the mask here.
252     C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
253     C-----
254 jmc 1.1
255     WRITE(suff,'(I10.10)') myIter
256     c CALL WRITE_FLD_XY_RS('hFac_surfC.',suff,hFac_surfC,
257     c & myIter,myThid)
258    
259     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
260     #endif /* NONLIN_FRSURF */
261    
262     RETURN
263     END
  ViewVC Help
Powered by ViewVC 1.1.22