pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.5 2005/02/11 19:33:59 jmc 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"
 
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)
      _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _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)

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

         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.
           evpAtm(i,j) = 0.
           flxSW (i,j) = 0.
           snowPrc(i,j,bi,bj) = 0. _d 0
           siceAlb(i,j,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
#ifdef ALLOW_BULK_FORCE
         IF ( useBulkforce ) THEN
           CALL THSICE_GET_PRECIP(
     I                  iceMask,
     O                  prcAtm, snowPrc(1-Olx,1-Oly,bi,bj), flxSW,
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
         ENDIF
#endif

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

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

c      ENDDO
c     ENDDO

c       IF ( .FALSE. ) THEN
        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.6	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.5 2005/02/11 19:33:59 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5
6	jmc	1.2	CBOP
7	jmc	1.1	C !ROUTINE: THSICE_MAIN
8			C !INTERFACE:
9			SUBROUTINE THSICE_MAIN(
10			I myTime, myIter, myThid )
11	jmc	1.2	C !DESCRIPTION: \bv
12	jmc	1.1	C ==========================================================
13	jmc	1.2	C \| S/R THSICE_MAIN
14	jmc	1.1	C \| o Therm_SeaIce main routine.
15			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	jmc	1.1
33			C !INPUT/OUTPUT PARAMETERS:
34			C === Routine arguments ===
35			C myIter :: iteration counter for this thread
36			C myTime :: time counter for this thread
37			C myThid :: thread number for this instance of the routine.
38			_RL myTime
39			INTEGER myIter
40			INTEGER myThid
41	jmc	1.2	CEOP
42	jmc	1.1
43			#ifdef ALLOW_THSICE
44			C !LOCAL VARIABLES:
45			C === Local variables ===
46			INTEGER i,j
47			INTEGER bi,bj
48			INTEGER iMin, iMax
49			INTEGER jMin, jMax
50	jmc	1.2	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51			_RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52			_RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53	jmc	1.1
54			_RL tauFac
55
56			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
57
58			IF ( stressReduction.GT. 0. _d 0 ) THEN
59	jmc	1.5	C- needs new Ice Fraction in halo region to apply wind-stress reduction
60	jmc	1.1	iMin = 1-Olx
61			iMax = sNx+Olx-1
62			jMin = 1-Oly
63			jMax = sNy+Oly-1
64	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
65			ELSEIF ( useRealFreshWaterFlux ) THEN
66			C- needs sea-ice loading in part of the halo regions for grad.Phi0surf
67			C to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
68			iMin = 0
69			iMax = sNx+1
70			jMin = 0
71			jMax = sNy+1
72			#endif /* ATMOSPHERIC_LOADING */
73	jmc	1.1	ELSE
74			iMin = 1
75			iMax = sNx
76			jMin = 1
77			jMax = sNy
78			ENDIF
79
80			DO bj=myByLo(myThid),myByHi(myThid)
81			DO bi=myBxLo(myThid),myBxHi(myThid)
82	jmc	1.2
83	jmc	1.3	C-- Mixed layer thickness: take the 1rst layer
84			#ifdef NONLIN_FRSURF
85			IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
86			IF ( select_rStar.GT.0 ) THEN
87			DO j = jMin, jMax
88			DO i = iMin, iMax
89			hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
90			& *rStarFacC(i,j,bi,bj)
91			ENDDO
92			ENDDO
93			ELSE
94			DO j = jMin, jMax
95			DO i = iMin, iMax
96			IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
97			hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
98			ELSE
99			hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
100			ENDIF
101			ENDDO
102			ENDDO
103			ENDIF
104			ELSE
105			#else /* ndef NONLIN_FRSURF */
106			IF (.TRUE.) THEN
107			#endif /* NONLIN_FRSURF */
108			DO j = jMin, jMax
109			DO i = iMin, iMax
110			hOceMxL(i,j,bi,bj) = drF(1)*hfacC(i,j,1,bi,bj)
111			ENDDO
112			ENDDO
113			ENDIF
114
115	jmc	1.2	DO j = jMin, jMax
116			DO i = iMin, iMax
117			tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
118			sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
119			v2ocMxL(i,j,bi,bj) =
120			& ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
121			& + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
122			& + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
123			& + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
124			& )*0.5 _d 0
125			prcAtm(i,j) = 0.
126			evpAtm(i,j) = 0.
127			flxSW (i,j) = 0.
128			snowPrc(i,j,bi,bj) = 0. _d 0
129	jmc	1.4	siceAlb(i,j,bi,bj) = 0. _d 0
130	jmc	1.6	ENDDO
131			ENDDO
132	jmc	1.2	#ifdef ALLOW_BULK_FORCE
133	jmc	1.6	IF ( useBulkforce ) THEN
134			CALL THSICE_GET_PRECIP(
135			I iceMask,
136			O prcAtm, snowPrc(1-Olx,1-Oly,bi,bj), flxSW,
137			I iMin,iMax,jMin,jMax, bi,bj, myThid )
138			ENDIF
139	jmc	1.2	#endif
140
141	jmc	1.1	CALL THSICE_STEP_FWD(
142			I bi, bj, iMin, iMax, jMin, jMax,
143	jmc	1.2	I prcAtm,
144			U evpAtm, flxSW,
145	jmc	1.1	I myTime, myIter, myThid )
146	jmc	1.2
147			CALL THSICE_AVE(
148			I evpAtm, flxSW,
149			I bi,bj, myTime, myIter, myThid )
150
151	jmc	1.1	c ENDDO
152			c ENDDO
153
154			c IF ( .FALSE. ) THEN
155			IF ( stressReduction.GT. 0. _d 0 ) THEN
156			DO j = jMin, jMax
157			DO i = iMin+1,iMax
158			tauFac = stressReduction
159			& (iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
160			fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
161			ENDDO
162			ENDDO
163			DO j = jMin+1, jMax
164			DO i = iMin, iMax
165			tauFac = stressReduction
166			& (iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
167			fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
168			ENDDO
169			ENDDO
170			ENDIF
171
172			C-- end bi,bj loop
173			ENDDO
174			ENDDO
175
176	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
177			c IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)
178			#endif
179
180	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
181			#endif /ALLOW_THSICE/
182
183			RETURN
184			END