/[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.8 by torge, Mon Oct 22 16:02:25 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 111  C     a_QSWbyATM_open   - short wave hea Line 115  C     a_QSWbyATM_open   - short wave hea
115  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2
116        _RL a_QSWbyATM_open     (1:sNx,1:sNy)        _RL a_QSWbyATM_open     (1:sNx,1:sNy)
117        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)
118  C     a_QbyOCN :: available heat (in in W/m^2) due to the  C     a_QbyOCN :: available heat (in W/m^2) due to the
119  C             interaction of the ice pack and the ocean surface  C             interaction of the ice pack and the ocean surface
120  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting
121  C             processes have been accounted for  C             processes have been accounted for
# Line 138  c     The change of mean ice thickness d Line 142  c     The change of mean ice thickness d
142  #endif  #endif
143    
144  c     The change of mean ice thickness due to out-of-bounds values following  c     The change of mean ice thickness due to out-of-bounds values following
145  c     sea ice dyhnamics  c     sea ice dynamics
146        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
147    
148  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
# Line 192  C     AREA_PRE :: hold sea-ice fraction Line 196  C     AREA_PRE :: hold sea-ice fraction
196        _RL AREApreTH           (1:sNx,1:sNy)        _RL AREApreTH           (1:sNx,1:sNy)
197        _RL HEFFpreTH           (1:sNx,1:sNy)        _RL HEFFpreTH           (1:sNx,1:sNy)
198        _RL HSNWpreTH           (1:sNx,1:sNy)        _RL HSNWpreTH           (1:sNx,1:sNy)
199    #ifdef SEAICE_ITD
200          _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
201          _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
202          _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
203          _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
204          _RL leadIceThickMin
205    #endif
206    
207  C     wind speed  C     wind speed
208        _RL UG                  (1:sNx,1:sNy)        _RL UG                  (1:sNx,1:sNy)
# Line 221  C temporary variables available for the Line 232  C temporary variables available for the
232        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
233        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
234        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
235    #ifdef SEAICE_ITD
236          _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
237          _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
238    #endif
239        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
240        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
241        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)        _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
242    #ifdef SEAICE_ITD
243          _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
244          _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
245    #endif
246  C     Helper variables: reciprocal of some constants  C     Helper variables: reciprocal of some constants
247        _RL recip_multDim        _RL recip_multDim
248        _RL recip_deltaTtherm        _RL recip_deltaTtherm
# Line 259  C ====================================== Line 278  C ======================================
278        ENDIF        ENDIF
279    
280  C     avoid unnecessary divisions in loops  C     avoid unnecessary divisions in loops
281    c#ifdef SEAICE_ITD
282    CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
283    C    (see SEAICE_SIZE.h and seaice_readparms.F)
284    c     SEAICE_multDim = nITD
285    c#endif
286        recip_multDim        = SEAICE_multDim        recip_multDim        = SEAICE_multDim
287        recip_multDim        = ONE / recip_multDim        recip_multDim        = ONE / recip_multDim
288  C     above/below: double/single precision calculation of recip_multDim  C     above/below: double/single precision calculation of recip_multDim
# Line 373  c Line 397  c
397  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
398              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
399  #endif  #endif
400    #ifdef SEAICE_ITD
401                r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
402                r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
403    #endif
404            ENDDO            ENDDO
405           ENDDO           ENDDO
406          ENDDO          ENDDO
# Line 406  C ====================================== Line 434  C ======================================
434  #endif  #endif
435           ENDDO           ENDDO
436          ENDDO          ENDDO
437    #ifdef SEAICE_ITD
438            DO IT=1,nITD
439             DO J=1,sNy
440              DO I=1,sNx
441               HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
442               HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
443               AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
444              ENDDO
445             ENDDO
446            ENDDO
447    #endif
448    
449  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
450    
# Line 433  C         d_HEFFbyRLX(i,j) = 0. _d 0 Line 472  C         d_HEFFbyRLX(i,j) = 0. _d 0
472  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)
473              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
474             ENDIF             ENDIF
475    #ifdef SEAICE_ITD
476                AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
477         &                           +  d_AREAbyRLX(i,j)
478                HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
479         &                           +  d_HEFFbyRLX(i,j)
480    #endif
481              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)
482              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)
483           ENDDO           ENDDO
# Line 449  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 494  CADJ STORE area(:,:,bi,bj) = comlev1_bib
494  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
495          DO J=1,sNy          DO J=1,sNy
496           DO I=1,sNx           DO I=1,sNx
497    #ifdef SEAICE_ITD
498              DO IT=1,nITD
499               tmpscal2=0. _d 0
500               tmpscal3=0. _d 0
501               tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
502               HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
503               d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
504               tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
505               HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
506               d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
507               AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
508              ENDDO
509    CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
510    C         by calling SEAICE_ITD_SUM
511    #else
512            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)
513            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)
514            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)
515            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)
516            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)
517    #endif
518           ENDDO           ENDDO
519          ENDDO          ENDDO
520    
# Line 464  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 525  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
525  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
526          DO J=1,sNy          DO J=1,sNy
527           DO I=1,sNx           DO I=1,sNx
528            tmpscal2=0. _d 0  #ifdef SEAICE_ITD
529            tmpscal3=0. _d 0            DO IT=1,nITD
530    #endif
531              tmpscal2=0. _d 0
532              tmpscal3=0. _d 0
533    #ifdef SEAICE_ITD
534               IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
535                tmpscal2=-HEFFITD(I,J,IT,bi,bj)
536                tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
537                TICES(I,J,IT,bi,bj)=celsius2K
538    CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
539               ENDIF
540               HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
541               HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
542    #else
543            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
544              tmpscal2=-HEFF(I,J,bi,bj)             tmpscal2=-HEFF(I,J,bi,bj)
545              tmpscal3=-HSNOW(I,J,bi,bj)             tmpscal3=-HSNOW(I,J,bi,bj)
546              TICE(I,J,bi,bj)=celsius2K             TICE(I,J,bi,bj)=celsius2K
547              DO IT=1,SEAICE_multDim             DO IT=1,SEAICE_multDim
548                TICES(I,J,IT,bi,bj)=celsius2K              TICES(I,J,IT,bi,bj)=celsius2K
549              ENDDO             ENDDO
550            ENDIF            ENDIF
551            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  
552            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
553    #endif
554              d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
555            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
556    #ifdef SEAICE_ITD
557              ENDDO
558    #endif
559           ENDDO           ENDDO
560          ENDDO          ENDDO
561    
# Line 489  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 567  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
567  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
568          DO J=1,sNy          DO J=1,sNy
569           DO I=1,sNx           DO I=1,sNx
570    #ifdef SEAICE_ITD
571              DO IT=1,nITD
572               IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
573         &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
574         &      AREAITD(I,J,IT,bi,bj)=0. _d 0
575              ENDDO
576    #else
577            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
578       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
579    #endif
580           ENDDO           ENDDO
581          ENDDO          ENDDO
582    
# Line 502  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 588  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
588  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
589          DO J=1,sNy          DO J=1,sNy
590           DO I=1,sNx           DO I=1,sNx
591    #ifdef SEAICE_ITD
592              DO IT=1,nITD
593               IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
594         &         (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
595    CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1
596    C           hence use SEAICE_area_floor devided by nITD
597    C           (or install a warning in e.g. seaice_readparms.F)
598                AREAITD(I,J,IT,bi,bj)=
599         &         MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
600               ENDIF
601              ENDDO
602    #else
603            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
604             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)
605            ENDIF            ENDIF
606    #endif
607           ENDDO           ENDDO
608          ENDDO          ENDDO
609  #endif /* DISABLE_AREA_FLOOR */  #endif /* DISABLE_AREA_FLOOR */
610    
611  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:
612    
613    CToM   for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
614    C      which is called at end of PART 1 below
615    #ifndef SEAICE_ITD
616  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
617  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
618  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 525  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 627  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
627            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)
628           ENDDO           ENDDO
629          ENDDO          ENDDO
630    #endif /* notSEAICE_ITD */
631    
632    #ifdef SEAICE_ITD
633    CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
634            DO J=1,sNy
635             DO I=1,sNx
636    C    TICES was changed above (item 1.25), now update TICE as ice volume
637    C     weighted average of TICES
638    C    also compute total of AREAITD (needed for finishing item 2.5, see below)
639              tmpscal1 = 0. _d 0
640              tmpscal2 = 0. _d 0
641              tmpscal3 = 0. _d 0
642              DO IT=1,nITD
643               tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
644               tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
645               tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
646              ENDDO
647              TICE(I,J,bi,bj)=tmpscal1/tmpscal2
648    C    lines of item 2.5 that were omitted:
649    C    in 2.5 these lines are executed before "ridging" is applied to AREA
650    C    hence we execute them here before SEAICE_ITD_REDIST is called
651    C    although this means that AREA has not been completely regularized
652    #ifdef ALLOW_DIAGNOSTICS
653              DIAGarrayA(I,J) = tmpscal3
654    #endif
655    #ifdef ALLOW_SITRACER
656              SItrAREA(I,J,bi,bj,1)=tmpscal3
657    #endif
658             ENDDO
659            ENDDO
660    
661    CToM finally make sure that all categories meet their thickness limits
662    C    which includes ridging as in item 2.5
663    C    and update AREA, HEFF, and HSNOW
664            CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
665            CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
666    
667    c ToM<<< debug seaice_growth
668            WRITE(msgBuf,'(A,7F8.4)')
669         &    ' SEAICE_GROWTH: Heff increments 0, HEFFITD = ',
670         &     HEFFITD(1,1,:,bi,bj)
671            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
672         &    SQUEEZE_RIGHT , myThid)
673            WRITE(msgBuf,'(A,7F8.4)')
674         &    ' SEAICE_GROWTH: Area increments 0, AREAITD = ',
675         &     AREAITD(1,1,:,bi,bj)
676            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
677         &    SQUEEZE_RIGHT , myThid)
678    #else
679            WRITE(msgBuf,'(A,7F8.4)')
680         &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
681         &     HEFF(1,1,bi,bj)
682            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
683         &    SQUEEZE_RIGHT , myThid)
684            WRITE(msgBuf,'(A,7F8.4)')
685         &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
686         &     AREA(1,1,bi,bj)
687            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
688         &    SQUEEZE_RIGHT , myThid)
689    c ToM>>>
690    #endif /* SEAICE_ITD */
691  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
692    C        end SEAICEadjMODE.EQ.0 statement:
693          ENDIF          ENDIF
694  #endif  #endif
695    
# Line 547  C 3) store regularized values of heff, h Line 710  C 3) store regularized values of heff, h
710  #endif  #endif
711           ENDDO           ENDDO
712          ENDDO          ENDDO
713    #ifdef SEAICE_ITD
714            DO IT=1,nITD
715             DO J=1,sNy
716              DO I=1,sNx
717               HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
718               HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
719               AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
720    
721    C memorize areal and volume fraction of each ITD category
722               IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
723                areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
724               ELSE
725    C           if there's no ice, potential growth starts in 1st category
726                IF (IT .EQ. 1) THEN
727                 areaFracFactor(I,J,IT)=ONE
728                ELSE
729                 areaFracFactor(I,J,IT)=ZERO
730                ENDIF
731               ENDIF
732              ENDDO
733             ENDDO
734            ENDDO
735    C prepare SItrHEFF to be computed as cumulative sum
736            DO iTr=2,5
737             DO J=1,sNy
738              DO I=1,sNx
739               SItrHEFF(I,J,bi,bj,iTr)=ZERO
740              ENDDO
741             ENDDO
742            ENDDO
743    C prepare SItrAREA to be computed as cumulative sum
744            DO J=1,sNy
745             DO I=1,sNx
746              SItrAREA(I,J,bi,bj,3)=ZERO
747             ENDDO
748            ENDDO
749    #endif
750    
751  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
752  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
# Line 601  Cgf no additional dependency of air-sea Line 801  Cgf no additional dependency of air-sea
801             AREApreTH(I,J) = 0. _d 0             AREApreTH(I,J) = 0. _d 0
802            ENDDO            ENDDO
803           ENDDO           ENDDO
804    #ifdef SEAICE_ITD
805             DO IT=1,nITD
806              DO J=1,sNy
807               DO I=1,sNx
808                HEFFITDpreTH(I,J,IT) = 0. _d 0
809                HSNWITDpreTH(I,J,IT) = 0. _d 0
810                AREAITDpreTH(I,J,IT) = 0. _d 0
811               ENDDO
812              ENDDO
813             ENDDO
814    #endif
815          ENDIF          ENDIF
816  #endif  #endif
817    
# Line 620  CADJ STORE AREApreTH = comlev1_bibj, key Line 831  CADJ STORE AREApreTH = comlev1_bibj, key
831  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
832  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
833  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
834    #ifdef SEAICE_ITD
835            DO IT=1,nITD
836    #endif
837          DO J=1,sNy          DO J=1,sNy
838           DO I=1,sNx           DO I=1,sNx
839    
840    #ifdef SEAICE_ITD
841              IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
842    #ifdef SEAICE_GROWTH_LEGACY
843                tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
844         &                              AREAITDpreTH(I,J,IT))
845                hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
846                tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
847                heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
848    #else /* SEAICE_GROWTH_LEGACY */
849    cif        regularize AREA with SEAICE_area_reg
850               tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
851         &                     + area_reg_sq)
852    cif        heffActual calculated with the regularized AREA
853               tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
854    cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
855               heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
856         &                                  + hice_reg_sq)
857    cif        hsnowActual calculated with the regularized AREA
858               hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
859    #endif /* SEAICE_GROWTH_LEGACY */
860    cif        regularize the inverse of heffActual by hice_reg
861               recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
862         &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
863         &                      + hice_reg_sq)
864    cif       Do not regularize when HEFFpreTH = 0
865              ELSE
866                heffActualMult(I,J,IT) = ZERO
867                hsnowActualMult(I,J,IT) = ZERO
868                recip_heffActualMult(I,J,IT)  = ZERO
869              ENDIF
870    #else /* SEAICE_ITD */
871            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
872  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
873              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 893  cif       Do not regularize when HEFFpre
893              hsnowActual(I,J) = ZERO              hsnowActual(I,J) = ZERO
894              recip_heffActual(I,J)  = ZERO              recip_heffActual(I,J)  = ZERO
895            ENDIF            ENDIF
896    #endif /* SEAICE_ITD */
897    
898           ENDDO           ENDDO
899          ENDDO          ENDDO
900    #ifdef SEAICE_ITD
901            ENDDO
902    #endif
903    
904  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
905          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 910  cif       Do not regularize when HEFFpre
910  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
911  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
912  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
913    #ifdef SEAICE_ITD
914            DO IT=1,nITD
915    #endif
916          DO J=1,sNy          DO J=1,sNy
917           DO I=1,sNx           DO I=1,sNx
918  c        The latent heat flux over the sea ice which  c        The latent heat flux over the sea ice which
919  c        will sublimate all of the snow and ice over one time  c        will sublimate all of the snow and ice over one time
920  c        step (W/m^2)  c        step (W/m^2)
921    #ifdef SEAICE_ITD
922              IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
923                latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
924         &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
925         &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
926         &         /AREAITDpreTH(I,J,IT)
927              ELSE
928                latentHeatFluxMaxMult(I,J,IT) = ZERO
929              ENDIF
930    #else
931            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
932              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *              latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
933       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +       &         (HEFFpreTH(I,J) * SEAICE_rhoIce +
# Line 673  c        step (W/m^2) Line 935  c        step (W/m^2)
935            ELSE            ELSE
936              latentHeatFluxMax(I,J) = ZERO              latentHeatFluxMax(I,J) = ZERO
937            ENDIF            ENDIF
938    #endif
939           ENDDO           ENDDO
940          ENDDO          ENDDO
941    #ifdef SEAICE_ITD
942            ENDDO
943    #endif
944  #endif /* SEAICE_CAP_SUBLIM */  #endif /* SEAICE_CAP_SUBLIM */
945    
946  C ===================================================================  C ===================================================================
# Line 705  cCADJ STORE TmixLoc = comlev1_bibj, key Line 971  cCADJ STORE TmixLoc = comlev1_bibj, key
971       I       TmixLoc,       I       TmixLoc,
972       O       a_QbyATM_open, a_QSWbyATM_open,       O       a_QbyATM_open, a_QSWbyATM_open,
973       I       bi, bj, myTime, myIter, myThid )       I       bi, bj, myTime, myIter, myThid )
974    c ToM<<< debugging
975            print*,' '
976            print*,'UG = ',UG(1,1)
977            print*,'Tsurf = ',TmixLoc(1,1)
978            print*,'a_QbyATM_open = ',a_QbyATM_open(1,1)
979            print*,' '
980    c ToM>>>
981    
982  C determine available heat due to the atmosphere -- for ice covered water  C determine available heat due to the atmosphere -- for ice covered water
983  C =======================================================================  C =======================================================================
# Line 746  CADJ &                       key = iicek Line 1019  CADJ &                       key = iicek
1019  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1020    
1021  C--   Start loop over multi-categories  C--   Start loop over multi-categories
1022    #ifdef SEAICE_ITD
1023    CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
1024    #endif
1025          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1026  c        homogeneous distribution between 0 and 2 x heffActual  c        homogeneous distribution between 0 and 2 x heffActual
1027    #ifndef SEAICE_ITD
1028           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
1029    #endif
1030           DO J=1,sNy           DO J=1,sNy
1031            DO I=1,sNx            DO I=1,sNx
1032    #ifndef SEAICE_ITD
1033    CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
1034    C    (instead of heffActual and latentHeatFluxMax)
1035             heffActualMult(I,J,IT)= heffActual(I,J)*pFac             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1036  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1037             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1038  #endif  #endif
1039    #endif
1040             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1041             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1042             TICE(I,J,bi,bj) = ZERO             TICE(I,J,bi,bj) = ZERO
# Line 780  CADJ &     comlev1_bibj, key = iicekey, Line 1062  CADJ &     comlev1_bibj, key = iicekey,
1062    
1063          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1064           CALL SEAICE_SOLVE4TEMP(           CALL SEAICE_SOLVE4TEMP(
1065    #ifdef SEAICE_ITD
1066         I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1067    #else
1068       I        UG, heffActualMult(1,1,IT), hsnowActual,       I        UG, heffActualMult(1,1,IT), hsnowActual,
1069    #endif
1070  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1071       I        latentHeatFluxMaxMult(1,1,IT),       I        latentHeatFluxMaxMult(1,1,IT),
1072  #endif  #endif
# Line 809  CADJ &     comlev1_bibj, key = iicekey, Line 1095  CADJ &     comlev1_bibj, key = iicekey,
1095           DO J=1,sNy           DO J=1,sNy
1096            DO I=1,sNx            DO I=1,sNx
1097  C     update TICE & TICES  C     update TICE & TICES
1098    #ifdef SEAICE_ITD
1099    C     calculate area weighted mean
1100    C     (although the ice's temperature relates to its energy content
1101    C      and hence should be averaged weighted by ice volume,
1102    C      the temperature here is a result of the fluxes through the ice surface
1103    C      computed individually for each single category in SEAICE_SOLVE4TEMP
1104    C      and hence is averaged area weighted [areaFracFactor])
1105               TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1106         &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1107    #else
1108             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1109       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
1110    #endif
1111             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1112  C     average over categories  C     average over categories
1113    #ifdef SEAICE_ITD
1114    C     calculate area weighted mean
1115    C     (fluxes are per unit (ice surface) area and are thus area weighted)
1116               a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1117         &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1118               a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1119         &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1120               a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1121         &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1122    #else
1123             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1124       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1125             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1126       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1127             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1128       &          + a_FWbySublimMult(I,J,IT)*recip_multDim       &          + a_FWbySublimMult(I,J,IT)*recip_multDim
1129    #endif
1130            ENDDO            ENDDO
1131           ENDDO           ENDDO
1132          ENDDO          ENDDO
1133    c ToM<<< debugging
1134            print*,' '
1135            print*,'after SOLVE4TEMP: '
1136            print*,'TICE  = ',TICE(1,1,bi,bj)
1137            print*,'TICES = ',TICES(1,1,:,bi,bj)
1138            print*,'a_QSWbyATM_cover     = ',a_QSWbyATM_cover(1,1)
1139            print*,'a_QSWbyATMmult_cover = ',a_QSWbyATMmult_cover(1,1,:)
1140            print*,' '
1141    c ToM>>>
1142    
1143  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1144  # ifdef ALLOW_DIAGNOSTICS  # ifdef ALLOW_DIAGNOSTICS
# Line 851  CADJ STORE a_FWbySublim    = comlev1_bib Line 1168  CADJ STORE a_FWbySublim    = comlev1_bib
1168  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1169    
1170  C switch heat fluxes from W/m2 to 'effective' ice meters  C switch heat fluxes from W/m2 to 'effective' ice meters
1171    #ifdef SEAICE_ITD
1172            DO IT=1,nITD
1173             DO J=1,sNy
1174              DO I=1,sNx
1175               a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
1176         &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1177               a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1178         &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1179    C and initialize r_QbyATMmult_cover
1180               r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
1181    C     Convert fresh water flux by sublimation to 'effective' ice meters.
1182    C     Negative sublimation is resublimation and will be added as snow.
1183    #ifdef SEAICE_DISABLE_SUBLIM
1184               a_FWbySublimMult(I,J,IT) = ZERO
1185    #endif
1186               a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
1187         &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
1188               r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
1189              ENDDO
1190             ENDDO
1191            ENDDO
1192            DO J=1,sNy
1193             DO I=1,sNx
1194              a_QbyATM_open(I,J)    = a_QbyATM_open(I,J)
1195         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1196              a_QSWbyATM_open(I,J)  = a_QSWbyATM_open(I,J)
1197         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1198    C and initialize r_QbyATM_open
1199              r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1200             ENDDO
1201            ENDDO
1202    #else /* SEAICE_ITD */
1203          DO J=1,sNy          DO J=1,sNy
1204           DO I=1,sNx           DO I=1,sNx
1205            a_QbyATM_cover(I,J)   = a_QbyATM_cover(I,J)            a_QbyATM_cover(I,J)   = a_QbyATM_cover(I,J)
# Line 869  C     Negative sublimation is resublimat Line 1218  C     Negative sublimation is resublimat
1218  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1219  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
1220            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1221  #endif /* SEAICE_CAP_SUBLIM */  #endif
1222            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1223       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1224            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1225           ENDDO           ENDDO
1226          ENDDO          ENDDO
1227    #endif /* SEAICE_ITD */
1228    
1229  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1230  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte
# Line 891  CADJ STORE r_FWbySublim    = comlev1_bib Line 1241  CADJ STORE r_FWbySublim    = comlev1_bib
1241  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1242  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1243          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1244    #ifdef SEAICE_ITD
1245             DO IT=1,nITD
1246              DO J=1,sNy
1247               DO I=1,sNx
1248                a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1249                r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1250                a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
1251               ENDDO
1252              ENDDO
1253             ENDDO
1254    #else
1255           DO J=1,sNy           DO J=1,sNy
1256            DO I=1,sNx            DO I=1,sNx
1257             a_QbyATM_cover(I,J)   = 0. _d 0             a_QbyATM_cover(I,J)   = 0. _d 0
# Line 898  Cgf no additional dependency through ice Line 1259  Cgf no additional dependency through ice
1259             a_QSWbyATM_cover(I,J) = 0. _d 0             a_QSWbyATM_cover(I,J) = 0. _d 0
1260            ENDDO            ENDDO
1261           ENDDO           ENDDO
1262    #endif
1263          ENDIF          ENDIF
1264  #endif  #endif
1265    
# Line 942  c available turbulent flux Line 1304  c available turbulent flux
1304            a_QbyOCN(i,j) =            a_QbyOCN(i,j) =
1305       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
1306            r_QbyOCN(i,j) = a_QbyOCN(i,j)            r_QbyOCN(i,j) = a_QbyOCN(i,j)
1307    c ToM<<< debugging
1308              if (i.eq.1 .and. j.eq.1) then
1309                print *, 'salt [psu]     = ',salt(i,j,kSurface,bi,bj)
1310                print *, 'theta [degC]   = ',theta(i,j,kSurface,bi,bj)
1311                print *, 'tempFrz [degC] = ',tempFrz
1312                print *, 'max turb flx [m]   = ',tmpscal2
1313                print *, 'avail trub flx [m] = ',a_QbyOCN(i,j)
1314              endif
1315    c ToM>>>
1316           ENDDO           ENDDO
1317          ENDDO          ENDDO
1318    
# Line 961  C ====================================== Line 1332  C ======================================
1332  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1333  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1334  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1335    #ifdef SEAICE_ITD
1336            DO IT=1,nITD
1337    #endif
1338          DO J=1,sNy          DO J=1,sNy
1339           DO I=1,sNx           DO I=1,sNx
1340  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1341            tmpscal2 =            tmpscal2 =
1342    #ifdef SEAICE_ITD
1343         &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1344         &             *SNOW2ICE),ZERO)
1345    C         accumulate change over ITD categories
1346              d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
1347         &                                                       *ICE2SNOW
1348              HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
1349         &                                                       *ICE2SNOW
1350              r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
1351    #else
1352       &     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)
1353            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1354            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW            HSNOW(I,J,bi,bj)    = HSNOW(I,J,bi,bj)  - tmpscal2*ICE2SNOW
1355            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1356    #endif
1357           ENDDO           ENDDO
1358          ENDDO          ENDDO
1359  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 979  CADJ STORE r_FWbySublim    = comlev1_bib Line 1364  CADJ STORE r_FWbySublim    = comlev1_bib
1364           DO I=1,sNx           DO I=1,sNx
1365  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1366            tmpscal2 =            tmpscal2 =
1367    #ifdef SEAICE_ITD
1368         &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1369    C         accumulate change over ITD categories
1370              d_HSNWbySublim(I,J)      = d_HSNWbySublim(I,J)      - tmpscal2
1371              HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
1372              r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
1373    #else
1374       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1375            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
1376            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2            HEFF(I,J,bi,bj)     = HEFF(I,J,bi,bj)   - tmpscal2
1377            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2            r_FWbySublim(I,J)   = r_FWbySublim(I,J) - tmpscal2
1378    #endif
1379           ENDDO           ENDDO
1380          ENDDO          ENDDO
1381          DO J=1,sNy          DO J=1,sNy
1382           DO I=1,sNx           DO I=1,sNx
1383  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.
1384  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
1385  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).
1386    #ifdef SEAICE_ITD
1387              a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1388         &                               - r_FWbySublimMult(I,J,IT)
1389              r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1390         &                               - r_FWbySublimMult(I,J,IT)
1391    #else
1392            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)
1393            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)
1394    #endif
1395           ENDDO           ENDDO
1396          ENDDO          ENDDO
1397    #ifdef SEAICE_ITD
1398    C       end IT loop
1399            ENDDO
1400    #endif
1401    c ToM<<< debug seaice_growth
1402            WRITE(msgBuf,'(A,7F8.4)')
1403    #ifdef SEAICE_ITD
1404         &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1405         &     HEFFITD(1,1,:,bi,bj)
1406            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1407         &    SQUEEZE_RIGHT , myThid)
1408            WRITE(msgBuf,'(A,7F8.4)')
1409         &    ' SEAICE_GROWTH: Area increments 1, AREAITD = ',
1410         &     AREAITD(1,1,:,bi,bj)
1411    #else
1412         &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1413         &     HEFF(1,1,bi,bj)
1414    #endif
1415            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1416         &    SQUEEZE_RIGHT , myThid)
1417    c ToM>>>
1418    
1419  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1420  C ===========================================================  C ===========================================================
# Line 1003  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1424  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1424  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1425  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1426    
1427    #ifdef SEAICE_ITD
1428            DO IT=1,nITD
1429             DO J=1,sNy
1430              DO I=1,sNx
1431    C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
1432    C          equally distributed under the ice and hence weighted by
1433    C          fractional area of each thickness category
1434               tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1435         &                               -HEFFITD(I,J,IT,bi,bj))
1436               d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
1437               HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal1
1438    #ifdef ALLOW_SITRACER
1439               SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1440         &                           + HEFFITD(I,J,IT,bi,bj)
1441    #endif
1442              ENDDO
1443             ENDDO
1444            ENDDO
1445            DO J=1,sNy
1446             DO I=1,sNx
1447              r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1448             ENDDO
1449            ENDDO
1450    #else /* SEAICE_ITD */
1451          DO J=1,sNy          DO J=1,sNy
1452           DO I=1,sNx           DO I=1,sNx
1453            d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))            d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
# Line 1013  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1458  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1458  #endif  #endif
1459           ENDDO           ENDDO
1460          ENDDO          ENDDO
1461    #endif /* SEAICE_ITD */
1462    c ToM<<< debug seaice_growth
1463            WRITE(msgBuf,'(A,7F8.4)')
1464    #ifdef SEAICE_ITD
1465         &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1466         &     HEFFITD(1,1,:,bi,bj)
1467            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1468         &    SQUEEZE_RIGHT , myThid)
1469            WRITE(msgBuf,'(A,7F8.4)')
1470         &    ' SEAICE_GROWTH: Area increments 2, AREAITD = ',
1471         &     AREAITD(1,1,:,bi,bj)
1472    #else
1473         &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1474         &     HEFF(1,1,bi,bj)
1475    #endif
1476            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1477         &    SQUEEZE_RIGHT , myThid)
1478    c ToM>>>
1479    
1480  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1481  C ==================================================================  C ==================================================================
# Line 1022  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1485  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1485  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1486  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1487    
1488    #ifdef SEAICE_ITD
1489            DO IT=1,nITD
1490             DO J=1,sNy
1491              DO I=1,sNx
1492    C     Convert to standard units (meters of ice) rather than to meters
1493    C     of snow. This appears to be more robust.
1494               tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1495         &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1496               tmpscal2=MIN(tmpscal1,0. _d 0)
1497    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1498    Cgf no additional dependency through snow
1499               IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1500    #endif
1501               d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1502         &                           + tmpscal2*ICE2SNOW
1503               HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
1504         &                           + tmpscal2*ICE2SNOW
1505               r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1506         &                           - tmpscal2
1507              ENDDO
1508             ENDDO
1509            ENDDO
1510    #else /* SEAICE_ITD */
1511          DO J=1,sNy          DO J=1,sNy
1512           DO I=1,sNx           DO I=1,sNx
1513  C     Convert to standard units (meters of ice) rather than to meters  C     Convert to standard units (meters of ice) rather than to meters
# Line 1037  Cgf no additional dependency through sno Line 1523  Cgf no additional dependency through sno
1523            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1524           ENDDO           ENDDO
1525          ENDDO          ENDDO
1526    #endif /* SEAICE_ITD */
1527    c ToM<<< debug seaice_growth
1528            WRITE(msgBuf,'(A,7F8.4)')
1529    #ifdef SEAICE_ITD
1530         &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1531         &     HEFFITD(1,1,:,bi,bj)
1532            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1533         &    SQUEEZE_RIGHT , myThid)
1534            WRITE(msgBuf,'(A,7F8.4)')
1535         &    ' SEAICE_GROWTH: Area increments 3, AREAITD = ',
1536         &     AREAITD(1,1,:,bi,bj)
1537    #else
1538         &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1539         &     HEFF(1,1,bi,bj)
1540    #endif
1541            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1542         &    SQUEEZE_RIGHT , myThid)
1543    c ToM>>>
1544    
1545  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1546  C ====================================================  C ====================================================
# Line 1051  Cgf where all experiments start in Janua Line 1555  Cgf where all experiments start in Janua
1555  Cgf the v1.81=>v1.82 revision would change results in  Cgf the v1.81=>v1.82 revision would change results in
1556  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1557    
1558    #ifdef SEAICE_ITD
1559            DO IT=1,nITD
1560             DO J=1,sNy
1561              DO I=1,sNx
1562    #ifdef SEAICE_GROWTH_LEGACY
1563               tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1564         &                     r_QbyATMmult_cover(I,J,IT))
1565    #else
1566               tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1567         &                     r_QbyATMmult_cover(I,J,IT)
1568    c         Limit ice growth by potential melt by ocean
1569         &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1570    #endif /* SEAICE_GROWTH_LEGACY */
1571               d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1572         &                                 + tmpscal2
1573               d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1574         &                                 + tmpscal2
1575               r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
1576         &                                 - tmpscal2
1577               HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
1578    
1579    #ifdef ALLOW_SITRACER
1580               SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1581         &                           + HEFFITD(I,J,IT,bi,bj)
1582    #endif
1583              ENDDO
1584             ENDDO
1585            ENDDO
1586    #else /* SEAICE_ITD */
1587          DO J=1,sNy          DO J=1,sNy
1588           DO I=1,sNx           DO I=1,sNx
1589    
# Line 1070  c         Limit ice growth by potential Line 1603  c         Limit ice growth by potential
1603  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1604            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1605  #endif  #endif
1606           ENDDO           ENDDO
1607          ENDDO          ENDDO
1608    #endif /* SEAICE_ITD */
1609    c ToM<<< debug seaice_growth
1610            WRITE(msgBuf,'(A,7F8.4)')
1611    #ifdef SEAICE_ITD
1612         &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1613         &     HEFFITD(1,1,:,bi,bj)
1614            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1615         &    SQUEEZE_RIGHT , myThid)
1616            WRITE(msgBuf,'(A,7F8.4)')
1617         &    ' SEAICE_GROWTH: Area increments 4, AREAITD = ',
1618         &     AREAITD(1,1,:,bi,bj)
1619    #else
1620         &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1621         &     HEFF(1,1,bi,bj)
1622    #endif
1623            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1624         &    SQUEEZE_RIGHT , myThid)
1625    c ToM>>>
1626    
1627  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1628  C depending on atmospheric conditions.  C depending on atmospheric conditions.
# Line 1098  C           add precip to the fresh wate Line 1649  C           add precip to the fresh wate
1649       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1650              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
1651            ENDIF            ENDIF
1652             ENDDO
1653            ENDDO
1654    #ifdef SEAICE_ITD
1655            DO IT=1,nITD
1656             DO J=1,sNy
1657              DO I=1,sNx
1658               HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1659         &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1660              ENDDO
1661             ENDDO
1662            ENDDO
1663    #else
1664            DO J=1,sNy
1665             DO I=1,sNx
1666            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)
1667           ENDDO           ENDDO
1668          ENDDO          ENDDO
1669    #endif
1670  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1671  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1672  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1673  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1674    c ToM<<< debug seaice_growth
1675            WRITE(msgBuf,'(A,7F8.4)')
1676    #ifdef SEAICE_ITD
1677         &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1678         &     HEFFITD(1,1,:,bi,bj)
1679            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1680         &    SQUEEZE_RIGHT , myThid)
1681            WRITE(msgBuf,'(A,7F8.4)')
1682         &    ' SEAICE_GROWTH: Area increments 5, AREAITD = ',
1683         &     AREAITD(1,1,:,bi,bj)
1684    #else
1685         &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1686         &     HEFF(1,1,bi,bj)
1687    #endif
1688            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1689         &    SQUEEZE_RIGHT , myThid)
1690    c ToM>>>
1691    
1692  C compute snow melt due to heat available from ocean.  C compute snow melt due to heat available from ocean.
1693  C =================================================================  C =================================================================
# Line 1116  Cph( very sensitive bit here by JZ Line 1699  Cph( very sensitive bit here by JZ
1699  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1700  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1701  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1702    
1703    #ifdef SEAICE_ITD
1704            DO IT=1,nITD
1705             DO J=1,sNy
1706              DO I=1,sNx
1707               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1708         &                  -HSNOWITD(I,J,IT,bi,bj))
1709               tmpscal2=MIN(tmpscal1,0. _d 0)
1710    #ifdef SEAICE_MODIFY_GROWTH_ADJ
1711    Cgf no additional dependency through snow
1712               if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1713    #endif
1714               d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1715               r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1716               HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) + tmpscal2
1717              ENDDO
1718             ENDDO
1719            ENDDO
1720    #else /* SEAICE_ITD */
1721          DO J=1,sNy          DO J=1,sNy
1722           DO I=1,sNx           DO I=1,sNx
1723            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
# Line 1130  Cgf no additional dependency through sno Line 1732  Cgf no additional dependency through sno
1732            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)
1733           ENDDO           ENDDO
1734          ENDDO          ENDDO
1735    #endif /* SEAICE_ITD */
1736  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1737  Cph)  Cph)
1738    c ToM<<< debug seaice_growth
1739            WRITE(msgBuf,'(A,7F8.4)')
1740    #ifdef SEAICE_ITD
1741         &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1742         &     HEFFITD(1,1,:,bi,bj)
1743            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1744         &    SQUEEZE_RIGHT , myThid)
1745            WRITE(msgBuf,'(A,7F8.4)')
1746         &    ' SEAICE_GROWTH: Area increments 6, AREAITD = ',
1747         &     AREAITD(1,1,:,bi,bj)
1748    #else
1749         &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1750         &     HEFF(1,1,bi,bj)
1751    #endif
1752            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1753         &    SQUEEZE_RIGHT , myThid)
1754    c ToM>>>
1755    
1756  C gain of new ice over open water  C gain of new ice over open water
1757  C ===============================  C ===============================
# Line 1156  C           impose -HEFF as the maxmum m Line 1776  C           impose -HEFF as the maxmum m
1776  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
1777              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
1778       &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)       &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
1779    c ToM<<< debugging
1780                if (I.eq.1 .and. J.eq.1) then
1781                 print*,'r_QbyATM_open(I,J) = ', r_QbyATM_open(I,J)
1782                 print*,'r_QbyOCN(i,j) = ', r_QbyOCN(i,j)
1783                 print*,'1 - AREApreTH = ', (1.0 _d 0 - AREApreTH(I,J))
1784                 print*,'tmpscal1 = ', tmpscal1
1785                 print*,' '
1786                 print*,'SWFracB = ', SWFracB
1787                 print*,'a_QSWbyATM_open(I,J) = ', a_QSWbyATM_open(I,J)
1788                 print*,'tmpscal2 = ', tmpscal2
1789                 print*,' '
1790                 print*,'facOpenGrow = ', facOpenGrow
1791                 print*,'HEFF(I,J,bi,bj) = ', HEFF(I,J,bi,bj)
1792                 print*,'facOpenMelt = ', facOpenMelt
1793                 print*,'MAX = ', MAX(tmpscal1-tmpscal2,
1794         &              -HEFF(I,J,bi,bj)*facOpenMelt)
1795                 print*,'HEFFM(I,J,bi,bj) = ', HEFFM(I,J,bi,bj)
1796                 print*,'tmpscal3 = ', tmpscal3
1797                 print*,' '
1798                endif
1799    c ToM>>>
1800            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1801            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1802            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1803    #ifdef SEAICE_ITD
1804    cC         open water area fraction
1805    c          tmpscal0 = ONE-AREApreTH(I,J)
1806    cC         determine thickness of new ice
1807    cctomC         considering the entire open water area to refreeze
1808    cctom          tmpscal1 = tmpscal3/tmpscal0
1809    cC         considering a minimum lead ice thickness of 10 cm
1810    cC         WATCH that leadIceThickMin is smaller that Hlimit(1)!
1811    c          leadIceThickMin = 0.1
1812    c          tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
1813    cC         adjust ice area fraction covered by new ice
1814    c         tmpscal0 = tmpscal3/tmpscal1
1815    cC         then add new ice volume to appropriate thickness category
1816    c          DO IT=1,nITD
1817    c          IF (tmpscal1.LT.Hlimit(IT)) THEN
1818    c            HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
1819    c           tmpscal3=ZERO
1820    cC not sure if AREAITD should be changed here since AREA is incremented
1821    cC   in PART 4 below in non-itd code
1822    cC in this scenario all open water is covered by new ice instantaneously,
1823    cC   i.e. no delayed lead closing is concidered such as is associated with
1824    cC   Hibler's h_0 parameter
1825    c           AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
1826    c     &                          + tmpscal0
1827    c           tmpscal0=ZERO
1828    c           ENDIF
1829    c          ENDDO
1830    ctom debugging: 1 category only
1831              HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj) + tmpscal3
1832    #else
1833            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1834    #endif
1835           ENDDO           ENDDO
1836          ENDDO          ENDDO
1837    
1838  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1839    #ifdef SEAICE_ITD
1840            DO IT=1,nITD
1841             DO J=1,sNy
1842              DO I=1,sNx
1843    c needs to be here to allow use also with LEGACY branch
1844               SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1845         &                           + HEFFITD(I,J,IT,bi,bj)
1846              ENDDO
1847             ENDDO
1848            ENDDO
1849    #else
1850          DO J=1,sNy          DO J=1,sNy
1851           DO I=1,sNx           DO I=1,sNx
1852  c needs to be here to allow use also with LEGACY branch  c needs to be here to allow use also with LEGACY branch
1853            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1854           ENDDO           ENDDO
1855          ENDDO          ENDDO
1856    #endif
1857  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1858    c ToM<<< debug seaice_growth
1859            WRITE(msgBuf,'(A,7F8.4)')
1860    #ifdef SEAICE_ITD
1861         &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1862         &     HEFFITD(1,1,:,bi,bj)
1863            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1864         &    SQUEEZE_RIGHT , myThid)
1865            WRITE(msgBuf,'(A,7F8.4)')
1866         &    ' SEAICE_GROWTH: Area increments 7, AREAITD = ',
1867         &     AREAITD(1,1,:,bi,bj)
1868    #else
1869         &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1870         &     HEFF(1,1,bi,bj)
1871            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1872         &    SQUEEZE_RIGHT , myThid)
1873            WRITE(msgBuf,'(A,7F8.4)')
1874         &    ' SEAICE_GROWTH: Area increments 7, AREA = ',
1875         &     AREA(1,1,bi,bj)
1876    #endif
1877            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1878         &    SQUEEZE_RIGHT , myThid)
1879    c ToM>>>
1880    
1881  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1882  C =================================  C =================================
# Line 1182  CADJ STORE heff(:,:,bi,bj) = comlev1_bib Line 1888  CADJ STORE heff(:,:,bi,bj) = comlev1_bib
1888  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1889  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1890          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1891    #ifdef SEAICE_ITD
1892             DO IT=1,nITD
1893              DO J=1,sNy
1894               DO I=1,sNx
1895                tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
1896         &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
1897         &                                        *recip_rhoConst
1898                tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
1899                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
1900                HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1
1901                HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
1902         &                            * ICE2SNOW
1903               ENDDO
1904              ENDDO
1905             ENDDO
1906    #else
1907           DO J=1,sNy           DO J=1,sNy
1908            DO I=1,sNx            DO I=1,sNx
1909             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow             tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
# Line 1191  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1913  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1913             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)
1914             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1915       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
1916            ENDDO            ENDDO
1917           ENDDO           ENDDO
1918    #endif
1919          ENDIF          ENDIF
1920  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1921    c ToM<<< debug seaice_growth
1922            WRITE(msgBuf,'(A,7F8.4)')
1923    #ifdef SEAICE_ITD
1924         &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1925         &     HEFFITD(1,1,:,bi,bj)
1926            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1927         &    SQUEEZE_RIGHT , myThid)
1928            WRITE(msgBuf,'(A,7F8.4)')
1929         &    ' SEAICE_GROWTH: Area increments 8, AREAITD = ',
1930         &     AREAITD(1,1,:,bi,bj)
1931    #else
1932         &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1933         &     HEFF(1,1,bi,bj)
1934            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1935         &    SQUEEZE_RIGHT , myThid)
1936            WRITE(msgBuf,'(A,7F8.4)')
1937         &    ' SEAICE_GROWTH: Area increments 8, AREA = ',
1938         &     AREA(1,1,bi,bj)
1939    #endif
1940            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1941         &    SQUEEZE_RIGHT , myThid)
1942    c ToM>>>
1943    
1944  C ===================================================================  C ===================================================================
1945  C ==========PART 4: determine ice cover fraction increments=========-  C ==========PART 4: determine ice cover fraction increments=========-
# Line 1220  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi Line 1965  CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
1965  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1966  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1967    
1968    c#ifdef SEAICE_ITD
1969    cC       replaces Hibler '79 scheme and lead closing parameter
1970    cC       because ITD accounts explicitly for lead openings and
1971    cC       different melt rates due to varying ice thickness
1972    cC
1973    cC       only consider ice area loss due to total ice thickness loss;
1974    cC       ice area gain due to freezing of open water is handled above
1975    cC       under "gain of new ice over open water"
1976    cC
1977    cC       does not account for lateral melt of ice floes
1978    cC
1979    cC        AREAITD is incremented in section "gain of new ice over open water" above
1980    cC
1981    c        DO IT=1,nITD
1982    c         DO J=1,sNy
1983    c          DO I=1,sNx
1984    c          IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1985    c           AREAITD(I,J,IT,bi,bj)=ZERO
1986    c          ENDIF
1987    c#ifdef ALLOW_SITRACER
1988    c           SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1989    c     &                           + AREAITD(I,J,IT,bi,bj)
1990    c#endif /* ALLOW_SITRACER */
1991    c         ENDDO
1992    c        ENDDO
1993    c       ENDDO
1994    c#else /* SEAICE_ITD */
1995          DO J=1,sNy          DO J=1,sNy
1996           DO I=1,sNx           DO I=1,sNx
1997    
1998    ctom<<< debugging
1999    #ifdef SEAICE_ITD
2000              HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
2001              AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
2002              HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
2003              HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
2004              AREApreTH(I,J)=AREAITDpreTH(I,J,1)
2005              recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
2006    #endif
2007    ctom>>> debugging
2008    
2009            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
2010             recip_HO=1. _d 0 / HO_south             recip_HO=1. _d 0 / HO_south
2011            ELSE            ELSE
# Line 1287  C apply tendency Line 2070  C apply tendency
2070            d_AREAbyOCN(I,J)=            d_AREAbyOCN(I,J)=
2071       &        HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )       &        HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
2072  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
2073    ctom<<< debugging
2074    #ifdef SEAICE_ITD
2075              HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
2076              AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
2077              HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
2078    #endif
2079    ctom>>> debugging
2080           ENDDO           ENDDO
2081          ENDDO          ENDDO
2082    c#endif /* SEAICE_ITD */
2083    
2084  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
2085  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
2086          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
2087    #ifdef SEAICE_ITD
2088             DO IT=1,nITD
2089              DO J=1,sNy
2090               DO I=1,sNx
2091                AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
2092         &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
2093               ENDDO
2094              ENDDO
2095             ENDDO
2096    #else
2097           DO J=1,sNy           DO J=1,sNy
2098            DO I=1,sNx            DO I=1,sNx
2099  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 2101  C            AREA(I,J,bi,bj) = 0.1 _d 0
2101       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )       &               ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
2102            ENDDO            ENDDO
2103           ENDDO           ENDDO
2104    #endif
2105          ENDIF          ENDIF
2106  #endif  #endif
2107    #ifdef SEAICE_ITD
2108    C check categories for consistency with limits after growth/melt
2109            CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
2110    C finally update total AREA, HEFF, HSNOW
2111            CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
2112    #endif
2113    
2114    c ToM<<< debugging
2115            DO J=1,sNy
2116             DO I=1,sNx
2117              if (I.eq.1 .and. J.eq.1) then
2118               print *, 'd_HEFFbyNEG(I,J)      = ', d_HEFFbyNEG(I,J)
2119               print *, 'd_HEFFbyOCNonICE(I,J) = ', d_HEFFbyOCNonICE(I,J)
2120               print *, 'd_HEFFbyATMonOCN(I,J) = ', d_HEFFbyATMonOCN(I,J)
2121               print *, 'd_HEFFbyATMonOCN_cover(I,J) = ',
2122         &               d_HEFFbyATMonOCN_cover(I,J)
2123               print *, 'd_HEFFbyATMonOCN_open(I,J) = ',
2124         &               d_HEFFbyATMonOCN_open(I,J)
2125               print *, 'd_HEFFbyFLOODING(I,J) = ', d_HEFFbyFLOODING(I,J)
2126               print *, 'd_HEFFbySublim(I,J)   = ', d_HEFFbySublim(I,J)
2127              endif
2128             ENDDO
2129            ENDDO
2130    c ToM>>>
2131  C ===================================================================  C ===================================================================
2132  C =============PART 5: determine ice salinity increments=============  C =============PART 5: determine ice salinity increments=============
2133  C ===================================================================  C ===================================================================
# Line 1516  C compute net heat flux leaving/entering Line 2341  C compute net heat flux leaving/entering
2341  C accounting for the part used in melt/freeze processes  C accounting for the part used in melt/freeze processes
2342  C =====================================================  C =====================================================
2343    
2344    #ifdef SEAICE_ITD
2345    C compute total of "mult" fluxes for ocean forcing
2346            DO J=1,sNy
2347             DO I=1,sNx
2348              a_QbyATM_cover(I,J)   = 0.0 _d 0
2349              r_QbyATM_cover(I,J)   = 0.0 _d 0
2350              a_QSWbyATM_cover(I,J) = 0.0 _d 0
2351              r_FWbySublim(I,J)     = 0.0 _d 0
2352             ENDDO
2353            ENDDO
2354            DO IT=1,nITD
2355             DO J=1,sNy
2356              DO I=1,sNx
2357    cToM if fluxes in W/m^2 then
2358    c           a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2359    c     &      + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2360    c           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2361    c     &      + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2362    c           a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2363    c     &      + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2364    c           r_FWbySublim(I,J)=r_FWbySublim(I,J)
2365    c     &      + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
2366    cToM if fluxes in effective ice meters, i.e. ice volume per area, then
2367               a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2368         &      + a_QbyATMmult_cover(I,J,IT)
2369               r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2370         &      + r_QbyATMmult_cover(I,J,IT)
2371               a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2372         &      + a_QSWbyATMmult_cover(I,J,IT)
2373               r_FWbySublim(I,J)=r_FWbySublim(I,J)
2374         &      + r_FWbySublimMult(I,J,IT)
2375              ENDDO
2376             ENDDO
2377            ENDDO
2378    #endif
2379    
2380  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
2381  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
2382  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
# Line 1540  C for backward compatibility it is left Line 2401  C for backward compatibility it is left
2401            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
2402           ENDDO           ENDDO
2403          ENDDO          ENDDO
2404    cToM<<< debugging
2405            print*,'------------------'
2406            print*,'OcnModFrc: QNET = ',QNET(1,1,bi,bj)
2407            print*,'OcnModFrc: QSW  = ',QSW(1,1,bi,bj)
2408            print*,' '
2409            print*,'r_QbyATM_cover  = ', r_QbyATM_cover(1,1)
2410            print*,'r_QbyATM_open   = ', r_QbyATM_open(1,1)
2411            print*,'a_QSWbyATM_cover = ', a_QSWbyATM_cover(1,1)
2412            print*,'d_HEFFbyOCNonICE = ', d_HEFFbyOCNonICE(1,1)
2413            print*,'d_HSNWbyOCNonSNW = ', d_HSNWbyOCNonSNW(1,1)
2414            print*,'d_HEFFbyNEG = ', d_HEFFbyNEG(1,1)
2415            print*,'d_HSNWbyNEG = ', d_HSNWbyNEG(1,1)
2416            print*,'SNOW2ICE = ',SNOW2ICE
2417            print*,'maskC       = ', maskC(1,1,kSurface,bi,bj)
2418            print*,'------------------'
2419    cToM>>>
2420    
2421  C switch heat fluxes from 'effective' ice meters to W/m2  C switch heat fluxes from 'effective' ice meters to W/m2
2422  C ======================================================  C ======================================================
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22