/[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.4 by torge, Fri Sep 28 19:01:17 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 141  c     The change of mean ice thickness d Line 143  c     The change of mean ice thickness d
143        _RL d_HEFFbyRLX         (1:sNx,1:sNy)        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
144  #endif  #endif
145    
 #ifdef SEAICE_ITD  
 c     The change of mean ice area due to out-of-bounds values following  
 c     sea ice dynamics  
       _RL d_AREAbyNEG         (1:sNx,1:sNy)  
 #endif  
146  c     The change of mean ice thickness due to out-of-bounds values following  c     The change of mean ice thickness due to out-of-bounds values following
147  c     sea ice dynamics  c     sea ice dynamics
148        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
# Line 179  C                     as EVAP (positive Line 176  C                     as EVAP (positive
176        _RL d_HEFFbySublim      (1:sNx,1:sNy)        _RL d_HEFFbySublim      (1:sNx,1:sNy)
177        _RL d_HSNWbySublim      (1:sNx,1:sNy)        _RL d_HSNWbySublim      (1:sNx,1:sNy)
178    
179  #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))  #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
182        _RL latentHeatFluxMax   (1:sNx,1:sNy)        _RL latentHeatFluxMax   (1:sNx,1:sNy)
# Line 206  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 233  C temporary variables available for the Line 230  C temporary variables available for the
230    
231        INTEGER ilockey        INTEGER ilockey
232        INTEGER it        INTEGER it
 #ifdef SEAICE_ITD  
       INTEGER K  
 #endif  
233        _RL pFac        _RL pFac
234        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
235        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
# Line 286  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 367  C ===================== Line 363  C =====================
363            d_HEFFbyRLX(I,J)           = 0.0 _d 0            d_HEFFbyRLX(I,J)           = 0.0 _d 0
364  #endif  #endif
365    
 #ifdef SEAICE_ITD  
           d_AREAbyNEG(I,J)           = 0.0 _d 0  
 #endif  
366            d_HEFFbyNEG(I,J)           = 0.0 _d 0            d_HEFFbyNEG(I,J)           = 0.0 _d 0
367            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
368            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
# Line 403  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
405  #endif  #endif
 #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))  
             latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0  
 #endif  
406            ENDDO            ENDDO
407           ENDDO           ENDDO
408          ENDDO          ENDDO
# Line 444  C ====================================== Line 437  C ======================================
437           ENDDO           ENDDO
438          ENDDO          ENDDO
439  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
440          DO K=1,nITD          DO IT=1,nITD
441           DO J=1,sNy           DO J=1,sNy
442            DO I=1,sNx            DO I=1,sNx
443             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
444             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
445             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
446            ENDDO            ENDDO
447           ENDDO           ENDDO
448          ENDDO          ENDDO
# Line 504  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 497  CADJ STORE area(:,:,bi,bj) = comlev1_bib
497          DO J=1,sNy          DO J=1,sNy
498           DO I=1,sNx           DO I=1,sNx
499  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
500            DO K=1,nITD            DO IT=1,nITD
            tmpscal1=0. _d 0  
501             tmpscal2=0. _d 0             tmpscal2=0. _d 0
502             tmpscal3=0. _d 0             tmpscal3=0. _d 0
503             tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)             tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
504             HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2             HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
505             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
506             tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)             tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
507             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,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             tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)             AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1  
            d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1  
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
513  #else  #else
514            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
           d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)  
           AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)  
515            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
516              d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
517            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
518              AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
519  #endif  #endif
520           ENDDO           ENDDO
521          ENDDO          ENDDO
# Line 538  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 528  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
528          DO J=1,sNy          DO J=1,sNy
529           DO I=1,sNx           DO I=1,sNx
530  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
531            DO K=1,nITD            DO IT=1,nITD
532  #endif  #endif
533            tmpscal2=0. _d 0            tmpscal2=0. _d 0
534            tmpscal3=0. _d 0            tmpscal3=0. _d 0
535  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
536             IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN             IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
537              tmpscal2=-HEFFITD(I,J,K,bi,bj)              tmpscal2=-HEFFITD(I,J,IT,bi,bj)
538              tmpscal3=-HSNOWITD(I,J,K,bi,bj)              tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
539              TICES(I,J,K,bi,bj)=celsius2K              TICES(I,J,IT,bi,bj)=celsius2K
540  CToM        TICE will be updated at end of Part 1 together with AREA and HEFF  CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
541             ENDIF             ENDIF
542             HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2             HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
543             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3             HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
544  #else  #else
545            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
546             tmpscal2=-HEFF(I,J,bi,bj)             tmpscal2=-HEFF(I,J,bi,bj)
# Line 580  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 570  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
570          DO J=1,sNy          DO J=1,sNy
571           DO I=1,sNx           DO I=1,sNx
572  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
573            DO K=1,nITD            DO IT=1,nITD
574             IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.             IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
575       &         (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))       &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
576       &      AREAITD(I,J,K,bi,bj)=0. _d 0       &      AREAITD(I,J,IT,bi,bj)=0. _d 0
577            ENDDO            ENDDO
578  #else  #else
579            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
# Line 601  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 591  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
591          DO J=1,sNy          DO J=1,sNy
592           DO I=1,sNx           DO I=1,sNx
593  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
594            DO K=1,nITD            DO IT=1,nITD
595             IF ((HEFFITD(i,j,k,bi,bj).GT.0).OR.             IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
596       &         (HSNOWITD(i,j,k,bi,bj).GT.0)) THEN       &         (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
597  CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1  CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1
598  C           hence use SEAICE_area_floor devided by nITD  C           hence use SEAICE_area_floor devided by nITD
599  C           (or install a warning in e.g. seaice_readparms.F)  C           (or install a warning in e.g. seaice_readparms.F)
600              AREAITD(I,J,K,bi,bj)=              AREAITD(I,J,IT,bi,bj)=
601       &         MAX(AREAITD(I,J,K,bi,bj),SEAICE_area_floor/float(nITD))       &         MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
602             ENDIF             ENDIF
603            ENDDO            ENDDO
604  #else  #else
# Line 639  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 629  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
629            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
630           ENDDO           ENDDO
631          ENDDO          ENDDO
632  #endif /* SEAICE_ITD */  #endif /* notSEAICE_ITD */
633    
634  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
635  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
 C    first, update AREA and HEFF:  
         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)  
 C  
636          DO J=1,sNy          DO J=1,sNy
637           DO I=1,sNx           DO I=1,sNx
638  C    TICES was changed above (item 1.25), now update TICE as ice volume  C    TICES was changed above (item 1.25), now update TICE as ice volume
639  C     weighted average of TICES  C     weighted average of TICES
640    C    also compute total of AREAITD (needed for finishing item 2.5, see below)
641            tmpscal1 = 0. _d 0            tmpscal1 = 0. _d 0
642            tmpscal2 = 0. _d 0            tmpscal2 = 0. _d 0
643            DO K=1,nITD            tmpscal3 = 0. _d 0
644             tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)            DO IT=1,nITD
645             tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)             tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
646               tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
647               tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
648            ENDDO            ENDDO
649            TICE(I,J,bi,bj)=tmpscal1/tmpscal2            TICE(I,J,bi,bj)=tmpscal1/tmpscal2
650  C    lines of item 2.5 that were omitted:  C    lines of item 2.5 that were omitted:
# Line 662  C    in 2.5 these lines are executed bef Line 652  C    in 2.5 these lines are executed bef
652  C    hence we execute them here before SEAICE_ITD_REDIST is called  C    hence we execute them here before SEAICE_ITD_REDIST is called
653  C    although this means that AREA has not been completely regularized  C    although this means that AREA has not been completely regularized
654  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
655            DIAGarrayA(I,J) = AREA(I,J,bi,bj)            DIAGarrayA(I,J) = tmpscal3
656  #endif  #endif
657  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
658            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)            SItrAREA(I,J,bi,bj,1)=tmpscal3
659  #endif  #endif
660           ENDDO           ENDDO
661          ENDDO          ENDDO
# Line 676  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
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  #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        ENDIF SEAICEadjMODE.EQ.0  C        end SEAICEadjMODE.EQ.0 statement:
699          ENDIF          ENDIF
700  #endif  #endif
701    
# Line 700  C 3) store regularized values of heff, h Line 717  C 3) store regularized values of heff, h
717           ENDDO           ENDDO
718          ENDDO          ENDDO
719  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
720          DO K=1,nITD          DO IT=1,nITD
721           DO J=1,sNy           DO J=1,sNy
722            DO I=1,sNx            DO I=1,sNx
723             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
724             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
725             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,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,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)              areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
            ELSE  
             areaFracFactor(I,J,K)=ZERO  
            ENDIF  
            IF (HEFF(I,J,bi,bj).GT.0) THEN  
             heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)  
730             ELSE             ELSE
731              heffFracFactor(I,J,K)=ZERO  C           if there's no ice, potential growth starts in 1st category
732                IF (IT .EQ. 1) THEN
733                 areaFracFactor(I,J,IT)=ONE
734                ELSE
735                 areaFracFactor(I,J,IT)=ZERO
736                ENDIF
737             ENDIF             ENDIF
738            ENDDO            ENDDO
739           ENDDO           ENDDO
740          ENDDO          ENDDO
741  C prepare SItrHEFF to be computed as cumulative sum  C prepare SItrHEFF to be computed as cumulative sum
742          DO K=2,5          DO iTr=2,5
743           DO J=1,sNy           DO J=1,sNy
744            DO I=1,sNx            DO I=1,sNx
745             SItrHEFF(I,J,bi,bj,K)=ZERO             SItrHEFF(I,J,bi,bj,iTr)=ZERO
746            ENDDO            ENDDO
747           ENDDO           ENDDO
748          ENDDO          ENDDO
# Line 791  Cgf no additional dependency of air-sea Line 808  Cgf no additional dependency of air-sea
808            ENDDO            ENDDO
809           ENDDO           ENDDO
810  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
811           DO K=1,nITD           DO IT=1,nITD
812            DO J=1,sNy            DO J=1,sNy
813             DO I=1,sNx             DO I=1,sNx
814              HEFFITDpreTH(I,J,K) = 0. _d 0              HEFFITDpreTH(I,J,IT) = 0. _d 0
815              HSNWITDpreTH(I,J,K) = 0. _d 0              HSNWITDpreTH(I,J,IT) = 0. _d 0
816              AREAITDpreTH(I,J,K) = 0. _d 0              AREAITDpreTH(I,J,IT) = 0. _d 0
817             ENDDO             ENDDO
818            ENDDO            ENDDO
819           ENDDO           ENDDO
# Line 821  CADJ STORE HEFFpreTH = comlev1_bibj, key Line 838  CADJ STORE HEFFpreTH = comlev1_bibj, key
838  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
839  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
840  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
841          DO K=1,nITD          DO IT=1,nITD
842  #endif  #endif
843          DO J=1,sNy          DO J=1,sNy
844           DO I=1,sNx           DO I=1,sNx
845    
846  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
847            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
848  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
849              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
850       &                              AREAITDpreTH(I,J,K))       &                              AREAITDpreTH(I,J,IT))
851              hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1              hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
852              tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1              tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
853              heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)              heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
854  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
855  cif        regularize AREA with SEAICE_area_reg  cif        regularize AREA with SEAICE_area_reg
856             tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
857       &                     + area_reg_sq)       &                     + area_reg_sq)
858  cif        heffActual calculated with the regularized AREA  cif        heffActual calculated with the regularized AREA
859             tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1             tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
860  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
861             heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2             heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
862       &                                  + hice_reg_sq)       &                                  + hice_reg_sq)
863  cif        hsnowActual calculated with the regularized AREA  cif        hsnowActual calculated with the regularized AREA
864             hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
865  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
866  cif        regularize the inverse of heffActual by hice_reg  cif        regularize the inverse of heffActual by hice_reg
867             recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /             recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
868       &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
869       &                      + hice_reg_sq)       &                      + hice_reg_sq)
870  cif       Do not regularize when HEFFpreTH = 0  cif       Do not regularize when HEFFpreTH = 0
871            ELSE            ELSE
872              heffActualMult(I,J,K) = ZERO              heffActualMult(I,J,IT) = ZERO
873              hsnowActualMult(I,J,K) = ZERO              hsnowActualMult(I,J,IT) = ZERO
874              recip_heffActualMult(I,J,K)  = ZERO              recip_heffActualMult(I,J,IT)  = ZERO
875            ENDIF            ENDIF
876  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
877            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
# Line 900  cif       Do not regularize when HEFFpre Line 917  cif       Do not regularize when HEFFpre
917  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
918  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
919  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
920          DO K=1,nITD          DO IT=1,nITD
921  #endif  #endif
922          DO J=1,sNy          DO J=1,sNy
923           DO I=1,sNx           DO I=1,sNx
# Line 908  c        The latent heat flux over the s Line 925  c        The latent heat flux over the s
925  c        will sublimate all of the snow and ice over one time  c        will sublimate all of the snow and ice over one time
926  c        step (W/m^2)  c        step (W/m^2)
927  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
928            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
929              latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *              latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
930       &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +       &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
931       &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)       &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
932         &         /AREAITDpreTH(I,J,IT)
933            ELSE            ELSE
934              latentHeatFluxMaxMult(I,J,K) = ZERO              latentHeatFluxMaxMult(I,J,IT) = ZERO
935            ENDIF            ENDIF
936  #else  #else
937            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
# Line 1079  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])
1104             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1105       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1106  #else  #else
1107             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1108       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
# Line 1095  C     average over categories Line 1113  C     average over categories
1113  C     calculate area weighted mean  C     calculate area weighted mean
1114  C     (fluxes are per unit (ice surface) area and are thus area weighted)  C     (fluxes are per unit (ice surface) area and are thus area weighted)
1115             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1116       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1117             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1118       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1119             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1120       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1121  #else  #else
1122             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1123       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
# Line 1141  CADJ STORE a_FWbySublim    = comlev1_bib Line 1159  CADJ STORE a_FWbySublim    = comlev1_bib
1159    
1160  C switch heat fluxes from W/m2 to 'effective' ice meters  C switch heat fluxes from W/m2 to 'effective' ice meters
1161  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1162          DO K=1,nITD          DO IT=1,nITD
1163           DO J=1,sNy           DO J=1,sNy
1164            DO I=1,sNx            DO I=1,sNx
1165             a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)             a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
1166       &          * convertQ2HI * AREAITDpreTH(I,J,K)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1167             a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)             a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1168       &          * convertQ2HI * AREAITDpreTH(I,J,K)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1169  C and initialize r_QbyATM_cover  C and initialize r_QbyATMmult_cover
1170             r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)             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.
1173  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1174             a_FWbySublimMult(I,J,K) = ZERO             a_FWbySublimMult(I,J,IT) = ZERO
1175  #endif  #endif
1176             a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce             a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
1177       &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)       &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
1178             r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)             r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
1179            ENDDO            ENDDO
1180           ENDDO           ENDDO
1181          ENDDO          ENDDO
# Line 1214  CADJ STORE r_FWbySublim    = comlev1_bib Line 1232  CADJ STORE r_FWbySublim    = comlev1_bib
1232  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1233          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1234  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1235           DO K=1,nITD           DO IT=1,nITD
1236            DO J=1,sNy            DO J=1,sNy
1237             DO I=1,sNx             DO I=1,sNx
1238              a_QbyATMmult_cover(I,J,K)   = 0. _d 0              a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1239              r_QbyATMmult_cover(I,J,K)   = 0. _d 0              r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1240              a_QSWbyATMmult_cover(I,J,K) = 0. _d 0              a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
1241             ENDDO             ENDDO
1242            ENDDO            ENDDO
1243           ENDDO           ENDDO
# Line 1296  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1314  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1314  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1315  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1316  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1317          DO K=1,nITD          DO IT=1,nITD
1318  #endif  #endif
1319          DO J=1,sNy          DO J=1,sNy
1320           DO I=1,sNx           DO I=1,sNx
1321  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1322            tmpscal2 =            tmpscal2 =
1323  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1324       &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1325       &             *SNOW2ICE),ZERO)       &             *SNOW2ICE),ZERO)
1326  C         accumulate change over ITD categories  C         accumulate change over ITD categories
1327            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
1328       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1329            HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2            HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
1330       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1331            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - 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 1330  CADJ STORE r_FWbySublim    = comlev1_bib Line 1346  CADJ STORE r_FWbySublim    = comlev1_bib
1346  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1347            tmpscal2 =            tmpscal2 =
1348  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1349       &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,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,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2            HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
1353            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - 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 1351  C     If anything is left, it will be ev Line 1365  C     If anything is left, it will be ev
1365  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
1366  C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).  C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
1367  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1368            a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)            a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1369       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1370            r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)            r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1371       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1372           ENDDO  #else
         ENDDO  
 C       end K loop  
         ENDDO  
 C       then update totals        
         DO J=1,sNy  
          DO I=1,sNx  
 #endif  
1373            a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)            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)            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
1380            ENDDO
1381    #endif
1382    #ifdef SEAICE_DEBUG
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 1389  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1408  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1408  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1409    
1410  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1411          DO K=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             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),  C          fractional area of each thickness category
1417       &                               -HEFFITD(I,J,K,bi,bj))             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1418             d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1       &                               -HEFFITD(I,J,IT,bi,bj))
1419             r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1             d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
1420             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,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,K,bi,bj)       &                           + HEFFITD(I,J,IT,bi,bj)
1424  #endif  #endif
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 1418  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 1440  CADJ STORE r_QbyATM_cover = comlev1_bibj Line 1471  CADJ STORE r_QbyATM_cover = comlev1_bibj
1471  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1472    
1473  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1474          DO K=1,nITD          DO IT=1,nITD
1475           DO J=1,sNy           DO J=1,sNy
1476            DO I=1,sNx            DO I=1,sNx
1477  C     Convert to standard units (meters of ice) rather than to meters  C     Convert to standard units (meters of ice) rather than to meters
1478  C     of snow. This appears to be more robust.  C     of snow. This appears to be more robust.
1479            tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)             tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1480       &                                           *SNOW2ICE)       &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1481            tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1482  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1483  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1484            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1485  #endif  #endif
1486            d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW             d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1487            HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW       &                           + tmpscal2*ICE2SNOW
1488            r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2             HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
1489  C         keep the total up to date, too       &                           + tmpscal2*ICE2SNOW
1490            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2             r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1491         &                           - tmpscal2
1492            ENDDO            ENDDO
1493           ENDDO           ENDDO
1494          ENDDO          ENDDO
# Line 1477  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 1504  Cgf the v1.81=>v1.82 revision would chan Line 1543  Cgf the v1.81=>v1.82 revision would chan
1543  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1544    
1545  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1546          DO K=1,nITD          DO IT=1,nITD
1547           DO J=1,sNy           DO J=1,sNy
1548            DO I=1,sNx            DO I=1,sNx
1549  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
1550             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))             tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1551         &                     r_QbyATMmult_cover(I,J,IT))
1552  #else  #else
1553             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)             tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1554         &                     r_QbyATMmult_cover(I,J,IT)
1555  c         Limit ice growth by potential melt by ocean  c         Limit ice growth by potential melt by ocean
1556       &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1557  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1558             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
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,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
1565    
1566  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1567             SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)             SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1568       &                           + HEFFITD(I,J,K,bi,bj)       &                           + HEFFITD(I,J,IT,bi,bj)
1569  #endif  #endif
1570            ENDDO            ENDDO
1571           ENDDO           ENDDO
# Line 1552  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 1593  C           add precip to the fresh wate Line 1641  C           add precip to the fresh wate
1641           ENDDO           ENDDO
1642          ENDDO          ENDDO
1643  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1644          DO K=1,nITD          DO IT=1,nITD
1645           DO J=1,sNy           DO J=1,sNy
1646            DO I=1,sNx            DO I=1,sNx
1647             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)             HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1648       &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)       &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1649            ENDDO            ENDDO
1650           ENDDO           ENDDO
1651          ENDDO          ENDDO
# Line 1612  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 1637  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1692  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1692  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1693    
1694  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1695          DO K=1,nITD          DO IT=1,nITD
1696           DO J=1,sNy           DO J=1,sNy
1697            DO I=1,sNx            DO I=1,sNx
1698             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1699       &                  -HSNOWITD(I,J,K,bi,bj))       &                  -HSNOWITD(I,J,IT,bi,bj))
1700             tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1701  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1702  Cgf no additional dependency through snow  Cgf no additional dependency through snow
# Line 1649  Cgf no additional dependency through sno Line 1704  Cgf no additional dependency through sno
1704  #endif  #endif
1705             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1706             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1707             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2             HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) + tmpscal2
1708            ENDDO            ENDDO
1709           ENDDO           ENDDO
1710          ENDDO          ENDDO
# Line 1671  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 1714  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 K=1,nITD            DO IT=1,nITD
1789             IF (tmpscal1.LT.Hlimit(K)) THEN             IF (tmpscal1.LT.Hlimit(IT)) THEN
1790              HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3              HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
1791              tmpscal3=ZERO              tmpscal3=ZERO
1792  C not sure if AREAITD should be changed here since AREA is incremented  C not sure if AREAITD should be changed here since AREA is incremented
1793  C   in PART 4 below in non-itd code  C   in PART 4 below in non-itd code
1794  C in this scenario all open water is covered by new ice instantaneously,  C in this scenario all open water is covered by new ice instantaneously,
1795  C   i.e. no delayed lead closing is concidered such as is associated with  C   i.e. no delayed lead closing is concidered such as is associated with
1796  C   Hibler's h_0 parameter  C   Hibler's h_0 parameter
1797              AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)              AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
1798       &                           + tmpscal0       &                           + tmpscal0
1799              tmpscal0=ZERO              tmpscal0=ZERO
1800             ENDIF             ENDIF
# Line 1739  C   Hibler's h_0 parameter Line 1807  C   Hibler's h_0 parameter
1807    
1808  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1809  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1810          DO K=1,nITD          DO IT=1,nITD
1811           DO J=1,sNy           DO J=1,sNy
1812            DO I=1,sNx            DO I=1,sNx
1813  c needs to be here to allow use also with LEGACY branch  c needs to be here to allow use also with LEGACY branch
1814             SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)             SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1815       &                           + HEFFITD(I,J,K,bi,bj)       &                           + HEFFITD(I,J,IT,bi,bj)
1816            ENDDO            ENDDO
1817           ENDDO           ENDDO
1818          ENDDO          ENDDO
# Line 1757  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 1781  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1861  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1861  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1862          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1863  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1864           DO K=1,nITD           DO IT=1,nITD
1865            DO J=1,sNy            DO J=1,sNy
1866             DO I=1,sNx             DO I=1,sNx
1867              tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow              tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
1868       &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
1869              tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))       &                                        *recip_rhoConst
1870              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1              tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
1871              HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
1872              HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1              HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1
1873                HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
1874       &                            * ICE2SNOW       &                            * ICE2SNOW
1875             ENDDO             ENDDO
1876            ENDDO            ENDDO
# Line 1809  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 1851  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
1952  C  C
1953  C        AREAITD is incremented in section "gain of new ice over open water" above  C        AREAITD is incremented in section "gain of new ice over open water" above
1954  C  C
1955          DO K=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,K,bi,bj).LE.ZERO) THEN             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1959              AREAITD(I,J,K,bi,bj)=ZERO              AREAITD(I,J,IT,bi,bj)=ZERO
1960             ENDIF             ENDIF
1961  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1962             SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)             SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1963       &                           + AREAITD(I,J,K,bi,bj)       &                           + AREAITD(I,J,IT,bi,bj)
1964  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1965            ENDDO            ENDDO
1966           ENDDO           ENDDO
# Line 1948  C apply tendency Line 2041  C apply tendency
2041  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
2042          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
2043  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
2044           DO K=1,nITD           DO IT=1,nITD
2045            DO J=1,sNy            DO J=1,sNy
2046             DO I=1,sNx             DO I=1,sNx
2047              AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *              AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
2048       &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )       &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
2049             ENDDO             ENDDO
2050            ENDDO            ENDDO
2051           ENDDO           ENDDO
# Line 2187  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.4  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22