/[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.19 by jmc, Tue Apr 28 18:21:18 2009 UTC revision 1.27 by jmc, Sun Jan 20 15:09:54 2013 UTC
# Line 28  C     === Global variables === Line 28  C     === Global variables ===
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_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
34  # include "tamc.h"  # include "tamc.h"
# Line 51  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  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57  c     _RL flxAtm(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)  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 60  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy Line 61  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy
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    C-    EXF does not provide valid fields in overlap
66           iMin = 1
67           iMax = sNx
68           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  C-     needs new Ice Fraction in halo region to apply wind-stress reduction
72         iMin = 1-OLx         iMin = 1-OLx
73         iMax = sNx+OLx-1         iMax = sNx+OLx-1
74         jMin = 1-OLy         jMin = 1-OLy
75         jMax = sNy+OLy-1         jMax = sNy+OLy-1
76  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
77        ELSEIF ( useRealFreshWaterFlux .AND. .NOT.useSEAICE ) THEN        ELSEIF ( useRealFreshWaterFlux ) THEN
78  C-     needs sea-ice loading in part of the halo regions for grad.Phi0surf  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 )  C      to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
80         iMin = 0         iMin = 0
# Line 93  C      to be valid at the boundaries ( d Line 100  C      to be valid at the boundaries ( d
100            act3 = myThid - 1            act3 = myThid - 1
101            max3 = nTx*nTy            max3 = nTx*nTy
102            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
103            iicekey = (act1 + 1) + act2*max1            ticekey = (act1 + 1) + act2*max1
104       &                         + act3*max1*max2       &                         + act3*max1*max2
105       &                         + act4*max1*max2*max3       &                         + act4*max1*max2*max3
106  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
107    
108  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
109  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
110  CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
111  CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
112  # ifdef ALLOW_EXF  # ifdef ALLOW_EXF
113  CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
114  # endif  # endif
115  #endif  #endif
116    
# Line 111  C--     Mixed layer thickness: take the Line 118  C--     Mixed layer thickness: take the
118  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
119          IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN          IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
120           IF ( select_rStar.GT.0 ) THEN           IF ( select_rStar.GT.0 ) THEN
121            DO j = jMin, jMax            DO j=1-OLy,sNy+OLy
122             DO i = iMin, iMax             DO i=1-OLx,sNx+OLx
123               hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)               hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
124       &                                  *rStarFacC(i,j,bi,bj)       &                                  *rStarFacC(i,j,bi,bj)
125             ENDDO             ENDDO
126            ENDDO            ENDDO
127           ELSE           ELSE
128            DO j = jMin, jMax            DO j=1-OLy,sNy+OLy
129             DO i = iMin, iMax             DO i=1-OLx,sNx+OLx
130              IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN              IF ( kSurfC(i,j,bi,bj).EQ.1 ) THEN
131               hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)               hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
132              ELSE              ELSE
133               hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)               hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
# Line 132  C--     Mixed layer thickness: take the Line 139  C--     Mixed layer thickness: take the
139  #else /* ndef NONLIN_FRSURF */  #else /* ndef NONLIN_FRSURF */
140          IF (.TRUE.) THEN          IF (.TRUE.) THEN
141  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
142            DO j = jMin, jMax            DO j=1-OLy,sNy+OLy
143             DO i = iMin, iMax             DO i=1-OLx,sNx+OLx
144               hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)               hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
145             ENDDO             ENDDO
146            ENDDO            ENDDO
147          ENDIF          ENDIF
148    
149  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
150  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
151  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
152  #endif  #endif
153    
154          DO j = jMin, jMax          DO j = jMin, jMax
# Line 154  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_ Line 161  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_
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            icFrwAtm(i,j,bi,bj) = 0. _d 0            icFrwAtm(i,j,bi,bj) = 0. _d 0
166            icFlxAtm(i,j,bi,bj) = 0. _d 0            icFlxAtm(i,j,bi,bj) = 0. _d 0
167            icFlxSW (i,j,bi,bj) = 0. _d 0            icFlxSW (i,j,bi,bj) = 0. _d 0
168            snowPrc(i,j,bi,bj) = 0. _d 0            snowPrc (i,j,bi,bj) = 0. _d 0
169            siceAlb(i,j,bi,bj) = 0. _d 0            siceAlb (i,j,bi,bj) = 0. _d 0
170           ENDDO           ENDDO
171          ENDDO          ENDDO
172    
173  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
174  CADJ STORE iceMask = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
175  CADJ STORE iceHeight  = comlev1, key = iicekey  CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
176  CADJ STORE snowHeight = comlev1, key = iicekey  CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
177  CADJ STORE Tsrf    = comlev1, key = iicekey  CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
178  CADJ STORE Qice1   = comlev1, key = iicekey  CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
179  CADJ STORE Qice2   = comlev1, key = iicekey  CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
180  CADJ STORE snowAge = comlev1, key = iicekey  CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
181  CADJ STORE snowPrc  = comlev1, key = iicekey  CADJ STORE snowPrc(:,:,bi,bj)  = comlev1_bibj, key = ticekey
182    
183  CADJ STORE hOceMxL = comlev1, key = iicekey  CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
184  CADJ STORE tOceMxL = comlev1, key = iicekey  CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
185  CADJ STORE sOceMxL = comlev1, key = iicekey  CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
186  CADJ STORE v2ocMxL = comlev1, key = iicekey  CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
   
 CADJ STORE empmr   = comlev1, key = iicekey  
 CADJ STORE qnet    = comlev1, key = iicekey  
187  #endif  #endif
188    
189  C-      do sea-ice advection before getting surface fluxes  C-      do sea-ice advection before getting surface fluxes
# Line 192  C Note: will inline this S/R once thSIce Line 196  C Note: will inline this S/R once thSIce
196          IF ( useBulkforce ) THEN          IF ( useBulkforce ) THEN
197           CALL THSICE_GET_PRECIP(           CALL THSICE_GET_PRECIP(
198       I                  iceMask,       I                  iceMask,
199       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
200         O                  snowPrc(1-OLx,1-OLy,bi,bj),
201       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
202       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
203          ENDIF          ENDIF
# Line 201  C Note: will inline this S/R once thSIce Line 206  C Note: will inline this S/R once thSIce
206          IF ( useEXF ) THEN          IF ( useEXF ) THEN
207           CALL THSICE_MAP_EXF(           CALL THSICE_MAP_EXF(
208       I                  iceMask,       I                  iceMask,
209       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       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),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
212       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
213          ENDIF          ENDIF
214  #endif  #endif
215    
216          CALL THSICE_STEP_TEMP(  #ifdef ALLOW_AUTODIFF_TAMC
217    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
218    CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
219    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,       I                     bi, bj, iMin, iMax, jMin, jMax,
225       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
226    
227  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
228  CADJ STORE empmr, qnet = comlev1, key = iicekey  CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
229  CADJ STORE iceMask = comlev1, key = iicekey  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
230  CADJ STORE iceHeight  = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
231  CADJ STORE snowHeight = comlev1, key = iicekey  CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
232  CADJ STORE Tsrf    = comlev1, key = iicekey  CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
233  CADJ STORE Qice1   = comlev1, key = iicekey  cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
234  CADJ STORE Qice2   = comlev1, key = iicekey  CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
235  CADJ STORE snowAge = comlev1, key = iicekey  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  #endif
239              CALL THSICE_STEP_FWD(
         CALL THSICE_STEP_FWD(  
240       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
241       I                     prcAtm,       I                     prcAtm(1-OLx,1-OLy,bi,bj),
242       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
243    #ifndef ALLOW_AUTODIFF_TAMC
244            ENDIF
245    #endif
246    
247    C--  end bi,bj loop
248           ENDDO
249          ENDDO
250    
251    #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(          CALL THSICE_AVE(
263       I                     bi,bj, myTime, myIter, myThid )       I                     bi,bj, myTime, myIter, myThid )
264           ENDDO
265          ENDDO
266    
267  c      ENDDO  C     add a small piece of code to check AddFluid implementation:
268  c     ENDDO  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.
282         &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
283         &  _EXCH_XY_RS( sIceLoad, myThid )
284    #endif
285    
286          DO bj=myByLo(myThid),myByHi(myThid)
287           DO bi=myBxLo(myThid),myBxHi(myThid)
288  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
289  C--   and stressReduction is always set to zero  C--   and stressReduction is always set to zero
290  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
291  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
292  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
293  #endif  #endif
294          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
295            DO j = jMin, jMax            DO j = jMin, jMax
# Line 260  C--  end bi,bj loop Line 312  C--  end bi,bj loop
312         ENDDO         ENDDO
313        ENDDO        ENDDO
314    
 C     add a small piece of code to check AddFluid implementation:  
 c#include "thsice_test_addfluid.h"  
   
       IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN  
 C--   Exchange fields that are advected by seaice dynamics  
         _EXCH_XY_RL( iceMask, myThid )  
         _EXCH_XY_RL( iceHeight, myThid )  
         _EXCH_XY_RL( snowHeight, myThid )  
         _EXCH_XY_RL( Qice1, myThid )  
         _EXCH_XY_RL( Qice2, myThid )  
 #ifdef ATMOSPHERIC_LOADING  
         IF (useRealFreshWaterFlux)  
      &  _EXCH_XY_RS( sIceLoad, myThid )  
 #endif  
       ENDIF  
   
315  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
317    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22