/[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.21 by heimbach, Tue Apr 17 23:42:33 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    #ifdef ALLOW_AUTODIFF_TAMC
34    # include "tamc.h"
35    # include "tamc_keys.h"
36    #endif
37    
38          INTEGER siLo, siHi, sjLo, sjHi
39          PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
40          PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
41    
42  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
43  C     === Routine arguments ===  C     === Routine arguments ===
44  C     myIter :: iteration counter for this thread  C- input:
45  C     myTime :: time counter for this thread  C     bi,bj   :: tile indices
46  C     myThid :: thread number for this instance of the routine.  C   iMin,iMax :: computation domain: 1rst index range
47    C   jMin,jMax :: computation domain: 2nd  index range
48    C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
49    C     myTime  :: current Time of simulation [s]
50    C     myIter  :: current Iteration number in simulation
51    C     myThid  :: my Thread Id number
52    C-- Use fluxes hold in commom blocks
53    C- input:
54    C     icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
55    C     icFlxAtm  :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
56    C     icFrwAtm  :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
57    C- output
58    C     icFlxAtm  :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
59    C     icFrwAtm  :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
60        INTEGER bi,bj        INTEGER bi,bj
61        INTEGER iMin, iMax        INTEGER iMin, iMax
62        INTEGER jMin, jMax        INTEGER jMin, jMax
63          _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64        _RL  myTime        _RL  myTime
65        INTEGER myIter        INTEGER myIter
66        INTEGER myThid        INTEGER myThid
67    CEOP
68    
69  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
70  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
71  C     === Local variables ===  C     === Local variables ===
72  C     snowAge       :: snow age (s)  C     iceFrac   :: fraction of grid area covered in ice
73  C     albedo        :: surface albedo [0-1]  C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
74  C     fSWabs        :: net Short-Wave (+ = down) at surface (W m-2)  C     frw2oc    :: fresh-water flux from the ice to the ocean
75  C     Fbot          :: the oceanic heat flux already incorporated (ice_therm)  C     fsalt     :: mass salt flux to the ocean
76  C     flxAtm        :: net heat flux from the atmosphere ( >0 downward)  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
77  C     evpAtm        :: evaporation to the atmosphere  C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
78  C     frwAtm        :: net fresh-water flux (E-P-R) to the atmosphere (m/s)  C     isIceFree :: true for ice-free grid-cell that remains ice-free
79  C     qleft         :: net heat flux from the ice to the ocean  C     ageFac    :: snow aging factor [1]
80  C     ffresh        :: fresh-water flux from the ice to the ocean  C     snowFac   :: snowing refreshing-age factor [units of 1/snowPr]
81  C     fsalt         :: mass salt flux to the ocean        LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL     iceFrac  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     flx2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL     frw2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL     fsalt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86          _RL     tFrzOce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL     frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL ageFac
89          _RL snowFac
90          _RL cphm
91          _RL opFrac, icFrac
92    #ifdef ALLOW_DIAGNOSTICS
93          _RL tmpFac
94    #endif
95        INTEGER i,j        INTEGER i,j
96        _RL snowAge        LOGICAL dBugFlag
       _RL albedo  
       _RL fSWabs  
       _RL qleft, qNewIce  
       _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)  
97    
98        LOGICAL dBug  C-    define grid-point location where to print debugging values
99    #include "THSICE_DEBUG.h"
100    
101        dBug = .FALSE.   1010 FORMAT(A,1P4E14.6)
  1010 FORMAT(A,1P4E11.3)  
102    
103    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104    
105    #ifdef ALLOW_AUTODIFF_TAMC
106          act1 = bi - myBxLo(myThid)
107          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
108          act2 = bj - myByLo(myThid)
109          max2 = myByHi(myThid) - myByLo(myThid) + 1
110          act3 = myThid - 1
111          max3 = nTx*nTy
112          act4 = ikey_dynamics - 1
113          iicekey = (act1 + 1) + act2*max1
114         &                     + act3*max1*max2
115         &                     + act4*max1*max2*max3
116    #endif /* ALLOW_AUTODIFF_TAMC */
117    
118    C-    Initialise
119          dBugFlag = debugLevel.GE.debLevB
120          DO j = 1-OLy, sNy+OLy
121            DO i = 1-OLx, sNx+OLx
122              isIceFree(i,j) = .FALSE.
123              saltFlux(i,j,bi,bj) = 0. _d 0
124    #ifdef ALLOW_AUTODIFF_TAMC
125              iceFrac(i,j) = 0.
126    #endif
127            ENDDO
128          ENDDO
129    
130          ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
131          snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
132    
133    #ifdef ALLOW_AUTODIFF_TAMC
134    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
135    #endif
136        DO j = jMin, jMax        DO j = jMin, jMax
137         DO i = iMin, iMax         DO i = iMin, iMax
138  c       dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 )          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
139    C--  Snow aging :
140          Tf     = -mu_Tf*salt(i,j,1,bi,bj)            snowAge(i,j,bi,bj) = thSIce_deltaT
141          cphm   = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)       &                       + snowAge(i,j,bi,bj)*ageFac
142          oceTs  = theta(i,j,1,bi,bj)            IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
143          frzmlt = (Tf-oceTs)*cphm/thSIce_deltaT       &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
144          compact= iceMask(i,j,bi,bj)       &          * EXP( - snowFac*snowPrc(i,j,bi,bj) )
145          hIce   = iceHeight(i,j,bi,bj)  c    &          * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
146          hSnow  = snowHeight(i,j,bi,bj)  c    &                  /hNewSnowAge )
147          Fbot   = 0. _d 0  C-------
148          snow(i,j,bi,bj)     = 0. _d 0  C note: Any flux of mass (here fresh water) that enter or leave the system
149          saltFlux(i,j,bi,bj) = 0. _d 0  C       with a non zero energy HAS TO be counted: add snow precip.
150              icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
151          IF (dBug .AND. (frzmlt.GT.0. .OR. compact.GT.0.) ) THEN       &                        - Lfresh*snowPrc(i,j,bi,bj)
152            WRITE(6,1010) 'ThSI_FWD:-0- iceMask,hIc,hSn,Qnet=',  C--
      &                  compact, hIce, hSnow, Qnet(i,j,bi,bj)  
           WRITE(6,1010) 'ThSI_FWD: ocTs,Tf,frzmlt=',  
      &                            oceTs,Tf,frzmlt  
153          ENDIF          ENDIF
154           ENDDO
155          ENDDO
156    
157    #ifdef ALLOW_DIAGNOSTICS
158          IF ( useDiagnostics ) THEN
159            tmpFac = 1. _d 0
160            CALL DIAGNOSTICS_FRACT_FILL(
161         I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
162         I                   0,1,1,bi,bj,myThid)
163            CALL DIAGNOSTICS_FRACT_FILL(
164         I                   siceAlb,   iceMask,tmpFac,1,'SIalbedo',
165         I                   0,1,1,bi,bj,myThid)
166          ENDIF
167    #endif /* ALLOW_DIAGNOSTICS */
168          DO j = jMin, jMax
169           DO i = iMin, iMax
170              siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
171           ENDDO
172          ENDDO
173    
174  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
175  C    part.1 : ice-covered fraction ;  C    part.2 : ice-covered fraction ;
176  C     can only reduce the ice-fraction but not increase it.  C     change in ice/snow thickness and ice-fraction
177    C     note: can only reduce the ice-fraction but not increase it.
178  C-------  C-------
179          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN        DO j = jMin, jMax
180            oceHs   = hfacC(i,j,1,bi,bj)*drF(1)         DO i = iMin, iMax
           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  
           snowAge = sage(i,j,bi,bj)  
 c         snowAge = thSIce_deltaT  
           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_ALBEDO(hIce,hSnow,Tsf,snowAge,albedo)  
           fSWabs = solar(i,j,bi,bj)*(1. _d 0 - albedo)  
           CALL THSICE_THERM(  
      I          fSWabs, 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)  
181    
182  C-- Diagnostic of Atmospheric Fluxes over sea-ice :          tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
183            frwAtm = evpAtm - snow(i,j,bi,bj)*rhos/rhofw          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
184  C note: Any flux of mass (here fresh water) that enter or leave the system          frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
185  C       with a non zero energy HAS TO be counted: add snow precip.       &                 * cphm/ocean_deltaT
186            flxAtm = flxAtm - Lfresh*snow(i,j,bi,bj)*rhos          iceFrac(i,j) = iceMask(i,j,bi,bj)
187            flx2oc(i,j)  = icFlxSW(i,j,bi,bj)
188    C-------
189    #ifdef ALLOW_DBUG_THSICE
190            IF ( dBug(i,j,bi,bj) ) THEN
191             IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
192              WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
193              WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
194         &                  iceFrac(i,j), iceHeight(i,j,bi,bj),
195         &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
196              WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
197         &                     tOceMxL(i,j,bi,bj), tFrzOce(i,j),
198         &                     frzmltMxL(i,j), Qnet(i,j,bi,bj)
199             ENDIF
200             IF (iceFrac(i,j).GT.0.)
201         &    WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
202         &      iceFrac(i,j), icFlxAtm(i,j,bi,bj),
203         &      icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
204            ENDIF
205    #endif
206           ENDDO
207          ENDDO
208    
209  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  #ifdef ALLOW_AUTODIFF_TAMC
210        IF (dBug) WRITE(6,1010) 'ThSI_FWD: iceFrac,flxAtm,evpAtm,flxSnw=',  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
211       &  iceMask(i,j,bi,bj),flxAtm,evpAtm,-Lfresh*snow(i,j,bi,bj)*rhos  #endif
       IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qleft,fsalt,ffresh=',  
      &   compact,qleft,fsalt,ffresh  
 #ifdef CHECK_ENERGY_CONSERV  
           iceFrac = iceMask(i,j,bi,bj)  
           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-|--+----|  
212    
213  C--    Update Sea-Ice state :        CALL THSICE_CALC_THICKN(
214  c         iceMask(i,j,bi,bj)=compact       I          bi, bj, siLo, siHi, sjLo, sjHi,
215            iceheight(i,j,bi,bj) = hIce       I          iMin,iMax, jMin,jMax, dBugFlag,
216            snowheight(i,j,bi,bj)= hSnow       I          iceMask(siLo,sjLo,bi,bj), tFrzOce,
217            Tsrf(i,j,bi,bj) =Tsf       I          tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
218            Tice1(i,j,bi,bj)=Tice(1)       I          snowPrc(siLo,sjLo,bi,bj), prcAtm,
219            Tice2(i,j,bi,bj)=Tice(2)       I          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
220            Qice1(i,j,bi,bj)=qicen(1)       U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
221            Qice2(i,j,bi,bj)=qicen(2)       U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
222         U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
223         U          icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
224         O          frw2oc, fsalt,
225         I          myTime, myIter, myThid )
226    
227  C--    Net fluxes :  C--    Net fluxes :
228            ffresh = ffresh/rhofw        DO j = jMin, jMax
229            ffresh = -ffresh-rain(i,j,bi,bj)-runoff(i,j,bi,bj)         DO i = iMin, iMax
230            frwAtm =  frwAtm-rain(i,j,bi,bj)-runoff(i,j,bi,bj)  #ifdef ALLOW_AUTODIFF_TAMC
231            iceFrac = iceMask(i,j,bi,bj)            ikey_1 = i
232            openFrac= 1. _d 0-iceFrac       &         + sNx*(j-1)
233  #ifdef ALLOW_TIMEAVE       &         + sNx*sNy*act1
234            ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)       &         + sNx*sNy*max1*act2
235       &          + ( -iceFrac*flxAtm + openFrac*Qnet(i,j,bi,bj)       &         + sNx*sNy*max1*max2*act3
236       &            )*thSIce_deltaT       &         + sNx*sNy*max1*max2*max3*act4
237            ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)  #endif /* ALLOW_AUTODIFF_TAMC */
238       &          + ( iceFrac*frwAtm + openFrac*EmPmR(i,j,bi,bj)  C--
239       &            )*thSIce_deltaT  #ifdef ALLOW_AUTODIFF_TAMC
240            ICE_albedo_AVE(i,j,bi,bj) = ICE_albedo_AVE(i,j,bi,bj)  CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
241       &          + iceFrac*albedo*thSIce_deltaT  #endif
242  #endif /*ALLOW_TIMEAVE*/          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
243            Qnet(i,j,bi,bj)=-iceFrac*qleft + openFrac*Qnet(i,j,bi,bj)  C-     weighted average net fluxes:
244            EmPmR(i,j,bi,bj)=iceFrac*ffresh+openFrac*EmPmR(i,j,bi,bj)  #ifdef ALLOW_AUTODIFF_TAMC
245            saltFlux(i,j,bi,bj)=-iceFrac*fsalt  CADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
246    CADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
247            IF (dBug) WRITE(6,1010)'ThSI_FWD:-1- compact,hIc,hSn,Qnet=',  CADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
248       &                      compact,hIce,hSnow,Qnet(i,j,bi,bj)  CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
249    #endif
250          ELSEIF (hFacC(i,j,1,bi,bj).gt.0. _d 0) THEN            icFrac = iceMask(i,j,bi,bj)
251              opFrac= 1. _d 0-icFrac
252  #ifdef ALLOW_TIMEAVE            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
253           ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)       &                        - opFrac*Qnet(i,j,bi,bj)
254       &                   +Qnet(i,j,bi,bj)*thSIce_deltaT            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
255           ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)       &                        + opFrac*rhofw*EmPmR(i,j,bi,bj)
256       &                   +EmPmR(i,j,bi,bj)*thSIce_deltaT            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
257  #endif /*ALLOW_TIMEAVE*/            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
258         &                    +  opFrac*EmPmR(i,j,bi,bj)
259              saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
260    
261    #ifdef ALLOW_DBUG_THSICE
262              IF (dBug(i,j,bi,bj)) WRITE(6,1010)
263         &          'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
264         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
265         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
266    #endif
267    
268            ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
269              icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
270              icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
271            ELSE
272              icFlxAtm(i,j,bi,bj) = 0. _d 0
273              icFrwAtm(i,j,bi,bj) = 0. _d 0
274          ENDIF          ENDIF
275           ENDDO
276          ENDDO
277    
278  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
279  C    part.2 : freezing of sea-water  C    part.3 : freezing of sea-water
280  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
281  C-------  C-------
282          esurp = frzmlt - Fbot*iceMask(i,j,bi,bj)        CALL THSICE_EXTEND(
283          IF (esurp.GT.0. _d 0) THEN       I          bi, bj, siLo, siHi, sjLo, sjHi,
284            iceFrac = compact       I          iMin,iMax, jMin,jMax, dBugFlag,
285            IF ( compact.GT.0. _d 0 ) THEN       I          frzmltMxL, tFrzOce,
286              qicen(1)= Qice1(i,j,bi,bj)       I          tOceMxL(siLo,sjLo,bi,bj),
287              qicen(2)= Qice2(i,j,bi,bj)       U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
288            ELSE       U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
289              qicen(1)= -cpwater*Tmlt1       U          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
290       &               + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)       U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
291              qicen(2)= -cpice *Tf + Lfresh       O          flx2oc, frw2oc, fsalt,
292            ENDIF       I          myTime, myIter, myThid )
293            qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0  
294            CALL THSICE_START( myThid,  #ifdef ALLOW_AUTODIFF_TAMC
295       I             esurp, qicAv, Tf,  CADJ STORE snowHeight(:,:,bi,bj) =
296       O             qNewIce, ffresh, fsalt,  CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
297       U             oceTs, compact, hIce, hSnow )  #endif
298  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|        DO j = jMin, jMax
299        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qNewIce,fsalt,ffresh='         DO i = iMin, iMax
300       &                 ,compact,qNewIce,fsalt,ffresh          IF (frzmltMxL(i,j).GT.0. _d 0) THEN
 #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)  
           ENDIF  
           iceheight(i,j,bi,bj) = hIce  
           snowheight(i,j,bi,bj)= hSnow  
301  C--    Net fluxes :  C--    Net fluxes :
302            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - qNewIce            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
303            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
304            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
305    
306            IF (dBug) WRITE(6,1010)'ThSI_FWD:-2- compact,hIc,hSn,Qnet=',  #ifdef ALLOW_DBUG_THSICE
307       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
308  C--   - if esurp > 0 : end       &         'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
309         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
310         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
311    #endif
312          ENDIF          ENDIF
313    
314          IF ( compact .GT. 0. _d 0 ) THEN          IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
315            iceMask(i,j,bi,bj)=compact       &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
316            IF ( hSnow .EQ. 0. _d 0 ) sage(i,j,bi,bj) = 0. _d 0       &                  .AND.   iceFrac(i,j) .LE.0. _d 0
317            IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
318              iceMask(i,j,bi,bj)=iceFrac(i,j)
319              IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
320         &     snowAge(i,j,bi,bj) = 0. _d 0
321          ELSE          ELSE
322            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
323            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
324            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
325            sage(i,j,bi,bj)     = 0. _d 0            snowAge(i,j,bi,bj)  = 0. _d 0
326            Tsrf(i,j,bi,bj)     = oceTs            Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
327            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
328            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
329            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = Lfresh
330            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = Lfresh
331          ENDIF          ENDIF
332           ENDDO
333          ENDDO
334    
335  #ifndef CHECK_ENERGY_CONSERV  #ifdef ATMOSPHERIC_LOADING
336  #ifdef ALLOW_TIMEAVE  # ifdef ALLOW_AUTODIFF_TAMC
337            ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)  CADJ STORE snowHeight(:,:,bi,bj) =
338       &         + ( Qnet(i,j,bi,bj)  CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
339       &            )*thSIce_deltaT  # endif
340            ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)        DO j = jMin, jMax
341       &         + ( EmPmR(i,j,bi,bj)         DO i = iMin, iMax
342       &            )*thSIce_deltaT  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
343            ICE_salFx_AVE(i,j,bi,bj)=ICE_salFx_AVE(i,j,bi,bj)          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
344       &            +saltFlux(i,j,bi,bj)*thSIce_deltaT       &                         + iceHeight(i,j,bi,bj)*rhoi
345  #endif /* ALLOW_TIMEAVE */       &                        )*iceMask(i,j,bi,bj)
 #endif /* CHECK_ENERGY_CONSERV */  
   
346         ENDDO         ENDDO
347        ENDDO        ENDDO
348    #endif
349    
350          IF ( thSIceAdvScheme.GT.0 ) THEN
351    C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
352           DO j = jMin, jMax
353            DO i = iMin, iMax
354             IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
355              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
356              EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
357              saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
358             ENDIF
359            ENDDO
360           ENDDO
361          ENDIF
362    
363    #ifdef ALLOW_BULK_FORCE
364          IF ( useBulkForce ) THEN
365            CALL BULKF_FLUX_ADJUST(
366         I                          bi, bj, iMin, iMax, jMin, jMax,
367         I                          isIceFree, myTime, myIter, myThid )
368          ENDIF
369    #endif /* ALLOW_BULK_FORCE */
370    
371  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
372  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22