/[MITgcm]/MITgcm/pkg/thsice/thsice_step_fwd.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_step_fwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.18 by jmc, Wed Apr 4 02:40:42 2007 UTC revision 1.23 by jscott, Fri Apr 20 19:24:47 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    #ifdef ALLOW_ATM2D
6    # include "ctrparam.h"
7    #endif
8    
9  CBOP  CBOP
10  C     !ROUTINE: THSICE_STEP_FWD  C     !ROUTINE: THSICE_STEP_FWD
# Line 25  C     === Global variables === Line 28  C     === Global variables ===
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "FFIELDS.h"  #include "FFIELDS.h"
31    #ifdef  ALLOW_ATM2D
32    # include "ATMSIZE.h"
33    # include "ATM2D_VARS.h"
34    #endif
35  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
36  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
37  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
38  #include "THSICE_TAVE.h"  #include "THSICE_TAVE.h"
39  #include "THSICE_2DYN.h"  #include "THSICE_2DYN.h"
40    #ifdef ALLOW_AUTODIFF_TAMC
41    # include "tamc.h"
42    # include "tamc_keys.h"
43    #endif
44    
45        INTEGER siLo, siHi, sjLo, sjHi        INTEGER siLo, siHi, sjLo, sjHi
46        PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )        PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
47        PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )        PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
# Line 97  C-    define grid-point location where t Line 109  C-    define grid-point location where t
109    
110  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
111    
112    #ifdef ALLOW_AUTODIFF_TAMC
113          act1 = bi - myBxLo(myThid)
114          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
115          act2 = bj - myByLo(myThid)
116          max2 = myByHi(myThid) - myByLo(myThid) + 1
117          act3 = myThid - 1
118          max3 = nTx*nTy
119          act4 = ikey_dynamics - 1
120          iicekey = (act1 + 1) + act2*max1
121         &                     + act3*max1*max2
122         &                     + act4*max1*max2*max3
123    #endif /* ALLOW_AUTODIFF_TAMC */
124    
125  C-    Initialise  C-    Initialise
126        dBugFlag = debugLevel.GE.debLevB        dBugFlag = debugLevel.GE.debLevB
127        DO j = 1-OLy, sNy+OLy        DO j = 1-OLy, sNy+OLy
128          DO i = 1-OLx, sNx+OLx          DO i = 1-OLx, sNx+OLx
129            isIceFree(i,j) = .FALSE.            isIceFree(i,j) = .FALSE.
130    #ifdef ALLOW_ATM2D
131              sFluxFromIce(i,j) = 0. _d 0
132    #else
133            saltFlux(i,j,bi,bj) = 0. _d 0            saltFlux(i,j,bi,bj) = 0. _d 0
134    #endif
135  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
136            iceFrac(i,j) = 0.            iceFrac(i,j) = 0.
137  #endif  #endif
# Line 111  C-    Initialise Line 140  C-    Initialise
140    
141        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
142        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
143    
144    #ifdef ALLOW_AUTODIFF_TAMC
145    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
146    #endif
147        DO j = jMin, jMax        DO j = jMin, jMax
148         DO i = iMin, iMax         DO i = iMin, iMax
149          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
# Line 135  C-- Line 168  C--
168  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
169        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
170          tmpFac = 1. _d 0          tmpFac = 1. _d 0
171            CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
172          CALL DIAGNOSTICS_FRACT_FILL(          CALL DIAGNOSTICS_FRACT_FILL(
173       I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',       I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
174       I                   0,1,1,bi,bj,myThid)       I                   0,1,1,bi,bj,myThid)
# Line 184  C------- Line 218  C-------
218         ENDDO         ENDDO
219        ENDDO        ENDDO
220    
221    #ifdef ALLOW_AUTODIFF_TAMC
222    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
223    #endif
224    
225        CALL THSICE_CALC_THICKN(        CALL THSICE_CALC_THICKN(
226       I          bi, bj, siLo, siHi, sjLo, sjHi,       I          bi, bj, siLo, siHi, sjLo, sjHi,
227       I          iMin,iMax, jMin,jMax, dBugFlag,       I          iMin,iMax, jMin,jMax, dBugFlag,
# Line 201  C------- Line 239  C-------
239  C--    Net fluxes :  C--    Net fluxes :
240        DO j = jMin, jMax        DO j = jMin, jMax
241         DO i = iMin, iMax         DO i = iMin, iMax
242    #ifdef ALLOW_AUTODIFF_TAMC
243              ikey_1 = i
244         &         + sNx*(j-1)
245         &         + sNx*sNy*act1
246         &         + sNx*sNy*max1*act2
247         &         + sNx*sNy*max1*max2*act3
248         &         + sNx*sNy*max1*max2*max3*act4
249    #endif /* ALLOW_AUTODIFF_TAMC */
250    C--
251    #ifdef ALLOW_AUTODIFF_TAMC
252    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
253    #endif
254          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
255  C-     weighted average net fluxes:  C-     weighted average net fluxes:
256    #ifdef ALLOW_AUTODIFF_TAMC
257    CADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
258    CADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
259    CADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
260    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
261    #endif
262            icFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
263            opFrac= 1. _d 0-icFrac            opFrac= 1. _d 0-icFrac
264    #ifdef ALLOW_ATM2D
265              pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
266              pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
267              sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
268    #else
269            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
270       &                        - opFrac*Qnet(i,j,bi,bj)       &                        - opFrac*Qnet(i,j,bi,bj)
271            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
# Line 213  C-     weighted average net fluxes: Line 274  C-     weighted average net fluxes:
274            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
275       &                    +  opFrac*EmPmR(i,j,bi,bj)       &                    +  opFrac*EmPmR(i,j,bi,bj)
276            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
277    #endif
278    
279  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
280            IF (dBug(i,j,bi,bj)) WRITE(6,1010)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
# Line 247  C------- Line 309  C-------
309       O          flx2oc, frw2oc, fsalt,       O          flx2oc, frw2oc, fsalt,
310       I          myTime, myIter, myThid )       I          myTime, myIter, myThid )
311    
312    #ifdef ALLOW_AUTODIFF_TAMC
313    CADJ STORE snowHeight(:,:,bi,bj) =
314    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
315    #endif
316        DO j = jMin, jMax        DO j = jMin, jMax
317         DO i = iMin, iMax         DO i = iMin, iMax
318          IF (frzmltMxL(i,j).GT.0. _d 0) THEN          IF (frzmltMxL(i,j).GT.0. _d 0) THEN
319  C--    Net fluxes :  C--    Net fluxes :
320    #ifdef ALLOW_ATM2D
321              pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
322              pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
323              sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
324    #else
325            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
326            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
327            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
328    #endif
329    
330  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
331            IF (dBug(i,j,bi,bj)) WRITE(6,1010)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
# Line 278  C--    Net fluxes : Line 350  C--    Net fluxes :
350            Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)            Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
351            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
352            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
353            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = Lfresh
354            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = Lfresh
355          ENDIF          ENDIF
356           ENDDO
357          ENDDO
358    
359  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
360    # ifdef ALLOW_AUTODIFF_TAMC
361    CADJ STORE snowHeight(:,:,bi,bj) =
362    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
363    # endif
364          DO j = jMin, jMax
365           DO i = iMin, iMax
366  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
367          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
368       &                         + iceHeight(i,j,bi,bj)*rhoi       &                         + iceHeight(i,j,bi,bj)*rhoi
369       &                        )*iceMask(i,j,bi,bj)       &                        )*iceMask(i,j,bi,bj)
 #endif  
   
370         ENDDO         ENDDO
371        ENDDO        ENDDO
372    #endif
373    
374        IF ( thSIceAdvScheme.GT.0 ) THEN        IF ( thSIceAdvScheme.GT.0 ) THEN
375  C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux  C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22