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