/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_growth.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by heimbach, Fri Dec 15 15:48:23 2006 UTC revision 1.10 by mlosch, Tue Jan 9 13:33:49 2007 UTC
# Line 40  C     i,j,bi,bj - Loop counters Line 40  C     i,j,bi,bj - Loop counters
40        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
41  C     number of surface interface layer  C     number of surface interface layer
42        INTEGER kSurface        INTEGER kSurface
43        _RL TBC, salinity_ice, SDF, ICE_DENS, Q0, QS  C     constants
44          _RL TBC, salinity_ice, SDF, ICE2WATR, ICE2SNOW
45          _RL QI, recip_QI, QS, recip_QS
46    C     auxillary variables
47          _RL availHeat, hEffOld, snowEnergy, snowAsIce
48          _RL growthHEFF, growthNeg
49  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
50        _RL hDraft, hFlood        _RL hDraft, hFlood
51  #endif /* ALLOW_SEAICE_FLOODING */  #endif /* ALLOW_SEAICE_FLOODING */
# Line 67  C     QSWI   - short wave heat flux unde Line 72  C     QSWI   - short wave heat flux unde
72        _RL QSWI   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL QSWI   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73  C      C    
74        _RL HCORR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL HCORR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75  C     SEAICE_SALT contains m of ice melted (<0) or created (>0)  C     frWtrIce contains m of ice melted (<0) or created (>0)
76        _RL SEAICE_SALT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL frWtrIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77  C     actual ice thickness  C     actual ice thickness with upper and lower limit
78        _RL HICE   (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL HICE   (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
79  C     actual snow thickness  C     actual snow thickness
80        _RL hSnwLoc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL hSnwLoc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
# Line 77  C     wind speed Line 82  C     wind speed
82        _RL UG     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL UG     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
83        _RL SPEED_SQ        _RL SPEED_SQ
84  C     local copy of AREA  C     local copy of AREA
85        _RL areaLoc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL areaLoc
86    
87  #ifdef SEAICE_MULTILEVEL  #ifdef SEAICE_MULTICATEGORY
88        INTEGER it        INTEGER it
89        INTEGER ilockey        INTEGER ilockey
90        _RL RK        _RL RK
91        _RL HICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL HICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
92        _RL FICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL FICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
93          _RL QSWIP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
94  #endif  #endif
95    
96        if ( buoyancyRelation .eq. 'OCEANICP' ) then        if ( buoyancyRelation .eq. 'OCEANICP' ) then
# Line 100  C     FREEZING TEMP. OF SEA WATER (deg C Line 106  C     FREEZING TEMP. OF SEA WATER (deg C
106  C     RATIO OF WATER DESITY TO SNOW DENSITY  C     RATIO OF WATER DESITY TO SNOW DENSITY
107        SDF          = 1000.0 _d 0/330.0 _d 0        SDF          = 1000.0 _d 0/330.0 _d 0
108  C     RATIO OF SEA ICE DESITY TO WATER DENSITY  C     RATIO OF SEA ICE DESITY TO WATER DENSITY
109        ICE_DENS     = 0.920 _d 0        ICE2WATR     = 0.920 _d 0
110  C     INVERSE HEAT OF FUSION OF ICE (m^3/J)  C     this makes more sense, but changes the results
111        Q0           = 1.0 _d -06 / 302.0 _d +00  C      ICE2WATR     = SEAICE_rhoIce * 1. _d -03
112    C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
113          ICE2SNOW     = SDF * ICE2WATR
114    C     HEAT OF FUSION OF ICE (m^3/J)
115          QI           = 302.0 _d +06
116          recip_QI     = 1.0 _d 0 / QI
117  C     HEAT OF FUSION OF SNOW (J/m^3)  C     HEAT OF FUSION OF SNOW (J/m^3)
118        QS           = 1.1 _d +08        QS           = 1.1 _d +08
119          recip_QS     = 1.1 _d 0 / QS
120    
121        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
122         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 134  CADJ &                         key = iic Line 146  CADJ &                         key = iic
146  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
147          DO J=1,sNy          DO J=1,sNy
148           DO I=1,sNx           DO I=1,sNx
149            areaLoc(I,J)     = MAX(A22,AREA(I,J,2,bi,bj))            FHEFF(I,J)      = 0.0 _d 0
150            FHEFF(I,J)       = 0.0 _d 0            FICE (I,J)      = 0.0 _d 0
151            FICE (I,J)       = 0.0 _d 0  #ifdef SEAICE_MULTICATEGORY
152  #ifdef SEAICE_MULTILEVEL            FICEP(I,J)      = 0.0 _d 0
153            FICEP(I,J)       = 0.0 _d 0            QSWIP(I,J)      = 0.0 _d 0
154  #endif  #endif
155            FHEFF(I,J)       = 0.0 _d 0            FHEFF(I,J)      = 0.0 _d 0
156            FICE (I,J)       = 0.0 _d 0            FICE (I,J)      = 0.0 _d 0
157            QNETO(I,J)       = 0.0 _d 0            QNETO(I,J)      = 0.0 _d 0
158            QNETI(I,J)       = 0.0 _d 0            QNETI(I,J)      = 0.0 _d 0
159            QSWO (I,J)       = 0.0 _d 0            QSWO (I,J)      = 0.0 _d 0
160            QSWI (I,J)       = 0.0 _d 0            QSWI (I,J)      = 0.0 _d 0
161            HCORR(I,J)       = 0.0 _d 0            HCORR(I,J)      = 0.0 _d 0
162            SEAICE_SALT(I,J) = 0.0 _d 0            frWtrIce(I,J)   = 0.0 _d 0
163            RESID_HEAT (I,J) = 0.0 _d 0            RESID_HEAT(I,J) = 0.0 _d 0
164           ENDDO           ENDDO
165          ENDDO          ENDDO
166  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 159  CADJ &                         key = iic Line 171  CADJ &                         key = iic
171  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
172          DO J=1,sNy          DO J=1,sNy
173           DO I=1,sNx           DO I=1,sNx
174  cph need to adjoint-store AREA again before using it in further init.  C     COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
175  cph (all these initialisations involving AREA are nasty "non-linear")  C     ON ICE THICKNESS FOR BUDGET COMPUTATION
176            HICE(I,J)    = HEFF(I,J,2,bi,bj)/areaLoc(I,J)  C     The default of A22 = 0.15 is a common threshold for defining
177            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc(I,J)  C     the ice edge. This ice concentration usually does not occur
178    C     due to thermodynamics but due to advection.
179              areaLoc      = MAX(A22,AREA(I,J,2,bi,bj))
180              HICE(I,J)    = HEFF(I,J,2,bi,bj)/areaLoc
181    C     Do we know what this is for?
182              HICE(I,J)    = MAX(HICE(I,J),0.05 _d +00)
183    C     Capping the actual ice thickness effectively enforces a
184    C     minimum of heat flux through the ice and helps getting rid of
185    C     very thick ice.
186              HICE(I,J)    = MIN(HICE(I,J),9.0 _d +00)
187              hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
188           ENDDO           ENDDO
189          ENDDO          ENDDO
190    
# Line 202  cphCADJ STORE uwind  = comlev1, key = ik Line 224  cphCADJ STORE uwind  = comlev1, key = ik
224  cphCADJ STORE vwind  = comlev1, key = ikey_dynamics  cphCADJ STORE vwind  = comlev1, key = ikey_dynamics
225  c  c
226  CADJ STORE tice   = comlev1, key = ikey_dynamics  CADJ STORE tice   = comlev1, key = ikey_dynamics
227  # ifdef SEAICE_MULTILEVEL  # ifdef SEAICE_MULTICATEGORY
228  CADJ STORE tices  = comlev1, key = ikey_dynamics  CADJ STORE tices  = comlev1, key = ikey_dynamics
229  # endif  # endif
230  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 216  C FIRST DO OPEN WATER Line 238  C FIRST DO OPEN WATER
238       I       bi, bj)       I       bi, bj)
239    
240  C NOW DO ICE  C NOW DO ICE
241  #ifdef SEAICE_MULTILEVEL  #ifdef SEAICE_MULTICATEGORY
242  C--  Start loop over muli-levels  C--  Start loop over muli-categories
243          DO IT=1,MULTDIM          DO IT=1,MULTDIM
244  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
245           ilockey = (iicekey-1)*MULTDIM + IT           ilockey = (iicekey-1)*MULTDIM + IT
246  CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,  CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
247  CADJ &                           key = ilockey, byte = isbyte  CADJ &                           key = ilockey, byte = isbyte
248  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
249             RK=REAL(IT)
250           DO J=1,sNy           DO J=1,sNy
251            DO I=1,sNx            DO I=1,sNx
252             RK=IT*1.0             HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)
            HICEP(I,J)=(HICE(I,J)/7.0 _d 0)*((2.0 _d 0*RK)-1.0 _d 0)  
253             TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)             TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
254            ENDDO            ENDDO
255           ENDDO           ENDDO
256           CALL SEAICE_BUDGET_ICE(           CALL SEAICE_BUDGET_ICE(
257       I        UG, HICE, hSnwLoc,       I        UG, HICEP, hSnwLoc,
258       U        TICE,       U        TICE,
259       O        FICE, QSWI,       O        FICEP, QSWIP,
260       I        bi, bj)       I        bi, bj)
261           DO J=1,sNy           DO J=1,sNy
262            DO I=1,sNx            DO I=1,sNx
263             FICEP(I,J)=(FICE(I,J)/7.0 _d 0)+FICEP(I,J)  C     average surface heat fluxes/growth rates
264             TICES(I,J,IT,bi,bj)=TICE(I,J,bi,bj)             FICE (I,J)          = FICE(I,J) + FICEP(I,J)/MULTDIM
265               QSWI (I,J)          = QSWI(I,J) + QSWIP(I,J)/MULTDIM
266               TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
267            ENDDO            ENDDO
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270  C--  End loop over muli-levels  C--  End loop over multi-categories
271          DO J=1,sNy  #else  /* SEAICE_MULTICATEGORY */
          DO I=1,sNx  
           FICE(I,J)=FICEP(I,J)  
          ENDDO  
         ENDDO  
 #else  /* SEAICE_MULTILEVEL */  
272          CALL SEAICE_BUDGET_ICE(          CALL SEAICE_BUDGET_ICE(
273       I       UG, HICE, hSnwLoc,       I       UG, HICE, hSnwLoc,
274       U       TICE,       U       TICE,
275       O       FICE, QSWI,       O       FICE, QSWI,
276       I       bi, bj)       I       bi, bj)
277  #endif /* SEAICE_MULTILEVEL */  #endif /* SEAICE_MULTICATEGORY */
278            
279  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
280  CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,  CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
# Line 263  CADJ &                         key = iic Line 282  CADJ &                         key = iic
282  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
283  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
284  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
285    C
286    C--   compute and apply ice growth due to oceanic heat flux from below
287    C
288          DO J=1,sNy          DO J=1,sNy
289           DO I=1,sNx           DO I=1,sNx
290  C--   Create or melt sea-ice so that first-level oceanic temperature  C--   Create or melt sea-ice so that first-level oceanic temperature
291  C     is approximately at the freezing point when there is sea-ice.  C     is approximately at the freezing point when there is sea-ice.
292  C     Initially the units of YNEG are m of sea-ice.  C     Initially the units of YNEG/availHeat are m of sea-ice.
293  C     The factor dRf(1)/72.0764, used to convert temperature  C     The factor dRf(1)/72.0764, used to convert temperature
294  C     change in deg K to m of sea-ice, is approximately:  C     change in deg K to m of sea-ice, is approximately:
295  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)
296  C        * (density of sea-water = 1026 kg/m^3)  C        * (density of sea-water = 1026 kg/m^3)
297  C        / (latent heat of fusion of sea-ice = 334000 J/kg)  C        / (latent heat of fusion of sea-ice = 334000 J/kg)
298  C        / (density of sea-ice = 910 kg/m^3)  C        / (density of sea-ice = 910 kg/m^3)
299  C     Negative YNEG leads to ice growth.  C     Negative YNEG/availHeat leads to ice growth.
300  C     Positive YNEG leads to ice melting.  C     Positive YNEG/availHeat leads to ice melting.
301            IF ( .NOT. inAdMode ) THEN            IF ( .NOT. inAdMode ) THEN
302  #ifdef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
303             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
304  #endif /* SEAICE_VARIABLE_FREEZING_POINT */  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
305             YNEG(I,J,bi,bj) = (theta(I,J,kSurface,bi,bj)-TBC)             availHeat = (theta(I,J,kSurface,bi,bj)-TBC)
306       &          *dRf(1)/72.0764 _d 0       &          *dRf(1)/72.0764 _d 0
307            ELSE            ELSE
308             YNEG(I,J,bi,bj)= 0.             availHeat = 0.
309            ENDIF            ENDIF
310            GHEFF(I,J)=HEFF(I,J,1,bi,bj)  C     local copy of old effective ice thickness
311  C     Melt (YNEG>0) or create (YNEG<0) sea ice            hEffOld           = HEFF(I,J,1,bi,bj)
312            HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj))  C     Melt (availHeat>0) or create (availHeat<0) sea ice
313            RESID_HEAT(I,J)  = YNEG(I,J,bi,bj)            HEFF(I,J,1,bi,bj) = MAX(ZERO,HEFF(I,J,1,bi,bj)-availHeat)
314            YNEG(I,J,bi,bj)  = GHEFF(I,J)-HEFF(I,J,1,bi,bj)  C    
315            SEAICE_SALT(I,J) = SEAICE_SALT(I,J)-YNEG(I,J,bi,bj)            YNEG(I,J,bi,bj)   = hEffOld       - HEFF(I,J,1,bi,bj)
316            RESID_HEAT(I,J)  =  RESID_HEAT(I,J)-YNEG(I,J,bi,bj)  C    
317              frWtrIce(I,J)     = frWtrIce(I,J) - YNEG(I,J,bi,bj)
318              RESID_HEAT(I,J)   = availHeat     - YNEG(I,J,bi,bj)
319  C     YNEG now contains m of ice melted (>0) or created (<0)  C     YNEG now contains m of ice melted (>0) or created (<0)
320  C     SEAICE_SALT contains m of ice melted (<0) or created (>0)  C     frWtrIce contains m of ice melted (<0) or created (>0)
321  C     RESID_HEAT is residual heat above freezing in equivalent m of ice  C     RESID_HEAT is residual heat above freezing in equivalent m of ice
322           ENDDO           ENDDO
323          ENDDO          ENDDO
# Line 314  CADJ STORE fice(:,:)         = comlev1_b Line 338  CADJ STORE fice(:,:)         = comlev1_b
338  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
339  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
340  cph)  cph)
341    C
342    C--   compute and apply ice growth due to atmospheric fluxes from above
343    C    
344          DO J=1,sNy          DO J=1,sNy
345           DO I=1,sNx           DO I=1,sNx
346  C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)  C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)
# Line 330  CADJ &                         key = iic Line 356  CADJ &                         key = iic
356          DO J=1,sNy          DO J=1,sNy
357           DO I=1,sNx           DO I=1,sNx
358            IF(FICE(I,J).LT.ZERO.AND.AREA(I,J,2,bi,bj).GT.ZERO) THEN            IF(FICE(I,J).LT.ZERO.AND.AREA(I,J,2,bi,bj).GT.ZERO) THEN
359  C use FICE to melt snow and CALCULATE CORRECTED GROWTH  C     use FICE to melt snow and CALCULATE CORRECTED GROWTH
360             GAREA(I,J)=HSNOW(I,J,bi,bj)*QS ! effective snow thickness in J/m^2  C     effective snow thickness in J/m^2
361             IF(GHEFF(I,J).LE.GAREA(I,J)) THEN             snowEnergy=HSNOW(I,J,bi,bj)*QS
362  C not enough heat to melt all snow; use up all heat flux FICE             IF(GHEFF(I,J).LE.snowEnergy) THEN
363    C     not enough heat to melt all snow; use up all heat flux FICE
364              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
365  C SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt  C     SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt
366  C The factor 1/SDF/ICE_DENS converts m of snow to m of sea-ice  C     The factor 1/ICE2SNOW converts m of snow to m of sea-ice
367              SEAICE_SALT(I,J)=SEAICE_SALT(I,J)              frWtrIce(I,J) = frWtrIce(I,J) - GHEFF(I,J)/(QS*ICE2SNOW)
368       &           -GHEFF(I,J)/QS/SDF/ICE_DENS              FICE    (I,J) = ZERO
             FICE(I,J)=ZERO  
369             ELSE             ELSE
370  C enought heat to melt snow completely;  C     enought heat to melt snow completely;
371  C compute remaining heat flux that will melt ice  C     compute remaining heat flux that will melt ice
372              FICE(I,J)=-(GHEFF(I,J)-GAREA(I,J))/              FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/
373       &           SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)       &           SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)
374  C     convert all snow to melt water (fresh water flux)  C     convert all snow to melt water (fresh water flux)
375              SEAICE_SALT(I,J)=SEAICE_SALT(I,J)              frWtrIce(I,J) = frWtrIce(I,J)
376       &           -HSNOW(I,J,bi,bj)/SDF/ICE_DENS       &           -HSNOW(I,J,bi,bj)/ICE2SNOW
377              HSNOW(I,J,bi,bj)=0.0              HSNOW(I,J,bi,bj)=0.0
378             END IF             END IF
379            END IF            END IF
# Line 361  CADJ &                         key = iic Line 387  CADJ &                         key = iic
387    
388          DO J=1,sNy          DO J=1,sNy
389           DO I=1,sNx           DO I=1,sNx
390  C NOW GET TOTAL GROWTH RATE in W/m^2, >0 causes ice growth  C     now get cell averaged growth rate in W/m^2, >0 causes ice growth
391            FHEFF(I,J)= FICE(I,J)  * AREA(I,J,2,bi,bj)            FHEFF(I,J)= FICE(I,J)  * AREA(I,J,2,bi,bj)
392       &              + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))       &              + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
393           ENDDO           ENDDO
# Line 384  CADJ STORE qswo(:,:)         = comlev1_b Line 410  CADJ STORE qswo(:,:)         = comlev1_b
410  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
411  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
412  cph)  cph)
         DO J=1,sNy  
          DO I=1,sNx  
 C NOW UPDATE AREA  
           GHEFF(I,J) = -SEAICE_deltaTtherm*FHEFF(I,J)*Q0  
           GAREA(I,J) =  SEAICE_deltaTtherm*QNETO(I,J)*Q0  
           GHEFF(I,J) = -ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))  
           GAREA(I,J) = MAX(ZERO,GAREA(I,J))  
           HCORR(I,J) = MIN(ZERO,GHEFF(I,J))  
          ENDDO  
         ENDDO  
413  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
414  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
415  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
416  #endif  #endif
417    C
418    C     First update (freeze or melt) ice area
419    C
420          DO J=1,sNy          DO J=1,sNy
421           DO I=1,sNx           DO I=1,sNx
422            GAREA(I,J)=(ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO  C     negative growth in meters of ice (>0 for melting)
423       &    +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
424       &    /(HEFF(I,J,1,bi,bj)+.00001 _d 0)  C     negative growth must not exceed effective ice thickness (=volume)
425            AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+GAREA(I,J)  C     (that is, cannot melt more than all the ice)
426              growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
427    C     growthHEFF < 0 means melting
428              HCORR(I,J) = MIN(ZERO,growthHEFF)
429    C     gain of new effective ice thickness over open water (>0 by definition)
430              GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)
431    CML   removed these loops and moved TAMC store directive up
432    CML         ENDDO
433    CML        ENDDO
434    CML#ifdef ALLOW_AUTODIFF_TAMC
435    CMLCADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
436    CMLCADJ &                         key = iicekey, byte = isbyte
437    CML#endif
438    CML        DO J=1,sNy
439    CML         DO I=1,sNx
440    C     Here we finally compute the new AREA
441              AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+
442         &         (ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
443         &         +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)
444         &         /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
445           ENDDO           ENDDO
446          ENDDO          ENDDO
447  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
448  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
449  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
450  #endif  #endif
451    C    
452    C     now update (freeze or melt) HEFF
453    C        
454          DO J=1,sNy          DO J=1,sNy
455           DO I=1,sNx           DO I=1,sNx
456    C     negative growth (>0 for melting) of existing ice in meters
457              growthNeg        = -SEAICE_deltaTtherm*
458         &                       FICE(I,J)*recip_QI*AREA(I,J,2,bi,bj)
459    C     negative growth must not exceed effective ice thickness (=volume)
460    C     (that is, cannot melt more than all the ice)
461              growthHEFF       = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
462    C     growthHEFF < 0 means melting
463              HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj) + growthHEFF
464    C     add effective growth to fresh water of ice
465              frWtrIce(I,J)    = frWtrIce(I,J)     + growthHEFF
466    
467    C     now calculate QNETI under ice (if any) as the difference between
468    C     the available "heat flux" growthNeg and the actual growthHEFF;
469    C     keep in mind that growthNeg and growthHEFF have different signs
470    C     by construction
471              QNETI(I,J) =  (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm
472    
473  C NOW UPDATE HEFF  C     now update other things
           GHEFF(I,J)       = -SEAICE_deltaTtherm*  
      &                       FICE(I,J)*Q0*AREA(I,J,2,bi,bj)  
           GHEFF(I,J)       = -ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))  
           HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj)+GHEFF(I,J)  
           SEAICE_SALT(I,J) = SEAICE_SALT(I,J)+GHEFF(I,J)  
   
 C NOW CALCULATE QNETI UNDER ICE IF ANY  
           QNETI(I,J)       = (GHEFF(I,J)-SEAICE_deltaTtherm*  
      &       FICE(I,J)*Q0*AREA(I,J,2,bi,bj))/Q0/SEAICE_deltaTtherm  
   
 C NOW UPDATE OTHER THINGS  
474    
475            IF(FICE(I,J).GT.ZERO) THEN            IF(FICE(I,J).GT.ZERO) THEN
476  C FREEZING, PRECIP ADDED AS SNOW  C     freezing, add precip as snow
477             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
478       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
479            ELSE            ELSE
480  C ADD PRECIP AS RAIN, WATER CONVERTED INTO equivalent m of ICE BY 1/ICE_DENS  C     add precip as rain, water converted into equivalent m of
481             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)  C     ice by 1/ICE2WATR.
482    C     Do not get confused by the sign:
483    C     precip   > 0 for downward flux of fresh water
484    C     frWtrIce > 0 for more ice (corresponds to an upward "fresh water flux"),
485    C     so that here the rain is added *as if* it is melted ice (which is not
486    C     true, but just a trick; physically the rain just runs as water
487    C     through the ice into the ocean)
488               frWtrIce(I,J) = frWtrIce(I,J)
489       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
490       &            SEAICE_deltaTtherm/ICE_DENS       &            SEAICE_deltaTtherm/ICE2WATR
491            ENDIF            ENDIF
492    
493  C Now add in precip over open water directly into ocean as negative salt  C     Now melt snow if there is residual heat left in surface level
494            SEAICE_SALT(I,J) = SEAICE_SALT(I,J)  C     Note that units of YNEG and frWtrIce are m of ice
      &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))  
      &         *SEAICE_deltaTtherm/ICE_DENS  
   
 C Now melt snow if there is residual heat left in surface level  
 C Note that units of YNEG and SEAICE_SALT are m of ice  
495  cph( very sensitive bit here by JZ  cph( very sensitive bit here by JZ
496            IF( RESID_HEAT(I,J) .GT. ZERO            IF( RESID_HEAT(I,J) .GT. ZERO .AND.
497       &        .AND. HSNOW(I,J,bi,bj) .GT. ZERO ) THEN       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
498             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE_DENS,             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,
499       &                             RESID_HEAT(I,J) )       &                             RESID_HEAT(I,J) )
500             YNEG(I,J,bi,bj)  = YNEG(I,J,bi,bj)+GHEFF(I,J)             YNEG(I,J,bi,bj)  = YNEG(I,J,bi,bj) +GHEFF(I,J)
501             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE_DENS             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR
502             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)-GHEFF(I,J)             frWtrIce(I,J)    = frWtrIce(I,J)   -GHEFF(I,J)
503            ENDIF            ENDIF
504  cph)  cph)
505    
506  C NOW GET FRESH WATER FLUX  C NOW GET FRESH WATER FLUX
507            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
508       &         EVAP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))       &         ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
509       &         -RUNOFF(I,J,bi,bj)       &         * ( ONE - AREA(I,J,2,bi,bj) )
510       &         +SEAICE_SALT(I,J)*ICE_DENS/SEAICE_deltaTtherm       &         - RUNOFF(I,J,bi,bj)
511         &         + frWtrIce(I,J)*ICE2WATR/SEAICE_deltaTtherm
512       &         )       &         )
513    
514  C NOW GET TOTAL QNET AND QSW  C NOW GET TOTAL QNET AND QSW
# Line 491  CML       CALL PLOT_FIELD_XYRL( FHEFF,'C Line 540  CML       CALL PLOT_FIELD_XYRL( FHEFF,'C
540         CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )         CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
541         CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )         CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
542         CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )         CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
543          DO j=1-OLy,sNy+OLy  CML       DO j=1-OLy,sNy+OLy
544           DO i=1-OLx,sNx+OLx  CML        DO i=1-OLx,sNx+OLx
545            GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2)  CML         GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2)
546            GAREA(I,J)=HEFF(I,J,1,bi,bj)  CML         GAREA(I,J)=HEFF(I,J,1,bi,bj)
547            print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj)  CML         print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj)
548           ENDDO  CML        ENDDO
549          ENDDO  CML       ENDDO
550         CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid )  CML       CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid )
551         CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid )  CML       CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid )
552          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
553           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
554            if(HEFF(i,j,1,bi,bj).gt.1.) then            if(HEFF(i,j,1,bi,bj).gt.1.) then
# Line 569  C     convert snow to ice if submerged Line 618  C     convert snow to ice if submerged
618       &              +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0       &              +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0
619             hFlood = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))             hFlood = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))
620             HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj) + hFlood             HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj) + hFlood
621             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,HSNOW(I,J,bi,bj)-hFlood/SDF)             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,
622         &          HSNOW(I,J,bi,bj)-hFlood*ICE2SNOW)
623            ENDDO            ENDDO
624           ENDDO           ENDDO
625          ENDIF          ENDIF
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22