/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_growth.F

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

revision 1.112 by gforget, Mon Feb 14 23:51:07 2011 UTC revision 1.113 by dimitri, Wed Feb 23 21:12:45 2011 UTC
# Line 117  C conversion factors to go from precip ( Line 117  C conversion factors to go from precip (
117        _RL convertPRECIP2HI, convertHI2PRECIP        _RL convertPRECIP2HI, convertHI2PRECIP
118    
119  C ICE/SNOW stocks tendencies associated with the various melt/freeze processes  C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
120        _RL d_AREA              (1:sNx,1:sNy)        _RL d_AREAbyATM         (1:sNx,1:sNy)
121    #ifdef FENTY_AREA_EXPANSION_CONTRACTION
122        _RL d_AREAbyOCN         (1:sNx,1:sNy)        _RL d_AREAbyOCN         (1:sNx,1:sNy)
123        _RL d_AREAbyATM_cover   (1:sNx,1:sNy)        _RL d_AREAbyICE         (1:sNx,1:sNy)
124        _RL d_AREAbyATM_open    (1:sNx,1:sNy)  #endif
125    
126    c     The change of mean ice thickness due to out-of-bounds values following
127    c     sea ice dyhnamics
128        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
       _RL d_HEFFbyOCN         (1:sNx,1:sNy)  
       _RL d_HEFFbyATM         (1:sNx,1:sNy)  
       _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)  
129    
130        _RL d_HEFFbyATM_open    (1:sNx,1:sNy)  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
131        _RL d_HEFFbyATM_cover   (1:sNx,1:sNy)        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
132    
133    c     The change of mean ice thickness due to atmospheric fluxes over the open water
134    c     fraction of the grid cell
135          _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
136    
137    c     The change of mean ice thickness due to atmospheric fluxes over the ice
138          _RL d_HEFFbyATMonICE    (1:sNx,1:sNy)
139    
140    c     The change of mean ice thickness due to flooding by snow
141          _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
142    
143        _RL d_HSNWbyNEG         (1:sNx,1:sNy)        _RL d_HSNWbyNEG         (1:sNx,1:sNy)
144        _RL d_HSNWbyATM         (1:sNx,1:sNy)        _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)
145        _RL d_HSNWbyOCN         (1:sNx,1:sNy)        _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)
146        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
147    
148        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
# Line 167  C     pathological cases thresholds Line 177  C     pathological cases thresholds
177        _RL heffTooThin, heffTooHeavy        _RL heffTooThin, heffTooHeavy
178    
179  C temporary variables available for the various computations  C temporary variables available for the various computations
180        _RL tmpscal0,tmpscal1, tmpscal2, tmpscal3, tmpscal4  #ifdef SEAICE_GROWTH_LEGACY
181          _RL tmpscal0
182    #endif
183          _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4
184        _RL tmparr1             (1:sNx,1:sNy)        _RL tmparr1             (1:sNx,1:sNy)
185    
186  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
# Line 192  C temporary variables available for the Line 205  C temporary variables available for the
205        PARAMETER ( nDim = 1 )        PARAMETER ( nDim = 1 )
206  #endif /* SEAICE_MULTICATEGORY */  #endif /* SEAICE_MULTICATEGORY */
207    
208    #ifdef FENTY_AREA_EXPANSION_CONTRACTION
209          _RL heff_star
210    #endif
211    
212  #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX  #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
213  C     Mixed Layer factor (dimensionless)  C     Factor by which we increase the upper ocean friction velocity (u*) when
214        _RL mlTurbFac  C     ice is absent in a grid cell  (dimensionless)
215  c     Stanton number (dimensionless)        _RL MixedLayerTurbulenceFactor
216    
217    c     The Stanton number for the McPhee
218    c     ocean-ice heat flux parameterization (dimensionless)
219        _RL STANTON_NUMBER        _RL STANTON_NUMBER
220  c     Friction velocity (m/s)  
221    c     A typical friction velocity beneath sea ice for the
222    c     McPhee heat flux parameterization (m/s)
223        _RL USTAR_BASE        _RL USTAR_BASE
224    
225          _RL surf_theta
226  #endif  #endif
227    
228  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
229    c     Helper variables for diagnostics
230        _RL DIAGarray     (1:sNx,1:sNy)        _RL DIAGarray     (1:sNx,1:sNy)
231          _RL DIAGarrayA    (1:sNx,1:sNy)
232          _RL DIAGarrayB    (1:sNx,1:sNy)
233          _RL DIAGarrayC    (1:sNx,1:sNy)
234          _RL DIAGarrayD    (1:sNx,1:sNy)
235  #endif  #endif
236    
237  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 225  C     Cutoff for iceload Line 254  C     Cutoff for iceload
254  C     FREEZING TEMP. OF SEA WATER (deg C)  C     FREEZING TEMP. OF SEA WATER (deg C)
255        TBC          = SEAICE_freeze        TBC          = SEAICE_freeze
256    
257    
258  C     RATIO OF SEA ICE DENSITY to SNOW DENSITY  C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
259        ICE2SNOW     = SEAICE_rhoIce/SEAICE_rhoSnow        ICE2SNOW     = SEAICE_rhoIce/SEAICE_rhoSnow
260    
# Line 235  C     HEAT OF FUSION OF SNOW (J/m^3) Line 265  C     HEAT OF FUSION OF SNOW (J/m^3)
265  C  C
266  C     note that QI/QS=ICE2SNOW -- except it wasnt in old code  C     note that QI/QS=ICE2SNOW -- except it wasnt in old code
267    
268    #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
269          STANTON_NUMBER = 0.0056 _d 0
270          USTAR_BASE = 0.0125 _d 0
271    #endif
272    
273  C conversion factors to go from Q (W/m2) to HEFF (ice meters)  C conversion factors to go from Q (W/m2) to HEFF (ice meters)
274        convertQ2HI=SEAICE_deltaTtherm/QI        convertQ2HI=SEAICE_deltaTtherm/QI
275        convertHI2Q=1/convertQ2HI        convertHI2Q=1/convertQ2HI
# Line 242  C conversion factors to go from precip ( Line 277  C conversion factors to go from precip (
277        convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce        convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
278        convertHI2PRECIP=1./convertPRECIP2HI        convertHI2PRECIP=1./convertPRECIP2HI
279    
 #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX  
       STANTON_NUMBER = 0.0056 _d 0  
       USTAR_BASE = 0.0125 _d 0  
 #endif  
   
280        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
281         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
282    
# Line 280  C ===================== Line 310  C =====================
310            a_QbyOCN (I,J)             = 0.0 _d 0            a_QbyOCN (I,J)             = 0.0 _d 0
311            r_QbyOCN (I,J)             = 0.0 _d 0            r_QbyOCN (I,J)             = 0.0 _d 0
312    
313            d_AREA(I,J)                = 0.0 _d 0            d_AREAbyATM(I,J)           = 0.0 _d 0
314    #ifdef FENTY_AREA_EXPANSION_CONTRACTION
315              d_AREAbyICE(I,J)           = 0.0 _d 0
316            d_AREAbyOCN(I,J)           = 0.0 _d 0            d_AREAbyOCN(I,J)           = 0.0 _d 0
317            d_AREAbyATM_cover(I,J)     = 0.0 _d 0  #endif
           d_AREAbyATM_open(I,J)      = 0.0 _d 0  
318    
319            d_HEFFbyOCN(I,J)           = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
320            d_HEFFbyATM(I,J)           = 0.0 _d 0  cif delta HEFF by atmospheric fluxes over the open ocean (ice-free)
321              d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
322    cif delta HEFF by atmospheric fluxes over the ice (not to the ocean)
323              d_HEFFbyATMonICE(I,J)      = 0.0 _d 0
324            d_HEFFbyFLOODING(I,J)      = 0.0 _d 0            d_HEFFbyFLOODING(I,J)      = 0.0 _d 0
325    
326            d_HEFFbyATM_open(I,J)      = 0.0 _d 0            d_HSNWbyATMonSNW(I,J)      = 0.0 _d 0
327            d_HEFFbyATM_cover(I,J)     = 0.0 _d 0            d_HSNWbyOCNonSNW(I,J)      = 0.0 _d 0
   
           d_HSNWbyATM(I,J)           = 0.0 _d 0  
           d_HSNWbyOCN(I,J)           = 0.0 _d 0  
328            d_HSNWbyRAIN(I,J)          = 0.0 _d 0            d_HSNWbyRAIN(I,J)          = 0.0 _d 0
329            a_FWbySublim(I,J)          = 0.0 _d 0            a_FWbySublim(I,J)          = 0.0 _d 0
330  #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET  #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
# Line 490  CADJ STORE HSNWpreTH = comlev1_bibj, key Line 521  CADJ STORE HSNWpreTH = comlev1_bibj, key
521            hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1            hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
522            tmpscal2         = HEFFpreTH(I,J)/tmpscal1            tmpscal2         = HEFFpreTH(I,J)/tmpscal1
523            heffActual(I,J)  = MAX(tmpscal2,hiceMin)            heffActual(I,J)  = MAX(tmpscal2,hiceMin)
524    Cgf do we need to keep this comment?
525    C     Capping the actual ice thickness effectively enforces a
526    C     minimum of heat flux through the ice and helps getting rid of
527    C     very thick ice.
528    Cdm actually, this does exactly the opposite, i.e., ice is thicker
529    Cdm when heffActual is capped, so I am commenting out
530    Cdm          heffActual(I,J)    = MIN(heffActual(I,J),9.0 _d +00)
531           ENDDO           ENDDO
532          ENDDO          ENDDO
533    
# Line 501  CADJ STORE HSNWpreTH = comlev1_bibj, key Line 539  CADJ STORE HSNWpreTH = comlev1_bibj, key
539  #endif  #endif
540    
541    
542    
543  C ===================================================================  C ===================================================================
544  C ================PART 2: determine heat fluxes/stocks===============  C ================PART 2: determine heat fluxes/stocks===============
545  C ===================================================================  C ===================================================================
# Line 690  C     Negative sublimation is resublimat Line 729  C     Negative sublimation is resublimat
729           ENDDO           ENDDO
730          ENDDO          ENDDO
731    
732    #ifdef ALLOW_DIAGNOSTICS
733            IF ( useDiagnostics ) THEN
734             IF ( DIAGNOSTICS_IS_ON('SIaQbATC',myThid) ) THEN
735              CALL DIAGNOSTICS_FILL(a_QbyATM_cover,
736         &      'SIaQbATC',0,1,3,bi,bj,myThid)
737             ENDIF
738             IF ( DIAGNOSTICS_IS_ON('SIaQbATO',myThid) ) THEN
739              CALL DIAGNOSTICS_FILL(a_QbyATM_open,
740         &      'SIaQbATO',0,1,3,bi,bj,myThid)
741             ENDIF
742           ENDIF
743    #endif
744  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
745  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
746  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
# Line 713  C ====================================== Line 764  C ======================================
764  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
765  #endif  #endif
766    
767    
768    
769          DO J=1,sNy          DO J=1,sNy
770           DO I=1,sNx           DO I=1,sNx
771    
772  #ifdef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
773             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
774  #endif /* SEAICE_VARIABLE_FREEZING_POINT */  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
775    
776  #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX  #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
777  #  ifdef MCPHEE_STEP_FUNCTION  
778             mlTurbFac = 12.5 _d 0  c     Bound the ocean temperature to be at or above the freezing point.
779             IF (AREApreTH(I,J) .GT. 0. _d 0) mlTurbFac = 1. _d 0             surf_theta = max(theta(I,J,kSurface,bi,bj), TBC)
780  #  else  #ifdef GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR
781             mlTurbFac = 12.5 _d 0 - 11.5 _d 0 *AREApreTH(I,J)             MixedLayerTurbulenceFactor = 12.5 _d 0 -
782  #  endif /* GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR */       &            11.5 _d 0 *AREApreTH(I,J)
783    #else
784    c     If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
785               IF (AREApreTH(I,J) .GT. 0. _d 0) THEN
786                   MixedLayerTurbulenceFactor = ONE
787                ELSE
788                   MixedLayerTurbulenceFactor = 12.5 _d 0
789               ENDIF
790    #endif /* GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR */
791    
792    c          The turbulent ocean-ice heat flux
793               a_QbyOCN(I,J) = -STANTON_NUMBER * USTAR_BASE * rhoConst *
794         &       HeatCapacity_Cp  *(surf_theta - TBC)*
795         &       MixedLayerTurbulenceFactor*hFacC(i,j,kSurface,bi,bj)
796    
797    c          The turbulent ocean-ice heat flux converted to meters
798    c          of potential ice melt
799               a_QbyOCN(I,J) = a_QbyOCN(I,J) * convertQ2HI
800    
801    #else
802    
803    c reproduces v1.109
804    #ifdef Reproduce_v1p109
805              IF ( .NOT. inAdMode ) THEN
806             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
807               tmpscal1 = STANTON_NUMBER * USTAR_BASE * mlTurbFac                a_QbyOCN(i,j) = -SEAICE_availHeatFrac
808         &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
809         &             * hFacC(i,j,kSurface,bi,bj) *
810         &             (HeatCapacity_Cp*rhoConst/QI)
811             ELSE             ELSE
812               tmpscal1 = 0. _d 0                a_QbyOCN(i,j) = -SEAICE_availHeatFracFrz
813         &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
814         &             * hFacC(i,j,kSurface,bi,bj) *
815         &             (HeatCapacity_Cp*rhoConst/QI)
816             ENDIF             ENDIF
817             tmpscal1 = tmpscal1            ELSE
818       &              * SEAICE_deltaTtherm * hFacC(i,j,kSurface,bi,bj)             a_QbyOCN(i,j) = 0.
819  #else /* MCPHEE_OCEAN_ICE_HEAT_FLUX */            ENDIF
820    #else  /* Reproduce_v1p109 */      
821    cif reproduces seaice_growth v1.111
822             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
823               tmpscal1 = SEAICE_availHeatFrac               tmpscal1 = SEAICE_availHeatFrac
824             ELSE             ELSE
# Line 741  CADJ STORE theta(:,:,:,bi,bj) = comlev1_ Line 826  CADJ STORE theta(:,:,:,bi,bj) = comlev1_
826             ENDIF             ENDIF
827             tmpscal1 = tmpscal1             tmpscal1 = tmpscal1
828       &              * dRf(kSurface) * hFacC(i,j,kSurface,bi,bj)       &              * dRf(kSurface) * hFacC(i,j,kSurface,bi,bj)
 #endif /* MCPHEE_OCEAN_ICE_HEAT_FLUX */  
829    
830             a_QbyOCN(i,j) = -tmpscal1 * (HeatCapacity_Cp*rhoConst/QI)             a_QbyOCN(i,j) = -tmpscal1 * (HeatCapacity_Cp*rhoConst/QI)
831       &                   * (theta(I,J,kSurface,bi,bj)-TBC)       &                   * (theta(I,J,kSurface,bi,bj)-TBC)
832             r_QbyOCN(i,j) = a_QbyOCN(i,j)  
833    #endif /* Reproduce_v1p109 */
834    
835    #endif /* MCPHEE_OCEAN_ICE_HEAT_FLUX */
836    
837    #ifdef ALLOW_DIAGNOSTICS
838              DIAGarrayA (I,J) = a_QbyOCN(I,J)
839    #endif
840              r_QbyOCN(i,j) = a_QbyOCN(i,j)
841    
842           ENDDO           ENDDO
843          ENDDO          ENDDO
844            
845    #ifdef ALLOW_DIAGNOSTICS
846            IF ( useDiagnostics ) THEN
847             IF ( DIAGNOSTICS_IS_ON('SIaQbOCN',myThid) ) THEN
848              CALL DIAGNOSTICS_FILL(DIAGarrayA,
849         &      'SIaQbOCN',0,1,3,bi,bj,myThid)
850             ENDIF
851            ENDIF
852    #endif
853    
854    
855    
856  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
857  #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ  #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ
# Line 770  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 874  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
874    
875          DO J=1,sNy          DO J=1,sNy
876           DO I=1,sNx           DO I=1,sNx
877            d_HEFFbyOCN(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))
878            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCN(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
879            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCN(I,J)            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
880    #ifdef ALLOW_DIAGNOSTICS
881              DIAGarrayA(I,J) = d_HEFFbyOCNonICE(i,j)
882    #endif
883           ENDDO           ENDDO
884          ENDDO          ENDDO
885    
886    #ifdef ALLOW_DIAGNOSTICS
887            IF ( useDiagnostics ) THEN
888             IF ( DIAGNOSTICS_IS_ON('SIdHbOCN',myThid) ) THEN
889              CALL DIAGNOSTICS_FILL(DIAGarrayA,
890         &      'SIdHbOCN',0,1,3,bi,bj,myThid)
891             ENDIF
892            ENDIF
893    #endif
894    
895  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
896  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
897          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
898           IF ( DIAGNOSTICS_IS_ON('SIyneg  ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIyneg  ',myThid) ) THEN
899            CALL DIAGNOSTICS_FILL(d_HEFFbyOCN,            CALL DIAGNOSTICS_FILL(d_HEFFbyOCNonICE,
900       &      'SIyneg  ',0,1,1,bi,bj,myThid)       &      'SIyneg  ',0,1,1,bi,bj,myThid)
901           ENDIF           ENDIF
902          ENDIF          ENDIF
# Line 829  CADJ STORE r_QbyATM_cover = comlev1_bibj Line 945  CADJ STORE r_QbyATM_cover = comlev1_bibj
945  Cgf no additional dependency through snow  Cgf no additional dependency through snow
946            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
947  #endif  #endif
948            d_HSNWbyATM(I,J)= tmpscal2            d_HSNWbyATMonSNW(I,J)= tmpscal2
949            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2
950            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2/ICE2SNOW            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2/ICE2SNOW
951           ENDDO           ENDDO
# Line 861  C     If anything is left, sublimate ice Line 977  C     If anything is left, sublimate ice
977  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
978  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
979  CADJ STORE r_QbyATM_cover  = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_cover  = comlev1_bibj,key=iicekey,byte=isbyte
 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  
980  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
981    
982  Cgf note: this block is not actually tested by lab_sea  Cgf note: this block is not actually tested by lab_sea
# Line 871  Cgf warming conditions, the lab_sea resu Line 986  Cgf warming conditions, the lab_sea resu
986    
987          DO J=1,sNy          DO J=1,sNy
988           DO I=1,sNx           DO I=1,sNx
989    cif This is a real mess because the default behavior of the code was changed between v1.190
990    cif v.1.111.  After #ifdef FENTY_DELTA_HEFF_OPEN_WATER_FLUXES I left the v1.111 code.
991    cif It seems that now undef SEAICE_GROWTH_LEGACY = define FENTY_DELTA_HEFF_OPEN_WATER FLUXES
992    #ifdef FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
993    c         Thickening of the existing ice pack is limited by the
994    c         residual capability of the ocean to melt (weighted by
995    c         the ice-covered area)
996              tmpscal2 = MAX(-HEFF(i,j,bi,bj) ,
997         &       r_QbyATM_cover(I,J) + AREApreTH(I,J) * r_QbyOCN(I,J))
998        
999    #else /* FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */
1000    
1001  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
1002            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1003  #else  #else
1004            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1005  c         Limit ice growth by potential melt by ocean  c         Limit ice growth by potential melt by ocean
1006       &        AREApreTH(I,J) * r_QbyOCN(I,J))       &        AREApreTH(I,J) * r_QbyOCN(I,J))
1007  #endif  #endif /* SEAICE_GROWTH_LEGACY */
1008            d_HEFFbyATM_cover(I,J)=tmpscal2  
1009            d_HEFFbyATM(I,J)=d_HEFFbyATM(I,J)+tmpscal2  #endif /*FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */
1010    
1011              d_HEFFbyATMonICE(I,J)=d_HEFFbyATMonICE(I,J)+tmpscal2
1012            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1013            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
1014    
1015    #ifdef ALLOW_DIAGNOSTICS
1016              DIAGarrayA(I,J) = d_HEFFbyATMonICE(I,J)
1017    #endif
1018           ENDDO           ENDDO
1019          ENDDO          ENDDO
1020    
1021    #ifdef ALLOW_DIAGNOSTICS
1022             IF ( useDiagnostics ) THEN
1023              IF ( DIAGNOSTICS_IS_ON('SIdHbATC',myThid) ) THEN
1024              CALL DIAGNOSTICS_FILL(DIAGarrayA,
1025         &      'SIdHbATC',0,1,3,bi,bj,myThid)
1026             ENDIF
1027            ENDIF
1028    #endif /* ALLOW DIAGNOSTICS */
1029    
1030    
1031  C attribute precip to fresh water or snow stock,  C attribute precip to fresh water or snow stock,
1032  C depending on atmospheric conditions.  C depending on atmospheric conditions.
1033  C =================================================  C =================================================
# Line 949  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1092  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1092  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1093            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1094  #endif  #endif
1095            d_HSNWbyOCN(I,J) = tmpscal2            d_HSNWbyOCNonSNW(I,J) = tmpscal2
1096            r_QbyOCN(I,J)=r_QbyOCN(I,J)            r_QbyOCN(I,J)=r_QbyOCN(I,J)
1097       &                               -d_HSNWbyOCN(I,J)/ICE2SNOW       &                               -d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
1098            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCN(I,J)            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1099           ENDDO           ENDDO
1100          ENDDO          ENDDO
1101  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
# Line 961  Cph) Line 1104  Cph)
1104  C gain of new ice over open water  C gain of new ice over open water
1105  C ===============================  C ===============================
1106  #ifndef SEAICE_GROWTH_LEGACY  #ifndef SEAICE_GROWTH_LEGACY
 #ifdef SEAICE_DO_OPEN_WATER_GROWTH  
1107  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1108  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1109  CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1110  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1111    CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1112  CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1113  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1114          DO J=1,sNy          DO J=1,sNy
1115           DO I=1,sNx           DO I=1,sNx
1116    #ifdef ALLOW_DIAGNOSTICS        
1117                DIAGarrayA(I,J) = ZERO
1118    #endif
1119    
1120    cif         another mess because the default behavior was changed between v1.109 and v1.111
1121    #ifdef FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
1122    c           The sum of ice growth tendencies from open water fluxes
1123    c           + any residual ocean-ice fluxes (weighted by the open
1124    c           water fraction).
1125    
1126    c           Using the McPhee parameterization, r_QbyOCN .LE ZERO
1127    c           Therefore, initial ice growth is triggered by open water
1128    c           heat flux divergence overcoming the melting tendency of
1129    c           the ocean-ice heat fluxes.
1130                tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1131         &            (1.0 _d 0 - AREApreTH(i,J))
1132    
1133    c           A fraction (SWFRACB) of the open-water heat flux includes
1134    c           shortwave radiative convergence at depths below the upper
1135    c           ocean grid cell.  These must be subtracted out.
1136                tmpscal2=SWFRACB * a_QSWbyATM_open(I,J)
1137    
1138    #ifdef FENTY_OPEN_WATER_FLUXES_MELT_ICE
1139    c           Convergent open water heat fluxes melt ice directly,
1140    c           bypassing the ocean
1141                tmpscal3 = tmpscal1-tmpscal2
1142    #else
1143    c           Convergent open water heat fluxes warm the ocean
1144    c           leading to an indirect ice melt.
1145                tmpscal3 = MAX(0.0 _d 0, tmpscal1-tmpscal2)
1146    #endif /* FENTY_OPEN_WATER_FLUXES_MELT_ICE */
1147    
1148                d_HEFFbyATMonOCN(I,J)=MAX(tmpscal3, -HEFF(I,J,bi,bj))
1149                r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-d_HEFFbyATMonOCN(I,J)
1150                HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyATMonOCN(I,J)
1151    
1152    #else /* FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */
1153    
1154    #ifdef SEAICE_DO_OPEN_WATER_GROWTH
1155            
1156  c           Initial ice growth is triggered by open water  c           Initial ice growth is triggered by open water
1157  c           heat flux overcoming potential melt by ocean  c           heat flux overcoming potential melt by ocean
1158              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
# Line 982  C           allow not only growth but al Line 1165  C           allow not only growth but al
1165  #else  #else
1166              tmpscal3=MAX(tmpscal1-tmpscal2, 0. _d 0)              tmpscal3=MAX(tmpscal1-tmpscal2, 0. _d 0)
1167  #endif /* SEAICE_DO_OPEN_WATER_MELT */  #endif /* SEAICE_DO_OPEN_WATER_MELT */
1168            d_HEFFbyATM_open(I,J)=tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
           d_HEFFbyATM(I,J)=d_HEFFbyATM(I,J)+tmpscal3  
1169            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1170            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1171              
1172    #endif /* SEAICE_DO_OPEN_WATER_GROWTH */
1173    #endif /* FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */
1174    
1175    #ifdef ALLOW_DIAGNOSTICS
1176                DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
1177         &           * d_HEFFbyATMonOCN(I,J)
1178    #endif
1179           ENDDO           ENDDO
1180          ENDDO          ENDDO
1181  #endif /* SEAICE_DO_OPEN_WATER_GROWTH */  
1182    #ifdef ALLOW_DIAGNOSTICS
1183            IF ( useDiagnostics ) THEN
1184             IF ( DIAGNOSTICS_IS_ON('SIdHbATO',myThid) ) THEN
1185              CALL DIAGNOSTICS_FILL(DIAGarrayA,
1186         &      'SIdHbATO',0,1,3,bi,bj,myThid)
1187             ENDIF
1188            ENDIF
1189    #endif /* ALLOW DIAGNOSTICS */
1190        
1191  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1192    
1193  C convert snow to ice if submerged.  C convert snow to ice if submerged.
# Line 1023  C ==========PART 4: determine ice cover Line 1222  C ==========PART 4: determine ice cover
1222  C ===================================================================  C ===================================================================
1223    
1224  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1225  CADJ STORE d_HEFFbyATM = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1226  CADJ STORE d_HEFFbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyATMonICE = comlev1_bibj,key=iicekey,byte=isbyte
1227  CADJ STORE d_HEFFbyATM_open=comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
 CADJ STORE d_HEFFbyATM_cover=comlev1_bibj,key=iicekey,byte=isbyte  
1228  CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1229  CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1230  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
# Line 1037  CADJ STORE AREA(:,:,bi,bj) = comlev1_bib Line 1235  CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
1235    
1236          DO J=1,sNy          DO J=1,sNy
1237           DO I=1,sNx           DO I=1,sNx
 C compute ice melt due to ATM (and OCN) heat stocks  
1238    
1239            IF (SEAICEareaFormula.EQ.1) THEN  C         compute ice melt due to ATM (and OCN) heat stocks
1240    #ifdef FENTY_AREA_EXPANSION_CONTRACTION
1241    
1242    #ifdef ALLOW_DIAGNOSTICS
1243              DIAGarrayA(I,J) = ZERO
1244              DIAGarrayB(I,J) = ZERO
1245              DIAGarrayC(I,J) = ZERO
1246              DIAGarray(I,J) = ZERO
1247    #endif
1248    
1249    c         The minimum effective thickness used in the ice contraction
1250    c         parameterization.
1251              heff_star = sqrt(HEFFpreTH(I,J)*HEFFPreTH(I,J) + 0.01 _d 0)
1252    
1253    c         the various thickness tendency terms
1254              tmpscal1 = d_HEFFbyATMOnOCN(I,J)
1255              tmpscal2 = d_HEFFbyATMOnICE(I,J)
1256              tmpscal3 = d_HEFFbyOCNonICE(I,J)
1257    
1258    c         Part 1: Treat expansion/contraction of ice cover in open-water areas
1259    
1260    C         All new ice cover is created from divergent air-sea heat fluxes.  Divergent
1261    c         air-sea heat fluxes must exceed the potential convergent ocean-ice heat fluxes
1262    c         for ice to form.  
1263              IF (tmpscal1 .GT. ZERO) then
1264                IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1265                 d_AREAbyATM(I,J)=tmpscal1/HO_south
1266                ELSE
1267                 d_AREAbyATM(I,J)=tmpscal1/HO
1268                ENDIF
1269    
1270    c         If there are convergent air-sea heat fluxes and convergent air-sea
1271    c         heat fluxes are permitted to melt ice, tmpscal1 < 0.
1272              ELSEIF (AREApreTH(I,J) .GT. ZERO) THEN
1273                 d_AREAbyATM(i,J) = HALF*tmpscal1*AREApreTH(I,J)/heff_star
1274              ENDIF
1275    
1276    c         Part 2: Reduce ice concentration from ice thinning from above
1277    
1278    c         Ice concentrations are reduced whenever existing ice thins from surface
1279    c         heat flux convergence.
1280              if (tmpscal2 .LE. ZERO) then
1281                 d_AREAbyICE(I,J) =
1282         &            HALF * tmpscal2 * AREApreTH(I,J)/heff_star
1283              endif
1284    
1285    c         Part 3: Reduce ice concentration from ice thinning from below
1286    
1287    c         Sensible heat transfer from the ocean to the sea ice thins it and
1288    c         reduces concentrations
1289              if (tmpscal3 .LE.ZERO) then
1290                 d_AREAbyOCN(I,J) =
1291         &         HALF * tmpscal3 * AREApreTH(I,J)/heff_star
1292              endif
1293    
1294    #ifdef ALLOW_DIAGNOSTICS
1295              DIAGarrayA(I,J) = d_AREAbyATM(I,J)
1296              DIAGarrayB(I,J) = d_AREAbyICE(I,J)
1297              DIAGarrayC(I,J) = d_AREAbyOCN(I,J)
1298              DIAGarray(I,J) = d_AREAbyICE(I,J) + d_AREAbyATM(I,J)
1299         &                     + d_AREAbyOCN(I,J)
1300    #endif
1301    
1302    #else /* FENTY_AREA_EXPANSION_CONTRACTION */
1303    
1304    #ifdef SEAICE_GROWTH_LEGACY
1305    
1306  C compute heff after ice melt by ocn:  C compute heff after ice melt by ocn:
1307            tmpscal0=HEFF(I,J,bi,bj)            tmpscal0=HEFF(I,J,bi,bj)
1308       &            - d_HEFFbyATM(I,J) - d_HEFFbyFLOODING(I,J)       &            - d_HEFFbyATMonOCN(I,J) - d_HEFFbyFLOODING(I,J)
1309         &            - d_HEFFbyATMonICE(I,J)
1310    
1311  C compute available heat left after snow melt by atm:  C compute available heat left after snow melt by atm:
1312            tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)            tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
1313       &            - d_HSNWbyATM(I,J)/ICE2SNOW       &            - d_HSNWbyATMonSNW(I,J)/ICE2SNOW
1314  C (cannot melt more than all the ice)  C (cannot melt more than all the ice)
1315            tmpscal2 = MAX(-tmpscal0,tmpscal1)            tmpscal2 = MAX(-tmpscal0,tmpscal1)
1316            tmpscal3 = MIN(ZERO,tmpscal2)            tmpscal3 = MIN(ZERO,tmpscal2)
# Line 1054  C (cannot melt more than all the ice) Line 1319  C (cannot melt more than all the ice)
1319  #endif  #endif
1320  C gain of new ice over open water  C gain of new ice over open water
1321            tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))            tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))
 C compute cover fraction tendency  
           IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN  
            d_AREA(I,J)=tmpscal4/HO_south  
           ELSE  
            d_AREA(I,J)=tmpscal4/HO  
           ENDIF  
           d_AREA(I,J)=d_AREA(I,J)  
      &         +HALF*tmpscal3*AREApreTH(I,J)  
      &         /(tmpscal0+.00001 _d 0)  
1322    
1323            ELSEIF (SEAICEareaFormula.EQ.2) THEN  #else /* SEAICE_GROWTH_LEGACY */
1324    
1325    # ifdef SEAICE_OCN_MELT_ACT_ON_AREA
1326  C ice cover reduction by joint OCN+ATM melt  C ice cover reduction by joint OCN+ATM melt
1327            tmpscal3 = MIN( 0. _d 0 ,d_HEFFbyATM(I,J)+d_HEFFbyOCN(I,J) )            tmpscal3 = MIN( 0. _d 0 ,
1328         &              d_HEFFbyATMonOCN(I,J)+d_HEFFbyOCNonICE(I,J)
1329         &              + d_HEFFbyATMonICE(I,J)                     )
1330    # else
1331    C ice cover reduction by ATM melt only -- as in legacy code
1332              tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN(I,J)
1333         &                            + d_HEFFbyATMonICE(I,J) )
1334    # endif
1335  C gain of new ice over open water  C gain of new ice over open water
1336    
1337  # ifdef SEAICE_DO_OPEN_WATER_GROWTH  # ifdef SEAICE_DO_OPEN_WATER_GROWTH
1338  C the one effectively used to increment HEFF  C the one effectively used to increment HEFF
1339            tmpscal4 = MAX(ZERO,d_HEFFbyATM_open(I,J))            tmpscal4 = d_HEFFbyATMonOCN(I,J)
1340  # else  # else
1341  C the virtual one -- as in legcy code  C the virtual one -- as in legcy code
1342            tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))            tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))
1343  # endif  # endif
1344    #endif /* SEAICE_GROWTH_LEGACY */
1345    
1346  C compute cover fraction tendency  C compute cover fraction tendency
1347            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1348             d_AREA(I,J)=tmpscal4/HO_south             d_AREAbyATM(I,J)=tmpscal4/HO_south
1349            ELSE            ELSE
1350             d_AREA(I,J)=tmpscal4/HO             d_AREAbyATM(I,J)=tmpscal4/HO
1351            ENDIF            ENDIF
1352            d_AREA(I,J)=d_AREA(I,J)            d_AREAbyATM(I,J)=d_AREAbyATM(I,J)
1353       &         +HALF*tmpscal3/heffActual(I,J)  #ifdef SEAICE_GROWTH_LEGACY
1354         &         +HALF*tmpscal3*AREApreTH(I,J)
1355            ELSEIF (SEAICEareaFormula.EQ.3) THEN       &         /(tmpscal0+.00001 _d 0)
 c         a) extend/reduce ice cover due to open-water ice growth/melt  
 # ifdef SEAICE_DO_OPEN_WATER_GROWTH  
 C the one effectively used to increment HEFF  
           tmpscal1 = d_HEFFbyATM_open(I,J)  
1356  #else  #else
1357  C the virtual one -- as in legcy code       &         +HALF*tmpscal3/heffActual(I,J)
           tmpscal1 = a_QbyATM_open(I,J)  
1358  #endif  #endif
           IF (tmpscal1 .GT. ZERO) then  
             IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN  
              d_AREAbyATM_open(I,J)=tmpscal1/HO_south  
             ELSE  
              d_AREAbyATM_open(I,J)=tmpscal1/HO  
             ENDIF  
           ELSE  
              d_AREAbyATM_open(i,J) = HALF*tmpscal1/heffActual(I,J)  
           ENDIF  
 c         b) Reduce ice concentration due to ice thinning by atmospheric heat  
           d_AREAbyATM_cover(I,J) = MIN( 0. _d 0 ,  
      &           HALF * d_HEFFbyATM_cover(I,J) /heffActual(I,J) )  
 c         c) Reduce ice concentration due to ice thinning by ocean heat  
           d_AREAbyOCN(I,J) = MIN( 0. _d 0 ,  
      &           HALF * d_HEFFbyOCN(I,J) /heffActual(I,J) )  
 c         d) add up the various terms  
           d_AREA(I,J) = d_AREAbyATM_open(I,J)  
      &                + d_AREAbyATM_cover(I,J)  
      &                + d_AREAbyOCN(I,J)  
1359    
1360            ENDIF  #endif /* FENTY_AREA_EXPANSION_CONTRACTION */
1361    
1362  C apply tendency  C apply tendency
1363            IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.            IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
1364       &        (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN       &        (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
1365             AREA(I,J,bi,bj)=max(0. _d 0 , min( 1. _d 0,             AREA(I,J,bi,bj)=max(0. _d 0 , min( 1. _d 0,
1366       &                     AREA(I,J,bi,bj)+d_AREA(I,J) ) )       &                     AREA(I,J,bi,bj)+d_AREAbyATM(I,J)
1367    #ifdef FENTY_AREA_EXPANSION_CONTRACTION
1368         &                   + d_AREAbyOCN(I,J) + d_AREAbyICE(I,J)
1369    #endif
1370         &                                       ))
1371            ELSE            ELSE
1372             AREA(I,J,bi,bj)=0. _d 0             AREA(I,J,bi,bj)=0. _d 0
1373            ENDIF            ENDIF
1374           ENDDO           ENDDO
1375          ENDDO          ENDDO
1376    
1377    #ifdef ALLOW_DIAGNOSTICS
1378            IF ( useDiagnostics ) THEN
1379             IF ( DIAGNOSTICS_IS_ON('SIdA    ',myThid) ) THEN
1380              CALL DIAGNOSTICS_FILL(DIAGarray,
1381         &      'SIdA    ',0,1,3,bi,bj,myThid)
1382             ENDIF
1383             IF ( DIAGNOSTICS_IS_ON('SIdAbATO',myThid) ) THEN
1384              CALL DIAGNOSTICS_FILL(DIAGarrayA,
1385         &      'SIdAbATO',0,1,3,bi,bj,myThid)
1386             ENDIF
1387             IF ( DIAGNOSTICS_IS_ON('SIdAbATC',myThid) ) THEN
1388              CALL DIAGNOSTICS_FILL(DIAGarrayB,
1389         &      'SIdAbATC',0,1,3,bi,bj,myThid)
1390             ENDIF
1391             IF ( DIAGNOSTICS_IS_ON('SIdAbOCN',myThid) ) THEN
1392              CALL DIAGNOSTICS_FILL(DIAGarrayC,
1393         &      'SIdAbOCN',0,1,3,bi,bj,myThid)
1394             ENDIF
1395            ENDIF
1396    #endif
1397    
1398  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1399  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1400  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
# Line 1160  C ====================================== Line 1429  C ======================================
1429  # ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
1430  #  ifdef ALLOW_SALT_PLUME  #  ifdef ALLOW_SALT_PLUME
1431  CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1432  CADJ STORE d_HEFFbyOCN = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1433  CADJ STORE d_HEFFbyATM = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1434    CADJ STORE d_HEFFbyATMonICE = comlev1_bibj,key=iicekey,byte=isbyte
1435  CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
1436  CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,  CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1437  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
# Line 1170  CADJ &                       key = iicek Line 1440  CADJ &                       key = iicek
1440          DO J=1,sNy          DO J=1,sNy
1441           DO I=1,sNx           DO I=1,sNx
1442  Cgf note: flooding does count negatively  Cgf note: flooding does count negatively
1443            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCN(I,J) +            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
1444       &               d_HEFFbyATM(I,J) - d_HEFFbyFLOODING(I,J)       &               d_HEFFbyATMonICE(I,J) +
1445         &               d_HEFFbyATMonOCN(I,J) - d_HEFFbyFLOODING(I,J)
1446            tmpscal2 = tmpscal1 * SIsal0 * HEFFM(I,J,bi,bj)            tmpscal2 = tmpscal1 * SIsal0 * HEFFM(I,J,bi,bj)
1447       &            /SEAICE_deltaTtherm * ICE2WATR * rhoConstFresh       &            /SEAICE_deltaTtherm * ICE2WATR * rhoConstFresh
1448            saltFlux(I,J,bi,bj) = tmpscal2            saltFlux(I,J,bi,bj) = tmpscal2
# Line 1210  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bi Line 1481  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bi
1481          DO J=1,sNy          DO J=1,sNy
1482           DO I=1,sNx           DO I=1,sNx
1483  C sum up the terms that affect the salt content of the ice pack  C sum up the terms that affect the salt content of the ice pack
1484           tmpscal1=d_HEFFbyOCN(I,J)+d_HEFFbyATM(I,J)           tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)+
1485         &       d_HEFFbyATMonICE(I,J)
1486    
1487  C recompute HEFF before thermodyncamic updates (which is not AREApreTH in legacy code)  C recompute HEFF before thermodyncamic updates (which is not AREApreTH in legacy code)
1488           tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)           tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
1489  C tmpscal1 > 0 : m of sea ice that is created  C tmpscal1 > 0 : m of sea ice that is created
# Line 1437  C ====================================== Line 1710  C ======================================
1710          DO J=1,sNy          DO J=1,sNy
1711           DO I=1,sNx           DO I=1,sNx
1712            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
1713       &         - ( d_HEFFbyOCN(I,J) +  #ifdef FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
1714       &             d_HSNWbyOCN(I,J)/ICE2SNOW +       &         +   a_QSWbyATM_cover(I,J)
1715    #endif /* FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */          
1716         &         - ( d_HEFFbyOCNonICE(I,J) +
1717         &             d_HSNWbyOCNonSNW(I,J)/ICE2SNOW +
1718       &             d_HEFFbyNEG(I,J) +       &             d_HEFFbyNEG(I,J) +
1719       &             d_HSNWbyNEG(I,J)/ICE2SNOW )       &             d_HSNWbyNEG(I,J)/ICE2SNOW )
1720       &         * maskC(I,J,kSurface,bi,bj)       &         * maskC(I,J,kSurface,bi,bj)
# Line 1484  C ====================================== Line 1760  C ======================================
1760  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
1761          DO J=1,sNy          DO J=1,sNy
1762           DO I=1,sNx           DO I=1,sNx
1763            tmpscal1= d_HSNWbyATM(I,J)/ICE2SNOW            tmpscal1= d_HSNWbyATMonSNW(I,J)/ICE2SNOW
1764       &             +d_HFRWbyRAIN(I,J)       &             +d_HFRWbyRAIN(I,J)
1765       &             +d_HSNWbyOCN(I,J)/ICE2SNOW       &             +d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
1766       &             +d_HEFFbyOCN(I,J)       &             +d_HEFFbyOCNonICE(I,J)
1767       &             +d_HEFFbyATM(I,J)       &             +d_HEFFbyATMonOCN(I,J)
1768         &             +d_HEFFbyATMonICE(I,J)    
1769       &             +d_HEFFbyNEG(I,J)       &             +d_HEFFbyNEG(I,J)
1770       &             +d_HSNWbyNEG(I,J)/ICE2SNOW       &             +d_HSNWbyNEG(I,J)/ICE2SNOW
1771  #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET  #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
Legend:
Removed from v.1.112  
changed lines
  Added in v.1.113
  ViewVC Help
Powered by ViewVC 1.1.22