/[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.1 by jmc, Sun Nov 23 01:20:13 2003 UTC revision 1.18 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
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,
12       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
13    C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | SUBROUTINE  THSICE_STEP_FWD              C     | S/R  THSICE_STEP_FWD
16  C     | o Step Forward Therm-SeaIce model.  C     | o Step Forward Therm-SeaIce model.
17  C     *==========================================================*  C     *==========================================================*
18    C     \ev
19    
20  C     !USES:  C     !USES:
21        IMPLICIT NONE        IMPLICIT NONE
22    
23  C     === Global variables ===  C     === Global variables ===
24  #include "SIZE.h"  #include "SIZE.h"
25  #include "EEPARAMS.h"  #include "EEPARAMS.h"
26  #include "PARAMS.h"  #include "PARAMS.h"
27  #include "FFIELDS.h"  #include "FFIELDS.h"
 #include "DYNVARS.h"  
 #include "GRID.h"  
 #ifdef ALLOW_BULK_FORCE  
 #include "BULKF.h"  
 #endif  
28  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
29  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
30  #include "THSICE.h"  #include "THSICE_VARS.h"
31  #include "THSICE_DIAGS.h"  #include "THSICE_TAVE.h"
32    #include "THSICE_2DYN.h"
33          INTEGER siLo, siHi, sjLo, sjHi
34          PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
35          PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
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- input:
40  C     myTime :: time counter for this thread  C     bi,bj   :: tile indices
41  C     myThid :: thread number for this instance of the routine.  C   iMin,iMax :: computation domain: 1rst index range
42    C   jMin,jMax :: computation domain: 2nd  index range
43    C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
44    C     myTime  :: current Time of simulation [s]
45    C     myIter  :: current Iteration number in simulation
46    C     myThid  :: my Thread Id number
47    C-- Use fluxes hold in commom blocks
48    C- input:
49    C     icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
50    C     icFlxAtm  :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
51    C     icFrwAtm  :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
52    C- output
53    C     icFlxAtm  :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
54    C     icFrwAtm  :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
55        INTEGER bi,bj        INTEGER bi,bj
56        INTEGER iMin, iMax        INTEGER iMin, iMax
57        INTEGER jMin, jMax        INTEGER jMin, jMax
58          _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59        _RL  myTime        _RL  myTime
60        INTEGER myIter        INTEGER myIter
61        INTEGER myThid        INTEGER myThid
62    CEOP
63    
64  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
65  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
66  C     === Local variables ===  C     === Local variables ===
67  C     Fbot          :: the oceanic heat flux already incorporated (ice_therm)  C     iceFrac   :: fraction of grid area covered in ice
68  C     flxAtm        :: net heat flux from the atmosphere ( >0 downward)  C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
69  C     evpAtm        :: evaporation to the atmosphere  C     frw2oc    :: fresh-water flux from the ice to the ocean
70  C     frwAtm        :: net fresh-water flux (E-P-R) to the atmosphere (m/s)  C     fsalt     :: mass salt flux to the ocean
71  C     qleft         :: net heat flux from the ice to the ocean  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
72  C     ffresh        :: fresh-water flux from the ice to the ocean  C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
73  C     fsalt         :: mass salt flux to the ocean  C     isIceFree :: true for ice-free grid-cell that remains ice-free
74    C     ageFac    :: snow aging factor [1]
75    C     snowFac   :: snowing refreshing-age factor [units of 1/snowPr]
76          LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RL     iceFrac  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RL     flx2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RL     frw2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RL     fsalt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL     tFrzOce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL     frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL ageFac
84          _RL snowFac
85          _RL cphm
86          _RL opFrac, icFrac
87    #ifdef ALLOW_DIAGNOSTICS
88          _RL tmpFac
89    #endif
90        INTEGER i,j        INTEGER i,j
91        _RL fswdown, qleft, qNewIce        LOGICAL dBugFlag
       _RL fsalt  
       _RL ffresh  
       _RL Tf, cphm, frzmlt  
       _RL Fbot, esurp  
       _RL flxAtm, evpAtm, frwAtm  
       _RL openFrac, iceFrac, qicAv  
       _RL oceHs, oceV2s, oceSs, oceTs  
       _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)  
92    
93        LOGICAL dBug  C-    define grid-point location where to print debugging values
94    #include "THSICE_DEBUG.h"
95    
96        dBug = .FALSE.   1010 FORMAT(A,1P4E14.6)
  1010 FORMAT(A,1P4E11.3)  
97    
98        DO j = jMin, jMax  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
        DO i = iMin, iMax  
 c       dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 )  
99    
100          Tf     = -mu_Tf*salt(i,j,1,bi,bj)  C-    Initialise
101          cphm   = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)        dBugFlag = debugLevel.GE.debLevB
102          frzmlt = (Tf-theta(i,j,1,bi,bj))*cphm/thSIce_deltaT        DO j = 1-OLy, sNy+OLy
103          Fbot   = 0. _d 0          DO i = 1-OLx, sNx+OLx
104          compact= 0. _d 0            isIceFree(i,j) = .FALSE.
105          snow(i,j,bi,bj)     = 0. _d 0            saltFlux(i,j,bi,bj) = 0. _d 0
106          saltFlux(i,j,bi,bj) = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
107              iceFrac(i,j) = 0.
108          IF (dBug.AND.(frzmlt.GT.0. .OR.iceMask(i,j,bi,bj).GT.0.)) THEN  #endif
109            WRITE(6,1010) 'ThSI_FWD:-0- iceMask,hIc,hSn,Qnet=',          ENDDO
110       &       iceMask(i,j,bi,bj),iceHeight(i,j,bi,bj),        ENDDO
      &       snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)  
           WRITE(6,1010) 'ThSI_FWD: ocTs,Tf,frzmlt=',  
      &              theta(i,j,1,bi,bj),Tf,frzmlt  
         ENDIF  
111    
112  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
113  C    part.1 : ice-covered fraction ;        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
114  C     can only reduce the ice-fraction but not increase it.        DO j = jMin, jMax
115           DO i = iMin, iMax
116            IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
117    C--  Snow aging :
118              snowAge(i,j,bi,bj) = thSIce_deltaT
119         &                       + snowAge(i,j,bi,bj)*ageFac
120              IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
121         &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
122         &          * EXP( - snowFac*snowPrc(i,j,bi,bj) )
123    c    &          * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
124    c    &                  /hNewSnowAge )
125  C-------  C-------
         IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN  
           fswdown = solar(i,j,bi,bj)  
           oceHs   = hfacC(i,j,1,bi,bj)*drF(1)  
           oceTs   = theta(i,j,1,bi,bj)  
           oceSs   = salt (i,j,1,bi,bj)  
           oceV2s  = ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)  
      &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)  
      &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)  
      &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj) )*0.5 _d 0  
           compact = iceMask(i,j,bi,bj)  
           hIce    = iceHeight(i,j,bi,bj)  
           hSnow   = snowHeight(i,j,bi,bj)  
           Tsf     = Tsrf(i,j,bi,bj)  
           Tice(1) = Tice1(i,j,bi,bj)  
           Tice(2) = Tice2(i,j,bi,bj)  
           qicen(1)= Qice1(i,j,bi,bj)  
           qicen(2)= Qice2(i,j,bi,bj)  
           CALL THSICE_THERM(  
      I          fswdown, oceHs, oceV2s, oceSs, oceTs,  
      U          compact, hIce, hSnow, Tsf, Tice, qicen,  
      O          qleft, ffresh, fsalt, Fbot,  
      O          flxAtm, evpAtm,  
      I          i,j, bi,bj, myThid)  
   
 C-- Diagnostic of Atmospheric Fluxes over sea-ice :  
           frwAtm = evpAtm - snow(i,j,bi,bj)*rhos/rhofw  
126  C note: Any flux of mass (here fresh water) that enter or leave the system  C note: Any flux of mass (here fresh water) that enter or leave the system
127  C       with a non zero energy HAS TO be counted: add snow precip.  C       with a non zero energy HAS TO be counted: add snow precip.
128            flxAtm = flxAtm - Lfresh*snow(i,j,bi,bj)*rhos            icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
129         &                        - Lfresh*snowPrc(i,j,bi,bj)
130    C--
131            ENDIF
132           ENDDO
133          ENDDO
134    
135    #ifdef ALLOW_DIAGNOSTICS
136          IF ( useDiagnostics ) THEN
137            tmpFac = 1. _d 0
138            CALL DIAGNOSTICS_FRACT_FILL(
139         I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
140         I                   0,1,1,bi,bj,myThid)
141            CALL DIAGNOSTICS_FRACT_FILL(
142         I                   siceAlb,   iceMask,tmpFac,1,'SIalbedo',
143         I                   0,1,1,bi,bj,myThid)
144          ENDIF
145    #endif /* ALLOW_DIAGNOSTICS */
146          DO j = jMin, jMax
147           DO i = iMin, iMax
148              siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
149           ENDDO
150          ENDDO
151    
152  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
153        IF (dBug) WRITE(6,1010) 'ThSI_FWD: iceFrac,flxAtm,evpAtm,flxSnw=',  C    part.2 : ice-covered fraction ;
154       &  iceMask(i,j,bi,bj),flxAtm,evpAtm,-Lfresh*snow(i,j,bi,bj)*rhos  C     change in ice/snow thickness and ice-fraction
155        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qleft,fsalt,ffresh=',  C     note: can only reduce the ice-fraction but not increase it.
156       &   compact,qleft,fsalt,ffresh  C-------
157  #ifdef CHECK_ENERGY_CONSERV        DO j = jMin, jMax
158            iceFrac = iceMask(i,j,bi,bj)         DO i = iMin, iMax
           CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,  
      I            iceFrac, compact, hIce, hSnow, qicen,  
      I            qleft, ffresh, fsalt, flxAtm, frwAtm,  
      I            myTime, myIter, myThid )  
 #endif /* CHECK_ENERGY_CONSERV */  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
159    
160  C--    Update Sea-Ice state :          tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
161  c         theta(i,j,1,bi,bj) = oceTs          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
162  c         iceMask(i,j,bi,bj)=compact          frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
163            iceheight(i,j,bi,bj) = hIce       &                 * cphm/ocean_deltaT
164            snowheight(i,j,bi,bj)= hSnow          iceFrac(i,j) = iceMask(i,j,bi,bj)
165            Tsrf(i,j,bi,bj) =Tsf          flx2oc(i,j)  = icFlxSW(i,j,bi,bj)
166            Tice1(i,j,bi,bj)=Tice(1)  C-------
167            Tice2(i,j,bi,bj)=Tice(2)  #ifdef ALLOW_DBUG_THSICE
168            Qice1(i,j,bi,bj)=qicen(1)          IF ( dBug(i,j,bi,bj) ) THEN
169            Qice2(i,j,bi,bj)=qicen(2)           IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
170              WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
171              WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
172         &                  iceFrac(i,j), iceHeight(i,j,bi,bj),
173         &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
174              WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
175         &                     tOceMxL(i,j,bi,bj), tFrzOce(i,j),
176         &                     frzmltMxL(i,j), Qnet(i,j,bi,bj)
177             ENDIF
178             IF (iceFrac(i,j).GT.0.)
179         &    WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
180         &      iceFrac(i,j), icFlxAtm(i,j,bi,bj),
181         &      icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
182            ENDIF
183    #endif
184           ENDDO
185          ENDDO
186    
187          CALL THSICE_CALC_THICKN(
188         I          bi, bj, siLo, siHi, sjLo, sjHi,
189         I          iMin,iMax, jMin,jMax, dBugFlag,
190         I          iceMask(siLo,sjLo,bi,bj), tFrzOce,
191         I          tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
192         I          snowPrc(siLo,sjLo,bi,bj), prcAtm,
193         I          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
194         U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
195         U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
196         U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
197         U          icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
198         O          frw2oc, fsalt,
199         I          myTime, myIter, myThid )
200    
201  C--    Net fluxes :  C--    Net fluxes :
202            ffresh = ffresh/rhofw        DO j = jMin, jMax
203            ffresh = -ffresh-rain(i,j,bi,bj)-runoff(i,j,bi,bj)         DO i = iMin, iMax
204            frwAtm =  frwAtm-rain(i,j,bi,bj)-runoff(i,j,bi,bj)          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
205            iceFrac = iceMask(i,j,bi,bj)  C-     weighted average net fluxes:
206            openFrac= 1. _d 0-iceFrac            icFrac = iceMask(i,j,bi,bj)
207  #ifdef ALLOW_TIMEAVE            opFrac= 1. _d 0-icFrac
208            ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
209       &          + ( -iceFrac*flxAtm + openFrac*Qnet(i,j,bi,bj)       &                        - opFrac*Qnet(i,j,bi,bj)
210       &            )*thSIce_deltaT            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
211            ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)       &                        + opFrac*rhofw*EmPmR(i,j,bi,bj)
212       &          + ( iceFrac*frwAtm + openFrac*EmPmR(i,j,bi,bj)            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
213       &            )*thSIce_deltaT            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
214  #endif /*ALLOW_TIMEAVE*/       &                    +  opFrac*EmPmR(i,j,bi,bj)
215            Qnet(i,j,bi,bj)=-iceFrac*qleft + openFrac*Qnet(i,j,bi,bj)            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
216            EmPmR(i,j,bi,bj)=iceFrac*ffresh+openFrac*EmPmR(i,j,bi,bj)  
217            saltFlux(i,j,bi,bj)=-iceFrac*fsalt  #ifdef ALLOW_DBUG_THSICE
218              IF (dBug(i,j,bi,bj)) WRITE(6,1010)
219            IF (dBug) WRITE(6,1010)'ThSI_FWD:-1- compact,hIc,hSn,Qnet=',       &          'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
220       &                      compact,hIce,hSnow,Qnet(i,j,bi,bj)       &           iceFrac(i,j), iceHeight(i,j,bi,bj),
221         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
222          ELSEIF (hFacC(i,j,1,bi,bj).gt.0. _d 0) THEN  #endif
   
 #ifdef ALLOW_TIMEAVE  
          ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)  
      &                   +Qnet(i,j,bi,bj)*thSIce_deltaT  
          ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)  
      &                   +EmPmR(i,j,bi,bj)*thSIce_deltaT  
 #endif /*ALLOW_TIMEAVE*/  
223    
224            ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
225              icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
226              icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
227            ELSE
228              icFlxAtm(i,j,bi,bj) = 0. _d 0
229              icFrwAtm(i,j,bi,bj) = 0. _d 0
230          ENDIF          ENDIF
231           ENDDO
232          ENDDO
233    
234  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
235  C    part.2 : freezing of sea-water  C    part.3 : freezing of sea-water
236  C     over ice-free fraction and what is left from ice-covered fraction  C     over ice-free fraction and what is left from ice-covered fraction
237  C-------  C-------
238          esurp = frzmlt - Fbot*iceMask(i,j,bi,bj)        CALL THSICE_EXTEND(
239          IF (esurp.GT.0. _d 0) THEN       I          bi, bj, siLo, siHi, sjLo, sjHi,
240            iceFrac = compact       I          iMin,iMax, jMin,jMax, dBugFlag,
241            IF ( compact.GT.0. _d 0 ) THEN       I          frzmltMxL, tFrzOce,
242              qicen(1)= Qice1(i,j,bi,bj)       I          tOceMxL(siLo,sjLo,bi,bj),
243              qicen(2)= Qice2(i,j,bi,bj)       U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
244            ELSE       U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
245              qicen(1)= -cpwater*Tmlt1       U          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
246       &               + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)       U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
247              qicen(2)= -cpice *Tf + Lfresh       O          flx2oc, frw2oc, fsalt,
248            ENDIF       I          myTime, myIter, myThid )
249            qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0  
250            oceTs = theta(i,j,1,bi,bj)        DO j = jMin, jMax
251            hIce  = iceHeight(i,j,bi,bj)         DO i = iMin, iMax
252            hSnow = snowHeight(i,j,bi,bj)          IF (frzmltMxL(i,j).GT.0. _d 0) THEN
           CALL THSICE_START( myThid,  
      I             esurp, qicAv, Tf,  
      O             qNewIce, ffresh, fsalt,  
      U             oceTs, compact, hIce, hSnow )  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
       IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qNewIce,fsalt,ffresh='  
      &                 ,compact,qNewIce,fsalt,ffresh  
 #ifdef CHECK_ENERGY_CONSERV  
           flxAtm = 0.  
           frwAtm = 0.  
           CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,  
      I            iceFrac, compact, hIce, hSnow, qicen,  
      I            qNewIce, ffresh, fsalt, flxAtm, frwAtm,  
      I            myTime, myIter, myThid )  
 #endif /* CHECK_ENERGY_CONSERV */  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C--    Update Sea-Ice state :  
           IF ( compact.GT.0. _d 0 .AND. iceFrac.EQ.0. _d 0) THEN  
              Tsrf(i,j,bi,bj)  = Tf  
              Tice1(i,j,bi,bj) = Tf  
              Tice2(i,j,bi,bj) = Tf  
              Qice1(i,j,bi,bj) = qicen(1)  
              Qice2(i,j,bi,bj) = qicen(2)  
 c            theta(i,j,1,bi,bj)= oceTs  
           ENDIF  
           iceheight(i,j,bi,bj) = hIce  
           snowheight(i,j,bi,bj)= hSnow  
253  C--    Net fluxes :  C--    Net fluxes :
254            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - qNewIce            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
255            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- ffresh/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
256            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
257    
258            IF (dBug) WRITE(6,1010)'ThSI_FWD:-2- compact,hIc,hSn,Qnet=',  #ifdef ALLOW_DBUG_THSICE
259       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
260  C--   - if esurp > 0 : end       &         'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
261         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
262         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
263    #endif
264          ENDIF          ENDIF
265    
266          IF ( compact .GT. 0. _d 0 ) THEN          IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
267            iceMask(i,j,bi,bj)=compact       &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
268         &                  .AND.   iceFrac(i,j) .LE.0. _d 0
269            IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
270              iceMask(i,j,bi,bj)=iceFrac(i,j)
271              IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
272         &     snowAge(i,j,bi,bj) = 0. _d 0
273          ELSE          ELSE
274            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
275            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
276            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
277            Tsrf(i,j,bi,bj)=theta(i,j,1,bi,bj)            snowAge(i,j,bi,bj)  = 0. _d 0
278              Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
279            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
280            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
281            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = 0. _d 0
282            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = 0. _d 0
283          ENDIF          ENDIF
284    
285  #ifndef CHECK_ENERGY_CONSERV  #ifdef ATMOSPHERIC_LOADING
286  #ifdef ALLOW_TIMEAVE  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
287            ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
288       &         + ( Qnet(i,j,bi,bj)       &                         + iceHeight(i,j,bi,bj)*rhoi
289       &            )*thSIce_deltaT       &                        )*iceMask(i,j,bi,bj)
290            ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)  #endif
      &         + ( EmPmR(i,j,bi,bj)  
      &            )*thSIce_deltaT  
           ICE_salFx_AVE(i,j,bi,bj)=ICE_salFx_AVE(i,j,bi,bj)  
      &            +saltFlux(i,j,bi,bj)*thSIce_deltaT  
 #endif /* ALLOW_TIMEAVE */  
 #endif /* CHECK_ENERGY_CONSERV */  
291    
292         ENDDO         ENDDO
293        ENDDO        ENDDO
294    
295          IF ( thSIceAdvScheme.GT.0 ) THEN
296    C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
297           DO j = jMin, jMax
298            DO i = iMin, iMax
299             IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
300              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
301              EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
302              saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
303             ENDIF
304            ENDDO
305           ENDDO
306          ENDIF
307    
308    #ifdef ALLOW_BULK_FORCE
309          IF ( useBulkForce ) THEN
310            CALL BULKF_FLUX_ADJUST(
311         I                          bi, bj, iMin, iMax, jMin, jMax,
312         I                          isIceFree, myTime, myIter, myThid )
313          ENDIF
314    #endif /* ALLOW_BULK_FORCE */
315    
316  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
317  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
318    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22