pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.7 2006/04/09 17:35:30 heimbach Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"

CBOP
C     !ROUTINE: THSICE_MAIN
C     !INTERFACE:
      SUBROUTINE THSICE_MAIN(
     I                        myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R  THSICE_MAIN
C     | o Therm_SeaIce main routine.
C     |   step forward Thermodynamic_SeaIce variables and modify
C     |    ocean surface forcing accordingly.
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myIter :: iteration counter for this thread
C     myTime :: time counter for this thread
C     myThid :: thread number for this instance of the routine.
      _RL  myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_THSICE
C     !LOCAL VARIABLES:
C     === Local variables ===
      INTEGER i,j
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

      _RL tauFac

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( stressReduction.GT. 0. _d 0 ) THEN
C-     needs new Ice Fraction in halo region to apply wind-stress reduction
       iMin = 1-OLx
       iMax = sNx+OLx-1
       jMin = 1-OLy
       jMax = sNy+OLy-1
#ifdef ATMOSPHERIC_LOADING
      ELSEIF ( useRealFreshWaterFlux ) THEN
C-     needs sea-ice loading in part of the halo regions for grad.Phi0surf
C      to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
       iMin = 0
       iMax = sNx+1
       jMin = 0
       jMax = sNy+1
#endif /* ATMOSPHERIC_LOADING */
      ELSE
       iMin = 1
       iMax = sNx
       jMin = 1
       jMax = sNy
      ENDIF

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

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          iicekey = (act1 + 1) + act2*max1
     &                         + act3*max1*max2
     &                         + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

C--     Mixed layer thickness: take the 1rst layer
#ifdef NONLIN_FRSURF
        IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
         IF ( select_rStar.GT.0 ) THEN
          DO j = jMin, jMax
           DO i = iMin, iMax
             hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
     &                                  *rStarFacC(i,j,bi,bj)
           ENDDO
          ENDDO
         ELSE
          DO j = jMin, jMax
           DO i = iMin, iMax
            IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
             hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
            ELSE
             hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
            ENDIF
           ENDDO
          ENDDO
         ENDIF
        ELSE
#else /* ndef NONLIN_FRSURF */
        IF (.TRUE.) THEN
#endif /* NONLIN_FRSURF */
          DO j = jMin, jMax
           DO i = iMin, iMax
             hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
           ENDDO
          ENDDO
        ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
#endif

         DO j = jMin, jMax
          DO i = iMin, iMax
           tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
           sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
           v2ocMxL(i,j,bi,bj) =
     &              ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
     &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
     &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
     &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
     &              )*0.5 _d 0
           prcAtm(i,j) = 0.
           icFrwAtm(i,j,bi,bj) = 0. _d 0
           icFlxAtm(i,j,bi,bj) = 0. _d 0
           icFlxSW (i,j,bi,bj) = 0. _d 0
           snowPrc(i,j,bi,bj) = 0. _d 0
           siceAlb(i,j,bi,bj) = 0. _d 0
          ENDDO
         ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask = comlev1, key = iicekey
CADJ STORE iceHeight  = comlev1, key = iicekey
CADJ STORE snowHeight = comlev1, key = iicekey
CADJ STORE Tsrf    = comlev1, key = iicekey
CADJ STORE Qice1   = comlev1, key = iicekey
CADJ STORE Qice2   = comlev1, key = iicekey
CADJ STORE snowAge = comlev1, key = iicekey

CADJ STORE sHeating = comlev1, key = iicekey
CADJ STORE flxCndBt = comlev1, key = iicekey
CADJ STORE snowPrc  = comlev1, key = iicekey

CADJ STORE hOceMxL = comlev1, key = iicekey
CADJ STORE tOceMxL = comlev1, key = iicekey
CADJ STORE sOceMxL = comlev1, key = iicekey
CADJ STORE v2ocMxL = comlev1, key = iicekey

CADJ STORE empmr   = comlev1, key = iicekey
CADJ STORE qnet    = comlev1, key = iicekey
#endif

#ifdef ALLOW_BULK_FORCE
         IF ( useBulkforce ) THEN
           CALL THSICE_GET_PRECIP(
     I                  iceMask,
     O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
     O                  icFlxSW(1-OLx,1-OLy,bi,bj),
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
         ENDIF
#endif


         CALL THSICE_STEP_TEMP(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     myTime, myIter, myThid )

         CALL THSICE_STEP_FWD(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     prcAtm,
     I                     myTime, myIter, myThid )

         CALL THSICE_AVE(
     I                     bi,bj, myTime, myIter, myThid )

c      ENDDO
c     ENDDO

c       IF ( .FALSE. ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
#endif
        IF ( stressReduction.GT. 0. _d 0 ) THEN
         DO j = jMin, jMax
          DO i = iMin+1,iMax
            tauFac = stressReduction
     &             *(iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
            fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
          ENDDO
         ENDDO
         DO j = jMin+1, jMax
          DO i = iMin, iMax
            tauFac = stressReduction
     &             *(iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
            fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF

C--  end bi,bj loop
       ENDDO
      ENDDO

#ifdef ATMOSPHERIC_LOADING
c     IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)
#endif

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif  /*ALLOW_THSICE*/

      RETURN
      END
1	jmc	1.8	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.7 2006/04/09 17:35:30 heimbach Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	jmc	1.8
6	jmc	1.2	CBOP
7	jmc	1.1	C !ROUTINE: THSICE_MAIN
8			C !INTERFACE:
9	jmc	1.8	SUBROUTINE THSICE_MAIN(
10	jmc	1.1	I myTime, myIter, myThid )
11	jmc	1.2	C !DESCRIPTION: \bv
12	jmc	1.1	C ==========================================================
13	jmc	1.8	C \| S/R THSICE_MAIN
14			C \| o Therm_SeaIce main routine.
15	jmc	1.1	C \| step forward Thermodynamic_SeaIce variables and modify
16			C \| ocean surface forcing accordingly.
17			C ==========================================================
18
19			C !USES:
20			IMPLICIT NONE
21	jmc	1.2
22	jmc	1.1	C === Global variables ===
23			#include "SIZE.h"
24			#include "EEPARAMS.h"
25			#include "PARAMS.h"
26	jmc	1.2	#include "GRID.h"
27	jmc	1.3	#include "SURFACE.h"
28	jmc	1.2	#include "DYNVARS.h"
29	jmc	1.1	#include "FFIELDS.h"
30			#include "THSICE_PARAMS.h"
31	jmc	1.2	#include "THSICE_VARS.h"
32	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
33			# include "tamc.h"
34			# include "tamc_keys.h"
35			#endif
36	jmc	1.8
37	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
38			C === Routine arguments ===
39			C myIter :: iteration counter for this thread
40			C myTime :: time counter for this thread
41			C myThid :: thread number for this instance of the routine.
42			_RL myTime
43			INTEGER myIter
44			INTEGER myThid
45	jmc	1.2	CEOP
46	jmc	1.1
47			#ifdef ALLOW_THSICE
48			C !LOCAL VARIABLES:
49			C === Local variables ===
50			INTEGER i,j
51			INTEGER bi,bj
52			INTEGER iMin, iMax
53			INTEGER jMin, jMax
54	jmc	1.2	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55	jmc	1.8	c _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56			c _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57			c _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58	jmc	1.1
59			_RL tauFac
60
61			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
62
63			IF ( stressReduction.GT. 0. _d 0 ) THEN
64	jmc	1.5	C- needs new Ice Fraction in halo region to apply wind-stress reduction
65	jmc	1.8	iMin = 1-OLx
66			iMax = sNx+OLx-1
67			jMin = 1-OLy
68			jMax = sNy+OLy-1
69	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
70			ELSEIF ( useRealFreshWaterFlux ) THEN
71			C- needs sea-ice loading in part of the halo regions for grad.Phi0surf
72			C to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
73			iMin = 0
74			iMax = sNx+1
75			jMin = 0
76			jMax = sNy+1
77			#endif /* ATMOSPHERIC_LOADING */
78	jmc	1.1	ELSE
79			iMin = 1
80			iMax = sNx
81			jMin = 1
82			jMax = sNy
83			ENDIF
84
85			DO bj=myByLo(myThid),myByHi(myThid)
86			DO bi=myBxLo(myThid),myBxHi(myThid)
87	jmc	1.2
88	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
89			act1 = bi - myBxLo(myThid)
90			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
91			act2 = bj - myByLo(myThid)
92			max2 = myByHi(myThid) - myByLo(myThid) + 1
93			act3 = myThid - 1
94			max3 = nTx*nTy
95			act4 = ikey_dynamics - 1
96			iicekey = (act1 + 1) + act2*max1
97			& + act3max1max2
98			& + act4max1max2*max3
99			#endif /* ALLOW_AUTODIFF_TAMC */
100
101	jmc	1.3	C-- Mixed layer thickness: take the 1rst layer
102			#ifdef NONLIN_FRSURF
103			IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
104			IF ( select_rStar.GT.0 ) THEN
105			DO j = jMin, jMax
106			DO i = iMin, iMax
107			hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
108	jmc	1.8	& *rStarFacC(i,j,bi,bj)
109	jmc	1.3	ENDDO
110			ENDDO
111			ELSE
112			DO j = jMin, jMax
113			DO i = iMin, iMax
114			IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
115			hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
116			ELSE
117			hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
118			ENDIF
119			ENDDO
120			ENDDO
121			ENDIF
122			ELSE
123			#else /* ndef NONLIN_FRSURF */
124			IF (.TRUE.) THEN
125			#endif /* NONLIN_FRSURF */
126			DO j = jMin, jMax
127			DO i = iMin, iMax
128			hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
129			ENDDO
130			ENDDO
131			ENDIF
132
133	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
134			CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
135			CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
136			#endif
137
138	jmc	1.2	DO j = jMin, jMax
139			DO i = iMin, iMax
140			tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
141			sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
142	jmc	1.8	v2ocMxL(i,j,bi,bj) =
143	jmc	1.2	& ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
144			& + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
145			& + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
146			& + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
147			& )*0.5 _d 0
148			prcAtm(i,j) = 0.
149	jmc	1.8	icFrwAtm(i,j,bi,bj) = 0. _d 0
150			icFlxAtm(i,j,bi,bj) = 0. _d 0
151			icFlxSW (i,j,bi,bj) = 0. _d 0
152	jmc	1.2	snowPrc(i,j,bi,bj) = 0. _d 0
153	jmc	1.4	siceAlb(i,j,bi,bj) = 0. _d 0
154	jmc	1.6	ENDDO
155			ENDDO
156	heimbach	1.7
157			#ifdef ALLOW_AUTODIFF_TAMC
158			CADJ STORE iceMask = comlev1, key = iicekey
159			CADJ STORE iceHeight = comlev1, key = iicekey
160			CADJ STORE snowHeight = comlev1, key = iicekey
161			CADJ STORE Tsrf = comlev1, key = iicekey
162			CADJ STORE Qice1 = comlev1, key = iicekey
163			CADJ STORE Qice2 = comlev1, key = iicekey
164			CADJ STORE snowAge = comlev1, key = iicekey
165
166			CADJ STORE sHeating = comlev1, key = iicekey
167			CADJ STORE flxCndBt = comlev1, key = iicekey
168			CADJ STORE snowPrc = comlev1, key = iicekey
169
170			CADJ STORE hOceMxL = comlev1, key = iicekey
171			CADJ STORE tOceMxL = comlev1, key = iicekey
172			CADJ STORE sOceMxL = comlev1, key = iicekey
173			CADJ STORE v2ocMxL = comlev1, key = iicekey
174
175			CADJ STORE empmr = comlev1, key = iicekey
176			CADJ STORE qnet = comlev1, key = iicekey
177			#endif
178
179	jmc	1.2	#ifdef ALLOW_BULK_FORCE
180	jmc	1.6	IF ( useBulkforce ) THEN
181			CALL THSICE_GET_PRECIP(
182			I iceMask,
183	jmc	1.8	O prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
184			O icFlxSW(1-OLx,1-OLy,bi,bj),
185	jmc	1.6	I iMin,iMax,jMin,jMax, bi,bj, myThid )
186			ENDIF
187	jmc	1.2	#endif
188
189	heimbach	1.7
190	jmc	1.8	CALL THSICE_STEP_TEMP(
191			I bi, bj, iMin, iMax, jMin, jMax,
192			I myTime, myIter, myThid )
193
194			CALL THSICE_STEP_FWD(
195			I bi, bj, iMin, iMax, jMin, jMax,
196			I prcAtm,
197	jmc	1.1	I myTime, myIter, myThid )
198	jmc	1.2
199	jmc	1.8	CALL THSICE_AVE(
200			I bi,bj, myTime, myIter, myThid )
201	jmc	1.2
202	jmc	1.1	c ENDDO
203			c ENDDO
204
205			c IF ( .FALSE. ) THEN
206	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
207			CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
208			CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
209			#endif
210	jmc	1.1	IF ( stressReduction.GT. 0. _d 0 ) THEN
211			DO j = jMin, jMax
212			DO i = iMin+1,iMax
213			tauFac = stressReduction
214			& (iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
215			fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
216			ENDDO
217			ENDDO
218			DO j = jMin+1, jMax
219			DO i = iMin, iMax
220			tauFac = stressReduction
221			& (iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
222			fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
223			ENDDO
224			ENDDO
225			ENDIF
226
227			C-- end bi,bj loop
228			ENDDO
229			ENDDO
230
231	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
232			c IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)
233			#endif
234
235	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
236			#endif /ALLOW_THSICE/
237
238			RETURN
239			END