model/src/calc_r_star.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.6 2005/06/22 00:18:00 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 ii,jj
      INTEGER 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

      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_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
#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.7	C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.6 2005/06/22 00:18:00 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.4	INTEGER ii,jj
52	jmc	1.7	INTEGER numbWrite, numbWrMax
53	jmc	1.1	_RL tmpfldW, tmpfldS
54			c CHARACTER*(MAX_LEN_MBUF) suff
55			CEOP
56	jmc	1.2	DATA numbWrite / 0 /
57			numbWrMax = Nx*Ny
58	jmc	1.1
59			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
60
61	jmc	1.4	#ifdef ALLOW_DEBUG
62			IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
63			#endif
64
65	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
66			DO bi=myBxLo(myThid), myBxHi(myThid)
67			C- 1rst bi,bj loop :
68
69			IF (myIter.EQ.-1) THEN
70			C-- Initialise arrays :
71			DO j=1-Oly,sNy+Oly
72			DO i=1-Olx,sNx+Olx
73			rStarFacC(i,j,bi,bj) = 1.
74			rStarFacW(i,j,bi,bj) = 1.
75			rStarFacS(i,j,bi,bj) = 1.
76			rStarExpC(i,j,bi,bj) = 1.
77			rStarExpW(i,j,bi,bj) = 1.
78			rStarExpS(i,j,bi,bj) = 1.
79			rStarDhCDt(i,j,bi,bj) = 0.
80			rStarDhWDt(i,j,bi,bj) = 0.
81			rStarDhSDt(i,j,bi,bj) = 0.
82			PmEpR(i,j,bi,bj) = 0.
83			ENDDO
84			ENDDO
85			DO k=1,Nr
86			DO j=1-Oly,sNy+Oly
87			DO i=1-Olx,sNx+Olx
88			h0FacC(i,j,k,bi,bj) = hFacC(i,j,k,bi,bj)
89			h0FacW(i,j,k,bi,bj) = hFacW(i,j,k,bi,bj)
90			h0FacS(i,j,k,bi,bj) = hFacS(i,j,k,bi,bj)
91			ENDDO
92			ENDDO
93			ENDDO
94			ELSE
95			C-- copy rStarFacX -> rStarExpX
96			DO j=1-Oly,sNy+Oly
97			DO i=1-Olx,sNx+Olx
98			rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
99			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
100			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
101			ENDDO
102			ENDDO
103			ENDIF
104
105			C-- Compute the new column thikness :
106			DO j=0,sNy+1
107			DO i=0,sNx+1
108			IF (maskH(i,j,bi,bj).EQ.1. ) THEN
109			rStarFacC(i,j,bi,bj) =
110			& (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
111			& *recip_Rcol(i,j,bi,bj)
112			ELSE
113			rStarFacC(i,j,bi,bj) = 1.
114			ENDIF
115			ENDDO
116			ENDDO
117	jmc	1.2	DO j=1,sNy
118	jmc	1.1	DO i=1,sNx+1
119	jmc	1.7	IF ( ksurfW(i,j,bi,bj).LE.Nr ) THEN
120	jmc	1.1	tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
121			& - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
122			rStarFacW(i,j,bi,bj) =
123			& ( 0.5 _d 0 ( etaFld(i-1,j,bi,bj)rA(i-1,j,bi,bj)
124			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
125			& )*recip_rAw(i,j,bi,bj)
126			& +tmpfldW )/tmpfldW
127			ELSE
128			rStarFacW(i,j,bi,bj) = 1.
129			ENDIF
130	jmc	1.2	ENDDO
131			ENDDO
132			DO j=1,sNy+1
133			DO i=1,sNx
134	jmc	1.7	IF ( ksurfS(i,j,bi,bj).LE.Nr ) THEN
135	jmc	1.1	tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
136			& - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
137			rStarFacS(i,j,bi,bj) =
138			& ( 0.5 _d 0 ( etaFld(i,j-1,bi,bj)rA(i,j-1,bi,bj)
139			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
140			& )*recip_rAs(i,j,bi,bj)
141			& +tmpfldS )/tmpfldS
142			ELSE
143			rStarFacS(i,j,bi,bj) = 1.
144	jmc	1.2	ENDIF
145			ENDDO
146			ENDDO
147
148			C- Needs to do something when r* ratio is too small ;
149			C for now, just stop
150			DO j=1,sNy+1
151			DO i=1,sNx+1
152			IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
153			numbWrite = numbWrite + 1
154			WRITE(errorMessageUnit,'(2A,5I4,I10)')
155			& 'WARNING: r*FacC < hFacInf at:',
156			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
157			WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
158			& 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
159			& Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
160			WRITE(errorMessageUnit,'(A)')
161			& 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
162			STOP 'ABNORMAL END: S/R CALC_SURF_DR'
163			ENDIF
164			IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
165			numbWrite = numbWrite + 1
166			tmpfldW = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
167			& - MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
168			WRITE(errorMessageUnit,'(2A,5I4,I10)')
169			& 'WARNING: r*FacW < hFacInf at:',
170			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
171			WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
172			& 'rStarFac,H,eta =', rStarFacW(i,j,bi,bj), tmpfldW,
173			& etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
174			WRITE(errorMessageUnit,'(A)')
175			& 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
176			STOP 'ABNORMAL END: S/R CALC_SURF_DR'
177			ENDIF
178			IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
179			numbWrite = numbWrite + 1
180			tmpfldS = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
181			& - MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
182			WRITE(errorMessageUnit,'(2A,5I4,I10)')
183			& 'WARNING: r*FacS < hFacInf at:',
184			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
185			WRITE(errorMessageUnit,'(A,1F10.6,1P3E14.6)')
186			& 'rStarFac,H,eta =', rStarFacS(i,j,bi,bj), tmpfldS,
187			& etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
188			WRITE(errorMessageUnit,'(A)')
189			& 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
190			STOP 'ABNORMAL END: S/R CALC_R_STAR'
191			ENDIF
192			C-- Usefull warning when r*Fac becomes very large:
193			IF ( numbWrite.LE.numbWrMax .AND.
194			& rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
195			numbWrite = numbWrite + 1
196			WRITE(errorMessageUnit,'(2A,5I4,I10)')
197			& 'WARNING: hFacC > hFacSup at:',
198			& ' i,j,bi,bj,Thid,Iter=',i,j,bi,bj,myThid,myIter
199			WRITE(errorMessageUnit,'(A,1F10.6,1P2E14.6)')
200			& 'rStarFac,H,eta =', rStarFacC(i,j,bi,bj),
201			& Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj), etaFld(i,j,bi,bj)
202	jmc	1.1	ENDIF
203			ENDDO
204			ENDDO
205
206			C- end 1rst bi,bj loop.
207			ENDDO
208	jmc	1.3	ENDDO
209	jmc	1.1
210			_EXCH_XY_RL( rStarFacC, myThid )
211			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.7	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
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			DO i=1-Olx,sNx+Olx
251			rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
252			& -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