model/src/calc_r_star.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.5 2005/06/19 23:26:25 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_R_STAR
C     !INTERFACE:
      SUBROUTINE CALC_R_STAR( etaFld,
     I                        myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE CALC_R_STAR                                   
C     | o Calculate new column thickness & scaling factor for r*
C     |    according to the surface r-position (Non-Lin Free-Surf)          
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
C     i,j,k,bi,bj  :: loop counter
C     numbWrite    :: count the Number of warning written on STD-ERR file
C     numbWrMax    ::  maximum  Number of warning written on STD-ERR file
      INTEGER i,j,k,bi,bj
      INTEGER ii,jj
      INTEGER km, numbWrite, numbWrMax
      _RL tmpfldW, tmpfldS
c     CHARACTER*(MAX_LEN_MBUF) suff
CEOP
      DATA numbWrite / 0 /
      numbWrMax = Nx*Ny

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

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
#endif

      IF (usingPCoords) THEN
       km = Nr
      ELSE
       km = 1
      ENDIF 

      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
C-    1rst bi,bj loop :

        IF (myIter.EQ.-1) THEN
C-- Initialise arrays :
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx 
            rStarFacC(i,j,bi,bj) = 1.
            rStarFacW(i,j,bi,bj) = 1.
            rStarFacS(i,j,bi,bj) = 1.
            rStarExpC(i,j,bi,bj) = 1.
            rStarExpW(i,j,bi,bj) = 1.
            rStarExpS(i,j,bi,bj) = 1.
            rStarDhCDt(i,j,bi,bj) = 0.
            rStarDhWDt(i,j,bi,bj) = 0.
            rStarDhSDt(i,j,bi,bj) = 0.
            PmEpR(i,j,bi,bj) = 0.
          ENDDO
         ENDDO
         DO k=1,Nr
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx 
            h0FacC(i,j,k,bi,bj) = hFacC(i,j,k,bi,bj)
            h0FacW(i,j,k,bi,bj) = hFacW(i,j,k,bi,bj)
            h0FacS(i,j,k,bi,bj) = hFacS(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ELSE
C-- copy rStarFacX -> rStarExpX
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx 
            rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
            rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
            rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF

C-- Compute the new column thikness :
        DO j=0,sNy+1
         DO i=0,sNx+1
          IF (maskH(i,j,bi,bj).EQ.1. ) THEN
           rStarFacC(i,j,bi,bj) = 
     &      (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
     &      *recip_Rcol(i,j,bi,bj)
          ELSE
           rStarFacC(i,j,bi,bj) = 1.
          ENDIF
         ENDDO
        ENDDO
        DO j=1,sNy
         DO i=1,sNx+1
          IF (maskW(i,j,km,bi,bj).EQ.1. ) THEN
           tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
     &             - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
           rStarFacW(i,j,bi,bj) = 
     &       ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
     &                    +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
     &                   )*recip_rAw(i,j,bi,bj)
     &        +tmpfldW )/tmpfldW
          ELSE
           rStarFacW(i,j,bi,bj) = 1.
          ENDIF
         ENDDO
        ENDDO
        DO j=1,sNy+1
         DO i=1,sNx
          IF (maskS(i,j,km,bi,bj).EQ.1. ) THEN
           tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
     &             - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
           rStarFacS(i,j,bi,bj) = 
     &       ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
     &                    +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
     &                   )*recip_rAs(i,j,bi,bj)
     &        +tmpfldS )/tmpfldS
          ELSE
           rStarFacS(i,j,bi,bj) = 1.
          ENDIF
         ENDDO
        ENDDO

C-    Needs to do something when r* ratio is too small ;
C     for now, just stop
        DO j=1,sNy+1
         DO i=1,sNx+1
          IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN  
            numbWrite = numbWrite + 1
            WRITE(errorMessageUnit,'(2A,5I4,I10)')
     &       'WARNING: r*FacC < hFacInf at:',
     &       ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
            WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
     &       'rStarFac,H,eta =', rStarFacC(i,j,bi,bj), 
     &       Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
            WRITE(errorMessageUnit,'(A)')
     &       'STOP in CALC_R_STAR : too SMALL rStarFacC !'
             STOP 'ABNORMAL END: S/R CALC_SURF_DR'  
          ENDIF
          IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN  
            numbWrite = numbWrite + 1
            tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
     &              - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
            WRITE(errorMessageUnit,'(2A,5I4,I10)')
     &       'WARNING: r*FacW < hFacInf at:',
     &       ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
            WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
     &       'rStarFac,H,eta =', rStarFacW(i,j,bi,bj), tmpfldW,
     &        etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
            WRITE(errorMessageUnit,'(A)')
     &       'STOP in CALC_R_STAR : too SMALL rStarFacW !'
             STOP 'ABNORMAL END: S/R CALC_SURF_DR'  
          ENDIF
          IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN  
            numbWrite = numbWrite + 1
            tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
     &              - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
            WRITE(errorMessageUnit,'(2A,5I4,I10)')
     &       'WARNING: r*FacS < hFacInf at:',
     &       ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
            WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
     &       'rStarFac,H,eta =', rStarFacS(i,j,bi,bj), tmpfldS,
     &        etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
            WRITE(errorMessageUnit,'(A)')
     &       'STOP in CALC_R_STAR : too SMALL rStarFacS !'
             STOP 'ABNORMAL END: S/R CALC_R_STAR'  
          ENDIF
C-- Usefull warning when r*Fac becomes very large:
          IF ( numbWrite.LE.numbWrMax .AND.
     &         rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
            numbWrite = numbWrite + 1
            WRITE(errorMessageUnit,'(2A,5I4,I10)')
     &       'WARNING: hFacC > hFacSup at:',
     &       ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
            WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
     &       'rStarFac,H,eta =', rStarFacC(i,j,bi,bj), 
     &       Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO

C-    end 1rst bi,bj loop.
       ENDDO
      ENDDO

       _EXCH_XY_RL( rStarFacC, myThid ) 
      CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)

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

      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
C-    2nd bi,bj loop :

#ifdef ALLOW_EXCH2
C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
C        in the corner regions of the tile (e.g.:[1-Olx:0,1-Oly:0])
C       => need to add those lines (or to fix exch2):
        DO j=1,Oly
         DO i=1,Olx
          ii = sNx+i
          jj = sNy+j

          IF (maskH(1-i,1-j,bi,bj).EQ.0.) rStarFacC(1-i,1-j,bi,bj) = 1.
          IF (maskH(ii, 1-j,bi,bj).EQ.0.) rStarFacC(ii, 1-j,bi,bj) = 1.
          IF (maskH(1-i,jj, bi,bj).EQ.0.) rStarFacC(1-i,jj, bi,bj) = 1.
          IF (maskH(ii, jj, bi,bj).EQ.0.) rStarFacC(ii, jj, bi,bj) = 1.

c         IF (ksurfW(1-i,1-j,bi,bj).LE.Nr) rStarFacW(1-i,1-j,bi,bj)=1.
          IF (maskW(1-i,1-j,km,bi,bj).EQ.0.) rStarFacW(1-i,1-j,bi,bj)=1.
          IF (maskW(ii, 1-j,km,bi,bj).EQ.0.) rStarFacW(ii, 1-j,bi,bj)=1.
          IF (maskW(1-i,jj, km,bi,bj).EQ.0.) rStarFacW(1-i,jj, bi,bj)=1.
          IF (maskW(ii, jj, km,bi,bj).EQ.0.) rStarFacW(ii, jj, bi,bj)=1.

          IF (maskS(1-i,1-j,km,bi,bj).EQ.0.) rStarFacS(1-i,1-j,bi,bj)=1.
          IF (maskS(ii, 1-j,km,bi,bj).EQ.0.) rStarFacS(ii, 1-j,bi,bj)=1.
          IF (maskS(1-i,jj, km,bi,bj).EQ.0.) rStarFacS(1-i,jj, bi,bj)=1.
          IF (maskS(ii, jj, km,bi,bj).EQ.0.) rStarFacS(ii, jj, bi,bj)=1.
         ENDDO
        ENDDO
#endif /* ALLOW_EXCH2 */

        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx 
           rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj) 
     &                           -rStarExpC(i,j,bi,bj))/deltaTfreesurf
           rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
     &                           -rStarExpW(i,j,bi,bj))/deltaTfreesurf
           rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
     &                           -rStarExpS(i,j,bi,bj))/deltaTfreesurf
           rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
     &                          / rStarExpC(i,j,bi,bj)
           rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
     &                          / rStarExpW(i,j,bi,bj)
           rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
     &                          / rStarExpS(i,j,bi,bj)
         ENDDO
        ENDDO

C-    end 2nd bi,bj loop.
        ENDDO
       ENDDO

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
#endif

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

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.5 2005/06/19 23:26:25 jmc Exp $
2	jmc	1.1	C $Name: $
3
4	jmc	1.4	#include "PACKAGES_CONFIG.h"
5	jmc	1.1	#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: CALC_R_STAR
9			C !INTERFACE:
10			SUBROUTINE CALC_R_STAR( etaFld,
11			I myTime, myIter, myThid )
12			C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE CALC_R_STAR
15			C \| o Calculate new column thickness & scaling factor for r*
16			C \| according to the surface r-position (Non-Lin Free-Surf)
17			C ==========================================================
18			C \ev
19
20			C !USES:
21			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			C !INPUT/OUTPUT PARAMETERS:
30			C == Routine arguments ==
31			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			_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			C !LOCAL VARIABLES:
43			C Local variables
44			C i,j,k,bi,bj :: loop counter
45	jmc	1.2	C numbWrite :: count the Number of warning written on STD-ERR file
46			C numbWrMax :: maximum Number of warning written on STD-ERR file
47	jmc	1.1	INTEGER i,j,k,bi,bj
48	jmc	1.4	INTEGER ii,jj
49	jmc	1.2	INTEGER km, numbWrite, numbWrMax
50	jmc	1.1	_RL tmpfldW, tmpfldS
51			c CHARACTER*(MAX_LEN_MBUF) suff
52			CEOP
53	jmc	1.2	DATA numbWrite / 0 /
54			numbWrMax = Nx*Ny
55	jmc	1.1
56			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
57
58	jmc	1.4	#ifdef ALLOW_DEBUG
59			IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
60			#endif
61
62	jmc	1.6	IF (usingPCoords) THEN
63	jmc	1.5	km = Nr
64			ELSE
65	jmc	1.1	km = 1
66			ENDIF
67
68			DO bj=myByLo(myThid), myByHi(myThid)
69			DO bi=myBxLo(myThid), myBxHi(myThid)
70			C- 1rst bi,bj loop :
71
72			IF (myIter.EQ.-1) THEN
73			C-- Initialise arrays :
74			DO j=1-Oly,sNy+Oly
75			DO i=1-Olx,sNx+Olx
76			rStarFacC(i,j,bi,bj) = 1.
77			rStarFacW(i,j,bi,bj) = 1.
78			rStarFacS(i,j,bi,bj) = 1.
79			rStarExpC(i,j,bi,bj) = 1.
80			rStarExpW(i,j,bi,bj) = 1.
81			rStarExpS(i,j,bi,bj) = 1.
82			rStarDhCDt(i,j,bi,bj) = 0.
83			rStarDhWDt(i,j,bi,bj) = 0.
84			rStarDhSDt(i,j,bi,bj) = 0.
85			PmEpR(i,j,bi,bj) = 0.
86			ENDDO
87			ENDDO
88			DO k=1,Nr
89			DO j=1-Oly,sNy+Oly
90			DO i=1-Olx,sNx+Olx
91			h0FacC(i,j,k,bi,bj) = hFacC(i,j,k,bi,bj)
92			h0FacW(i,j,k,bi,bj) = hFacW(i,j,k,bi,bj)
93			h0FacS(i,j,k,bi,bj) = hFacS(i,j,k,bi,bj)
94			ENDDO
95			ENDDO
96			ENDDO
97			ELSE
98			C-- copy rStarFacX -> rStarExpX
99			DO j=1-Oly,sNy+Oly
100			DO i=1-Olx,sNx+Olx
101			rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
102			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
103			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
104			ENDDO
105			ENDDO
106			ENDIF
107
108			C-- Compute the new column thikness :
109			DO j=0,sNy+1
110			DO i=0,sNx+1
111			IF (maskH(i,j,bi,bj).EQ.1. ) THEN
112			rStarFacC(i,j,bi,bj) =
113			& (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
114			& *recip_Rcol(i,j,bi,bj)
115			ELSE
116			rStarFacC(i,j,bi,bj) = 1.
117			ENDIF
118			ENDDO
119			ENDDO
120	jmc	1.2	DO j=1,sNy
121	jmc	1.1	DO i=1,sNx+1
122			IF (maskW(i,j,km,bi,bj).EQ.1. ) THEN
123			tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
124			& - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
125			rStarFacW(i,j,bi,bj) =
126			& ( 0.5 _d 0 ( etaFld(i-1,j,bi,bj)rA(i-1,j,bi,bj)
127			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
128			& )*recip_rAw(i,j,bi,bj)
129			& +tmpfldW )/tmpfldW
130			ELSE
131			rStarFacW(i,j,bi,bj) = 1.
132			ENDIF
133	jmc	1.2	ENDDO
134			ENDDO
135			DO j=1,sNy+1
136			DO i=1,sNx
137	jmc	1.1	IF (maskS(i,j,km,bi,bj).EQ.1. ) THEN
138			tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
139			& - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
140			rStarFacS(i,j,bi,bj) =
141			& ( 0.5 _d 0 ( etaFld(i,j-1,bi,bj)rA(i,j-1,bi,bj)
142			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
143			& )*recip_rAs(i,j,bi,bj)
144			& +tmpfldS )/tmpfldS
145			ELSE
146			rStarFacS(i,j,bi,bj) = 1.
147	jmc	1.2	ENDIF
148			ENDDO
149			ENDDO
150
151			C- Needs to do something when r* ratio is too small ;
152			C for now, just stop
153			DO j=1,sNy+1
154			DO i=1,sNx+1
155			IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
156			numbWrite = numbWrite + 1
157			WRITE(errorMessageUnit,'(2A,5I4,I10)')
158			& 'WARNING: r*FacC < hFacInf at:',
159			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
160			WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
161			& 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
162			& Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
163			WRITE(errorMessageUnit,'(A)')
164			& 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
165			STOP 'ABNORMAL END: S/R CALC_SURF_DR'
166			ENDIF
167			IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
168			numbWrite = numbWrite + 1
169			tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
170			& - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
171			WRITE(errorMessageUnit,'(2A,5I4,I10)')
172			& 'WARNING: r*FacW < hFacInf at:',
173			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
174			WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
175			& 'rStarFac,H,eta =', rStarFacW(i,j,bi,bj), tmpfldW,
176			& etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
177			WRITE(errorMessageUnit,'(A)')
178			& 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
179			STOP 'ABNORMAL END: S/R CALC_SURF_DR'
180			ENDIF
181			IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
182			numbWrite = numbWrite + 1
183			tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
184			& - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
185			WRITE(errorMessageUnit,'(2A,5I4,I10)')
186			& 'WARNING: r*FacS < hFacInf at:',
187			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
188			WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
189			& 'rStarFac,H,eta =', rStarFacS(i,j,bi,bj), tmpfldS,
190			& etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
191			WRITE(errorMessageUnit,'(A)')
192			& 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
193			STOP 'ABNORMAL END: S/R CALC_R_STAR'
194			ENDIF
195			C-- Usefull warning when r*Fac becomes very large:
196			IF ( numbWrite.LE.numbWrMax .AND.
197			& rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
198			numbWrite = numbWrite + 1
199			WRITE(errorMessageUnit,'(2A,5I4,I10)')
200			& 'WARNING: hFacC > hFacSup at:',
201			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
202			WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
203			& 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
204			& Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
205	jmc	1.1	ENDIF
206			ENDDO
207			ENDDO
208
209			C- end 1rst bi,bj loop.
210			ENDDO
211	jmc	1.3	ENDDO
212	jmc	1.1
213			_EXCH_XY_RL( rStarFacC, myThid )
214			CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
215
216			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
217
218			DO bj=myByLo(myThid), myByHi(myThid)
219			DO bi=myBxLo(myThid), myBxHi(myThid)
220			C- 2nd bi,bj loop :
221
222	jmc	1.4	#ifdef ALLOW_EXCH2
223			C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
224			C in the corner regions of the tile (e.g.:[1-Olx:0,1-Oly:0])
225			C => need to add those lines (or to fix exch2):
226			DO j=1,Oly
227			DO i=1,Olx
228			ii = sNx+i
229			jj = sNy+j
230
231			IF (maskH(1-i,1-j,bi,bj).EQ.0.) rStarFacC(1-i,1-j,bi,bj) = 1.
232			IF (maskH(ii, 1-j,bi,bj).EQ.0.) rStarFacC(ii, 1-j,bi,bj) = 1.
233			IF (maskH(1-i,jj, bi,bj).EQ.0.) rStarFacC(1-i,jj, bi,bj) = 1.
234			IF (maskH(ii, jj, bi,bj).EQ.0.) rStarFacC(ii, jj, bi,bj) = 1.
235
236			c IF (ksurfW(1-i,1-j,bi,bj).LE.Nr) rStarFacW(1-i,1-j,bi,bj)=1.
237			IF (maskW(1-i,1-j,km,bi,bj).EQ.0.) rStarFacW(1-i,1-j,bi,bj)=1.
238			IF (maskW(ii, 1-j,km,bi,bj).EQ.0.) rStarFacW(ii, 1-j,bi,bj)=1.
239			IF (maskW(1-i,jj, km,bi,bj).EQ.0.) rStarFacW(1-i,jj, bi,bj)=1.
240			IF (maskW(ii, jj, km,bi,bj).EQ.0.) rStarFacW(ii, jj, bi,bj)=1.
241
242			IF (maskS(1-i,1-j,km,bi,bj).EQ.0.) rStarFacS(1-i,1-j,bi,bj)=1.
243			IF (maskS(ii, 1-j,km,bi,bj).EQ.0.) rStarFacS(ii, 1-j,bi,bj)=1.
244			IF (maskS(1-i,jj, km,bi,bj).EQ.0.) rStarFacS(1-i,jj, bi,bj)=1.
245			IF (maskS(ii, jj, km,bi,bj).EQ.0.) rStarFacS(ii, jj, bi,bj)=1.
246			ENDDO
247			ENDDO
248			#endif /* ALLOW_EXCH2 */
249
250	jmc	1.1	DO j=1-Oly,sNy+Oly
251			DO i=1-Olx,sNx+Olx
252			rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
253			& -rStarExpC(i,j,bi,bj))/deltaTfreesurf
254			rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
255			& -rStarExpW(i,j,bi,bj))/deltaTfreesurf
256			rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
257			& -rStarExpS(i,j,bi,bj))/deltaTfreesurf
258			rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
259			& / rStarExpC(i,j,bi,bj)
260			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
261			& / rStarExpW(i,j,bi,bj)
262			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
263			& / rStarExpS(i,j,bi,bj)
264			ENDDO
265			ENDDO
266
267			C- end 2nd bi,bj loop.
268			ENDDO
269			ENDDO
270
271	jmc	1.4	#ifdef ALLOW_DEBUG
272			IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
273			#endif
274
275	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
276			#endif /* NONLIN_FRSURF */
277
278			RETURN
279			END