model/src/calc_surf_dr.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_surf_dr.F,v 1.1 2001/08/27 18:50:10 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

      SUBROUTINE CALC_SURF_DR(etaFld,
     I                        myTime, myIter, myThid )
C     /==========================================================\
C     | SUBROUTINE CALC_SURF_DR                                  |
C     | o Calculate the new surface level thickness according to |
C     |   the surface r-position  (Non-Linear Free-Surf)         |
C     | o take decision if grid box becomes too thin or too thick|
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
c #include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"

C     == Routine arguments ==
C     myTime - Current time in simulation
C     myIter - Current iteration number in simulation
C     myThid - Thread number for this instance of the routine.
C     etaFld - current eta field used to update the hFactor
      _RL myTime
      INTEGER myIter
      INTEGER myThid
      _RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)

#ifdef NONLIN_FRSURF

C     Local variables in common block
C     Rmin_surf : minimum r_value of the free surface position 
C                  that satisfy  the hFacInf criteria
      COMMON /LOCAL_CALC_SURF_DR/ Rmin_surf
      _RL Rmin_surf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C     Local variables
C     i,j,k,bi,bj - loop counter
C     rSurftmp : free surface r-position that is used to compute hFac_surf
C     adjust_nb_pt : Nb of grid points where rSurf is adjusted (hFactInf)
C     adjust_volum : adjustment effect on the volume (domain size)

      INTEGER i,j,k,bi,bj
      INTEGER ks
      _RL hFacInfMOM, Rmin_tmp, hFactmp, adjust_nb_pt, adjust_volum
      _RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS hhm, hhp
      CHARACTER*(MAX_LEN_MBUF) suff

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF (myIter.EQ.nIter0) THEN

       hFacInfMOM = hFacInf

       DO bj=myByLo(myThid), myByHi(myThid)
        DO bi=myBxLo(myThid), myBxHi(myThid)

C-- Initialise arrays :
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx 
           hFac_surfC(i,j,bi,bj) = 0.
           hFac_surfW(i,j,bi,bj) = 0.
           hFac_surfS(i,j,bi,bj) = 0.
           Rmin_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
          ENDDO
         ENDDO

C-- Compute the mimimum value of r_surf (used for computing hFac_surfC)
         DO j=1,sNy
          DO i=1,sNx
           ks = ksurfC(i,j,bi,bj)
           IF (ks.LE.Nr) THEN
             Rmin_tmp = rF(ks+1)
             IF ( ks.EQ.ksurfW(i,j,bi,bj))
     &          Rmin_tmp = MAX(Rmin_tmp, R_low(i-1,j,bi,bj))
             IF ( ks.EQ.ksurfW(i+1,j,bi,bj)) 
     &          Rmin_tmp = MAX(Rmin_tmp, R_low(i+1,j,bi,bj))
             IF ( ks.EQ.ksurfS(i,j,bi,bj))
     &          Rmin_tmp = MAX(Rmin_tmp, R_low(i,j-1,bi,bj))
             IF ( ks.EQ.ksurfS(i,j+1,bi,bj))
     &          Rmin_tmp = MAX(Rmin_tmp, R_low(i,j+1,bi,bj))

             Rmin_surf(i,j,bi,bj) =
     &        MAX( MAX(rF(ks+1),R_low(i,j,bi,bj)) + hFacInf*drF(ks),  
     &                                Rmin_tmp + hFacInfMOM*drF(ks)
     &           )
           ENDIF
          ENDDO
         ENDDO

C-    end bi,bj loop.
        ENDDO
       ENDDO

       _EXCH_XY_R8( Rmin_surf, myThid ) 

C-    end of initialization block
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      adjust_nb_pt = 0.
      adjust_volum = 0.

      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Compute the new fractional thickness of surface level (ksurfC):

        DO j=0,sNy+1
         DO i=0,sNx+1
          rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
          ks = ksurfC(i,j,bi,bj)
          IF (ks.LE.Nr) THEN
           IF (rSurftmp(i,j) .LT. Rmin_surf(i,j,bi,bj)) THEN
C-- Needs to do something :
             hFactmp = ( rSurftmp(i,j)-MAX(rF(ks+1),R_low(i,j,bi,bj))
     &                 )*recip_drF(ks)
             IF (hFactmp.LT.hFacInf) THEN
              write(0,'(2A,6I4,I10)') 'WARNING: hFacC < hFacInf at:',
     &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
             ELSE
              write(0,'(2A,6I4,I10)') 'WARNING: hFac < hFacInf near:',
     &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
             ENDIF
              write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
     &          hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
C-- Decide to STOP :
c             write(0,'(2A)') 'STOP in CALC_SURF_DR :',
c    &                        ' too SMALL hFac !'
c             STOP 'ABNORMAL END: S/R CALC_SURF_DR' 
C----------

C-- Continue with Rmin_surf:
             IF ( i.GE.1.AND.i.LE.sNx .AND.
     &            j.GE.1.AND.j.LE.sNy ) THEN
               adjust_nb_pt = adjust_nb_pt + 1.
               adjust_volum = adjust_volum 
     &          + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
             ENDIF
             rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
C----------
           ENDIF

C-- Set hFac_surfC :
           hFac_surfC(i,j,bi,bj) = 
     &         ( rSurftmp(i,j) - MAX(rF(ks+1), R_low(i,j,bi,bj))
     &         )*recip_drF(ks)*maskC(i,j,ks,bi,bj)

           IF (hFac_surfC(i,j,bi,bj).GT.hFacSup) THEN
C-- Usefull warning when hFac becomes very large:
              write(0,'(2A,6I4,I10)') 'WARNING: hFacC > hFacSup at:',
     &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
              write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
     &          hfacC(i,j,ks,bi,bj), hFac_surfC(i,j,bi,bj), 
     &          etaFld(i,j,bi,bj)
C-- Decide to STOP :
c             write(0,'(2A)') 'STOP in CALC_SURF_DR :',
c    &                        ' too LARGE hFac !'
c             STOP 'ABNORMAL END: S/R CALC_SURF_DR' 
C----------
           ENDIF
          ENDIF

         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Compute fractional thickness of surface level, for U & V point:

        DO j=1,sNy
         DO i=1,sNx+1
          ks = ksurfW(i,j,bi,bj)
          IF (ks.LE.Nr) THEN
            hhm = rF(ks)
            IF(ks.EQ.ksurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
            hhp = rF(ks)
            IF(ks.EQ.ksurfC(i,j,bi,bj))   hhp = rSurftmp(i,j)  
            hFac_surfW(i,j,bi,bj) = 
     &         ( MIN(hhm,hhp) 
     &          - MAX(rF(ks+1),R_low(i-1,j,bi,bj),R_low(i,j,bi,bj)) 
     &         )*recip_drF(ks)*maskW(i,j,ks,bi,bj)
          ENDIF
         ENDDO
        ENDDO

        DO j=1,sNy+1
         DO i=1,sNx
          ks = ksurfS(i,j,bi,bj)
          IF (ks.LE.Nr) THEN
            hhm = rF(ks)
            IF(ks.EQ.ksurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
            hhp = rF(ks)
            IF(ks.EQ.ksurfC(i,j,bi,bj))   hhp = rSurftmp(i,j)
            hFac_surfS(i,j,bi,bj) = 
     &         ( MIN(hhm,hhp) 
     &          - MAX(rF(ks+1),R_low(i,j-1,bi,bj),R_low(i,j,bi,bj)) 
     &         )*recip_drF(ks)*maskS(i,j,ks,bi,bj)
          ENDIF
         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C-    end bi,bj loop.
       ENDDO
      ENDDO

C-- Global diagnostic :
      _GLOBAL_SUM_R8( adjust_nb_pt , myThid ) 
      _GLOBAL_SUM_R8( adjust_volum , myThid ) 
      IF (adjust_nb_pt .GE.1.) THEN
        _BEGIN_MASTER( myThid ) 
        write(*,'(2(A,I10),1PE16.8)') ' SURF_ADJUSTMENT: Iter=',
     &   myIter, ' Nb_pts,Vol=', nint(adjust_nb_pt), adjust_volum
        _END_MASTER( )
      ENDIF

      _EXCH_XY_R4(hFac_surfC, myThid ) 
      CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)

C-----
C Note: testing ksurfW,S is equivalent to a full height mask 
C   ==> no need for applying the mask here.
C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
C-----

      WRITE(suff,'(I10.10)') myIter
c     CALL WRITE_FLD_XY_RS('hFac_surfC.',suff,hFac_surfC,
c    &                     myIter,myThid)                    

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* NONLIN_FRSURF */

      RETURN
      END
1	jmc	1.2	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_surf_dr.F,v 1.1 2001/08/27 18:50:10 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			SUBROUTINE CALC_SURF_DR(etaFld,
7			I myTime, myIter, myThid )
8			C /==========================================================\
9			C \| SUBROUTINE CALC_SURF_DR \|
10			C \| o Calculate the new surface level thickness according to \|
11			C \| the surface r-position (Non-Linear Free-Surf) \|
12			C \| o take decision if grid box becomes too thin or too thick\|
13			C \==========================================================/
14			IMPLICIT NONE
15
16			C == Global variables
17			#include "SIZE.h"
18			#include "EEPARAMS.h"
19			#include "PARAMS.h"
20			c #include "DYNVARS.h"
21			#include "GRID.h"
22			#include "SURFACE.h"
23
24			C == Routine arguments ==
25			C myTime - Current time in simulation
26			C myIter - Current iteration number in simulation
27			C myThid - Thread number for this instance of the routine.
28			C etaFld - current eta field used to update the hFactor
29			_RL myTime
30			INTEGER myIter
31			INTEGER myThid
32			_RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
33
34			#ifdef NONLIN_FRSURF
35
36	jmc	1.2	C Local variables in common block
37			C Rmin_surf : minimum r_value of the free surface position
38			C that satisfy the hFacInf criteria
39			COMMON /LOCAL_CALC_SURF_DR/ Rmin_surf
40			_RL Rmin_surf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
41
42	jmc	1.1	C Local variables
43			C i,j,k,bi,bj - loop counter
44	jmc	1.2	C rSurftmp : free surface r-position that is used to compute hFac_surf
45			C adjust_nb_pt : Nb of grid points where rSurf is adjusted (hFactInf)
46			C adjust_volum : adjustment effect on the volume (domain size)
47
48	jmc	1.1	INTEGER i,j,k,bi,bj
49			INTEGER ks
50	jmc	1.2	_RL hFacInfMOM, Rmin_tmp, hFactmp, adjust_nb_pt, adjust_volum
51			_RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52	jmc	1.1	_RS hhm, hhp
53			CHARACTER*(MAX_LEN_MBUF) suff
54
55			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
56
57			IF (myIter.EQ.nIter0) THEN
58	jmc	1.2
59			hFacInfMOM = hFacInf
60
61	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
62			DO bi=myBxLo(myThid), myBxHi(myThid)
63	jmc	1.2
64			C-- Initialise arrays :
65	jmc	1.1	DO j=1-Oly,sNy+Oly
66			DO i=1-Olx,sNx+Olx
67			hFac_surfC(i,j,bi,bj) = 0.
68			hFac_surfW(i,j,bi,bj) = 0.
69			hFac_surfS(i,j,bi,bj) = 0.
70	jmc	1.2	Rmin_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
71	jmc	1.1	ENDDO
72			ENDDO
73	jmc	1.2
74			C-- Compute the mimimum value of r_surf (used for computing hFac_surfC)
75			DO j=1,sNy
76			DO i=1,sNx
77			ks = ksurfC(i,j,bi,bj)
78			IF (ks.LE.Nr) THEN
79			Rmin_tmp = rF(ks+1)
80			IF ( ks.EQ.ksurfW(i,j,bi,bj))
81			& Rmin_tmp = MAX(Rmin_tmp, R_low(i-1,j,bi,bj))
82			IF ( ks.EQ.ksurfW(i+1,j,bi,bj))
83			& Rmin_tmp = MAX(Rmin_tmp, R_low(i+1,j,bi,bj))
84			IF ( ks.EQ.ksurfS(i,j,bi,bj))
85			& Rmin_tmp = MAX(Rmin_tmp, R_low(i,j-1,bi,bj))
86			IF ( ks.EQ.ksurfS(i,j+1,bi,bj))
87			& Rmin_tmp = MAX(Rmin_tmp, R_low(i,j+1,bi,bj))
88
89			Rmin_surf(i,j,bi,bj) =
90			& MAX( MAX(rF(ks+1),R_low(i,j,bi,bj)) + hFacInf*drF(ks),
91			& Rmin_tmp + hFacInfMOM*drF(ks)
92			& )
93			ENDIF
94			ENDDO
95			ENDDO
96
97			C- end bi,bj loop.
98	jmc	1.1	ENDDO
99			ENDDO
100	jmc	1.2
101			_EXCH_XY_R8( Rmin_surf, myThid )
102
103			C- end of initialization block
104	jmc	1.1	ENDIF
105
106			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
107
108	jmc	1.2	adjust_nb_pt = 0.
109			adjust_volum = 0.
110
111	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
112			DO bi=myBxLo(myThid), myBxHi(myThid)
113	jmc	1.2
114			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
115			C-- Compute the new fractional thickness of surface level (ksurfC):
116
117			DO j=0,sNy+1
118			DO i=0,sNx+1
119			rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
120	jmc	1.1	ks = ksurfC(i,j,bi,bj)
121			IF (ks.LE.Nr) THEN
122	jmc	1.2	IF (rSurftmp(i,j) .LT. Rmin_surf(i,j,bi,bj)) THEN
123			C-- Needs to do something :
124			hFactmp = ( rSurftmp(i,j)-MAX(rF(ks+1),R_low(i,j,bi,bj))
125			& )*recip_drF(ks)
126	jmc	1.1	IF (hFactmp.LT.hFacInf) THEN
127			write(0,'(2A,6I4,I10)') 'WARNING: hFacC < hFacInf at:',
128			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
129	jmc	1.2	ELSE
130			write(0,'(2A,6I4,I10)') 'WARNING: hFac < hFacInf near:',
131			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
132	jmc	1.1	ENDIF
133	jmc	1.2	write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
134			& hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
135			C-- Decide to STOP :
136			c write(0,'(2A)') 'STOP in CALC_SURF_DR :',
137			c & ' too SMALL hFac !'
138			c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
139			C----------
140
141			C-- Continue with Rmin_surf:
142			IF ( i.GE.1.AND.i.LE.sNx .AND.
143			& j.GE.1.AND.j.LE.sNy ) THEN
144			adjust_nb_pt = adjust_nb_pt + 1.
145			adjust_volum = adjust_volum
146			& + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
147			ENDIF
148			rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
149			C----------
150			ENDIF
151
152			C-- Set hFac_surfC :
153			hFac_surfC(i,j,bi,bj) =
154			& ( rSurftmp(i,j) - MAX(rF(ks+1), R_low(i,j,bi,bj))
155			& )recip_drF(ks)maskC(i,j,ks,bi,bj)
156
157			IF (hFac_surfC(i,j,bi,bj).GT.hFacSup) THEN
158			C-- Usefull warning when hFac becomes very large:
159	jmc	1.1	write(0,'(2A,6I4,I10)') 'WARNING: hFacC > hFacSup at:',
160			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
161			write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
162	jmc	1.2	& hfacC(i,j,ks,bi,bj), hFac_surfC(i,j,bi,bj),
163			& etaFld(i,j,bi,bj)
164			C-- Decide to STOP :
165			c write(0,'(2A)') 'STOP in CALC_SURF_DR :',
166			c & ' too LARGE hFac !'
167			c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
168	jmc	1.1	C----------
169			ENDIF
170			ENDIF
171	jmc	1.2
172	jmc	1.1	ENDDO
173			ENDDO
174
175			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
176			C-- Compute fractional thickness of surface level, for U & V point:
177
178			DO j=1,sNy
179			DO i=1,sNx+1
180			ks = ksurfW(i,j,bi,bj)
181			IF (ks.LE.Nr) THEN
182	jmc	1.2	hhm = rF(ks)
183			IF(ks.EQ.ksurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
184			hhp = rF(ks)
185			IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
186			hFac_surfW(i,j,bi,bj) =
187			& ( MIN(hhm,hhp)
188			& - MAX(rF(ks+1),R_low(i-1,j,bi,bj),R_low(i,j,bi,bj))
189			& )recip_drF(ks)maskW(i,j,ks,bi,bj)
190	jmc	1.1	ENDIF
191			ENDDO
192			ENDDO
193	jmc	1.2
194	jmc	1.1	DO j=1,sNy+1
195			DO i=1,sNx
196			ks = ksurfS(i,j,bi,bj)
197			IF (ks.LE.Nr) THEN
198	jmc	1.2	hhm = rF(ks)
199			IF(ks.EQ.ksurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
200			hhp = rF(ks)
201			IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
202			hFac_surfS(i,j,bi,bj) =
203			& ( MIN(hhm,hhp)
204			& - MAX(rF(ks+1),R_low(i,j-1,bi,bj),R_low(i,j,bi,bj))
205			& )recip_drF(ks)maskS(i,j,ks,bi,bj)
206	jmc	1.1	ENDIF
207			ENDDO
208			ENDDO
209	jmc	1.2
210			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
211
212			C- end bi,bj loop.
213	jmc	1.1	ENDDO
214			ENDDO
215
216	jmc	1.2	C-- Global diagnostic :
217			_GLOBAL_SUM_R8( adjust_nb_pt , myThid )
218			_GLOBAL_SUM_R8( adjust_volum , myThid )
219			IF (adjust_nb_pt .GE.1.) THEN
220			_BEGIN_MASTER( myThid )
221			write(*,'(2(A,I10),1PE16.8)') ' SURF_ADJUSTMENT: Iter=',
222			& myIter, ' Nb_pts,Vol=', nint(adjust_nb_pt), adjust_volum
223			_END_MASTER( )
224			ENDIF
225
226			_EXCH_XY_R4(hFac_surfC, myThid )
227	jmc	1.1	CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)
228
229	jmc	1.2	C-----
230			C Note: testing ksurfW,S is equivalent to a full height mask
231			C ==> no need for applying the mask here.
232			C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
233			C-----
234	jmc	1.1
235			WRITE(suff,'(I10.10)') myIter
236			c CALL WRITE_FLD_XY_RS('hFac_surfC.',suff,hFac_surfC,
237			c & myIter,myThid)
238
239			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
240			#endif /* NONLIN_FRSURF */
241
242			RETURN
243			END