/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_growth.F
ViewVC logotype

Diff of /MITgcm_contrib/torge/itd/code/seaice_growth.F

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

revision 1.15 by torge, Mon Dec 10 22:41:54 2012 UTC revision 1.16 by torge, Wed Mar 27 18:59:52 2013 UTC
# Line 458  C ====================================== Line 458  C ======================================
458  C ===========PART 1: treat pathological cases (post advdiff)===========  C ===========PART 1: treat pathological cases (post advdiff)===========
459  C =====================================================================  C =====================================================================
460    
 #ifdef SEAICE_GROWTH_LEGACY  
   
         DO J=1,sNy  
          DO I=1,sNx  
           HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)  
           HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)  
           AREApreTH(I,J)=AREANM1(I,J,bi,bj)  
           d_HEFFbyNEG(I,J)=0. _d 0  
           d_HSNWbyNEG(I,J)=0. _d 0  
 #ifdef ALLOW_DIAGNOSTICS  
           DIAGarrayA(I,J) = AREANM1(I,J,bi,bj)  
           DIAGarrayB(I,J) = AREANM1(I,J,bi,bj)  
           DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj)  
           DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)  
 #endif  
          ENDDO  
         ENDDO  
 #ifdef SEAICE_ITD  
         DO IT=1,nITD  
          DO J=1,sNy  
           DO I=1,sNx  
            HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)  
            HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)  
            AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
 #endif  
   
 #else /* SEAICE_GROWTH_LEGACY */  
   
461  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
462  Cgf no dependency through pathological cases treatment  Cgf no dependency through pathological cases treatment
463          IF ( SEAICEadjMODE.EQ.0 ) THEN          IF ( SEAICEadjMODE.EQ.0 ) THEN
# Line 829  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 798  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
798          ENDDO          ENDDO
799  #endif /* SEAICE_VARIABLE_SALINITY */  #endif /* SEAICE_VARIABLE_SALINITY */
800    
 #endif /* SEAICE_GROWTH_LEGACY */  
   
801  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
802          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
803           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
# Line 902  CADJ STORE HSNWpreTH = comlev1_bibj, key Line 869  CADJ STORE HSNWpreTH = comlev1_bibj, key
869    
870  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
871            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
 #ifdef SEAICE_GROWTH_LEGACY  
             tmpscal1          = MAX(SEAICE_area_reg/float(nITD),  
      &                              AREAITDpreTH(I,J,IT))  
             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1  
             tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1  
             heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)  
 #else /* SEAICE_GROWTH_LEGACY */  
872  cif        regularize AREA with SEAICE_area_reg  cif        regularize AREA with SEAICE_area_reg
873             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
874       &                     + area_reg_sq)       &                     + area_reg_sq)
# Line 919  cif        regularize heffActual with SE Line 879  cif        regularize heffActual with SE
879       &                                  + hice_reg_sq)       &                                  + hice_reg_sq)
880  cif        hsnowActual calculated with the regularized AREA  cif        hsnowActual calculated with the regularized AREA
881             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
 #endif /* SEAICE_GROWTH_LEGACY */  
882  cif        regularize the inverse of heffActual by hice_reg  cif        regularize the inverse of heffActual by hice_reg
883             recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /             recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
884       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
# Line 932  cif       Do not regularize when HEFFpre Line 891  cif       Do not regularize when HEFFpre
891            ENDIF            ENDIF
892  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
893            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
 #ifdef SEAICE_GROWTH_LEGACY  
             tmpscal1         = MAX(SEAICE_area_reg,AREApreTH(I,J))  
             hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1  
             tmpscal2         = HEFFpreTH(I,J)/tmpscal1  
             heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)  
 #else /* SEAICE_GROWTH_LEGACY */  
894  Cif        regularize AREA with SEAICE_area_reg  Cif        regularize AREA with SEAICE_area_reg
895             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
896  Cif        heffActual calculated with the regularized AREA  Cif        heffActual calculated with the regularized AREA
# Line 946  Cif        regularize heffActual with SE Line 899  Cif        regularize heffActual with SE
899             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
900  Cif        hsnowActual calculated with the regularized AREA  Cif        hsnowActual calculated with the regularized AREA
901             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
 #endif /* SEAICE_GROWTH_LEGACY */  
902  Cif        regularize the inverse of heffActual by hice_reg  Cif        regularize the inverse of heffActual by hice_reg
903             recip_heffActual(I,J)  = AREApreTH(I,J) /             recip_heffActual(I,J)  = AREApreTH(I,J) /
904       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
# Line 1605  Cgf warming conditions, the lab_sea resu Line 1557  Cgf warming conditions, the lab_sea resu
1557             tmpscal1 = HEFFITDpreTH(I,J,IT)             tmpscal1 = HEFFITDpreTH(I,J,IT)
1558       &              + d_HEFFbySublim_ITD(I,J,IT)       &              + d_HEFFbySublim_ITD(I,J,IT)
1559       &              + d_HEFFbyOCNonICE_ITD(I,J,IT)       &              + d_HEFFbyOCNonICE_ITD(I,J,IT)
 #ifdef SEAICE_GROWTH_LEGACY  
            tmpscal2 = MAX(-tmpscal1,  
      &                     r_QbyATMmult_cover(I,J,IT))  
 #else  
1560             tmpscal2 = MAX(-tmpscal1,             tmpscal2 = MAX(-tmpscal1,
1561       &                     r_QbyATMmult_cover(I,J,IT)       &                     r_QbyATMmult_cover(I,J,IT)
1562  c         Limit ice growth by potential melt by ocean  c         Limit ice growth by potential melt by ocean
1563       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
 #endif /* SEAICE_GROWTH_LEGACY */  
1564             d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2             d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2
1565             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1566       &                                 + tmpscal2       &                                 + tmpscal2
# Line 1638  c         Limit ice growth by potential Line 1585  c         Limit ice growth by potential
1585          DO J=1,sNy          DO J=1,sNy
1586           DO I=1,sNx           DO I=1,sNx
1587    
 #ifdef SEAICE_GROWTH_LEGACY  
           tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))  
 #else  
1588            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1589  C         Limit ice growth by potential melt by ocean  C         Limit ice growth by potential melt by ocean
1590       &        AREApreTH(I,J) * r_QbyOCN(I,J))       &        AREApreTH(I,J) * r_QbyOCN(I,J))
 #endif /* SEAICE_GROWTH_LEGACY */  
1591    
1592            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1593            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
# Line 1907  c ToM>>> Line 1850  c ToM>>>
1850  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1851  C =================================  C =================================
1852    
 #ifndef SEAICE_GROWTH_LEGACY  
1853  C note: in legacy, this process is done at the end  C note: in legacy, this process is done at the end
1854  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1855  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
# Line 1949  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1891  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1891           ENDDO           ENDDO
1892  #endif  #endif
1893          ENDIF          ENDIF
1894  #endif /* SEAICE_GROWTH_LEGACY */  
1895  #ifdef SEAICE_DEBUG  #ifdef SEAICE_DEBUG
1896  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1897          WRITE(msgBuf,msgBufForm)          WRITE(msgBuf,msgBufForm)
# Line 2068  C       i.e. change sin AREA only occur Line 2010  C       i.e. change sin AREA only occur
2010            ELSE            ELSE
2011             recip_HO=1. _d 0 / HO             recip_HO=1. _d 0 / HO
2012            ENDIF            ENDIF
2013  #ifdef SEAICE_GROWTH_LEGACY            recip_HH = recip_heffActual(I,J)
            tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)  
            recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0)  
 #else  
            recip_HH = recip_heffActual(I,J)  
 #endif  
2014    
2015  C gain of ice over open water : computed from  C gain of ice over open water : computed from
2016  C   (SEAICE_areaGainFormula.EQ.1) from growth by ATM  C   (SEAICE_areaGainFormula.EQ.1) from growth by ATM
# Line 2209  CADJ &                       key = iicek Line 2146  CADJ &                       key = iicek
2146    
2147  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
2148    
 #ifdef SEAICE_GROWTH_LEGACY  
 # ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  
 # endif /* ALLOW_AUTODIFF_TAMC */  
         DO J=1,sNy  
          DO I=1,sNx  
 C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean  
           IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN  
              saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *  
      &            HSALT(I,J,bi,bj) * recip_deltaTtherm  
              HSALT(I,J,bi,bj) = 0.0 _d 0  
           ENDIF  
          ENDDO  
         ENDDO  
 #endif /* SEAICE_GROWTH_LEGACY */  
   
2149  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
2150  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2151  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 2271  C update HSALT based on surface saltFlux Line 2192  C update HSALT based on surface saltFlux
2192       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
2193            saltFlux(I,J,bi,bj) =            saltFlux(I,J,bi,bj) =
2194       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
 #ifdef SEAICE_GROWTH_LEGACY  
 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean  
           IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN  
            saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -  
      &          HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm  
            HSALT(I,J,bi,bj) = 0.0 _d 0  
 #ifdef ALLOW_SALT_PLUME  
            saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0  
 #endif /* ALLOW_SALT_PLUME */  
           ENDIF  
 #endif /* SEAICE_GROWTH_LEGACY */  
2195           ENDDO           ENDDO
2196          ENDDO          ENDDO
2197  #endif /* SEAICE_VARIABLE_SALINITY */  #endif /* SEAICE_VARIABLE_SALINITY */
2198    
 C =======================================================================  
 C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===  
 C =======================================================================  
   
 #ifdef SEAICE_GROWTH_LEGACY  
   
 C treat values of ice cover fraction oustide  
 C the [0 1] range, and other such issues.  
 C ===========================================  
   
 Cgf note: this part cannot be heat and water conserving  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE area(:,:,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     NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE  
 CML   replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably  
 CML   meant to be something like a minimum thickness  
           AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4)  
          ENDDO  
         ENDDO  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  
 CADJ &                       key = iicekey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         DO J=1,sNy  
          DO I=1,sNx  
 C NOW TRUNCATE AREA  
           AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))  
          ENDDO  
         ENDDO  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,  
 CADJ &                        key = iicekey, byte = isbyte  
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  
 CADJ &                         key = iicekey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         DO J=1,sNy  
          DO I=1,sNx  
           AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))  
           HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))  
           AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)  
           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)  
 #ifdef SEAICE_CAP_HEFF  
 C     This is not energy conserving, but at least it conserves fresh water  
           tmpscal0         = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)  
           d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0  
           HEFF(I,J,bi,bj)  = HEFF(I,J,bi,bj)  + tmpscal0  
 #endif /* SEAICE_CAP_HEFF */  
           HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)  
          ENDDO  
         ENDDO  
   
 C convert snow to ice if submerged.  
 C =================================  
   
         IF ( SEAICEuseFlooding ) THEN  
          DO J=1,sNy  
           DO I=1,sNx  
            tmpscal0     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow  
      &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst  
            tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))  
            d_HEFFbyFLOODING(I,J)=tmpscal1  
            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)  
            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-  
      &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW  
           ENDDO  
          ENDDO  
         ENDIF  
   
 #endif /* SEAICE_GROWTH_LEGACY */  
   
2199  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
2200          DO J=1,sNy          DO J=1,sNy
2201           DO I=1,sNx           DO I=1,sNx
# Line 2426  CADJ STORE d_hsnwbyocnonsnw = comlev1_bi Line 2257  CADJ STORE d_hsnwbyocnonsnw = comlev1_bi
2257          DO J=1,sNy          DO J=1,sNy
2258           DO I=1,sNx           DO I=1,sNx
2259            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
 #ifndef SEAICE_GROWTH_LEGACY  
 C in principle a_QSWbyATM_cover should always be included here, however  
 C for backward compatibility it is left out of the LEGACY branch  
2260       &         +   a_QSWbyATM_cover(I,J)       &         +   a_QSWbyATM_cover(I,J)
 #endif /* SEAICE_GROWTH_LEGACY */  
2261       &         - ( d_HEFFbyOCNonICE(I,J)       &         - ( d_HEFFbyOCNonICE(I,J)
2262       &           + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE       &           + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2263       &           + d_HEFFbyNEG(I,J)       &           + d_HEFFbyNEG(I,J)
# Line 2536  CML So this diagnostic is only good for Line 2363  CML So this diagnostic is only good for
2363  CML the ice-ocean system.  CML the ice-ocean system.
2364             SIatmQnt(I,J,bi,bj) =             SIatmQnt(I,J,bi,bj) =
2365       &            maskC(I,J,kSurface,bi,bj)*convertHI2Q*(       &            maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
 #ifndef SEAICE_GROWTH_LEGACY  
2366       &            a_QSWbyATM_cover(I,J) +       &            a_QSWbyATM_cover(I,J) +
 #endif /* SEAICE_GROWTH_LEGACY */  
2367       &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )       &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2368  cgf 2) SItflux (analogous to tflux; includes advection by water  cgf 2) SItflux (analogous to tflux; includes advection by water
2369  C             exchanged between atmosphere and ocean+ice)  C             exchanged between atmosphere and ocean+ice)
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22