/[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.1 by jmc, Sun Nov 23 01:20:13 2003 UTC revision 1.13 by heimbach, Mon Apr 16 22:38:24 2007 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    
6    CBOP
7  C     !ROUTINE: THSICE_MAIN  C     !ROUTINE: THSICE_MAIN
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE THSICE_MAIN(        SUBROUTINE THSICE_MAIN(
10       I                        myTime, myIter, myThid )       I                        myTime, myIter, myThid )
11    C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | SUBROUTINE  THSICE_MAIN              C     | S/R  THSICE_MAIN
14  C     | o Therm_SeaIce main routine.  C     | o Therm_SeaIce main routine.
15  C     |   step forward Thermodynamic_SeaIce variables and modify  C     |   step forward Thermodynamic_SeaIce variables and modify
16  C     |    ocean surface forcing accordingly.  C     |    ocean surface forcing accordingly.
17  C     *==========================================================*  C     *==========================================================*
18    
19  C     !USES:  C     !USES:
20        IMPLICIT NONE        IMPLICIT NONE
21    
22  C     === Global variables ===  C     === Global variables ===
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26    #include "GRID.h"
27    #include "SURFACE.h"
28    #include "DYNVARS.h"
29  #include "FFIELDS.h"  #include "FFIELDS.h"
30  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
31  #include "THSICE.h"  #include "THSICE_VARS.h"
32    #ifdef ALLOW_AUTODIFF_TAMC
33    # include "tamc.h"
34    # include "tamc_keys.h"
35    C--
36    # include "THSICE_2DYN.h"
37    C--
38    #endif
39    
40  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
41  C     === Routine arguments ===  C     === Routine arguments ===
42  C     myIter :: iteration counter for this thread  C     myTime    :: Current time in simulation (s)
43  C     myTime :: time counter for this thread  C     myIter    :: Current iteration number
44  C     myThid :: thread number for this instance of the routine.  C     myThid    :: My Thread Id. number
45        _RL  myTime        _RL     myTime
46        INTEGER myIter        INTEGER myIter
47        INTEGER myThid        INTEGER myThid
48    CEOP
49    
50  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
51  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
# Line 40  C     === Local variables === Line 54  C     === Local variables ===
54        INTEGER bi,bj        INTEGER bi,bj
55        INTEGER iMin, iMax        INTEGER iMin, iMax
56        INTEGER jMin, jMax        INTEGER jMin, jMax
57          _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58    c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59    c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60    c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61    
62        _RL tauFac        _RL tauFac
63    
64  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
65    
66        IF ( stressReduction.GT. 0. _d 0 ) THEN        IF ( stressReduction.GT. 0. _d 0 ) THEN
67         iMin = 1-Olx  C-     needs new Ice Fraction in halo region to apply wind-stress reduction
68         iMax = sNx+Olx-1         iMin = 1-OLx
69         jMin = 1-Oly         iMax = sNx+OLx-1
70         jMax = sNy+Oly-1         jMin = 1-OLy
71           jMax = sNy+OLy-1
72    #ifdef ATMOSPHERIC_LOADING
73          ELSEIF ( useRealFreshWaterFlux .AND. .NOT.useSEAICE ) THEN
74    C-     needs sea-ice loading in part of the halo regions for grad.Phi0surf
75    C      to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
76           iMin = 0
77           iMax = sNx+1
78           jMin = 0
79           jMax = sNy+1
80    #endif /* ATMOSPHERIC_LOADING */
81        ELSE        ELSE
82         iMin = 1         iMin = 1
83         iMax = sNx         iMax = sNx
# Line 59  C---+----1----+----2----+----3----+----4 Line 87  C---+----1----+----2----+----3----+----4
87    
88        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
89         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
90           CALL THSICE_STEP_FWD(  
91       I                     bi, bj, iMin, iMax, jMin, jMax,  #ifdef ALLOW_AUTODIFF_TAMC
92              act1 = bi - myBxLo(myThid)
93              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
94              act2 = bj - myByLo(myThid)
95              max2 = myByHi(myThid) - myByLo(myThid) + 1
96              act3 = myThid - 1
97              max3 = nTx*nTy
98              act4 = ikey_dynamics - 1
99              iicekey = (act1 + 1) + act2*max1
100         &                         + act3*max1*max2
101         &                         + act4*max1*max2*max3
102    #endif /* ALLOW_AUTODIFF_TAMC */
103    
104    #ifdef ALLOW_AUTODIFF_TAMC
105    CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
106    CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
107    CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
108    #endif
109    
110    C--     Mixed layer thickness: take the 1rst layer
111    #ifdef NONLIN_FRSURF
112            IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
113             IF ( select_rStar.GT.0 ) THEN
114              DO j = jMin, jMax
115               DO i = iMin, iMax
116                 hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
117         &                                  *rStarFacC(i,j,bi,bj)
118               ENDDO
119              ENDDO
120             ELSE
121              DO j = jMin, jMax
122               DO i = iMin, iMax
123                IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
124                 hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
125                ELSE
126                 hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
127                ENDIF
128               ENDDO
129              ENDDO
130             ENDIF
131            ELSE
132    #else /* ndef NONLIN_FRSURF */
133            IF (.TRUE.) THEN
134    #endif /* NONLIN_FRSURF */
135              DO j = jMin, jMax
136               DO i = iMin, iMax
137                 hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
138               ENDDO
139              ENDDO
140            ENDIF
141    
142    #ifdef ALLOW_AUTODIFF_TAMC
143    CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
144    CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
145    #endif
146    
147            DO j = jMin, jMax
148             DO i = iMin, iMax
149              tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
150              sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
151              v2ocMxL(i,j,bi,bj) =
152         &              ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
153         &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
154         &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
155         &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
156         &              )*0.5 _d 0
157              prcAtm(i,j) = 0.
158              icFrwAtm(i,j,bi,bj) = 0. _d 0
159              icFlxAtm(i,j,bi,bj) = 0. _d 0
160              icFlxSW (i,j,bi,bj) = 0. _d 0
161              snowPrc(i,j,bi,bj) = 0. _d 0
162              siceAlb(i,j,bi,bj) = 0. _d 0
163             ENDDO
164            ENDDO
165    
166    #ifdef ALLOW_AUTODIFF_TAMC
167    CADJ STORE iceMask = comlev1, key = iicekey
168    CADJ STORE iceHeight  = comlev1, key = iicekey
169    CADJ STORE snowHeight = comlev1, key = iicekey
170    CADJ STORE Tsrf    = comlev1, key = iicekey
171    CADJ STORE Qice1   = comlev1, key = iicekey
172    CADJ STORE Qice2   = comlev1, key = iicekey
173    CADJ STORE snowAge = comlev1, key = iicekey
174    CADJ STORE snowPrc  = comlev1, key = iicekey
175    
176    CADJ STORE hOceMxL = comlev1, key = iicekey
177    CADJ STORE tOceMxL = comlev1, key = iicekey
178    CADJ STORE sOceMxL = comlev1, key = iicekey
179    CADJ STORE v2ocMxL = comlev1, key = iicekey
180    
181    CADJ STORE empmr   = comlev1, key = iicekey
182    CADJ STORE qnet    = comlev1, key = iicekey
183    #endif
184    
185    C-      do sea-ice advection before getting surface fluxes
186    C Note: will inline this S/R once thSIce in Atmos. set-up is settled
187            IF ( thSIceAdvScheme.GT.0 )
188         &   CALL THSICE_DO_ADVECT(
189         I                   bi,bj, myTime, myIter, myThid )
190    
191    #ifdef ALLOW_BULK_FORCE
192            IF ( useBulkforce ) THEN
193             CALL THSICE_GET_PRECIP(
194         I                  iceMask,
195         O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
196         O                  icFlxSW(1-OLx,1-OLy,bi,bj),
197         I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
198            ENDIF
199    #endif
200    #ifdef ALLOW_EXF
201            IF ( useEXF ) THEN
202             CALL THSICE_MAP_EXF(
203         I                  iceMask,
204         O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
205         O                  icFlxSW(1-OLx,1-OLy,bi,bj),
206         I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
207            ENDIF
208    #endif
209    
210    
211            CALL THSICE_STEP_TEMP(
212         I                     bi, bj, iMin, iMax, jMin, jMax,
213         I                     myTime, myIter, myThid )
214    
215    #ifdef ALLOW_AUTODIFF_TAMC
216    CADJ STORE empmr, qnet = comlev1, key = iicekey
217    CADJ STORE iceMask = comlev1, key = iicekey
218    CADJ STORE iceHeight  = comlev1, key = iicekey
219    CADJ STORE snowHeight = comlev1, key = iicekey
220    CADJ STORE Tsrf    = comlev1, key = iicekey
221    CADJ STORE Qice1   = comlev1, key = iicekey
222    CADJ STORE Qice2   = comlev1, key = iicekey
223    CADJ STORE snowAge = comlev1, key = iicekey
224    #endif
225    
226            CALL THSICE_STEP_FWD(
227         I                     bi, bj, iMin, iMax, jMin, jMax,
228         I                     prcAtm,
229       I                     myTime, myIter, myThid )       I                     myTime, myIter, myThid )
230           IF ( thSIce_taveFreq .GT. 0. _d 0 )  
231       &   CALL THSICE_AVE( bi,bj, mythid )          CALL THSICE_AVE(
232         I                     bi,bj, myTime, myIter, myThid )
233    
234  c      ENDDO  c      ENDDO
235  c     ENDDO  c     ENDDO
236    
237  c       IF ( .FALSE. ) THEN  C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
238    C--   and stressReduction is always set to zero
239    #ifdef ALLOW_AUTODIFF_TAMC
240    CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
241    CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
242    #endif
243          IF ( stressReduction.GT. 0. _d 0 ) THEN          IF ( stressReduction.GT. 0. _d 0 ) THEN
244           DO j = jMin, jMax            DO j = jMin, jMax
245            DO i = iMin+1,iMax             DO i = iMin+1,iMax
246              tauFac = stressReduction              tauFac = stressReduction
247       &             *(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
248              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)
249               ENDDO
250            ENDDO            ENDDO
251           ENDDO            DO j = jMin+1, jMax
252           DO j = jMin+1, jMax             DO i = iMin, iMax
           DO i = iMin, iMax  
253              tauFac = stressReduction              tauFac = stressReduction
254       &             *(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
255              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)
256               ENDDO
257            ENDDO            ENDDO
          ENDDO  
258          ENDIF          ENDIF
259    
260  C--  end bi,bj loop  C--  end bi,bj loop
261         ENDDO         ENDDO
262        ENDDO        ENDDO
263    
264    
265          IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN
266    C--   Exchange fields that are advected by seaice dynamics
267            _EXCH_XY_R8( iceMask, myThid )
268            _EXCH_XY_R8( iceHeight, myThid )
269            _EXCH_XY_R8( snowHeight, myThid )
270            _EXCH_XY_R8( Qice1, myThid )
271            _EXCH_XY_R8( Qice2, myThid )
272    #ifdef ATMOSPHERIC_LOADING
273            IF (useRealFreshWaterFlux)
274         &  _EXCH_XY_RS( sIceLoad, myThid )
275    #endif
276          ENDIF
277    
278  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
279  #endif  /*ALLOW_THSICE*/  #endif  /*ALLOW_THSICE*/
280    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22