/[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.16 by heimbach, Thu Apr 19 13:22:36 2007 UTC revision 1.29 by jmc, Sat Mar 30 18:38:50 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"
 C--  
 # include "THSICE_2DYN.h"  
 C--  
34  #endif  #endif
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 54  C     === Local variables === Line 50  C     === Local variables ===
50        INTEGER bi,bj        INTEGER bi,bj
51        INTEGER iMin, iMax        INTEGER iMin, iMax
52        INTEGER jMin, jMax        INTEGER jMin, jMax
53        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55  c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 63  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy Line 59  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy
59    
60  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
61    
62        IF ( stressReduction.GT. 0. _d 0 ) THEN        IF ( useEXF .OR. useSEAICE ) THEN
63    C-    EXF does not provide valid fields in overlap
64           iMin = 1
65           iMax = sNx
66           jMin = 1
67           jMax = sNy
68          ELSEIF ( stressReduction.GT. 0. _d 0 ) THEN
69  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
70         iMin = 1-OLx         iMin = 1-OLx
71         iMax = sNx+OLx-1         iMax = sNx+OLx-1
72         jMin = 1-OLy         jMin = 1-OLy
73         jMax = sNy+OLy-1         jMax = sNy+OLy-1
74  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
75        ELSEIF ( useRealFreshWaterFlux .AND. .NOT.useSEAICE ) THEN        ELSEIF ( useRealFreshWaterFlux ) THEN
76  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
77  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 )
78         iMin = 0         iMin = 0
# Line 96  C      to be valid at the boundaries ( d Line 98  C      to be valid at the boundaries ( d
98            act3 = myThid - 1            act3 = myThid - 1
99            max3 = nTx*nTy            max3 = nTx*nTy
100            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
101            iicekey = (act1 + 1) + act2*max1            ticekey = (act1 + 1) + act2*max1
102       &                         + act3*max1*max2       &                         + act3*max1*max2
103       &                         + act4*max1*max2*max3       &                         + act4*max1*max2*max3
104  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
105    
106  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
107  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
108  CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
109  CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
110  # ifdef ALLOW_EXF  # ifdef ALLOW_EXF
111  CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
112  # endif  # endif
113  #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 = 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  
   
114  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
115  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
116  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
117  #endif  #endif
118    
119          DO j = jMin, jMax          DO j=1-OLy,sNy+OLy
120           DO i = iMin, iMax           DO i=1-OLx,sNx+OLx
121            tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)            prcAtm  (i,j,bi,bj) = 0. _d 0
           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.  
           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  
122           ENDDO           ENDDO
123          ENDDO          ENDDO
124    
125  #ifdef ALLOW_AUTODIFF_TAMC          CALL THSICE_GET_OCEAN(
126  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  
127    
128  CADJ STORE empmr   = comlev1, key = iicekey  #ifdef ALLOW_AUTODIFF_TAMC
129  CADJ STORE qnet    = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
130    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
131    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
132    CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
133    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
134    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
135    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
136    CADJ STORE snowPrc(:,:,bi,bj)  = comlev1_bibj, key = ticekey
137    
138    CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
139    CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
140    CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
141    CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
142  #endif  #endif
143    
144    #ifdef OLD_THSICE_CALL_SEQUENCE
145  C-      do sea-ice advection before getting surface fluxes  C-      do sea-ice advection before getting surface fluxes
146  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
147          IF ( thSIceAdvScheme.GT.0 )          IF ( thSIceAdvScheme.GT.0 )
148       &   CALL THSICE_DO_ADVECT(       &   CALL THSICE_DO_ADVECT(
149       I                   bi,bj, myTime, myIter, myThid )       I                   bi,bj, myTime, myIter, myThid )
150    #endif /* OLD_THSICE_CALL_SEQUENCE */
151    
152  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
153          IF ( useBulkforce ) THEN          IF ( useBulkforce ) THEN
154           CALL THSICE_GET_PRECIP(           CALL THSICE_GET_PRECIP(
155       I                  iceMask,       I                  iceMask,
156       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
157         O                  snowPrc(1-OLx,1-OLy,bi,bj),
158       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
159       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
160          ENDIF          ENDIF
# Line 204  C Note: will inline this S/R once thSIce Line 163  C Note: will inline this S/R once thSIce
163          IF ( useEXF ) THEN          IF ( useEXF ) THEN
164           CALL THSICE_MAP_EXF(           CALL THSICE_MAP_EXF(
165       I                  iceMask,       I                  iceMask,
166       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
167         O                  snowPrc(1-OLx,1-OLy,bi,bj),
168       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
169       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
170          ENDIF          ENDIF
171  #endif  #endif
172    
173          CALL THSICE_STEP_TEMP(  #ifdef ALLOW_AUTODIFF_TAMC
174    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
175    CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
176    CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
177    #else
178            IF ( .NOT.thSIce_skipThermo ) THEN
179    #endif
180              CALL THSICE_STEP_TEMP(
181       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
182       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  CADJ STORE empmr, qnet = comlev1, key = iicekey  CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
186  CADJ STORE iceMask = comlev1, key = iicekey  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
187  CADJ STORE iceHeight  = comlev1, key = iicekey  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
188  CADJ STORE snowHeight = comlev1, key = iicekey  CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
189  CADJ STORE Tsrf    = comlev1, key = iicekey  CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
190  CADJ STORE Qice1   = comlev1, key = iicekey  cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
191  CADJ STORE Qice2   = comlev1, key = iicekey  CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
192  CADJ STORE snowAge = comlev1, key = iicekey  CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
193    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
194    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
195    #else
196            ENDIF
197            IF ( .NOT.thSIce_skipThermo ) THEN
198  #endif  #endif
199              CALL THSICE_STEP_FWD(
         CALL THSICE_STEP_FWD(  
200       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
201       I                     prcAtm,       I                     prcAtm(1-OLx,1-OLy,bi,bj),
202       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
203    #ifndef ALLOW_AUTODIFF_TAMC
204            ENDIF
205    #endif
206    
207    C--  end bi,bj loop
208           ENDDO
209          ENDDO
210    
211    #ifdef ALLOW_BALANCE_FLUXES
212    C--   Balance net Fresh-Water flux from Atm+Land
213          IF ( thSIceBalanceAtmFW.NE.0 ) THEN
214            CALL THSICE_BALANCE_FRW(
215         I                      iMin, iMax, jMin, jMax,
216         I                      prcAtm, myTime, myIter, myThid )
217          ENDIF
218    #endif
219    
220          DO bj=myByLo(myThid),myByHi(myThid)
221           DO bi=myBxLo(myThid),myBxHi(myThid)
222          CALL THSICE_AVE(          CALL THSICE_AVE(
223       I                     bi,bj, myTime, myIter, myThid )       I                     bi,bj, myTime, myIter, myThid )
224           ENDDO
225          ENDDO
226    
227    C     add a small piece of code to check AddFluid implementation:
228    c#include "thsice_test_addfluid.h"
229    
230    C--   Exchange fields that are advected by seaice dynamics
231          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
232         &       .OR. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
233            CALL THSICE_DO_EXCH( myThid )
234          ENDIF
235    #ifdef OLD_THSICE_CALL_SEQUENCE
236    #ifdef ATMOSPHERIC_LOADING
237          IF ( useRealFreshWaterFlux .AND.
238         &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
239         &  _EXCH_XY_RS( sIceLoad, myThid )
240    #endif
241    #else /* OLD_THSICE_CALL_SEQUENCE */
242    #ifdef ATMOSPHERIC_LOADING
243          IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
244         &                           .AND. thSIceAdvScheme.LE.0 )
245         &  _EXCH_XY_RS( sIceLoad, myThid )
246    #endif
247    
248  c      ENDDO  C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
249  c     ENDDO  C     otherwise, call it from here, after thsice-thermodynamics is done
250          IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
251             CALL THSICE_DO_ADVECT(
252         I                          0, 0, myTime, myIter, myThid )
253          ENDIF
254    #endif /* OLD_THSICE_CALL_SEQUENCE */
255    
256          DO bj=myByLo(myThid),myByHi(myThid)
257           DO bi=myBxLo(myThid),myBxHi(myThid)
258  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
259  C--   and stressReduction is always set to zero  C--   and stressReduction is always set to zero
260  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
261  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
262  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
263  #endif  #endif
264          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
265            DO j = jMin, jMax            DO j = jMin, jMax
# Line 263  C--  end bi,bj loop Line 282  C--  end bi,bj loop
282         ENDDO         ENDDO
283        ENDDO        ENDDO
284    
       IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN  
 C--   Exchange fields that are advected by seaice dynamics  
         _EXCH_XY_R8( iceMask, myThid )  
         _EXCH_XY_R8( iceHeight, myThid )  
         _EXCH_XY_R8( snowHeight, myThid )  
         _EXCH_XY_R8( Qice1, myThid )  
         _EXCH_XY_R8( Qice2, myThid )  
 #ifdef ATMOSPHERIC_LOADING  
         IF (useRealFreshWaterFlux)  
      &  _EXCH_XY_RS( sIceLoad, myThid )  
 #endif  
       ENDIF  
   
285  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
286  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
287    
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22