model/src/calc_surf_dr.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_surf_dr.F,v 1.3 2001/09/26 18:09:14 cnh Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_SURF_DR
C     !INTERFACE:
      SUBROUTINE CALC_SURF_DR(etaFld,
     I                        myTime, myIter, myThid )
C     !DESCRIPTION: \bv
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     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C     !INPUT/OUTPUT PARAMETERS:
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:
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
CEOP

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

      IF (myIter.EQ.-1) 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.4	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_surf_dr.F,v 1.3 2001/09/26 18:09:14 cnh 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.1	INTEGER i,j,k,bi,bj
54			INTEGER ks
55	jmc	1.2	_RL hFacInfMOM, Rmin_tmp, hFactmp, adjust_nb_pt, adjust_volum
56			_RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57	jmc	1.1	_RS hhm, hhp
58			CHARACTER*(MAX_LEN_MBUF) suff
59	cnh	1.3	CEOP
60	jmc	1.1
61			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
62
63	jmc	1.4	IF (myIter.EQ.-1) THEN
64	jmc	1.2
65			hFacInfMOM = hFacInf
66
67	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
68			DO bi=myBxLo(myThid), myBxHi(myThid)
69	jmc	1.2
70			C-- Initialise arrays :
71	jmc	1.1	DO j=1-Oly,sNy+Oly
72			DO i=1-Olx,sNx+Olx
73			hFac_surfC(i,j,bi,bj) = 0.
74			hFac_surfW(i,j,bi,bj) = 0.
75			hFac_surfS(i,j,bi,bj) = 0.
76	jmc	1.2	Rmin_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
77	jmc	1.1	ENDDO
78			ENDDO
79	jmc	1.2
80			C-- Compute the mimimum value of r_surf (used for computing hFac_surfC)
81			DO j=1,sNy
82			DO i=1,sNx
83			ks = ksurfC(i,j,bi,bj)
84			IF (ks.LE.Nr) THEN
85			Rmin_tmp = rF(ks+1)
86			IF ( ks.EQ.ksurfW(i,j,bi,bj))
87			& Rmin_tmp = MAX(Rmin_tmp, R_low(i-1,j,bi,bj))
88			IF ( ks.EQ.ksurfW(i+1,j,bi,bj))
89			& Rmin_tmp = MAX(Rmin_tmp, R_low(i+1,j,bi,bj))
90			IF ( ks.EQ.ksurfS(i,j,bi,bj))
91			& Rmin_tmp = MAX(Rmin_tmp, R_low(i,j-1,bi,bj))
92			IF ( ks.EQ.ksurfS(i,j+1,bi,bj))
93			& Rmin_tmp = MAX(Rmin_tmp, R_low(i,j+1,bi,bj))
94
95			Rmin_surf(i,j,bi,bj) =
96			& MAX( MAX(rF(ks+1),R_low(i,j,bi,bj)) + hFacInf*drF(ks),
97			& Rmin_tmp + hFacInfMOM*drF(ks)
98			& )
99			ENDIF
100			ENDDO
101			ENDDO
102
103			C- end bi,bj loop.
104	jmc	1.1	ENDDO
105			ENDDO
106	jmc	1.2
107			_EXCH_XY_R8( Rmin_surf, myThid )
108
109			C- end of initialization block
110	jmc	1.1	ENDIF
111
112			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
113
114	jmc	1.2	adjust_nb_pt = 0.
115			adjust_volum = 0.
116
117	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
118			DO bi=myBxLo(myThid), myBxHi(myThid)
119	jmc	1.2
120			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
121			C-- Compute the new fractional thickness of surface level (ksurfC):
122
123			DO j=0,sNy+1
124			DO i=0,sNx+1
125			rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
126	jmc	1.1	ks = ksurfC(i,j,bi,bj)
127			IF (ks.LE.Nr) THEN
128	jmc	1.2	IF (rSurftmp(i,j) .LT. Rmin_surf(i,j,bi,bj)) THEN
129			C-- Needs to do something :
130			hFactmp = ( rSurftmp(i,j)-MAX(rF(ks+1),R_low(i,j,bi,bj))
131			& )*recip_drF(ks)
132	jmc	1.1	IF (hFactmp.LT.hFacInf) THEN
133			write(0,'(2A,6I4,I10)') 'WARNING: hFacC < hFacInf at:',
134			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
135	jmc	1.2	ELSE
136			write(0,'(2A,6I4,I10)') 'WARNING: hFac < hFacInf near:',
137			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
138	jmc	1.1	ENDIF
139	jmc	1.2	write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
140			& hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
141			C-- Decide to STOP :
142			c write(0,'(2A)') 'STOP in CALC_SURF_DR :',
143			c & ' too SMALL hFac !'
144			c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
145			C----------
146
147			C-- Continue with Rmin_surf:
148			IF ( i.GE.1.AND.i.LE.sNx .AND.
149			& j.GE.1.AND.j.LE.sNy ) THEN
150			adjust_nb_pt = adjust_nb_pt + 1.
151			adjust_volum = adjust_volum
152			& + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
153			ENDIF
154			rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
155			C----------
156			ENDIF
157
158			C-- Set hFac_surfC :
159			hFac_surfC(i,j,bi,bj) =
160			& ( rSurftmp(i,j) - MAX(rF(ks+1), R_low(i,j,bi,bj))
161			& )recip_drF(ks)maskC(i,j,ks,bi,bj)
162
163			IF (hFac_surfC(i,j,bi,bj).GT.hFacSup) THEN
164			C-- Usefull warning when hFac becomes very large:
165	jmc	1.1	write(0,'(2A,6I4,I10)') 'WARNING: hFacC > hFacSup at:',
166			& ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
167			write(0,'(A,2F10.6,1PE14.6)') 'hFac_n-1,hFac_n,eta =',
168	jmc	1.2	& hfacC(i,j,ks,bi,bj), hFac_surfC(i,j,bi,bj),
169			& etaFld(i,j,bi,bj)
170			C-- Decide to STOP :
171			c write(0,'(2A)') 'STOP in CALC_SURF_DR :',
172			c & ' too LARGE hFac !'
173			c STOP 'ABNORMAL END: S/R CALC_SURF_DR'
174	jmc	1.1	C----------
175			ENDIF
176			ENDIF
177	jmc	1.2
178	jmc	1.1	ENDDO
179			ENDDO
180
181			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
182			C-- Compute fractional thickness of surface level, for U & V point:
183
184			DO j=1,sNy
185			DO i=1,sNx+1
186			ks = ksurfW(i,j,bi,bj)
187			IF (ks.LE.Nr) THEN
188	jmc	1.2	hhm = rF(ks)
189			IF(ks.EQ.ksurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
190			hhp = rF(ks)
191			IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
192			hFac_surfW(i,j,bi,bj) =
193			& ( MIN(hhm,hhp)
194			& - MAX(rF(ks+1),R_low(i-1,j,bi,bj),R_low(i,j,bi,bj))
195			& )recip_drF(ks)maskW(i,j,ks,bi,bj)
196	jmc	1.1	ENDIF
197			ENDDO
198			ENDDO
199	jmc	1.2
200	jmc	1.1	DO j=1,sNy+1
201			DO i=1,sNx
202			ks = ksurfS(i,j,bi,bj)
203			IF (ks.LE.Nr) THEN
204	jmc	1.2	hhm = rF(ks)
205			IF(ks.EQ.ksurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
206			hhp = rF(ks)
207			IF(ks.EQ.ksurfC(i,j,bi,bj)) hhp = rSurftmp(i,j)
208			hFac_surfS(i,j,bi,bj) =
209			& ( MIN(hhm,hhp)
210			& - MAX(rF(ks+1),R_low(i,j-1,bi,bj),R_low(i,j,bi,bj))
211			& )recip_drF(ks)maskS(i,j,ks,bi,bj)
212	jmc	1.1	ENDIF
213			ENDDO
214			ENDDO
215	jmc	1.2
216			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
217
218			C- end bi,bj loop.
219	jmc	1.1	ENDDO
220			ENDDO
221
222	jmc	1.2	C-- Global diagnostic :
223			_GLOBAL_SUM_R8( adjust_nb_pt , myThid )
224			_GLOBAL_SUM_R8( adjust_volum , myThid )
225			IF (adjust_nb_pt .GE.1.) THEN
226			_BEGIN_MASTER( myThid )
227			write(*,'(2(A,I10),1PE16.8)') ' SURF_ADJUSTMENT: Iter=',
228			& myIter, ' Nb_pts,Vol=', nint(adjust_nb_pt), adjust_volum
229			_END_MASTER( )
230			ENDIF
231
232			_EXCH_XY_R4(hFac_surfC, myThid )
233	jmc	1.1	CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)
234
235	jmc	1.2	C-----
236			C Note: testing ksurfW,S is equivalent to a full height mask
237			C ==> no need for applying the mask here.
238			C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
239			C-----
240	jmc	1.1
241			WRITE(suff,'(I10.10)') myIter
242			c CALL WRITE_FLD_XY_RS('hFac_surfC.',suff,hFac_surfC,
243			c & myIter,myThid)
244
245			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
246			#endif /* NONLIN_FRSURF */
247
248			RETURN
249			END