pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.4 2005/01/31 19:37:06 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"
#ifdef ALLOW_BULK_FORCE
#include "BULKF.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)
      _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
#ifdef ALLOW_BULK_FORCE
           prcAtm(i,j) = ( rain(i,j,bi,bj)+runoff(i,j,bi,bj) )*rhofw
           flxSW (i,j) = solar(i,j,bi,bj)
           IF ( iceMask(i,j,bi,bj).GT.0. _d 0
     &       .AND. Tair(i,j,bi,bj).LE.Tf0kel )  THEN
             snowPrc(i,j,bi,bj) = rain(i,j,bi,bj)*rhofw
           ENDIF
#endif
          ENDDO
         ENDDO

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