model/src/calc_r_star.F

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