/[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.19 by jmc, Sun Apr 15 21:58:33 2007 UTC revision 1.26 by jmc, Mon Oct 1 13:38:08 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"  #ifdef ALLOW_AUTODIFF_TAMC
40    # include "tamc.h"
41    # include "tamc_keys.h"
42    #endif
43    
44        INTEGER siLo, siHi, sjLo, sjHi        INTEGER siLo, siHi, sjLo, sjHi
45        PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )        PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
46        PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )        PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
# Line 97  C-    define grid-point location where t Line 108  C-    define grid-point location where t
108    
109  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
110    
111    #ifdef ALLOW_AUTODIFF_TAMC
112          act1 = bi - myBxLo(myThid)
113          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
114          act2 = bj - myByLo(myThid)
115          max2 = myByHi(myThid) - myByLo(myThid) + 1
116          act3 = myThid - 1
117          max3 = nTx*nTy
118          act4 = ikey_dynamics - 1
119          iicekey = (act1 + 1) + act2*max1
120         &                     + act3*max1*max2
121         &                     + act4*max1*max2*max3
122    #endif /* ALLOW_AUTODIFF_TAMC */
123    
124  C-    Initialise  C-    Initialise
125        dBugFlag = debugLevel.GE.debLevB        dBugFlag = debugLevel.GE.debLevB
126        DO j = 1-OLy, sNy+OLy        DO j = 1-OLy, sNy+OLy
127          DO i = 1-OLx, sNx+OLx          DO i = 1-OLx, sNx+OLx
128            isIceFree(i,j) = .FALSE.            isIceFree(i,j) = .FALSE.
129    #ifdef ALLOW_ATM2D
130              sFluxFromIce(i,j) = 0. _d 0
131    #else
132            saltFlux(i,j,bi,bj) = 0. _d 0            saltFlux(i,j,bi,bj) = 0. _d 0
133    #endif
134  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
135            iceFrac(i,j) = 0.            iceFrac(i,j) = 0.
136  #endif  #endif
# Line 111  C-    Initialise Line 139  C-    Initialise
139    
140        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
141        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
142    
143    #ifdef ALLOW_AUTODIFF_TAMC
144    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
145    #endif
146        DO j = jMin, jMax        DO j = jMin, jMax
147         DO i = iMin, iMax         DO i = iMin, iMax
148          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 167  C--
167  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
168        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
169          tmpFac = 1. _d 0          tmpFac = 1. _d 0
170            CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
171          CALL DIAGNOSTICS_FRACT_FILL(          CALL DIAGNOSTICS_FRACT_FILL(
172       I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',       I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
173       I                   0,1,1,bi,bj,myThid)       I                   0,1,1,bi,bj,myThid)
# Line 184  C------- Line 217  C-------
217         ENDDO         ENDDO
218        ENDDO        ENDDO
219    
220    #ifdef ALLOW_AUTODIFF_TAMC
221    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
222    #endif
223    
224        CALL THSICE_CALC_THICKN(        CALL THSICE_CALC_THICKN(
225       I          bi, bj, siLo, siHi, sjLo, sjHi,       I          bi, bj, siLo, siHi, sjLo, sjHi,
226       I          iMin,iMax, jMin,jMax, dBugFlag,       I          iMin,iMax, jMin,jMax, dBugFlag,
# Line 201  C------- Line 238  C-------
238  C--    Net fluxes :  C--    Net fluxes :
239        DO j = jMin, jMax        DO j = jMin, jMax
240         DO i = iMin, iMax         DO i = iMin, iMax
241    #ifdef ALLOW_AUTODIFF_TAMC
242              ikey_1 = i
243         &         + sNx*(j-1)
244         &         + sNx*sNy*act1
245         &         + sNx*sNy*max1*act2
246         &         + sNx*sNy*max1*max2*act3
247         &         + sNx*sNy*max1*max2*max3*act4
248    #endif /* ALLOW_AUTODIFF_TAMC */
249    C--
250    #ifdef ALLOW_AUTODIFF_TAMC
251    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
252    #endif
253          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
254  C-     weighted average net fluxes:  C-     weighted average net fluxes:
255    #ifdef ALLOW_AUTODIFF_TAMC
256    CADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
257    CADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
258    CADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
259    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
260    #endif
261            icFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
262            opFrac= 1. _d 0-icFrac            opFrac= 1. _d 0-icFrac
263    #ifdef ALLOW_ATM2D
264              pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
265              pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
266              sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
267    #else
268            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
269       &                        - opFrac*Qnet(i,j,bi,bj)       &                        - opFrac*Qnet(i,j,bi,bj)
270            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
271       &                        + opFrac*rhofw*EmPmR(i,j,bi,bj)       &                        + opFrac*EmPmR(i,j,bi,bj)
272            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
273            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
274       &                    +  opFrac*EmPmR(i,j,bi,bj)       &                    +  opFrac*EmPmR(i,j,bi,bj)
275            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
276    #endif
277    
278  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
279            IF (dBug(i,j,bi,bj)) WRITE(6,1010)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
# Line 221  C-     weighted average net fluxes: Line 282  C-     weighted average net fluxes:
282       &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)       &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
283  #endif  #endif
284    
285          ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN          ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
286            icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
287            icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
288          ELSE          ELSE
289            icFlxAtm(i,j,bi,bj) = 0. _d 0            icFlxAtm(i,j,bi,bj) = 0. _d 0
290            icFrwAtm(i,j,bi,bj) = 0. _d 0            icFrwAtm(i,j,bi,bj) = 0. _d 0
# Line 247  C------- Line 308  C-------
308       O          flx2oc, frw2oc, fsalt,       O          flx2oc, frw2oc, fsalt,
309       I          myTime, myIter, myThid )       I          myTime, myIter, myThid )
310    
311    #ifdef ALLOW_AUTODIFF_TAMC
312    CADJ STORE snowHeight(:,:,bi,bj) =
313    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
314    #endif
315        DO j = jMin, jMax        DO j = jMin, jMax
316         DO i = iMin, iMax         DO i = iMin, iMax
317          IF (frzmltMxL(i,j).GT.0. _d 0) THEN          IF (frzmltMxL(i,j).GT.0. _d 0) THEN
318  C--    Net fluxes :  C--    Net fluxes :
319    #ifdef ALLOW_ATM2D
320              pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
321              pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
322              sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
323    #else
324            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)
325            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)
326            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)
327    #endif
328    
329  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
330            IF (dBug(i,j,bi,bj)) WRITE(6,1010)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
# Line 281  C--    Net fluxes : Line 352  C--    Net fluxes :
352            Qice1(i,j,bi,bj)    = Lfresh            Qice1(i,j,bi,bj)    = Lfresh
353            Qice2(i,j,bi,bj)    = Lfresh            Qice2(i,j,bi,bj)    = Lfresh
354          ENDIF          ENDIF
355           ENDDO
356          ENDDO
357    
358  #ifdef ATMOSPHERIC_LOADING  # ifdef ALLOW_AUTODIFF_TAMC
359    CADJ STORE snowHeight(:,:,bi,bj) =
360    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
361    # endif
362          DO j = jMin, jMax
363           DO i = iMin, iMax
364  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
365          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
366       &                         + iceHeight(i,j,bi,bj)*rhoi       &                         + iceHeight(i,j,bi,bj)*rhoi
367       &                        )*iceMask(i,j,bi,bj)       &                        )*iceMask(i,j,bi,bj)
 #endif  
   
368         ENDDO         ENDDO
369        ENDDO        ENDDO
370    
# Line 298  C--   note: those fluxes should to be ad Line 374  C--   note: those fluxes should to be ad
374          DO i = iMin, iMax          DO i = iMin, iMax
375           IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN           IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
376            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
377            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
378            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
379           ENDIF           ENDIF
380          ENDDO          ENDDO
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22