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