/[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.2 by dimitri, Fri Apr 27 22:25:23 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 175  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 228  C temporary variables available for the Line 232  C temporary variables available for the
232  #endif  #endif
233    
234        INTEGER ilockey        INTEGER ilockey
235  CToM<<<        INTEGER it
236  C      INTEGER it  #ifdef SEAICE_ITD
237        INTEGER IT, K        INTEGER K
238  C>>>ToM  #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)
# Line 282  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 400  c Line 407  c
407              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
408              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
409  #endif  #endif
410  #ifdef SEAICE_CAP_SUBLIM  #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 443  C ====================================== Line 450  C ======================================
450             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
451             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
452             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
            IF (AREA(I,J,bi,bj).GT.0) THEN  
             areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)  
            ELSE  
             areaFracFactor(I,J,K)=ZERO  
            ENDIF  
            IF (HEFF(I,J,bi,bj).GT.0) THEN  
             heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)  
            ELSE  
             heffFracFactor(I,J,K)=ZERO  
            ENDIF  
453            ENDDO            ENDDO
454           ENDDO           ENDDO
455          ENDDO          ENDDO
# Line 508  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 505  CADJ STORE area(:,:,bi,bj) = comlev1_bib
505           DO I=1,sNx           DO I=1,sNx
506  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
507            DO K=1,nITD            DO K=1,nITD
508               tmpscal1=0. _d 0
509             tmpscal2=0. _d 0             tmpscal2=0. _d 0
510             tmpscal3=0. _d 0             tmpscal3=0. _d 0
            tmpscal4=0. _d 0  
511             tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)             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             HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
513             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
514             tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)             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             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
516             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
517             tmpscal4=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)             tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
518             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1
519             d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal4             d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1
520            ENDDO            ENDDO
521            AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyNEG(I,J)  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
522    C         by calling SEAICE_ITD_SUM
523  #else  #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)
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)
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)
 #endif  
527            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)
528            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)
529    #endif
530           ENDDO           ENDDO
531          ENDDO          ENDDO
532    
# Line 541  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 539  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
539           DO I=1,sNx           DO I=1,sNx
540  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
541            DO K=1,nITD            DO K=1,nITD
542             tmpscal2=0. _d 0  #endif
543             tmpscal3=0. _d 0            tmpscal2=0. _d 0
544              tmpscal3=0. _d 0
545    #ifdef SEAICE_ITD
546             IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN             IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
547              tmpscal2=-HEFFITD(I,J,K,bi,bj)              tmpscal2=-HEFFITD(I,J,K,bi,bj)
548              tmpscal3=-HSNOWITD(I,J,K,bi,bj)              tmpscal3=-HSNOWITD(I,J,K,bi,bj)
549              TICES(I,J,K,bi,bj)=celsius2K              TICES(I,J,K,bi,bj)=celsius2K
550              HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2  CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3  
 c            
             TICE(I,J,bi,bj)=celsius2K  
 c  
             HEFF(I,J,bi,bj) =HEFF(I,J,bi,bj) +tmpscal2  
             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2  
             HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3  
             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3  
551             ENDIF             ENDIF
552            ENDDO             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  #else
           tmpscal2=0. _d 0  
           tmpscal3=0. _d 0  
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  #endif
571           ENDDO           ENDDO
572          ENDDO          ENDDO
# Line 585  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
 CToM<<<  
           IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.  
 C     &        (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)) THEN  
            AREA(I,J,bi,bj)=0. _d 0  
582  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
583             DO K=1,nITD            DO K=1,nITD
584              AREAITD(I,J,K,bi,bj)=0. _d 0             IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.
585              HEFFITD(I,J,K,bi,bj)=0. _d 0       &         (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))
586              HSNOWITD(I,J,K,bi,bj)=0. _d 0       &      AREAITD(I,J,K,bi,bj)=0. _d 0
587             ENDDO            ENDDO
588    #else
589              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
591  #endif  #endif
           ENDIF  
 C>>>ToM  
592           ENDDO           ENDDO
593          ENDDO          ENDDO
594    
# Line 610  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
           IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN  
603  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
604             tmpscal2=AREA(I,J,bi,bj)            DO K=1,nITD
605  #endif             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
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)
 #ifdef SEAICE_ITD  
 c          ice area added (tmpscal3 is .ge.0):  
            tmpscal3=AREA(I,J,bi,bj)-tmpscal2  
 c          distribute this gain proportionally over categories  
            DO K=1,nITD  
             tmpscal4=AREAITD(I,J,K,bi,bj)/tmpscal2*tmpscal3  
             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4  
            ENDDO  
 #endif  
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 */
631          DO J=1,sNy          DO J=1,sNy
632           DO I=1,sNx           DO I=1,sNx
 #ifdef SEAICE_ITD  
           tmpscal2=AREA(I,J,bi,bj)  
 #endif  
633  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
634            DIAGarrayA(I,J) = AREA(I,J,bi,bj)            DIAGarrayA(I,J) = AREA(I,J,bi,bj)
635  #endif  #endif
# Line 646  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 637  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
637            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
638  #endif  #endif
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
641            ENDDO
642    #endif /* SEAICE_ITD */
643    
644  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
645  c         ice area subtracted (tmpscal3 is .ge.0):  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
646            tmpscal3=tmpscal2-AREA(I,J,bi,bj)  C    first, update AREA and HEFF:
647  c         distribute this loss proportionally over categories          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            DO K=1,nITD
656             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)             tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)
657       &                         -tmpscal3*areaFracFactor(I,J,K)             tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)
658            ENDDO            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  #endif
670           ENDDO           ENDDO
671          ENDDO          ENDDO
672    
673  #ifdef SEAICE_ITD  CToM finally make sure that all categories meet their thickness limits
674  C If AREAITD is changed due to regularization (but HEFFITD not) then the  C    which includes ridging as in item 2.5
675  C actual ice thickness (HEFFITD/AREAITD) in a category can be changed so  C    and update AREA, HEFF, and HSNOW
676  C that it does not fit its category limits anymore and redistribution is necessary          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
677          CALL SEAICE_ITD_REDIST(myTime, myIter, myThid)          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
678  C this should not affect the respective sums (AREA, HEFF, ...)  
 C ... except a non-conserving redistribution scheme is used; then call:  
 c        CALL SEAICE_ITD_SUM(myTime, myIter, myThid)  
679  #endif  #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 695  C 3) store regularized values of heff, h Line 706  C 3) store regularized values of heff, h
706             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
707             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
708             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             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            ENDDO
730           ENDDO           ENDDO
731          ENDDO          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  #endif
739    
740  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
# Line 792  CADJ STORE HSNWpreTH = comlev1_bibj, key Line 829  CADJ STORE HSNWpreTH = comlev1_bibj, key
829  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
830            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
831  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
832              tmpscal1          = MAX(SEAICE_area_reg,AREAITDpreTH(I,J,K))              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
833         &                              AREAITDpreTH(I,J,K))
834              hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1              hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
835              tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1              tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1
836              heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)              heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)
# Line 818  cif       Do not regularize when HEFFpre Line 856  cif       Do not regularize when HEFFpre
856              hsnowActualMult(I,J,K) = ZERO              hsnowActualMult(I,J,K) = ZERO
857              recip_heffActualMult(I,J,K)  = ZERO              recip_heffActualMult(I,J,K)  = ZERO
858            ENDIF            ENDIF
859  #else  #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 844  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  #endif /* SEAICE_ITD */
886    
887           ENDDO           ENDDO
888          ENDDO          ENDDO
# Line 962  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  #ifndef SEAICE_ITD
# Line 1037  CADJ &     comlev1_bibj, key = iicekey, Line 1078  CADJ &     comlev1_bibj, key = iicekey,
1078  C     update TICE & TICES  C     update TICE & TICES
1079  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1080  C     calculate area weighted mean  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)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1087       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)
1088  #else  #else
# Line 1047  C     calculate area weighted mean Line 1093  C     calculate area weighted mean
1093  C     average over categories  C     average over categories
1094  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1095  C     calculate area weighted mean  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)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1098       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1099             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
# Line 1093  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 1117  cgf just for those who may need to omit Line 1192  cgf just for those who may need to omit
1192            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1193           ENDDO           ENDDO
1194          ENDDO          ENDDO
1195  #ifdef SEAICE_ITD  #endif /* SEAICE_ITD */
         DO K=1,nITD  
          DO J=1,sNy  
           DO I=1,sNx  
            a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)  
      &          * convertQ2HI * AREAITDpreTH(I,J,K)  
            a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)  
      &          * convertQ2HI * AREAITDpreTH(I,J,K)  
            r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)  
 #ifdef SEAICE_DISABLE_SUBLIM  
            a_FWbySublimMult(I,J,K) = ZERO  
 #endif  
            a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce  
      &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)  
            r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)  
           ENDDO  
          ENDDO  
         ENDDO  
 #endif  
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 1152  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
          DO J=1,sNy  
           DO I=1,sNx  
            a_QbyATM_cover(I,J)   = 0. _d 0  
            r_QbyATM_cover(I,J)   = 0. _d 0  
            a_QSWbyATM_cover(I,J) = 0. _d 0  
           ENDDO  
          ENDDO  
1212  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1213           DO K=1,nITD           DO K=1,nITD
1214            DO J=1,sNy            DO J=1,sNy
# Line 1169  Cgf no additional dependency through ice Line 1219  Cgf no additional dependency through ice
1219             ENDDO             ENDDO
1220            ENDDO            ENDDO
1221           ENDDO           ENDDO
1222    #else
1223             DO J=1,sNy
1224              DO I=1,sNx
1225               a_QbyATM_cover(I,J)   = 0. _d 0
1226               r_QbyATM_cover(I,J)   = 0. _d 0
1227               a_QSWbyATM_cover(I,J) = 0. _d 0
1228              ENDDO
1229             ENDDO
1230  #endif  #endif
1231          ENDIF          ENDIF
1232  #endif  #endif
# Line 1214  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 1269  C     If anything is left, sublimate ice Line 1340  C     If anything is left, sublimate ice
1340            tmpscal2 =            tmpscal2 =
1341  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1342       &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)       &     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            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1345            HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2            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            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
# Line 1304  C       then update totals Line 1376  C       then update totals
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 1321  C          ice growth/melt due to ocean Line 1400  C          ice growth/melt due to ocean
1400  C          and hence weighted by fractional area of each thickness category  C          and hence weighted by fractional area of each thickness category
1401             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
1402       &                               -HEFFITD(I,J,K,bi,bj))       &                               -HEFFITD(I,J,K,bi,bj))
1403             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1             d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
1404             d_HEFFbyOCNonICE(I,J)=d_HEFFbyOCNonICE(I,J) + tmpscal1             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            ENDDO
1411           ENDDO           ENDDO
1412          ENDDO          ENDDO
1413  #endif  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
 #ifndef SEAICE_ITD  
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))
 #endif  
1424            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1425            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
1426  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
# Line 1339  C          and hence weighted by fractio Line 1428  C          and hence weighted by fractio
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 1362  Cgf no additional dependency through sno Line 1466  Cgf no additional dependency through sno
1466            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1467  #endif  #endif
1468            d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW            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            r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
1471  C         keep the total up to date, too  C         keep the total up to date, too
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          ENDDO          ENDDO
1476  #else  #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 1380  Cgf no additional dependency through sno Line 1485  Cgf no additional dependency through sno
1485            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1486  #endif  #endif
1487            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
1488              HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
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 */  #endif /* SEAICE_ITD */
1493          DO J=1,sNy  c ToM<<< debug seaice_growth
1494           DO I=1,sNx          WRITE(msgBuf,'(A,7F6.2)')
1495            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyATMonSNW(I,J)       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1496           ENDDO       &     HEFFITD(20,20,:,bi,bj)
1497          ENDDO          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 1412  Cgf warming conditions, the lab_sea resu Line 1520  Cgf warming conditions, the lab_sea resu
1520  #else  #else
1521             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
1522  c         Limit ice growth by potential melt by ocean  c         Limit ice growth by potential melt by ocean
1523       &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J)*areaFracFactor(I,J,K))       &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))
1524  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1525             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1526       &                                 + tmpscal2       &                                 + tmpscal2
# Line 1421  c         Limit ice growth by potential Line 1529  c         Limit ice growth by potential
1529             r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)             r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
1530       &                                 - tmpscal2       &                                 - tmpscal2
1531             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2             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            ENDDO
1538           ENDDO           ENDDO
1539          ENDDO          ENDDO
1540  #else  #else /* SEAICE_ITD */
1541          DO J=1,sNy          DO J=1,sNy
1542           DO I=1,sNx           DO I=1,sNx
1543    
# Line 1439  c         Limit ice growth by potential Line 1552  c         Limit ice growth by potential
1552            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1553            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1554            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1555           ENDDO            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
         ENDDO  
 #endif /* SEAICE_ITD */  
         DO J=1,sNy  
          DO I=1,sNx  
           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +d_HEFFbyATMonOCN_cover(I,J)  
1556    
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 1477  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
           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)  
1596           ENDDO           ENDDO
1597          ENDDO          ENDDO
1598  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
# Line 1489  C           add precip to the fresh wate Line 1604  C           add precip to the fresh wate
1604            ENDDO            ENDDO
1605           ENDDO           ENDDO
1606          ENDDO          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)
1611             ENDDO
1612            ENDDO
1613  #endif  #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 1511  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1639  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1639           DO J=1,sNy           DO J=1,sNy
1640            DO I=1,sNx            DO I=1,sNx
1641             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
1642       &                                        -HSNOW(I,J,bi,bj))       &                  -HSNOWITD(I,J,K,bi,bj))
1643             tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1644  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1645  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1646             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1647  #endif  #endif
            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2  
1648             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2             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            ENDDO
1652           ENDDO           ENDDO
1653          ENDDO          ENDDO
1654  #endif  #else /* SEAICE_ITD */
1655          DO J=1,sNy          DO J=1,sNy
1656           DO I=1,sNx           DO I=1,sNx
 #ifndef SEAICE_ITD  
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))
1658            tmpscal2=MIN(tmpscal1,0. _d 0)            tmpscal2=MIN(tmpscal1,0. _d 0)
1659  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
# Line 1533  Cgf no additional dependency through sno Line 1661  Cgf no additional dependency through sno
1661            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1662  #endif  #endif
1663            d_HSNWbyOCNonSNW(I,J) = tmpscal2            d_HSNWbyOCNonSNW(I,J) = tmpscal2
 #endif  
1664            r_QbyOCN(I,J)=r_QbyOCN(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)
1665       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
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 1569  C           or 0. otherwise (no melting Line 1704  C           or 0. otherwise (no melting
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  #ifdef SEAICE_ITD
1707    C         open water area fraction
1708              tmpscal0 = ONE-AREApreTH(I,J)
1709  C         determine thickness of new ice  C         determine thickness of new ice
1710  C         considering the entire open water area to refreeze  C         considering the entire open water area to refreeze
1711            tmpscal4 = tmpscal3/(ONE-AREApreTH(I,J))            tmpscal1 = tmpscal3/tmpscal0
1712  C         then add new ice volume to appropriate thickness category  C         then add new ice volume to appropriate thickness category
1713            DO K=1,nITD            DO K=1,nITD
1714             IF (tmpscal4.LT.Hlimit(K)) THEN             IF (tmpscal1.LT.Hlimit(K)) THEN
1715              HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3              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)              AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
1723       &                           + ONE-AREApreTH(I,J)       &                           + tmpscal0
1724                tmpscal0=ZERO
1725             ENDIF             ENDIF
1726            ENDDO            ENDDO
1727  C         in this case no open water is left after this step  #else
           AREA(I,J,bi,bj) = ONE  
 #endif  
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 1627  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1789  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1789       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1790             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1791             d_HEFFbyFLOODING(I,J)=tmpscal1             d_HEFFbyFLOODING(I,J)=tmpscal1
           ENDDO  
          ENDDO  
 #endif  
          DO J=1,sNy  
           DO I=1,sNx  
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 1675  C       under "gain of new ice over open Line 1840  C       under "gain of new ice over open
1840  C  C
1841  C       does not account for lateral melt of ice floes  C       does not account for lateral melt of ice floes
1842  C  C
1843    C        AREAITD is incremented in section "gain of new ice over open water" above
1844    C
1845          DO K=1,nITD          DO K=1,nITD
1846           DO J=1,sNy           DO J=1,sNy
1847            DO I=1,sNx            DO I=1,sNx
1848             IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN             IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
1849              AREAITD(I,J,K,bi,bj)=ZERO              AREAITD(I,J,K,bi,bj)=ZERO
1850             ENDIF             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            ENDDO
1856           ENDDO           ENDDO
1857          ENDDO          ENDDO
1858  C       update total AREA, HEFF, HSNOW  #else /* SEAICE_ITD */
         CALL SEAICE_ITD_SUM(myTime,myIter,myThid)  
 #else  
1859          DO J=1,sNy          DO J=1,sNy
1860           DO I=1,sNx           DO I=1,sNx
1861    
# Line 1770  Cgf 'bulk' linearization of area=f(HEFF) Line 1939  Cgf 'bulk' linearization of area=f(HEFF)
1939             ENDDO             ENDDO
1940            ENDDO            ENDDO
1941           ENDDO           ENDDO
 C        update total AREA, HEFF, HSNOW  
          CALL SEAICE_ITD_SUM(myTime,myIter,myThid)  
1942  #else  #else
1943           DO J=1,sNy           DO J=1,sNy
1944            DO I=1,sNx            DO I=1,sNx
# Line 1783  C            AREA(I,J,bi,bj) = 0.1 _d 0 Line 1950  C            AREA(I,J,bi,bj) = 0.1 _d 0
1950  #endif  #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.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22