/[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.10 by jmc, Sun Jul 10 01:24:25 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     Fbot          :: the oceanic heat flux already incorporated (ice_therm)  C     snowPr    :: snow precipitation [kg/m2/s]
65  C     flxAtm        :: net heat flux from the atmosphere ( >0 downward)  C     agingTime :: aging time scale (s)
66  C     evpAtm        :: evaporation to the atmosphere  C     ageFac    :: snow aging factor [1]
67  C     frwAtm        :: net fresh-water flux (E-P-R) to the atmosphere (m/s)  C     albedo    :: surface albedo [0-1]
68  C     qleft         :: net heat flux from the ice to the ocean  C     flxAtm    :: net heat flux from the atmosphere (+=down) [W/m2]
69  C     ffresh        :: fresh-water flux from the ice to the ocean  C     frwAtm    :: net fresh-water flux (E-P) to the atmosphere  [kg/m2/s]
70  C     fsalt         :: mass salt flux to the ocean  C     Fbot      :: the oceanic heat flux already incorporated (ice_therm)
71    C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
72    C     frw2oc    :: fresh-water flux from the ice to the ocean
73    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 fswdown, qleft, qNewIce        _RL snowPr
78          _RL agingTime, ageFac
79          _RL albedo
80          _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _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 tmpFac
93    #endif
94    
95        LOGICAL dBug        LOGICAL dBug
96    
       dBug = .FALSE.  
97   1010 FORMAT(A,1P4E11.3)   1010 FORMAT(A,1P4E11.3)
98          dBug = .FALSE.
99    C-    Initialise flxAtm
100           DO j = 1-Oly, sNy+Oly
101            DO i = 1-Olx, sNx+Olx
102              flxAtm(i,j) = 0.
103            ENDDO
104           ENDDO
105    
106          IF ( fluidIsWater ) THEN
107           DO j = jMin, jMax
108            DO i = iMin, iMax
109    c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
110    
111    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
112    C    part.1 : ice-covered fraction ;
113    C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes
114    C-------
115             IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
116              icFrac  = iceMask(i,j,bi,bj)
117              TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
118              hIce    = iceHeight(i,j,bi,bj)
119              hSnow   = snowHeight(i,j,bi,bj)
120              Tsf     = Tsrf(i,j,bi,bj)
121              qicen(1)= Qice1(i,j,bi,bj)
122              qicen(2)= Qice2(i,j,bi,bj)
123              IF ( dBug ) THEN
124               WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
125               WRITE(6,1010) 'ThSI_FWD:-0- iceMask, hIc, hSn, Tsf  =',
126         &                                 icFrac, hIce,hSnow,Tsf
127              ENDIF
128    
129              CALL THSICE_ALBEDO(
130         I               hIce, hSnow, Tsf, snowAge(i,j,bi,bj),
131         O               albedo,
132         I               myThid )
133              flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)
134    
135              CALL THSICE_SOLVE4TEMP(
136         I          useBulkforce, tmpflx, TFrzOce, hIce, hSnow,
137         U          flxSW(i,j), Tsf, qicen,
138         O          Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
139         O          tmpdTs, flxAtm(i,j), evpAtm(i,j),
140         I          i,j, bi,bj, myThid)
141    
142    #ifdef SHORTWAVE_HEATING
143    C--    Update Fluxes :
144              opFrac= 1. _d 0-icFrac
145              Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj)
146    #endif
147    C--    Update Sea-Ice state :
148              Tsrf(i,j,bi,bj) =Tsf
149              Tice1(i,j,bi,bj)=Tice(1)
150              Tice2(i,j,bi,bj)=Tice(2)
151              Qice1(i,j,bi,bj)=qicen(1)
152              Qice2(i,j,bi,bj)=qicen(2)
153              siceAlb(i,j,bi,bj) = icFrac*albedo
154              IF ( dBug ) THEN
155               WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =',
156         &                              Tsf, Tice, frzmltMxL
157               WRITE(6,1010) 'ThSI_FWD: sHeat,fxCndBt, fxAtm,evAtm=',
158         &                  sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
159         &                  flxAtm(i,j), evpAtm(i,j)
160              ENDIF
161             ENDIF
162            ENDDO
163           ENDDO
164          ENDIF
165          dBug = .FALSE.
166    
167    #ifdef ALLOW_DIAGNOSTICS
168          IF ( useDiagnostics ) THEN
169            tmpFac = 1. _d 0
170            CALL DIAGNOSTICS_FRACT_FILL(
171         I                   snowPrc,   iceMask,tmpFac,'SIsnwPrc',
172         I                   0,1,1,bi,bj,myThid)
173          ENDIF
174    #endif /* ALLOW_DIAGNOSTICS */
175    
176    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
177    C    part.2 : ice-covered fraction ;
178    C     change in ice/snow thickness and ice-fraction
179    C     note: can only reduce the ice-fraction but not increase it.
180    C-------
181          agingTime = 50. _d 0 * 86400. _d 0
182          ageFac = 1. _d 0 - thSIce_deltaT/agingTime
183        DO j = jMin, jMax        DO j = jMin, jMax
184         DO i = iMin, iMax         DO i = iMin, iMax
185  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 )
186    
187            TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
188            oceTs   = tOceMxL(i,j,bi,bj)
189            cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
190            frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT
191    
         Tf     = -mu_Tf*salt(i,j,1,bi,bj)  
         cphm   = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)  
         frzmlt = (Tf-theta(i,j,1,bi,bj))*cphm/thSIce_deltaT  
192          Fbot   = 0. _d 0          Fbot   = 0. _d 0
         compact= 0. _d 0  
         snow(i,j,bi,bj)     = 0. _d 0  
193          saltFlux(i,j,bi,bj) = 0. _d 0          saltFlux(i,j,bi,bj) = 0. _d 0
194            compact= iceMask(i,j,bi,bj)
195          IF (dBug.AND.(frzmlt.GT.0. .OR.iceMask(i,j,bi,bj).GT.0.)) THEN  C-------
196            WRITE(6,1010) 'ThSI_FWD:-0- iceMask,hIc,hSn,Qnet=',          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN
197       &       iceMask(i,j,bi,bj),iceHeight(i,j,bi,bj),            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
198       &       snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
199            WRITE(6,1010) 'ThSI_FWD: ocTs,Tf,frzmlt=',       &                  compact, iceHeight(i,j,bi,bj),
200       &              theta(i,j,1,bi,bj),Tf,frzmlt       &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
201              WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',
202         &                     oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)
203          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.  
204  C-------  C-------
205          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
206            fswdown = solar(i,j,bi,bj)  
207            oceHs   = hfacC(i,j,1,bi,bj)*drF(1)            oceV2s  = v2ocMxL(i,j,bi,bj)
208            oceTs   = theta(i,j,1,bi,bj)            snowPr  = snowPrc(i,j,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)  
209            hIce    = iceHeight(i,j,bi,bj)            hIce    = iceHeight(i,j,bi,bj)
210            hSnow   = snowHeight(i,j,bi,bj)            hSnow   = snowHeight(i,j,bi,bj)
211            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)  
212            qicen(1)= Qice1(i,j,bi,bj)            qicen(1)= Qice1(i,j,bi,bj)
213            qicen(2)= Qice2(i,j,bi,bj)            qicen(2)= Qice2(i,j,bi,bj)
214            CALL THSICE_THERM(            flx2oc  = flxSW(i,j)
215       I          fswdown, oceHs, oceV2s, oceSs, oceTs,  
216       U          compact, hIce, hSnow, Tsf, Tice, qicen,            CALL THSICE_CALC_THICKN(
217       O          qleft, ffresh, fsalt, Fbot,       I          frzmltMxL, TFrzOce, oceTs, oceV2s, snowPr,
218       O          flxAtm, evpAtm,       I          sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), evpAtm(i,j),
219       I          i,j, bi,bj, myThid)       U          compact, hIce, hSnow, Tsf, qicen, flx2oc,
220         O          frw2oc, fsalt, Fbot,
221         I          dBug, myThid)
222    
223    C- note : snowPr was not supposed to be modified in THSICE_THERM ;
224    C         but to reproduce old results, is reset to zero if Tsf >= 0
225              snowPrc(i,j,bi,bj) = snowPr
226    
227    C--  Snow aging :
228              snowAge(i,j,bi,bj) = thSIce_deltaT
229         &                       + snowAge(i,j,bi,bj)*ageFac
230              IF ( snowPr.GT.0. _d 0 )
231         &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
232         &          * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )
233    C--
234    
235  C-- Diagnostic of Atmospheric Fluxes over sea-ice :  C-- Diagnostic of Atmospheric Fluxes over sea-ice :
236            frwAtm = evpAtm - snow(i,j,bi,bj)*rhos/rhofw            frwAtm = evpAtm(i,j) - prcAtm(i,j)
237  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
238  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.
239            flxAtm = flxAtm - Lfresh*snow(i,j,bi,bj)*rhos            flxAtm(i,j) = flxAtm(i,j) - Lfresh*snowPrc(i,j,bi,bj)
240    
241  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
242        IF (dBug) WRITE(6,1010) 'ThSI_FWD: iceFrac,flxAtm,evpAtm,flxSnw=',        IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
243       &  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
244        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qleft,fsalt,ffresh=',        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=',
245       &   compact,qleft,fsalt,ffresh       &   compact,flx2oc,fsalt,frw2oc
246  #ifdef CHECK_ENERGY_CONSERV  #ifdef CHECK_ENERGY_CONSERV
247            iceFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
248            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,
249       I            iceFrac, compact, hIce, hSnow, qicen,       I            icFrac, compact, hIce, hSnow, qicen,
250       I            qleft, ffresh, fsalt, flxAtm, frwAtm,       I            flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm,
251       I            myTime, myIter, myThid )       I            myTime, myIter, myThid )
252  #endif /* CHECK_ENERGY_CONSERV */  #endif /* CHECK_ENERGY_CONSERV */
253  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
254    
255  C--    Update Sea-Ice state :  C--    Update Sea-Ice state :
 c         theta(i,j,1,bi,bj) = oceTs  
256  c         iceMask(i,j,bi,bj)=compact  c         iceMask(i,j,bi,bj)=compact
257            iceheight(i,j,bi,bj) = hIce            iceHeight(i,j,bi,bj) = hIce
258            snowheight(i,j,bi,bj)= hSnow            snowHeight(i,j,bi,bj)= hSnow
259            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)  
260            Qice1(i,j,bi,bj)=qicen(1)            Qice1(i,j,bi,bj)=qicen(1)
261            Qice2(i,j,bi,bj)=qicen(2)            Qice2(i,j,bi,bj)=qicen(2)
262    
263  C--    Net fluxes :  C--    Net fluxes :
264            ffresh = ffresh/rhofw            frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj))
265            ffresh = -ffresh-rain(i,j,bi,bj)-runoff(i,j,bi,bj)  C-     weighted average net fluxes:
266            frwAtm =  frwAtm-rain(i,j,bi,bj)-runoff(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
267            iceFrac = iceMask(i,j,bi,bj)            opFrac= 1. _d 0-icFrac
268            openFrac= 1. _d 0-iceFrac            flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj)
269  #ifdef ALLOW_TIMEAVE            frwAtm =     icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj)
270            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)
271       &          + ( -iceFrac*flxAtm + openFrac*Qnet(i,j,bi,bj)            EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj)
272       &            )*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  
 #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*/  
273    
274              IF (dBug) WRITE(6,1010)
275         &          'ThSI_FWD:-3- compact, hIc, hSn, Qnet =',
276         &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)
277    
278            ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
279              flxAtm(i,j) =  -Qnet(i,j,bi,bj)
280              frwAtm = rhofw*EmPmR(i,j,bi,bj)
281            ELSE
282              flxAtm(i,j) = 0. _d 0
283              frwAtm      = 0. _d 0
284          ENDIF          ENDIF
285    
286  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
287  C    part.2 : freezing of sea-water  C    part.3 : freezing of sea-water
288  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
289  C-------  C-------
290          esurp = frzmlt - Fbot*iceMask(i,j,bi,bj)  c       compact= iceMask(i,j,bi,bj)
291            hIce   = iceHeight(i,j,bi,bj)
292            hSnow  = snowHeight(i,j,bi,bj)
293    
294            esurp  = frzmltMxL - Fbot*iceMask(i,j,bi,bj)
295          IF (esurp.GT.0. _d 0) THEN          IF (esurp.GT.0. _d 0) THEN
296            iceFrac = compact            icFrac = compact
297            IF ( compact.GT.0. _d 0 ) THEN            qicen(1)= Qice1(i,j,bi,bj)
298              qicen(1)= Qice1(i,j,bi,bj)            qicen(2)= Qice2(i,j,bi,bj)
299              qicen(2)= Qice2(i,j,bi,bj)            CALL THSICE_EXTEND(
300            ELSE       I               esurp, TFrzOce,
301              qicen(1)= -cpwater*Tmlt1       U               oceTs, compact, hIce, hSnow, qicen,
302       &               + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)       O               flx2oc, frw2oc, fsalt,
303              qicen(2)= -cpice *Tf + Lfresh       I               dBug, myThid )
           ENDIF  
           qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0  
           oceTs = theta(i,j,1,bi,bj)  
           hIce  = iceHeight(i,j,bi,bj)  
           hSnow = snowHeight(i,j,bi,bj)  
           CALL THSICE_START( myThid,  
      I             esurp, qicAv, Tf,  
      O             qNewIce, ffresh, fsalt,  
      U             oceTs, compact, hIce, hSnow )  
304  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
305        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qNewIce,fsalt,ffresh='        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='
306       &                 ,compact,qNewIce,fsalt,ffresh       &                 ,compact,flx2oc,fsalt,frw2oc
307  #ifdef CHECK_ENERGY_CONSERV  #ifdef CHECK_ENERGY_CONSERV
308            flxAtm = 0.            tmpflx(1) = 0.
309            frwAtm = 0.            tmpflx(2) = 0.
310            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,            CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,
311       I            iceFrac, compact, hIce, hSnow, qicen,       I            icFrac, compact, hIce, hSnow, qicen,
312       I            qNewIce, ffresh, fsalt, flxAtm, frwAtm,       I            flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),
313       I            myTime, myIter, myThid )       I            myTime, myIter, myThid )
314  #endif /* CHECK_ENERGY_CONSERV */  #endif /* CHECK_ENERGY_CONSERV */
315  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316  C--    Update Sea-Ice state :  C--    Update Sea-Ice state :
317            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
318               Tsrf(i,j,bi,bj)  = Tf               Tsrf(i,j,bi,bj)  = TFrzOce
319               Tice1(i,j,bi,bj) = Tf               Tice1(i,j,bi,bj) = TFrzOce
320               Tice2(i,j,bi,bj) = Tf               Tice2(i,j,bi,bj) = TFrzOce
321               Qice1(i,j,bi,bj) = qicen(1)               Qice1(i,j,bi,bj) = qicen(1)
322               Qice2(i,j,bi,bj) = qicen(2)               Qice2(i,j,bi,bj) = qicen(2)
 c            theta(i,j,1,bi,bj)= oceTs  
323            ENDIF            ENDIF
324            iceheight(i,j,bi,bj) = hIce            iceHeight(i,j,bi,bj) = hIce
325            snowheight(i,j,bi,bj)= hSnow            snowHeight(i,j,bi,bj)= hSnow
326  C--    Net fluxes :  C--    Net fluxes :
327            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - qNewIce            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc
328            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- ffresh/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw
329            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt
330    
331            IF (dBug) WRITE(6,1010)'ThSI_FWD:-2- compact,hIc,hSn,Qnet=',            IF (dBug) WRITE(6,1010)
332         &          'ThSI_FWD:-4- compact, hIc, hSn, Qnet =',
333       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)
334  C--   - if esurp > 0 : end  C--   - if esurp > 0 : end
335          ENDIF          ENDIF
336    
337          IF ( compact .GT. 0. _d 0 ) THEN          IF ( compact .GT. 0. _d 0 ) THEN
338            iceMask(i,j,bi,bj)=compact            iceMask(i,j,bi,bj)=compact
339              IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0
340          ELSE          ELSE
341            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
342            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
343            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
344            Tsrf(i,j,bi,bj)=theta(i,j,1,bi,bj)            snowAge(i,j,bi,bj)  = 0. _d 0
345              Tsrf(i,j,bi,bj)     = oceTs
346            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
347            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
348            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = 0. _d 0
349            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = 0. _d 0
350          ENDIF          ENDIF
351    
352  #ifndef CHECK_ENERGY_CONSERV  C--     Return atmospheric fluxes in evpAtm & flxSW (same sign and units):
353  #ifdef ALLOW_TIMEAVE          evpAtm(i,j) = frwAtm
354            ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)          flxSW (i,j) = flxAtm(i,j)
355       &         + ( Qnet(i,j,bi,bj)  
356       &            )*thSIce_deltaT  #ifdef ATMOSPHERIC_LOADING
357            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)
358       &         + ( EmPmR(i,j,bi,bj)          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
359       &            )*thSIce_deltaT       &                         + iceHeight(i,j,bi,bj)*rhoi
360            ICE_salFx_AVE(i,j,bi,bj)=ICE_salFx_AVE(i,j,bi,bj)       &                        )*iceMask(i,j,bi,bj)
361       &            +saltFlux(i,j,bi,bj)*thSIce_deltaT  #endif
 #endif /* ALLOW_TIMEAVE */  
 #endif /* CHECK_ENERGY_CONSERV */  
362    
363         ENDDO         ENDDO
364        ENDDO        ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22