model/src/calc_r_star.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.11 2007/05/22 11:34:18 mlosch 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,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 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 :

        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 ( ksurfW(i,j,bi,bj).LE.Nr ) 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 ( ksurfS(i,j,bi,bj).LE.Nr ) 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
        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 =
     &         MIN( Ro_surf(iiw-1,jjw,bi,bj),Ro_surf(iiw,jjw,bi,bj) )
     &       - MAX( R_low(iiw-1,jjw,bi,bj), R_low(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 =
     &         MIN( Ro_surf(iis,jjs-1,bi,bj),Ro_surf(iis,jjs,bi,bj) )
     &       - MAX( R_low(iis,jjs-1,bi,bj), R_low(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.12	C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.11 2007/05/22 11:34:18 mlosch 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			C i,j,k,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.1	INTEGER i,j,k,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			IF (myIter.EQ.-1) THEN
74			C-- Initialise arrays :
75			DO j=1-Oly,sNy+Oly
76	jmc	1.12	DO i=1-Olx,sNx+Olx
77	jmc	1.1	rStarFacC(i,j,bi,bj) = 1.
78			rStarFacW(i,j,bi,bj) = 1.
79			rStarFacS(i,j,bi,bj) = 1.
80			rStarExpC(i,j,bi,bj) = 1.
81			rStarExpW(i,j,bi,bj) = 1.
82			rStarExpS(i,j,bi,bj) = 1.
83			rStarDhCDt(i,j,bi,bj) = 0.
84			rStarDhWDt(i,j,bi,bj) = 0.
85			rStarDhSDt(i,j,bi,bj) = 0.
86			PmEpR(i,j,bi,bj) = 0.
87			ENDDO
88			ENDDO
89			DO k=1,Nr
90			DO j=1-Oly,sNy+Oly
91	jmc	1.12	DO i=1-Olx,sNx+Olx
92	heimbach	1.10	h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
93			h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
94			h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
95	jmc	1.1	ENDDO
96			ENDDO
97			ENDDO
98			ELSE
99			C-- copy rStarFacX -> rStarExpX
100			DO j=1-Oly,sNy+Oly
101	jmc	1.12	DO i=1-Olx,sNx+Olx
102	jmc	1.1	rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
103			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
104			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
105			ENDDO
106			ENDDO
107			ENDIF
108
109			C-- Compute the new column thikness :
110			DO j=0,sNy+1
111			DO i=0,sNx+1
112			IF (maskH(i,j,bi,bj).EQ.1. ) THEN
113	jmc	1.12	rStarFacC(i,j,bi,bj) =
114	jmc	1.1	& (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
115			& *recip_Rcol(i,j,bi,bj)
116			ELSE
117			rStarFacC(i,j,bi,bj) = 1.
118			ENDIF
119			ENDDO
120			ENDDO
121	jmc	1.2	DO j=1,sNy
122	jmc	1.1	DO i=1,sNx+1
123	jmc	1.7	IF ( ksurfW(i,j,bi,bj).LE.Nr ) THEN
124	jmc	1.1	tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
125			& - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
126	jmc	1.12	rStarFacW(i,j,bi,bj) =
127	jmc	1.1	& ( 0.5 _d 0 ( etaFld(i-1,j,bi,bj)rA(i-1,j,bi,bj)
128			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
129			& )*recip_rAw(i,j,bi,bj)
130			& +tmpfldW )/tmpfldW
131			ELSE
132			rStarFacW(i,j,bi,bj) = 1.
133			ENDIF
134	jmc	1.2	ENDDO
135			ENDDO
136			DO j=1,sNy+1
137			DO i=1,sNx
138	jmc	1.7	IF ( ksurfS(i,j,bi,bj).LE.Nr ) THEN
139	jmc	1.1	tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
140			& - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
141	jmc	1.12	rStarFacS(i,j,bi,bj) =
142	jmc	1.1	& ( 0.5 _d 0 ( etaFld(i,j-1,bi,bj)rA(i,j-1,bi,bj)
143			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
144			& )*recip_rAs(i,j,bi,bj)
145			& +tmpfldS )/tmpfldS
146			ELSE
147			rStarFacS(i,j,bi,bj) = 1.
148	jmc	1.2	ENDIF
149			ENDDO
150			ENDDO
151
152			C- Needs to do something when r* ratio is too small ;
153			C for now, just stop
154	mlosch	1.11	icntc1 = 0
155			icntc2 = 0
156			icntw = 0
157			icnts = 0
158	jmc	1.2	DO j=1,sNy+1
159			DO i=1,sNx+1
160	jmc	1.12	IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
161	mlosch	1.11	icntc1 = icntc1 + 1
162			iic = i
163			jjc = j
164	jmc	1.2	ENDIF
165	jmc	1.12	IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
166	mlosch	1.11	icntw = icntw + 1
167			iiw = i
168			jjw = j
169	jmc	1.2	ENDIF
170	jmc	1.12	IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
171	mlosch	1.11	icnts = icnts + 1
172			iis = i
173			jjs = j
174	jmc	1.2	ENDIF
175	mlosch	1.11	IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
176			icntc2 = icntc2 + 1
177			iic2 = i
178			jjc2 = j
179	jmc	1.1	ENDIF
180			ENDDO
181			ENDDO
182	jmc	1.12	IF ( icntc1+icnts+icntw .GT. 0 ) THEN
183			C- Print an error msg and then stop:
184			IF ( icntw .GT. 0 ) THEN
185			tmpfldW =
186			& MIN( Ro_surf(iiw-1,jjw,bi,bj),Ro_surf(iiw,jjw,bi,bj) )
187			& - MAX( R_low(iiw-1,jjw,bi,bj), R_low(iiw,jjw,bi,bj) )
188			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
189			& 'WARNING: r*FacW < hFacInf at',icntw,
190			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
191			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
192			& ' e.g. at i,j=',iiw,jjw,' ; rStarFac,H,eta =',
193			& rStarFacW(iiw,jjw,bi,bj), tmpfldW,
194			& etaFld(iiw-1,jjw,bi,bj), etaFld(iiw,jjw,bi,bj)
195			WRITE(errorMessageUnit,'(A)')
196			& 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
197			ENDIF
198			IF ( icnts .GT. 0 ) THEN
199			tmpfldS =
200			& MIN( Ro_surf(iis,jjs-1,bi,bj),Ro_surf(iis,jjs,bi,bj) )
201			& - MAX( R_low(iis,jjs-1,bi,bj), R_low(iis,jjs,bi,bj) )
202			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
203			& 'WARNING: r*FacS < hFacInf at',icnts,
204			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
205			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
206			& ' e.g. at i,j=',iis,jjs,' ; rStarFac,H,eta =',
207			& rStarFacS(iis,jjs,bi,bj), tmpfldS,
208			& etaFld(iis,jjs-1,bi,bj), etaFld(iis,jjs,bi,bj)
209			WRITE(errorMessageUnit,'(A)')
210			& 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
211			ENDIF
212			IF ( icntc1 .GT. 0 ) THEN
213			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
214			& 'WARNING: r*FacC < hFacInf at',icntc1,
215			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
216			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
217			& ' e.g. at i,j=',iic,jjc,' ; rStarFac,H,eta =',
218			& rStarFacC(iic,jjc,bi,bj),
219			& Ro_surf(iic,jjc,bi,bj)-R_low(iic,jjc,bi,bj),
220			& etaFld(iic,jjc,bi,bj)
221			WRITE(errorMessageUnit,'(A)')
222			& 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
223			ENDIF
224	mlosch	1.11	STOP 'ABNORMAL END: S/R CALC_R_STAR'
225	jmc	1.12	ENDIF
226			C-- Usefull warning when r*Fac becomes very large:
227			IF ( icntc2.GT.0 .AND. numbWrite.LE.numbWrMax ) THEN
228			numbWrite = numbWrite + 1
229			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
230			& 'WARNING: r*FacC > hFacSup at',icntc2,
231			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
232			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
233			& ' e.g. at i,j=',iic2,jjc2,' ; rStarFac,H,eta =',
234			& rStarFacC(iic2,jjc2,bi,bj),
235			& Ro_surf(iic2,jjc2,bi,bj)-R_low(iic2,jjc2,bi,bj),
236	mlosch	1.11	& etaFld(iic2,jjc2,bi,bj)
237	jmc	1.12	ENDIF
238	jmc	1.1
239			C- end 1rst bi,bj loop.
240			ENDDO
241	jmc	1.3	ENDDO
242	jmc	1.1
243	jmc	1.12	_EXCH_XY_RL( rStarFacC, myThid )
244	jmc	1.1	CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
245
246			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
247
248			DO bj=myByLo(myThid), myByHi(myThid)
249			DO bi=myBxLo(myThid), myBxHi(myThid)
250			C- 2nd bi,bj loop :
251
252	jmc	1.4	#ifdef ALLOW_EXCH2
253	jmc	1.7	#ifdef W2_FILL_NULL_REGIONS
254	jmc	1.4	C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
255			C in the corner regions of the tile (e.g.:[1-Olx:0,1-Oly:0])
256			C => need to add those lines (or to fix exch2):
257			DO j=1,Oly
258			DO i=1,Olx
259			ii = sNx+i
260			jj = sNy+j
261
262	jmc	1.7	IF (ksurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
263			IF (ksurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
264			IF (ksurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
265			IF (ksurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.
266
267			IF (ksurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
268			IF (ksurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
269			IF (ksurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
270			IF (ksurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.
271
272			IF (ksurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
273			IF (ksurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
274			IF (ksurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
275			IF (ksurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.
276
277	jmc	1.4	ENDDO
278			ENDDO
279	jmc	1.7	#endif /* W2_FILL_NULL_REGIONS */
280	jmc	1.4	#endif /* ALLOW_EXCH2 */
281
282	jmc	1.1	DO j=1-Oly,sNy+Oly
283	jmc	1.12	DO i=1-Olx,sNx+Olx
284			rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
285	jmc	1.1	& -rStarExpC(i,j,bi,bj))/deltaTfreesurf
286			rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
287			& -rStarExpW(i,j,bi,bj))/deltaTfreesurf
288			rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
289			& -rStarExpS(i,j,bi,bj))/deltaTfreesurf
290			rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
291			& / rStarExpC(i,j,bi,bj)
292			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
293			& / rStarExpW(i,j,bi,bj)
294			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
295			& / rStarExpS(i,j,bi,bj)
296			ENDDO
297			ENDDO
298
299			C- end 2nd bi,bj loop.
300			ENDDO
301			ENDDO
302
303	jmc	1.4	#ifdef ALLOW_DEBUG
304			IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
305			#endif
306
307	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
308			#endif /* NONLIN_FRSURF */
309
310			RETURN
311			END