/[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.2 by jmc, Wed Dec 31 17:44:32 2003 UTC revision 1.9 by jmc, Mon Jan 31 19:37:06 2005 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         U             evpAtm, flxSW,
13       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
14    C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | SUBROUTINE  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
20    
21  C     !USES:  C     !USES:
22        IMPLICIT NONE        IMPLICIT NONE
23    
24  C     === Global variables ===  C     === Global variables ===
25  #include "SIZE.h"  #include "SIZE.h"
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "FFIELDS.h"  #include "FFIELDS.h"
 #include "DYNVARS.h"  
 #include "GRID.h"  
 #ifdef ALLOW_BULK_FORCE  
 #include "BULKF.h"  
 #endif  
29  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
30  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
31  #include "THSICE.h"  #include "THSICE_VARS.h"
32  #include "THSICE_DIAGS.h"  #include "THSICE_TAVE.h"
33    
34  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
35  C     === Routine arguments ===  C     === Routine arguments ===
36  C     myIter :: iteration counter for this thread  C     bi,bj   :: tile indices
37  C     myTime :: time counter for this thread  C   iMin,iMax :: computation domain: 1rst index range
38  C     myThid :: thread number for this instance of the routine.  C   jMin,jMax :: computation domain: 2nd  index range
39    C- input:
40    C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
41    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
43    C                      (part.1), becomes net SW flux into ocean (part.2).
44    C- output
45    C     evpAtm  :: (Out) net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
46    C     flxSW   :: (Out) net surf. heat flux from the atmosphere [W/m2], (+=down)
47    C     myTime  :: time counter for this thread
48    C     myIter  :: iteration counter for this thread
49    C     myThid  :: thread number for this instance of the routine.
50        INTEGER bi,bj        INTEGER bi,bj
51        INTEGER iMin, iMax        INTEGER iMin, iMax
52        INTEGER jMin, jMax        INTEGER jMin, jMax
53          _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54          _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55          _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56        _RL  myTime        _RL  myTime
57        INTEGER myIter        INTEGER myIter
58        INTEGER myThid        INTEGER myThid
59    CEOP
60    
61  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
62  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
63  C     === Local variables ===  C     === Local variables ===
64  C     snowAge       :: snow age (s)  C     snowPr    :: snow precipitation [kg/m2/s]
65  C     albedo        :: surface albedo [0-1]  C     agingTime :: aging time scale (s)
66  C     fSWabs        :: net Short-Wave (+ = down) at surface (W m-2)  C     ageFac    :: snow aging factor [1]
67  C     Fbot          :: the oceanic heat flux already incorporated (ice_therm)  C     albedo    :: surface albedo [0-1]
68  C     flxAtm        :: net heat flux from the atmosphere ( >0 downward)  C     flxAtm    :: net heat flux from the atmosphere (+=down) [W/m2]
69  C     evpAtm        :: evaporation to the atmosphere  C     frwAtm    :: net fresh-water flux (E-P) to the atmosphere  [kg/m2/s]
70  C     frwAtm        :: net fresh-water flux (E-P-R) to the atmosphere (m/s)  C     Fbot      :: the oceanic heat flux already incorporated (ice_therm)
71  C     qleft         :: net heat flux from the ice to the ocean  C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
72  C     ffresh        :: fresh-water flux from the ice to the ocean  C     frw2oc    :: fresh-water flux from the ice to the ocean
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]
75    C     TFrzOce   :: sea-water freezing temperature [oC] (function of S)
76        INTEGER i,j        INTEGER i,j
77        _RL snowAge        _RL snowPr
78          _RL agingTime, ageFac
79        _RL albedo        _RL albedo
80        _RL fSWabs        _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL qleft, qNewIce        _RL frwAtm
82          _RL flx2oc
83          _RL frw2oc
84        _RL fsalt        _RL fsalt
85        _RL ffresh        _RL TFrzOce, cphm, frzmltMxL
       _RL Tf, cphm, frzmlt  
86        _RL Fbot, esurp        _RL Fbot, esurp
87        _RL flxAtm, evpAtm, frwAtm        _RL opFrac, icFrac
88        _RL openFrac, iceFrac, qicAv        _RL oceV2s, oceTs
       _RL oceHs, oceV2s, oceSs, oceTs  
89        _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)        _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)
90          _RL tmpflx(0:2), tmpdTs
91    #ifdef ALLOW_DIAGNOSTICS
92          _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          LOGICAL  DIAGNOSTICS_IS_ON
94          EXTERNAL DIAGNOSTICS_IS_ON      
95    #endif
96    
97        LOGICAL dBug        LOGICAL dBug
98    
       dBug = .FALSE.  
99   1010 FORMAT(A,1P4E11.3)   1010 FORMAT(A,1P4E11.3)
100          dBug = .FALSE.
101    C-    Initialise flxAtm
102           DO j = 1-Oly, sNy+Oly
103            DO i = 1-Olx, sNx+Olx
104              flxAtm(i,j) = 0.
105            ENDDO
106           ENDDO
107    
108          IF ( fluidIsWater ) THEN
109           DO j = jMin, jMax
110            DO i = iMin, iMax
111    c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
112    
113    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
114    C    part.1 : ice-covered fraction ;
115    C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes
116    C-------
117             IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
118              icFrac  = iceMask(i,j,bi,bj)
119              TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
120              hIce    = iceHeight(i,j,bi,bj)
121              hSnow   = snowHeight(i,j,bi,bj)
122              Tsf     = Tsrf(i,j,bi,bj)
123              qicen(1)= Qice1(i,j,bi,bj)
124              qicen(2)= Qice2(i,j,bi,bj)
125              IF ( dBug ) THEN
126               WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
127               WRITE(6,1010) 'ThSI_FWD:-0- iceMask, hIc, hSn, Tsf  =',
128         &                                 icFrac, hIce,hSnow,Tsf
129              ENDIF
130    
131              CALL THSICE_ALBEDO(
132         I               hIce, hSnow, Tsf, snowAge(i,j,bi,bj),
133         O               albedo,
134         I               myThid )
135              flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)
136    
137              CALL THSICE_SOLVE4TEMP(
138         I          useBulkforce, tmpflx, TFrzOce, hIce, hSnow,
139         U          flxSW(i,j), Tsf, qicen,
140         O          Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
141         O          tmpdTs, flxAtm(i,j), evpAtm(i,j),
142         I          i,j, bi,bj, myThid)
143    
144    #ifdef SHORTWAVE_HEATING
145    C--    Update Fluxes :
146              opFrac= 1. _d 0-icFrac
147              Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj)
148    #endif
149    C--    Update Sea-Ice state :
150              Tsrf(i,j,bi,bj) =Tsf
151              Tice1(i,j,bi,bj)=Tice(1)
152              Tice2(i,j,bi,bj)=Tice(2)
153              Qice1(i,j,bi,bj)=qicen(1)
154              Qice2(i,j,bi,bj)=qicen(2)
155              siceAlb(i,j,bi,bj) = icFrac*albedo
156              IF ( dBug ) THEN
157               WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =',
158         &                              Tsf, Tice, frzmltMxL
159               WRITE(6,1010) 'ThSI_FWD: sHeat,fxCndBt, fxAtm,evAtm=',
160         &                  sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
161         &                  flxAtm(i,j), evpAtm(i,j)
162              ENDIF
163             ENDIF
164            ENDDO
165           ENDDO
166          ENDIF
167          dBug = .FALSE.
168    
169    #ifdef ALLOW_DIAGNOSTICS
170          IF ( useDiagnostics ) THEN
171    
172           IF ( DIAGNOSTICS_IS_ON('SIsnwPrc',myThid) ) THEN
173            DO j=1,sNy
174             DO i=1,sNx
175              tmpFld(i,j) = iceMask(i,j,bi,bj)*snowPrc(i,j,bi,bj)
176             ENDDO
177            ENDDO
178            CALL DIAGNOSTICS_FILL(tmpFld,'SIsnwPrc',0,1,2,bi,bj,myThid)
179           ENDIF
180    
181          ENDIF
182    #endif /* ALLOW_DIAGNOSTICS */
183    
184    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
185    C    part.2 : ice-covered fraction ;
186    C     change in ice/snow thickness and ice-fraction
187    C     note: can only reduce the ice-fraction but not increase it.
188    C-------
189          agingTime = 50. _d 0 * 86400. _d 0
190          ageFac = 1. _d 0 - thSIce_deltaT/agingTime
191        DO j = jMin, jMax        DO j = jMin, jMax
192         DO i = iMin, iMax         DO i = iMin, iMax
193  c       dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 )  c       dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
194    
195            TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
196            oceTs   = tOceMxL(i,j,bi,bj)
197            cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
198            frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT
199    
         Tf     = -mu_Tf*salt(i,j,1,bi,bj)  
         cphm   = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)  
         oceTs  = theta(i,j,1,bi,bj)  
         frzmlt = (Tf-oceTs)*cphm/thSIce_deltaT  
         compact= iceMask(i,j,bi,bj)  
         hIce   = iceHeight(i,j,bi,bj)  
         hSnow  = snowHeight(i,j,bi,bj)  
200          Fbot   = 0. _d 0          Fbot   = 0. _d 0
         snow(i,j,bi,bj)     = 0. _d 0  
201          saltFlux(i,j,bi,bj) = 0. _d 0          saltFlux(i,j,bi,bj) = 0. _d 0
202            compact= iceMask(i,j,bi,bj)
203          IF (dBug .AND. (frzmlt.GT.0. .OR. compact.GT.0.) ) THEN  C-------
204            WRITE(6,1010) 'ThSI_FWD:-0- iceMask,hIc,hSn,Qnet=',          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN
205       &                  compact, hIce, hSnow, Qnet(i,j,bi,bj)            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
206            WRITE(6,1010) 'ThSI_FWD: ocTs,Tf,frzmlt=',            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
207       &                            oceTs,Tf,frzmlt       &                  compact, iceHeight(i,j,bi,bj),
208         &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
209              WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',
210         &                     oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)
211          ENDIF          ENDIF
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C    part.1 : ice-covered fraction ;  
 C     can only reduce the ice-fraction but not increase it.  
212  C-------  C-------
213          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
214            oceHs   = hfacC(i,j,1,bi,bj)*drF(1)  
215            oceSs   = salt (i,j,1,bi,bj)            oceV2s  = v2ocMxL(i,j,bi,bj)
216            oceV2s  = ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)            snowPr  = snowPrc(i,j,bi,bj)
217       &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)            hIce    = iceHeight(i,j,bi,bj)
218       &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)            hSnow   = snowHeight(i,j,bi,bj)
      &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj) )*0.5 _d 0  
           snowAge = sage(i,j,bi,bj)  
 c         snowAge = thSIce_deltaT  
219            Tsf     = Tsrf(i,j,bi,bj)            Tsf     = Tsrf(i,j,bi,bj)
           Tice(1) = Tice1(i,j,bi,bj)  
           Tice(2) = Tice2(i,j,bi,bj)  
220            qicen(1)= Qice1(i,j,bi,bj)            qicen(1)= Qice1(i,j,bi,bj)
221            qicen(2)= Qice2(i,j,bi,bj)            qicen(2)= Qice2(i,j,bi,bj)
222            CALL THSICE_ALBEDO(hIce,hSnow,Tsf,snowAge,albedo)            flx2oc  = flxSW(i,j)
223            fSWabs = solar(i,j,bi,bj)*(1. _d 0 - albedo)  
224            CALL THSICE_THERM(            CALL THSICE_CALC_THICKN(
225       I          fSWabs, oceHs, oceV2s, oceSs, oceTs,       I          frzmltMxL, TFrzOce, oceTs, oceV2s, snowPr,
226       U          compact, hIce, hSnow, Tsf, Tice, qicen,       I          sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), evpAtm(i,j),
227       O          qleft, ffresh, fsalt, Fbot,       U          compact, hIce, hSnow, Tsf, qicen, flx2oc,
228       O          flxAtm, evpAtm,       O          frw2oc, fsalt, Fbot,
229       I          i,j, bi,bj, myThid)       I          dBug, myThid)
230    
231    C- note : snowPr was not supposed to be modified in THSICE_THERM ;
232    C         but to reproduce old results, is reset to zero if Tsf >= 0
233              snowPrc(i,j,bi,bj) = snowPr
234    
235    C--  Snow aging :
236              snowAge(i,j,bi,bj) = thSIce_deltaT
237         &                       + snowAge(i,j,bi,bj)*ageFac
238              IF ( snowPr.GT.0. _d 0 )
239         &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
240         &          * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )
241    C--
242    
243  C-- Diagnostic of Atmospheric Fluxes over sea-ice :  C-- Diagnostic of Atmospheric Fluxes over sea-ice :
244            frwAtm = evpAtm - snow(i,j,bi,bj)*rhos/rhofw            frwAtm = evpAtm(i,j) - prcAtm(i,j)
245  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
246  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.
247            flxAtm = flxAtm - Lfresh*snow(i,j,bi,bj)*rhos            flxAtm(i,j) = flxAtm(i,j) - Lfresh*snowPrc(i,j,bi,bj)
248    
249  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
250        IF (dBug) WRITE(6,1010) 'ThSI_FWD: iceFrac,flxAtm,evpAtm,flxSnw=',        IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
251       &  iceMask(i,j,bi,bj),flxAtm,evpAtm,-Lfresh*snow(i,j,bi,bj)*rhos       &  iceMask(i,j,bi,bj),flxAtm(i,j),evpAtm(i,j),-Lfresh*snowPr
252        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qleft,fsalt,ffresh=',        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=',
253       &   compact,qleft,fsalt,ffresh       &   compact,flx2oc,fsalt,frw2oc
254  #ifdef CHECK_ENERGY_CONSERV  #ifdef CHECK_ENERGY_CONSERV
255            iceFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
256            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,
257       I            iceFrac, compact, hIce, hSnow, qicen,       I            icFrac, compact, hIce, hSnow, qicen,
258       I            qleft, ffresh, fsalt, flxAtm, frwAtm,       I            flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm,
259       I            myTime, myIter, myThid )       I            myTime, myIter, myThid )
260  #endif /* CHECK_ENERGY_CONSERV */  #endif /* CHECK_ENERGY_CONSERV */
261  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
262    
263  C--    Update Sea-Ice state :  C--    Update Sea-Ice state :
264  c         iceMask(i,j,bi,bj)=compact  c         iceMask(i,j,bi,bj)=compact
265            iceheight(i,j,bi,bj) = hIce            iceHeight(i,j,bi,bj) = hIce
266            snowheight(i,j,bi,bj)= hSnow            snowHeight(i,j,bi,bj)= hSnow
267            Tsrf(i,j,bi,bj) =Tsf            Tsrf(i,j,bi,bj) =Tsf
           Tice1(i,j,bi,bj)=Tice(1)  
           Tice2(i,j,bi,bj)=Tice(2)  
268            Qice1(i,j,bi,bj)=qicen(1)            Qice1(i,j,bi,bj)=qicen(1)
269            Qice2(i,j,bi,bj)=qicen(2)            Qice2(i,j,bi,bj)=qicen(2)
270    
271  C--    Net fluxes :  C--    Net fluxes :
272            ffresh = ffresh/rhofw            frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj))
273            ffresh = -ffresh-rain(i,j,bi,bj)-runoff(i,j,bi,bj)  C-     weighted average net fluxes:
274            frwAtm =  frwAtm-rain(i,j,bi,bj)-runoff(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
275            iceFrac = iceMask(i,j,bi,bj)            opFrac= 1. _d 0-icFrac
276            openFrac= 1. _d 0-iceFrac            flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj)
277  #ifdef ALLOW_TIMEAVE            frwAtm =     icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj)
278            ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)            Qnet(i,j,bi,bj)=-icFrac*flx2oc +opFrac*Qnet(i,j,bi,bj)
279       &          + ( -iceFrac*flxAtm + openFrac*Qnet(i,j,bi,bj)            EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj)
280       &            )*thSIce_deltaT            saltFlux(i,j,bi,bj)=-icFrac*fsalt
           ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)  
      &          + ( iceFrac*frwAtm + openFrac*EmPmR(i,j,bi,bj)  
      &            )*thSIce_deltaT  
           ICE_albedo_AVE(i,j,bi,bj) = ICE_albedo_AVE(i,j,bi,bj)  
      &          + iceFrac*albedo*thSIce_deltaT  
 #endif /*ALLOW_TIMEAVE*/  
           Qnet(i,j,bi,bj)=-iceFrac*qleft + openFrac*Qnet(i,j,bi,bj)  
           EmPmR(i,j,bi,bj)=iceFrac*ffresh+openFrac*EmPmR(i,j,bi,bj)  
           saltFlux(i,j,bi,bj)=-iceFrac*fsalt  
   
           IF (dBug) WRITE(6,1010)'ThSI_FWD:-1- compact,hIc,hSn,Qnet=',  
      &                      compact,hIce,hSnow,Qnet(i,j,bi,bj)  
   
         ELSEIF (hFacC(i,j,1,bi,bj).gt.0. _d 0) THEN  
   
 #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*/  
281    
282              IF (dBug) WRITE(6,1010)
283         &          'ThSI_FWD:-3- compact, hIc, hSn, Qnet =',
284         &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)
285    
286            ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
287              flxAtm(i,j) =  -Qnet(i,j,bi,bj)
288              frwAtm = rhofw*EmPmR(i,j,bi,bj)
289            ELSE
290              flxAtm(i,j) = 0. _d 0
291              frwAtm      = 0. _d 0
292          ENDIF          ENDIF
293    
294  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
295  C    part.2 : freezing of sea-water  C    part.3 : freezing of sea-water
296  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
297  C-------  C-------
298          esurp = frzmlt - Fbot*iceMask(i,j,bi,bj)  c       compact= iceMask(i,j,bi,bj)
299            hIce   = iceHeight(i,j,bi,bj)
300            hSnow  = snowHeight(i,j,bi,bj)
301    
302            esurp  = frzmltMxL - Fbot*iceMask(i,j,bi,bj)
303          IF (esurp.GT.0. _d 0) THEN          IF (esurp.GT.0. _d 0) THEN
304            iceFrac = compact            icFrac = compact
305            IF ( compact.GT.0. _d 0 ) THEN            qicen(1)= Qice1(i,j,bi,bj)
306              qicen(1)= Qice1(i,j,bi,bj)            qicen(2)= Qice2(i,j,bi,bj)
307              qicen(2)= Qice2(i,j,bi,bj)            CALL THSICE_EXTEND(
308            ELSE       I               esurp, TFrzOce,
309              qicen(1)= -cpwater*Tmlt1       U               oceTs, compact, hIce, hSnow, qicen,
310       &               + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)       O               flx2oc, frw2oc, fsalt,
311              qicen(2)= -cpice *Tf + Lfresh       I               dBug, myThid )
           ENDIF  
           qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0  
           CALL THSICE_START( myThid,  
      I             esurp, qicAv, Tf,  
      O             qNewIce, ffresh, fsalt,  
      U             oceTs, compact, hIce, hSnow )  
312  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qNewIce,fsalt,ffresh='        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='
314       &                 ,compact,qNewIce,fsalt,ffresh       &                 ,compact,flx2oc,fsalt,frw2oc
315  #ifdef CHECK_ENERGY_CONSERV  #ifdef CHECK_ENERGY_CONSERV
316            flxAtm = 0.            tmpflx(1) = 0.
317            frwAtm = 0.            tmpflx(2) = 0.
318            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,
319       I            iceFrac, compact, hIce, hSnow, qicen,       I            icFrac, compact, hIce, hSnow, qicen,
320       I            qNewIce, ffresh, fsalt, flxAtm, frwAtm,       I            flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),
321       I            myTime, myIter, myThid )       I            myTime, myIter, myThid )
322  #endif /* CHECK_ENERGY_CONSERV */  #endif /* CHECK_ENERGY_CONSERV */
323  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
324  C--    Update Sea-Ice state :  C--    Update Sea-Ice state :
325            IF ( compact.GT.0. _d 0 .AND. iceFrac.EQ.0. _d 0) THEN            IF ( compact.GT.0. _d 0 .AND. icFrac.EQ.0. _d 0) THEN
326               Tsrf(i,j,bi,bj)  = Tf               Tsrf(i,j,bi,bj)  = TFrzOce
327               Tice1(i,j,bi,bj) = Tf               Tice1(i,j,bi,bj) = TFrzOce
328               Tice2(i,j,bi,bj) = Tf               Tice2(i,j,bi,bj) = TFrzOce
329               Qice1(i,j,bi,bj) = qicen(1)               Qice1(i,j,bi,bj) = qicen(1)
330               Qice2(i,j,bi,bj) = qicen(2)               Qice2(i,j,bi,bj) = qicen(2)
331            ENDIF            ENDIF
332            iceheight(i,j,bi,bj) = hIce            iceHeight(i,j,bi,bj) = hIce
333            snowheight(i,j,bi,bj)= hSnow            snowHeight(i,j,bi,bj)= hSnow
334  C--    Net fluxes :  C--    Net fluxes :
335            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - qNewIce            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc
336            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- ffresh/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw
337            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt
338    
339            IF (dBug) WRITE(6,1010)'ThSI_FWD:-2- compact,hIc,hSn,Qnet=',            IF (dBug) WRITE(6,1010)
340         &          'ThSI_FWD:-4- compact, hIc, hSn, Qnet =',
341       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)
342  C--   - if esurp > 0 : end  C--   - if esurp > 0 : end
343          ENDIF          ENDIF
344    
345          IF ( compact .GT. 0. _d 0 ) THEN          IF ( compact .GT. 0. _d 0 ) THEN
346            iceMask(i,j,bi,bj)=compact            iceMask(i,j,bi,bj)=compact
347            IF ( hSnow .EQ. 0. _d 0 ) sage(i,j,bi,bj) = 0. _d 0            IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0
348          ELSE          ELSE
349            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
350            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
351            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
352            sage(i,j,bi,bj)     = 0. _d 0            snowAge(i,j,bi,bj)  = 0. _d 0
353            Tsrf(i,j,bi,bj)     = oceTs            Tsrf(i,j,bi,bj)     = oceTs
354            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
355            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
# Line 254  C--   - if esurp > 0 : end Line 357  C--   - if esurp > 0 : end
357            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = 0. _d 0
358          ENDIF          ENDIF
359    
360  #ifndef CHECK_ENERGY_CONSERV  C--     Return atmospheric fluxes in evpAtm & flxSW (same sign and units):
361  #ifdef ALLOW_TIMEAVE          evpAtm(i,j) = frwAtm
362            ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)          flxSW (i,j) = flxAtm(i,j)
363       &         + ( Qnet(i,j,bi,bj)  
364       &            )*thSIce_deltaT  #ifdef ATMOSPHERIC_LOADING
365            ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
366       &         + ( EmPmR(i,j,bi,bj)          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
367       &            )*thSIce_deltaT       &                         + iceHeight(i,j,bi,bj)*rhoi
368            ICE_salFx_AVE(i,j,bi,bj)=ICE_salFx_AVE(i,j,bi,bj)       &                        )*iceMask(i,j,bi,bj)
369       &            +saltFlux(i,j,bi,bj)*thSIce_deltaT  #endif
 #endif /* ALLOW_TIMEAVE */  
 #endif /* CHECK_ENERGY_CONSERV */  
370    
371         ENDDO         ENDDO
372        ENDDO        ENDDO
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22