/[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.7 by mlosch, Wed Dec 20 12:25:15 2006 UTC revision 1.8 by mlosch, Wed Dec 20 20:49:12 2006 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, ICE_DENS, 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 with upper and lower limit  C     actual ice thickness with upper and lower limit
       _RL hIceLoc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
   
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 104  C     RATIO OF WATER DESITY TO SNOW DENS Line 107  C     RATIO OF WATER DESITY TO SNOW DENS
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        ICE_DENS     = 0.920 _d 0
110  C     INVERSE HEAT OF FUSION OF ICE (m^3/J)  C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
111        Q0           = 1.0 _d -06 / 302.0 _d +00        ICE2SNOW     = SDF * ICE_DENS
112    C     HEAT OF FUSION OF ICE (m^3/J)
113          QI           = 302.0 _d +06
114          recip_QI     = 1.0 _d 0 / QI
115  C     HEAT OF FUSION OF SNOW (J/m^3)  C     HEAT OF FUSION OF SNOW (J/m^3)
116        QS           = 1.1 _d +08        QS           = 1.1 _d +08
117          recip_QS     = 1.1 _d 0 / QS
118    
119        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
120         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 137  CADJ &                         key = iic Line 144  CADJ &                         key = iic
144  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
145          DO J=1,sNy          DO J=1,sNy
146           DO I=1,sNx           DO I=1,sNx
147            FHEFF(I,J)       = 0.0 _d 0            FHEFF(I,J)      = 0.0 _d 0
148            FICE (I,J)       = 0.0 _d 0            FICE (I,J)      = 0.0 _d 0
149  #ifdef SEAICE_MULTICATEGORY  #ifdef SEAICE_MULTICATEGORY
150            FICEP(I,J)       = 0.0 _d 0            FICEP(I,J)      = 0.0 _d 0
151            QSWIP(I,J)       = 0.0 _d 0            QSWIP(I,J)      = 0.0 _d 0
152  #endif  #endif
153            FHEFF(I,J)       = 0.0 _d 0            FHEFF(I,J)      = 0.0 _d 0
154            FICE (I,J)       = 0.0 _d 0            FICE (I,J)      = 0.0 _d 0
155            QNETO(I,J)       = 0.0 _d 0            QNETO(I,J)      = 0.0 _d 0
156            QNETI(I,J)       = 0.0 _d 0            QNETI(I,J)      = 0.0 _d 0
157            QSWO (I,J)       = 0.0 _d 0            QSWO (I,J)      = 0.0 _d 0
158            QSWI (I,J)       = 0.0 _d 0            QSWI (I,J)      = 0.0 _d 0
159            HCORR(I,J)       = 0.0 _d 0            HCORR(I,J)      = 0.0 _d 0
160            SEAICE_SALT(I,J) = 0.0 _d 0            frWtrIce(I,J)   = 0.0 _d 0
161            RESID_HEAT (I,J) = 0.0 _d 0            RESID_HEAT(I,J) = 0.0 _d 0
162           ENDDO           ENDDO
163          ENDDO          ENDDO
164  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 164  CADJ &                         key = iic Line 171  CADJ &                         key = iic
171           DO I=1,sNx           DO I=1,sNx
172  C     COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM  C     COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
173  C     ON ICE THICKNESS FOR BUDGET COMPUTATION  C     ON ICE THICKNESS FOR BUDGET COMPUTATION
174    C     The default of A22 = 0.15 is a common threshold for defining
175    C     the ice edge. This ice concentration usually does not occur
176    C     due to thermodynamics but due to advection.
177            areaLoc      = MAX(A22,AREA(I,J,2,bi,bj))            areaLoc      = MAX(A22,AREA(I,J,2,bi,bj))
178            HICE(I,J)    = HEFF(I,J,2,bi,bj)/areaLoc            HICE(I,J)    = HEFF(I,J,2,bi,bj)/areaLoc
179    C     Do we know what this is for?
180            HICE(I,J)    = MAX(HICE(I,J),0.05 _d +00)            HICE(I,J)    = MAX(HICE(I,J),0.05 _d +00)
181    C     Capping the actual ice thickness effectively enforces a
182    C     minimum of heat flux through the ice and helps getting rid of
183    C     very thick ice.
184            HICE(I,J)    = MIN(HICE(I,J),9.0 _d +00)            HICE(I,J)    = MIN(HICE(I,J),9.0 _d +00)
185            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
186           ENDDO           ENDDO
# Line 266  CADJ &                         key = iic Line 280  CADJ &                         key = iic
280  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
281  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
282  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
283    C
284    C--   compute and apply ice growth due to oceanic heat flux from below
285    C
286          DO J=1,sNy          DO J=1,sNy
287           DO I=1,sNx           DO I=1,sNx
288  C--   Create or melt sea-ice so that first-level oceanic temperature  C--   Create or melt sea-ice so that first-level oceanic temperature
289  C     is approximately at the freezing point when there is sea-ice.  C     is approximately at the freezing point when there is sea-ice.
290  C     Initially the units of YNEG are m of sea-ice.  C     Initially the units of YNEG/availHeat are m of sea-ice.
291  C     The factor dRf(1)/72.0764, used to convert temperature  C     The factor dRf(1)/72.0764, used to convert temperature
292  C     change in deg K to m of sea-ice, is approximately:  C     change in deg K to m of sea-ice, is approximately:
293  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)
294  C        * (density of sea-water = 1026 kg/m^3)  C        * (density of sea-water = 1026 kg/m^3)
295  C        / (latent heat of fusion of sea-ice = 334000 J/kg)  C        / (latent heat of fusion of sea-ice = 334000 J/kg)
296  C        / (density of sea-ice = 910 kg/m^3)  C        / (density of sea-ice = 910 kg/m^3)
297  C     Negative YNEG leads to ice growth.  C     Negative YNEG/availHeat leads to ice growth.
298  C     Positive YNEG leads to ice melting.  C     Positive YNEG/availHeat leads to ice melting.
299            IF ( .NOT. inAdMode ) THEN            IF ( .NOT. inAdMode ) THEN
300  #ifdef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
301             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
302  #endif /* SEAICE_VARIABLE_FREEZING_POINT */  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
303             YNEG(I,J,bi,bj) = (theta(I,J,kSurface,bi,bj)-TBC)             availHeat = (theta(I,J,kSurface,bi,bj)-TBC)
304       &          *dRf(1)/72.0764 _d 0       &          *dRf(1)/72.0764 _d 0
305            ELSE            ELSE
306             YNEG(I,J,bi,bj)= 0.             availHeat = 0.
307            ENDIF            ENDIF
308            GHEFF(I,J)=HEFF(I,J,1,bi,bj)  C     local copy of old effective ice thickness
309  C     Melt (YNEG>0) or create (YNEG<0) sea ice            hEffOld           = HEFF(I,J,1,bi,bj)
310            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
311            RESID_HEAT(I,J)  = YNEG(I,J,bi,bj)            HEFF(I,J,1,bi,bj) = MAX(ZERO,HEFF(I,J,1,bi,bj)-availHeat)
312            YNEG(I,J,bi,bj)  = GHEFF(I,J)-HEFF(I,J,1,bi,bj)  C    
313            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)
314            RESID_HEAT(I,J)  =  RESID_HEAT(I,J)-YNEG(I,J,bi,bj)  C    
315              frWtrIce(I,J)     = frWtrIce(I,J) - YNEG(I,J,bi,bj)
316              RESID_HEAT(I,J)   = availHeat     - YNEG(I,J,bi,bj)
317  C     YNEG now contains m of ice melted (>0) or created (<0)  C     YNEG now contains m of ice melted (>0) or created (<0)
318  C     SEAICE_SALT contains m of ice melted (<0) or created (>0)  C     frWtrIce contains m of ice melted (<0) or created (>0)
319  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
320           ENDDO           ENDDO
321          ENDDO          ENDDO
# Line 317  CADJ STORE fice(:,:)         = comlev1_b Line 336  CADJ STORE fice(:,:)         = comlev1_b
336  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
337  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
338  cph)  cph)
339    C
340    C--   compute and apply ice growth due to atmospheric fluxes from above
341    C    
342          DO J=1,sNy          DO J=1,sNy
343           DO I=1,sNx           DO I=1,sNx
344  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 333  CADJ &                         key = iic Line 354  CADJ &                         key = iic
354          DO J=1,sNy          DO J=1,sNy
355           DO I=1,sNx           DO I=1,sNx
356            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
357  C use FICE to melt snow and CALCULATE CORRECTED GROWTH  C     use FICE to melt snow and CALCULATE CORRECTED GROWTH
358             GAREA(I,J)=HSNOW(I,J,bi,bj)*QS ! effective snow thickness in J/m^2  C     effective snow thickness in J/m^2
359             IF(GHEFF(I,J).LE.GAREA(I,J)) THEN             snowEnergy=HSNOW(I,J,bi,bj)*QS
360  C not enough heat to melt all snow; use up all heat flux FICE             IF(GHEFF(I,J).LE.snowEnergy) THEN
361    C     not enough heat to melt all snow; use up all heat flux FICE
362              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
363  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
364  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
365              SEAICE_SALT(I,J)=SEAICE_SALT(I,J)              frWtrIce(I,J) = frWtrIce(I,J) - GHEFF(I,J)/(QS*ICE2SNOW)
366       &           -GHEFF(I,J)/QS/SDF/ICE_DENS              FICE    (I,J) = ZERO
             FICE(I,J)=ZERO  
367             ELSE             ELSE
368  C enought heat to melt snow completely;  C     enought heat to melt snow completely;
369  C compute remaining heat flux that will melt ice  C     compute remaining heat flux that will melt ice
370              FICE(I,J)=-(GHEFF(I,J)-GAREA(I,J))/              FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/
371       &           SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)       &           SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)
372  C     convert all snow to melt water (fresh water flux)  C     convert all snow to melt water (fresh water flux)
373              SEAICE_SALT(I,J)=SEAICE_SALT(I,J)              frWtrIce(I,J) = frWtrIce(I,J)
374       &           -HSNOW(I,J,bi,bj)/SDF/ICE_DENS       &           -HSNOW(I,J,bi,bj)/ICE2SNOW
375              HSNOW(I,J,bi,bj)=0.0              HSNOW(I,J,bi,bj)=0.0
376             END IF             END IF
377            END IF            END IF
# Line 364  CADJ &                         key = iic Line 385  CADJ &                         key = iic
385    
386          DO J=1,sNy          DO J=1,sNy
387           DO I=1,sNx           DO I=1,sNx
388  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
389            FHEFF(I,J)= FICE(I,J)  * AREA(I,J,2,bi,bj)            FHEFF(I,J)= FICE(I,J)  * AREA(I,J,2,bi,bj)
390       &              + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))       &              + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
391           ENDDO           ENDDO
# Line 387  CADJ STORE qswo(:,:)         = comlev1_b Line 408  CADJ STORE qswo(:,:)         = comlev1_b
408  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
409  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
410  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  
411  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
412  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
413  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
414  #endif  #endif
415    C
416    C     First update (freeze or melt) ice area
417    C
418          DO J=1,sNy          DO J=1,sNy
419           DO I=1,sNx           DO I=1,sNx
420            GAREA(I,J)=(ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO  C     negative growth in meters of ice (>0 for melting)
421       &    +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
422       &    /(HEFF(I,J,1,bi,bj)+.00001 _d 0)  C     negative growth must not exceed effective ice thickness (=volume)
423            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)
424              growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
425    C     growthHEFF < 0 means melting
426              HCORR(I,J) = MIN(ZERO,growthHEFF)
427    C     gain of new effective ice thickness over open water (>0 by definition)
428              GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)
429    CML   removed these loops and moved TAMC store directive up
430    CML         ENDDO
431    CML        ENDDO
432    CML#ifdef ALLOW_AUTODIFF_TAMC
433    CMLCADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
434    CMLCADJ &                         key = iicekey, byte = isbyte
435    CML#endif
436    CML        DO J=1,sNy
437    CML         DO I=1,sNx
438    C     Here we finally compute the new AREA
439              AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+
440         &         (ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
441         &         +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)
442         &         /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
443           ENDDO           ENDDO
444          ENDDO          ENDDO
445  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
446  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
447  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
448  #endif  #endif
449    C    
450    C     now update (freeze or melt) HEFF
451    C        
452          DO J=1,sNy          DO J=1,sNy
453           DO I=1,sNx           DO I=1,sNx
454    C     negative growth (>0 for melting) of existing ice in meters
455              growthNeg        = -SEAICE_deltaTtherm*
456         &                       FICE(I,J)*recip_QI*AREA(I,J,2,bi,bj)
457    C     negative growth must not exceed effective ice thickness (=volume)
458    C     (that is, cannot melt more than all the ice)
459              growthHEFF       = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
460    C     growthHEFF < 0 means melting
461              HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj) + growthHEFF
462    C     add effective growth to fresh water of ice
463              frWtrIce(I,J)    = frWtrIce(I,J)     + growthHEFF
464    
465    C     now calculate QNETI under ice (if any) as the difference between
466    C     the available "heat flux" growthNeg and the actual growthHEFF;
467    C     keep in mind that growthNeg and growthHEFF have different signs
468    C     by construction
469              QNETI(I,J) =  (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm
470    
471  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  
472    
473            IF(FICE(I,J).GT.ZERO) THEN            IF(FICE(I,J).GT.ZERO) THEN
474  C FREEZING, PRECIP ADDED AS SNOW  C     freezing, add precip as snow
475             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
476       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
477            ELSE            ELSE
478  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
479             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)  C     ice by 1/ICE_DENS.
480    C     Sign issue? frWtrIce (>0 for more ice)
481    C     has the opposite sign of precip (>0 for downward flux of water)
482               frWtrIce(I,J) = frWtrIce(I,J)
483       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
484       &            SEAICE_deltaTtherm/ICE_DENS       &            SEAICE_deltaTtherm/ICE_DENS
485            ENDIF            ENDIF
486    
487  C Now add in precip over open water directly into ocean as negative salt  C     now add in precip over open water directly into ocean as negative salt
488            SEAICE_SALT(I,J) = SEAICE_SALT(I,J)  C     Sign issue? frWtrIce (>0 for more ice)
489    C     has the opposite sign of precip (>0 for downward flux of water)
490              frWtrIce(I,J) = frWtrIce(I,J)
491       &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))       &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
492       &         *SEAICE_deltaTtherm/ICE_DENS       &         *SEAICE_deltaTtherm/ICE_DENS
493    
494  C Now melt snow if there is residual heat left in surface level  C     Now melt snow if there is residual heat left in surface level
495  C Note that units of YNEG and SEAICE_SALT are m of ice  C     Note that units of YNEG and frWtrIce are m of ice
496  cph( very sensitive bit here by JZ  cph( very sensitive bit here by JZ
497            IF( RESID_HEAT(I,J) .GT. ZERO            IF( RESID_HEAT(I,J) .GT. ZERO .AND.
498       &        .AND. HSNOW(I,J,bi,bj) .GT. ZERO ) THEN       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
499             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE_DENS,             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE_DENS,
500       &                             RESID_HEAT(I,J) )       &                             RESID_HEAT(I,J) )
501             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)
502             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*ICE_DENS
503             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)-GHEFF(I,J)             frWtrIce(I,J)    = frWtrIce(I,J)   -GHEFF(I,J)
504            ENDIF            ENDIF
505  cph)  cph)
506    
# Line 462  C NOW GET FRESH WATER FLUX Line 508  C NOW GET FRESH WATER FLUX
508            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
509       &         EVAP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))       &         EVAP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
510       &         -RUNOFF(I,J,bi,bj)       &         -RUNOFF(I,J,bi,bj)
511       &         +SEAICE_SALT(I,J)*ICE_DENS/SEAICE_deltaTtherm       &         +frWtrIce(I,J)*ICE_DENS/SEAICE_deltaTtherm
512       &         )       &         )
513    
514  C NOW GET TOTAL QNET AND QSW  C NOW GET TOTAL QNET AND QSW
# Line 494  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
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22