/[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.30 by jmc, Thu Apr 4 00:42:06 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 51  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 60  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
# Line 93  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    
 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  
   
106  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
107  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
108  CADJ STORE vvel (:,:,1,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=ticekey, byte=isbyte
110    # ifdef ALLOW_EXF
111    CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
112    # endif
113    #endif
114    #ifdef ALLOW_AUTODIFF_TAMC
115    CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
116    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  
           ENDDO  
122           ENDDO           ENDDO
123            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 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  
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
146    C Note: will inline this S/R once thSIce in Atmos. set-up is settled
147            IF ( thSIceAdvScheme.GT.0 )
148         &   CALL THSICE_DO_ADVECT(
149         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
161  #endif  #endif
162  #ifdef ALLOW_EXF  #ifdef ALLOW_EXF
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    #ifdef ALLOW_AUTODIFF_TAMC
174           CALL THSICE_STEP_TEMP(  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           CALL THSICE_STEP_FWD(  #ifdef ALLOW_AUTODIFF_TAMC
185    CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
186    CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
187    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
188    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
189    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
190    cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
191    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
192    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
199              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           CALL THSICE_AVE(  C--  end bi,bj loop
208       I                     bi,bj, myTime, myIter, myThid )         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    C     add a small piece of code to check AddFluid implementation:
221    c#include "thsice_test_addfluid.h"
222    
223    C--   Exchange fields that are advected by seaice dynamics
224          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
225         &       .OR. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
226            CALL THSICE_DO_EXCH( myThid )
227          ENDIF
228    #ifdef OLD_THSICE_CALL_SEQUENCE
229    #ifdef ATMOSPHERIC_LOADING
230          IF ( useRealFreshWaterFlux .AND.
231         &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
232         &  _EXCH_XY_RS( sIceLoad, myThid )
233    #endif
234    #else /* OLD_THSICE_CALL_SEQUENCE */
235    #ifdef ATMOSPHERIC_LOADING
236          IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
237         &                           .AND. thSIceAdvScheme.LE.0 )
238         &  _EXCH_XY_RS( sIceLoad, myThid )
239    #endif
240    
241  c      ENDDO  C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
242  c     ENDDO  C     otherwise, call it from here, after thsice-thermodynamics is done
243          IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
244             CALL THSICE_DO_ADVECT(
245         I                          0, 0, myTime, myIter, myThid )
246          ENDIF
247    #endif /* OLD_THSICE_CALL_SEQUENCE */
248    
249  c       IF ( .FALSE. ) THEN        DO bj=myByLo(myThid),myByHi(myThid)
250           DO bi=myBxLo(myThid),myBxHi(myThid)
251    C--   Cumulate time-averaged fields and also fill-up flux diagnostics
252    C     (if not done in THSICE_DO_ADVECT call)
253    #ifdef OLD_THSICE_CALL_SEQUENCE
254            IF ( .TRUE. ) THEN
255    #else /* OLD_THSICE_CALL_SEQUENCE */
256            IF ( thSIceAdvScheme.LE.0 ) THEN
257    #endif /* OLD_THSICE_CALL_SEQUENCE */
258             CALL THSICE_AVE(
259         I                     bi,bj, myTime, myIter, myThid )
260            ENDIF
261    C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
262    C--   and stressReduction is always set to zero
263  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
264  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
265  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
266  #endif  #endif
267          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
268           DO j = jMin, jMax            DO j = jMin, jMax
269            DO i = iMin+1,iMax             DO i = iMin+1,iMax
270              tauFac = stressReduction              tauFac = stressReduction
271       &             *(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
272              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)
273               ENDDO
274            ENDDO            ENDDO
275           ENDDO            DO j = jMin+1, jMax
276           DO j = jMin+1, jMax             DO i = iMin, iMax
           DO i = iMin, iMax  
277              tauFac = stressReduction              tauFac = stressReduction
278       &             *(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
279              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)
280               ENDDO
281            ENDDO            ENDDO
          ENDDO  
282          ENDIF          ENDIF
283    
284  C--  end bi,bj loop  C--  end bi,bj loop
285         ENDDO         ENDDO
286        ENDDO        ENDDO
287    
 #ifdef ATMOSPHERIC_LOADING  
 c     IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)  
 #endif  
   
288  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
289  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
290    
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22