/[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.5 by torge, Tue Oct 2 16:47:41 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 138  c     The change of mean ice thickness d Line 142  c     The change of mean ice thickness d
142  #endif  #endif
143    
144  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
145  c     sea ice dyhnamics  c     sea ice dynamics
146        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
147    
148  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 173  C                     as EVAP (positive Line 177  C                     as EVAP (positive
177  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
178  C     The latent heat flux which will sublimate all snow and ice  C     The latent heat flux which will sublimate all snow and ice
179  C     over one time step  C     over one time step
180        _RL latentHeatFluxMax   (1:sNx,1:sNy)  #ifdef SEAICE_ITD
181        _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)        _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
182    #else
183          _RL latentHeatFluxMax   (1:sNx,1:sNy)
184    #endif
185  #endif  #endif
186    
187  C     actual ice thickness (with upper and lower limit)  C     actual ice thickness (with upper and lower limit)
# Line 192  C     AREA_PRE :: hold sea-ice fraction Line 199  C     AREA_PRE :: hold sea-ice fraction
199        _RL AREApreTH           (1:sNx,1:sNy)        _RL AREApreTH           (1:sNx,1:sNy)
200        _RL HEFFpreTH           (1:sNx,1:sNy)        _RL HEFFpreTH           (1:sNx,1:sNy)
201        _RL HSNWpreTH           (1:sNx,1:sNy)        _RL HSNWpreTH           (1:sNx,1:sNy)
202    #ifdef SEAICE_ITD
203          _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
204          _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
205          _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
206          _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
207          _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)
208    #endif
209    
210  C     wind speed  C     wind speed
211        _RL UG                  (1:sNx,1:sNy)        _RL UG                  (1:sNx,1:sNy)
# Line 221  C temporary variables available for the Line 235  C temporary variables available for the
235        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
236        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
237        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
238    #ifdef SEAICE_ITD
239          _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
240          _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
241    #endif
242        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
243        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
244        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
245    #ifdef SEAICE_ITD
246          _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
247          _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
248    #endif
249  C     Helper variables: reciprocal of some constants  C     Helper variables: reciprocal of some constants
250        _RL recip_multDim        _RL recip_multDim
251        _RL recip_deltaTtherm        _RL recip_deltaTtherm
# Line 259  C ====================================== Line 281  C ======================================
281        ENDIF        ENDIF
282    
283  C     avoid unnecessary divisions in loops  C     avoid unnecessary divisions in loops
284    #ifdef SEAICE_ITD
285    CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
286    #endif
287        recip_multDim        = SEAICE_multDim        recip_multDim        = SEAICE_multDim
288        recip_multDim        = ONE / recip_multDim        recip_multDim        = ONE / recip_multDim
289  C     above/below: double/single precision calculation of recip_multDim  C     above/below: double/single precision calculation of recip_multDim
# Line 354  C ===================== Line 379  C =====================
379            d_HEFFbySublim(I,J)        = 0.0 _d 0            d_HEFFbySublim(I,J)        = 0.0 _d 0
380            d_HSNWbySublim(I,J)        = 0.0 _d 0            d_HSNWbySublim(I,J)        = 0.0 _d 0
381  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
382    #ifdef SEAICE_ITD
383              DO IT=1,SEAICE_multDim
384               latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
385              ENDDO
386    #else
387            latentHeatFluxMax(I,J)     = 0.0 _d 0            latentHeatFluxMax(I,J)     = 0.0 _d 0
388  #endif  #endif
389    #endif
390  c  c
391            d_HFRWbyRAIN(I,J)          = 0.0 _d 0            d_HFRWbyRAIN(I,J)          = 0.0 _d 0
392    
# Line 370  c Line 401  c
401              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
402              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
403              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
404  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_ITD
405              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
406                r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
407  #endif  #endif
408            ENDDO            ENDDO
409           ENDDO           ENDDO
# Line 406  C ====================================== Line 438  C ======================================
438  #endif  #endif
439           ENDDO           ENDDO
440          ENDDO          ENDDO
441    #ifdef SEAICE_ITD
442            DO IT=1,nITD
443             DO J=1,sNy
444              DO I=1,sNx
445               HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
446               HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
447               AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
448              ENDDO
449             ENDDO
450            ENDDO
451    #endif
452    
453  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
454    
# Line 433  C         d_HEFFbyRLX(i,j) = 0. _d 0 Line 476  C         d_HEFFbyRLX(i,j) = 0. _d 0
476  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)
477              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
478             ENDIF             ENDIF
479    #ifdef SEAICE_ITD
480                AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
481         &                           +  d_AREAbyRLX(i,j)
482                HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
483         &                           +  d_HEFFbyRLX(i,j)
484    #endif
485              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)
486              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)
487           ENDDO           ENDDO
# Line 449  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 498  CADJ STORE area(:,:,bi,bj) = comlev1_bib
498  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
499          DO J=1,sNy          DO J=1,sNy
500           DO I=1,sNx           DO I=1,sNx
501    #ifdef SEAICE_ITD
502              DO IT=1,nITD
503               tmpscal2=0. _d 0
504               tmpscal3=0. _d 0
505               tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
506               HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
507               d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
508               tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
509               HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
510               d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
511               AREAITD(I,J,IT,bi,bj)=MAX(-AREAITD(I,J,IT,bi,bj),0. _d 0)
512              ENDDO
513    CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
514    C         by calling SEAICE_ITD_SUM
515    #else
516            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)
517            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
518            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)
519            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
520            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)
521    #endif
522           ENDDO           ENDDO
523          ENDDO          ENDDO
524    
# Line 464  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 529  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
529  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
530          DO J=1,sNy          DO J=1,sNy
531           DO I=1,sNx           DO I=1,sNx
532            tmpscal2=0. _d 0  #ifdef SEAICE_ITD
533            tmpscal3=0. _d 0            DO IT=1,nITD
534    #endif
535              tmpscal2=0. _d 0
536              tmpscal3=0. _d 0
537    #ifdef SEAICE_ITD
538               IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
539                tmpscal2=-HEFFITD(I,J,IT,bi,bj)
540                tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
541                TICES(I,J,IT,bi,bj)=celsius2K
542    CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
543               ENDIF
544               HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
545               HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
546    #else
547            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
548              tmpscal2=-HEFF(I,J,bi,bj)             tmpscal2=-HEFF(I,J,bi,bj)
549              tmpscal3=-HSNOW(I,J,bi,bj)             tmpscal3=-HSNOW(I,J,bi,bj)
550              TICE(I,J,bi,bj)=celsius2K             TICE(I,J,bi,bj)=celsius2K
551              DO IT=1,SEAICE_multDim             DO IT=1,SEAICE_multDim
552                TICES(I,J,IT,bi,bj)=celsius2K              TICES(I,J,IT,bi,bj)=celsius2K
553              ENDDO             ENDDO
554            ENDIF            ENDIF
555            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  
556            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
557    #endif
558              d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
559            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
560    #ifdef SEAICE_ITD
561              ENDDO
562    #endif
563           ENDDO           ENDDO
564          ENDDO          ENDDO
565    
# Line 489  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 571  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
571  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
572          DO J=1,sNy          DO J=1,sNy
573           DO I=1,sNx           DO I=1,sNx
574    #ifdef SEAICE_ITD
575              DO IT=1,nITD
576               IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
577         &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
578         &      AREAITD(I,J,IT,bi,bj)=0. _d 0
579              ENDDO
580    #else
581            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
582       &        (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
583    #endif
584           ENDDO           ENDDO
585          ENDDO          ENDDO
586    
# Line 502  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 592  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
592  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
593          DO J=1,sNy          DO J=1,sNy
594           DO I=1,sNx           DO I=1,sNx
595    #ifdef SEAICE_ITD
596              DO IT=1,nITD
597               IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
598         &         (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
599    CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1
600    C           hence use SEAICE_area_floor devided by nITD
601    C           (or install a warning in e.g. seaice_readparms.F)
602                AREAITD(I,J,IT,bi,bj)=
603         &         MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
604               ENDIF
605              ENDDO
606    #else
607            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
608             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)
609            ENDIF            ENDIF
610    #endif
611           ENDDO           ENDDO
612          ENDDO          ENDDO
613  #endif /* DISABLE_AREA_FLOOR */  #endif /* DISABLE_AREA_FLOOR */
614    
615  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:
616    
617    CToM   for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
618    C      which is called at end of PART 1 below
619    #ifndef SEAICE_ITD
620  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
621  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
622  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 525  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 631  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
631            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)
632           ENDDO           ENDDO
633          ENDDO          ENDDO
634    #endif /* notSEAICE_ITD */
635    
636    #ifdef SEAICE_ITD
637    CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
638            DO J=1,sNy
639             DO I=1,sNx
640    C    TICES was changed above (item 1.25), now update TICE as ice volume
641    C     weighted average of TICES
642    C    also compute total of AREAITD (needed for finishing item 2.5, see below)
643              tmpscal1 = 0. _d 0
644              tmpscal2 = 0. _d 0
645              tmpscal3 = 0. _d 0
646              DO IT=1,nITD
647               tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
648               tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
649               tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
650              ENDDO
651              TICE(I,J,bi,bj)=tmpscal1/tmpscal2
652    C    lines of item 2.5 that were omitted:
653    C    in 2.5 these lines are executed before "ridging" is applied to AREA
654    C    hence we execute them here before SEAICE_ITD_REDIST is called
655    C    although this means that AREA has not been completely regularized
656    #ifdef ALLOW_DIAGNOSTICS
657              DIAGarrayA(I,J) = tmpscal3
658    #endif
659    #ifdef ALLOW_SITRACER
660              SItrAREA(I,J,bi,bj,1)=tmpscal3
661    #endif
662             ENDDO
663            ENDDO
664    
665    CToM finally make sure that all categories meet their thickness limits
666    C    which includes ridging as in item 2.5
667    C    and update AREA, HEFF, and HSNOW
668            CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
669            CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
670    
671    #endif
672  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
673    C        end SEAICEadjMODE.EQ.0 statement:
674          ENDIF          ENDIF
675  #endif  #endif
676    
# Line 547  C 3) store regularized values of heff, h Line 691  C 3) store regularized values of heff, h
691  #endif  #endif
692           ENDDO           ENDDO
693          ENDDO          ENDDO
694    #ifdef SEAICE_ITD
695            DO IT=1,nITD
696             DO J=1,sNy
697              DO I=1,sNx
698               HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
699               HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
700               AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
701    
702    C memorize areal and volume fraction of each ITD category
703               IF (AREA(I,J,bi,bj).GT.0) THEN
704                areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
705               ELSE
706                areaFracFactor(I,J,IT)=ZERO
707               ENDIF
708               IF (HEFF(I,J,bi,bj).GT.0) THEN
709                heffFracFactor(I,J,IT)=HEFFITD(I,J,IT,bi,bj)/HEFF(I,J,bi,bj)
710               ELSE
711                heffFracFactor(I,J,IT)=ZERO
712               ENDIF
713              ENDDO
714             ENDDO
715            ENDDO
716    C prepare SItrHEFF to be computed as cumulative sum
717            DO iTr=2,5
718             DO J=1,sNy
719              DO I=1,sNx
720               SItrHEFF(I,J,bi,bj,iTr)=ZERO
721              ENDDO
722             ENDDO
723            ENDDO
724    C prepare SItrAREA to be computed as cumulative sum
725            DO J=1,sNy
726             DO I=1,sNx
727              SItrAREA(I,J,bi,bj,3)=ZERO
728             ENDDO
729            ENDDO
730    #endif
731    
732  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
733  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
# Line 601  Cgf no additional dependency of air-sea Line 782  Cgf no additional dependency of air-sea
782             AREApreTH(I,J) = 0. _d 0             AREApreTH(I,J) = 0. _d 0
783            ENDDO            ENDDO
784           ENDDO           ENDDO
785    #ifdef SEAICE_ITD
786             DO IT=1,nITD
787              DO J=1,sNy
788               DO I=1,sNx
789                HEFFITDpreTH(I,J,IT) = 0. _d 0
790                HSNWITDpreTH(I,J,IT) = 0. _d 0
791                AREAITDpreTH(I,J,IT) = 0. _d 0
792               ENDDO
793              ENDDO
794             ENDDO
795    #endif
796          ENDIF          ENDIF
797  #endif  #endif
798    
# Line 620  CADJ STORE AREApreTH = comlev1_bibj, key Line 812  CADJ STORE AREApreTH = comlev1_bibj, key
812  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
813  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
814  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
815    #ifdef SEAICE_ITD
816            DO IT=1,nITD
817    #endif
818          DO J=1,sNy          DO J=1,sNy
819           DO I=1,sNx           DO I=1,sNx
820    
821    #ifdef SEAICE_ITD
822              IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
823    #ifdef SEAICE_GROWTH_LEGACY
824                tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
825         &                              AREAITDpreTH(I,J,IT))
826                hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
827                tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
828                heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
829    #else /* SEAICE_GROWTH_LEGACY */
830    cif        regularize AREA with SEAICE_area_reg
831               tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
832         &                     + area_reg_sq)
833    cif        heffActual calculated with the regularized AREA
834               tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
835    cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
836               heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
837         &                                  + hice_reg_sq)
838    cif        hsnowActual calculated with the regularized AREA
839               hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
840    #endif /* SEAICE_GROWTH_LEGACY */
841    cif        regularize the inverse of heffActual by hice_reg
842               recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
843         &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
844         &                      + hice_reg_sq)
845    cif       Do not regularize when HEFFpreTH = 0
846              ELSE
847                heffActualMult(I,J,IT) = ZERO
848                hsnowActualMult(I,J,IT) = ZERO
849                recip_heffActualMult(I,J,IT)  = ZERO
850              ENDIF
851    #else /* SEAICE_ITD */
852            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
853  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
854              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 874  cif       Do not regularize when HEFFpre
874              hsnowActual(I,J) = ZERO              hsnowActual(I,J) = ZERO
875              recip_heffActual(I,J)  = ZERO              recip_heffActual(I,J)  = ZERO
876            ENDIF            ENDIF
877    #endif /* SEAICE_ITD */
878    
879           ENDDO           ENDDO
880          ENDDO          ENDDO
881    #ifdef SEAICE_ITD
882            ENDDO
883    #endif
884    
885  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
886          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 891  cif       Do not regularize when HEFFpre
891  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
892  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
893  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
894    #ifdef SEAICE_ITD
895            DO IT=1,nITD
896    #endif
897          DO J=1,sNy          DO J=1,sNy
898           DO I=1,sNx           DO I=1,sNx
899  c        The latent heat flux over the sea ice which  c        The latent heat flux over the sea ice which
900  c        will sublimate all of the snow and ice over one time  c        will sublimate all of the snow and ice over one time
901  c        step (W/m^2)  c        step (W/m^2)
902    #ifdef SEAICE_ITD
903              IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
904                latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
905         &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
906         &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
907         &         /AREAITDpreTH(I,J,IT)
908              ELSE
909                latentHeatFluxMaxMult(I,J,IT) = ZERO
910              ENDIF
911    #else
912            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
913              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
914       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +
# Line 673  c        step (W/m^2) Line 916  c        step (W/m^2)
916            ELSE            ELSE
917              latentHeatFluxMax(I,J) = ZERO              latentHeatFluxMax(I,J) = ZERO
918            ENDIF            ENDIF
919    #endif
920           ENDDO           ENDDO
921          ENDDO          ENDDO
922    #ifdef SEAICE_ITD
923            ENDDO
924    #endif
925  #endif /* SEAICE_CAP_SUBLIM */  #endif /* SEAICE_CAP_SUBLIM */
926    
927  C ===================================================================  C ===================================================================
# Line 746  CADJ &                       key = iicek Line 993  CADJ &                       key = iicek
993  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
994    
995  C--   Start loop over multi-categories  C--   Start loop over multi-categories
996    #ifdef SEAICE_ITD
997    CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
998    #endif
999          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1000  c        homogeneous distribution between 0 and 2 x heffActual  c        homogeneous distribution between 0 and 2 x heffActual
1001    #ifndef SEAICE_ITD
1002           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
1003    #endif
1004           DO J=1,sNy           DO J=1,sNy
1005            DO I=1,sNx            DO I=1,sNx
1006    #ifndef SEAICE_ITD
1007    CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
1008    C    (instead of heffActual and latentHeatFluxMax)
1009             heffActualMult(I,J,IT)= heffActual(I,J)*pFac             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1010  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1011             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1012  #endif  #endif
1013    #endif
1014             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1015             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1016             TICE(I,J,bi,bj) = ZERO             TICE(I,J,bi,bj) = ZERO
# Line 780  CADJ &     comlev1_bibj, key = iicekey, Line 1036  CADJ &     comlev1_bibj, key = iicekey,
1036    
1037          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1038           CALL SEAICE_SOLVE4TEMP(           CALL SEAICE_SOLVE4TEMP(
1039    #ifdef SEAICE_ITD
1040         I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1041    #else
1042       I        UG, heffActualMult(1,1,IT), hsnowActual,       I        UG, heffActualMult(1,1,IT), hsnowActual,
1043    #endif
1044  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1045       I        latentHeatFluxMaxMult(1,1,IT),       I        latentHeatFluxMaxMult(1,1,IT),
1046  #endif  #endif
# Line 809  CADJ &     comlev1_bibj, key = iicekey, Line 1069  CADJ &     comlev1_bibj, key = iicekey,
1069           DO J=1,sNy           DO J=1,sNy
1070            DO I=1,sNx            DO I=1,sNx
1071  C     update TICE & TICES  C     update TICE & TICES
1072    #ifdef SEAICE_ITD
1073    C     calculate area weighted mean
1074    C     (although the ice's temperature relates to its energy content
1075    C      and hence should be averaged weighted by ice volume [heffFracFactor],
1076    C      the temperature here is a result of the fluxes through the ice surface
1077    C      computed individually for each single category in SEAICE_SOLVE4TEMP
1078    C      and hence is averaged area weighted [areaFracFactor])
1079               TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1080         &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1081    #else
1082             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1083       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
1084    #endif
1085             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1086  C     average over categories  C     average over categories
1087    #ifdef SEAICE_ITD
1088    C     calculate area weighted mean
1089    C     (fluxes are per unit (ice surface) area and are thus area weighted)
1090               a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1091         &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1092               a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1093         &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1094               a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1095         &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1096    #else
1097             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1098       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1099             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1100       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1101             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1102       &          + a_FWbySublimMult(I,J,IT)*recip_multDim       &          + a_FWbySublimMult(I,J,IT)*recip_multDim
1103    #endif
1104            ENDDO            ENDDO
1105           ENDDO           ENDDO
1106          ENDDO          ENDDO
# Line 851  CADJ STORE a_FWbySublim    = comlev1_bib Line 1133  CADJ STORE a_FWbySublim    = comlev1_bib
1133  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1134    
1135  C switch heat fluxes from W/m2 to 'effective' ice meters  C switch heat fluxes from W/m2 to 'effective' ice meters
1136    #ifdef SEAICE_ITD
1137            DO IT=1,nITD
1138             DO J=1,sNy
1139              DO I=1,sNx
1140               a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
1141         &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1142               a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1143         &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1144    C and initialize r_QbyATM_cover
1145               r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
1146    C     Convert fresh water flux by sublimation to 'effective' ice meters.
1147    C     Negative sublimation is resublimation and will be added as snow.
1148    #ifdef SEAICE_DISABLE_SUBLIM
1149               a_FWbySublimMult(I,J,IT) = ZERO
1150    #endif
1151               a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
1152         &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
1153               r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
1154              ENDDO
1155             ENDDO
1156            ENDDO
1157            DO J=1,sNy
1158             DO I=1,sNx
1159              a_QbyATM_open(I,J)    = a_QbyATM_open(I,J)
1160         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1161              a_QSWbyATM_open(I,J)  = a_QSWbyATM_open(I,J)
1162         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1163    C and initialize r_QbyATM_open
1164              r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1165             ENDDO
1166            ENDDO
1167    #else /* SEAICE_ITD */
1168          DO J=1,sNy          DO J=1,sNy
1169           DO I=1,sNx           DO I=1,sNx
1170            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 1183  C     Negative sublimation is resublimat
1183  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1184  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
1185            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1186  #endif /* SEAICE_CAP_SUBLIM */  #endif
1187            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1188       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1189            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1190           ENDDO           ENDDO
1191          ENDDO          ENDDO
1192    #endif /* SEAICE_ITD */
1193    
1194  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1195  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 1206  CADJ STORE r_FWbySublim    = comlev1_bib
1206  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1207  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1208          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1209    #ifdef SEAICE_ITD
1210             DO IT=1,nITD
1211              DO J=1,sNy
1212               DO I=1,sNx
1213                a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1214                r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1215                a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
1216               ENDDO
1217              ENDDO
1218             ENDDO
1219    #else
1220           DO J=1,sNy           DO J=1,sNy
1221            DO I=1,sNx            DO I=1,sNx
1222             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 1224  Cgf no additional dependency through ice
1224             a_QSWbyATM_cover(I,J) = 0. _d 0             a_QSWbyATM_cover(I,J) = 0. _d 0
1225            ENDDO            ENDDO
1226           ENDDO           ENDDO
1227    #endif
1228          ENDIF          ENDIF
1229  #endif  #endif
1230    
# Line 961  C ====================================== Line 1288  C ======================================
1288  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1289  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1290  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1291    #ifdef SEAICE_ITD
1292            DO IT=1,nITD
1293    #endif
1294          DO J=1,sNy          DO J=1,sNy
1295           DO I=1,sNx           DO I=1,sNx
1296  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1297            tmpscal2 =            tmpscal2 =
1298    #ifdef SEAICE_ITD
1299         &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1300         &             *SNOW2ICE),ZERO)
1301    C         accumulate change over ITD categories
1302              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
1303         &                                                       *ICE2SNOW
1304              HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
1305         &                                                       *ICE2SNOW
1306              r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
1307    C         keep total up to date, too
1308              r_FWbySublim(I,J)       = r_FWbySublim(I,J)        - tmpscal2
1309    #else
1310       &     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)
1311            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1312            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW
1313            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1314    #endif
1315           ENDDO           ENDDO
1316          ENDDO          ENDDO
1317  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 979  CADJ STORE r_FWbySublim    = comlev1_bib Line 1322  CADJ STORE r_FWbySublim    = comlev1_bib
1322           DO I=1,sNx           DO I=1,sNx
1323  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1324            tmpscal2 =            tmpscal2 =
1325    #ifdef SEAICE_ITD
1326         &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1327    C         accumulate change over ITD categories
1328              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1329              HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
1330              r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
1331    C         keep total up to date, too
1332              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1333    #else
1334       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1335            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
1336            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2
1337            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1338    #endif
1339           ENDDO           ENDDO
1340          ENDDO          ENDDO
1341          DO J=1,sNy          DO J=1,sNy
1342           DO I=1,sNx           DO I=1,sNx
1343  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.
1344  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
1345  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).
1346    #ifdef SEAICE_ITD
1347              a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1348         &                               - r_FWbySublimMult(I,J,IT)
1349              r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1350         &                               - r_FWbySublimMult(I,J,IT)
1351             ENDDO
1352            ENDDO
1353    C       end IT loop
1354            ENDDO
1355    C       then update totals      
1356            DO J=1,sNy
1357             DO I=1,sNx
1358    #endif
1359            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)
1360            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)
1361           ENDDO           ENDDO
1362          ENDDO          ENDDO
1363    c ToM<<< debug seaice_growth
1364            WRITE(msgBuf,'(A,7F6.2)')
1365    #ifdef SEAICE_ITD
1366         &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1367         &     HEFFITD(20,20,:,bi,bj)
1368    #else
1369         &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1370         &     HEFF(20,20,bi,bj)
1371    #endif
1372            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1373         &    SQUEEZE_RIGHT , myThid)
1374    c ToM>>>
1375    
1376  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1377  C ===========================================================  C ===========================================================
# Line 1003  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1381  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1381  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1382  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1383    
1384    #ifdef SEAICE_ITD
1385            DO IT=1,nITD
1386             DO J=1,sNy
1387              DO I=1,sNx
1388    C          ice growth/melt due to ocean heat is equally distributed under the ice
1389    C          and hence weighted by fractional area of each thickness category
1390               tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1391         &                               -HEFFITD(I,J,IT,bi,bj))
1392               d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
1393               r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1
1394               HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1
1395    #ifdef ALLOW_SITRACER
1396               SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1397         &                           + HEFFITD(I,J,IT,bi,bj)
1398    #endif
1399              ENDDO
1400             ENDDO
1401            ENDDO
1402    #else /* SEAICE_ITD */
1403          DO J=1,sNy          DO J=1,sNy
1404           DO I=1,sNx           DO I=1,sNx
1405            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 1410  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1410  #endif  #endif
1411           ENDDO           ENDDO
1412          ENDDO          ENDDO
1413    #endif /* SEAICE_ITD */
1414    c ToM<<< debug seaice_growth
1415            WRITE(msgBuf,'(A,7F6.2)')
1416    #ifdef SEAICE_ITD
1417         &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1418         &     HEFFITD(20,20,:,bi,bj)
1419    #else
1420         &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1421         &     HEFF(20,20,bi,bj)
1422    #endif
1423            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1424         &    SQUEEZE_RIGHT , myThid)
1425    c ToM>>>
1426    
1427  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1428  C ==================================================================  C ==================================================================
# Line 1022  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1432  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1432  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1433  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1434    
1435    #ifdef SEAICE_ITD
1436            DO IT=1,nITD
1437             DO J=1,sNy
1438              DO I=1,sNx
1439    C     Convert to standard units (meters of ice) rather than to meters
1440    C     of snow. This appears to be more robust.
1441               tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1442         &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1443               tmpscal2=MIN(tmpscal1,0. _d 0)
1444    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1445    Cgf no additional dependency through snow
1446               IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1447    #endif
1448               d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1449         &                           + tmpscal2*ICE2SNOW
1450               HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
1451         &                           + tmpscal2*ICE2SNOW
1452               r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1453         &                           - tmpscal2
1454    C          keep the total up to date, too
1455               r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1456              ENDDO
1457             ENDDO
1458            ENDDO
1459    #else /* SEAICE_ITD */
1460          DO J=1,sNy          DO J=1,sNy
1461           DO I=1,sNx           DO I=1,sNx
1462  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 1472  Cgf no additional dependency through sno
1472            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1473           ENDDO           ENDDO
1474          ENDDO          ENDDO
1475    #endif /* SEAICE_ITD */
1476    c ToM<<< debug seaice_growth
1477            WRITE(msgBuf,'(A,7F6.2)')
1478    #ifdef SEAICE_ITD
1479         &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1480         &     HEFFITD(20,20,:,bi,bj)
1481    #else
1482         &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1483         &     HEFF(20,20,bi,bj)
1484    #endif
1485            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1486         &    SQUEEZE_RIGHT , myThid)
1487    c ToM>>>
1488    
1489  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1490  C ====================================================  C ====================================================
# Line 1051  Cgf where all experiments start in Janua Line 1499  Cgf where all experiments start in Janua
1499  Cgf the v1.81=>v1.82 revision would change results in  Cgf the v1.81=>v1.82 revision would change results in
1500  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1501    
1502    #ifdef SEAICE_ITD
1503            DO IT=1,nITD
1504             DO J=1,sNy
1505              DO I=1,sNx
1506    #ifdef SEAICE_GROWTH_LEGACY
1507               tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1508         &                     r_QbyATMmult_cover(I,J,IT))
1509    #else
1510               tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1511         &                     r_QbyATMmult_cover(I,J,IT)
1512    c         Limit ice growth by potential melt by ocean
1513         &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1514    #endif /* SEAICE_GROWTH_LEGACY */
1515               d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1516         &                                 + tmpscal2
1517               d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1518         &                                 + tmpscal2
1519               r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
1520         &                                 - tmpscal2
1521               HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
1522    
1523    #ifdef ALLOW_SITRACER
1524               SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1525         &                           + HEFFITD(I,J,IT,bi,bj)
1526    #endif
1527              ENDDO
1528             ENDDO
1529            ENDDO
1530    #else /* SEAICE_ITD */
1531          DO J=1,sNy          DO J=1,sNy
1532           DO I=1,sNx           DO I=1,sNx
1533    
# Line 1070  c         Limit ice growth by potential Line 1547  c         Limit ice growth by potential
1547  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1548            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1549  #endif  #endif
1550           ENDDO           ENDDO
1551          ENDDO          ENDDO
1552    #endif /* SEAICE_ITD */
1553    c ToM<<< debug seaice_growth
1554            WRITE(msgBuf,'(A,7F6.2)')
1555    #ifdef SEAICE_ITD
1556         &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1557         &     HEFFITD(20,20,:,bi,bj)
1558    #else
1559         &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1560         &     HEFF(20,20,bi,bj)
1561    #endif
1562            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1563         &    SQUEEZE_RIGHT , myThid)
1564    c ToM>>>
1565    
1566  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1567  C depending on atmospheric conditions.  C depending on atmospheric conditions.
# Line 1098  C           add precip to the fresh wate Line 1588  C           add precip to the fresh wate
1588       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1589              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
1590            ENDIF            ENDIF
1591             ENDDO
1592            ENDDO
1593    #ifdef SEAICE_ITD
1594            DO IT=1,nITD
1595             DO J=1,sNy
1596              DO I=1,sNx
1597               HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1598         &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1599              ENDDO
1600             ENDDO
1601            ENDDO
1602    #else
1603            DO J=1,sNy
1604             DO I=1,sNx
1605            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)
1606           ENDDO           ENDDO
1607          ENDDO          ENDDO
1608    #endif
1609  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1610  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1611  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1612  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1613    c ToM<<< debug seaice_growth
1614            WRITE(msgBuf,'(A,7F6.2)')
1615    #ifdef SEAICE_ITD
1616         &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1617         &     HEFFITD(20,20,:,bi,bj)
1618    #else
1619         &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1620         &     HEFF(20,20,bi,bj)
1621    #endif
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 IT=1,nITD
1639             DO J=1,sNy
1640              DO I=1,sNx
1641               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1642         &                  -HSNOWITD(I,J,IT,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,IT,bi,bj) = HSNOWITD(I,J,IT,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    #ifdef SEAICE_ITD
1675         &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1676         &     HEFFITD(20,20,:,bi,bj)
1677    #else
1678         &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1679         &     HEFF(20,20,bi,bj)
1680    #endif
1681            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1682         &    SQUEEZE_RIGHT , myThid)
1683    c ToM>>>
1684    
1685  C gain of new ice over open water  C gain of new ice over open water
1686  C ===============================  C ===============================
# Line 1159  C           or 0. otherwise (no melting Line 1708  C           or 0. otherwise (no melting
1708            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1709            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1710            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1711    #ifdef SEAICE_ITD
1712    C         open water area fraction
1713              tmpscal0 = ONE-AREApreTH(I,J)
1714    C         determine thickness of new ice
1715    C         considering the entire open water area to refreeze
1716              tmpscal1 = tmpscal3/tmpscal0
1717    C         then add new ice volume to appropriate thickness category
1718              DO IT=1,nITD
1719               IF (tmpscal1.LT.Hlimit(IT)) THEN
1720                HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
1721                tmpscal3=ZERO
1722    C not sure if AREAITD should be changed here since AREA is incremented
1723    C   in PART 4 below in non-itd code
1724    C in this scenario all open water is covered by new ice instantaneously,
1725    C   i.e. no delayed lead closing is concidered such as is associated with
1726    C   Hibler's h_0 parameter
1727                AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
1728         &                           + tmpscal0
1729                tmpscal0=ZERO
1730               ENDIF
1731              ENDDO
1732    #else
1733            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1734    #endif
1735           ENDDO           ENDDO
1736          ENDDO          ENDDO
1737    
1738  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1739    #ifdef SEAICE_ITD
1740            DO IT=1,nITD
1741             DO J=1,sNy
1742              DO I=1,sNx
1743    c needs to be here to allow use also with LEGACY branch
1744               SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1745         &                           + HEFFITD(I,J,IT,bi,bj)
1746              ENDDO
1747             ENDDO
1748            ENDDO
1749    #else
1750          DO J=1,sNy          DO J=1,sNy
1751           DO I=1,sNx           DO I=1,sNx
1752  c needs to be here to allow use also with LEGACY branch  c needs to be here to allow use also with LEGACY branch
1753            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1754           ENDDO           ENDDO
1755          ENDDO          ENDDO
1756    #endif
1757  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1758    c ToM<<< debug seaice_growth
1759            WRITE(msgBuf,'(A,7F6.2)')
1760    #ifdef SEAICE_ITD
1761         &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1762         &     HEFFITD(20,20,:,bi,bj)
1763    #else
1764         &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1765         &     HEFF(20,20,bi,bj)
1766    #endif
1767            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1768         &    SQUEEZE_RIGHT , myThid)
1769    c ToM>>>
1770    
1771  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1772  C =================================  C =================================
# Line 1182  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1778  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1778  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1779  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1780          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1781    #ifdef SEAICE_ITD
1782             DO IT=1,nITD
1783              DO J=1,sNy
1784               DO I=1,sNx
1785                tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
1786         &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
1787         &                                        *recip_rhoConst
1788                tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
1789                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
1790                HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1
1791                HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
1792         &                            * ICE2SNOW
1793               ENDDO
1794              ENDDO
1795             ENDDO
1796    #else
1797           DO J=1,sNy           DO J=1,sNy
1798            DO I=1,sNx            DO I=1,sNx
1799             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 1803  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1803             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)
1804             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1805       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
1806            ENDDO            ENDDO
1807           ENDDO           ENDDO
1808    #endif
1809          ENDIF          ENDIF
1810  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1811    c ToM<<< debug seaice_growth
1812            WRITE(msgBuf,'(A,7F6.2)')
1813    #ifdef SEAICE_ITD
1814         &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1815         &     HEFFITD(20,20,:,bi,bj)
1816    #else
1817         &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1818         &     HEFF(20,20,bi,bj)
1819    #endif
1820            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1821         &    SQUEEZE_RIGHT , myThid)
1822    c ToM>>>
1823    
1824  C ===================================================================  C ===================================================================
1825  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 1845  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
1845  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1846  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1847    
1848    #ifdef SEAICE_ITD
1849    C       replaces Hibler '79 scheme and lead closing parameter
1850    C       because ITD accounts explicitly for lead openings and
1851    C       different melt rates due to varying ice thickness
1852    C
1853    C       only consider ice area loss due to total ice thickness loss
1854    C       ice area gain due to freezing of open water as handled above
1855    C       under "gain of new ice over open water"
1856    C
1857    C       does not account for lateral melt of ice floes
1858    C
1859    C        AREAITD is incremented in section "gain of new ice over open water" above
1860    C
1861            DO IT=1,nITD
1862             DO J=1,sNy
1863              DO I=1,sNx
1864               IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1865                AREAITD(I,J,IT,bi,bj)=ZERO
1866               ENDIF
1867    #ifdef ALLOW_SITRACER
1868               SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1869         &                           + AREAITD(I,J,IT,bi,bj)
1870    #endif /* ALLOW_SITRACER */
1871              ENDDO
1872             ENDDO
1873            ENDDO
1874    #else /* SEAICE_ITD */
1875          DO J=1,sNy          DO J=1,sNy
1876           DO I=1,sNx           DO I=1,sNx
1877    
# Line 1289  C apply tendency Line 1941  C apply tendency
1941  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
1942           ENDDO           ENDDO
1943          ENDDO          ENDDO
1944    #endif /* SEAICE_ITD */
1945    
1946  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1947  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
1948          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
1949    #ifdef SEAICE_ITD
1950             DO IT=1,nITD
1951              DO J=1,sNy
1952               DO I=1,sNx
1953                AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
1954         &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
1955               ENDDO
1956              ENDDO
1957             ENDDO
1958    #else
1959           DO J=1,sNy           DO J=1,sNy
1960            DO I=1,sNx            DO I=1,sNx
1961  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 1963  C            AREA(I,J,bi,bj) = 0.1 _d 0
1963       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1964            ENDDO            ENDDO
1965           ENDDO           ENDDO
1966    #endif
1967          ENDIF          ENDIF
1968  #endif  #endif
1969    #ifdef SEAICE_ITD
1970    C check categories for consistency with limits after growth/melt
1971            CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
1972    C finally update total AREA, HEFF, HSNOW
1973            CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
1974    #endif
1975    
1976  C ===================================================================  C ===================================================================
1977  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.5
  ViewVC Help
Powered by ViewVC 1.1.22