/[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.20 by jmc, Thu Jun 11 19:23:54 2009 UTC revision 1.31 by jmc, Thu Apr 4 16:44:34 2013 UTC
# Line 23  C     === Global variables === Line 23  C     === Global variables ===
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
 #include "SURFACE.h"  
 #include "DYNVARS.h"  
26  #include "FFIELDS.h"  #include "FFIELDS.h"
27  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
28    #include "THSICE_SIZE.h"
29  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
30  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
31    # include "DYNVARS.h"
32  # include "tamc.h"  # include "tamc.h"
33  # include "tamc_keys.h"  # include "tamc_keys.h"
34  #endif  #endif
# Line 47  CEOP Line 46  CEOP
46  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
47  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
48  C     === Local variables ===  C     === Local variables ===
49    C     qPrcRn    :: Energy content of Precip+RunOff (+=down) [W/m2]
50        INTEGER i,j        INTEGER i,j
51        INTEGER bi,bj        INTEGER bi,bj
52        INTEGER iMin, iMax        INTEGER iMin, iMax
53        INTEGER jMin, jMax        INTEGER jMin, jMax
54        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55          _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
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
   
 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=1-OLy,sNy+OLy  
            DO i=1-OLx,sNx+OLx  
              hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)  
      &                                  *rStarFacC(i,j,bi,bj)  
            ENDDO  
           ENDDO  
          ELSE  
           DO j=1-OLy,sNy+OLy  
            DO i=1-OLx,sNx+OLx  
             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=1-OLy,sNy+OLy  
            DO i=1-OLx,sNx+OLx  
              hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)  
            ENDDO  
           ENDDO  
         ENDIF  
   
116  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
117  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
118  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
119  #endif  #endif
120    
121          DO j = jMin, jMax          DO j=1-OLy,sNy+OLy
122           DO i = iMin, iMax           DO i=1-OLx,sNx+OLx
123            tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)            prcAtm  (i,j,bi,bj) = 0. _d 0
124            sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)            qPrcRn  (i,j) = 0. _d 0
           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  
125           ENDDO           ENDDO
126          ENDDO          ENDDO
127    
128  #ifdef ALLOW_AUTODIFF_TAMC          CALL THSICE_GET_OCEAN(
129  CADJ STORE iceMask = comlev1, key = iicekey       I                         bi, bj, myTime, myIter, myThid )
 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 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  
130    
131  CADJ STORE empmr   = comlev1, key = iicekey  #ifdef ALLOW_AUTODIFF_TAMC
132  CADJ STORE qnet    = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
133    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
134    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
135    CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
136    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
137    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
138    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
139    CADJ STORE snowPrc(:,:,bi,bj)  = comlev1_bibj, key = ticekey
140    
141    CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
142    CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
143    CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
144    CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
145  #endif  #endif
146    
147    #ifdef OLD_THSICE_CALL_SEQUENCE
148  C-      do sea-ice advection before getting surface fluxes  C-      do sea-ice advection before getting surface fluxes
149  C Note: will inline this S/R once thSIce in Atmos. set-up is settled  C Note: will inline this S/R once thSIce in Atmos. set-up is settled
150          IF ( thSIceAdvScheme.GT.0 )          IF ( thSIceAdvScheme.GT.0 )
151       &   CALL THSICE_DO_ADVECT(       &   CALL THSICE_DO_ADVECT(
152       I                   bi,bj, myTime, myIter, myThid )       I                   bi,bj, myTime, myIter, myThid )
153    #endif /* OLD_THSICE_CALL_SEQUENCE */
154    
155  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
156          IF ( useBulkforce ) THEN          IF ( useBulkforce ) THEN
157           CALL THSICE_GET_PRECIP(           CALL THSICE_GET_PRECIP(
158       I                  iceMask,       I                  iceMask,
159       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
160         O                  snowPrc(1-OLx,1-OLy,bi,bj),
161       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
162       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
163          ENDIF          ENDIF
# Line 200  C Note: will inline this S/R once thSIce Line 165  C Note: will inline this S/R once thSIce
165  #ifdef ALLOW_EXF  #ifdef ALLOW_EXF
166          IF ( useEXF ) THEN          IF ( useEXF ) THEN
167           CALL THSICE_MAP_EXF(           CALL THSICE_MAP_EXF(
168       I                  iceMask,       I                  iceMask, tOceMxL,
169       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
170         O                  snowPrc(1-OLx,1-OLy,bi,bj), qPrcRn,
171       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
172       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
173          ENDIF          ENDIF
174  #endif  #endif
175    
176          CALL THSICE_STEP_TEMP(  #ifdef ALLOW_AUTODIFF_TAMC
177    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
178    CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
179    CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
180    #else
181            IF ( .NOT.thSIce_skipThermo ) THEN
182    #endif
183              CALL THSICE_STEP_TEMP(
184       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
185       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
186    
187  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
188  CADJ STORE empmr, qnet = comlev1, key = iicekey  CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
189  CADJ STORE iceMask = comlev1, key = iicekey  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
190  CADJ STORE iceHeight  = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
191  CADJ STORE snowHeight = comlev1, key = iicekey  CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
192  CADJ STORE Tsrf    = comlev1, key = iicekey  CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
193  CADJ STORE Qice1   = comlev1, key = iicekey  cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
194  CADJ STORE Qice2   = comlev1, key = iicekey  CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
195  CADJ STORE snowAge = comlev1, key = iicekey  CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
196    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
197    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
198    #else
199            ENDIF
200            IF ( .NOT.thSIce_skipThermo ) THEN
201  #endif  #endif
202              CALL THSICE_STEP_FWD(
         CALL THSICE_STEP_FWD(  
203       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
204       I                     prcAtm,       I                     prcAtm(1-OLx,1-OLy,bi,bj), qPrcRn,
205       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
206    #ifndef ALLOW_AUTODIFF_TAMC
207            ENDIF
208    #endif
209    
210          CALL THSICE_AVE(  C--  end bi,bj loop
211       I                     bi,bj, myTime, myIter, myThid )         ENDDO
212          ENDDO
213    
214  c      ENDDO  #ifdef ALLOW_BALANCE_FLUXES
215  c     ENDDO  C--   Balance net Fresh-Water flux from Atm+Land
216          IF ( thSIceBalanceAtmFW.NE.0 ) THEN
217            CALL THSICE_BALANCE_FRW(
218         I                      iMin, iMax, jMin, jMax,
219         I                      prcAtm, myTime, myIter, myThid )
220          ENDIF
221    #endif
222    
223    C     add a small piece of code to check AddFluid implementation:
224    c#include "thsice_test_addfluid.h"
225    
226    C--   Exchange fields that are advected by seaice dynamics
227          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
228         &       .OR. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
229            CALL THSICE_DO_EXCH( myThid )
230          ENDIF
231    #ifdef OLD_THSICE_CALL_SEQUENCE
232    #ifdef ATMOSPHERIC_LOADING
233          IF ( useRealFreshWaterFlux .AND.
234         &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
235         &  _EXCH_XY_RS( sIceLoad, myThid )
236    #endif
237    #else /* OLD_THSICE_CALL_SEQUENCE */
238    #ifdef ATMOSPHERIC_LOADING
239          IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
240         &                           .AND. thSIceAdvScheme.LE.0 )
241         &  _EXCH_XY_RS( sIceLoad, myThid )
242    #endif
243    
244    C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
245    C     otherwise, call it from here, after thsice-thermodynamics is done
246          IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
247             CALL THSICE_DO_ADVECT(
248         I                          0, 0, myTime, myIter, myThid )
249          ENDIF
250    #endif /* OLD_THSICE_CALL_SEQUENCE */
251    
252          DO bj=myByLo(myThid),myByHi(myThid)
253           DO bi=myBxLo(myThid),myBxHi(myThid)
254    C--   Cumulate time-averaged fields and also fill-up flux diagnostics
255    C     (if not done in THSICE_DO_ADVECT call)
256    #ifdef OLD_THSICE_CALL_SEQUENCE
257            IF ( .TRUE. ) THEN
258    #else /* OLD_THSICE_CALL_SEQUENCE */
259            IF ( thSIceAdvScheme.LE.0 ) THEN
260    #endif /* OLD_THSICE_CALL_SEQUENCE */
261             CALL THSICE_AVE(
262         I                     bi,bj, myTime, myIter, myThid )
263            ENDIF
264  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
265  C--   and stressReduction is always set to zero  C--   and stressReduction is always set to zero
266  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
267  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
268  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
269  #endif  #endif
270          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
271            DO j = jMin, jMax            DO j = jMin, jMax
# Line 260  C--  end bi,bj loop Line 288  C--  end bi,bj loop
288         ENDDO         ENDDO
289        ENDDO        ENDDO
290    
 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  
   
291  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
292  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
293    
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22