/[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.33 by jmc, Sat Apr 20 23:59:07 2013 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "THSICE_OPTIONS.h"  #include "THSICE_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF_TAMC
6    # ifdef ALLOW_EXF
7    #  include "EXF_OPTIONS.h"
8    # endif
9    #endif
10    
11  CBOP  CBOP
12  C     !ROUTINE: THSICE_MAIN  C     !ROUTINE: THSICE_MAIN
# Line 23  C     === Global variables === Line 28  C     === Global variables ===
28  #include "SIZE.h"  #include "SIZE.h"
29  #include "EEPARAMS.h"  #include "EEPARAMS.h"
30  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
 #include "SURFACE.h"  
 #include "DYNVARS.h"  
31  #include "FFIELDS.h"  #include "FFIELDS.h"
32  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
33    #include "THSICE_SIZE.h"
34  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
35  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
36    # include "THSICE_TAVE.h"
37    # include "THSICE_COST.h"
38    # include "DYNVARS.h"
39  # include "tamc.h"  # include "tamc.h"
40  # include "tamc_keys.h"  # include "tamc_keys.h"
41    # ifdef ALLOW_EXF
42    #  include "EXF_FIELDS.h"
43    #  include "EXF_PARAM.h"
44    #  include "EXF_CONSTANTS.h"
45    # endif /* ALLOW_EXF */
46  #endif  #endif
47    
48  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
49  C     === Routine arguments ===  C     === Routine arguments ===
50  C     myIter :: iteration counter for this thread  C     myTime    :: Current time in simulation (s)
51  C     myTime :: time counter for this thread  C     myIter    :: Current iteration number
52  C     myThid :: thread number for this instance of the routine.  C     myThid    :: My Thread Id. number
53        _RL  myTime        _RL     myTime
54        INTEGER myIter        INTEGER myIter
55        INTEGER myThid        INTEGER myThid
56  CEOP  CEOP
# Line 47  CEOP Line 58  CEOP
58  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
59  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
60  C     === Local variables ===  C     === Local variables ===
61    C     qPrcRn    :: Energy content of Precip+RunOff (+=down) [W/m2]
62        INTEGER i,j        INTEGER i,j
63        INTEGER bi,bj        INTEGER bi,bj
64        INTEGER iMin, iMax        INTEGER iMin, iMax
65        INTEGER jMin, jMax        INTEGER jMin, jMax
66        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67          _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69  c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70  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 73  c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy
73    
74  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
75    
76        IF ( stressReduction.GT. 0. _d 0 ) THEN        IF ( useEXF .OR. useSEAICE ) THEN
77    C-    EXF does not provide valid fields in overlap
78           iMin = 1
79           iMax = sNx
80           jMin = 1
81           jMax = sNy
82          ELSEIF ( stressReduction.GT. 0. _d 0 ) THEN
83  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
84         iMin = 1-OLx         iMin = 1-OLx
85         iMax = sNx+OLx-1         iMax = sNx+OLx-1
# Line 93  C      to be valid at the boundaries ( d Line 112  C      to be valid at the boundaries ( d
112            act3 = myThid - 1            act3 = myThid - 1
113            max3 = nTx*nTy            max3 = nTx*nTy
114            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
115            iicekey = (act1 + 1) + act2*max1            ticekey = (act1 + 1) + act2*max1
116       &                         + act3*max1*max2       &                         + act3*max1*max2
117       &                         + act4*max1*max2*max3       &                         + act4*max1*max2*max3
118  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
119    
 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  
   
120  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
121  CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
122  CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
123    CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
124    # ifdef ALLOW_EXF
125    CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
126    # endif
127    #endif
128    #ifdef ALLOW_AUTODIFF_TAMC
129    CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
130    CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
131  #endif  #endif
132    
133           DO j = jMin, jMax          DO j=1-OLy,sNy+OLy
134            DO i = iMin, iMax           DO i=1-OLx,sNx+OLx
135             tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)            prcAtm  (i,j,bi,bj) = 0. _d 0
136             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  
137           ENDDO           ENDDO
138            ENDDO
139    
140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_CHEAPAML
141  CADJ STORE iceMask = comlev1, key = iicekey          IF ( .NOT.useCheapAML ) THEN
142  CADJ STORE iceHeight  = comlev1, key = iicekey  #endif
 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  
143    
144  CADJ STORE empmr   = comlev1, key = iicekey          CALL THSICE_GET_OCEAN(
145  CADJ STORE qnet    = comlev1, key = iicekey       I                         bi, bj, myTime, myIter, myThid )
146    
147    #ifdef ALLOW_AUTODIFF_TAMC
148    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
149    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
150    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
151    CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
152    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
153    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
154    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
155    CADJ STORE snowPrc(:,:,bi,bj)  = comlev1_bibj, key = ticekey
156    
157    CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
158    CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
159    CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
160    CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
161  #endif  #endif
162    
163    #ifdef OLD_THSICE_CALL_SEQUENCE
164    C-      do sea-ice advection before getting surface fluxes
165    C Note: will inline this S/R once thSIce in Atmos. set-up is settled
166            IF ( thSIceAdvScheme.GT.0 )
167         &   CALL THSICE_DO_ADVECT(
168         I                   bi,bj, myTime, myIter, myThid )
169    #endif /* OLD_THSICE_CALL_SEQUENCE */
170    
171  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
172           IF ( useBulkforce ) THEN          IF ( useBulkforce ) THEN
173             CALL THSICE_GET_PRECIP(           CALL THSICE_GET_PRECIP(
174       I                  iceMask,       I                  iceMask,
175       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
176         O                  snowPrc(1-OLx,1-OLy,bi,bj),
177       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
178       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
179           ENDIF          ENDIF
180  #endif  #endif
181  #ifdef ALLOW_EXF  #ifdef ALLOW_EXF
182           IF ( useEXF ) THEN          IF ( useEXF ) THEN
183             CALL THSICE_MAP_EXF(           CALL THSICE_MAP_EXF(
184       I                  iceMask,       I                  iceMask, tOceMxL,
185       O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),       O                  prcAtm(1-OLx,1-OLy,bi,bj),
186         O                  snowPrc(1-OLx,1-OLy,bi,bj), qPrcRn,
187       O                  icFlxSW(1-OLx,1-OLy,bi,bj),       O                  icFlxSW(1-OLx,1-OLy,bi,bj),
188       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )       I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
189           ENDIF          ENDIF
190  #endif  #endif
191    
192    #ifdef ALLOW_AUTODIFF_TAMC
193           CALL THSICE_STEP_TEMP(  CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
194    CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
195    CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
196    #else
197            IF ( .NOT.thSIce_skipThermo ) THEN
198    #endif
199              CALL THSICE_STEP_TEMP(
200       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
201       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
202    
203           CALL THSICE_STEP_FWD(  #ifdef ALLOW_CHEAPAML
204            ENDIF
205    #endif
206    #ifdef ALLOW_AUTODIFF_TAMC
207    CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
208    CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
209    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
210    CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
211    CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
212    cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
213    CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
214    CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
215    CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
216    CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
217    #else
218            ENDIF
219            IF ( .NOT.thSIce_skipThermo ) THEN
220    #endif
221              CALL THSICE_STEP_FWD(
222       I                     bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, iMin, iMax, jMin, jMax,
223       I                     prcAtm,       I                     prcAtm(1-OLx,1-OLy,bi,bj), qPrcRn,
224       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
225    #ifndef ALLOW_AUTODIFF_TAMC
226            ENDIF
227    #endif
228    
229           CALL THSICE_AVE(  C--  end bi,bj loop
230       I                     bi,bj, myTime, myIter, myThid )         ENDDO
231          ENDDO
232    
233  c      ENDDO  #ifdef ALLOW_BALANCE_FLUXES
234  c     ENDDO  C--   Balance net Fresh-Water flux from Atm+Land
235          IF ( thSIceBalanceAtmFW.NE.0 ) THEN
236            CALL THSICE_BALANCE_FRW(
237         I                      iMin, iMax, jMin, jMax,
238         I                      prcAtm, myTime, myIter, myThid )
239          ENDIF
240    #endif
241    
242  c       IF ( .FALSE. ) THEN  C     add a small piece of code to check AddFluid implementation:
243    c#include "thsice_test_addfluid.h"
244    
245    C--   Exchange fields that are advected by seaice dynamics
246          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
247         &       .OR. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
248            CALL THSICE_DO_EXCH( myThid )
249          ENDIF
250    #ifdef OLD_THSICE_CALL_SEQUENCE
251    #ifdef ATMOSPHERIC_LOADING
252          IF ( useRealFreshWaterFlux .AND.
253         &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
254         &  _EXCH_XY_RS( sIceLoad, myThid )
255    #endif
256    #else /* OLD_THSICE_CALL_SEQUENCE */
257    #ifdef ATMOSPHERIC_LOADING
258          IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
259         &                           .AND. thSIceAdvScheme.LE.0 )
260         &  _EXCH_XY_RS( sIceLoad, myThid )
261    #endif
262    
263    C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
264    C     otherwise, call it from here, after thsice-thermodynamics is done
265          IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
266             CALL THSICE_DO_ADVECT(
267         I                          0, 0, myTime, myIter, myThid )
268          ENDIF
269    #endif /* OLD_THSICE_CALL_SEQUENCE */
270    
271          DO bj=myByLo(myThid),myByHi(myThid)
272           DO bi=myBxLo(myThid),myBxHi(myThid)
273    C--   Cumulate time-averaged fields and also fill-up flux diagnostics
274    C     (if not done in THSICE_DO_ADVECT call)
275    #ifdef OLD_THSICE_CALL_SEQUENCE
276            IF ( .TRUE. ) THEN
277    #else /* OLD_THSICE_CALL_SEQUENCE */
278            IF ( thSIceAdvScheme.LE.0 ) THEN
279    #endif /* OLD_THSICE_CALL_SEQUENCE */
280             CALL THSICE_AVE(
281         I                     bi,bj, myTime, myIter, myThid )
282            ENDIF
283    C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
284    C--   and stressReduction is always set to zero
285  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
286  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
287  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte  CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
288  #endif  #endif
289          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
290           DO j = jMin, jMax            DO j = jMin, jMax
291            DO i = iMin+1,iMax             DO i = iMin+1,iMax
292              tauFac = stressReduction              tauFac = stressReduction
293       &             *(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
294              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)
295               ENDDO
296            ENDDO            ENDDO
297           ENDDO            DO j = jMin+1, jMax
298           DO j = jMin+1, jMax             DO i = iMin, iMax
           DO i = iMin, iMax  
299              tauFac = stressReduction              tauFac = stressReduction
300       &             *(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
301              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)
302               ENDDO
303            ENDDO            ENDDO
          ENDDO  
304          ENDIF          ENDIF
305    
306  C--  end bi,bj loop  C--  end bi,bj loop
307         ENDDO         ENDDO
308        ENDDO        ENDDO
309    
 #ifdef ATMOSPHERIC_LOADING  
 c     IF (useRealFreshWaterFlux) _EXCH_XY_RS(sIceLoad, myThid)  
 #endif  
   
310  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
311  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
312    
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22