/[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.5 by torge, Tue Oct 2 16:47:41 2012 UTC revision 1.9 by torge, Mon Oct 22 16:36:45 2012 UTC
# Line 58  C     !FUNCTIONS: Line 58  C     !FUNCTIONS:
58    
59  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
60  C     === Local variables ===  C     === Local variables ===
61    #ifdef SEAICE_DEBUG
62  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
63  C     msgBuf      :: Informational/error message buffer  C     msgBuf      :: Informational/error message buffer
64        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
65  c ToM>>>  c ToM>>>
66    #endif
67  C  C
68  C unit/sign convention:  C unit/sign convention:
69  C    Within the thermodynamic computation all stocks, except HSNOW,  C    Within the thermodynamic computation all stocks, except HSNOW,
# Line 115  C     a_QSWbyATM_open   - short wave hea Line 117  C     a_QSWbyATM_open   - short wave hea
117  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2
118        _RL a_QSWbyATM_open     (1:sNx,1:sNy)        _RL a_QSWbyATM_open     (1:sNx,1:sNy)
119        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)
120  C     a_QbyOCN :: available heat (in in W/m^2) due to the  C     a_QbyOCN :: available heat (in W/m^2) due to the
121  C             interaction of the ice pack and the ocean surface  C             interaction of the ice pack and the ocean surface
122  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting
123  C             processes have been accounted for  C             processes have been accounted for
# Line 177  C                     as EVAP (positive Line 179  C                     as EVAP (positive
179  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
180  C     The latent heat flux which will sublimate all snow and ice  C     The latent heat flux which will sublimate all snow and ice
181  C     over one time step  C     over one time step
 #ifdef SEAICE_ITD  
       _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)  
 #else  
182        _RL latentHeatFluxMax   (1:sNx,1:sNy)        _RL latentHeatFluxMax   (1:sNx,1:sNy)
183  #endif        _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
184  #endif  #endif
185    
186  C     actual ice thickness (with upper and lower limit)  C     actual ice thickness (with upper and lower limit)
# Line 204  C     AREA_PRE :: hold sea-ice fraction Line 203  C     AREA_PRE :: hold sea-ice fraction
203        _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)        _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
204        _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)        _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
205        _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)        _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
206        _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)        _RL leadIceThickMin
207  #endif  #endif
208    
209  C     wind speed  C     wind speed
# Line 281  C ====================================== Line 280  C ======================================
280        ENDIF        ENDIF
281    
282  C     avoid unnecessary divisions in loops  C     avoid unnecessary divisions in loops
283  #ifdef SEAICE_ITD  c#ifdef SEAICE_ITD
284  CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)  CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
285  #endif  C    (see SEAICE_SIZE.h and seaice_readparms.F)
286    c     SEAICE_multDim = nITD
287    c#endif
288        recip_multDim        = SEAICE_multDim        recip_multDim        = SEAICE_multDim
289        recip_multDim        = ONE / recip_multDim        recip_multDim        = ONE / recip_multDim
290  C     above/below: double/single precision calculation of recip_multDim  C     above/below: double/single precision calculation of recip_multDim
# Line 379  C ===================== Line 380  C =====================
380            d_HEFFbySublim(I,J)        = 0.0 _d 0            d_HEFFbySublim(I,J)        = 0.0 _d 0
381            d_HSNWbySublim(I,J)        = 0.0 _d 0            d_HSNWbySublim(I,J)        = 0.0 _d 0
382  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
 #ifdef SEAICE_ITD  
           DO IT=1,SEAICE_multDim  
            latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0  
           ENDDO  
 #else  
383            latentHeatFluxMax(I,J)     = 0.0 _d 0            latentHeatFluxMax(I,J)     = 0.0 _d 0
384  #endif  #endif
 #endif  
385  c  c
386            d_HFRWbyRAIN(I,J)          = 0.0 _d 0            d_HFRWbyRAIN(I,J)          = 0.0 _d 0
387    
# Line 401  c Line 396  c
396              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
397              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
398              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
399    #ifdef SEAICE_CAP_SUBLIM
400                latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
401    #endif
402  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
403              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
404              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
# Line 508  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 506  CADJ STORE area(:,:,bi,bj) = comlev1_bib
506             tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)             tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
507             HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3             HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
508             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
509             AREAITD(I,J,IT,bi,bj)=MAX(-AREAITD(I,J,IT,bi,bj),0. _d 0)             AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
510            ENDDO            ENDDO
511  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
512  C         by calling SEAICE_ITD_SUM  C         by calling SEAICE_ITD_SUM
# Line 668  C    and update AREA, HEFF, and HSNOW Line 666  C    and update AREA, HEFF, and HSNOW
666          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
667          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
668    
669    #ifdef SEAICE_DEBUG
670    c ToM<<< debug seaice_growth
671            WRITE(msgBuf,'(A,7F8.4)')
672         &    ' SEAICE_GROWTH: Heff increments 0, HEFFITD = ',
673         &     HEFFITD(1,1,:,bi,bj)
674            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
675         &    SQUEEZE_RIGHT , myThid)
676            WRITE(msgBuf,'(A,7F8.4)')
677         &    ' SEAICE_GROWTH: Area increments 0, AREAITD = ',
678         &     AREAITD(1,1,:,bi,bj)
679            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
680         &    SQUEEZE_RIGHT , myThid)
681  #endif  #endif
682    #else
683    #ifdef SEAICE_DEBUG
684            WRITE(msgBuf,'(A,7F8.4)')
685         &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
686         &     HEFF(1,1,bi,bj)
687            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
688         &    SQUEEZE_RIGHT , myThid)
689            WRITE(msgBuf,'(A,7F8.4)')
690         &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
691         &     AREA(1,1,bi,bj)
692            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
693         &    SQUEEZE_RIGHT , myThid)
694    c ToM>>>
695    #endif
696    #endif /* SEAICE_ITD */
697  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
698  C        end SEAICEadjMODE.EQ.0 statement:  C        end SEAICEadjMODE.EQ.0 statement:
699          ENDIF          ENDIF
# Line 700  C 3) store regularized values of heff, h Line 725  C 3) store regularized values of heff, h
725             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
726    
727  C memorize areal and volume fraction of each ITD category  C memorize areal and volume fraction of each ITD category
728             IF (AREA(I,J,bi,bj).GT.0) THEN             IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
729              areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)              areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
730             ELSE             ELSE
731              areaFracFactor(I,J,IT)=ZERO  C           if there's no ice, potential growth starts in 1st category
732             ENDIF              IF (IT .EQ. 1) THEN
733             IF (HEFF(I,J,bi,bj).GT.0) THEN               areaFracFactor(I,J,IT)=ONE
734              heffFracFactor(I,J,IT)=HEFFITD(I,J,IT,bi,bj)/HEFF(I,J,bi,bj)              ELSE
735             ELSE               areaFracFactor(I,J,IT)=ZERO
736              heffFracFactor(I,J,IT)=ZERO              ENDIF
737             ENDIF             ENDIF
738            ENDDO            ENDDO
739           ENDDO           ENDDO
# Line 1072  C     update TICE & TICES Line 1097  C     update TICE & TICES
1097  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1098  C     calculate area weighted mean  C     calculate area weighted mean
1099  C     (although the ice's temperature relates to its energy content  C     (although the ice's temperature relates to its energy content
1100  C      and hence should be averaged weighted by ice volume [heffFracFactor],  C      and hence should be averaged weighted by ice volume,
1101  C      the temperature here is a result of the fluxes through the ice surface  C      the temperature here is a result of the fluxes through the ice surface
1102  C      computed individually for each single category in SEAICE_SOLVE4TEMP  C      computed individually for each single category in SEAICE_SOLVE4TEMP
1103  C      and hence is averaged area weighted [areaFracFactor])  C      and hence is averaged area weighted [areaFracFactor])
# Line 1141  C switch heat fluxes from W/m2 to 'effec Line 1166  C switch heat fluxes from W/m2 to 'effec
1166       &          * convertQ2HI * AREAITDpreTH(I,J,IT)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1167             a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)             a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1168       &          * convertQ2HI * AREAITDpreTH(I,J,IT)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1169  C and initialize r_QbyATM_cover  C and initialize r_QbyATMmult_cover
1170             r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)             r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
1171  C     Convert fresh water flux by sublimation to 'effective' ice meters.  C     Convert fresh water flux by sublimation to 'effective' ice meters.
1172  C     Negative sublimation is resublimation and will be added as snow.  C     Negative sublimation is resublimation and will be added as snow.
# Line 1304  C         accumulate change over ITD cat Line 1329  C         accumulate change over ITD cat
1329            HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2            HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
1330       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1331            r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2            r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
 C         keep total up to date, too  
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)        - tmpscal2  
1332  #else  #else
1333       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
1334            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
# Line 1325  C     If anything is left, sublimate ice Line 1348  C     If anything is left, sublimate ice
1348  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1349       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1350  C         accumulate change over ITD categories  C         accumulate change over ITD categories
1351            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2            d_HSNWbySublim(I,J)      = d_HSNWbySublim(I,J)      - tmpscal2
1352            HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2            HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
1353            r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2            r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
 C         keep total up to date, too  
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2  
1354  #else  #else
1355       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1356            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
# Line 1348  C     remove the fusion part for the res Line 1369  C     remove the fusion part for the res
1369       &                               - r_FWbySublimMult(I,J,IT)       &                               - r_FWbySublimMult(I,J,IT)
1370            r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)            r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1371       &                               - r_FWbySublimMult(I,J,IT)       &                               - r_FWbySublimMult(I,J,IT)
1372    #else
1373              a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1374              r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1375    #endif
1376           ENDDO           ENDDO
1377          ENDDO          ENDDO
1378    #ifdef SEAICE_ITD
1379  C       end IT loop  C       end IT loop
1380          ENDDO          ENDDO
 C       then update totals        
         DO J=1,sNy  
          DO I=1,sNx  
1381  #endif  #endif
1382            a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)  #ifdef SEAICE_DEBUG
           r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)  
          ENDDO  
         ENDDO  
1383  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1384          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1385  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1386       &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1387       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1388            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1389         &    SQUEEZE_RIGHT , myThid)
1390            WRITE(msgBuf,'(A,7F8.4)')
1391         &    ' SEAICE_GROWTH: Area increments 1, AREAITD = ',
1392         &     AREAITD(1,1,:,bi,bj)
1393  #else  #else
1394       &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1395       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1396  #endif  #endif
1397          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1398       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1399  c ToM>>>  c ToM>>>
1400    #endif
1401    
1402  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1403  C ===========================================================  C ===========================================================
# Line 1385  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1411  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1411          DO IT=1,nITD          DO IT=1,nITD
1412           DO J=1,sNy           DO J=1,sNy
1413            DO I=1,sNx            DO I=1,sNx
1414  C          ice growth/melt due to ocean heat is equally distributed under the ice  C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
1415  C          and hence weighted by fractional area of each thickness category  C          equally distributed under the ice and hence weighted by
1416    C          fractional area of each thickness category
1417             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1418       &                               -HEFFITD(I,J,IT,bi,bj))       &                               -HEFFITD(I,J,IT,bi,bj))
1419             d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1             d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
1420             r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal1
            HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1  
1421  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1422             SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)             SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1423       &                           + HEFFITD(I,J,IT,bi,bj)       &                           + HEFFITD(I,J,IT,bi,bj)
# Line 1399  C          and hence weighted by fractio Line 1425  C          and hence weighted by fractio
1425            ENDDO            ENDDO
1426           ENDDO           ENDDO
1427          ENDDO          ENDDO
1428            DO J=1,sNy
1429             DO I=1,sNx
1430              r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1431             ENDDO
1432            ENDDO
1433  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
1434          DO J=1,sNy          DO J=1,sNy
1435           DO I=1,sNx           DO I=1,sNx
# Line 1411  C          and hence weighted by fractio Line 1442  C          and hence weighted by fractio
1442           ENDDO           ENDDO
1443          ENDDO          ENDDO
1444  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1445    #ifdef SEAICE_DEBUG
1446  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1447          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1448  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1449       &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1450       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1451            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1452         &    SQUEEZE_RIGHT , myThid)
1453            WRITE(msgBuf,'(A,7F8.4)')
1454         &    ' SEAICE_GROWTH: Area increments 2, AREAITD = ',
1455         &     AREAITD(1,1,:,bi,bj)
1456  #else  #else
1457       &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1458       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1459  #endif  #endif
1460          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1461       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1462  c ToM>>>  c ToM>>>
1463    #endif
1464    
1465  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1466  C ==================================================================  C ==================================================================
# Line 1451  Cgf no additional dependency through sno Line 1489  Cgf no additional dependency through sno
1489       &                           + tmpscal2*ICE2SNOW       &                           + tmpscal2*ICE2SNOW
1490             r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)             r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1491       &                           - tmpscal2       &                           - tmpscal2
 C          keep the total up to date, too  
            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2  
1492            ENDDO            ENDDO
1493           ENDDO           ENDDO
1494          ENDDO          ENDDO
# Line 1473  Cgf no additional dependency through sno Line 1509  Cgf no additional dependency through sno
1509           ENDDO           ENDDO
1510          ENDDO          ENDDO
1511  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1512    #ifdef SEAICE_DEBUG
1513  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1514          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1515  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1516       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1517       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1518            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1519         &    SQUEEZE_RIGHT , myThid)
1520            WRITE(msgBuf,'(A,7F8.4)')
1521         &    ' SEAICE_GROWTH: Area increments 3, AREAITD = ',
1522         &     AREAITD(1,1,:,bi,bj)
1523  #else  #else
1524       &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1525       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1526  #endif  #endif
1527          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1528       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1529  c ToM>>>  c ToM>>>
1530    #endif
1531    
1532  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1533  C ====================================================  C ====================================================
# Line 1516  c         Limit ice growth by potential Line 1559  c         Limit ice growth by potential
1559       &                                 + tmpscal2       &                                 + tmpscal2
1560             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1561       &                                 + tmpscal2       &                                 + tmpscal2
1562             r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)             r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
1563       &                                 - tmpscal2       &                                 - tmpscal2
1564             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
1565    
# Line 1550  c         Limit ice growth by potential Line 1593  c         Limit ice growth by potential
1593           ENDDO           ENDDO
1594          ENDDO          ENDDO
1595  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1596    #ifdef SEAICE_DEBUG
1597  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1598          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1599  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1600       &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1601       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1602            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1603         &    SQUEEZE_RIGHT , myThid)
1604            WRITE(msgBuf,'(A,7F8.4)')
1605         &    ' SEAICE_GROWTH: Area increments 4, AREAITD = ',
1606         &     AREAITD(1,1,:,bi,bj)
1607  #else  #else
1608       &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1609       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1610  #endif  #endif
1611          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1612       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1613  c ToM>>>  c ToM>>>
1614    #endif
1615    
1616  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1617  C depending on atmospheric conditions.  C depending on atmospheric conditions.
# Line 1610  Cgf note: this does not affect air-sea h Line 1660  Cgf note: this does not affect air-sea h
1660  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1661  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1662  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1663    #ifdef SEAICE_DEBUG
1664  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1665          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1666  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1667       &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1668       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1669            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1670         &    SQUEEZE_RIGHT , myThid)
1671            WRITE(msgBuf,'(A,7F8.4)')
1672         &    ' SEAICE_GROWTH: Area increments 5, AREAITD = ',
1673         &     AREAITD(1,1,:,bi,bj)
1674  #else  #else
1675       &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1676       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1677  #endif  #endif
1678          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1679       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1680  c ToM>>>  c ToM>>>
1681    #endif
1682    
1683  C compute snow melt due to heat available from ocean.  C compute snow melt due to heat available from ocean.
1684  C =================================================================  C =================================================================
# Line 1669  Cgf no additional dependency through sno Line 1726  Cgf no additional dependency through sno
1726  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1727  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1728  Cph)  Cph)
1729    #ifdef SEAICE_DEBUG
1730  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1731          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1732  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1733       &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1734       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1735            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1736         &    SQUEEZE_RIGHT , myThid)
1737            WRITE(msgBuf,'(A,7F8.4)')
1738         &    ' SEAICE_GROWTH: Area increments 6, AREAITD = ',
1739         &     AREAITD(1,1,:,bi,bj)
1740  #else  #else
1741       &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1742       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1743  #endif  #endif
1744          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1745       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1746  c ToM>>>  c ToM>>>
1747    #endif
1748    
1749  C gain of new ice over open water  C gain of new ice over open water
1750  C ===============================  C ===============================
# Line 1712  C           or 0. otherwise (no melting Line 1776  C           or 0. otherwise (no melting
1776  C         open water area fraction  C         open water area fraction
1777            tmpscal0 = ONE-AREApreTH(I,J)            tmpscal0 = ONE-AREApreTH(I,J)
1778  C         determine thickness of new ice  C         determine thickness of new ice
1779  C         considering the entire open water area to refreeze  ctomC         considering the entire open water area to refreeze
1780            tmpscal1 = tmpscal3/tmpscal0  ctom          tmpscal1 = tmpscal3/tmpscal0
1781    C         considering a minimum lead ice thickness of 10 cm
1782    C         WATCH that leadIceThickMin is smaller that Hlimit(1)!
1783              leadIceThickMin = 0.1
1784              tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
1785    C         adjust ice area fraction covered by new ice
1786              tmpscal0 = tmpscal3/tmpscal1
1787  C         then add new ice volume to appropriate thickness category  C         then add new ice volume to appropriate thickness category
1788            DO IT=1,nITD            DO IT=1,nITD
1789             IF (tmpscal1.LT.Hlimit(IT)) THEN             IF (tmpscal1.LT.Hlimit(IT)) THEN
# Line 1755  c needs to be here to allow use also wit Line 1825  c needs to be here to allow use also wit
1825          ENDDO          ENDDO
1826  #endif  #endif
1827  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1828    #ifdef SEAICE_DEBUG
1829  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1830          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1831  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1832       &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1833       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1834            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1835         &    SQUEEZE_RIGHT , myThid)
1836            WRITE(msgBuf,'(A,7F8.4)')
1837         &    ' SEAICE_GROWTH: Area increments 7, AREAITD = ',
1838         &     AREAITD(1,1,:,bi,bj)
1839  #else  #else
1840       &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1841       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1842            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1843         &    SQUEEZE_RIGHT , myThid)
1844            WRITE(msgBuf,'(A,7F8.4)')
1845         &    ' SEAICE_GROWTH: Area increments 7, AREA = ',
1846         &     AREA(1,1,bi,bj)
1847  #endif  #endif
1848          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1849       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1850  c ToM>>>  c ToM>>>
1851    #endif
1852    
1853  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1854  C =================================  C =================================
# Line 1808  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1890  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1890  #endif  #endif
1891          ENDIF          ENDIF
1892  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1893    #ifdef SEAICE_DEBUG
1894  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1895          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,'(A,7F8.4)')
1896  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1897       &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1898       &     HEFFITD(20,20,:,bi,bj)       &     HEFFITD(1,1,:,bi,bj)
1899            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1900         &    SQUEEZE_RIGHT , myThid)
1901            WRITE(msgBuf,'(A,7F8.4)')
1902         &    ' SEAICE_GROWTH: Area increments 8, AREAITD = ',
1903         &     AREAITD(1,1,:,bi,bj)
1904  #else  #else
1905       &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',       &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1906       &     HEFF(20,20,bi,bj)       &     HEFF(1,1,bi,bj)
1907            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1908         &    SQUEEZE_RIGHT , myThid)
1909            WRITE(msgBuf,'(A,7F8.4)')
1910         &    ' SEAICE_GROWTH: Area increments 8, AREA = ',
1911         &     AREA(1,1,bi,bj)
1912  #endif  #endif
1913          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1914       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1915  c ToM>>>  c ToM>>>
1916    #endif
1917    
1918  C ===================================================================  C ===================================================================
1919  C ==========PART 4: determine ice cover fraction increments=========-  C ==========PART 4: determine ice cover fraction increments=========-
# Line 1850  C       replaces Hibler '79 scheme and l Line 1944  C       replaces Hibler '79 scheme and l
1944  C       because ITD accounts explicitly for lead openings and  C       because ITD accounts explicitly for lead openings and
1945  C       different melt rates due to varying ice thickness  C       different melt rates due to varying ice thickness
1946  C  C
1947  C       only consider ice area loss due to total ice thickness loss  C       only consider ice area loss due to total ice thickness loss;
1948  C       ice area gain due to freezing of open water as handled above  C       ice area gain due to freezing of open water is handled above
1949  C       under "gain of new ice over open water"  C       under "gain of new ice over open water"
1950  C  C
1951  C       does not account for lateral melt of ice floes  C       does not account for lateral melt of ice floes
# Line 1860  C        AREAITD is incremented in secti Line 1954  C        AREAITD is incremented in secti
1954  C  C
1955          DO IT=1,nITD          DO IT=1,nITD
1956           DO J=1,sNy           DO J=1,sNy
1957            DO I=1,sNx           DO I=1,sNx
1958             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1959              AREAITD(I,J,IT,bi,bj)=ZERO              AREAITD(I,J,IT,bi,bj)=ZERO
1960             ENDIF             ENDIF
# Line 2186  C compute net heat flux leaving/entering Line 2280  C compute net heat flux leaving/entering
2280  C accounting for the part used in melt/freeze processes  C accounting for the part used in melt/freeze processes
2281  C =====================================================  C =====================================================
2282    
2283    #ifdef SEAICE_ITD
2284    C compute total of "mult" fluxes for ocean forcing
2285            DO J=1,sNy
2286             DO I=1,sNx
2287              a_QbyATM_cover(I,J)   = 0.0 _d 0
2288              r_QbyATM_cover(I,J)   = 0.0 _d 0
2289              a_QSWbyATM_cover(I,J) = 0.0 _d 0
2290              r_FWbySublim(I,J)     = 0.0 _d 0
2291             ENDDO
2292            ENDDO
2293            DO IT=1,nITD
2294             DO J=1,sNy
2295              DO I=1,sNx
2296    cToM if fluxes in W/m^2 then
2297    c           a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2298    c     &      + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2299    c           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2300    c     &      + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2301    c           a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2302    c     &      + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2303    c           r_FWbySublim(I,J)=r_FWbySublim(I,J)
2304    c     &      + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
2305    cToM if fluxes in effective ice meters, i.e. ice volume per area, then
2306               a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2307         &      + a_QbyATMmult_cover(I,J,IT)
2308               r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2309         &      + r_QbyATMmult_cover(I,J,IT)
2310               a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2311         &      + a_QSWbyATMmult_cover(I,J,IT)
2312               r_FWbySublim(I,J)=r_FWbySublim(I,J)
2313         &      + r_FWbySublimMult(I,J,IT)
2314              ENDDO
2315             ENDDO
2316            ENDDO
2317    #endif
2318    
2319  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
2320  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
2321  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22