/[MITgcm]/MITgcm/pkg/thsice/thsice_main.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_main.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by jmc, Wed Apr 7 23:40:34 2004 UTC revision 1.26 by jmc, Wed Nov 21 01:53:34 2012 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "THSICE_OPTIONS.h"  #include "THSICE_OPTIONS.h"
5    
6  CBOP  CBOP
7  C     !ROUTINE: THSICE_MAIN  C     !ROUTINE: THSICE_MAIN
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE THSICE_MAIN(        SUBROUTINE THSICE_MAIN(
10       I                        myTime, myIter, myThid )       I                        myTime, myIter, myThid )
11  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | S/R  THSICE_MAIN              C     | S/R  THSICE_MAIN
14  C     | o Therm_SeaIce main routine.  C     | o Therm_SeaIce main routine.
15  C     |   step forward Thermodynamic_SeaIce variables and modify  C     |   step forward Thermodynamic_SeaIce variables and modify
16  C     |    ocean surface forcing accordingly.  C     |    ocean surface forcing accordingly.
17  C     *==========================================================*  C     *==========================================================*
# Line 24  C     === Global variables === Line 24  C     === Global variables ===
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "GRID.h"  #include "GRID.h"
27    #include "SURFACE.h"
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "FFIELDS.h"  #include "FFIELDS.h"
30  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
31    #include "THSICE_SIZE.h"
32  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
33  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_AUTODIFF_TAMC
34  #include "BULKF.h"  # include "tamc.h"
35    # include "tamc_keys.h"
36  #endif  #endif
37    
38  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
39  C     === Routine arguments ===  C     === Routine arguments ===
40  C     myIter :: iteration counter for this thread  C     myTime    :: Current time in simulation (s)
41  C     myTime :: time counter for this thread  C     myIter    :: Current iteration number
42  C     myThid :: thread number for this instance of the routine.  C     myThid    :: My Thread Id. number
43        _RL  myTime        _RL     myTime
44        INTEGER myIter        INTEGER myIter
45        INTEGER myThid        INTEGER myThid
46  CEOP  CEOP
# Line 49  C     === Local variables === Line 52  C     === Local variables ===
52        INTEGER bi,bj        INTEGER bi,bj
53        INTEGER iMin, iMax        INTEGER iMin, iMax
54        INTEGER jMin, jMax        INTEGER jMin, jMax
55        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56        _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57        _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58    c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59    
60        _RL tauFac        _RL tauFac
61    
62  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
63    
64        IF ( stressReduction.GT. 0. _d 0 ) THEN        IF ( useEXF .OR. useSEAICE ) THEN
65         iMin = 1-Olx  C-    EXF does not provide valid fields in overlap
66         iMax = sNx+Olx-1         iMin = 1
67         jMin = 1-Oly         iMax = sNx
68         jMax = sNy+Oly-1         jMin = 1
69           jMax = sNy
70          ELSEIF ( stressReduction.GT. 0. _d 0 ) THEN
71    C-     needs new Ice Fraction in halo region to apply wind-stress reduction
72           iMin = 1-OLx
73           iMax = sNx+OLx-1
74           jMin = 1-OLy
75           jMax = sNy+OLy-1
76    #ifdef ATMOSPHERIC_LOADING
77          ELSEIF ( useRealFreshWaterFlux ) THEN
78    C-     needs sea-ice loading in part of the halo regions for grad.Phi0surf
79    C      to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
80           iMin = 0
81           iMax = sNx+1
82           jMin = 0
83           jMax = sNy+1
84    #endif /* ATMOSPHERIC_LOADING */
85        ELSE        ELSE
86         iMin = 1         iMin = 1
87         iMax = sNx         iMax = sNx
# Line 72  C---+----1----+----2----+----3----+----4 Line 92  C---+----1----+----2----+----3----+----4
92        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
93         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
94    
95           DO j = jMin, jMax  #ifdef ALLOW_AUTODIFF_TAMC
96            DO i = iMin, iMax            act1 = bi - myBxLo(myThid)
97             hOceMxL(i,j,bi,bj) = hfacC(i,j,1,bi,bj)*drF(1)            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
98             tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)            act2 = bj - myByLo(myThid)
99             sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)            max2 = myByHi(myThid) - myByLo(myThid) + 1
100             v2ocMxL(i,j,bi,bj) =            act3 = myThid - 1
101              max3 = nTx*nTy
102              act4 = ikey_dynamics - 1
103              ticekey = (act1 + 1) + act2*max1
104         &                         + act3*max1*max2
105         &                         + act4*max1*max2*max3
106    #endif /* ALLOW_AUTODIFF_TAMC */
107    
108    #ifdef ALLOW_AUTODIFF_TAMC
109    CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
110    CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
111    CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
112    # ifdef ALLOW_EXF
113    CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
114    # endif
115    #endif
116    
117    C--     Mixed layer thickness: take the 1rst layer
118    #ifdef NONLIN_FRSURF
119            IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
120             IF ( select_rStar.GT.0 ) THEN
121              DO j=1-OLy,sNy+OLy
122               DO i=1-OLx,sNx+OLx
123                 hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
124         &                                  *rStarFacC(i,j,bi,bj)
125               ENDDO
126              ENDDO
127             ELSE
128              DO j=1-OLy,sNy+OLy
129               DO i=1-OLx,sNx+OLx
130                IF ( kSurfC(i,j,bi,bj).EQ.1 ) THEN
131                 hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
132                ELSE
133                 hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
134                ENDIF
135               ENDDO
136              ENDDO
137             ENDIF
138            ELSE
139    #else /* ndef NONLIN_FRSURF */
140            IF (.TRUE.) THEN
141    #endif /* NONLIN_FRSURF */
142              DO j=1-OLy,sNy+OLy
143               DO i=1-OLx,sNx+OLx
144                 hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
145               ENDDO
146              ENDDO
147            ENDIF
148    
149    #ifdef ALLOW_AUTODIFF_TAMC
150    CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
151    CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
152    #endif
153    
154            DO j = jMin, jMax
155             DO i = iMin, iMax
156              tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
157              sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
158              v2ocMxL(i,j,bi,bj) =
159       &              ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)       &              ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
160       &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)       &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
161       &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)       &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
162       &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)       &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
163       &              )*0.5 _d 0       &              )*0.5 _d 0
164             prcAtm(i,j) = 0.            prcAtm  (i,j,bi,bj) = 0. _d 0
165             evpAtm(i,j) = 0.            icFrwAtm(i,j,bi,bj) = 0. _d 0
166             flxSW (i,j) = 0.            icFlxAtm(i,j,bi,bj) = 0. _d 0
167             snowPrc(i,j,bi,bj) = 0. _d 0            icFlxSW (i,j,bi,bj) = 0. _d 0
168              snowPrc (i,j,bi,bj) = 0. _d 0
169              siceAlb (i,j,bi,bj) = 0. _d 0
170             ENDDO
171            ENDDO
172    
173    #ifdef ALLOW_AUTODIFF_TAMC
174    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
175    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
176    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
177    CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
178    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
179    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
180    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
181    CADJ STORE snowPrc(:,:,bi,bj)  = comlev1_bibj, key = ticekey
182    
183    CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
184    CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
185    CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
186    CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
187    #endif
188    
189    C-      do sea-ice advection before getting surface fluxes
190    C Note: will inline this S/R once thSIce in Atmos. set-up is settled
191            IF ( thSIceAdvScheme.GT.0 )
192         &   CALL THSICE_DO_ADVECT(
193         I                   bi,bj, myTime, myIter, myThid )
194    
195  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
196             prcAtm(i,j) = ( rain(i,j,bi,bj)+runoff(i,j,bi,bj) )*rhofw          IF ( useBulkforce ) THEN
197             flxSW (i,j) = solar(i,j,bi,bj)           CALL THSICE_GET_PRECIP(
198             IF ( iceMask(i,j,bi,bj).GT.0. _d 0       I                  iceMask,
199       &       .AND. Tair(i,j,bi,bj).LE.Tf0kel )  THEN       O                  prcAtm(1-OLx,1-OLy,bi,bj),
200               snowPrc(i,j,bi,bj) = rain(i,j,bi,bj)*rhofw       O                  snowPrc(1-OLx,1-OLy,bi,bj),
201             ENDIF       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
202         I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
203            ENDIF
204    #endif
205    #ifdef ALLOW_EXF
206            IF ( useEXF ) THEN
207             CALL THSICE_MAP_EXF(
208         I                  iceMask,
209         O                  prcAtm(1-OLx,1-OLy,bi,bj),
210         O                  snowPrc(1-OLx,1-OLy,bi,bj),
211         O                  icFlxSW(1-OLx,1-OLy,bi,bj),
212         I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
213            ENDIF
214  #endif  #endif
           ENDDO  
          ENDDO  
215    
216           CALL THSICE_STEP_FWD(  #ifdef ALLOW_AUTODIFF_TAMC
217       I                     bi, bj, iMin, iMax, jMin, jMax,  CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
218       I                     prcAtm,  CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
219       U                     evpAtm, flxSW,  CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
220    #else
221            IF ( .NOT.thSIce_skipThermo ) THEN
222    #endif
223              CALL THSICE_STEP_TEMP(
224         I                     bi, bj, iMin, iMax, jMin, jMax,
225       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
226    
227           CALL THSICE_AVE(  #ifdef ALLOW_AUTODIFF_TAMC
228       I                   evpAtm, flxSW,  CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
229       I                   bi,bj, myTime, myIter, myThid )  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
230    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
231    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
232    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
233    cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
234    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
235    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
236    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
237    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
238    #endif
239              CALL THSICE_STEP_FWD(
240         I                     bi, bj, iMin, iMax, jMin, jMax,
241         I                     prcAtm(1-OLx,1-OLy,bi,bj),
242         I                     myTime, myIter, myThid )
243    #ifndef ALLOW_AUTODIFF_TAMC
244            ENDIF
245    #endif
246    
247  c      ENDDO  C--  end bi,bj loop
248  c     ENDDO         ENDDO
249          ENDDO
250    
251  c       IF ( .FALSE. ) THEN  #ifdef ALLOW_BALANCE_FLUXES
252    C--   Balance net Fresh-Water flux from Atm+Land
253          IF ( thSIceBalanceAtmFW.NE.0 ) THEN
254            CALL THSICE_BALANCE_FRW(
255         I                      iMin, iMax, jMin, jMax,
256         I                      prcAtm, myTime, myIter, myThid )
257          ENDIF
258    #endif
259    
260          DO bj=myByLo(myThid),myByHi(myThid)
261           DO bi=myBxLo(myThid),myBxHi(myThid)
262            CALL THSICE_AVE(
263         I                     bi,bj, myTime, myIter, myThid )
264           ENDDO
265          ENDDO
266    
267    C     add a small piece of code to check AddFluid implementation:
268    c#include "thsice_test_addfluid.h"
269    
270          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN
271    C--   Exchange fields that are advected by seaice dynamics
272            _EXCH_XY_RL( iceMask, myThid )
273            _EXCH_XY_RL( iceHeight, myThid )
274            _EXCH_XY_RL( snowHeight, myThid )
275            _EXCH_XY_RL( Qice1, myThid )
276            _EXCH_XY_RL( Qice2, myThid )
277          ELSEIF ( useEXF .AND. stressReduction.GT. 0. _d 0 ) THEN
278            _EXCH_XY_RL( iceMask, myThid )
279          ENDIF
280    #ifdef ATMOSPHERIC_LOADING
281          IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE ) )
282         &  _EXCH_XY_RS( sIceLoad, myThid )
283    #endif
284    
285          DO bj=myByLo(myThid),myByHi(myThid)
286           DO bi=myBxLo(myThid),myBxHi(myThid)
287    C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
288    C--   and stressReduction is always set to zero
289    #ifdef ALLOW_AUTODIFF_TAMC
290    CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
291    CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
292    #endif
293          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
294           DO j = jMin, jMax            DO j = jMin, jMax
295            DO i = iMin+1,iMax             DO i = iMin+1,iMax
296              tauFac = stressReduction              tauFac = stressReduction
297       &             *(iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0       &             *(iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
298              fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)              fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
299               ENDDO
300            ENDDO            ENDDO
301           ENDDO            DO j = jMin+1, jMax
302           DO j = jMin+1, jMax             DO i = iMin, iMax
           DO i = iMin, iMax  
303              tauFac = stressReduction              tauFac = stressReduction
304       &             *(iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0       &             *(iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
305              fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)              fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
306               ENDDO
307            ENDDO            ENDDO
          ENDDO  
308          ENDIF          ENDIF
309    
310  C--  end bi,bj loop  C--  end bi,bj loop
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22