/[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.9 by mlosch, Tue May 30 22:49:00 2006 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
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
37  C     === Routine arguments ===  C     === Routine arguments ===
38  C     myIter :: iteration counter for this thread  C     myTime    :: Current time in simulation (s)
39  C     myTime :: time counter for this thread  C     myIter    :: Current iteration number
40  C     myThid :: thread number for this instance of the routine.  C     myThid    :: My Thread Id. number
41        _RL  myTime        _RL     myTime
42        INTEGER myIter        INTEGER myIter
43        INTEGER myThid        INTEGER myThid
44  CEOP  CEOP
# 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
# 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    
 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  
   
108  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
109  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
110  CADJ STORE vvel (:,:,1,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=ticekey, byte=isbyte
112    # ifdef ALLOW_EXF
113    CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
114    # endif
115    #endif
116    #ifdef ALLOW_AUTODIFF_TAMC
117    CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
118    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  
           ENDDO  
125           ENDDO           ENDDO
126            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 sHeating = comlev1, key = iicekey  
 CADJ STORE flxCndBt = 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
149    C Note: will inline this S/R once thSIce in Atmos. set-up is settled
150            IF ( thSIceAdvScheme.GT.0 )
151         &   CALL THSICE_DO_ADVECT(
152         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
164  #endif  #endif
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    #ifdef ALLOW_AUTODIFF_TAMC
177           CALL THSICE_STEP_TEMP(  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           CALL THSICE_STEP_FWD(  #ifdef ALLOW_AUTODIFF_TAMC
188    CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
189    CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
190    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
191    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
192    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
193    cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
194    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
195    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
202              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    #ifdef ALLOW_BALANCE_FLUXES
215    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      ENDDO  C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
245  c     ENDDO  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  c       IF ( .FALSE. ) THEN        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,
265    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
272            DO i = iMin+1,iMax             DO i = iMin+1,iMax
273              tauFac = stressReduction              tauFac = stressReduction
274       &             *(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
275              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)
276               ENDDO
277            ENDDO            ENDDO
278           ENDDO            DO j = jMin+1, jMax
279           DO j = jMin+1, jMax             DO i = iMin, iMax
           DO i = iMin, iMax  
280              tauFac = stressReduction              tauFac = stressReduction
281       &             *(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
282              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)
283               ENDDO
284            ENDDO            ENDDO
          ENDDO  
285          ENDIF          ENDIF
286    
287  C--  end bi,bj loop  C--  end bi,bj loop
288         ENDDO         ENDDO
289        ENDDO        ENDDO
290    
 #ifdef ATMOSPHERIC_LOADING  
 c     IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)  
 #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.9  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22