/[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.71 by gforget, Thu Sep 30 02:22:25 2010 UTC revision 1.72 by gforget, Fri Oct 1 01:45:15 2010 UTC
# Line 61  C     auxillary variables Line 61  C     auxillary variables
61  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
62        _RL hDraft        _RL hDraft
63  #endif /* ALLOW_SEAICE_FLOODING */  #endif /* ALLOW_SEAICE_FLOODING */
       _RL hEffOld       (1:sNx,1:sNy)  
64        _RL GHEFF         (1:sNx,1:sNy)        _RL GHEFF         (1:sNx,1:sNy)
65  C     heat above freezing in equivalent m of ice        
       _RL a_HEFFbyICEonOCN     (1:sNx,1:sNy)  
       _RL r_HEFFbyICEonOCN    (1:sNx,1:sNy)  
 cgf =====  
 cgf YNEG is defined in common block -- for now I redefine/initialize it in this routine  
       _RL d_HEFFbyICEonOCN    (1:sNx,1:sNy,nSx,nSy)  
 cgf =====  
66  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
67        _RL saltFluxAdjust(1:sNx,1:sNy)        _RL saltFluxAdjust(1:sNx,1:sNy)
68  #endif  #endif
# Line 82  c Line 75  c
75        _RL d_AREAbyATM    (1:sNx,1:sNy)        _RL d_AREAbyATM    (1:sNx,1:sNy)
76  c  c
77        _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4        _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4
78    c
79          _RL a_QbyICE    (1:sNx,1:sNy)
80          _RL r_QbyICE    (1:sNx,1:sNy)
81          _RL d_QbyICE    (1:sNx,1:sNy)
82          _RL d_QbySNW    (1:sNx,1:sNy)
83          _RL d_HEFFbyICEonOCN    (1:sNx,1:sNy)
84          _RL d_HEFFbySNWonOCN    (1:sNx,1:sNy)
85          
86    
87    
88  C     a_QbyATM_cover  - thermodynamic ice growth rate over sea ice in W/m^2  C     a_QbyATM_cover  - thermodynamic ice growth rate over sea ice in W/m^2
# Line 154  C     HEAT OF FUSION OF ICE (J/m^3) Line 155  C     HEAT OF FUSION OF ICE (J/m^3)
155  C     HEAT OF FUSION OF SNOW (J/m^3)  C     HEAT OF FUSION OF SNOW (J/m^3)
156        QS           = 1.1 _d +08        QS           = 1.1 _d +08
157    
 cgf =====  
 cgf YNEG is defined in common block -- for now I redefine/initialize it in this routine  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           d_HEFFbyICEonOCN(I,J,bi,bj)      = 0.0 _d 0  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
 cgf =====  
158    
159        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
160         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 201  CADJ &                       key = iicek Line 190  CADJ &                       key = iicek
190            QNETI(I,J)      = 0.0 _d 0            QNETI(I,J)      = 0.0 _d 0
191            a_QSWbyATM_open (I,J)      = 0.0 _d 0            a_QSWbyATM_open (I,J)      = 0.0 _d 0
192            a_QSWbyATM_cover (I,J)      = 0.0 _d 0            a_QSWbyATM_cover (I,J)      = 0.0 _d 0
           r_HEFFbyICEonOCN(I,J) = 0.0 _d 0  
193  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
194            saltFluxAdjust(I,J) = 0.0 _d 0            saltFluxAdjust(I,J) = 0.0 _d 0
195  #endif  #endif
# Line 215  cgf new variables: Line 203  cgf new variables:
203            d_QbyATMonSNW(I,J)      = 0.0 _d 0            d_QbyATMonSNW(I,J)      = 0.0 _d 0
204  c  c
205            d_AREAbyATM(I,J)      = 0.0 _d 0            d_AREAbyATM(I,J)      = 0.0 _d 0
206    c
207              d_HEFFbyICEonOCN(I,J)      = 0.0 _d 0
208              d_HEFFbySNWonOCN(I,J)      = 0.0 _d 0
209           ENDDO           ENDDO
210          ENDDO          ENDDO
211          DO J=1-oLy,sNy+oLy          DO J=1-oLy,sNy+oLy
# Line 399  CADJ &                          key = ii Line 390  CADJ &                          key = ii
390  C  C
391  C--   compute and apply ice growth due to oceanic heat flux from below  C--   compute and apply ice growth due to oceanic heat flux from below
392  C  C
393    
394    c a_QbyICE: available heat that may be extracted by the ICE acting
395    c on the OCN surface/mixed layer to bring it back to freezig point
396    c (in m, negative out of the ocean -- different convention than a_QbyATM)
397    cgf unit change that breaks lab_sea (in W/m2, positive out of the ocean -- same convention as a_QbyATM)
398          DO J=1,sNy          DO J=1,sNy
399           DO I=1,sNx           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 d_HEFFbyICEonOCN/a_HEFFbyICEonOCN 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 d_HEFFbyICEonOCN/a_HEFFbyICEonOCN leads to ice growth.  
 C     Positive d_HEFFbyICEonOCN/a_HEFFbyICEonOCN leads to ice melting.  
400            IF ( .NOT. inAdMode ) THEN            IF ( .NOT. inAdMode ) THEN
401  #ifdef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
402             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
403  #endif /* SEAICE_VARIABLE_FREEZING_POINT */  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
404             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
405                a_HEFFbyICEonOCN(i,j) = SEAICE_availHeatFrac                a_QbyICE(i,j) = SEAICE_availHeatFrac
406       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
407       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
408    #ifdef SEAICE_QBYICE_IS_WPERM2
409    cgf using W/m2 as the Q unit changes lab_sea results...
410         &             * ( - QI / SEAICE_deltaTtherm )
411    #endif
412             ELSE             ELSE
413                a_HEFFbyICEonOCN(i,j) = SEAICE_availHeatFracFrz                a_QbyICE(i,j) = SEAICE_availHeatFracFrz
414       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
415       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
416    #ifdef SEAICE_QBYICE_IS_WPERM2
417         &             * ( - QI / SEAICE_deltaTtherm )
418    #endif
419             ENDIF             ENDIF
420            ELSE            ELSE
421             a_HEFFbyICEonOCN(i,j) = 0.             a_QbyICE(i,j) = 0.
422            ENDIF            ENDIF
423  C     local copy of old effective ice thickness           ENDDO
424            hEffOld(I,J)      = HEFF(I,J,bi,bj)          ENDDO
425  C     Melt (a_HEFFbyICEonOCN>0) or create (a_HEFFbyICEonOCN<0) sea ice  
426            HEFF(I,J,bi,bj) = MAX(ZERO,          DO J=1,sNy
427       &       HEFF(I,J,bi,bj)-a_HEFFbyICEonOCN(I,J))           DO I=1,sNx
428  C           tmpscal1=a_QbyICE(i,j)
429            d_HEFFbyICEonOCN(I,J,bi,bj) = hEffOld(I,J) - HEFF(I,J,bi,bj)  #ifdef SEAICE_QBYICE_IS_WPERM2
430  C       &             / ( - QI / SEAICE_deltaTtherm )
431    #endif
432              d_HEFFbyICEonOCN(I,J) =
433         &     MAX(ZERO, HEFF(I,J,bi,bj)-tmpscal1)- HEFF(I,J,bi,bj)
434              d_QbyICE(I,J)=-d_HEFFbyICEonOCN(I,J)
435    #ifdef SEAICE_QBYICE_IS_WPERM2
436         &             * ( - QI / SEAICE_deltaTtherm )
437    #endif
438              r_QbyICE(I,J)=a_QbyICE(I,J)-d_QbyICE(I,J)
439    c
440              HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyICEonOCN(I,J)
441            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)
442       &                              - d_HEFFbyICEonOCN(I,J,bi,bj)       &                              + d_HEFFbyICEonOCN(I,J)
           r_HEFFbyICEonOCN(I,J)   = a_HEFFbyICEonOCN(I,J)    
      &                              - d_HEFFbyICEonOCN(I,J,bi,bj)  
443  C     d_HEFFbyICEonOCN now contains m of ice melted (>0) or created (<0)  C     d_HEFFbyICEonOCN now contains m of ice melted (>0) or created (<0)
444  C     saltWtrIce contains m of ice melted (<0) or created (>0)  C     saltWtrIce contains m of ice melted (<0) or created (>0)
445  C     r_HEFFbyICEonOCN is residual heat above freezing in equivalent m of ice  C     r_QbyICE is residual heat above freezing in equivalent m of ice
446           ENDDO           ENDDO
447          ENDDO          ENDDO
448    
# Line 504  C     The factor 1/ICE2SNOW converts m o Line 504  C     The factor 1/ICE2SNOW converts m o
504             ELSE             ELSE
505  C     enought heat to melt snow completely;  C     enought heat to melt snow completely;
506  C     compute remaining heat flux that will melt ice  C     compute remaining heat flux that will melt ice
507              a_QbyATM_cover(I,J)=snowEnergy/              d_QbyATMonSNW(I,J)=-(GHEFF(I,J)-snowEnergy)/
508       &           SEAICE_deltaTtherm/AREANm1(I,J,bi,bj)       &        SEAICE_deltaTtherm/AREANm1(I,J,bi,bj)-a_QbyATM_cover(I,J)
509  C     convert all snow to melt water (fresh water flux)  C     convert all snow to melt water (fresh water flux)
510              d_HFRWbyATMonSNW(I,J)=-HSNOW(I,J,bi,bj)/ICE2SNOW              d_HFRWbyATMonSNW(I,J)=-HSNOW(I,J,bi,bj)/ICE2SNOW
511              d_HSNWbyATMonSNW(I,J)=-HSNOW(I,J,bi,bj)              d_HSNWbyATMonSNW(I,J)=-HSNOW(I,J,bi,bj)
# Line 659  cph( very sensitive bit here by JZ Line 659  cph( very sensitive bit here by JZ
659           DO I=1,sNx           DO I=1,sNx
660  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
661  C     Note that units of d_HEFFbyICEonOCN and frWtrIce are m of ice  C     Note that units of d_HEFFbyICEonOCN and frWtrIce are m of ice
           IF( r_HEFFbyICEonOCN(I,J) .GT. ZERO .AND.  
      &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN  
            GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,  
      &                             r_HEFFbyICEonOCN(I,J) )  
            d_HEFFbyICEonOCN(I,J,bi,bj)  =  
      &        d_HEFFbyICEonOCN(I,J,bi,bj) +GHEFF(I,J)  
           ENDIF  
          ENDDO  
         ENDDO  
662    
663  #ifdef ALLOW_AUTODIFF_TAMC           tmpscal1=r_QbyICE(i,j)
664  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  #ifdef SEAICE_QBYICE_IS_WPERM2
665  CADJ &                        key = iicekey, byte = isbyte       &             / ( - QI / SEAICE_deltaTtherm )
 CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,  
 CADJ &                        key = iicekey, byte = isbyte  
666  #endif  #endif
667          DO J=1,sNy            IF( tmpscal1 .GT. ZERO .AND.
          DO I=1,sNx  
           IF( r_HEFFbyICEonOCN(I,J) .GT. ZERO .AND.  
668       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
669             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR             d_HEFFbySNWonOCN(I,J) =
670             frWtrIce(I,J,bi,bj)    = frWtrIce(I,J,bi,bj)-GHEFF(I,J)       &       - MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR , tmpscal1 )
671              ELSE
672               d_HEFFbySNWonOCN(I,J)       = 0. _d 0
673            ENDIF            ENDIF
674              d_QbySNW(I,J)=-d_HEFFbySNWonOCN(I,J)
675    #ifdef SEAICE_QBYICE_IS_WPERM2
676         &             * ( - QI / SEAICE_deltaTtherm )
677    #endif
678              r_QbyICE(I,J)=a_QbyICE(I,J)-d_QbySNW(I,J)
679              HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)
680         &        +d_HEFFbySNWonOCN(I,J)*SDF*ICE2WATR
681              frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
682         &        +d_HEFFbySNWonOCN(I,J)
683           ENDDO           ENDDO
684          ENDDO          ENDDO
685  #endif  #endif
# Line 846  C NOW GET TOTAL QNET AND QSW Line 844  C NOW GET TOTAL QNET AND QSW
844          ENDIF          ENDIF
845  #endif  #endif
846    
 #ifdef ALLOW_AUTODIFF_TAMC  
 cphCADJ STORE d_HEFFbyICEonOCN(:,:,bi,bj) = comlev1_bibj,  
 cphCADJ &                       key = iicekey, byte = isbyte  
 #endif  
         DO J=1,sNy  
          DO I=1,sNx  
 C Now convert d_HEFFbyICEonOCN back to deg K.  
           d_HEFFbyICEonOCN(I,J,bi,bj) =  
      &       d_HEFFbyICEonOCN(I,J,bi,bj)*recip_dRf(kSurface) *  
      &       recip_hFacC(i,j,kSurface,bi,bj)*72.0764 _d 0  
          ENDDO  
         ENDDO  
847    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE d_HEFFbyICEonOCN(:,:,bi,bj) = comlev1_bibj,  
 CADJ &                       key = iicekey, byte = isbyte  
 #endif  
848          DO J=1,sNy          DO J=1,sNy
849           DO I=1,sNx           DO I=1,sNx
850    
851  C Add d_HEFFbyICEonOCN contribution to QNET  C Add d_HEFFbyICEonOCN contribution to QNET
852    #ifdef SEAICE_QBYICE_IS_WPERM2
853              tmpscal1=
854         &         -maskC(I,J,kSurface,bi,bj)
855         &         *(HeatCapacity_Cp*rUnit2mass/QI)
856         &         *72.0764 _d 0
857              QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
858         &         + ( d_QbyICE(I,J) + d_QbySNW(I,J) )
859         &         * tmpscal1
860    #else
861              tmpscal1 =
862         &       - ( d_HEFFbyICEonOCN(I,J)+d_HEFFbySNWonOCN(I,J) )
863         &       *recip_dRf(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
864         &       *72.0764 _d 0
865            QNET(I,J,bi,bj) = QNET(I,J,bi,bj)            QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
866       &         +d_HEFFbyICEonOCN(I,J,bi,bj)/SEAICE_deltaTtherm       &         +tmpscal1
867       &         *maskC(I,J,kSurface,bi,bj)       &         /SEAICE_deltaTtherm*maskC(I,J,kSurface,bi,bj)
868       &         *HeatCapacity_Cp*rUnit2mass       &         *HeatCapacity_Cp*rUnit2mass
869       &         *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)       &         *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
870    #endif
871    
872           ENDDO           ENDDO
873          ENDDO          ENDDO
874    
Legend:
Removed from v.1.71  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22