/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_growth.F
ViewVC logotype

Diff of /MITgcm_contrib/torge/itd/code/seaice_growth.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by dimitri, Fri Apr 27 22:25:23 2012 UTC revision 1.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 137  c     The change of mean ice thickness d Line 141  c     The change of mean ice thickness d
141        _RL d_HEFFbyRLX         (1:sNx,1:sNy)        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
142  #endif  #endif
143    
 #ifdef SEAICE_ITD  
 c     The change of mean ice area due to out-of-bounds values following  
 c     sea ice dynamics  
       _RL d_AREAbyNEG         (1:sNx,1:sNy)  
 #endif  
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 dynamics  c     sea ice dynamics
146        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
# Line 202  C     AREA_PRE :: hold sea-ice fraction Line 201  C     AREA_PRE :: hold sea-ice fraction
201        _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)        _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
202        _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)        _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
203        _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)        _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
204        _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)        _RL leadIceThickMin
205  #endif  #endif
206    
207  C     wind speed  C     wind speed
# Line 228  C temporary variables available for the Line 227  C temporary variables available for the
227  #endif  #endif
228    
229        INTEGER ilockey        INTEGER ilockey
230  CToM<<<        INTEGER it
 C      INTEGER it  
       INTEGER IT, K  
 C>>>ToM  
231        _RL pFac        _RL pFac
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)
# Line 282  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 360  C ===================== Line 361  C =====================
361            d_HEFFbyRLX(I,J)           = 0.0 _d 0            d_HEFFbyRLX(I,J)           = 0.0 _d 0
362  #endif  #endif
363    
 #ifdef SEAICE_ITD  
           d_AREAbyNEG(I,J)           = 0.0 _d 0  
 #endif  
364            d_HEFFbyNEG(I,J)           = 0.0 _d 0            d_HEFFbyNEG(I,J)           = 0.0 _d 0
365            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
366            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
# Line 396  c Line 394  c
394              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
395              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
396              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
397    #ifdef SEAICE_CAP_SUBLIM
398                latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
399    #endif
400  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
401              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
402              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
403  #endif  #endif
 #ifdef SEAICE_CAP_SUBLIM  
             latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0  
 #endif  
404            ENDDO            ENDDO
405           ENDDO           ENDDO
406          ENDDO          ENDDO
# Line 437  C ====================================== Line 435  C ======================================
435           ENDDO           ENDDO
436          ENDDO          ENDDO
437  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
438          DO K=1,nITD          DO IT=1,nITD
439           DO J=1,sNy           DO J=1,sNy
440            DO I=1,sNx            DO I=1,sNx
441             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
442             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
443             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
            IF (AREA(I,J,bi,bj).GT.0) THEN  
             areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)  
            ELSE  
             areaFracFactor(I,J,K)=ZERO  
            ENDIF  
            IF (HEFF(I,J,bi,bj).GT.0) THEN  
             heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)  
            ELSE  
             heffFracFactor(I,J,K)=ZERO  
            ENDIF  
444            ENDDO            ENDDO
445           ENDDO           ENDDO
446          ENDDO          ENDDO
# Line 507  CADJ STORE area(:,:,bi,bj) = comlev1_bib Line 495  CADJ STORE area(:,:,bi,bj) = comlev1_bib
495          DO J=1,sNy          DO J=1,sNy
496           DO I=1,sNx           DO I=1,sNx
497  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
498            DO K=1,nITD            DO IT=1,nITD
499             tmpscal2=0. _d 0             tmpscal2=0. _d 0
500             tmpscal3=0. _d 0             tmpscal3=0. _d 0
501             tmpscal4=0. _d 0             tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
502             tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)             HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
            HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2  
503             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
504             tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)             tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
505             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3             HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
506             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
507             tmpscal4=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)             AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4  
            d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal4  
508            ENDDO            ENDDO
509            AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyNEG(I,J)  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
510    C         by calling SEAICE_ITD_SUM
511  #else  #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)
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)
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  #endif
           HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)  
           HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)  
518           ENDDO           ENDDO
519          ENDDO          ENDDO
520    
# Line 540  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 526  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
526          DO J=1,sNy          DO J=1,sNy
527           DO I=1,sNx           DO I=1,sNx
528  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
529            DO K=1,nITD            DO IT=1,nITD
530             tmpscal2=0. _d 0  #endif
531             tmpscal3=0. _d 0            tmpscal2=0. _d 0
532             IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN            tmpscal3=0. _d 0
533              tmpscal2=-HEFFITD(I,J,K,bi,bj)  #ifdef SEAICE_ITD
534              tmpscal3=-HSNOWITD(I,J,K,bi,bj)             IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
535              TICES(I,J,K,bi,bj)=celsius2K              tmpscal2=-HEFFITD(I,J,IT,bi,bj)
536              HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2              tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
537              HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3              TICES(I,J,IT,bi,bj)=celsius2K
538  c            CToM        TICE will be updated at end of Part 1 together with AREA and HEFF
             TICE(I,J,bi,bj)=celsius2K  
 c  
             HEFF(I,J,bi,bj) =HEFF(I,J,bi,bj) +tmpscal2  
             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2  
             HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3  
             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3  
539             ENDIF             ENDIF
540            ENDDO             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  #else
           tmpscal2=0. _d 0  
           tmpscal3=0. _d 0  
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  #endif
559           ENDDO           ENDDO
560          ENDDO          ENDDO
# Line 585  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
 CToM<<<  
           IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.  
 C     &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0  
      &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) THEN  
            AREA(I,J,bi,bj)=0. _d 0  
570  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
571             DO K=1,nITD            DO IT=1,nITD
572              AREAITD(I,J,K,bi,bj)=0. _d 0             IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
573              HEFFITD(I,J,K,bi,bj)=0. _d 0       &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
574              HSNOWITD(I,J,K,bi,bj)=0. _d 0       &      AREAITD(I,J,IT,bi,bj)=0. _d 0
575             ENDDO            ENDDO
576    #else
577              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
579  #endif  #endif
           ENDIF  
 C>>>ToM  
580           ENDDO           ENDDO
581          ENDDO          ENDDO
582    
# Line 610  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
           IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN  
591  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
592             tmpscal2=AREA(I,J,bi,bj)            DO IT=1,nITD
593  #endif             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
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)
 #ifdef SEAICE_ITD  
 c          ice area added (tmpscal3 is .ge.0):  
            tmpscal3=AREA(I,J,bi,bj)-tmpscal2  
 c          distribute this gain proportionally over categories  
            DO K=1,nITD  
             tmpscal4=AREAITD(I,J,K,bi,bj)/tmpscal2*tmpscal3  
             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4  
            ENDDO  
 #endif  
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 */
619          DO J=1,sNy          DO J=1,sNy
620           DO I=1,sNx           DO I=1,sNx
 #ifdef SEAICE_ITD  
           tmpscal2=AREA(I,J,bi,bj)  
 #endif  
621  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
622            DIAGarrayA(I,J) = AREA(I,J,bi,bj)            DIAGarrayA(I,J) = AREA(I,J,bi,bj)
623  #endif  #endif
# Line 646  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 625  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
625            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
626  #endif  #endif
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
629            ENDDO
630    #endif /* notSEAICE_ITD */
631    
632  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
633  c         ice area subtracted (tmpscal3 is .ge.0):  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
634            tmpscal3=tmpscal2-AREA(I,J,bi,bj)          DO J=1,sNy
635  c         distribute this loss proportionally over categories           DO I=1,sNx
636            DO K=1,nITD  C    TICES was changed above (item 1.25), now update TICE as ice volume
637             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)  C     weighted average of TICES
638       &                         -tmpscal3*areaFracFactor(I,J,K)  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            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  #endif
658           ENDDO           ENDDO
659          ENDDO          ENDDO
660    
661  #ifdef SEAICE_ITD  CToM finally make sure that all categories meet their thickness limits
662  C If AREAITD is changed due to regularization (but HEFFITD not) then the  C    which includes ridging as in item 2.5
663  C actual ice thickness (HEFFITD/AREAITD) in a category can be changed so  C    and update AREA, HEFF, and HSNOW
664  C that it does not fit its category limits anymore and redistribution is necessary          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
665          CALL SEAICE_ITD_REDIST(myTime, myIter, myThid)          CALL SEAICE_ITD_SUM(bi, bj, 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:  c ToM<<< debug seaice_growth
668  c        CALL SEAICE_ITD_SUM(myTime, myIter, myThid)          WRITE(msgBuf,'(A,7F8.4)')
669  #endif       &    ' 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 689  C 3) store regularized values of heff, h Line 711  C 3) store regularized values of heff, h
711           ENDDO           ENDDO
712          ENDDO          ENDDO
713  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
714          DO K=1,nITD          DO IT=1,nITD
715           DO J=1,sNy           DO J=1,sNy
716            DO I=1,sNx            DO I=1,sNx
717             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
718             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
719             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             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            ENDDO
733           ENDDO           ENDDO
734          ENDDO          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  #endif
750    
751  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
# Line 754  Cgf no additional dependency of air-sea Line 802  Cgf no additional dependency of air-sea
802            ENDDO            ENDDO
803           ENDDO           ENDDO
804  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
805           DO K=1,nITD           DO IT=1,nITD
806            DO J=1,sNy            DO J=1,sNy
807             DO I=1,sNx             DO I=1,sNx
808              HEFFITDpreTH(I,J,K) = 0. _d 0              HEFFITDpreTH(I,J,IT) = 0. _d 0
809              HSNWITDpreTH(I,J,K) = 0. _d 0              HSNWITDpreTH(I,J,IT) = 0. _d 0
810              AREAITDpreTH(I,J,K) = 0. _d 0              AREAITDpreTH(I,J,IT) = 0. _d 0
811             ENDDO             ENDDO
812            ENDDO            ENDDO
813           ENDDO           ENDDO
# Line 784  CADJ STORE HEFFpreTH = comlev1_bibj, key Line 832  CADJ STORE HEFFpreTH = comlev1_bibj, key
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  #ifdef SEAICE_ITD
835          DO K=1,nITD          DO IT=1,nITD
836  #endif  #endif
837          DO J=1,sNy          DO J=1,sNy
838           DO I=1,sNx           DO I=1,sNx
839    
840  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
841            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
842  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
843              tmpscal1          = MAX(SEAICE_area_reg,AREAITDpreTH(I,J,K))              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
844              hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1       &                              AREAITDpreTH(I,J,IT))
845              tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1              hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
846              heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)              tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
847                heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
848  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
849  cif        regularize AREA with SEAICE_area_reg  cif        regularize AREA with SEAICE_area_reg
850             tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
851       &                     + area_reg_sq)       &                     + area_reg_sq)
852  cif        heffActual calculated with the regularized AREA  cif        heffActual calculated with the regularized AREA
853             tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1             tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
854  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
855             heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2             heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
856       &                                  + hice_reg_sq)       &                                  + hice_reg_sq)
857  cif        hsnowActual calculated with the regularized AREA  cif        hsnowActual calculated with the regularized AREA
858             hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
859  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
860  cif        regularize the inverse of heffActual by hice_reg  cif        regularize the inverse of heffActual by hice_reg
861             recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /             recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
862       &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
863       &                      + hice_reg_sq)       &                      + hice_reg_sq)
864  cif       Do not regularize when HEFFpreTH = 0  cif       Do not regularize when HEFFpreTH = 0
865            ELSE            ELSE
866              heffActualMult(I,J,K) = ZERO              heffActualMult(I,J,IT) = ZERO
867              hsnowActualMult(I,J,K) = ZERO              hsnowActualMult(I,J,IT) = ZERO
868              recip_heffActualMult(I,J,K)  = ZERO              recip_heffActualMult(I,J,IT)  = ZERO
869            ENDIF            ENDIF
870  #else  #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 844  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  #endif /* SEAICE_ITD */
897    
898           ENDDO           ENDDO
899          ENDDO          ENDDO
# Line 862  cif       Do not regularize when HEFFpre Line 911  cif       Do not regularize when HEFFpre
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  #ifdef SEAICE_ITD
914          DO K=1,nITD          DO IT=1,nITD
915  #endif  #endif
916          DO J=1,sNy          DO J=1,sNy
917           DO I=1,sNx           DO I=1,sNx
# Line 870  c        The latent heat flux over the s Line 919  c        The latent heat flux over the s
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  #ifdef SEAICE_ITD
922            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
923              latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *              latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
924       &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +       &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
925       &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)       &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
926         &         /AREAITDpreTH(I,J,IT)
927            ELSE            ELSE
928              latentHeatFluxMaxMult(I,J,K) = ZERO              latentHeatFluxMaxMult(I,J,IT) = ZERO
929            ENDIF            ENDIF
930  #else  #else
931            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
# Line 921  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 962  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  #ifndef SEAICE_ITD
# Line 1037  CADJ &     comlev1_bibj, key = iicekey, Line 1097  CADJ &     comlev1_bibj, key = iicekey,
1097  C     update TICE & TICES  C     update TICE & TICES
1098  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1099  C     calculate area weighted mean  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)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1106       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1107  #else  #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
# Line 1047  C     calculate area weighted mean Line 1112  C     calculate area weighted mean
1112  C     average over categories  C     average over categories
1113  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1114  C     calculate area weighted mean  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)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1117       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1118             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1119       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1120             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1121       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1122  #else  #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
# Line 1064  C     calculate area weighted mean Line 1130  C     calculate area weighted mean
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 1093  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 1117  cgf just for those who may need to omit Line 1224  cgf just for those who may need to omit
1224            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
1225           ENDDO           ENDDO
1226          ENDDO          ENDDO
1227  #ifdef SEAICE_ITD  #endif /* SEAICE_ITD */
         DO K=1,nITD  
          DO J=1,sNy  
           DO I=1,sNx  
            a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)  
      &          * convertQ2HI * AREAITDpreTH(I,J,K)  
            a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)  
      &          * convertQ2HI * AREAITDpreTH(I,J,K)  
            r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)  
 #ifdef SEAICE_DISABLE_SUBLIM  
            a_FWbySublimMult(I,J,K) = ZERO  
 #endif  
            a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce  
      &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)  
            r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)  
           ENDDO  
          ENDDO  
         ENDDO  
 #endif  
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 1152  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 1159  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
 #ifdef SEAICE_ITD  
          DO K=1,nITD  
           DO J=1,sNy  
            DO I=1,sNx  
             a_QbyATMmult_cover(I,J,K)   = 0. _d 0  
             r_QbyATMmult_cover(I,J,K)   = 0. _d 0  
             a_QSWbyATMmult_cover(I,J,K) = 0. _d 0  
            ENDDO  
           ENDDO  
          ENDDO  
1262  #endif  #endif
1263          ENDIF          ENDIF
1264  #endif  #endif
# Line 1214  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 1234  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1333  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
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  #ifdef SEAICE_ITD
1336          DO K=1,nITD          DO IT=1,nITD
1337  #endif  #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  #ifdef SEAICE_ITD
1343       &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1344       &             *SNOW2ICE),ZERO)       &             *SNOW2ICE),ZERO)
1345  C         accumulate change over ITD categories  C         accumulate change over ITD categories
1346            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
1347       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1348            HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2            HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
1349       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1350            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2            r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
 C         keep total up to date, too  
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2  
1351  #else  #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
# Line 1268  CADJ STORE r_FWbySublim    = comlev1_bib Line 1365  CADJ STORE r_FWbySublim    = comlev1_bib
1365  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1366            tmpscal2 =            tmpscal2 =
1367  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1368       &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1369            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2  C         accumulate change over ITD categories
1370            HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2            d_HSNWbySublim(I,J)      = d_HSNWbySublim(I,J)      - tmpscal2
1371            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2            HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
1372  C         keep total up to date, too            r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2  
1373  #else  #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
# Line 1288  C     If anything is left, it will be ev Line 1384  C     If anything is left, it will be ev
1384  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, 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  #ifdef SEAICE_ITD
1387            a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)            a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1388       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1389            r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)            r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1390       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1391           ENDDO  #else
         ENDDO  
 C       end K loop  
         ENDDO  
 C       then update totals        
         DO J=1,sNy  
          DO I=1,sNx  
 #endif  
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 1314  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1425  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1425  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1426    
1427  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1428          DO K=1,nITD          DO IT=1,nITD
1429           DO J=1,sNy           DO J=1,sNy
1430            DO I=1,sNx            DO I=1,sNx
1431  C          ice growth/melt due to ocean heat is equally distributed under the ice  C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
1432  C          and hence weighted by fractional area of each thickness category  C          equally distributed under the ice and hence weighted by
1433             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),  C          fractional area of each thickness category
1434       &                               -HEFFITD(I,J,K,bi,bj))             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1435             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1       &                               -HEFFITD(I,J,IT,bi,bj))
1436             d_HEFFbyOCNonICE(I,J)=d_HEFFbyOCNonICE(I,J) + tmpscal1             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            ENDDO
1443           ENDDO           ENDDO
1444          ENDDO          ENDDO
 #endif  
1445          DO J=1,sNy          DO J=1,sNy
1446           DO I=1,sNx           DO I=1,sNx
1447  #ifndef SEAICE_ITD            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
1452             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))
 #endif  
1454            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1455            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)
1456  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
# Line 1339  C          and hence weighted by fractio Line 1458  C          and hence weighted by fractio
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 1349  CADJ STORE r_QbyATM_cover = comlev1_bibj Line 1486  CADJ STORE r_QbyATM_cover = comlev1_bibj
1486  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1487    
1488  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1489          DO K=1,nITD          DO IT=1,nITD
1490           DO J=1,sNy           DO J=1,sNy
1491            DO I=1,sNx            DO I=1,sNx
1492  C     Convert to standard units (meters of ice) rather than to meters  C     Convert to standard units (meters of ice) rather than to meters
1493  C     of snow. This appears to be more robust.  C     of snow. This appears to be more robust.
1494            tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)             tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1495       &                                           *SNOW2ICE)       &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1496            tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1497  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1498  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1499            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1500  #endif  #endif
1501            d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW             d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1502            r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2       &                           + tmpscal2*ICE2SNOW
1503  C         keep the total up to date, too             HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
1504            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2       &                           + tmpscal2*ICE2SNOW
1505               r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1506         &                           - tmpscal2
1507            ENDDO            ENDDO
1508           ENDDO           ENDDO
1509          ENDDO          ENDDO
1510  #else  #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 1380  Cgf no additional dependency through sno Line 1519  Cgf no additional dependency through sno
1519            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1520  #endif  #endif
1521            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
1522              HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
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 */  #endif /* SEAICE_ITD */
1527          DO J=1,sNy  c ToM<<< debug seaice_growth
1528           DO I=1,sNx          WRITE(msgBuf,'(A,7F8.4)')
1529            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyATMonSNW(I,J)  #ifdef SEAICE_ITD
1530           ENDDO       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1531          ENDDO       &     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 1404  Cgf the v1.81=>v1.82 revision would chan Line 1556  Cgf the v1.81=>v1.82 revision would chan
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  #ifdef SEAICE_ITD
1559          DO K=1,nITD          DO IT=1,nITD
1560           DO J=1,sNy           DO J=1,sNy
1561            DO I=1,sNx            DO I=1,sNx
1562  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
1563             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))             tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1564         &                     r_QbyATMmult_cover(I,J,IT))
1565  #else  #else
1566             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)             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  c         Limit ice growth by potential melt by ocean
1569       &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J)*areaFracFactor(I,J,K))       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1570  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1571             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1572       &                                 + tmpscal2       &                                 + tmpscal2
1573             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1574       &                                 + tmpscal2       &                                 + tmpscal2
1575             r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)             r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
1576       &                                 - tmpscal2       &                                 - tmpscal2
1577             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2             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            ENDDO
1584           ENDDO           ENDDO
1585          ENDDO          ENDDO
1586  #else  #else /* SEAICE_ITD */
1587          DO J=1,sNy          DO J=1,sNy
1588           DO I=1,sNx           DO I=1,sNx
1589    
# Line 1439  c         Limit ice growth by potential Line 1598  c         Limit ice growth by potential
1598            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1599            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1600            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1601           ENDDO            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
         ENDDO  
 #endif /* SEAICE_ITD */  
         DO J=1,sNy  
          DO I=1,sNx  
           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +d_HEFFbyATMonOCN_cover(I,J)  
1602    
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 1477  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
           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)  
1652           ENDDO           ENDDO
1653          ENDDO          ENDDO
1654  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1655          DO K=1,nITD          DO IT=1,nITD
1656           DO J=1,sNy           DO J=1,sNy
1657            DO I=1,sNx            DO I=1,sNx
1658             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)             HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1659       &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)       &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1660            ENDDO            ENDDO
1661           ENDDO           ENDDO
1662          ENDDO          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)
1667             ENDDO
1668            ENDDO
1669  #endif  #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 1507  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1701  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1701  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1702    
1703  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1704          DO K=1,nITD          DO IT=1,nITD
1705           DO J=1,sNy           DO J=1,sNy
1706            DO I=1,sNx            DO I=1,sNx
1707             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1708       &                                        -HSNOW(I,J,bi,bj))       &                  -HSNOWITD(I,J,IT,bi,bj))
1709             tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1710  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1711  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1712             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1713  #endif  #endif
            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2  
1714             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2             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            ENDDO
1718           ENDDO           ENDDO
1719          ENDDO          ENDDO
1720  #endif  #else /* SEAICE_ITD */
1721          DO J=1,sNy          DO J=1,sNy
1722           DO I=1,sNx           DO I=1,sNx
 #ifndef SEAICE_ITD  
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))
1724            tmpscal2=MIN(tmpscal1,0. _d 0)            tmpscal2=MIN(tmpscal1,0. _d 0)
1725  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
# Line 1533  Cgf no additional dependency through sno Line 1727  Cgf no additional dependency through sno
1727            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1728  #endif  #endif
1729            d_HSNWbyOCNonSNW(I,J) = tmpscal2            d_HSNWbyOCNonSNW(I,J) = tmpscal2
 #endif  
1730            r_QbyOCN(I,J)=r_QbyOCN(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)
1731       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
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 1565  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  #ifdef SEAICE_ITD
1804  C         determine thickness of new ice  cC         open water area fraction
1805  C         considering the entire open water area to refreeze  c          tmpscal0 = ONE-AREApreTH(I,J)
1806            tmpscal4 = tmpscal3/(ONE-AREApreTH(I,J))  cC         determine thickness of new ice
1807  C         then add new ice volume to appropriate thickness category  cctomC         considering the entire open water area to refreeze
1808            DO K=1,nITD  cctom          tmpscal1 = tmpscal3/tmpscal0
1809             IF (tmpscal4.LT.Hlimit(K)) THEN  cC         considering a minimum lead ice thickness of 10 cm
1810              HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3  cC         WATCH that leadIceThickMin is smaller that Hlimit(1)!
1811              AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)  c          leadIceThickMin = 0.1
1812       &                           + ONE-AREApreTH(I,J)  c          tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
1813             ENDIF  cC         adjust ice area fraction covered by new ice
1814            ENDDO  c         tmpscal0 = tmpscal3/tmpscal1
1815  C         in this case no open water is left after this step  cC         then add new ice volume to appropriate thickness category
1816            AREA(I,J,bi,bj) = ONE  c          DO IT=1,nITD
1817  #endif  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 1607  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1889  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1889  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1890          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1891  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1892           DO K=1,nITD           DO IT=1,nITD
1893            DO J=1,sNy            DO J=1,sNy
1894             DO I=1,sNx             DO I=1,sNx
1895              tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow              tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
1896       &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
1897              tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))       &                                        *recip_rhoConst
1898              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1              tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
1899              HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
1900              HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1              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       &                            * ICE2SNOW
1903             ENDDO             ENDDO
1904            ENDDO            ENDDO
# Line 1627  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1910  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1910       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1911             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1912             d_HEFFbyFLOODING(I,J)=tmpscal1             d_HEFFbyFLOODING(I,J)=tmpscal1
           ENDDO  
          ENDDO  
 #endif  
          DO J=1,sNy  
           DO I=1,sNx  
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 1664  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  #ifdef SEAICE_ITD  c#ifdef SEAICE_ITD
1969  C       replaces Hibler '79 scheme and lead closing parameter  cC       replaces Hibler '79 scheme and lead closing parameter
1970  C       because ITD accounts explicitly for lead openings and  cC       because ITD accounts explicitly for lead openings and
1971  C       different melt rates due to varying ice thickness  cC       different melt rates due to varying ice thickness
1972  C  cC
1973  C       only consider ice area loss due to total ice thickness loss  cC       only consider ice area loss due to total ice thickness loss;
1974  C       ice area gain due to freezing of open water as handled above  cC       ice area gain due to freezing of open water is handled above
1975  C       under "gain of new ice over open water"  cC       under "gain of new ice over open water"
1976  C  cC
1977  C       does not account for lateral melt of ice floes  cC       does not account for lateral melt of ice floes
1978  C  cC
1979          DO K=1,nITD  cC        AREAITD is incremented in section "gain of new ice over open water" above
1980           DO J=1,sNy  cC
1981            DO I=1,sNx  c        DO IT=1,nITD
1982             IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN  c         DO J=1,sNy
1983              AREAITD(I,J,K,bi,bj)=ZERO  c          DO I=1,sNx
1984             ENDIF  c          IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1985            ENDDO  c           AREAITD(I,J,IT,bi,bj)=ZERO
1986           ENDDO  c          ENDIF
1987          ENDDO  c#ifdef ALLOW_SITRACER
1988  C       update total AREA, HEFF, HSNOW  c           SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1989          CALL SEAICE_ITD_SUM(myTime,myIter,myThid)  c     &                           + AREAITD(I,J,IT,bi,bj)
1990  #else  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 1754  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  #endif /* SEAICE_ITD */  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  #ifdef SEAICE_ITD
2088           DO K=1,nITD           DO IT=1,nITD
2089            DO J=1,sNy            DO J=1,sNy
2090             DO I=1,sNx             DO I=1,sNx
2091              AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *              AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
2092       &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )       &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
2093             ENDDO             ENDDO
2094            ENDDO            ENDDO
2095           ENDDO           ENDDO
 C        update total AREA, HEFF, HSNOW  
          CALL SEAICE_ITD_SUM(myTime,myIter,myThid)  
2096  #else  #else
2097           DO J=1,sNy           DO J=1,sNy
2098            DO I=1,sNx            DO I=1,sNx
# Line 1783  C            AREA(I,J,bi,bj) = 0.1 _d 0 Line 2104  C            AREA(I,J,bi,bj) = 0.1 _d 0
2104  #endif  #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 1997  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 2021  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.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22