/[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.7 by jmc, Thu Jul 22 22:52:59 2004 UTC revision 1.14 by jmc, Tue Mar 14 15:58:27 2006 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: THSICE_STEP_FWD  C     !ROUTINE: THSICE_STEP_FWD
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE THSICE_STEP_FWD(        SUBROUTINE THSICE_STEP_FWD(
10       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
11       I             prcAtm,       I             prcAtm,
12       U             evpAtm, flxSW,       U             evpAtm, flxSW,
13       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | S/R  THSICE_STEP_FWD              C     | S/R  THSICE_STEP_FWD
17  C     | o Step Forward Therm-SeaIce model.  C     | o Step Forward Therm-SeaIce model.
18  C     *==========================================================*  C     *==========================================================*
19  C     \ev  C     \ev
# Line 30  C     === Global variables === Line 30  C     === Global variables ===
30  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
31  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
32  #include "THSICE_TAVE.h"  #include "THSICE_TAVE.h"
33    
34  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
35  C     === Routine arguments ===  C     === Routine arguments ===
36  C     bi,bj   :: tile indices  C     bi,bj   :: tile indices
37  C   iMin,iMax :: computation domain: 1rst index range  C   iMin,iMax :: computation domain: 1rst index range
38  C   jMin,jMax :: computation domain: 2nd  index range  C   jMin,jMax :: computation domain: 2nd  index range
39  C- input:  C- input:
40  C     prcAtm  :: total precip from the atmosphere [kg/m2/s]  C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
41  C     evpAtm  :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)  C     evpAtm  :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
42  C     flxSW   :: (Inp) short-wave heat flux (+=down): downward comp. only  C     flxSW   :: (Inp) short-wave heat flux (+=down): downward comp. only
43  C                      (part.1), becomes net SW flux into ocean (part.2).  C                      (part.1), becomes net SW flux into ocean (part.2).
# Line 73  C     frw2oc    :: fresh-water flux from Line 73  C     frw2oc    :: fresh-water flux from
73  C     fsalt     :: mass salt flux to the ocean  C     fsalt     :: mass salt flux to the ocean
74  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
75  C     TFrzOce   :: sea-water freezing temperature [oC] (function of S)  C     TFrzOce   :: sea-water freezing temperature [oC] (function of S)
76    C     isIceFree :: true for ice-free grid-cell that remains ice-free
77        INTEGER i,j        INTEGER i,j
78        _RL snowPr        _RL snowPr
79        _RL agingTime, ageFac        _RL agingTime, ageFac
# Line 88  C     TFrzOce   :: sea-water freezing te Line 89  C     TFrzOce   :: sea-water freezing te
89        _RL oceV2s, oceTs        _RL oceV2s, oceTs
90        _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)        _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)
91        _RL tmpflx(0:2), tmpdTs        _RL tmpflx(0:2), tmpdTs
92          LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93    #ifdef ALLOW_DIAGNOSTICS
94          _RL tmpFac
95    #endif
96    
97        LOGICAL dBug        LOGICAL dBug
98    
99   1010 FORMAT(A,1P4E11.3)   1010 FORMAT(A,1P4E14.6)
100        dBug = .FALSE.        dBug = .FALSE.
101  C-    Initialise flxAtm  C-    Initialise flxAtm
102         DO j = 1-Oly, sNy+Oly         DO j = 1-Oly, sNy+Oly
103          DO i = 1-Olx, sNx+Olx          DO i = 1-Olx, sNx+Olx
104            flxAtm(i,j) = 0.            flxAtm(i,j) = 0.
105              isIceFree(i,j) = .FALSE.
106          ENDDO          ENDDO
107         ENDDO         ENDDO
108    
109        IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN        IF ( fluidIsWater ) THEN
110         DO j = jMin, jMax         DO j = jMin, jMax
111          DO i = iMin, iMax          DO i = iMin, iMax
112  c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
113    
114  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
115  C    part.1 : ice-covered fraction ;  C    part.1 : ice-covered fraction ;
116  C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes  C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes
117  C-------  C-------
118           IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN           IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
119            icFrac  = iceMask(i,j,bi,bj)            icFrac  = iceMask(i,j,bi,bj)
120            TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)            TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
121            hIce    = iceHeight(i,j,bi,bj)            hIce    = iceHeight(i,j,bi,bj)
# Line 128  C------- Line 134  C-------
134       O               albedo,       O               albedo,
135       I               myThid )       I               myThid )
136            flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)            flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)
137              siceAlb(i,j,bi,bj) = albedo
138    
139            CALL THSICE_SOLVE4TEMP(            CALL THSICE_SOLVE4TEMP(
140       I          useBulkforce, tmpflx, TFrzOce, hIce, hSnow,       I          useBulkForce, tmpflx, TFrzOce, hIce, hSnow,
141       U          flxSW(i,j), Tsf, qicen,       U          flxSW(i,j), Tsf, qicen,
142       O          Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),       O          Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
143       O          tmpdTs, flxAtm(i,j), evpAtm(i,j),       O          tmpdTs, flxAtm(i,j), evpAtm(i,j),
# Line 147  C--    Update Sea-Ice state : Line 154  C--    Update Sea-Ice state :
154            Tice2(i,j,bi,bj)=Tice(2)            Tice2(i,j,bi,bj)=Tice(2)
155            Qice1(i,j,bi,bj)=qicen(1)            Qice1(i,j,bi,bj)=qicen(1)
156            Qice2(i,j,bi,bj)=qicen(2)            Qice2(i,j,bi,bj)=qicen(2)
 #ifdef ALLOW_TIMEAVE  
           ice_albedo_Ave(i,j,bi,bj) = ice_albedo_Ave(i,j,bi,bj)  
      &                              + icFrac*albedo*thSIce_deltaT  
 #endif /*ALLOW_TIMEAVE*/  
157            IF ( dBug ) THEN            IF ( dBug ) THEN
158             WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =',             WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =',
159       &                              Tsf, Tice, frzmltMxL       &                              Tsf, Tice, frzmltMxL
# Line 164  C--    Update Sea-Ice state : Line 167  C--    Update Sea-Ice state :
167        ENDIF        ENDIF
168        dBug = .FALSE.        dBug = .FALSE.
169    
170    #ifdef ALLOW_DIAGNOSTICS
171          IF ( useDiagnostics ) THEN
172            tmpFac = 1. _d 0
173            CALL DIAGNOSTICS_FRACT_FILL(
174         I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
175         I                   0,1,1,bi,bj,myThid)
176            CALL DIAGNOSTICS_FRACT_FILL(
177         I                   siceAlb,   iceMask,tmpFac,1,'SIalbedo',
178         I                   0,1,1,bi,bj,myThid)
179          ENDIF
180    #endif /* ALLOW_DIAGNOSTICS */
181          DO j = jMin, jMax
182           DO i = iMin, iMax
183              siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
184           ENDDO
185          ENDDO
186    
187  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
188  C    part.2 : ice-covered fraction ;  C    part.2 : ice-covered fraction ;
189  C     change in ice/snow thickness and ice-fraction  C     change in ice/snow thickness and ice-fraction
190  C     note: can only reduce the ice-fraction but not increase it.  C     note: can only reduce the ice-fraction but not increase it.
191  C-------  C-------
192        agingTime = 50. _d 0 * 86400. _d 0        agingTime = 50. _d 0 * 86400. _d 0
# Line 187  C------- Line 207  C-------
207          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN
208            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
209            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
210       &                  compact, iceHeight(i,j,bi,bj),       &                  compact, iceHeight(i,j,bi,bj),
211       &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)       &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
212            WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',            WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',
213       &                     oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)       &                     oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)
214          ENDIF          ENDIF
215  C-------  C-------
216          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
217    
218            oceV2s  = v2ocMxL(i,j,bi,bj)            oceV2s  = v2ocMxL(i,j,bi,bj)
219            snowPr  = snowPrc(i,j,bi,bj)            snowPr  = snowPrc(i,j,bi,bj)
# Line 216  C         but to reproduce old results, Line 236  C         but to reproduce old results,
236            snowPrc(i,j,bi,bj) = snowPr            snowPrc(i,j,bi,bj) = snowPr
237    
238  C--  Snow aging :  C--  Snow aging :
239            snowAge(i,j,bi,bj) = thSIce_deltaT            snowAge(i,j,bi,bj) = thSIce_deltaT
240       &                       + snowAge(i,j,bi,bj)*ageFac       &                       + snowAge(i,j,bi,bj)*ageFac
241            IF ( snowPr.GT.0. _d 0 )            IF ( snowPr.GT.0. _d 0 )
242       &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)       &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
243       &          * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )       &          * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )
244  C--  C--
245    
# Line 295  c       compact= iceMask(i,j,bi,bj) Line 315  c       compact= iceMask(i,j,bi,bj)
315  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='
317       &                 ,compact,flx2oc,fsalt,frw2oc       &                 ,compact,flx2oc,fsalt,frw2oc
318  #ifdef CHECK_ENERGY_CONSERV  #ifdef CHECK_ENERGY_CONSERV
319            tmpflx(1) = 0.            tmpflx(1) = 0.
320            tmpflx(2) = 0.            tmpflx(2) = 0.
321            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,
322       I            icFrac, compact, hIce, hSnow, qicen,       I            icFrac, compact, hIce, hSnow, qicen,
323       I            flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),       I            flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),
# Line 316  C--    Update Sea-Ice state : Line 336  C--    Update Sea-Ice state :
336            snowHeight(i,j,bi,bj)= hSnow            snowHeight(i,j,bi,bj)= hSnow
337  C--    Net fluxes :  C--    Net fluxes :
338            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc
339            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw
340            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt
341    
342            IF (dBug) WRITE(6,1010)            IF (dBug) WRITE(6,1010)
# Line 325  C--    Net fluxes : Line 345  C--    Net fluxes :
345  C--   - if esurp > 0 : end  C--   - if esurp > 0 : end
346          ENDIF          ENDIF
347    
348            IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
349         &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
350         &                       .AND.   compact .LE.0. _d 0
351          IF ( compact .GT. 0. _d 0 ) THEN          IF ( compact .GT. 0. _d 0 ) THEN
352            iceMask(i,j,bi,bj)=compact            iceMask(i,j,bi,bj)=compact
353            IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0            IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0
# Line 354  C--     Compute Sea-Ice Loading (= mass Line 377  C--     Compute Sea-Ice Loading (= mass
377         ENDDO         ENDDO
378        ENDDO        ENDDO
379    
380    #ifdef ALLOW_BULK_FORCE
381          IF ( useBulkForce ) THEN
382            CALL BULKF_FLUX_ADJUST(
383         I                          bi, bj, iMin, iMax, jMin, jMax,
384         I                          isIceFree, myTime, myIter, myThid )
385          ENDIF
386    #endif /* ALLOW_BULK_FORCE */
387    
388  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
389  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
390    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22