model/src/calc_r_star.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.19 2011/01/14 01:28:31 gforget 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     etaFld :: current eta field used to update the hFactor
C     myTime :: Current time in simulation
C     myIter :: Current iteration number in simulation
C     myThid :: Thread number for this instance of the routine.
      _RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
      _RL myTime
      INTEGER myIter
      INTEGER myThid

#ifdef NONLIN_FRSURF

C     !LOCAL VARIABLES:
C     Local variables
C     rStarAreaWeight :: use area weighted average for rStarFac at U & V point
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
      LOGICAL rStarAreaWeight
      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
CEOP

#ifdef W2_FILL_NULL_REGIONS
      INTEGER ii,jj
#endif
      DATA numbWrite / 0 /
      numbWrMax = Nx*Ny

      rStarAreaWeight = .TRUE.
C-    Area-weighted average consistent with KE (& vert. advection):
      IF ( vectorInvariantMomentum .AND.
     &     (selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3)
     &   ) rStarAreaWeight =.FALSE.

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--   before updating rStarFacC/S/W save current fields
         DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
             rStarFacNm1C(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
             rStarFacNm1S(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
             rStarFacNm1W(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
           ENDDO
         ENDDO

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
       IF ( rStarAreaWeight ) THEN
C-     Area weighted average
        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
       ELSE
C-     Simple average
        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) + etaFld(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) + etaFld(i,j,bi,bj) )
     &        +tmpfldS )/tmpfldS
          ELSE
           rStarFacS(i,j,bi,bj) = 1.
          ENDIF
         ENDDO
        ENDDO
       ENDIF

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.20	C $Header: /u/gcmpack/MITgcm/model/src/calc_r_star.F,v 1.19 2011/01/14 01:28:31 gforget 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.20	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	jmc	1.20	C etaFld :: current eta field used to update the hFactor
35	jmc	1.1	C myTime :: Current time in simulation
36			C myIter :: Current iteration number in simulation
37			C myThid :: Thread number for this instance of the routine.
38	jmc	1.20	_RL etaFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
39	jmc	1.1	_RL myTime
40			INTEGER myIter
41			INTEGER myThid
42
43			#ifdef NONLIN_FRSURF
44
45			C !LOCAL VARIABLES:
46			C Local variables
47	jmc	1.17	C rStarAreaWeight :: use area weighted average for rStarFac at U & V point
48			C i,j, bi,bj :: loop counter
49			C numbWrite :: count the Number of warning written on STD-ERR file
50			C numbWrMax :: maximum Number of warning written on STD-ERR file
51			LOGICAL rStarAreaWeight
52	jmc	1.14	INTEGER i,j,bi,bj
53	jmc	1.7	INTEGER numbWrite, numbWrMax
54	mlosch	1.11	INTEGER icntc1, icntc2, icntw, icnts
55			INTEGER iic, jjc, iic2, jjc2, iiw, jjw, iis, jjs
56	jmc	1.1	_RL tmpfldW, tmpfldS
57			CEOP
58	jmc	1.17
59	jmc	1.20	#ifdef W2_FILL_NULL_REGIONS
60			INTEGER ii,jj
61			#endif
62			DATA numbWrite / 0 /
63			numbWrMax = Nx*Ny
64
65	jmc	1.17	rStarAreaWeight = .TRUE.
66			C- Area-weighted average consistent with KE (& vert. advection):
67	jmc	1.18	IF ( vectorInvariantMomentum .AND.
68			& (selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3)
69			& ) rStarAreaWeight =.FALSE.
70	jmc	1.17
71	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
72
73	jmc	1.4	#ifdef ALLOW_DEBUG
74			IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
75			#endif
76
77	jmc	1.1	DO bj=myByLo(myThid), myByHi(myThid)
78			DO bi=myBxLo(myThid), myBxHi(myThid)
79	jmc	1.20
80	gforget	1.19	C-- before updating rStarFacC/S/W save current fields
81			DO j=1-Oly,sNy+Oly
82			DO i=1-Olx,sNx+Olx
83			rStarFacNm1C(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
84			rStarFacNm1S(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
85			rStarFacNm1W(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
86			ENDDO
87			ENDDO
88	jmc	1.20
89	jmc	1.1	C- 1rst bi,bj loop :
90
91			C-- copy rStarFacX -> rStarExpX
92	jmc	1.13	DO j=1-Oly,sNy+Oly
93	jmc	1.12	DO i=1-Olx,sNx+Olx
94	jmc	1.1	rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
95			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
96			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
97			ENDDO
98	jmc	1.13	ENDDO
99	jmc	1.1
100			C-- Compute the new column thikness :
101			DO j=0,sNy+1
102			DO i=0,sNx+1
103	jmc	1.15	IF (kSurfC(i,j,bi,bj).LE.Nr ) THEN
104	jmc	1.12	rStarFacC(i,j,bi,bj) =
105	jmc	1.1	& (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
106			& *recip_Rcol(i,j,bi,bj)
107			ELSE
108			rStarFacC(i,j,bi,bj) = 1.
109			ENDIF
110			ENDDO
111			ENDDO
112	jmc	1.17	IF ( rStarAreaWeight ) THEN
113			C- Area weighted average
114	jmc	1.2	DO j=1,sNy
115	jmc	1.1	DO i=1,sNx+1
116	jmc	1.15	IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
117	jmc	1.16	tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
118	jmc	1.12	rStarFacW(i,j,bi,bj) =
119	jmc	1.1	& ( 0.5 _d 0 ( etaFld(i-1,j,bi,bj)rA(i-1,j,bi,bj)
120			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
121			& )*recip_rAw(i,j,bi,bj)
122			& +tmpfldW )/tmpfldW
123			ELSE
124			rStarFacW(i,j,bi,bj) = 1.
125			ENDIF
126	jmc	1.2	ENDDO
127			ENDDO
128			DO j=1,sNy+1
129			DO i=1,sNx
130	jmc	1.15	IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
131	jmc	1.16	tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
132	jmc	1.12	rStarFacS(i,j,bi,bj) =
133	jmc	1.1	& ( 0.5 _d 0 ( etaFld(i,j-1,bi,bj)rA(i,j-1,bi,bj)
134			& +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
135			& )*recip_rAs(i,j,bi,bj)
136			& +tmpfldS )/tmpfldS
137			ELSE
138			rStarFacS(i,j,bi,bj) = 1.
139	jmc	1.2	ENDIF
140			ENDDO
141			ENDDO
142	jmc	1.17	ELSE
143			C- Simple average
144			DO j=1,sNy
145			DO i=1,sNx+1
146			IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
147			tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
148			rStarFacW(i,j,bi,bj) =
149			& ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj) + etaFld(i,j,bi,bj) )
150			& +tmpfldW )/tmpfldW
151			ELSE
152			rStarFacW(i,j,bi,bj) = 1.
153			ENDIF
154			ENDDO
155			ENDDO
156			DO j=1,sNy+1
157			DO i=1,sNx
158			IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
159			tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
160			rStarFacS(i,j,bi,bj) =
161			& ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj) + etaFld(i,j,bi,bj) )
162			& +tmpfldS )/tmpfldS
163			ELSE
164			rStarFacS(i,j,bi,bj) = 1.
165			ENDIF
166			ENDDO
167			ENDDO
168			ENDIF
169	jmc	1.2
170			C- Needs to do something when r* ratio is too small ;
171			C for now, just stop
172	mlosch	1.11	icntc1 = 0
173			icntc2 = 0
174			icntw = 0
175			icnts = 0
176	jmc	1.2	DO j=1,sNy+1
177			DO i=1,sNx+1
178	jmc	1.12	IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
179	mlosch	1.11	icntc1 = icntc1 + 1
180			iic = i
181			jjc = j
182	jmc	1.2	ENDIF
183	jmc	1.12	IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
184	mlosch	1.11	icntw = icntw + 1
185			iiw = i
186			jjw = j
187	jmc	1.2	ENDIF
188	jmc	1.12	IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
189	mlosch	1.11	icnts = icnts + 1
190			iis = i
191			jjs = j
192	jmc	1.2	ENDIF
193	mlosch	1.11	IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
194			icntc2 = icntc2 + 1
195			iic2 = i
196			jjc2 = j
197	jmc	1.1	ENDIF
198			ENDDO
199			ENDDO
200	jmc	1.12	IF ( icntc1+icnts+icntw .GT. 0 ) THEN
201			C- Print an error msg and then stop:
202			IF ( icntw .GT. 0 ) THEN
203	jmc	1.16	tmpfldW = rSurfW(iiw,jjw,bi,bj) - rLowW(iiw,jjw,bi,bj)
204	jmc	1.12	WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
205			& 'WARNING: r*FacW < hFacInf at',icntw,
206			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
207			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
208			& ' e.g. at i,j=',iiw,jjw,' ; rStarFac,H,eta =',
209			& rStarFacW(iiw,jjw,bi,bj), tmpfldW,
210			& etaFld(iiw-1,jjw,bi,bj), etaFld(iiw,jjw,bi,bj)
211			WRITE(errorMessageUnit,'(A)')
212			& 'STOP in CALC_R_STAR : too SMALL rStarFacW !'
213			ENDIF
214			IF ( icnts .GT. 0 ) THEN
215	jmc	1.16	tmpfldS = rSurfS(iis,jjs,bi,bj) - rLowS(iis,jjs,bi,bj)
216	jmc	1.12	WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
217			& 'WARNING: r*FacS < hFacInf at',icnts,
218			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
219			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
220			& ' e.g. at i,j=',iis,jjs,' ; rStarFac,H,eta =',
221			& rStarFacS(iis,jjs,bi,bj), tmpfldS,
222			& etaFld(iis,jjs-1,bi,bj), etaFld(iis,jjs,bi,bj)
223			WRITE(errorMessageUnit,'(A)')
224			& 'STOP in CALC_R_STAR : too SMALL rStarFacS !'
225			ENDIF
226			IF ( icntc1 .GT. 0 ) THEN
227			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
228			& 'WARNING: r*FacC < hFacInf at',icntc1,
229			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
230			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
231			& ' e.g. at i,j=',iic,jjc,' ; rStarFac,H,eta =',
232			& rStarFacC(iic,jjc,bi,bj),
233			& Ro_surf(iic,jjc,bi,bj)-R_low(iic,jjc,bi,bj),
234			& etaFld(iic,jjc,bi,bj)
235			WRITE(errorMessageUnit,'(A)')
236			& 'STOP in CALC_R_STAR : too SMALL rStarFacC !'
237			ENDIF
238	mlosch	1.11	STOP 'ABNORMAL END: S/R CALC_R_STAR'
239	jmc	1.12	ENDIF
240			C-- Usefull warning when r*Fac becomes very large:
241			IF ( icntc2.GT.0 .AND. numbWrite.LE.numbWrMax ) THEN
242			numbWrite = numbWrite + 1
243			WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
244			& 'WARNING: r*FacC > hFacSup at',icntc2,
245			& ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
246			WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
247			& ' e.g. at i,j=',iic2,jjc2,' ; rStarFac,H,eta =',
248			& rStarFacC(iic2,jjc2,bi,bj),
249			& Ro_surf(iic2,jjc2,bi,bj)-R_low(iic2,jjc2,bi,bj),
250	mlosch	1.11	& etaFld(iic2,jjc2,bi,bj)
251	jmc	1.12	ENDIF
252	jmc	1.1
253			C- end 1rst bi,bj loop.
254			ENDDO
255	jmc	1.3	ENDDO
256	jmc	1.1
257	jmc	1.12	_EXCH_XY_RL( rStarFacC, myThid )
258	jmc	1.1	CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
259
260			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
261
262			DO bj=myByLo(myThid), myByHi(myThid)
263			DO bi=myBxLo(myThid), myBxHi(myThid)
264			C- 2nd bi,bj loop :
265
266	jmc	1.4	#ifdef ALLOW_EXCH2
267	jmc	1.7	#ifdef W2_FILL_NULL_REGIONS
268	jmc	1.4	C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
269			C in the corner regions of the tile (e.g.:[1-Olx:0,1-Oly:0])
270			C => need to add those lines (or to fix exch2):
271			DO j=1,Oly
272			DO i=1,Olx
273			ii = sNx+i
274			jj = sNy+j
275
276	jmc	1.15	IF (kSurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
277			IF (kSurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
278			IF (kSurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
279			IF (kSurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.
280
281			IF (kSurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
282			IF (kSurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
283			IF (kSurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
284			IF (kSurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.
285
286			IF (kSurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
287			IF (kSurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
288			IF (kSurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
289			IF (kSurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.
290	jmc	1.7
291	jmc	1.4	ENDDO
292			ENDDO
293	jmc	1.7	#endif /* W2_FILL_NULL_REGIONS */
294	jmc	1.4	#endif /* ALLOW_EXCH2 */
295
296	jmc	1.1	DO j=1-Oly,sNy+Oly
297	jmc	1.12	DO i=1-Olx,sNx+Olx
298			rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
299	jmc	1.1	& -rStarExpC(i,j,bi,bj))/deltaTfreesurf
300			rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
301			& -rStarExpW(i,j,bi,bj))/deltaTfreesurf
302			rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
303			& -rStarExpS(i,j,bi,bj))/deltaTfreesurf
304			rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
305			& / rStarExpC(i,j,bi,bj)
306			rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
307			& / rStarExpW(i,j,bi,bj)
308			rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
309			& / rStarExpS(i,j,bi,bj)
310			ENDDO
311			ENDDO
312
313			C- end 2nd bi,bj loop.
314			ENDDO
315			ENDDO
316
317	jmc	1.4	#ifdef ALLOW_DEBUG
318			IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
319			#endif
320
321	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
322			#endif /* NONLIN_FRSURF */
323
324			RETURN
325			END