/[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.1 by dimitri, Fri Apr 27 22:22:17 2012 UTC revision 1.4 by torge, Fri Sep 28 19:01:17 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    c ToM<<< debug seaice_growth
62    C     msgBuf      :: Informational/error message buffer
63          CHARACTER*(MAX_LEN_MBUF) msgBuf
64    c ToM>>>
65  C  C
66  C unit/sign convention:  C unit/sign convention:
67  C    Within the thermodynamic computation all stocks, except HSNOW,  C    Within the thermodynamic computation all stocks, except HSNOW,
# Line 137  c     The change of mean ice thickness d Line 141  c     The change of mean ice thickness d
141        _RL d_HEFFbyRLX         (1:sNx,1:sNy)        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
142  #endif  #endif
143    
144    #ifdef SEAICE_ITD
145    c     The change of mean ice area due to out-of-bounds values following
146    c     sea ice dynamics
147          _RL d_AREAbyNEG         (1:sNx,1:sNy)
148    #endif
149  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
150  c     sea ice dyhnamics  c     sea ice dynamics
151        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
152    
153  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
# Line 170  C                     as EVAP (positive Line 179  C                     as EVAP (positive
179        _RL d_HEFFbySublim      (1:sNx,1:sNy)        _RL d_HEFFbySublim      (1:sNx,1:sNy)
180        _RL d_HSNWbySublim      (1:sNx,1:sNy)        _RL d_HSNWbySublim      (1:sNx,1:sNy)
181    
182  #ifdef SEAICE_CAP_SUBLIM  #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
183  C     The latent heat flux which will sublimate all snow and ice  C     The latent heat flux which will sublimate all snow and ice
184  C     over one time step  C     over one time step
185        _RL latentHeatFluxMax   (1:sNx,1:sNy)        _RL latentHeatFluxMax   (1:sNx,1:sNy)
# Line 192  C     AREA_PRE :: hold sea-ice fraction Line 201  C     AREA_PRE :: hold sea-ice fraction
201        _RL AREApreTH           (1:sNx,1:sNy)        _RL AREApreTH           (1:sNx,1:sNy)
202        _RL HEFFpreTH           (1:sNx,1:sNy)        _RL HEFFpreTH           (1:sNx,1:sNy)
203        _RL HSNWpreTH           (1:sNx,1:sNy)        _RL HSNWpreTH           (1:sNx,1:sNy)
204    #ifdef SEAICE_ITD
205          _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
206          _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
207          _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
208          _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
209          _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)
210    #endif
211    
212  C     wind speed  C     wind speed
213        _RL UG                  (1:sNx,1:sNy)        _RL UG                  (1:sNx,1:sNy)
# Line 217  C temporary variables available for the Line 233  C temporary variables available for the
233    
234        INTEGER ilockey        INTEGER ilockey
235        INTEGER it        INTEGER it
236    #ifdef SEAICE_ITD
237          INTEGER K
238    #endif
239        _RL pFac        _RL pFac
240        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
241        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
242        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
243    #ifdef SEAICE_ITD
244          _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
245          _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
246    #endif
247        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
248        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
249        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
250    #ifdef SEAICE_ITD
251          _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
252          _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
253    #endif
254  C     Helper variables: reciprocal of some constants  C     Helper variables: reciprocal of some constants
255        _RL recip_multDim        _RL recip_multDim
256        _RL recip_deltaTtherm        _RL recip_deltaTtherm
# Line 259  C ====================================== Line 286  C ======================================
286        ENDIF        ENDIF
287    
288  C     avoid unnecessary divisions in loops  C     avoid unnecessary divisions in loops
289    #ifdef SEAICE_ITD
290    CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
291    #endif
292        recip_multDim        = SEAICE_multDim        recip_multDim        = SEAICE_multDim
293        recip_multDim        = ONE / recip_multDim        recip_multDim        = ONE / recip_multDim
294  C     above/below: double/single precision calculation of recip_multDim  C     above/below: double/single precision calculation of recip_multDim
# Line 337  C ===================== Line 367  C =====================
367            d_HEFFbyRLX(I,J)           = 0.0 _d 0            d_HEFFbyRLX(I,J)           = 0.0 _d 0
368  #endif  #endif
369    
370    #ifdef SEAICE_ITD
371              d_AREAbyNEG(I,J)           = 0.0 _d 0
372    #endif
373            d_HEFFbyNEG(I,J)           = 0.0 _d 0            d_HEFFbyNEG(I,J)           = 0.0 _d 0
374            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
375            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
# Line 370  c Line 403  c
403              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
404              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
405              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
406  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_ITD
407                r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
408                r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
409    #endif
410    #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
411              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
412  #endif  #endif
413            ENDDO            ENDDO
# Line 406  C ====================================== Line 443  C ======================================
443  #endif  #endif
444           ENDDO           ENDDO
445          ENDDO          ENDDO
446    #ifdef SEAICE_ITD
447            DO K=1,nITD
448             DO J=1,sNy
449              DO I=1,sNx
450               HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
451               HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
452               AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
453              ENDDO
454             ENDDO
455            ENDDO
456    #endif
457    
458  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
459    
# Line 433  C         d_HEFFbyRLX(i,j) = 0. _d 0 Line 481  C         d_HEFFbyRLX(i,j) = 0. _d 0
481  C           d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j)  C           d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j)
482              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
483             ENDIF             ENDIF
484    #ifdef SEAICE_ITD
485                AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
486         &                           +  d_AREAbyRLX(i,j)
487                HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
488         &                           +  d_HEFFbyRLX(i,j)
489    #endif
490              AREA(I,J,bi,bj) = AREA(I,J,bi,bj) +  d_AREAbyRLX(i,j)              AREA(I,J,bi,bj) = AREA(I,J,bi,bj) +  d_AREAbyRLX(i,j)
491              HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +  d_HEFFbyRLX(i,j)              HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +  d_HEFFbyRLX(i,j)
492           ENDDO           ENDDO
# Line 449  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 503  CADJ STORE area(:,:,bi,bj) = comlev1_bib
503  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
504          DO J=1,sNy          DO J=1,sNy
505           DO I=1,sNx           DO I=1,sNx
506    #ifdef SEAICE_ITD
507              DO K=1,nITD
508               tmpscal1=0. _d 0
509               tmpscal2=0. _d 0
510               tmpscal3=0. _d 0
511               tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)
512               HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
513               d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
514               tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)
515               HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
516               d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
517               tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
518               AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1
519               d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1
520              ENDDO
521    CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
522    C         by calling SEAICE_ITD_SUM
523    #else
524            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)
           HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)  
525            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
           HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)  
526            AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)            AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
527              HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
528              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
529    #endif
530           ENDDO           ENDDO
531          ENDDO          ENDDO
532    
# Line 464  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 537  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
537  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
538          DO J=1,sNy          DO J=1,sNy
539           DO I=1,sNx           DO I=1,sNx
540            tmpscal2=0. _d 0  #ifdef SEAICE_ITD
541            tmpscal3=0. _d 0            DO K=1,nITD
542    #endif
543              tmpscal2=0. _d 0
544              tmpscal3=0. _d 0
545    #ifdef SEAICE_ITD
546               IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
547                tmpscal2=-HEFFITD(I,J,K,bi,bj)
548                tmpscal3=-HSNOWITD(I,J,K,bi,bj)
549                TICES(I,J,K,bi,bj)=celsius2K
550    CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
551               ENDIF
552               HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2
553               HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
554    #else
555            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
556              tmpscal2=-HEFF(I,J,bi,bj)             tmpscal2=-HEFF(I,J,bi,bj)
557              tmpscal3=-HSNOW(I,J,bi,bj)             tmpscal3=-HSNOW(I,J,bi,bj)
558              TICE(I,J,bi,bj)=celsius2K             TICE(I,J,bi,bj)=celsius2K
559              DO IT=1,SEAICE_multDim             DO IT=1,SEAICE_multDim
560                TICES(I,J,IT,bi,bj)=celsius2K              TICES(I,J,IT,bi,bj)=celsius2K
561              ENDDO             ENDDO
562            ENDIF            ENDIF
563            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
           d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2  
564            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
565    #endif
566              d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
567            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
568    #ifdef SEAICE_ITD
569              ENDDO
570    #endif
571           ENDDO           ENDDO
572          ENDDO          ENDDO
573    
# Line 489  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 579  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
579  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
580          DO J=1,sNy          DO J=1,sNy
581           DO I=1,sNx           DO I=1,sNx
582    #ifdef SEAICE_ITD
583              DO K=1,nITD
584               IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.
585         &         (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))
586         &      AREAITD(I,J,K,bi,bj)=0. _d 0
587              ENDDO
588    #else
589            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
590       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
591    #endif
592           ENDDO           ENDDO
593          ENDDO          ENDDO
594    
# Line 502  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 600  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
600  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
601          DO J=1,sNy          DO J=1,sNy
602           DO I=1,sNx           DO I=1,sNx
603    #ifdef SEAICE_ITD
604              DO K=1,nITD
605               IF ((HEFFITD(i,j,k,bi,bj).GT.0).OR.
606         &         (HSNOWITD(i,j,k,bi,bj).GT.0)) THEN
607    CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1
608    C           hence use SEAICE_area_floor devided by nITD
609    C           (or install a warning in e.g. seaice_readparms.F)
610                AREAITD(I,J,K,bi,bj)=
611         &         MAX(AREAITD(I,J,K,bi,bj),SEAICE_area_floor/float(nITD))
612               ENDIF
613              ENDDO
614    #else
615            IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN            IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
616             AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)             AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
617            ENDIF            ENDIF
618    #endif
619           ENDDO           ENDDO
620          ENDDO          ENDDO
621  #endif /* DISABLE_AREA_FLOOR */  #endif /* DISABLE_AREA_FLOOR */
622    
623  C 2.5) treat case of excessive ice cover, e.g., due to ridging:  C 2.5) treat case of excessive ice cover, e.g., due to ridging:
624    
625    CToM   for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
626    C      which is called at end of PART 1 below
627    #ifndef SEAICE_ITD
628  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
629  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
630  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 525  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 639  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
639            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)
640           ENDDO           ENDDO
641          ENDDO          ENDDO
642    #endif /* SEAICE_ITD */
643    
644    #ifdef SEAICE_ITD
645    CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
646    C    first, update AREA and HEFF:
647            CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
648    C
649            DO J=1,sNy
650             DO I=1,sNx
651    C    TICES was changed above (item 1.25), now update TICE as ice volume
652    C     weighted average of TICES
653              tmpscal1 = 0. _d 0
654              tmpscal2 = 0. _d 0
655              DO K=1,nITD
656               tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)
657               tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)
658              ENDDO
659              TICE(I,J,bi,bj)=tmpscal1/tmpscal2
660    C    lines of item 2.5 that were omitted:
661    C    in 2.5 these lines are executed before "ridging" is applied to AREA
662    C    hence we execute them here before SEAICE_ITD_REDIST is called
663    C    although this means that AREA has not been completely regularized
664    #ifdef ALLOW_DIAGNOSTICS
665              DIAGarrayA(I,J) = AREA(I,J,bi,bj)
666    #endif
667    #ifdef ALLOW_SITRACER
668              SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
669    #endif
670             ENDDO
671            ENDDO
672    
673    CToM finally make sure that all categories meet their thickness limits
674    C    which includes ridging as in item 2.5
675    C    and update AREA, HEFF, and HSNOW
676            CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
677            CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
678    
679    #endif
680  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
681    C        ENDIF SEAICEadjMODE.EQ.0
682          ENDIF          ENDIF
683  #endif  #endif
684    
# Line 547  C 3) store regularized values of heff, h Line 699  C 3) store regularized values of heff, h
699  #endif  #endif
700           ENDDO           ENDDO
701          ENDDO          ENDDO
702    #ifdef SEAICE_ITD
703            DO K=1,nITD
704             DO J=1,sNy
705              DO I=1,sNx
706               HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
707               HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
708               AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
709    
710    C memorize areal and volume fraction of each ITD category
711               IF (AREA(I,J,bi,bj).GT.0) THEN
712                areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
713               ELSE
714                areaFracFactor(I,J,K)=ZERO
715               ENDIF
716               IF (HEFF(I,J,bi,bj).GT.0) THEN
717                heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)
718               ELSE
719                heffFracFactor(I,J,K)=ZERO
720               ENDIF
721              ENDDO
722             ENDDO
723            ENDDO
724    C prepare SItrHEFF to be computed as cumulative sum
725            DO K=2,5
726             DO J=1,sNy
727              DO I=1,sNx
728               SItrHEFF(I,J,bi,bj,K)=ZERO
729              ENDDO
730             ENDDO
731            ENDDO
732    C prepare SItrAREA to be computed as cumulative sum
733            DO J=1,sNy
734             DO I=1,sNx
735              SItrAREA(I,J,bi,bj,3)=ZERO
736             ENDDO
737            ENDDO
738    #endif
739    
740  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
741  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
# Line 601  Cgf no additional dependency of air-sea Line 790  Cgf no additional dependency of air-sea
790             AREApreTH(I,J) = 0. _d 0             AREApreTH(I,J) = 0. _d 0
791            ENDDO            ENDDO
792           ENDDO           ENDDO
793    #ifdef SEAICE_ITD
794             DO K=1,nITD
795              DO J=1,sNy
796               DO I=1,sNx
797                HEFFITDpreTH(I,J,K) = 0. _d 0
798                HSNWITDpreTH(I,J,K) = 0. _d 0
799                AREAITDpreTH(I,J,K) = 0. _d 0
800               ENDDO
801              ENDDO
802             ENDDO
803    #endif
804          ENDIF          ENDIF
805  #endif  #endif
806    
# Line 620  CADJ STORE AREApreTH = comlev1_bibj, key Line 820  CADJ STORE AREApreTH = comlev1_bibj, key
820  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
821  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
822  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
823    #ifdef SEAICE_ITD
824            DO K=1,nITD
825    #endif
826          DO J=1,sNy          DO J=1,sNy
827           DO I=1,sNx           DO I=1,sNx
828    
829    #ifdef SEAICE_ITD
830              IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
831    #ifdef SEAICE_GROWTH_LEGACY
832                tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
833         &                              AREAITDpreTH(I,J,K))
834                hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
835                tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1
836                heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)
837    #else /* SEAICE_GROWTH_LEGACY */
838    cif        regularize AREA with SEAICE_area_reg
839               tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)
840         &                     + area_reg_sq)
841    cif        heffActual calculated with the regularized AREA
842               tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1
843    cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
844               heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2
845         &                                  + hice_reg_sq)
846    cif        hsnowActual calculated with the regularized AREA
847               hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1
848    #endif /* SEAICE_GROWTH_LEGACY */
849    cif        regularize the inverse of heffActual by hice_reg
850               recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /
851         &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)
852         &                      + hice_reg_sq)
853    cif       Do not regularize when HEFFpreTH = 0
854              ELSE
855                heffActualMult(I,J,K) = ZERO
856                hsnowActualMult(I,J,K) = ZERO
857                recip_heffActualMult(I,J,K)  = ZERO
858              ENDIF
859    #else /* SEAICE_ITD */
860            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
861  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
862              tmpscal1         = MAX(SEAICE_area_reg,AREApreTH(I,J))              tmpscal1         = MAX(SEAICE_area_reg,AREApreTH(I,J))
# Line 648  cif       Do not regularize when HEFFpre Line 882  cif       Do not regularize when HEFFpre
882              hsnowActual(I,J) = ZERO              hsnowActual(I,J) = ZERO
883              recip_heffActual(I,J)  = ZERO              recip_heffActual(I,J)  = ZERO
884            ENDIF            ENDIF
885    #endif /* SEAICE_ITD */
886    
887           ENDDO           ENDDO
888          ENDDO          ENDDO
889    #ifdef SEAICE_ITD
890            ENDDO
891    #endif
892    
893  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
894          CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)          CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
# Line 661  cif       Do not regularize when HEFFpre Line 899  cif       Do not regularize when HEFFpre
899  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
900  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
901  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
902    #ifdef SEAICE_ITD
903            DO K=1,nITD
904    #endif
905          DO J=1,sNy          DO J=1,sNy
906           DO I=1,sNx           DO I=1,sNx
907  c        The latent heat flux over the sea ice which  c        The latent heat flux over the sea ice which
908  c        will sublimate all of the snow and ice over one time  c        will sublimate all of the snow and ice over one time
909  c        step (W/m^2)  c        step (W/m^2)
910    #ifdef SEAICE_ITD
911              IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
912                latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *
913         &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +
914         &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)
915              ELSE
916                latentHeatFluxMaxMult(I,J,K) = ZERO
917              ENDIF
918    #else
919            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
920              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
921       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +
# Line 673  c        step (W/m^2) Line 923  c        step (W/m^2)
923            ELSE            ELSE
924              latentHeatFluxMax(I,J) = ZERO              latentHeatFluxMax(I,J) = ZERO
925            ENDIF            ENDIF
926    #endif
927           ENDDO           ENDDO
928          ENDDO          ENDDO
929    #ifdef SEAICE_ITD
930            ENDDO
931    #endif
932  #endif /* SEAICE_CAP_SUBLIM */  #endif /* SEAICE_CAP_SUBLIM */
933    
934  C ===================================================================  C ===================================================================
# Line 746  CADJ &                       key = iicek Line 1000  CADJ &                       key = iicek
1000  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1001    
1002  C--   Start loop over multi-categories  C--   Start loop over multi-categories
1003    #ifdef SEAICE_ITD
1004    CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
1005    #endif
1006          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1007  c        homogeneous distribution between 0 and 2 x heffActual  c        homogeneous distribution between 0 and 2 x heffActual
1008    #ifndef SEAICE_ITD
1009           pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim           pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
1010    #endif
1011           DO J=1,sNy           DO J=1,sNy
1012            DO I=1,sNx            DO I=1,sNx
1013    #ifndef SEAICE_ITD
1014    CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
1015    C    (instead of heffActual and latentHeatFluxMax)
1016             heffActualMult(I,J,IT)= heffActual(I,J)*pFac             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1017  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1018             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1019  #endif  #endif
1020    #endif
1021             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1022             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1023             TICE(I,J,bi,bj) = ZERO             TICE(I,J,bi,bj) = ZERO
# Line 780  CADJ &     comlev1_bibj, key = iicekey, Line 1043  CADJ &     comlev1_bibj, key = iicekey,
1043    
1044          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1045           CALL SEAICE_SOLVE4TEMP(           CALL SEAICE_SOLVE4TEMP(
1046    #ifdef SEAICE_ITD
1047         I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1048    #else
1049       I        UG, heffActualMult(1,1,IT), hsnowActual,       I        UG, heffActualMult(1,1,IT), hsnowActual,
1050    #endif
1051  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1052       I        latentHeatFluxMaxMult(1,1,IT),       I        latentHeatFluxMaxMult(1,1,IT),
1053  #endif  #endif
# Line 809  CADJ &     comlev1_bibj, key = iicekey, Line 1076  CADJ &     comlev1_bibj, key = iicekey,
1076           DO J=1,sNy           DO J=1,sNy
1077            DO I=1,sNx            DO I=1,sNx
1078  C     update TICE & TICES  C     update TICE & TICES
1079    #ifdef SEAICE_ITD
1080    C     calculate area weighted mean
1081    C     (although the ice's temperature relates to its energy content
1082    C      and hence should be averaged weighted by ice volume [heffFracFactor],
1083    C      the temperature here is a result of the fluxes through the ice surface
1084    C      computed individually for each single category in SEAICE_SOLVE4TEMP
1085    C      and hence is averaged area weighted [areaFracFactor])
1086               TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1087         &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)
1088    #else
1089             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1090       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
1091    #endif
1092             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1093  C     average over categories  C     average over categories
1094    #ifdef SEAICE_ITD
1095    C     calculate area weighted mean
1096    C     (fluxes are per unit (ice surface) area and are thus area weighted)
1097               a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1098         &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1099               a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1100         &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1101               a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1102         &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)
1103    #else
1104             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1105       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1106             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1107       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1108             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1109       &          + a_FWbySublimMult(I,J,IT)*recip_multDim       &          + a_FWbySublimMult(I,J,IT)*recip_multDim
1110    #endif
1111            ENDDO            ENDDO
1112           ENDDO           ENDDO
1113          ENDDO          ENDDO
# Line 851  CADJ STORE a_FWbySublim    = comlev1_bib Line 1140  CADJ STORE a_FWbySublim    = comlev1_bib
1140  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1141    
1142  C switch heat fluxes from W/m2 to 'effective' ice meters  C switch heat fluxes from W/m2 to 'effective' ice meters
1143    #ifdef SEAICE_ITD
1144            DO K=1,nITD
1145             DO J=1,sNy
1146              DO I=1,sNx
1147               a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)
1148         &          * convertQ2HI * AREAITDpreTH(I,J,K)
1149               a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)
1150         &          * convertQ2HI * AREAITDpreTH(I,J,K)
1151    C and initialize r_QbyATM_cover
1152               r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)
1153    C     Convert fresh water flux by sublimation to 'effective' ice meters.
1154    C     Negative sublimation is resublimation and will be added as snow.
1155    #ifdef SEAICE_DISABLE_SUBLIM
1156               a_FWbySublimMult(I,J,K) = ZERO
1157    #endif
1158               a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce
1159         &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)
1160               r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)
1161              ENDDO
1162             ENDDO
1163            ENDDO
1164            DO J=1,sNy
1165             DO I=1,sNx
1166              a_QbyATM_open(I,J)    = a_QbyATM_open(I,J)
1167         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1168              a_QSWbyATM_open(I,J)  = a_QSWbyATM_open(I,J)
1169         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1170    C and initialize r_QbyATM_open
1171              r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1172             ENDDO
1173            ENDDO
1174    #else /* SEAICE_ITD */
1175          DO J=1,sNy          DO J=1,sNy
1176           DO I=1,sNx           DO I=1,sNx
1177            a_QbyATM_cover(I,J)   = a_QbyATM_cover(I,J)            a_QbyATM_cover(I,J)   = a_QbyATM_cover(I,J)
# Line 869  C     Negative sublimation is resublimat Line 1190  C     Negative sublimation is resublimat
1190  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1191  cgf just for those who may need to omit this term to reproduce old results  cgf just for those who may need to omit this term to reproduce old results
1192            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1193  #endif /* SEAICE_CAP_SUBLIM */  #endif
1194            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1195       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1196            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1197           ENDDO           ENDDO
1198          ENDDO          ENDDO
1199    #endif /* SEAICE_ITD */
1200    
1201  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1202  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte
# Line 891  CADJ STORE r_FWbySublim    = comlev1_bib Line 1213  CADJ STORE r_FWbySublim    = comlev1_bib
1213  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1214  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1215          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1216    #ifdef SEAICE_ITD
1217             DO K=1,nITD
1218              DO J=1,sNy
1219               DO I=1,sNx
1220                a_QbyATMmult_cover(I,J,K)   = 0. _d 0
1221                r_QbyATMmult_cover(I,J,K)   = 0. _d 0
1222                a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
1223               ENDDO
1224              ENDDO
1225             ENDDO
1226    #else
1227           DO J=1,sNy           DO J=1,sNy
1228            DO I=1,sNx            DO I=1,sNx
1229             a_QbyATM_cover(I,J)   = 0. _d 0             a_QbyATM_cover(I,J)   = 0. _d 0
# Line 898  Cgf no additional dependency through ice Line 1231  Cgf no additional dependency through ice
1231             a_QSWbyATM_cover(I,J) = 0. _d 0             a_QSWbyATM_cover(I,J) = 0. _d 0
1232            ENDDO            ENDDO
1233           ENDDO           ENDDO
1234    #endif
1235          ENDIF          ENDIF
1236  #endif  #endif
1237    
# Line 961  C ====================================== Line 1295  C ======================================
1295  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1296  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1297  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1298    #ifdef SEAICE_ITD
1299            DO K=1,nITD
1300    #endif
1301          DO J=1,sNy          DO J=1,sNy
1302           DO I=1,sNx           DO I=1,sNx
1303  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1304            tmpscal2 =            tmpscal2 =
1305    #ifdef SEAICE_ITD
1306         &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)
1307         &             *SNOW2ICE),ZERO)
1308    C         accumulate change over ITD categories
1309              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1310         &                                                       *ICE2SNOW
1311              HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2
1312         &                                                       *ICE2SNOW
1313              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1314    C         keep total up to date, too
1315              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1316    #else
1317       &     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)
1318            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1319            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW
1320            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1321    #endif
1322           ENDDO           ENDDO
1323          ENDDO          ENDDO
1324  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 979  CADJ STORE r_FWbySublim    = comlev1_bib Line 1329  CADJ STORE r_FWbySublim    = comlev1_bib
1329           DO I=1,sNx           DO I=1,sNx
1330  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1331            tmpscal2 =            tmpscal2 =
1332    #ifdef SEAICE_ITD
1333         &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)
1334    C         accumulate change over ITD categories
1335              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1336              HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2
1337              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1338    C         keep total up to date, too
1339              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1340    #else
1341       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1342            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
1343            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2
1344            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1345    #endif
1346           ENDDO           ENDDO
1347          ENDDO          ENDDO
1348          DO J=1,sNy          DO J=1,sNy
1349           DO I=1,sNx           DO I=1,sNx
1350  C     If anything is left, it will be evaporated from the ocean rather than sublimated.  C     If anything is left, it will be evaporated from the ocean rather than sublimated.
1351  C     Since a_QbyATM_cover was computed for sublimation, not simple evapation, we need to  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
1352  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).
1353    #ifdef SEAICE_ITD
1354              a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
1355         &                              - r_FWbySublimMult(I,J,K)
1356              r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)
1357         &                              - r_FWbySublimMult(I,J,K)
1358             ENDDO
1359            ENDDO
1360    C       end K loop
1361            ENDDO
1362    C       then update totals      
1363            DO J=1,sNy
1364             DO I=1,sNx
1365    #endif
1366            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)
1367            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)
1368           ENDDO           ENDDO
1369          ENDDO          ENDDO
1370    c ToM<<< debug seaice_growth
1371            WRITE(msgBuf,'(A,7F6.2)')
1372    #ifdef SEAICE_ITD
1373         &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1374         &     HEFFITD(20,20,:,bi,bj)
1375    #else
1376         &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1377         &     HEFF(20,20,bi,bj)
1378    #endif
1379            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1380         &    SQUEEZE_RIGHT , myThid)
1381    c ToM>>>
1382    
1383  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1384  C ===========================================================  C ===========================================================
# Line 1003  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1388  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1388  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1389  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1390    
1391    #ifdef SEAICE_ITD
1392            DO K=1,nITD
1393             DO J=1,sNy
1394              DO I=1,sNx
1395    C          ice growth/melt due to ocean heat is equally distributed under the ice
1396    C          and hence weighted by fractional area of each thickness category
1397               tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
1398         &                               -HEFFITD(I,J,K,bi,bj))
1399               d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
1400               r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1
1401               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj)  + tmpscal1
1402    #ifdef ALLOW_SITRACER
1403               SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1404         &                           + HEFFITD(I,J,K,bi,bj)
1405    #endif
1406              ENDDO
1407             ENDDO
1408            ENDDO
1409    #else /* SEAICE_ITD */
1410          DO J=1,sNy          DO J=1,sNy
1411           DO I=1,sNx           DO I=1,sNx
1412            d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))            d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
# Line 1013  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1417  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1417  #endif  #endif
1418           ENDDO           ENDDO
1419          ENDDO          ENDDO
1420    #endif /* SEAICE_ITD */
1421    c ToM<<< debug seaice_growth
1422            WRITE(msgBuf,'(A,7F6.2)')
1423    #ifdef SEAICE_ITD
1424         &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1425         &     HEFFITD(20,20,:,bi,bj)
1426    #else
1427         &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1428         &     HEFF(20,20,bi,bj)
1429    #endif
1430            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1431         &    SQUEEZE_RIGHT , myThid)
1432    c ToM>>>
1433    
1434  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1435  C ==================================================================  C ==================================================================
# Line 1022  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1439  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1439  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1440  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1441    
1442    #ifdef SEAICE_ITD
1443            DO K=1,nITD
1444             DO J=1,sNy
1445              DO I=1,sNx
1446    C     Convert to standard units (meters of ice) rather than to meters
1447    C     of snow. This appears to be more robust.
1448              tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)
1449         &                                           *SNOW2ICE)
1450              tmpscal2=MIN(tmpscal1,0. _d 0)
1451    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1452    Cgf no additional dependency through snow
1453              IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1454    #endif
1455              d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW
1456              HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW
1457              r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
1458    C         keep the total up to date, too
1459              r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1460              ENDDO
1461             ENDDO
1462            ENDDO
1463    #else /* SEAICE_ITD */
1464          DO J=1,sNy          DO J=1,sNy
1465           DO I=1,sNx           DO I=1,sNx
1466  C     Convert to standard units (meters of ice) rather than to meters  C     Convert to standard units (meters of ice) rather than to meters
# Line 1037  Cgf no additional dependency through sno Line 1476  Cgf no additional dependency through sno
1476            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1477           ENDDO           ENDDO
1478          ENDDO          ENDDO
1479    #endif /* SEAICE_ITD */
1480    c ToM<<< debug seaice_growth
1481            WRITE(msgBuf,'(A,7F6.2)')
1482    #ifdef SEAICE_ITD
1483         &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1484         &     HEFFITD(20,20,:,bi,bj)
1485    #else
1486         &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1487         &     HEFF(20,20,bi,bj)
1488    #endif
1489            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1490         &    SQUEEZE_RIGHT , myThid)
1491    c ToM>>>
1492    
1493  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1494  C ====================================================  C ====================================================
# Line 1051  Cgf where all experiments start in Janua Line 1503  Cgf where all experiments start in Janua
1503  Cgf the v1.81=>v1.82 revision would change results in  Cgf the v1.81=>v1.82 revision would change results in
1504  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1505    
1506    #ifdef SEAICE_ITD
1507            DO K=1,nITD
1508             DO J=1,sNy
1509              DO I=1,sNx
1510    #ifdef SEAICE_GROWTH_LEGACY
1511               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
1512    #else
1513               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
1514    c         Limit ice growth by potential melt by ocean
1515         &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))
1516    #endif /* SEAICE_GROWTH_LEGACY */
1517               d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1518         &                                 + tmpscal2
1519               d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1520         &                                 + tmpscal2
1521               r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
1522         &                                 - tmpscal2
1523               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2
1524    
1525    #ifdef ALLOW_SITRACER
1526               SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1527         &                           + HEFFITD(I,J,K,bi,bj)
1528    #endif
1529              ENDDO
1530             ENDDO
1531            ENDDO
1532    #else /* SEAICE_ITD */
1533          DO J=1,sNy          DO J=1,sNy
1534           DO I=1,sNx           DO I=1,sNx
1535    
# Line 1070  c         Limit ice growth by potential Line 1549  c         Limit ice growth by potential
1549  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1550            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1551  #endif  #endif
1552           ENDDO           ENDDO
1553          ENDDO          ENDDO
1554    #endif /* SEAICE_ITD */
1555    c ToM<<< debug seaice_growth
1556            WRITE(msgBuf,'(A,7F6.2)')
1557    #ifdef SEAICE_ITD
1558         &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1559         &     HEFFITD(20,20,:,bi,bj)
1560    #else
1561         &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1562         &     HEFF(20,20,bi,bj)
1563    #endif
1564            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1565         &    SQUEEZE_RIGHT , myThid)
1566    c ToM>>>
1567    
1568  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1569  C depending on atmospheric conditions.  C depending on atmospheric conditions.
# Line 1098  C           add precip to the fresh wate Line 1590  C           add precip to the fresh wate
1590       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1591              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
1592            ENDIF            ENDIF
1593             ENDDO
1594            ENDDO
1595    #ifdef SEAICE_ITD
1596            DO K=1,nITD
1597             DO J=1,sNy
1598              DO I=1,sNx
1599               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)
1600         &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
1601              ENDDO
1602             ENDDO
1603            ENDDO
1604    #else
1605            DO J=1,sNy
1606             DO I=1,sNx
1607            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1608           ENDDO           ENDDO
1609          ENDDO          ENDDO
1610    #endif
1611  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1612  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1613  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1614  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1615    c ToM<<< debug seaice_growth
1616            WRITE(msgBuf,'(A,7F6.2)')
1617    #ifdef SEAICE_ITD
1618         &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1619         &     HEFFITD(20,20,:,bi,bj)
1620    #else
1621         &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1622         &     HEFF(20,20,bi,bj)
1623    #endif
1624            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1625         &    SQUEEZE_RIGHT , myThid)
1626    c ToM>>>
1627    
1628  C compute snow melt due to heat available from ocean.  C compute snow melt due to heat available from ocean.
1629  C =================================================================  C =================================================================
# Line 1116  Cph( very sensitive bit here by JZ Line 1635  Cph( very sensitive bit here by JZ
1635  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1636  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1637  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1638    
1639    #ifdef SEAICE_ITD
1640            DO K=1,nITD
1641             DO J=1,sNy
1642              DO I=1,sNx
1643               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
1644         &                  -HSNOWITD(I,J,K,bi,bj))
1645               tmpscal2=MIN(tmpscal1,0. _d 0)
1646    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1647    Cgf no additional dependency through snow
1648               if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1649    #endif
1650               d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1651               r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1652               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
1653              ENDDO
1654             ENDDO
1655            ENDDO
1656    #else /* SEAICE_ITD */
1657          DO J=1,sNy          DO J=1,sNy
1658           DO I=1,sNx           DO I=1,sNx
1659            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
# Line 1130  Cgf no additional dependency through sno Line 1668  Cgf no additional dependency through sno
1668            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1669           ENDDO           ENDDO
1670          ENDDO          ENDDO
1671    #endif /* SEAICE_ITD */
1672  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1673  Cph)  Cph)
1674    c ToM<<< debug seaice_growth
1675            WRITE(msgBuf,'(A,7F6.2)')
1676    #ifdef SEAICE_ITD
1677         &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1678         &     HEFFITD(20,20,:,bi,bj)
1679    #else
1680         &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1681         &     HEFF(20,20,bi,bj)
1682    #endif
1683            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1684         &    SQUEEZE_RIGHT , myThid)
1685    c ToM>>>
1686    
1687  C gain of new ice over open water  C gain of new ice over open water
1688  C ===============================  C ===============================
# Line 1159  C           or 0. otherwise (no melting Line 1710  C           or 0. otherwise (no melting
1710            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1711            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1712            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1713    #ifdef SEAICE_ITD
1714    C         open water area fraction
1715              tmpscal0 = ONE-AREApreTH(I,J)
1716    C         determine thickness of new ice
1717    C         considering the entire open water area to refreeze
1718              tmpscal1 = tmpscal3/tmpscal0
1719    C         then add new ice volume to appropriate thickness category
1720              DO K=1,nITD
1721               IF (tmpscal1.LT.Hlimit(K)) THEN
1722                HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
1723                tmpscal3=ZERO
1724    C not sure if AREAITD should be changed here since AREA is incremented
1725    C   in PART 4 below in non-itd code
1726    C in this scenario all open water is covered by new ice instantaneously,
1727    C   i.e. no delayed lead closing is concidered such as is associated with
1728    C   Hibler's h_0 parameter
1729                AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
1730         &                           + tmpscal0
1731                tmpscal0=ZERO
1732               ENDIF
1733              ENDDO
1734    #else
1735            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1736    #endif
1737           ENDDO           ENDDO
1738          ENDDO          ENDDO
1739    
1740  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1741    #ifdef SEAICE_ITD
1742            DO K=1,nITD
1743             DO J=1,sNy
1744              DO I=1,sNx
1745    c needs to be here to allow use also with LEGACY branch
1746               SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1747         &                           + HEFFITD(I,J,K,bi,bj)
1748              ENDDO
1749             ENDDO
1750            ENDDO
1751    #else
1752          DO J=1,sNy          DO J=1,sNy
1753           DO I=1,sNx           DO I=1,sNx
1754  c needs to be here to allow use also with LEGACY branch  c needs to be here to allow use also with LEGACY branch
1755            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1756           ENDDO           ENDDO
1757          ENDDO          ENDDO
1758    #endif
1759  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1760    c ToM<<< debug seaice_growth
1761            WRITE(msgBuf,'(A,7F6.2)')
1762    #ifdef SEAICE_ITD
1763         &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1764         &     HEFFITD(20,20,:,bi,bj)
1765    #else
1766         &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1767         &     HEFF(20,20,bi,bj)
1768    #endif
1769            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1770         &    SQUEEZE_RIGHT , myThid)
1771    c ToM>>>
1772    
1773  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1774  C =================================  C =================================
# Line 1182  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1780  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1780  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1781  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1782          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1783    #ifdef SEAICE_ITD
1784             DO K=1,nITD
1785              DO J=1,sNy
1786               DO I=1,sNx
1787                tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow
1788         &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1789                tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
1790                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
1791                HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1
1792                HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
1793         &                            * ICE2SNOW
1794               ENDDO
1795              ENDDO
1796             ENDDO
1797    #else
1798           DO J=1,sNy           DO J=1,sNy
1799            DO I=1,sNx            DO I=1,sNx
1800             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
# Line 1191  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1804  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1804             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1805             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1806       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
1807            ENDDO            ENDDO
1808           ENDDO           ENDDO
1809    #endif
1810          ENDIF          ENDIF
1811  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1812    c ToM<<< debug seaice_growth
1813            WRITE(msgBuf,'(A,7F6.2)')
1814    #ifdef SEAICE_ITD
1815         &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1816         &     HEFFITD(20,20,:,bi,bj)
1817    #else
1818         &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1819         &     HEFF(20,20,bi,bj)
1820    #endif
1821            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1822         &    SQUEEZE_RIGHT , myThid)
1823    c ToM>>>
1824    
1825  C ===================================================================  C ===================================================================
1826  C ==========PART 4: determine ice cover fraction increments=========-  C ==========PART 4: determine ice cover fraction increments=========-
# Line 1220  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi Line 1846  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
1846  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1847  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1848    
1849    #ifdef SEAICE_ITD
1850    C       replaces Hibler '79 scheme and lead closing parameter
1851    C       because ITD accounts explicitly for lead openings and
1852    C       different melt rates due to varying ice thickness
1853    C
1854    C       only consider ice area loss due to total ice thickness loss
1855    C       ice area gain due to freezing of open water as handled above
1856    C       under "gain of new ice over open water"
1857    C
1858    C       does not account for lateral melt of ice floes
1859    C
1860    C        AREAITD is incremented in section "gain of new ice over open water" above
1861    C
1862            DO K=1,nITD
1863             DO J=1,sNy
1864              DO I=1,sNx
1865               IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
1866                AREAITD(I,J,K,bi,bj)=ZERO
1867               ENDIF
1868    #ifdef ALLOW_SITRACER
1869               SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1870         &                           + AREAITD(I,J,K,bi,bj)
1871    #endif /* ALLOW_SITRACER */
1872              ENDDO
1873             ENDDO
1874            ENDDO
1875    #else /* SEAICE_ITD */
1876          DO J=1,sNy          DO J=1,sNy
1877           DO I=1,sNx           DO I=1,sNx
1878    
# Line 1289  C apply tendency Line 1942  C apply tendency
1942  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
1943           ENDDO           ENDDO
1944          ENDDO          ENDDO
1945    #endif /* SEAICE_ITD */
1946    
1947  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1948  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
1949          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
1950    #ifdef SEAICE_ITD
1951             DO K=1,nITD
1952              DO J=1,sNy
1953               DO I=1,sNx
1954                AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *
1955         &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )
1956               ENDDO
1957              ENDDO
1958             ENDDO
1959    #else
1960           DO J=1,sNy           DO J=1,sNy
1961            DO I=1,sNx            DO I=1,sNx
1962  C            AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)  C            AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
# Line 1300  C            AREA(I,J,bi,bj) = 0.1 _d 0 Line 1964  C            AREA(I,J,bi,bj) = 0.1 _d 0
1964       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1965            ENDDO            ENDDO
1966           ENDDO           ENDDO
1967    #endif
1968          ENDIF          ENDIF
1969  #endif  #endif
1970    #ifdef SEAICE_ITD
1971    C check categories for consistency with limits after growth/melt
1972            CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
1973    C finally update total AREA, HEFF, HSNOW
1974            CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
1975    #endif
1976    
1977  C ===================================================================  C ===================================================================
1978  C =============PART 5: determine ice salinity increments=============  C =============PART 5: determine ice salinity increments=============
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22