model/src/calc_r_star.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.15 2009/12/21 00:21:35 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_EXCH2
# include "W2_OPTIONS.h"
#endif

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, 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,bi,bj
      INTEGER numbWrite, numbWrMax
      INTEGER icntc1, icntc2, icntw, icnts
      INTEGER iic, jjc, iic2, jjc2, iiw, jjw, iis, jjs
      _RL tmpfldW, tmpfldS
c     CHARACTER*(MAX_LEN_MBUF) suff
CEOP
#ifdef W2_FILL_NULL_REGIONS
      INTEGER ii,jj
#endif
      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

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

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

C-- Compute the new column thikness :
        DO j=0,sNy+1
         DO i=0,sNx+1
          IF (kSurfC(i,j,bi,bj).LE.Nr ) 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 ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
           tmpfldW = rSurfW(i,j,bi,bj) - rLowW(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 ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
           tmpfldS = rSurfS(i,j,bi,bj) - rLowS(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
        icntc1 = 0
        icntc2 = 0
        icntw  = 0
        icnts  = 0
        DO j=1,sNy+1
         DO i=1,sNx+1
          IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
            icntc1 = icntc1 + 1
            iic = i
            jjc = j
          ENDIF
          IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
            icntw = icntw + 1
            iiw = i
            jjw = j
          ENDIF
          IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
            icnts = icnts + 1
            iis = i
            jjs = j
          ENDIF
          IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
            icntc2 = icntc2 + 1
            iic2 = i
            jjc2 = j
          ENDIF
         ENDDO
        ENDDO
        IF ( icntc1+icnts+icntw .GT. 0 ) THEN
C-    Print an error msg and then stop:
         IF ( icntw  .GT. 0 ) THEN
          tmpfldW = rSurfW(iiw,jjw,bi,bj) - rLowW(iiw,jjw,bi,bj)
          WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
     &     'WARNING: r*FacW < hFacInf at',icntw,
     &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
          WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
     &     ' e.g. at i,j=',iiw,jjw,' ; rStarFac,H,eta =',
     &      rStarFacW(iiw,jjw,bi,bj), tmpfldW,
     &      etaFld(iiw-1,jjw,bi,bj), etaFld(iiw,jjw,bi,bj)
          WRITE(errorMessageUnit,'(A)')
     &     'STOP in CALC_R_STAR : too SMALL rStarFacW !'
         ENDIF
         IF ( icnts  .GT. 0 ) THEN
          tmpfldS = rSurfS(iis,jjs,bi,bj) - rLowS(iis,jjs,bi,bj)
          WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
     &     'WARNING: r*FacS < hFacInf at',icnts,
     &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
          WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
     &     ' e.g. at i,j=',iis,jjs,' ; rStarFac,H,eta =',
     &      rStarFacS(iis,jjs,bi,bj), tmpfldS,
     &      etaFld(iis,jjs-1,bi,bj), etaFld(iis,jjs,bi,bj)
          WRITE(errorMessageUnit,'(A)')
     &     'STOP in CALC_R_STAR : too SMALL rStarFacS !'
         ENDIF
         IF ( icntc1 .GT. 0 ) THEN
          WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
     &     'WARNING: r*FacC < hFacInf at',icntc1,
     &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
          WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
     &     ' e.g. at i,j=',iic,jjc,' ; rStarFac,H,eta =',
     &     rStarFacC(iic,jjc,bi,bj),
     &     Ro_surf(iic,jjc,bi,bj)-R_low(iic,jjc,bi,bj),
     &     etaFld(iic,jjc,bi,bj)
          WRITE(errorMessageUnit,'(A)')
     &     'STOP in CALC_R_STAR : too SMALL rStarFacC !'
         ENDIF
         STOP 'ABNORMAL END: S/R CALC_R_STAR'
        ENDIF
C-- Usefull warning when r*Fac becomes very large:
        IF ( icntc2.GT.0 .AND. numbWrite.LE.numbWrMax ) THEN
         numbWrite = numbWrite + 1
         WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
     &    'WARNING: r*FacC > hFacSup at',icntc2,
     &    ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
         WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
     &    ' e.g. at i,j=',iic2,jjc2,' ; rStarFac,H,eta =',
     &    rStarFacC(iic2,jjc2,bi,bj),
     &    Ro_surf(iic2,jjc2,bi,bj)-R_low(iic2,jjc2,bi,bj),
     &    etaFld(iic2,jjc2,bi,bj)
        ENDIF

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
#ifdef W2_FILL_NULL_REGIONS
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 (kSurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
          IF (kSurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
          IF (kSurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
          IF (kSurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.

          IF (kSurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
          IF (kSurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
          IF (kSurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
          IF (kSurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.

          IF (kSurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
          IF (kSurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
          IF (kSurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
          IF (kSurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.

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