/[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.3 by torge, Wed Sep 26 17:50: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    C and initialize r_QbyATM_open
1167              r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1168             ENDDO
1169            ENDDO
1170    #else /* SEAICE_ITD */
1171          DO J=1,sNy          DO J=1,sNy
1172           DO I=1,sNx           DO I=1,sNx
1173            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 1186  C     Negative sublimation is resublimat
1186  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1187  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
1188            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1189  #endif /* SEAICE_CAP_SUBLIM */  #endif
1190            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1191       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1192            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1193           ENDDO           ENDDO
1194          ENDDO          ENDDO
1195    #endif /* SEAICE_ITD */
1196    
1197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1198  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 1209  CADJ STORE r_FWbySublim    = comlev1_bib
1209  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1210  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1211          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1212    #ifdef SEAICE_ITD
1213             DO K=1,nITD
1214              DO J=1,sNy
1215               DO I=1,sNx
1216                a_QbyATMmult_cover(I,J,K)   = 0. _d 0
1217                r_QbyATMmult_cover(I,J,K)   = 0. _d 0
1218                a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
1219               ENDDO
1220              ENDDO
1221             ENDDO
1222    #else
1223           DO J=1,sNy           DO J=1,sNy
1224            DO I=1,sNx            DO I=1,sNx
1225             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 1227  Cgf no additional dependency through ice
1227             a_QSWbyATM_cover(I,J) = 0. _d 0             a_QSWbyATM_cover(I,J) = 0. _d 0
1228            ENDDO            ENDDO
1229           ENDDO           ENDDO
1230    #endif
1231          ENDIF          ENDIF
1232  #endif  #endif
1233    
# Line 942  c available turbulent flux Line 1272  c available turbulent flux
1272            a_QbyOCN(i,j) =            a_QbyOCN(i,j) =
1273       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
1274            r_QbyOCN(i,j) = a_QbyOCN(i,j)            r_QbyOCN(i,j) = a_QbyOCN(i,j)
1275    ctm
1276              if (i.eq.20 .and. j.eq.20) then
1277                print *, HeatCapacity_Cp
1278                print *, rhoConst
1279                print *, recip_QI
1280                print *, theta(20,20,kSurface,bi,bj)
1281                print *, tempFrz
1282                print *, SEAICE_deltaTtherm
1283                print *, maskC(20,20,kSurface,bi,bj)
1284                print *, tmpscal2
1285                print *, a_QbyOCN(20,20)
1286              endif
1287    ctm
1288           ENDDO           ENDDO
1289          ENDDO          ENDDO
1290    
# Line 961  C ====================================== Line 1304  C ======================================
1304  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1305  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1307    #ifdef SEAICE_ITD
1308            DO K=1,nITD
1309    #endif
1310          DO J=1,sNy          DO J=1,sNy
1311           DO I=1,sNx           DO I=1,sNx
1312  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1313            tmpscal2 =            tmpscal2 =
1314    #ifdef SEAICE_ITD
1315         &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)
1316         &             *SNOW2ICE),ZERO)
1317    C         accumulate change over ITD categories
1318              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1319         &                                                       *ICE2SNOW
1320              HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2
1321         &                                                       *ICE2SNOW
1322              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1323    C         keep total up to date, too
1324              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1325    #else
1326       &     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)
1327            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1328            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW
1329            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1330    #endif
1331           ENDDO           ENDDO
1332          ENDDO          ENDDO
1333  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 979  CADJ STORE r_FWbySublim    = comlev1_bib Line 1338  CADJ STORE r_FWbySublim    = comlev1_bib
1338           DO I=1,sNx           DO I=1,sNx
1339  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1340            tmpscal2 =            tmpscal2 =
1341    #ifdef SEAICE_ITD
1342         &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)
1343    C         accumulate change over ITD categories
1344              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1345              HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2
1346              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1347    C         keep total up to date, too
1348              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1349    #else
1350       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1351            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
1352            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2
1353            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1354    #endif
1355           ENDDO           ENDDO
1356          ENDDO          ENDDO
1357          DO J=1,sNy          DO J=1,sNy
1358           DO I=1,sNx           DO I=1,sNx
1359  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.
1360  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
1361  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).
1362    #ifdef SEAICE_ITD
1363              a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
1364         &                              - r_FWbySublimMult(I,J,K)
1365              r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)
1366         &                              - r_FWbySublimMult(I,J,K)
1367             ENDDO
1368            ENDDO
1369    C       end K loop
1370            ENDDO
1371    C       then update totals      
1372            DO J=1,sNy
1373             DO I=1,sNx
1374    #endif
1375            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)
1376            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)
1377           ENDDO           ENDDO
1378          ENDDO          ENDDO
1379    c ToM<<< debug seaice_growth
1380            WRITE(msgBuf,'(A,7F6.2)')
1381         &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1382         &     HEFFITD(20,20,:,bi,bj)
1383            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1384         &    SQUEEZE_RIGHT , myThid)
1385    c ToM>>>
1386    
1387  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1388  C ===========================================================  C ===========================================================
# Line 1003  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1392  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1392  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1393  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1394    
1395    #ifdef SEAICE_ITD
1396            DO K=1,nITD
1397             DO J=1,sNy
1398              DO I=1,sNx
1399    C          ice growth/melt due to ocean heat is equally distributed under the ice
1400    C          and hence weighted by fractional area of each thickness category
1401               tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
1402         &                               -HEFFITD(I,J,K,bi,bj))
1403               d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
1404               r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1
1405               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj)  + tmpscal1
1406    #ifdef ALLOW_SITRACER
1407               SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1408         &                           + HEFFITD(I,J,K,bi,bj)
1409    #endif
1410              ENDDO
1411             ENDDO
1412            ENDDO
1413    c ToM<<< debug seaice_growth
1414            WRITE(msgBuf,'(A,7F9.6)')
1415         &    ' SEAICE_GROWTH: d_HEFFbyOCNonICE w/ITD: ',
1416         &     d_HEFFbyOCNonICE(20,20)
1417            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1418         &    SQUEEZE_RIGHT , myThid)
1419    c ToM>>>
1420    #else /* SEAICE_ITD */
1421          DO J=1,sNy          DO J=1,sNy
1422           DO I=1,sNx           DO I=1,sNx
1423            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 1428  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1428  #endif  #endif
1429           ENDDO           ENDDO
1430          ENDDO          ENDDO
1431    c ToM<<< debug seaice_growth
1432            WRITE(msgBuf,'(A,7F9.6)')
1433         &    ' SEAICE_GROWTH: d_HEFFbyOCNonICE w/o ITD: ',
1434         &     d_HEFFbyOCNonICE(20,20)
1435            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1436         &    SQUEEZE_RIGHT , myThid)
1437    c ToM>>>
1438    #endif /* SEAICE_ITD */
1439    c ToM<<< debug seaice_growth
1440            WRITE(msgBuf,'(A,7F6.2)')
1441         &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1442         &     HEFFITD(20,20,:,bi,bj)
1443            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1444         &    SQUEEZE_RIGHT , myThid)
1445    c ToM>>>
1446    
1447  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1448  C ==================================================================  C ==================================================================
# Line 1022  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1452  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1452  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1453  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1454    
1455    #ifdef SEAICE_ITD
1456            DO K=1,nITD
1457             DO J=1,sNy
1458              DO I=1,sNx
1459    C     Convert to standard units (meters of ice) rather than to meters
1460    C     of snow. This appears to be more robust.
1461              tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)
1462         &                                           *SNOW2ICE)
1463              tmpscal2=MIN(tmpscal1,0. _d 0)
1464    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1465    Cgf no additional dependency through snow
1466              IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1467    #endif
1468              d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW
1469              HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW
1470              r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
1471    C         keep the total up to date, too
1472              r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1473              ENDDO
1474             ENDDO
1475            ENDDO
1476    #else /* SEAICE_ITD */
1477          DO J=1,sNy          DO J=1,sNy
1478           DO I=1,sNx           DO I=1,sNx
1479  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 1489  Cgf no additional dependency through sno
1489            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1490           ENDDO           ENDDO
1491          ENDDO          ENDDO
1492    #endif /* SEAICE_ITD */
1493    c ToM<<< debug seaice_growth
1494            WRITE(msgBuf,'(A,7F6.2)')
1495         &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1496         &     HEFFITD(20,20,:,bi,bj)
1497            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1498         &    SQUEEZE_RIGHT , myThid)
1499    c ToM>>>
1500    
1501  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1502  C ====================================================  C ====================================================
# Line 1051  Cgf where all experiments start in Janua Line 1511  Cgf where all experiments start in Janua
1511  Cgf the v1.81=>v1.82 revision would change results in  Cgf the v1.81=>v1.82 revision would change results in
1512  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1513    
1514    #ifdef SEAICE_ITD
1515            DO K=1,nITD
1516             DO J=1,sNy
1517              DO I=1,sNx
1518    #ifdef SEAICE_GROWTH_LEGACY
1519               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
1520    #else
1521               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
1522    c         Limit ice growth by potential melt by ocean
1523         &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))
1524    #endif /* SEAICE_GROWTH_LEGACY */
1525               d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1526         &                                 + tmpscal2
1527               d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1528         &                                 + tmpscal2
1529               r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
1530         &                                 - tmpscal2
1531               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2
1532    
1533    #ifdef ALLOW_SITRACER
1534               SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1535         &                           + HEFFITD(I,J,K,bi,bj)
1536    #endif
1537              ENDDO
1538             ENDDO
1539            ENDDO
1540    #else /* SEAICE_ITD */
1541          DO J=1,sNy          DO J=1,sNy
1542           DO I=1,sNx           DO I=1,sNx
1543    
# Line 1070  c         Limit ice growth by potential Line 1557  c         Limit ice growth by potential
1557  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1558            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1559  #endif  #endif
1560           ENDDO           ENDDO
1561          ENDDO          ENDDO
1562    #endif /* SEAICE_ITD */
1563    c ToM<<< debug seaice_growth
1564            WRITE(msgBuf,'(A,7F6.2)')
1565         &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1566         &     HEFFITD(20,20,:,bi,bj)
1567            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1568         &    SQUEEZE_RIGHT , myThid)
1569    c ToM>>>
1570    
1571  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1572  C depending on atmospheric conditions.  C depending on atmospheric conditions.
# Line 1098  C           add precip to the fresh wate Line 1593  C           add precip to the fresh wate
1593       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1594              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
1595            ENDIF            ENDIF
1596             ENDDO
1597            ENDDO
1598    #ifdef SEAICE_ITD
1599            DO K=1,nITD
1600             DO J=1,sNy
1601              DO I=1,sNx
1602               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)
1603         &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
1604              ENDDO
1605             ENDDO
1606            ENDDO
1607    #else
1608            DO J=1,sNy
1609             DO I=1,sNx
1610            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)
1611           ENDDO           ENDDO
1612          ENDDO          ENDDO
1613    #endif
1614  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1615  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1616  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1617  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1618    c ToM<<< debug seaice_growth
1619            WRITE(msgBuf,'(A,7F6.2)')
1620         &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1621         &     HEFFITD(20,20,:,bi,bj)
1622            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1623         &    SQUEEZE_RIGHT , myThid)
1624    c ToM>>>
1625    
1626  C compute snow melt due to heat available from ocean.  C compute snow melt due to heat available from ocean.
1627  C =================================================================  C =================================================================
# Line 1116  Cph( very sensitive bit here by JZ Line 1633  Cph( very sensitive bit here by JZ
1633  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1634  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1635  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1636    
1637    #ifdef SEAICE_ITD
1638            DO K=1,nITD
1639             DO J=1,sNy
1640              DO I=1,sNx
1641               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
1642         &                  -HSNOWITD(I,J,K,bi,bj))
1643               tmpscal2=MIN(tmpscal1,0. _d 0)
1644    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1645    Cgf no additional dependency through snow
1646               if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1647    #endif
1648               d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1649               r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1650               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
1651              ENDDO
1652             ENDDO
1653            ENDDO
1654    #else /* SEAICE_ITD */
1655          DO J=1,sNy          DO J=1,sNy
1656           DO I=1,sNx           DO I=1,sNx
1657            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 1666  Cgf no additional dependency through sno
1666            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)
1667           ENDDO           ENDDO
1668          ENDDO          ENDDO
1669    #endif /* SEAICE_ITD */
1670  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1671  Cph)  Cph)
1672    c ToM<<< debug seaice_growth
1673            WRITE(msgBuf,'(A,7F6.2)')
1674         &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1675         &     HEFFITD(20,20,:,bi,bj)
1676            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1677         &    SQUEEZE_RIGHT , myThid)
1678    c ToM>>>
1679    
1680  C gain of new ice over open water  C gain of new ice over open water
1681  C ===============================  C ===============================
# Line 1159  C           or 0. otherwise (no melting Line 1703  C           or 0. otherwise (no melting
1703            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1704            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1705            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1706    #ifdef SEAICE_ITD
1707    C         open water area fraction
1708              tmpscal0 = ONE-AREApreTH(I,J)
1709    C         determine thickness of new ice
1710    C         considering the entire open water area to refreeze
1711              tmpscal1 = tmpscal3/tmpscal0
1712    C         then add new ice volume to appropriate thickness category
1713              DO K=1,nITD
1714               IF (tmpscal1.LT.Hlimit(K)) THEN
1715                HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
1716                tmpscal3=ZERO
1717    C not sure if AREAITD should be changed here since AREA is incremented
1718    C   in PART 4 below in non-itd code
1719    C in this scenario all open water is covered by new ice instantaneously,
1720    C   i.e. no delayed lead closing is concidered such as is associated with
1721    C   Hibler's h_0 parameter
1722                AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
1723         &                           + tmpscal0
1724                tmpscal0=ZERO
1725               ENDIF
1726              ENDDO
1727    #else
1728            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1729    #endif
1730           ENDDO           ENDDO
1731          ENDDO          ENDDO
1732    
1733  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1734    #ifdef SEAICE_ITD
1735            DO K=1,nITD
1736             DO J=1,sNy
1737              DO I=1,sNx
1738    c needs to be here to allow use also with LEGACY branch
1739               SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1740         &                           + HEFFITD(I,J,K,bi,bj)
1741              ENDDO
1742             ENDDO
1743            ENDDO
1744    #else
1745          DO J=1,sNy          DO J=1,sNy
1746           DO I=1,sNx           DO I=1,sNx
1747  c needs to be here to allow use also with LEGACY branch  c needs to be here to allow use also with LEGACY branch
1748            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1749           ENDDO           ENDDO
1750          ENDDO          ENDDO
1751    #endif
1752  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1753    c ToM<<< debug seaice_growth
1754            WRITE(msgBuf,'(A,7F6.2)')
1755         &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1756         &     HEFFITD(20,20,:,bi,bj)
1757            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1758         &    SQUEEZE_RIGHT , myThid)
1759    c ToM>>>
1760    
1761  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1762  C =================================  C =================================
# Line 1182  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1768  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1768  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1769  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1770          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1771    #ifdef SEAICE_ITD
1772             DO K=1,nITD
1773              DO J=1,sNy
1774               DO I=1,sNx
1775                tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow
1776         &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1777                tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
1778                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
1779                HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1
1780                HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
1781         &                            * ICE2SNOW
1782               ENDDO
1783              ENDDO
1784             ENDDO
1785    #else
1786           DO J=1,sNy           DO J=1,sNy
1787            DO I=1,sNx            DO I=1,sNx
1788             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 1792  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1792             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)
1793             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1794       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
1795            ENDDO            ENDDO
1796           ENDDO           ENDDO
1797    #endif
1798          ENDIF          ENDIF
1799  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1800    c ToM<<< debug seaice_growth
1801            WRITE(msgBuf,'(A,7F6.2)')
1802         &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1803         &     HEFFITD(20,20,:,bi,bj)
1804            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1805         &    SQUEEZE_RIGHT , myThid)
1806    c ToM>>>
1807    
1808  C ===================================================================  C ===================================================================
1809  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 1829  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
1829  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1830  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1831    
1832    #ifdef SEAICE_ITD
1833    C       replaces Hibler '79 scheme and lead closing parameter
1834    C       because ITD accounts explicitly for lead openings and
1835    C       different melt rates due to varying ice thickness
1836    C
1837    C       only consider ice area loss due to total ice thickness loss
1838    C       ice area gain due to freezing of open water as handled above
1839    C       under "gain of new ice over open water"
1840    C
1841    C       does not account for lateral melt of ice floes
1842    C
1843    C        AREAITD is incremented in section "gain of new ice over open water" above
1844    C
1845            DO K=1,nITD
1846             DO J=1,sNy
1847              DO I=1,sNx
1848               IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
1849                AREAITD(I,J,K,bi,bj)=ZERO
1850               ENDIF
1851    #ifdef ALLOW_SITRACER
1852               SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1853         &                           + AREAITD(I,J,K,bi,bj)
1854    #endif /* ALLOW_SITRACER */
1855              ENDDO
1856             ENDDO
1857            ENDDO
1858    #else /* SEAICE_ITD */
1859          DO J=1,sNy          DO J=1,sNy
1860           DO I=1,sNx           DO I=1,sNx
1861    
# Line 1289  C apply tendency Line 1925  C apply tendency
1925  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
1926           ENDDO           ENDDO
1927          ENDDO          ENDDO
1928    #endif /* SEAICE_ITD */
1929    
1930  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1931  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
1932          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
1933    #ifdef SEAICE_ITD
1934             DO K=1,nITD
1935              DO J=1,sNy
1936               DO I=1,sNx
1937                AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *
1938         &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )
1939               ENDDO
1940              ENDDO
1941             ENDDO
1942    #else
1943           DO J=1,sNy           DO J=1,sNy
1944            DO I=1,sNx            DO I=1,sNx
1945  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 1947  C            AREA(I,J,bi,bj) = 0.1 _d 0
1947       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1948            ENDDO            ENDDO
1949           ENDDO           ENDDO
1950    #endif
1951          ENDIF          ENDIF
1952  #endif  #endif
1953    #ifdef SEAICE_ITD
1954    C check categories for consistency with limits after growth/melt
1955            CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
1956    C finally update total AREA, HEFF, HSNOW
1957            CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
1958    #endif
1959    
1960  C ===================================================================  C ===================================================================
1961  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.3
  ViewVC Help
Powered by ViewVC 1.1.22