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

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

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

-revision 1.9 by torge,
Mon Oct 22 16:36:45 2012 UTC
+revision 1.10 by torge,
Mon Oct 22 21:12:47 2012 UTC
 Line 97 
 C     i,j,bi,bj :: Loop counters
        INTEGER i, j, bi, bj
  C     number of surface interface layer
        INTEGER kSurface
+ C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
+       INTEGER IT
+       _RL pFac
  C     constants
        _RL tempFrz, ICE2SNOW, SNOW2ICE
        _RL QI, QS, recip_QI
+       _RL lhSublim
+ C conversion factors to go from Q (W/m2) to HEFF (ice meters)
+       _RL convertQ2HI, convertHI2Q
+ C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
+       _RL convertPRECIP2HI, convertHI2PRECIP
+ C     Factor by which we increase the upper ocean friction velocity (u*) when
+ C     ice is absent in a grid cell  (dimensionless)
+       _RL MixedLayerTurbulenceFactor
+ C     wind speed square
+       _RL SPEED_SQ
+ C     Regularization values squared
+       _RL area_reg_sq, hice_reg_sq
+ C     pathological cases thresholds
+       _RL heffTooHeavy
+ C     Helper variables: reciprocal of some constants
+       _RL recip_multDim
+       _RL recip_deltaTtherm
+       _RL recip_rhoIce
+ C     local value (=1/HO or 1/HO_south)
+       _RL recip_HO
+ C     local value (=1/ice thickness)
+       _RL recip_HH
+ C     facilitate multi-category snow implementation
+       _RL pFacSnow
+ C     temporary variables available for the various computations
+       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
+ #ifdef ALLOW_SITRACER
+       INTEGER iTr
+ #ifdef ALLOW_DIAGNOSTICS
+       CHARACTER*8   diagName
+ #endif
+ #endif /* ALLOW_SITRACER */
+ #ifdef ALLOW_AUTODIFF_TAMC
+       INTEGER ilockey
+ #endif
+ C==   local arrays ==
  C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
        _RL TmixLoc       (1:sNx,1:sNy)
+ C     actual ice thickness (with upper and lower limit)
+       _RL heffActual          (1:sNx,1:sNy)
+ C     actual snow thickness
+       _RL hsnowActual         (1:sNx,1:sNy)
+ C     actual ice thickness (with lower limit only) Reciprocal
+       _RL recip_heffActual    (1:sNx,1:sNy)
+ C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
+       _RL AREApreTH           (1:sNx,1:sNy)
+       _RL HEFFpreTH           (1:sNx,1:sNy)
+       _RL HSNWpreTH           (1:sNx,1:sNy)
+ #ifdef SEAICE_ITD
+       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
+       _RL leadIceThickMin
+ #endif
+ C     wind speed
+       _RL UG                  (1:sNx,1:sNy)
+ C     temporary variables available for the various computations
+       _RL tmparr1             (1:sNx,1:sNy)
+ #ifdef SEAICE_VARIABLE_SALINITY
+       _RL saltFluxAdjust      (1:sNx,1:sNy)
+ #endif
+       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
+       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
+       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
+       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
+ #endif
+       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
+       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #endif
  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
  C             the atmosphere and the ocean surface - for ice covered water
  C     a_QbyATM_open  :: same but for open water
-Line 124 
 C             processes have been accoun
+Line 212 
 C             processes have been accoun
        _RL a_QbyOCN            (1:sNx,1:sNy)
        _RL r_QbyOCN            (1:sNx,1:sNy)
- C conversion factors to go from Q (W/m2) to HEFF (ice meters)
+ C     The change of mean ice thickness due to out-of-bounds values following
-       _RL convertQ2HI, convertHI2Q
+ C     sea ice dyhnamics
- C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
-       _RL convertPRECIP2HI, convertHI2PRECIP
- #ifdef ALLOW_DIAGNOSTICS
- C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
-       _RL d_AREAbyATM         (1:sNx,1:sNy)
-       _RL d_AREAbyOCN         (1:sNx,1:sNy)
-       _RL d_AREAbyICE         (1:sNx,1:sNy)
- #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
- C     ICE/SNOW stocks tendency associated with relaxation towards observation
-       _RL d_AREAbyRLX         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to relaxation
-       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
- #endif
- c     The change of mean ice thickness due to out-of-bounds values following
- c     sea ice dynamics
        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
+ C     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
- c     The sum of mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The sum of mean ice thickness increments due to atmospheric fluxes over
- c     fraction and ice-covered fractions of the grid cell
+ C     the open water fraction and ice-covered fractions of the grid cell
        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
- c     The change of mean ice thickness due to flooding by snow
+ C     The change of mean ice thickness due to flooding by snow
        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
- c     The mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The mean ice thickness increments due to atmospheric fluxes over the open
- c     fraction and ice-covered fractions of the grid cell, respectively
+ C     water fraction and ice-covered fractions of the grid cell, respectively
        _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
        _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
-Line 167 
 c     fraction and ice-covered fractions
+Line 236 
 c     fraction and ice-covered fractions
        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
- C
  C     a_FWbySublim :: fresh water flux implied by latent heat of
  C                     sublimation to atmosphere, same sign convention
  C                     as EVAP (positive upward)
-Line 180 
 C                     as EVAP (positive
+Line 249 
 C                     as EVAP (positive
  C     The latent heat flux which will sublimate all snow and ice
  C     over one time step
        _RL latentHeatFluxMax   (1:sNx,1:sNy)
-       _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
+       _RL latentHeatFluxMaxMult(1:sNx,1:sNy,MULTDIM)
  #endif
- C     actual ice thickness (with upper and lower limit)
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-       _RL heffActual          (1:sNx,1:sNy)
+ C     ICE/SNOW stocks tendency associated with relaxation towards observation
- C     actual snow thickness
+       _RL d_AREAbyRLX         (1:sNx,1:sNy)
-       _RL hsnowActual         (1:sNx,1:sNy)
+ C     The change of mean ice thickness due to relaxation
- C     actual ice thickness (with lower limit only) Reciprocal
+       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
-       _RL recip_heffActual    (1:sNx,1:sNy)
- C     local value (=1/HO or 1/HO_south)
-       _RL recip_HO
- C     local value (=1/ice thickness)
-       _RL recip_HH
- C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
-       _RL AREApreTH           (1:sNx,1:sNy)
-       _RL HEFFpreTH           (1:sNx,1:sNy)
-       _RL HSNWpreTH           (1:sNx,1:sNy)
- #ifdef SEAICE_ITD
-       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
-       _RL leadIceThickMin
- #endif
- C     wind speed
-       _RL UG                  (1:sNx,1:sNy)
- #ifdef ALLOW_ATM_WIND
-       _RL SPEED_SQ
- #endif
- C     Regularization values squared
-       _RL area_reg_sq, hice_reg_sq
- C     pathological cases thresholds
-       _RL heffTooHeavy
-       _RL lhSublim
- C temporary variables available for the various computations
-       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
-       _RL tmparr1             (1:sNx,1:sNy)
- #ifdef SEAICE_VARIABLE_SALINITY
-       _RL saltFluxAdjust      (1:sNx,1:sNy)
- #endif
-       INTEGER ilockey
-       INTEGER it
-       _RL pFac
-       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
-       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
-       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
- #ifdef SEAICE_ITD
-       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
-       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
- #endif
-       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
-       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
-       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
- #ifdef SEAICE_ITD
-       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
-       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
  #endif
- C     Helper variables: reciprocal of some constants
-       _RL recip_multDim
-       _RL recip_deltaTtherm
-       _RL recip_rhoIce
- C     Factor by which we increase the upper ocean friction velocity (u*) when
- C     ice is absent in a grid cell  (dimensionless)
-       _RL MixedLayerTurbulenceFactor
- #ifdef ALLOW_SITRACER
-       INTEGER iTr
-       CHARACTER*8   diagName
- #endif
  #ifdef ALLOW_DIAGNOSTICS
- c     Helper variables for diagnostics
+ C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
+       _RL d_AREAbyATM         (1:sNx,1:sNy)
+       _RL d_AREAbyOCN         (1:sNx,1:sNy)
+       _RL d_AREAbyICE         (1:sNx,1:sNy)
+ C     Helper variables for diagnostics
        _RL DIAGarrayA    (1:sNx,1:sNy)
        _RL DIAGarrayB    (1:sNx,1:sNy)
        _RL DIAGarrayC    (1:sNx,1:sNy)
        _RL DIAGarrayD    (1:sNx,1:sNy)
- #endif
+ #endif /* ALLOW_DIAGNOSTICS */
+       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmQnt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+ #ifdef ALLOW_BALANCE_FLUXES
+       _RL FWFsiTile(nSx,nSy)
+       _RL FWFsiGlob
+       _RL HFsiTile(nSx,nSy)
+       _RL HFsiGlob
+       _RL FWF2HFsiTile(nSx,nSy)
+       _RL FWF2HFsiGlob
+       CHARACTER*(max_len_mbuf) msgbuf
+ #endif
  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-Line 335 
 C conversion factors to go from precip (
+Line 352 
 C conversion factors to go from precip (
       &                       + act4*max1*max2*max3
  #endif /* ALLOW_AUTODIFF_TAMC */
  C array initializations
  C =====================
-Line 358 
 C =====================
+Line 374 
 C =====================
            d_AREAbyOCN(I,J)           = 0.0 _d 0
  #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
            d_AREAbyRLX(I,J)           = 0.0 _d 0
            d_HEFFbyRLX(I,J)           = 0.0 _d 0
  #endif
-Line 382 
 C =====================
+Line 398 
 C =====================
  #ifdef SEAICE_CAP_SUBLIM
            latentHeatFluxMax(I,J)     = 0.0 _d 0
  #endif
- c
            d_HFRWbyRAIN(I,J)          = 0.0 _d 0
            tmparr1(I,J)               = 0.0 _d 0
  #ifdef SEAICE_VARIABLE_SALINITY
            saltFluxAdjust(I,J)        = 0.0 _d 0
  #endif
-Line 414 
 c
+Line 427 
 c
          ENDDO
  #endif
  C =====================================================================
  C ===========PART 1: treat pathological cases (post advdiff)===========
  C =====================================================================
-Line 455 
 Cgf no dependency through pathological c
+Line 467 
 Cgf no dependency through pathological c
          IF ( SEAICEadjMODE.EQ.0 ) THEN
  #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  C 0) relax sea ice concentration towards observation
-Line 485 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 497 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
           ENDDO
          ENDDO
          ENDIF
- #endif /* SEAICE_ALLOW_AREA_RELAXATION */
+ #endif /* EXF_ALLOW_SEAICE_RELAX */
  C 1) treat the case of negative values:
-Line 881 
 cif       Do not regularize when HEFFpre
+Line 893 
 cif       Do not regularize when HEFFpre
              tmpscal2         = HEFFpreTH(I,J)/tmpscal1
              heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)
  #else /* SEAICE_GROWTH_LEGACY */
- cif        regularize AREA with SEAICE_area_reg
+ Cif        regularize AREA with SEAICE_area_reg
             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
- cif        heffActual calculated with the regularized AREA
+ Cif        heffActual calculated with the regularized AREA
             tmpscal2 = HEFFpreTH(I,J) / tmpscal1
- cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
+ Cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
- cif        hsnowActual calculated with the regularized AREA
+ Cif        hsnowActual calculated with the regularized AREA
             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
  #endif /* SEAICE_GROWTH_LEGACY */
- cif        regularize the inverse of heffActual by hice_reg
+ Cif        regularize the inverse of heffActual by hice_reg
             recip_heffActual(I,J)  = AREApreTH(I,J) /
       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
- cif       Do not regularize when HEFFpreTH = 0
+ Cif       Do not regularize when HEFFpreTH = 0
            ELSE
              heffActual(I,J) = ZERO
              hsnowActual(I,J) = ZERO
-Line 921 
 C    AND SNOW THICKNESS
+Line 933 
 C    AND SNOW THICKNESS
  #endif
          DO J=1,sNy
           DO I=1,sNx
- c        The latent heat flux over the sea ice which
+ C        The latent heat flux over the sea ice which
- c        will sublimate all of the snow and ice over one time
+ C        will sublimate all of the snow and ice over one time
- c        step (W/m^2)
+ C        step (W/m^2)
  #ifdef SEAICE_ITD
            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
              latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
-Line 981 
 cCADJ STORE TmixLoc = comlev1_bibj, key
+Line 993 
 cCADJ STORE TmixLoc = comlev1_bibj, key
  C determine available heat due to the atmosphere -- for ice covered water
  C =======================================================================
- #ifdef ALLOW_ATM_WIND
+         IF (useRelativeWind.AND.useAtmWind) THEN
-         IF (useRelativeWind) THEN
  C     Compute relative wind speed over sea ice.
           DO J=1,sNy
            DO I=1,sNx
-Line 1003 
 C     Compute relative wind speed over s
+Line 1014 
 C     Compute relative wind speed over s
            ENDDO
           ENDDO
          ENDIF
- #endif /* ALLOW_ATM_WIND */
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE tice(:,:,bi,bj)
-Line 1022 
 C--   Start loop over multi-categories
+Line 1032 
 C--   Start loop over multi-categories
  CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
  #endif
          DO IT=1,SEAICE_multDim
- c        homogeneous distribution between 0 and 2 x heffActual
+ C        homogeneous distribution between 0 and 2 x heffActual
  #ifndef SEAICE_ITD
-          pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
+          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
+          pFacSnow = 1. _d 0
+          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
  #endif
           DO J=1,sNy
            DO I=1,sNx
-Line 1032 
 c        homogeneous distribution betwee
+Line 1044 
 c        homogeneous distribution betwee
  CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
  C    (instead of heffActual and latentHeatFluxMax)
             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
+            hsnowActualMult(I,J,IT)=hsnowActual(I,J)*pFacSnow
  #ifdef SEAICE_CAP_SUBLIM
             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
  #endif
-Line 1046 
 C    (instead of heffActual and latentHe
+Line 1059 
 C    (instead of heffActual and latentHe
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
+ CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
  CADJ STORE ticeInMult     = comlev1_bibj, key = iicekey, byte = isbyte
  # ifdef SEAICE_CAP_SUBLIM
  CADJ STORE latentHeatFluxMaxMult
-Line 1061 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1075 
 CADJ &     comlev1_bibj, key = iicekey,
          DO IT=1,SEAICE_multDim
           CALL SEAICE_SOLVE4TEMP(
- #ifdef SEAICE_ITD
       I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
- #else
-      I        UG, heffActualMult(1,1,IT), hsnowActual,
- #endif
  #ifdef SEAICE_CAP_SUBLIM
       I        latentHeatFluxMaxMult(1,1,IT),
  #endif
       U        ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
-      O        a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),
+      O        a_QbyATMmult_cover(1,1,IT),
+      O        a_QSWbyATMmult_cover(1,1,IT),
       O        a_FWbySublimMult(1,1,IT),
       I        bi, bj, myTime, myIter, myThid )
          ENDDO
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
+ CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
  CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte
  # ifdef SEAICE_CAP_SUBLIM
  CADJ STORE latentHeatFluxMaxMult
-Line 1134 
 C     (fluxes are per unit (ice surface)
+Line 1146 
 C     (fluxes are per unit (ice surface)
  # ifdef ALLOW_DIAGNOSTICS
          DO J=1,sNy
           DO I=1,sNx
- c          The actual latent heat flux realized by SOLVE4TEMP
+ C          The actual latent heat flux realized by SOLVE4TEMP
             DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
           ENDDO
          ENDDO
- cif     The actual vs. maximum latent heat flux
+ Cif     The actual vs. maximum latent heat flux
          IF ( useDiagnostics ) THEN
            CALL DIAGNOSTICS_FILL(DIAGarrayA,
       &     'SIactLHF',0,1,3,bi,bj,myThid)
-Line 1206 
 C and initialize r_QbyATM_cover/r_QbyATM
+Line 1218 
 C and initialize r_QbyATM_cover/r_QbyATM
  C     Convert fresh water flux by sublimation to 'effective' ice meters.
  C     Negative sublimation is resublimation and will be added as snow.
  #ifdef SEAICE_DISABLE_SUBLIM
- cgf just for those who may need to omit this term to reproduce old results
+ Cgf just for those who may need to omit this term to reproduce old results
            a_FWbySublim(I,J) = ZERO
- #endif
+ #endif /* SEAICE_DISABLE_SUBLIM */
            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
            r_FWbySublim(I,J)=a_FWbySublim(I,J)
-Line 1266 
 CADJ &                       key = iicek
+Line 1278 
 CADJ &                       key = iicek
          DO J=1,sNy
           DO I=1,sNx
- c         FREEZING TEMP. OF SEA WATER (deg C)
+ C         FREEZING TEMP. OF SEA WATER (deg C)
            tempFrz = SEAICE_tempFrz0 +
       &              SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
- c efficiency of turbulent fluxes : dependency to sign of THETA-TBC
+ C efficiency of turbulent fluxes : dependency to sign of THETA-TBC
            IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
             tmpscal1 = SEAICE_mcPheePiston
            ELSE
             tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
            ENDIF
- c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
+ C efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
            IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
       &         (.NOT.SEAICE_mcPheeStepFunc) ) THEN
             MixedLayerTurbulenceFactor = ONE -
-Line 1286 
 c efficiency of turbulent fluxes : depen
+Line 1298 
 c efficiency of turbulent fluxes : depen
            ELSE
             MixedLayerTurbulenceFactor = ONE
            ENDIF
- c maximum turbulent flux, in ice meters
+ C maximum turbulent flux, in ice meters
            tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
       &         * (theta(I,J,kSurface,bi,bj)-tempFrz)
       &         * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
- c available turbulent flux
+ C available turbulent flux
            a_QbyOCN(i,j) =
       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
            r_QbyOCN(i,j) = a_QbyOCN(i,j)
-Line 1301 
 c available turbulent flux
+Line 1313 
 c available turbulent flux
          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
  #endif
  C ===================================================================
  C =========PART 3: determine effective thicknesses increments========
  C ===================================================================
-Line 1578 
 c         Limit ice growth by potential
+Line 1589 
 c         Limit ice growth by potential
            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
  #else
            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
- c         Limit ice growth by potential melt by ocean
+ C         Limit ice growth by potential melt by ocean
       &        AREApreTH(I,J) * r_QbyOCN(I,J))
  #endif /* SEAICE_GROWTH_LEGACY */
-Line 1590 
 c         Limit ice growth by potential
+Line 1601 
 c         Limit ice growth by potential
  #ifdef ALLOW_SITRACER
            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
  #endif
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
  #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-Line 1613 
 c ToM<<< debug seaice_growth
+Line 1624 
 c ToM<<< debug seaice_growth
  c ToM>>>
  #endif
- C attribute precip to fresh water or snow stock,
+ C add snow precipitation to HSNOW.
- C depending on atmospheric conditions.
  C =================================================
  #ifdef ALLOW_ATM_TEMP
- #ifdef ALLOW_AUTODIFF_TAMC
+ # ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
  CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
- #endif /* ALLOW_AUTODIFF_TAMC */
+ # endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
+         IF ( snowPrecipFile .NE. ' ' ) THEN
-          DO I=1,sNx
+ C add snowPrecip to HSNOW
+          DO J=1,sNy
+           DO I=1,sNx
+            d_HSNWbyRAIN(I,J) = convertPRECIP2HI * ICE2SNOW *
+      &          snowPrecip(i,j,bi,bj) * AREApreTH(I,J)
+            d_HFRWbyRAIN(I,J) = -convertPRECIP2HI *
+      &          ( PRECIP(I,J,bi,bj) - snowPrecip(I,J,bi,bj) ) *
+      &          AREApreTH(I,J)
+            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
+           ENDDO
+          ENDDO
+         ELSE
+ C attribute precip to fresh water or snow stock,
+ C depending on atmospheric conditions.
+          DO J=1,sNy
+           DO I=1,sNx
  C possible alternatives to the a_QbyATM_cover criterium
  c          IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
  c          IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
-           IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
+            IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
  C           add precip as snow
              d_HFRWbyRAIN(I,J)=0. _d 0
              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
-           ELSE
+            ELSE
  C           add precip to the fresh water bucket
              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
              d_HSNWbyRAIN(I,J)=0. _d 0
-           ENDIF
+            ENDIF
           ENDDO
          ENDDO
  #ifdef SEAICE_ITD
-Line 1659 
 C           add precip to the fresh wate
+Line 1684 
 C           add precip to the fresh wate
  Cgf note: this does not affect air-sea heat flux,
  Cgf since the implied air heat gain to turn
  Cgf rain to snow is not a surface process.
+ C end of IF statement snowPrecipFile:
+         ENDIF
  #endif /* ALLOW_ATM_TEMP */
  #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-Line 1758 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
+Line 1785 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
- c           Initial ice growth is triggered by open water
+ C           Initial ice growth is triggered by open water
- c           heat flux overcoming potential melt by ocean
+ C           heat flux overcoming potential melt by ocean
              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
       &            (1.0 _d 0 - AREApreTH(I,J))
- c           Penetrative shortwave flux beyond first layer
+ C           Penetrative shortwave flux beyond first layer
- c           that is therefore not available to ice growth/melt
+ C           that is therefore not available to ice growth/melt
              tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
-Line 1819 
 c needs to be here to allow use also wit
+Line 1846 
 c needs to be here to allow use also wit
  #else
          DO J=1,sNy
           DO I=1,sNx
- c needs to be here to allow use also with LEGACY branch
+ C needs to be here to allow use also with LEGACY branch
            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
           ENDDO
          ENDDO
-Line 1885 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1912 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
            ENDDO
           ENDDO
  #endif
          ENDIF
  #endif /* SEAICE_GROWTH_LEGACY */
-Line 2088 
 CADJ &                       key = iicek
+Line 2115 
 CADJ &                       key = iicek
            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
       &               d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
       &             + d_HEFFbySublim(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-                    + d_HEFFbyRLX(I,J)
+      &             + d_HEFFbyRLX(I,J)
  #endif
            tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
       &            * recip_deltaTtherm * SEAICE_rhoIce
-Line 2183 
 C set HSALT = 0 if HEFF = 0 and compute
+Line 2210 
 C set HSALT = 0 if HEFF = 0 and compute
          ENDDO
  #endif /* SEAICE_VARIABLE_SALINITY */
  C =======================================================================
  C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
  C =======================================================================
-Line 2266 
 C =================================
+Line 2292 
 C =================================
  #ifdef ALLOW_SITRACER
          DO J=1,sNy
           DO I=1,sNx
- c needs to be here to allow use also with LEGACY branch
+ C needs to be here to allow use also with LEGACY branch
            SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
           ENDDO
          ENDDO
-Line 2329 
 C in principle a_QSWbyATM_cover should a
+Line 2355 
 C in principle a_QSWbyATM_cover should a
  C for backward compatibility it is left out of the LEGACY branch
       &         +   a_QSWbyATM_cover(I,J)
  #endif /* SEAICE_GROWTH_LEGACY */
-      &         - ( d_HEFFbyOCNonICE(I,J) +
+      &         - ( d_HEFFbyOCNonICE(I,J)
-      &             d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +
+      &           + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
-      &             d_HEFFbyNEG(I,J) +
+      &           + d_HEFFbyNEG(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-      &             d_HEFFbyRLX(I,J) +
+      &           + d_HEFFbyRLX(I,J)
  #endif
-      &             d_HSNWbyNEG(I,J)*SNOW2ICE )
+      &           + d_HSNWbyNEG(I,J)*SNOW2ICE
-      &         * maskC(I,J,kSurface,bi,bj)
+      &           - convertPRECIP2HI *
+      &             snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J))
+      &           ) * maskC(I,J,kSurface,bi,bj)
+          ENDDO
+         ENDDO
+         DO J=1,sNy
+          DO I=1,sNx
            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
           ENDDO
          ENDDO
-Line 2364 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
+Line 2396 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
  CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
  CADJ &                       key = iicekey, byte = isbyte
  # endif /* ALLOW_AUTODIFF_TAMC */
-       IF ( SEAICEheatConsFix ) THEN
+ cgf Unlike for evap and precip, the temperature of gained/lost
- c Unlike for evap and precip, the temperature of gained/lost
+ C ocean liquid water due to melt/freeze of solid water cannot be chosen
- c ocean liquid water due to melt/freeze of solid water cannot be chosen
+ C arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
- c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates
+ C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
- c on external_forcing_surf.F and applies the correction to QNET.
+ C at a different temperature, it leads to a loss of conservation in the
-       IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN
+ C ocean+ice system. While this is mostly a serious issue in the
- c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and
+ C real fresh water + non linear free surface framework, a mismatch
- c the small error of the synchronous time stepping case (see external_forcing_surf.F).
+ C between ice and ocean boundary condition can result in all cases.
+ C Below we therefore anticipate on external_forcing_surf.F
+ C to diagnoze and/or apply the correction to QNET.
          DO J=1,sNy
           DO I=1,sNx
- #ifdef ALLOW_DIAGNOSTICS
+ C ocean water going to ice/snow, in precip units
- c store unaltered QNET for diagnostic purposes
+            tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*(
-            DIAGarrayA(I,J)=QNET(I,J,bi,bj)
- #endif
- c compute the ocean water going to ice/snow, in precip units
-            tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
       &       ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
       &       + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
       &       + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
       &       + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
       &       * convertHI2PRECIP
- c factor in the heat content that external_forcing_surf.F
+      &       - snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J)) )
- c will associate with EMPMR, and remove it from QNET, so that
+ C factor in the heat content as done in external_forcing_surf.F
- c melt/freez water is in effect consistently gained/lost at 0degC
+       IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
-            IF (temp_EvPrRn.NE.UNSET_RL) THEN
+      &     .AND.(nonlinFreeSurf.NE.0) ) THEN
-              QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
+              tmpscal1 = - tmpscal3*
       &         HeatCapacity_Cp * temp_EvPrRn
-            ELSE
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
-              QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
+      &         .AND.(nonlinFreeSurf.NE.0) ) THEN
+              tmpscal1 = - tmpscal3*
       &         HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
-            ENDIF
+       ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+            tmpscal1 = - tmpscal3*HeatCapacity_Cp*
+      &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
+            tmpscal1 = ZERO
+       ENDIF
  #ifdef ALLOW_DIAGNOSTICS
- c back out the eventual TFLUX adjustement and fill diag
+ C in all cases, diagnoze the boundary condition mismatch to SIaaflux
-            DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)
+            DIAGarrayA(I,J)=tmpscal1
  #endif
+ C remove the mismatch when real fresh water is exchanged (at 0degC here)
+       IF ( useRealFreshWaterFlux.AND.(nonlinFreeSurf.GT.0).AND.
+      &   SEAICEheatConsFix ) QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+tmpscal1
           ENDDO
          ENDDO
-       ENDIF
  #ifdef ALLOW_DIAGNOSTICS
          CALL DIAGNOSTICS_FILL(DIAGarrayA,
       &      'SIaaflux',0,1,3,bi,bj,myThid)
  #endif
-       ENDIF
  #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
+ C compute the net heat flux, incl. adv. by water, entering ocean+ice
+ C ===================================================================
+          DO J=1,sNy
+           DO I=1,sNx
+ cgf 1) SIatmQnt (analogous to qnet; excl. adv. by water exch.)
+ CML If I consider the atmosphere above the ice, the surface flux
+ CML which is relevant for the air temperature dT/dt Eq
+ CML accounts for sensible and radiation (with different treatment
+ CML according to wave-length) fluxes but not for "latent heat flux",
+ CML since it does not contribute to heating the air.
+ CML So this diagnostic is only good for heat budget calculations within
+ CML the ice-ocean system.
+            SIatmQnt(I,J,bi,bj) =
+      &            maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
+ #ifndef SEAICE_GROWTH_LEGACY
+      &            a_QSWbyATM_cover(I,J) +
+ #endif /* SEAICE_GROWTH_LEGACY */
+      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
+ cgf 2) SItflux (analogous to tflux; includes advection by water
+ C             exchanged between atmosphere and ocean+ice)
+ C solid water going to atm, in precip units
+            tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
+      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
+      &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
+ C liquid water going to atm, in precip units
+            tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
+      &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
+      &         * ( ONE - AREApreTH(I,J) )
+ #ifdef ALLOW_RUNOFF
+      &         - RUNOFF(I,J,bi,bj)
+ #endif /* ALLOW_RUNOFF */
+      &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
+      &         *convertHI2PRECIP )
+ C In real fresh water flux + nonlinFS, we factor in the advected specific
+ C energy (referenced to 0 for 0deC liquid water). In virtual salt flux or
+ C linFS, rain/evap get a special treatment (see external_forcing_surf.F).
+            tmpscal1= - tmpscal1*
+      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
+       IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
+      &     .AND.(nonlinFreeSurf.NE.0) ) THEN
+            tmpscal2= - tmpscal2*
+      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
+      &        .AND.(nonlinFreeSurf.NE.0) ) THEN
+            tmpscal2= - tmpscal2*
+      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+            tmpscal2= - tmpscal2*HeatCapacity_Cp*
+      &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
+            tmpscal2= ZERO
+       ENDIF
+            SItflux(I,J,bi,bj)=
+      &            SIatmQnt(I,J,bi,bj)-tmpscal1-tmpscal2
+           ENDDO
+          ENDDO
  C compute net fresh water flux leaving/entering
  C the ocean, accounting for fresh/salt water stocks.
  C ==================================================
-Line 2422 
 C ======================================
+Line 2516 
 C ======================================
       &             +d_HEFFbyOCNonICE(I,J)
       &             +d_HEFFbyATMonOCN(I,J)
       &             +d_HEFFbyNEG(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
       &             +d_HEFFbyRLX(I,J)
  #endif
       &             +d_HSNWbyNEG(I,J)*SNOW2ICE
-Line 2436 
 C     If r_FWbySublim>0, then it is evap
+Line 2530 
 C     If r_FWbySublim>0, then it is evap
  #endif /* ALLOW_RUNOFF */
       &         + tmpscal1*convertHI2PRECIP
       &         )*rhoConstFresh
+ c and the flux leaving/entering the ocean+ice
+            SIatmFW(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
+      &          EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
+      &          - PRECIP(I,J,bi,bj)
+ #ifdef ALLOW_RUNOFF
+      &          - RUNOFF(I,J,bi,bj)
+ #endif /* ALLOW_RUNOFF */
+      &           )*rhoConstFresh
+      &     + a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
           ENDDO
          ENDDO
  #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
  C--
-Line 2508 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 2612 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
           ENDDO
          ENDIF
+ #ifdef ALLOW_BALANCE_FLUXES
+ C Compute tile integrals of heat/fresh water fluxes to/from atm.
+ C ==============================================================
+        FWFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceEmPmR ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           FWFsiTile(bi,bj) =
+      &      FWFsiTile(bi,bj) + SIatmFW(i,j,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+ c to translate global mean FWF adjustements (see below) we may need :
+        FWF2HFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           FWF2HFsiTile(bi,bj) = FWF2HFsiTile(bi,bj) +
+      &      HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+        HFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceQnet ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           HFsiTile(bi,bj) =
+      &      HFsiTile(bi,bj) + SItflux(i,j,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+ #endif
  C ===================================================================
  C ======================PART 8: diagnostics==========================
  C ===================================================================
-Line 2558 
 C three that actually need intermediate
+Line 2698 
 C three that actually need intermediate
  #ifdef ALLOW_ATM_TEMP
           DO J=1,sNy
            DO I=1,sNx
- CML If I consider the atmosphere above the ice, the surface flux
- CML which is relevant for the air temperature dT/dt Eq
- CML accounts for sensible and radiation (with different treatment
- CML according to wave-length) fluxes but not for "latent heat flux",
- CML since it does not contribute to heating the air.
- CML So this diagnostic is only good for heat budget calculations within
- CML the ice-ocean system.
-            DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
- #ifndef SEAICE_GROWTH_LEGACY
-      &            a_QSWbyATM_cover(I,J) +
- #endif /* SEAICE_GROWTH_LEGACY */
-      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
- C
             DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
       &       a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
- C
-            DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(
-      &          PRECIP(I,J,bi,bj)
-      &          - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
- #ifdef ALLOW_RUNOFF
-      &          + RUNOFF(I,J,bi,bj)
- #endif /* ALLOW_RUNOFF */
-      &           )*rhoConstFresh
-      &     - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
            ENDDO
           ENDDO
-          CALL DIAGNOSTICS_FILL(DIAGarrayA,
-      &      'SIatmQnt',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayB,
       &      'SIfwSubl',0,1,3,bi,bj,myThid)
-          CALL DIAGNOSTICS_FILL(DIAGarrayC,
-      &      'SIatmFW ',0,1,3,bi,bj,myThid)
  C
           DO J=1,sNy
            DO I=1,sNx
-Line 2597 
 C the actual Freshwater flux of sublimat
+Line 2711 
 C the actual Freshwater flux of sublimat
             DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
       &       * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
       &       * SEAICE_rhoIce * recip_deltaTtherm
- c the residual Freshwater flux of sublimated ice
+ C the residual Freshwater flux of sublimated ice
             DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
       &       * r_FWbySublim(I,J)
       &       * SEAICE_rhoIce * recip_deltaTtherm
-Line 2614 
 C the latent heat flux
+Line 2728 
 C the latent heat flux
           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl    ',0,1,3,bi,bj,myThid)
-          DO J=1,sNy
-           DO I=1,sNx
- c compute ice/snow water going to atm, in precip units
-            tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
-      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
-      &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
- c compute ocean water going to atm, in precip units
-            tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
-      &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
-      &         * ( ONE - AREApreTH(I,J) )
- #ifdef ALLOW_RUNOFF
-      &         - RUNOFF(I,J,bi,bj)
- #endif /* ALLOW_RUNOFF */
-      &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
-      &         *convertHI2PRECIP )
- c factor in the advected specific energy (referenced to 0 for 0deC liquid water)
-            tmpscal1= - tmpscal1*
-      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
-            IF (temp_EvPrRn.NE.UNSET_RL) THEN
-            tmpscal2= - tmpscal2*
-      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
-            ELSE
-            tmpscal2= - tmpscal2*
-      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
-            ENDIF
- c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX
-            DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
- #ifndef SEAICE_GROWTH_LEGACY
-      &            a_QSWbyATM_cover(I,J) +
- #endif
-      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
-      &            -tmpscal1-tmpscal2
-           ENDDO
-          ENDDO
-          CALL DIAGNOSTICS_FILL(DIAGarrayA,
-      &      'SItflux ',0,1,3,bi,bj,myThid)
  #endif /* ALLOW_ATM_TEMP */
          ENDIF
-Line 2660 
 C close bi,bj loops
+Line 2737 
 C close bi,bj loops
         ENDDO
        ENDDO
+ C ===================================================================
+ C =========PART 9: HF/FWF global integrals and balancing=============
+ C ===================================================================
+ #ifdef ALLOW_BALANCE_FLUXES
+ c 1)  global sums
+ # ifdef ALLOW_AUTODIFF_TAMC
+ CADJ STORE FWFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ CADJ STORE HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ CADJ STORE FWF2HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ # endif /* ALLOW_AUTODIFF_TAMC */
+       FWFsiGlob=0. _d 0
+       IF ( balanceEmPmR )
+      &   CALL GLOBAL_SUM_TILE_RL( FWFsiTile, FWFsiGlob, myThid )
+       FWF2HFsiGlob=0. _d 0
+       IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
+          CALL GLOBAL_SUM_TILE_RL(FWF2HFsiTile, FWF2HFsiGlob, myThid)
+       ELSEIF ( balanceEmPmR ) THEN
+          FWF2HFsiGlob=HeatCapacity_Cp * temp_EvPrRn * globalArea
+       ENDIF
+       HFsiGlob=0. _d 0
+       IF ( balanceQnet )
+      &   CALL GLOBAL_SUM_TILE_RL( HFsiTile, HFsiGlob, myThid )
+ c 2) global means
+ c mean SIatmFW
+       tmpscal0=FWFsiGlob / globalArea
+ c corresponding mean advection by atm to ocean+ice water exchange
+ c        (if mean SIatmFW was removed uniformely from ocean)
+       tmpscal1=FWFsiGlob / globalArea * FWF2HFsiGlob / globalArea
+ c mean SItflux (before potential adjustement due to SIatmFW)
+       tmpscal2=HFsiGlob / globalArea
+ c mean SItflux (after potential adjustement due to SIatmFW)
+       IF ( balanceEmPmR ) tmpscal2=tmpscal2-tmpscal1
+ c 3) balancing adjustments
+       IF ( balanceEmPmR ) THEN
+       DO bj=myByLo(myThid),myByHi(myThid)
+        DO bi=myBxLo(myThid),myBxHi(myThid)
+         DO j=1-OLy,sNy+OLy
+          DO i=1-OLx,sNx+OLx
+             empmr(i,j,bi,bj) = empmr(i,j,bi,bj) - tmpscal0
+             SIatmFW(i,j,bi,bj) = SIatmFW(i,j,bi,bj) - tmpscal0
+ c           adjust SItflux consistently
+             IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
+      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
+             tmpscal1=
+      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
+             ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
+      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
+             tmpscal1=
+      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
+             ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+             tmpscal1=
+      &       HeatCapacity_Cp*(temp_EvPrRn - theta(I,J,kSurface,bi,bj))
+             ELSE
+             tmpscal1=ZERO
+             ENDIF
+             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal0*tmpscal1
+ c           no qnet or tflux adjustement is needed
+          ENDDO
+         ENDDO
+        ENDDO
+       ENDDO
+       IF ( balancePrintMean ) THEN
+        _BEGIN_MASTER( myThid )
+        WRITE(msgbuf,'(a,a,e24.17)') 'rm Global mean of ',
+      &      'SIatmFW = ', tmpscal0
+        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &      SQUEEZE_RIGHT , myThid)
+        _END_MASTER( myThid )
+       ENDIF
+       ENDIF
+       IF ( balanceQnet ) THEN
+       DO bj=myByLo(myThid),myByHi(myThid)
+        DO bi=myBxLo(myThid),myBxHi(myThid)
+         DO j=1-OLy,sNy+OLy
+          DO i=1-OLx,sNx+OLx
+             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal2
+             qnet(i,j,bi,bj) = qnet(i,j,bi,bj) - tmpscal2
+             SIatmQnt(i,j,bi,bj) = SIatmQnt(i,j,bi,bj) - tmpscal2
+          ENDDO
+         ENDDO
+        ENDDO
+       ENDDO
+       IF ( balancePrintMean ) THEN
+        _BEGIN_MASTER( myThid )
+        WRITE(msgbuf,'(a,a,e24.17)') 'rm Global mean of ',
+      &      'SItflux = ', tmpscal2
+        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &      SQUEEZE_RIGHT , myThid)
+        _END_MASTER( myThid )
+       ENDIF
+       ENDIF
+ #endif /* */
+ #ifdef ALLOW_DIAGNOSTICS
+ c these diags need to be done outside of the bi,bj loop so that
+ c we may do potential global mean adjustement to them consistently.
+          CALL DIAGNOSTICS_FILL(SItflux,
+      &      'SItflux ',0,1,0,1,1,myThid)
+          CALL DIAGNOSTICS_FILL(SIatmQnt,
+      &      'SIatmQnt',0,1,0,1,1,myThid)
+ c SIatmFW follows the same convention as empmr -- SIatmFW diag does not
+          tmpscal1= - 1. _d 0
+          CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
+      &      tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
+ #endif /* ALLOW_DIAGNOSTICS */
        RETURN
        END

 Legend:



Removed from v.1.9
 


changed lines


 
Added in v.1.10
 Legend:



Removed from v.1.9
 


changed lines


 
Added in v.1.10
-Removed from v.1.9
+Added in v.1.10

	ViewVC Help
Powered by ViewVC 1.1.22