/[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.2 by dimitri, Fri Apr 27 22:25:23 2012 UTC
# Line 137  c     The change of mean ice thickness d Line 137  c     The change of mean ice thickness d
137        _RL d_HEFFbyRLX         (1:sNx,1:sNy)        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
138  #endif  #endif
139    
140    #ifdef SEAICE_ITD
141    c     The change of mean ice area due to out-of-bounds values following
142    c     sea ice dynamics
143          _RL d_AREAbyNEG         (1:sNx,1:sNy)
144    #endif
145  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
146  c     sea ice dyhnamics  c     sea ice dynamics
147        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
148    
149  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 192  C     AREA_PRE :: hold sea-ice fraction Line 197  C     AREA_PRE :: hold sea-ice fraction
197        _RL AREApreTH           (1:sNx,1:sNy)        _RL AREApreTH           (1:sNx,1:sNy)
198        _RL HEFFpreTH           (1:sNx,1:sNy)        _RL HEFFpreTH           (1:sNx,1:sNy)
199        _RL HSNWpreTH           (1:sNx,1:sNy)        _RL HSNWpreTH           (1:sNx,1:sNy)
200    #ifdef SEAICE_ITD
201          _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
202          _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
203          _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
204          _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
205          _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)
206    #endif
207    
208  C     wind speed  C     wind speed
209        _RL UG                  (1:sNx,1:sNy)        _RL UG                  (1:sNx,1:sNy)
# Line 216  C temporary variables available for the Line 228  C temporary variables available for the
228  #endif  #endif
229    
230        INTEGER ilockey        INTEGER ilockey
231        INTEGER it  CToM<<<
232    C      INTEGER it
233          INTEGER IT, K
234    C>>>ToM
235        _RL pFac        _RL pFac
236        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
237        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
238        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
239    #ifdef SEAICE_ITD
240          _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
241          _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
242    #endif
243        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
244        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
245        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
246    #ifdef SEAICE_ITD
247          _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
248          _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
249    #endif
250  C     Helper variables: reciprocal of some constants  C     Helper variables: reciprocal of some constants
251        _RL recip_multDim        _RL recip_multDim
252        _RL recip_deltaTtherm        _RL recip_deltaTtherm
# Line 337  C ===================== Line 360  C =====================
360            d_HEFFbyRLX(I,J)           = 0.0 _d 0            d_HEFFbyRLX(I,J)           = 0.0 _d 0
361  #endif  #endif
362    
363    #ifdef SEAICE_ITD
364              d_AREAbyNEG(I,J)           = 0.0 _d 0
365    #endif
366            d_HEFFbyNEG(I,J)           = 0.0 _d 0            d_HEFFbyNEG(I,J)           = 0.0 _d 0
367            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
368            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
# Line 370  c Line 396  c
396              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
397              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
398              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
399    #ifdef SEAICE_ITD
400                r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
401                r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
402    #endif
403  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
404              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
405  #endif  #endif
# Line 406  C ====================================== Line 436  C ======================================
436  #endif  #endif
437           ENDDO           ENDDO
438          ENDDO          ENDDO
439    #ifdef SEAICE_ITD
440            DO K=1,nITD
441             DO J=1,sNy
442              DO I=1,sNx
443               HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
444               HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
445               AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
446               IF (AREA(I,J,bi,bj).GT.0) THEN
447                areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
448               ELSE
449                areaFracFactor(I,J,K)=ZERO
450               ENDIF
451               IF (HEFF(I,J,bi,bj).GT.0) THEN
452                heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)
453               ELSE
454                heffFracFactor(I,J,K)=ZERO
455               ENDIF
456              ENDDO
457             ENDDO
458            ENDDO
459    #endif
460    
461  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
462    
# Line 433  C         d_HEFFbyRLX(i,j) = 0. _d 0 Line 484  C         d_HEFFbyRLX(i,j) = 0. _d 0
484  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)
485              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
486             ENDIF             ENDIF
487    #ifdef SEAICE_ITD
488                AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
489         &                           +  d_AREAbyRLX(i,j)
490                HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
491         &                           +  d_HEFFbyRLX(i,j)
492    #endif
493              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)
494              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)
495           ENDDO           ENDDO
# Line 449  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 506  CADJ STORE area(:,:,bi,bj) = comlev1_bib
506  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
507          DO J=1,sNy          DO J=1,sNy
508           DO I=1,sNx           DO I=1,sNx
509    #ifdef SEAICE_ITD
510              DO K=1,nITD
511               tmpscal2=0. _d 0
512               tmpscal3=0. _d 0
513               tmpscal4=0. _d 0
514               tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)
515               HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
516               d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
517               tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)
518               HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
519               d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
520               tmpscal4=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
521               AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4
522               d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal4
523              ENDDO
524              AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyNEG(I,J)
525    #else
526            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
           HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)  
527            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
           HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)  
528            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)
529    #endif
530              HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
531              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
532           ENDDO           ENDDO
533          ENDDO          ENDDO
534    
# Line 464  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 539  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
539  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
540          DO J=1,sNy          DO J=1,sNy
541           DO I=1,sNx           DO I=1,sNx
542    #ifdef SEAICE_ITD
543              DO K=1,nITD
544               tmpscal2=0. _d 0
545               tmpscal3=0. _d 0
546               IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
547                tmpscal2=-HEFFITD(I,J,K,bi,bj)
548                tmpscal3=-HSNOWITD(I,J,K,bi,bj)
549                TICES(I,J,K,bi,bj)=celsius2K
550                HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2
551                HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
552    c          
553                TICE(I,J,bi,bj)=celsius2K
554    c
555                HEFF(I,J,bi,bj) =HEFF(I,J,bi,bj) +tmpscal2
556                d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
557                HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
558                d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
559               ENDIF
560              ENDDO
561    #else
562            tmpscal2=0. _d 0            tmpscal2=0. _d 0
563            tmpscal3=0. _d 0            tmpscal3=0. _d 0
564            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
# Line 478  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 573  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
573            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
574            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
575            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
576    #endif
577           ENDDO           ENDDO
578          ENDDO          ENDDO
579    
# Line 489  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 585  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
585  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
586          DO J=1,sNy          DO J=1,sNy
587           DO I=1,sNx           DO I=1,sNx
588    CToM<<<
589            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
590       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0  C     &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
591         &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) THEN
592               AREA(I,J,bi,bj)=0. _d 0
593    #ifdef SEAICE_ITD
594               DO K=1,nITD
595                AREAITD(I,J,K,bi,bj)=0. _d 0
596                HEFFITD(I,J,K,bi,bj)=0. _d 0
597                HSNOWITD(I,J,K,bi,bj)=0. _d 0
598               ENDDO
599    #endif
600              ENDIF
601    C>>>ToM
602           ENDDO           ENDDO
603          ENDDO          ENDDO
604    
# Line 503  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 611  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
611          DO J=1,sNy          DO J=1,sNy
612           DO I=1,sNx           DO I=1,sNx
613            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
614    #ifdef SEAICE_ITD
615               tmpscal2=AREA(I,J,bi,bj)
616    #endif
617             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)
618    #ifdef SEAICE_ITD
619    c          ice area added (tmpscal3 is .ge.0):
620               tmpscal3=AREA(I,J,bi,bj)-tmpscal2
621    c          distribute this gain proportionally over categories
622               DO K=1,nITD
623                tmpscal4=AREAITD(I,J,K,bi,bj)/tmpscal2*tmpscal3
624                AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4
625               ENDDO
626    #endif
627            ENDIF            ENDIF
628           ENDDO           ENDDO
629          ENDDO          ENDDO
# Line 516  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 636  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
636  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
637          DO J=1,sNy          DO J=1,sNy
638           DO I=1,sNx           DO I=1,sNx
639    #ifdef SEAICE_ITD
640              tmpscal2=AREA(I,J,bi,bj)
641    #endif
642  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
643            DIAGarrayA(I,J) = AREA(I,J,bi,bj)            DIAGarrayA(I,J) = AREA(I,J,bi,bj)
644  #endif  #endif
# Line 523  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 646  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
646            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
647  #endif  #endif
648            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)
649    #ifdef SEAICE_ITD
650    c         ice area subtracted (tmpscal3 is .ge.0):
651              tmpscal3=tmpscal2-AREA(I,J,bi,bj)
652    c         distribute this loss proportionally over categories
653              DO K=1,nITD
654               AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)
655         &                         -tmpscal3*areaFracFactor(I,J,K)
656              ENDDO
657    #endif
658           ENDDO           ENDDO
659          ENDDO          ENDDO
660    
661    #ifdef SEAICE_ITD
662    C If AREAITD is changed due to regularization (but HEFFITD not) then the
663    C actual ice thickness (HEFFITD/AREAITD) in a category can be changed so
664    C that it does not fit its category limits anymore and redistribution is necessary
665            CALL SEAICE_ITD_REDIST(myTime, myIter, myThid)
666    C this should not affect the respective sums (AREA, HEFF, ...)
667    C ... except a non-conserving redistribution scheme is used; then call:
668    c        CALL SEAICE_ITD_SUM(myTime, myIter, myThid)
669    #endif
670  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
671          ENDIF          ENDIF
672  #endif  #endif
# Line 547  C 3) store regularized values of heff, h Line 688  C 3) store regularized values of heff, h
688  #endif  #endif
689           ENDDO           ENDDO
690          ENDDO          ENDDO
691    #ifdef SEAICE_ITD
692            DO K=1,nITD
693             DO J=1,sNy
694              DO I=1,sNx
695               HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
696               HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
697               AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
698              ENDDO
699             ENDDO
700            ENDDO
701    #endif
702    
703  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
704  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
# Line 601  Cgf no additional dependency of air-sea Line 753  Cgf no additional dependency of air-sea
753             AREApreTH(I,J) = 0. _d 0             AREApreTH(I,J) = 0. _d 0
754            ENDDO            ENDDO
755           ENDDO           ENDDO
756    #ifdef SEAICE_ITD
757             DO K=1,nITD
758              DO J=1,sNy
759               DO I=1,sNx
760                HEFFITDpreTH(I,J,K) = 0. _d 0
761                HSNWITDpreTH(I,J,K) = 0. _d 0
762                AREAITDpreTH(I,J,K) = 0. _d 0
763               ENDDO
764              ENDDO
765             ENDDO
766    #endif
767          ENDIF          ENDIF
768  #endif  #endif
769    
# Line 620  CADJ STORE AREApreTH = comlev1_bibj, key Line 783  CADJ STORE AREApreTH = comlev1_bibj, key
783  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
784  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
785  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
786    #ifdef SEAICE_ITD
787            DO K=1,nITD
788    #endif
789          DO J=1,sNy          DO J=1,sNy
790           DO I=1,sNx           DO I=1,sNx
791    
792    #ifdef SEAICE_ITD
793              IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
794    #ifdef SEAICE_GROWTH_LEGACY
795                tmpscal1          = MAX(SEAICE_area_reg,AREAITDpreTH(I,J,K))
796                hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
797                tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1
798                heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)
799    #else /* SEAICE_GROWTH_LEGACY */
800    cif        regularize AREA with SEAICE_area_reg
801               tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)
802         &                     + area_reg_sq)
803    cif        heffActual calculated with the regularized AREA
804               tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1
805    cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
806               heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2
807         &                                  + hice_reg_sq)
808    cif        hsnowActual calculated with the regularized AREA
809               hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1
810    #endif /* SEAICE_GROWTH_LEGACY */
811    cif        regularize the inverse of heffActual by hice_reg
812               recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /
813         &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)
814         &                      + hice_reg_sq)
815    cif       Do not regularize when HEFFpreTH = 0
816              ELSE
817                heffActualMult(I,J,K) = ZERO
818                hsnowActualMult(I,J,K) = ZERO
819                recip_heffActualMult(I,J,K)  = ZERO
820              ENDIF
821    #else
822            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
823  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
824              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 844  cif       Do not regularize when HEFFpre
844              hsnowActual(I,J) = ZERO              hsnowActual(I,J) = ZERO
845              recip_heffActual(I,J)  = ZERO              recip_heffActual(I,J)  = ZERO
846            ENDIF            ENDIF
847    #endif
848    
849           ENDDO           ENDDO
850          ENDDO          ENDDO
851    #ifdef SEAICE_ITD
852            ENDDO
853    #endif
854    
855  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
856          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 861  cif       Do not regularize when HEFFpre
861  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
862  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
863  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
864    #ifdef SEAICE_ITD
865            DO K=1,nITD
866    #endif
867          DO J=1,sNy          DO J=1,sNy
868           DO I=1,sNx           DO I=1,sNx
869  c        The latent heat flux over the sea ice which  c        The latent heat flux over the sea ice which
870  c        will sublimate all of the snow and ice over one time  c        will sublimate all of the snow and ice over one time
871  c        step (W/m^2)  c        step (W/m^2)
872    #ifdef SEAICE_ITD
873              IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
874                latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *
875         &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +
876         &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)
877              ELSE
878                latentHeatFluxMaxMult(I,J,K) = ZERO
879              ENDIF
880    #else
881            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
882              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
883       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +
# Line 673  c        step (W/m^2) Line 885  c        step (W/m^2)
885            ELSE            ELSE
886              latentHeatFluxMax(I,J) = ZERO              latentHeatFluxMax(I,J) = ZERO
887            ENDIF            ENDIF
888    #endif
889           ENDDO           ENDDO
890          ENDDO          ENDDO
891    #ifdef SEAICE_ITD
892            ENDDO
893    #endif
894  #endif /* SEAICE_CAP_SUBLIM */  #endif /* SEAICE_CAP_SUBLIM */
895    
896  C ===================================================================  C ===================================================================
# Line 748  CADJ &                       key = iicek Line 964  CADJ &                       key = iicek
964  C--   Start loop over multi-categories  C--   Start loop over multi-categories
965          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
966  c        homogeneous distribution between 0 and 2 x heffActual  c        homogeneous distribution between 0 and 2 x heffActual
967    #ifndef SEAICE_ITD
968           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
969    #endif
970           DO J=1,sNy           DO J=1,sNy
971            DO I=1,sNx            DO I=1,sNx
972    #ifndef SEAICE_ITD
973    CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
974    C    (instead of heffActual and latentHeatFluxMax)
975             heffActualMult(I,J,IT)= heffActual(I,J)*pFac             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
976  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
977             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
978  #endif  #endif
979    #endif
980             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
981             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
982             TICE(I,J,bi,bj) = ZERO             TICE(I,J,bi,bj) = ZERO
# Line 780  CADJ &     comlev1_bibj, key = iicekey, Line 1002  CADJ &     comlev1_bibj, key = iicekey,
1002    
1003          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1004           CALL SEAICE_SOLVE4TEMP(           CALL SEAICE_SOLVE4TEMP(
1005    #ifdef SEAICE_ITD
1006         I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1007    #else
1008       I        UG, heffActualMult(1,1,IT), hsnowActual,       I        UG, heffActualMult(1,1,IT), hsnowActual,
1009    #endif
1010  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1011       I        latentHeatFluxMaxMult(1,1,IT),       I        latentHeatFluxMaxMult(1,1,IT),
1012  #endif  #endif
# Line 809  CADJ &     comlev1_bibj, key = iicekey, Line 1035  CADJ &     comlev1_bibj, key = iicekey,
1035           DO J=1,sNy           DO J=1,sNy
1036            DO I=1,sNx            DO I=1,sNx
1037  C     update TICE & TICES  C     update TICE & TICES
1038    #ifdef SEAICE_ITD
1039    C     calculate area weighted mean
1040               TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1041         &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)
1042    #else
1043             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1044       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
1045    #endif
1046             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1047  C     average over categories  C     average over categories
1048    #ifdef SEAICE_ITD
1049    C     calculate area weighted mean
1050               a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1051         &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1052               a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1053         &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1054               a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1055         &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)
1056    #else
1057             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1058       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1059             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1060       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1061             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1062       &          + a_FWbySublimMult(I,J,IT)*recip_multDim       &          + a_FWbySublimMult(I,J,IT)*recip_multDim
1063    #endif
1064            ENDDO            ENDDO
1065           ENDDO           ENDDO
1066          ENDDO          ENDDO
# Line 869  C     Negative sublimation is resublimat Line 1111  C     Negative sublimation is resublimat
1111  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1112  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
1113            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1114  #endif /* SEAICE_CAP_SUBLIM */  #endif
1115            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1116       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1117            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1118           ENDDO           ENDDO
1119          ENDDO          ENDDO
1120    #ifdef SEAICE_ITD
1121            DO K=1,nITD
1122             DO J=1,sNy
1123              DO I=1,sNx
1124               a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)
1125         &          * convertQ2HI * AREAITDpreTH(I,J,K)
1126               a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)
1127         &          * convertQ2HI * AREAITDpreTH(I,J,K)
1128               r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)
1129    #ifdef SEAICE_DISABLE_SUBLIM
1130               a_FWbySublimMult(I,J,K) = ZERO
1131    #endif
1132               a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce
1133         &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)
1134               r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)
1135              ENDDO
1136             ENDDO
1137            ENDDO
1138    #endif
1139    
1140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1141  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte
# Line 898  Cgf no additional dependency through ice Line 1159  Cgf no additional dependency through ice
1159             a_QSWbyATM_cover(I,J) = 0. _d 0             a_QSWbyATM_cover(I,J) = 0. _d 0
1160            ENDDO            ENDDO
1161           ENDDO           ENDDO
1162    #ifdef SEAICE_ITD
1163             DO K=1,nITD
1164              DO J=1,sNy
1165               DO I=1,sNx
1166                a_QbyATMmult_cover(I,J,K)   = 0. _d 0
1167                r_QbyATMmult_cover(I,J,K)   = 0. _d 0
1168                a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
1169               ENDDO
1170              ENDDO
1171             ENDDO
1172    #endif
1173          ENDIF          ENDIF
1174  #endif  #endif
1175    
# Line 961  C ====================================== Line 1233  C ======================================
1233  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1234  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1235  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1236    #ifdef SEAICE_ITD
1237            DO K=1,nITD
1238    #endif
1239          DO J=1,sNy          DO J=1,sNy
1240           DO I=1,sNx           DO I=1,sNx
1241  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1242            tmpscal2 =            tmpscal2 =
1243    #ifdef SEAICE_ITD
1244         &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)
1245         &             *SNOW2ICE),ZERO)
1246    C         accumulate change over ITD categories
1247              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1248         &                                                       *ICE2SNOW
1249              HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2
1250         &                                                       *ICE2SNOW
1251              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1252    C         keep total up to date, too
1253              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1254    #else
1255       &     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)
1256            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1257            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW
1258            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1259    #endif
1260           ENDDO           ENDDO
1261          ENDDO          ENDDO
1262  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 979  CADJ STORE r_FWbySublim    = comlev1_bib Line 1267  CADJ STORE r_FWbySublim    = comlev1_bib
1267           DO I=1,sNx           DO I=1,sNx
1268  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1269            tmpscal2 =            tmpscal2 =
1270    #ifdef SEAICE_ITD
1271         &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)
1272              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
1273              HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2
1274              r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1275    C         keep total up to date, too
1276              r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
1277    #else
1278       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1279            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
1280            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2
1281            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1282    #endif
1283           ENDDO           ENDDO
1284          ENDDO          ENDDO
1285          DO J=1,sNy          DO J=1,sNy
1286           DO I=1,sNx           DO I=1,sNx
1287  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.
1288  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
1289  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).
1290    #ifdef SEAICE_ITD
1291              a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
1292         &                              - r_FWbySublimMult(I,J,K)
1293              r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)
1294         &                              - r_FWbySublimMult(I,J,K)
1295             ENDDO
1296            ENDDO
1297    C       end K loop
1298            ENDDO
1299    C       then update totals      
1300            DO J=1,sNy
1301             DO I=1,sNx
1302    #endif
1303            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)
1304            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)
1305           ENDDO           ENDDO
# Line 1003  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1313  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1313  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1314  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1315    
1316    #ifdef SEAICE_ITD
1317            DO K=1,nITD
1318             DO J=1,sNy
1319              DO I=1,sNx
1320    C          ice growth/melt due to ocean heat is equally distributed under the ice
1321    C          and hence weighted by fractional area of each thickness category
1322               tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
1323         &                               -HEFFITD(I,J,K,bi,bj))
1324               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1
1325               d_HEFFbyOCNonICE(I,J)=d_HEFFbyOCNonICE(I,J) + tmpscal1
1326              ENDDO
1327             ENDDO
1328            ENDDO
1329    #endif
1330          DO J=1,sNy          DO J=1,sNy
1331           DO I=1,sNx           DO I=1,sNx
1332    #ifndef SEAICE_ITD
1333            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))
1334    #endif
1335            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1336            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)
1337  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
# Line 1022  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1348  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1348  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1349  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1350    
1351    #ifdef SEAICE_ITD
1352            DO K=1,nITD
1353             DO J=1,sNy
1354              DO I=1,sNx
1355    C     Convert to standard units (meters of ice) rather than to meters
1356    C     of snow. This appears to be more robust.
1357              tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)
1358         &                                           *SNOW2ICE)
1359              tmpscal2=MIN(tmpscal1,0. _d 0)
1360    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1361    Cgf no additional dependency through snow
1362              IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1363    #endif
1364              d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW
1365              r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
1366    C         keep the total up to date, too
1367              r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1368              ENDDO
1369             ENDDO
1370            ENDDO
1371    #else
1372          DO J=1,sNy          DO J=1,sNy
1373           DO I=1,sNx           DO I=1,sNx
1374  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 1033  Cgf no additional dependency through sno Line 1380  Cgf no additional dependency through sno
1380            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1381  #endif  #endif
1382            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW  
1383            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1384           ENDDO           ENDDO
1385          ENDDO          ENDDO
1386    #endif /* SEAICE_ITD */
1387            DO J=1,sNy
1388             DO I=1,sNx
1389              HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyATMonSNW(I,J)
1390             ENDDO
1391            ENDDO
1392    
1393  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1394  C ====================================================  C ====================================================
# Line 1051  Cgf where all experiments start in Janua Line 1403  Cgf where all experiments start in Janua
1403  Cgf the v1.81=>v1.82 revision would change results in  Cgf the v1.81=>v1.82 revision would change results in
1404  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1405    
1406    #ifdef SEAICE_ITD
1407            DO K=1,nITD
1408             DO J=1,sNy
1409              DO I=1,sNx
1410    #ifdef SEAICE_GROWTH_LEGACY
1411               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
1412    #else
1413               tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
1414    c         Limit ice growth by potential melt by ocean
1415         &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J)*areaFracFactor(I,J,K))
1416    #endif /* SEAICE_GROWTH_LEGACY */
1417               d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1418         &                                 + tmpscal2
1419               d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1420         &                                 + tmpscal2
1421               r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
1422         &                                 - tmpscal2
1423               HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2
1424              ENDDO
1425             ENDDO
1426            ENDDO
1427    #else
1428          DO J=1,sNy          DO J=1,sNy
1429           DO I=1,sNx           DO I=1,sNx
1430    
# Line 1065  c         Limit ice growth by potential Line 1439  c         Limit ice growth by potential
1439            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1440            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1441            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1442            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2           ENDDO
1443            ENDDO
1444    #endif /* SEAICE_ITD */
1445            DO J=1,sNy
1446             DO I=1,sNx
1447              HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +d_HEFFbyATMonOCN_cover(I,J)
1448    
1449  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1450            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
# Line 1101  C           add precip to the fresh wate Line 1480  C           add precip to the fresh wate
1480            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)
1481           ENDDO           ENDDO
1482          ENDDO          ENDDO
1483    #ifdef SEAICE_ITD
1484            DO K=1,nITD
1485             DO J=1,sNy
1486              DO I=1,sNx
1487               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)
1488         &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
1489              ENDDO
1490             ENDDO
1491            ENDDO
1492    #endif
1493  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1494  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1495  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
# Line 1116  Cph( very sensitive bit here by JZ Line 1505  Cph( very sensitive bit here by JZ
1505  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1506  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1507  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1508    
1509    #ifdef SEAICE_ITD
1510            DO K=1,nITD
1511             DO J=1,sNy
1512              DO I=1,sNx
1513               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
1514         &                                        -HSNOW(I,J,bi,bj))
1515               tmpscal2=MIN(tmpscal1,0. _d 0)
1516    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1517    Cgf no additional dependency through snow
1518               if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1519    #endif
1520               HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
1521               d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1522              ENDDO
1523             ENDDO
1524            ENDDO
1525    #endif
1526          DO J=1,sNy          DO J=1,sNy
1527           DO I=1,sNx           DO I=1,sNx
1528    #ifndef SEAICE_ITD
1529            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
1530            tmpscal2=MIN(tmpscal1,0. _d 0)            tmpscal2=MIN(tmpscal1,0. _d 0)
1531  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
# Line 1125  Cgf no additional dependency through sno Line 1533  Cgf no additional dependency through sno
1533            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1534  #endif  #endif
1535            d_HSNWbyOCNonSNW(I,J) = tmpscal2            d_HSNWbyOCNonSNW(I,J) = tmpscal2
1536    #endif
1537            r_QbyOCN(I,J)=r_QbyOCN(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)
1538       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1539            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)
# Line 1159  C           or 0. otherwise (no melting Line 1568  C           or 0. otherwise (no melting
1568            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1569            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1570            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1571    #ifdef SEAICE_ITD
1572    C         determine thickness of new ice
1573    C         considering the entire open water area to refreeze
1574              tmpscal4 = tmpscal3/(ONE-AREApreTH(I,J))
1575    C         then add new ice volume to appropriate thickness category
1576              DO K=1,nITD
1577               IF (tmpscal4.LT.Hlimit(K)) THEN
1578                HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
1579                AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
1580         &                           + ONE-AREApreTH(I,J)
1581               ENDIF
1582              ENDDO
1583    C         in this case no open water is left after this step
1584              AREA(I,J,bi,bj) = ONE
1585    #endif
1586            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1587           ENDDO           ENDDO
1588          ENDDO          ENDDO
# Line 1182  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1606  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1606  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1607  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1608          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1609    #ifdef SEAICE_ITD
1610             DO K=1,nITD
1611              DO J=1,sNy
1612               DO I=1,sNx
1613                tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow
1614         &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1615                tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
1616                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
1617                HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1
1618                HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
1619         &                            * ICE2SNOW
1620               ENDDO
1621              ENDDO
1622             ENDDO
1623    #else
1624           DO J=1,sNy           DO J=1,sNy
1625            DO I=1,sNx            DO I=1,sNx
1626             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1627       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1628             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1629             d_HEFFbyFLOODING(I,J)=tmpscal1             d_HEFFbyFLOODING(I,J)=tmpscal1
1630              ENDDO
1631             ENDDO
1632    #endif
1633             DO J=1,sNy
1634              DO I=1,sNx
1635             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)
1636             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1637       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
# Line 1220  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi Line 1664  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
1664  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1665  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1666    
1667    #ifdef SEAICE_ITD
1668    C       replaces Hibler '79 scheme and lead closing parameter
1669    C       because ITD accounts explicitly for lead openings and
1670    C       different melt rates due to varying ice thickness
1671    C
1672    C       only consider ice area loss due to total ice thickness loss
1673    C       ice area gain due to freezing of open water as handled above
1674    C       under "gain of new ice over open water"
1675    C
1676    C       does not account for lateral melt of ice floes
1677    C
1678            DO K=1,nITD
1679             DO J=1,sNy
1680              DO I=1,sNx
1681               IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
1682                AREAITD(I,J,K,bi,bj)=ZERO
1683               ENDIF
1684              ENDDO
1685             ENDDO
1686            ENDDO
1687    C       update total AREA, HEFF, HSNOW
1688            CALL SEAICE_ITD_SUM(myTime,myIter,myThid)
1689    #else
1690          DO J=1,sNy          DO J=1,sNy
1691           DO I=1,sNx           DO I=1,sNx
1692    
# Line 1289  C apply tendency Line 1756  C apply tendency
1756  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
1757           ENDDO           ENDDO
1758          ENDDO          ENDDO
1759    #endif /* SEAICE_ITD */
1760    
1761  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1762  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
1763          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
1764    #ifdef SEAICE_ITD
1765             DO K=1,nITD
1766              DO J=1,sNy
1767               DO I=1,sNx
1768                AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *
1769         &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )
1770               ENDDO
1771              ENDDO
1772             ENDDO
1773    C        update total AREA, HEFF, HSNOW
1774             CALL SEAICE_ITD_SUM(myTime,myIter,myThid)
1775    #else
1776           DO J=1,sNy           DO J=1,sNy
1777            DO I=1,sNx            DO I=1,sNx
1778  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 1780  C            AREA(I,J,bi,bj) = 0.1 _d 0
1780       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1781            ENDDO            ENDDO
1782           ENDDO           ENDDO
1783    #endif
1784          ENDIF          ENDIF
1785  #endif  #endif
1786    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22