/[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.2 by heimbach, Thu Dec 14 22:31:18 2006 UTC revision 1.3 by mlosch, Fri Dec 15 15:04:53 2006 UTC
# Line 38  C     === Local variables === Line 38  C     === Local variables ===
38  C     i,j,bi,bj - Loop counters  C     i,j,bi,bj - Loop counters
39    
40        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
41    C     number of surface interface layer
42          INTEGER kSurface
43        _RL TBC, salinity_ice, SDF, ICE_DENS, Q0, QS        _RL TBC, salinity_ice, SDF, ICE_DENS, Q0, QS
44  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
45        _RL hDraft, hFlood        _RL hDraft, hFlood
# Line 74  C     actual snow thickness Line 76  C     actual snow thickness
76  C     wind speed  C     wind speed
77        _RL UG     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL UG     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
78        _RL SPEED_SQ        _RL SPEED_SQ
79    C     local copy of AREA
80          _RL areaLoc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
81    
82  #ifdef SEAICE_MULTILEVEL  #ifdef SEAICE_MULTILEVEL
83        INTEGER it        INTEGER it
# Line 83  C     wind speed Line 87  C     wind speed
87        _RL FICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)        _RL FICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
88  #endif  #endif
89    
 C     number of surface interface layer  
       INTEGER kSurface  
   
90        if ( buoyancyRelation .eq. 'OCEANICP' ) then        if ( buoyancyRelation .eq. 'OCEANICP' ) then
91         kSurface        = Nr         kSurface        = Nr
92        else        else
93         kSurface        = 1         kSurface        = 1
94        endif        endif
95    
96        salinity_ice=4.0 _d 0      ! ICE SALINITY (g/kg)  C     ICE SALINITY (g/kg)
97        TBC=SEAICE_freeze          ! FREEZING TEMP. OF SEA WATER (deg C)        salinity_ice = 4.0 _d 0
98        SDF=1000.0 _d 0/330.0 _d 0 ! RATIO OF WATER DESITY TO SNOW DENSITY  C     FREEZING TEMP. OF SEA WATER (deg C)
99        ICE_DENS=0.920 _d 0        ! RATIO OF SEA ICE DESITY TO WATER DENSITY        TBC          = SEAICE_freeze
100        Q0=1.0D-06/302.0 _d +00    ! INVERSE HEAT OF FUSION OF ICE (m^3/J)  C     RATIO OF WATER DESITY TO SNOW DENSITY
101        QS=1.1 _d +08              ! HEAT OF FUSION OF SNOW (J/m^3)        SDF          = 1000.0 _d 0/330.0 _d 0
102    C     RATIO OF SEA ICE DESITY TO WATER DENSITY
103          ICE_DENS     = 0.920 _d 0
104    C     INVERSE HEAT OF FUSION OF ICE (m^3/J)
105          Q0           = 1.0 _d -06 / 302.0 _d +00
106    C     HEAT OF FUSION OF SNOW (J/m^3)
107          QS           = 1.1 _d +08
108    
109        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
110         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 128  CADJ &                         key = iic Line 135  CADJ &                         key = iic
135  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
136          DO J=1,sNy          DO J=1,sNy
137           DO I=1,sNx           DO I=1,sNx
138            AREA(I,J,3,bi,bj) = MAX(A22,AREA(I,J,2,bi,bj))            areaLoc(I,J)     = MAX(A22,AREA(I,J,2,bi,bj))
139            FHEFF(I,J)        = 0.0 _d 0            FHEFF(I,J)       = 0.0 _d 0
140            FICE (I,J)        = 0.0 _d 0            FICE (I,J)       = 0.0 _d 0
141  #ifdef SEAICE_MULTILEVEL  #ifdef SEAICE_MULTILEVEL
142            FICEP(I,J)        = 0.0 _d 0            FICEP(I,J)       = 0.0 _d 0
143  #endif  #endif
144            FHEFF(I,J)        = 0.0 _d 0            FHEFF(I,J)       = 0.0 _d 0
145            FICE (I,J)        = 0.0 _d 0            FICE (I,J)       = 0.0 _d 0
146            QNETO(I,J)        = 0.0 _d 0            QNETO(I,J)       = 0.0 _d 0
147            QNETI(I,J)        = 0.0 _d 0            QNETI(I,J)       = 0.0 _d 0
148            QSWO (I,J)        = 0.0 _d 0            QSWO (I,J)       = 0.0 _d 0
149            QSWI (I,J)        = 0.0 _d 0            QSWI (I,J)       = 0.0 _d 0
150            HCORR(I,J)        = 0.0 _d 0            HCORR(I,J)       = 0.0 _d 0
151            SEAICE_SALT(I,J)  = 0.0 _d 0            SEAICE_SALT(I,J) = 0.0 _d 0
152            RESID_HEAT (I,J)  = 0.0 _d 0            RESID_HEAT (I,J) = 0.0 _d 0
153           ENDDO           ENDDO
154          ENDDO          ENDDO
155  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 157  CADJ &                         key = iic Line 164  CADJ &                         key = iic
164           DO I=1,sNx           DO I=1,sNx
165  cph need to adjoint-store AREA again before using it in further init.  cph need to adjoint-store AREA again before using it in further init.
166  cph (all these initialisations involving AREA are nasty "non-linear")  cph (all these initialisations involving AREA are nasty "non-linear")
167            HICE(I,J)         = HEFF(I,J,2,bi,bj)/AREA(I,J,3,bi,bj)            HICE(I,J)    = HEFF(I,J,2,bi,bj)/areaLoc(I,J)
168            hSnwLoc(I,J)      = HSNOW(I,J,bi,bj)/AREA(I,J,3,bi,bj)            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc(I,J)
169           ENDDO           ENDDO
170          ENDDO          ENDDO
171    
# Line 190  CML#endif /* SEAICE_EXTERNAL_FORCING */ Line 197  CML#endif /* SEAICE_EXTERNAL_FORCING */
197           ENDDO           ENDDO
198          ENDDO          ENDDO
199    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,  
 CADJ &                         key = iicekey, byte = isbyte  
 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  
 CADJ &                         key = iicekey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         DO J=1,sNy  
          DO I=1,sNx  
 C--   Create or melt sea-ice so that first-level oceanic temperature  
 C     is approximately at the freezing point when there is sea-ice.  
 C     Initially the units of YNEG are m of sea-ice.  
 C     The factor dRf(1)/72.0764, used to convert temperature  
 C     change in deg K to m of sea-ice, is approximately:  
 C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)  
 C        * (density of sea-water = 1026 kg/m^3)  
 C        / (latent heat of fusion of sea-ice = 334000 J/kg)  
 C        / (density of sea-ice = 910 kg/m^3)  
 C     Negative YNEG leads to ice growth.  
 C     Positive YNEG leads to ice melting.  
           if ( .NOT. inAdMode ) then  
 #ifdef SEAICE_VARIABLE_FREEZING_POINT  
            TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0  
 #endif /* SEAICE_VARIABLE_FREEZING_POINT */  
            YNEG(I,J,bi,bj)=(theta(I,J,kSurface,bi,bj)-TBC)  
      &          *dRf(1)/72.0764 _d 0  
           else  
            YNEG(I,J,bi,bj)= 0.  
           endif  
           GHEFF(I,J)=HEFF(I,J,1,bi,bj)  
 C     Melt (YNEG>0) or create (YNEG<0) sea ice  
           HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj))  
           RESID_HEAT(I,J)  =YNEG(I,J,bi,bj)  
           YNEG(I,J,bi,bj)  =GHEFF(I,J)-HEFF(I,J,1,bi,bj)  
           SEAICE_SALT(I,J) =SEAICE_SALT(I,J)-YNEG(I,J,bi,bj)  
           RESID_HEAT(I,J)  =RESID_HEAT(I,J)-YNEG(I,J,bi,bj)  
 C     YNEG now contains m of ice melted (>0) or created (<0)  
 C     SEAICE_SALT contains m of ice melted (<0) or created (>0)  
 C     RESID_HEAT is residual heat above freezing in equivalent m of ice  
          ENDDO  
         ENDDO  
200    
201  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
202  cphCADJ STORE heff   = comlev1, key = ikey_dynamics  cphCADJ STORE heff   = comlev1, key = ikey_dynamics
# Line 243  CADJ STORE tices  = comlev1, key = ikey_ Line 210  CADJ STORE tices  = comlev1, key = ikey_
210  # endif  # endif
211  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
212    
 C GROWTH SUBROUTINE CALCULATES TOTAL GROWTH TENDENCIES,  
 C INCLUDING SNOWFALL  
 CML  beginning of groatb code  
   
213  C NOW DETERMINE GROWTH RATES  C NOW DETERMINE GROWTH RATES
214  C FIRST DO OPEN WATER  C FIRST DO OPEN WATER
215          CALL SEAICE_BUDGET_OCEAN(          CALL SEAICE_BUDGET_OCEAN(
# Line 296  C--  End loop over muli-levels Line 259  C--  End loop over muli-levels
259       O       FICE, QSWI,       O       FICE, QSWI,
260       I       bi, bj)       I       bi, bj)
261  #endif /* SEAICE_MULTILEVEL */  #endif /* SEAICE_MULTILEVEL */
 CML   end of groatb code  
262            
263    #ifdef ALLOW_AUTODIFF_TAMC
264    CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
265    CADJ &                         key = iicekey, byte = isbyte
266    CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
267    CADJ &                         key = iicekey, byte = isbyte
268    #endif /* ALLOW_AUTODIFF_TAMC */
269            DO J=1,sNy
270             DO I=1,sNx
271    C--   Create or melt sea-ice so that first-level oceanic temperature
272    C     is approximately at the freezing point when there is sea-ice.
273    C     Initially the units of YNEG are m of sea-ice.
274    C     The factor dRf(1)/72.0764, used to convert temperature
275    C     change in deg K to m of sea-ice, is approximately:
276    C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)
277    C        * (density of sea-water = 1026 kg/m^3)
278    C        / (latent heat of fusion of sea-ice = 334000 J/kg)
279    C        / (density of sea-ice = 910 kg/m^3)
280    C     Negative YNEG leads to ice growth.
281    C     Positive YNEG leads to ice melting.
282              IF ( .NOT. inAdMode ) THEN
283    #ifdef SEAICE_VARIABLE_FREEZING_POINT
284               TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
285    #endif /* SEAICE_VARIABLE_FREEZING_POINT */
286               YNEG(I,J,bi,bj) = (theta(I,J,kSurface,bi,bj)-TBC)
287         &          *dRf(1)/72.0764 _d 0
288              ELSE
289               YNEG(I,J,bi,bj)= 0.
290              ENDIF
291              GHEFF(I,J)=HEFF(I,J,1,bi,bj)
292    C     Melt (YNEG>0) or create (YNEG<0) sea ice
293              HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj))
294              RESID_HEAT(I,J)  = YNEG(I,J,bi,bj)
295              YNEG(I,J,bi,bj)  = GHEFF(I,J)-HEFF(I,J,1,bi,bj)
296              SEAICE_SALT(I,J) = SEAICE_SALT(I,J)-YNEG(I,J,bi,bj)
297              RESID_HEAT(I,J)  =  RESID_HEAT(I,J)-YNEG(I,J,bi,bj)
298    C     YNEG now contains m of ice melted (>0) or created (<0)
299    C     SEAICE_SALT contains m of ice melted (<0) or created (>0)
300    C     RESID_HEAT is residual heat above freezing in equivalent m of ice
301             ENDDO
302            ENDDO
303    
304  cph(  cph(
305  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
306  cphCADJ STORE heff   = comlev1, key = ikey_dynamics  cphCADJ STORE heff   = comlev1, key = ikey_dynamics
# Line 317  cph) Line 320  cph)
320    
321          DO J=1,sNy          DO J=1,sNy
322           DO I=1,sNx           DO I=1,sNx
323  C NOW CALCULATE CORRECTED GROWTH  C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)
324            GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)            GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)*AREA(I,J,2,bi,bj)
      &         *AREA(I,J,2,bi,bj)        ! effective growth in J/m^2 (>0=melt)  
325           ENDDO           ENDDO
326          ENDDO          ENDDO
327    
328  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
329  CADJ STORE fice(:,:)   = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
330  CADJ &     key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
331  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
332    
333          DO J=1,sNy          DO J=1,sNy
# Line 356  C     convert all snow to melt water (fr Line 358  C     convert all snow to melt water (fr
358          ENDDO          ENDDO
359    
360  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
361  CADJ STORE fice(:,:)   = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
362  CADJ &     key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
363  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
364    
365          DO J=1,sNy          DO J=1,sNy
# Line 373  CADJ STORE heff(:,:,:,bi,bj) = comlev1_b Line 375  CADJ STORE heff(:,:,:,bi,bj) = comlev1_b
375  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
376  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,
377  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
378  CADJ STORE fice(:,:)   = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
379  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
380  CADJ STORE fheff(:,:)  = comlev1_bibj,  CADJ STORE fheff(:,:)        = comlev1_bibj,
381  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
382  CADJ STORE qneto(:,:)  = comlev1_bibj,  CADJ STORE qneto(:,:)        = comlev1_bibj,
383  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
384  CADJ STORE qswi(:,:)   = comlev1_bibj,  CADJ STORE qswi(:,:)         = comlev1_bibj,
385  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
386  CADJ STORE qswo(:,:)   = comlev1_bibj,  CADJ STORE qswo(:,:)         = comlev1_bibj,
387  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
388  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
389  cph)  cph)
390          DO J=1,sNy          DO J=1,sNy
391           DO I=1,sNx           DO I=1,sNx
392  C NOW UPDATE AREA  C NOW UPDATE AREA
393            GHEFF(I,J)=-SEAICE_deltaTtherm*FHEFF(I,J)*Q0            GHEFF(I,J) = -SEAICE_deltaTtherm*FHEFF(I,J)*Q0
394            GAREA(I,J)=SEAICE_deltaTtherm*QNETO(I,J)*Q0            GAREA(I,J) =  SEAICE_deltaTtherm*QNETO(I,J)*Q0
395            GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))            GHEFF(I,J) = -ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
396            GAREA(I,J)=MAX(ZERO,GAREA(I,J))            GAREA(I,J) = MAX(ZERO,GAREA(I,J))
397            HCORR(I,J)=MIN(ZERO,GHEFF(I,J))            HCORR(I,J) = MIN(ZERO,GHEFF(I,J))
398           ENDDO           ENDDO
399          ENDDO          ENDDO
400  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 415  CADJ &                         key = iic Line 417  CADJ &                         key = iic
417           DO I=1,sNx           DO I=1,sNx
418    
419  C NOW UPDATE HEFF  C NOW UPDATE HEFF
420            GHEFF(I,J)=-SEAICE_deltaTtherm*            GHEFF(I,J)       = -SEAICE_deltaTtherm*
421       &               FICE(I,J)*Q0*AREA(I,J,2,bi,bj)       &                       FICE(I,J)*Q0*AREA(I,J,2,bi,bj)
422            GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))            GHEFF(I,J)       = -ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
423            HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)+GHEFF(I,J)            HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj)+GHEFF(I,J)
424            SEAICE_SALT(I,J)=SEAICE_SALT(I,J)+GHEFF(I,J)            SEAICE_SALT(I,J) = SEAICE_SALT(I,J)+GHEFF(I,J)
425    
426  C NOW CALCULATE QNETI UNDER ICE IF ANY  C NOW CALCULATE QNETI UNDER ICE IF ANY
427            QNETI(I,J)=(GHEFF(I,J)-SEAICE_deltaTtherm*            QNETI(I,J)       = (GHEFF(I,J)-SEAICE_deltaTtherm*
428       &       FICE(I,J)*Q0*AREA(I,J,2,bi,bj))/Q0/SEAICE_deltaTtherm       &       FICE(I,J)*Q0*AREA(I,J,2,bi,bj))/Q0/SEAICE_deltaTtherm
429    
430  C NOW UPDATE OTHER THINGS  C NOW UPDATE OTHER THINGS
431    
432            IF(FICE(I,J).GT.ZERO) THEN            IF(FICE(I,J).GT.ZERO) THEN
433  C FREEZING, PRECIP ADDED AS SNOW  C FREEZING, PRECIP ADDED AS SNOW
434             HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
435       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF       &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
436            ELSE            ELSE
437  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 ICE BY 1/ICE_DENS
438               SEAICE_SALT(I,J)=SEAICE_SALT(I,J)             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)
439       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*       &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
440       &            SEAICE_deltaTtherm/ICE_DENS       &            SEAICE_deltaTtherm/ICE_DENS
441            ENDIF            ENDIF
442    
443  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
444            SEAICE_SALT(I,J)=SEAICE_SALT(I,J)            SEAICE_SALT(I,J) = SEAICE_SALT(I,J)
445       &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))       &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
446       &         *SEAICE_deltaTtherm/ICE_DENS       &         *SEAICE_deltaTtherm/ICE_DENS
447    
448  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
449  C Note that units of YNEG and SEAICE_SALT are m of ice  C Note that units of YNEG and SEAICE_SALT are m of ice
450            IF(RESID_HEAT(I,J).GT.ZERO.AND.            IF( RESID_HEAT(I,J)      .GT.ZERO
451       &         HSNOW(I,J,bi,bj).GT.ZERO) THEN       &        .AND.HSNOW(I,J,bi,bj).GT.ZERO ) THEN
452             GHEFF(I,J)=MIN(HSNOW(I,J,bi,bj)/SDF/ICE_DENS,             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE_DENS,
453       &         RESID_HEAT(I,J))       &                             RESID_HEAT(I,J) )
454             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)
455             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
456             SEAICE_SALT(I,J)=SEAICE_SALT(I,J)-GHEFF(I,J)             SEAICE_SALT(I,J) = SEAICE_SALT(I,J)-GHEFF(I,J)
457            ENDIF            ENDIF
458    
459  C NOW GET FRESH WATER FLUX  C NOW GET FRESH WATER FLUX
460            EmPmR(I,J,bi,bj)= maskC(I,J,kSurface,bi,bj)*(            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
461       &         EVAP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))       &         EVAP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
462       &         -RUNOFF(I,J,bi,bj)       &         -RUNOFF(I,J,bi,bj)
463       &         +SEAICE_SALT(I,J)*ICE_DENS/SEAICE_deltaTtherm       &         +SEAICE_SALT(I,J)*ICE_DENS/SEAICE_deltaTtherm
464       &         )       &         )
465    
466  C NOW GET TOTAL QNET AND QSW  C NOW GET TOTAL QNET AND QSW
467            QNET(I,J,bi,bj)=QNETI(I,J) *AREA(I,J,2,bi,bj)            QNET(I,J,bi,bj) = QNETI(I,J) * AREA(I,J,2,bi,bj)
468       &                   +QNETO(I,J) *(ONE-AREA(I,J,2,bi,bj))       &                     +QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
469            QSW(I,J,bi,bj) =QSWI(I,J)  *AREA(I,J,2,bi,bj)            QSW(I,J,bi,bj)  = QSWI(I,J)  * AREA(I,J,2,bi,bj)
470       &                   +QSWO(I,J)  *(ONE-AREA(I,J,2,bi,bj))       &                     +QSWO(I,J)  * (ONE-AREA(I,J,2,bi,bj))
 c #ifndef SHORTWAVE_HEATING  
 c         QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+QSW(I,J,bi,bj)  
 c #endif  
471    
472  C Now convert YNEG back to deg K.  C Now convert YNEG back to deg K.
473            YNEG(I,J,bi,bj)=YNEG(I,J,bi,bj)*recip_dRf(1)*72.0764 _d 0            YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj)*recip_dRf(1)*72.0764 _d 0
474    
475  C Add YNEG contribution to QNET  C Add YNEG contribution to QNET
476            QNET(I,J,bi,bj)=QNET(I,J,bi,bj)            QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
477       &         +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm       &         +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm
478       &         *maskC(I,J,kSurface,bi,bj)       &         *maskC(I,J,kSurface,bi,bj)
479       &         *HeatCapacity_Cp*recip_horiVertRatio*rhoConst       &         *HeatCapacity_Cp*recip_horiVertRatio*rhoConst
# Line 550  CADJ &                         key = iic Line 549  CADJ &                         key = iic
549  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
550          DO J=1,sNy          DO J=1,sNy
551           DO I=1,sNx           DO I=1,sNx
552            AREA(I,J,1,bi,bj)=MAX(ZERO,AREA(I,J,1,bi,bj))            AREA(I,J,1,bi,bj) = MAX(ZERO,AREA(I,J,1,bi,bj))
553            HSNOW(I,J,bi,bj)=MAX(ZERO,HSNOW(I,J,bi,bj))            HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))
554  #endif  #endif
555            AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)            AREA(I,J,1,bi,bj) = AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
556            HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)            HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
557  #ifdef DO_WE_NEED_THIS  #ifdef DO_WE_NEED_THIS
558  c          HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))  c          HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))
559  #endif  #endif
560            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)            HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
561           ENDDO           ENDDO
562          ENDDO          ENDDO
563    
# Line 571  C     convert snow to ice if submerged Line 570  C     convert snow to ice if submerged
570       &              +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0       &              +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0
571             hFlood = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))             hFlood = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))
572             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
573             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,HSNOW(I,J,bi,bj)-hFlood/SDF)
574            ENDDO            ENDDO
575           ENDDO           ENDDO
576          ENDIF          ENDIF
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22