cg/code/set_ddtvars.F

C $Header: /u/gcmpack/MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F,v 1.6 2008/05/04 20:39:32 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#ifdef ALLOW_MYPACKAGE
# include "MYPACKAGE_OPTIONS.h"
#else
# include "CPP_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: SET_DDTVARS
C     !INTERFACE:
      SUBROUTINE SET_DDTVARS( gUVelC, gVVelC, gThetaC, gSaltC,
     I                        cnx, cny, cnr,
     I                        myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SET_DDTVARS
C     | o Set tendencies to be used as additional forcing
C     *==========================================================*
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "ORIENTATION.h"
#include "EXT_MOM_TEND.h"
#ifdef ALLOW_MYPACKAGE
# include "MYPACKAGE.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime    :: Current time in simulation
C     myIter    :: Current iteration number
C     myThid    :: Thread Id number
      INTEGER cnx, cny, cnr
      REAL*8 gUVelC( cnx, cny, cnr)
      REAL*8 gVVelC( cnx, cny, cnr)
      REAL*8 gThetaC(cnx, cny, cnr)
      REAL*8 gSaltC( cnx, cny, cnr)
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C     !FUNCTIONS:
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !LOCAL VARIABLES:
C     == Local variables ==
C     i,j,k     :: Loop counters
      INTEGER i,j,k
      INTEGER bi,bj
C     Variables used for I/O
      CHARACTER*(10) suff
C     Variables used for smoothing momentum tendencies:
      _RL  viscAhDivDt, viscAhVortDt
      _RL  locDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  locVort(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  locFx  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  locFy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

      viscAhDivDt  = myPa_param1 * deltaTmom
      viscAhVortDt = myPa_param2 * deltaTmom

      IF ( cnx .NE. sNx ) THEN
       STOP 'SET_DTTVARS cnx NE sNx'
      ENDIF
      IF ( cny .NE. sNy ) THEN
       STOP 'SET_DTTVARS cny NE sNy'
      ENDIF
      IF ( cnr .NE.  Nr ) THEN
       STOP 'SET_DTTVARS cnr NE  Nr'
      ENDIF

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C---  Temp. & Salt Tendency from CG : just do a copy
        DO k=1,Nr
         DO j=1,sNy
          DO i=1,sNx
           myPa_TendScal1(i,j,k,bi,bj) = gThetaC(i,j,k)
           myPa_TendScal2(i,j,k,bi,bj) = gsaltC (i,j,k)
          ENDDO
         ENDDO
        ENDDO

C---  Momentum Tendency from CG
C--   1) rotate to get u,v aline with CG axes
        DO k=1,Nr
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            ext_gu(i,j,k,bi,bj) = 0. _d 0
            ext_gv(i,j,k,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
         IF ( myPa_doSwitch1 ) THEN
          DO j=1,sNy
           DO i=1,sNx
            IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN
             ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj)
             ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj)
            ENDIF
           ENDDO
          ENDDO
         ELSE
          DO j=1,sNy
           DO i=1,sNx
            IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN
             ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj)
     &                           - gVVelC(i,j,k)*sAngleFG(i,j,bi,bj)
             ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj)
     &                           + gVVelC(i,j,k)*cAngleFG(i,j,bi,bj)
            ENDIF
           ENDDO
          ENDDO
         ENDIF
        ENDDO

C-    end bi,bj loops
       ENDDO
      ENDDO

C--   2) Fill the overlap : ext_gu & ext_gv are both on A-grid
      CALL EXCH_UV_AGRID_3D_RL( ext_gu, ext_gv,
     &                          .TRUE., Nr, myThid )

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
C--   3) average to C-grid
         DO j=1-OLy,sNy+OLy
          DO i=2-OLx,sNx+OLx
            myPa_TendVelU(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
     &          *( ext_gu(i-1,j,k,bi,bj)
     &            +ext_gu( i ,j,k,bi,bj) )*0.5 _d 0
          ENDDO
         ENDDO
         DO j=2-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
            myPa_TendVelV(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
     &          *( ext_gv(i,j-1,k,bi,bj)
     &            +ext_gv(i, j ,k,bi,bj) )*0.5 _d 0
          ENDDO
         ENDDO
C--   4) Smooth Momentum tendency :
C       apply horizontal viscosity (*time-step = viscAhDivDt) on Divergence
C       tendency and on vorticity tend. (viscosity*time-step = viscAhVortDt)
         IF ( viscAhDivDt.GT.0. .OR. viscAhVortDt.GT.0. ) THEN
C     Fluxes
          DO j=2-OLy,sNy+OLy
           DO i=2-OLx,sNx+OLx
             locFx(i,j) = myPa_TendVelU(i,j,k,bi,bj)
     &                  *dyG(i,j,bi,bj)*hFacW(i,j,k,bi,bj)
             locFy(i,j) = myPa_TendVelV(i,j,k,bi,bj)
     &                  *dxG(i,j,bi,bj)*hFacS(i,j,k,bi,bj)
           ENDDO
          ENDDO
C     Divergence
          DO j=2-OLy,sNy+OLy-1
           DO i=2-OLx,sNx+OLx-1
            locDiv(i,j) =
     &                 (  ( locFx(i+1,j)-locFx(i,j) )
     &                   +( locFy(i,j+1)-locFy(i,j) )
     &                 )*recip_rA(i,j,bi,bj)*recip_hFacC(i,j,k,bi,bj)
           ENDDO
          ENDDO
C     Vorticity
          DO j=2-OLy,sNy+OLy
           DO i=2-OLx,sNx+OLx
            locVort(i,j) =
     &            ( ( myPa_TendVelV(i,j,k,bi,bj)*dyC(i,j,bi,bj)
     &               -myPa_TendVelV(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj) )
     &             -( myPa_TendVelU(i,j,k,bi,bj)*dxC(i,j,bi,bj)
     &               -myPa_TendVelU(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj) )
     &            )*recip_rAz(I,J,bi,bj)
           ENDDO
          ENDDO
C     Apply to C-grid tendencies:
          DO j=3-OLy,sNy+OLy-1
           DO i=3-OLx,sNx+OLx-1
            myPa_TendVelU(i,j,k,bi,bj) = myPa_TendVelU(i,j,k,bi,bj)
     &        + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i-1,j) )
     &                        * recip_dxC(i,j,bi,bj)
     &          - viscAhVortDt* ( locVort(i,j+1)-locVort(i,j) )
     &                        * recip_dyG(i,j,bi,bj)
     &          )*maskW(i,j,k,bi,bj)
            myPa_TendVelV(i,j,k,bi,bj) = myPa_TendVelV(i,j,k,bi,bj)
     &        + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i,j-1) )
     &                        * recip_dyC(i,j,bi,bj)
     &          + viscAhVortDt* ( locVort(i+1,j)-locVort(i,j) )
     &                        * recip_dxG(i,j,bi,bj)
     &          )*maskS(i,j,k,bi,bj)
           ENDDO
          ENDDO
C-    end smoothing & end of k loop
         ENDIF
        ENDDO

C-    end bi,bj loops
       ENDDO
      ENDDO

C--   Diagnostics
#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        CALL DIAGNOSTICS_FILL ( myPa_TendScal1,'MYPadTdt',
     &                          0, Nr, 0, 1, 1, myThid )
        CALL DIAGNOSTICS_FILL ( myPa_TendScal2,'MYPadSdt',
     &                          0, Nr, 0, 1, 1, myThid )
        CALL DIAGNOSTICS_FILL ( myPa_TendVelU, 'MYPadUdt',
     &                          0, Nr, 0, 1, 1, myThid )
        CALL DIAGNOSTICS_FILL ( myPa_TendVelV, 'MYPadVdt',
     &                          0, Nr, 0, 1, 1, myThid )
C-   This is a hack (to avoid changing mypackage_diagnostics_init.F)
        CALL DIAGNOSTICS_FILL ( ext_gu, 'MYPaStaU',
     &                          0, Nr, 0, 1, 1, myThid )
        CALL DIAGNOSTICS_FILL ( ext_gv, 'MYPaStaV',
     &                          0, Nr, 0, 1, 1, myThid )
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C--   Write snap-shot output:
      IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime,deltaTClock)
     &     .OR. dumpInitAndLast.AND.( myTime.EQ.endTime .OR.
     &                                myTime.EQ.startTime  )
     &   ) THEN
        WRITE(suff,'(I10.10)') myIter
        CALL WRITE_FLD_XYZ_RL( 'ext_gT.', suff, myPa_TendScal1,
     &                         myIter, myThid )
        CALL WRITE_FLD_XYZ_RL( 'ext_gS.', suff, myPa_TendScal2,
     &                         myIter, myThid )
        CALL WRITE_FLD_XYZ_RL( 'ext_gU.', suff, ext_gu,
     &                         myIter, myThid )
        CALL WRITE_FLD_XYZ_RL( 'ext_gV.', suff, ext_gv,
     &                         myIter, myThid )
      ENDIF

      RETURN
      END
1	jmc	1.7	C $Header: /u/gcmpack/MITgcm_contrib/PRM/multi_comp_setup/cg/code/set_ddtvars.F,v 1.6 2008/05/04 20:39:32 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5	jmc	1.5	#ifdef ALLOW_MYPACKAGE
6	jmc	1.3	# include "MYPACKAGE_OPTIONS.h"
7			#else
8			# include "CPP_OPTIONS.h"
9			#endif
10	jmc	1.1
11			CBOP
12			C !ROUTINE: SET_DDTVARS
13			C !INTERFACE:
14	jmc	1.6	SUBROUTINE SET_DDTVARS( gUVelC, gVVelC, gThetaC, gSaltC,
15	jmc	1.5	I cnx, cny, cnr,
16	jmc	1.1	I myTime, myIter, myThid )
17			C !DESCRIPTION: \bv
18			C ==========================================================
19			C \| S/R SET_DDTVARS
20	jmc	1.5	C \| o Set tendencies to be used as additional forcing
21	jmc	1.1	C ==========================================================
22			C ==========================================================
23			C \ev
24
25			C !USES:
26			IMPLICIT NONE
27
28			C == Global data ==
29			#include "SIZE.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32	jmc	1.5	#include "GRID.h"
33	jmc	1.1	#include "ORIENTATION.h"
34	jmc	1.5	#include "EXT_MOM_TEND.h"
35			#ifdef ALLOW_MYPACKAGE
36	jmc	1.3	# include "MYPACKAGE.h"
37			#endif
38	jmc	1.1
39			C !INPUT/OUTPUT PARAMETERS:
40			C == Routine arguments ==
41			C myTime :: Current time in simulation
42			C myIter :: Current iteration number
43			C myThid :: Thread Id number
44			INTEGER cnx, cny, cnr
45			REAL*8 gUVelC( cnx, cny, cnr)
46			REAL*8 gVVelC( cnx, cny, cnr)
47			REAL*8 gThetaC(cnx, cny, cnr)
48	jmc	1.6	REAL*8 gSaltC( cnx, cny, cnr)
49	jmc	1.1	_RL myTime
50			INTEGER myIter
51			INTEGER myThid
52
53			C !FUNCTIONS:
54			LOGICAL DIFFERENT_MULTIPLE
55			EXTERNAL DIFFERENT_MULTIPLE
56
57			C !LOCAL VARIABLES:
58			C == Local variables ==
59			C i,j,k :: Loop counters
60			INTEGER i,j,k
61	jmc	1.5	INTEGER bi,bj
62	jmc	1.1	C Variables used for I/O
63			CHARACTER*(10) suff
64	jmc	1.5	C Variables used for smoothing momentum tendencies:
65			_RL viscAhDivDt, viscAhVortDt
66			_RL locDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67			_RL locVort(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68			_RL locFx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69			_RL locFy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70	jmc	1.1	CEOP
71
72	jmc	1.5	viscAhDivDt = myPa_param1 * deltaTmom
73			viscAhVortDt = myPa_param2 * deltaTmom
74
75	jmc	1.1	IF ( cnx .NE. sNx ) THEN
76			STOP 'SET_DTTVARS cnx NE sNx'
77			ENDIF
78			IF ( cny .NE. sNy ) THEN
79			STOP 'SET_DTTVARS cny NE sNy'
80			ENDIF
81			IF ( cnr .NE. Nr ) THEN
82			STOP 'SET_DTTVARS cnr NE Nr'
83			ENDIF
84
85	jmc	1.5	DO bj = myByLo(myThid), myByHi(myThid)
86			DO bi = myBxLo(myThid), myBxHi(myThid)
87
88	jmc	1.6	C--- Temp. & Salt Tendency from CG : just do a copy
89	jmc	1.5	DO k=1,Nr
90			DO j=1,sNy
91			DO i=1,sNx
92			myPa_TendScal1(i,j,k,bi,bj) = gThetaC(i,j,k)
93	jmc	1.6	myPa_TendScal2(i,j,k,bi,bj) = gsaltC (i,j,k)
94	jmc	1.5	ENDDO
95			ENDDO
96			ENDDO
97
98			C--- Momentum Tendency from CG
99			C-- 1) rotate to get u,v aline with CG axes
100			DO k=1,Nr
101			DO j=1-Oly,sNy+Oly
102			DO i=1-Olx,sNx+Olx
103			ext_gu(i,j,k,bi,bj) = 0. _d 0
104			ext_gv(i,j,k,bi,bj) = 0. _d 0
105			ENDDO
106			ENDDO
107	jmc	1.7	IF ( myPa_doSwitch1 ) THEN
108			DO j=1,sNy
109			DO i=1,sNx
110			IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN
111			ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj)
112			ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj)
113			ENDIF
114			ENDDO
115			ENDDO
116			ELSE
117			DO j=1,sNy
118			DO i=1,sNx
119			IF ( velDvRMS(i,j,bi,bj).GT.0. _d 0 ) THEN
120	jmc	1.5	ext_gu(i,j,k,bi,bj) = gUVelC(i,j,k)*cAngleFG(i,j,bi,bj)
121			& - gVVelC(i,j,k)*sAngleFG(i,j,bi,bj)
122	jmc	1.7	ext_gv(i,j,k,bi,bj) = gUVelC(i,j,k)*sAngleFG(i,j,bi,bj)
123			& + gVVelC(i,j,k)*cAngleFG(i,j,bi,bj)
124			ENDIF
125			ENDDO
126	jmc	1.5	ENDDO
127	jmc	1.7	ENDIF
128	jmc	1.1	ENDDO
129	jmc	1.5
130			C- end bi,bj loops
131	jmc	1.1	ENDDO
132			ENDDO
133
134	jmc	1.5	C-- 2) Fill the overlap : ext_gu & ext_gv are both on A-grid
135			CALL EXCH_UV_AGRID_3D_RL( ext_gu, ext_gv,
136			& .TRUE., Nr, myThid )
137
138			DO bj = myByLo(myThid), myByHi(myThid)
139			DO bi = myBxLo(myThid), myBxHi(myThid)
140			DO k=1,Nr
141			C-- 3) average to C-grid
142			DO j=1-OLy,sNy+OLy
143			DO i=2-OLx,sNx+OLx
144			myPa_TendVelU(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
145			& *( ext_gu(i-1,j,k,bi,bj)
146			& +ext_gu( i ,j,k,bi,bj) )*0.5 _d 0
147			ENDDO
148			ENDDO
149			DO j=2-OLy,sNy+OLy
150			DO i=1-OLx,sNx+OLx
151			myPa_TendVelV(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
152			& *( ext_gv(i,j-1,k,bi,bj)
153			& +ext_gv(i, j ,k,bi,bj) )*0.5 _d 0
154			ENDDO
155			ENDDO
156			C-- 4) Smooth Momentum tendency :
157			C apply horizontal viscosity (*time-step = viscAhDivDt) on Divergence
158			C tendency and on vorticity tend. (viscosity*time-step = viscAhVortDt)
159			IF ( viscAhDivDt.GT.0. .OR. viscAhVortDt.GT.0. ) THEN
160			C Fluxes
161			DO j=2-OLy,sNy+OLy
162			DO i=2-OLx,sNx+OLx
163			locFx(i,j) = myPa_TendVelU(i,j,k,bi,bj)
164			& dyG(i,j,bi,bj)hFacW(i,j,k,bi,bj)
165			locFy(i,j) = myPa_TendVelV(i,j,k,bi,bj)
166			& dxG(i,j,bi,bj)hFacS(i,j,k,bi,bj)
167			ENDDO
168			ENDDO
169			C Divergence
170			DO j=2-OLy,sNy+OLy-1
171			DO i=2-OLx,sNx+OLx-1
172			locDiv(i,j) =
173			& ( ( locFx(i+1,j)-locFx(i,j) )
174			& +( locFy(i,j+1)-locFy(i,j) )
175			& )recip_rA(i,j,bi,bj)recip_hFacC(i,j,k,bi,bj)
176			ENDDO
177			ENDDO
178			C Vorticity
179			DO j=2-OLy,sNy+OLy
180			DO i=2-OLx,sNx+OLx
181			locVort(i,j) =
182			& ( ( myPa_TendVelV(i,j,k,bi,bj)*dyC(i,j,bi,bj)
183			& -myPa_TendVelV(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj) )
184			& -( myPa_TendVelU(i,j,k,bi,bj)*dxC(i,j,bi,bj)
185			& -myPa_TendVelU(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj) )
186			& )*recip_rAz(I,J,bi,bj)
187			ENDDO
188			ENDDO
189			C Apply to C-grid tendencies:
190			DO j=3-OLy,sNy+OLy-1
191			DO i=3-OLx,sNx+OLx-1
192			myPa_TendVelU(i,j,k,bi,bj) = myPa_TendVelU(i,j,k,bi,bj)
193			& + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i-1,j) )
194			& * recip_dxC(i,j,bi,bj)
195			& - viscAhVortDt* ( locVort(i,j+1)-locVort(i,j) )
196			& * recip_dyG(i,j,bi,bj)
197			& )*maskW(i,j,k,bi,bj)
198			myPa_TendVelV(i,j,k,bi,bj) = myPa_TendVelV(i,j,k,bi,bj)
199			& + ( viscAhDivDt * ( locDiv(i,j) - locDiv(i,j-1) )
200			& * recip_dyC(i,j,bi,bj)
201			& + viscAhVortDt* ( locVort(i+1,j)-locVort(i,j) )
202			& * recip_dxG(i,j,bi,bj)
203			& )*maskS(i,j,k,bi,bj)
204			ENDDO
205			ENDDO
206			C- end smoothing & end of k loop
207			ENDIF
208	jmc	1.1	ENDDO
209	jmc	1.5
210			C- end bi,bj loops
211	jmc	1.1	ENDDO
212			ENDDO
213
214	jmc	1.5	C-- Diagnostics
215	jmc	1.4	#ifdef ALLOW_DIAGNOSTICS
216			IF ( useDiagnostics ) THEN
217	jmc	1.5	CALL DIAGNOSTICS_FILL ( myPa_TendScal1,'MYPadTdt',
218			& 0, Nr, 0, 1, 1, myThid )
219	jmc	1.6	CALL DIAGNOSTICS_FILL ( myPa_TendScal2,'MYPadSdt',
220			& 0, Nr, 0, 1, 1, myThid )
221	jmc	1.4	CALL DIAGNOSTICS_FILL ( myPa_TendVelU, 'MYPadUdt',
222			& 0, Nr, 0, 1, 1, myThid )
223			CALL DIAGNOSTICS_FILL ( myPa_TendVelV, 'MYPadVdt',
224			& 0, Nr, 0, 1, 1, myThid )
225	jmc	1.5	C- This is a hack (to avoid changing mypackage_diagnostics_init.F)
226			CALL DIAGNOSTICS_FILL ( ext_gu, 'MYPaStaU',
227			& 0, Nr, 0, 1, 1, myThid )
228			CALL DIAGNOSTICS_FILL ( ext_gv, 'MYPaStaV',
229	jmc	1.4	& 0, Nr, 0, 1, 1, myThid )
230			ENDIF
231			#endif /* ALLOW_DIAGNOSTICS */
232
233	jmc	1.5	C-- Write snap-shot output:
234	jmc	1.1	IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime,deltaTClock)
235			& .OR. dumpInitAndLast.AND.( myTime.EQ.endTime .OR.
236			& myTime.EQ.startTime )
237			& ) THEN
238			WRITE(suff,'(I10.10)') myIter
239	jmc	1.5	CALL WRITE_FLD_XYZ_RL( 'ext_gT.', suff, myPa_TendScal1,
240	jmc	1.3	& myIter, myThid )
241	jmc	1.6	CALL WRITE_FLD_XYZ_RL( 'ext_gS.', suff, myPa_TendScal2,
242			& myIter, myThid )
243	jmc	1.5	CALL WRITE_FLD_XYZ_RL( 'ext_gU.', suff, ext_gu,
244	jmc	1.3	& myIter, myThid )
245	jmc	1.5	CALL WRITE_FLD_XYZ_RL( 'ext_gV.', suff, ext_gv,
246	jmc	1.3	& myIter, myThid )
247	jmc	1.1	ENDIF
248
249			RETURN
250			END