/[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.14 - (hide annotations) (download)
Tue Apr 28 18:01:14 2009 UTC (15 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62q, checkpoint62p, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62, checkpoint62b, checkpoint61n, checkpoint61q, checkpoint61o, checkpoint61m, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.13: +4 -4 lines
change macros (EXCH & GLOBAL_SUM/MAX) sufix _R4/_R8 to _RS/_RL
 when applied to _RS/_RL variable

1 jmc 1.14 C $Header: /u/gcmpack/MITgcm/model/src/calc_surf_dr.F,v 1.13 2007/09/04 16:49:44 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 jmc 1.11 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 cnh 1.3 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.1 C Local variables
44 cnh 1.3 C i,j,k,bi,bj :: loop counter
45     C rSurftmp :: free surface r-position that is used to compute hFac_surf
46     C adjust_nb_pt :: Nb of grid points where rSurf is adjusted (hFactInf)
47     C adjust_volum :: adjustment effect on the volume (domain size)
48 jmc 1.7 C numbWrite :: count the Number of warning written on STD-ERR file
49     C numbWrMax :: maximum Number of warning written on STD-ERR file
50 jmc 1.6 INTEGER i,j,bi,bj
51 jmc 1.7 INTEGER ks, numbWrite, numbWrMax
52 jmc 1.11 _RL hFactmp, adjust_nb_pt, adjust_volum
53 jmc 1.2 _RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54 jmc 1.1 _RS hhm, hhp
55 jmc 1.13 c CHARACTER*(MAX_LEN_MBUF) suff
56 cnh 1.3 CEOP
57 jmc 1.7 DATA numbWrite / 0 /
58     numbWrMax = Nx*Ny
59 jmc 1.1
60     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
61    
62 jmc 1.2 adjust_nb_pt = 0.
63     adjust_volum = 0.
64    
65 jmc 1.1 DO bj=myByLo(myThid), myByHi(myThid)
66     DO bi=myBxLo(myThid), myBxHi(myThid)
67 jmc 1.2
68     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
69     C-- Compute the new fractional thickness of surface level (ksurfC):
70    
71     DO j=0,sNy+1
72     DO i=0,sNx+1
73     rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
74 jmc 1.1 ks = ksurfC(i,j,bi,bj)
75     IF (ks.LE.Nr) THEN
76 jmc 1.2 IF (rSurftmp(i,j) .LT. Rmin_surf(i,j,bi,bj)) THEN
77     C-- Needs to do something :
78 jmc 1.7 IF (numbWrite.LE.numbWrMax) THEN
79 jmc 1.13 numbWrite = numbWrite + 1
80 jmc 1.2 hFactmp = ( rSurftmp(i,j)-MAX(rF(ks+1),R_low(i,j,bi,bj))
81     & )*recip_drF(ks)
82 jmc 1.1 IF (hFactmp.LT.hFacInf) THEN
83 jmc 1.13 WRITE(errorMessageUnit,'(2A,6I4,I10)')
84 jmc 1.7 & 'WARNING: hFacC < hFacInf at:',
85 jmc 1.1 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
86 jmc 1.2 ELSE
87 jmc 1.13 WRITE(errorMessageUnit,'(2A,6I4,I10)')
88 jmc 1.7 & 'WARNING: hFac < hFacInf near:',
89 jmc 1.2 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
90 jmc 1.1 ENDIF
91 jmc 1.7 WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
92     & 'hFac_n-1,hFac_n,eta =',
93 jmc 1.2 & hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
94 jmc 1.7 ENDIF
95 jmc 1.2 C-- Decide to STOP :
96 jmc 1.13 c WRITE(errorMessageUnit,'(A)')
97 jmc 1.7 c & 'STOP in CALC_SURF_DR : too SMALL hFac !'
98 jmc 1.13 c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
99 jmc 1.2 C----------
100    
101     C-- Continue with Rmin_surf:
102     IF ( i.GE.1.AND.i.LE.sNx .AND.
103     & j.GE.1.AND.j.LE.sNy ) THEN
104     adjust_nb_pt = adjust_nb_pt + 1.
105 jmc 1.13 adjust_volum = adjust_volum
106 jmc 1.2 & + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
107     ENDIF
108     rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
109     C----------
110     ENDIF
111    
112     C-- Set hFac_surfC :
113 jmc 1.13 hFac_surfC(i,j,bi,bj) =
114 jmc 1.2 & ( rSurftmp(i,j) - MAX(rF(ks+1), R_low(i,j,bi,bj))
115     & )*recip_drF(ks)*maskC(i,j,ks,bi,bj)
116    
117     C-- Usefull warning when hFac becomes very large:
118 jmc 1.8 IF ( numbWrite.LE.numbWrMax .AND.
119     & hFac_surfC(i,j,bi,bj).GT.hFacSup ) THEN
120 jmc 1.13 numbWrite = numbWrite + 1
121     WRITE(errorMessageUnit,'(2A,6I4,I10)')
122 jmc 1.7 & 'WARNING: hFacC > hFacSup at:',
123 jmc 1.1 & ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
124 jmc 1.13 WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
125 jmc 1.7 & 'hFac_n-1,hFac_n,eta =', hfacC(i,j,ks,bi,bj),
126     & hFac_surfC(i,j,bi,bj), etaFld(i,j,bi,bj)
127 jmc 1.8 ENDIF
128 jmc 1.1 C----------
129     ENDIF
130 jmc 1.2
131 jmc 1.1 ENDDO
132     ENDDO
133    
134     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
135     C-- Compute fractional thickness of surface level, for U & V point:
136    
137     DO j=1,sNy
138     DO i=1,sNx+1
139     ks = ksurfW(i,j,bi,bj)
140     IF (ks.LE.Nr) THEN
141 jmc 1.2 hhm = rF(ks)
142     IF(ks.EQ.ksurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
143     hhp = rF(ks)
144 jmc 1.13 IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
145     hFac_surfW(i,j,bi,bj) =
146     & ( MIN(hhm,hhp)
147     & - MAX(rF(ks+1),R_low(i-1,j,bi,bj),R_low(i,j,bi,bj))
148 jmc 1.2 & )*recip_drF(ks)*maskW(i,j,ks,bi,bj)
149 jmc 1.1 ENDIF
150     ENDDO
151     ENDDO
152 jmc 1.2
153 jmc 1.1 DO j=1,sNy+1
154     DO i=1,sNx
155     ks = ksurfS(i,j,bi,bj)
156     IF (ks.LE.Nr) THEN
157 jmc 1.2 hhm = rF(ks)
158     IF(ks.EQ.ksurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
159     hhp = rF(ks)
160     IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
161 jmc 1.13 hFac_surfS(i,j,bi,bj) =
162     & ( MIN(hhm,hhp)
163     & - MAX(rF(ks+1),R_low(i,j-1,bi,bj),R_low(i,j,bi,bj))
164 jmc 1.2 & )*recip_drF(ks)*maskS(i,j,ks,bi,bj)
165 jmc 1.1 ENDIF
166     ENDDO
167     ENDDO
168 jmc 1.2
169     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
170    
171     C- end bi,bj loop.
172 jmc 1.1 ENDDO
173     ENDDO
174    
175 jmc 1.2 C-- Global diagnostic :
176 jmc 1.14 _GLOBAL_SUM_RL( adjust_nb_pt , myThid )
177     _GLOBAL_SUM_RL( adjust_volum , myThid )
178 jmc 1.2 IF (adjust_nb_pt .GE.1.) THEN
179 jmc 1.13 _BEGIN_MASTER( myThid )
180     WRITE(standardMessageUnit,'(2(A,I10),1PE16.8)')
181 jmc 1.7 & ' SURF_ADJUSTMENT: Iter=', myIter,
182     & ' Nb_pts,Vol=', nint(adjust_nb_pt), adjust_volum
183     _END_MASTER( myThid )
184 jmc 1.2 ENDIF
185    
186 jmc 1.14 _EXCH_XY_RS(hFac_surfC, myThid )
187 jmc 1.1 CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)
188 jmc 1.5
189 jmc 1.2 C-----
190 jmc 1.13 C Note: testing ksurfW,S is equivalent to a full height mask
191 jmc 1.2 C ==> no need for applying the mask here.
192     C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
193     C-----
194 jmc 1.1
195 jmc 1.12 c IF ( myIter.GE.0 ) THEN
196     c WRITE(suff,'(I10.10)') myIter
197     c CALL WRITE_FLD_XY_RS( 'hFac_surfC.', suff, hFac_surfC,
198 jmc 1.13 c & myIter, myThid )
199 jmc 1.12 c ENDIF
200 jmc 1.1
201     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
202     #endif /* NONLIN_FRSURF */
203    
204     RETURN
205     END

  ViewVC Help
Powered by ViewVC 1.1.22