C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.212 2017/04/04 23:31:27 jmc Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif
#ifdef ALLOW_SALT_PLUME
# include "SALT_PLUME_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: SEAICE_GROWTH
C !INTERFACE:
SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE seaice_growth
C | o Updata ice thickness and snow depth
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "FFIELDS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
#ifdef ALLOW_EXF
# include "EXF_PARAM.h"
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_SALT_PLUME
# include "SALT_PLUME.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C myTime :: Simulation time
C myIter :: Simulation timestep number
C myThid :: Thread no. that called this routine.
_RL myTime
INTEGER myIter, myThid
CEOP
#if (defined ALLOW_EXF) (defined ALLOW_ATM_TEMP)
C !FUNCTIONS:
#ifdef ALLOW_DIAGNOSTICS
LOGICAL DIAGNOSTICS_IS_ON
EXTERNAL
#endif
C !LOCAL VARIABLES:
C === Local variables ===
C
C unit/sign convention:
C Within the thermodynamic computation all stocks, except HSNOW,
C are in 'effective ice meters' units, and >0 implies more ice.
C This holds for stocks due to ocean and atmosphere heat,
C at the outset of 'PART 2: determine heat fluxes/stocks'
C and until 'PART 7: determine ocean model forcing'
C This strategy minimizes the need for multiplications/divisions
C by ice fraction, heat capacity, etc. The only conversions that
C occurs are for the HSNOW (in effective snow meters) and
C PRECIP (fresh water m/s).
C
C HEFF is effective Hice thickness (m3/m2)
C HSNOW is Heffective snow thickness (m3/m2)
C HSALT is Heffective salt content (g/m2)
C AREA is the seaice cover fraction (0<=AREA<=1)
C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
C
C For all other stocks/increments, such as d_HEFFbyATMonOCN
C or a_QbyATM_cover, the naming convention is as follows:
C The prefix 'a_' means available, the prefix 'd_' means delta
C (i.e. increment), and the prefix 'r_' means residual.
C The suffix '_cover' denotes a value for the ice covered fraction
C of the grid cell, whereas '_open' is for the open water fraction.
C The main part of the name states what ice/snow stock is concerned
C (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
C is the increment of HEFF due to the ATMosphere extracting heat from the
C OCeaN surface, or providing heat to the OCeaN surface).
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
C msgBuf :: Informational/error message buffer
#ifdef ALLOW_BALANCE_FLUXES
CHARACTER*(MAX_LEN_MBUF) msgBuf
#elif (defined (SEAICE_DEBUG))
CHARACTER*(MAX_LEN_MBUF) msgBuf
CHARACTER*12 msgBufForm
#endif
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
#ifndef SEAICE_ITD
C facilitate multi-category snow implementation
_RL pFac, pFacSnow
#endif
C additional factors accounting for a non-uniform sea-ice PDF
_RL denominator, recip_denominator, areaPDFfac
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
#ifdef SEAICE_GREASE
INTEGER iTrGrease
_RL greaseDecayTime
_RL greaseNewFrazil
_RL THIRD
PARAMETER (THIRD = 1.0 _d 0 / 3.0 _d 0)
#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)
#ifndef SEAICE_ITD
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)
#endif
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)
#endif
C wind speed
_RL UG (1:sNx,1:sNy)
C temporary variables available for the various computations
_RL tmparr1 (1:sNx,1:sNy)
_RL ticeInMult (1:sNx,1:sNy,nITD)
_RL ticeOutMult (1:sNx,1:sNy,nITD)
_RL heffActualMult (1:sNx,1:sNy,nITD)
_RL hsnowActualMult (1:sNx,1:sNy,nITD)
#ifdef SEAICE_ITD
_RL recip_heffActualMult(1:sNx,1:sNy,nITD)
#endif
_RL a_QbyATMmult_cover (1:sNx,1:sNy,nITD)
_RL a_QSWbyATMmult_cover(1:sNx,1:sNy,nITD)
_RL a_FWbySublimMult (1:sNx,1:sNy,nITD)
#ifdef SEAICE_GREASE
_RL greaseLayerThick (1:sNx,1:sNy)
_RL d_HEFFbyGREASE (1:sNx,1:sNy)
#endif
#ifdef SEAICE_ITD
_RL r_QbyATMmult_cover (1:sNx,1:sNy,nITD)
_RL r_FWbySublimMult (1:sNx,1:sNy,nITD)
C for lateral melt parameterization:
_RL latMeltFrac (1:sNx,1:sNy,nITD)
_RL latMeltRate (1:sNx,1:sNy,nITD)
_RL floeAlpha
_RL floeDiameter
_RL floeDiameterMin
_RL floeDiameterMax
#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
C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
C r_QbyATM_open :: same but for open water
_RL a_QbyATM_cover (1:sNx,1:sNy)
_RL a_QbyATM_open (1:sNx,1:sNy)
_RL r_QbyATM_cover (1:sNx,1:sNy)
_RL r_QbyATM_open (1:sNx,1:sNy)
C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
_RL a_QSWbyATM_open (1:sNx,1:sNy)
_RL a_QSWbyATM_cover (1:sNx,1:sNy)
C a_QbyOCN :: available heat (in W/m^2) due to the
C interaction of the ice pack and the ocean surface
C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
C processes have been accounted for
_RL a_QbyOCN (1:sNx,1:sNy)
_RL r_QbyOCN (1:sNx,1:sNy)
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
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
_RL d_HEFFbyFLOODING (1:sNx,1:sNy)
C The mean ice thickness increments due to atmospheric fluxes over the open
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)
_RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
_RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
_RL d_HSNWbyRAIN (1:sNx,1:sNy)
_RL d_HFRWbyRAIN (1:sNx,1:sNy)
C a_FWbySublim :: fresh water flux implied by latent heat of
C sublimation to atmosphere, same sign convention
C as EVAP (positive upward)
_RL a_FWbySublim (1:sNx,1:sNy)
_RL r_FWbySublim (1:sNx,1:sNy)
_RL d_HEFFbySublim (1:sNx,1:sNy)
_RL d_HSNWbySublim (1:sNx,1:sNy)
#ifdef SEAICE_CAP_SUBLIM
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,nITD)
#endif
#ifdef SEAICE_ITD
_RL d_HEFFbySublim_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HSNWbySublim_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HEFFbyOCNonICE_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HSNWbyATMonSNW_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HEFFbyATMonOCN_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HEFFbyATMonOCN_cover_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HEFFbyATMonOCN_open_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HSNWbyRAIN_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HSNWbyOCNonSNW_ITD (1:sNx,1:sNy,1:nITD)
_RL d_HEFFbyFLOODING_ITD (1:sNx,1:sNy,1:nITD)
#endif
#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)
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 /* 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
#endif
#ifdef ALLOW_SALT_PLUME
_RL localSPfrac (1:sNx,1:sNy)
#ifdef SALT_PLUME_IN_LEADS
_RL leadPlumeFraction (1:sNx,1:sNy)
_RL IceGrowthRateInLeads (1:sNx,1:sNy)
#endif /* SALT_PLUME_IN_LEADS */
#endif /* ALLOW_SALT_PLUME */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C ===================================================================
C =================PART 0: constants and initializations=============
C ===================================================================
IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
C avoid unnecessary divisions in loops
recip_multDim = SEAICE_multDim
recip_multDim = ONE / recip_multDim
C above/below: double/single precision calculation of recip_multDim
c recip_multDim = 1./float(SEAICE_multDim)
recip_deltaTtherm = ONE / SEAICE_deltaTtherm
recip_rhoIce = ONE / SEAICE_rhoIce
C Cutoff for iceload
heffTooHeavy=drF(kSurface) / 5. _d 0
C RATIO OF SEA ICE DENSITY to SNOW DENSITY
ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
SNOW2ICE = ONE / ICE2SNOW
C HEAT OF FUSION OF ICE (J/m^3)
QI = SEAICE_rhoIce*SEAICE_lhFusion
recip_QI = ONE / QI
C HEAT OF FUSION OF SNOW (J/m^3)
QS = SEAICE_rhoSnow*SEAICE_lhFusion
C ICE LATENT HEAT CONSTANT
lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
C regularization constants
area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
C conversion factors to go from Q (W/m2) to HEFF (ice meters)
convertQ2HI=SEAICE_deltaTtherm/QI
convertHI2Q = ONE/convertQ2HI
C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
convertHI2PRECIP = ONE/convertPRECIP2HI
C compute parameters for thickness pdf (cheap enough to do it here
C and not in seaice_readparms and store in common block)
denominator = 0. _d 0
DO IT=1,SEAICE_multDim
denominator = denominator + IT * SEAICE_pdf(IT)
ENDDO
denominator = (2.0 _d 0 * denominator) - 1.0 _d 0
recip_denominator = 1. _d 0 / denominator
#ifdef SEAICE_ITD
areaPDFfac = 1. _d 0
#else
areaPDFfac = denominator * recip_multDim
#endif /* SEAICE_ITD */
#ifdef SEAICE_ITD
C constants for lateral melt parameterization:
C following Steele (1992), Equ. 2
floeAlpha = 0.66 _d 0
C typical mean diameter used in CICE 4.1:
C (this is currently computed as a function of ice concentration
C following a suggestion by Luepkes at al. (2012))
C floeDiameter = 300. _d 0
C parameters needed for variable floe diameter following Luepkes et al. (2012):
floeDiameterMin = 8. _d 0
floeDiameterMax = 300. _d 0
#endif
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
iicekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef SEAICE_GREASE
C time scale of grease ice decline by solidification,
C with 50% grease ice becoming solid pancake ice within 1 day:
greaseDecayTime=1.44*86400. _d 0
C store position of 'grease' in tracer array
iTrGrease=-1
DO iTr = 1, SItrNumInUse
if (SItrName(iTr).EQ.'grease') iTrGrease=iTr
ENDDO
#endif
C array initializations
C =====================
DO J=1,sNy
DO I=1,sNx
a_QbyATM_cover (I,J) = 0.0 _d 0
a_QbyATM_open(I,J) = 0.0 _d 0
r_QbyATM_cover (I,J) = 0.0 _d 0
r_QbyATM_open (I,J) = 0.0 _d 0
a_QSWbyATM_open (I,J) = 0.0 _d 0
a_QSWbyATM_cover (I,J) = 0.0 _d 0
a_QbyOCN (I,J) = 0.0 _d 0
r_QbyOCN (I,J) = 0.0 _d 0
#ifdef ALLOW_DIAGNOSTICS
d_AREAbyATM(I,J) = 0.0 _d 0
d_AREAbyICE(I,J) = 0.0 _d 0
d_AREAbyOCN(I,J) = 0.0 _d 0
#endif
d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
d_HEFFbyFLOODING(I,J) = 0.0 _d 0
d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
d_HEFFbyATMonOCN_cover(I,J) = 0.0 _d 0
d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
d_HSNWbyRAIN(I,J) = 0.0 _d 0
a_FWbySublim(I,J) = 0.0 _d 0
r_FWbySublim(I,J) = 0.0 _d 0
d_HEFFbySublim(I,J) = 0.0 _d 0
d_HSNWbySublim(I,J) = 0.0 _d 0
#ifdef SEAICE_CAP_SUBLIM
latentHeatFluxMax(I,J) = 0.0 _d 0
#endif
d_HFRWbyRAIN(I,J) = 0.0 _d 0
tmparr1(I,J) = 0.0 _d 0
#ifdef SEAICE_GREASE
greaseLayerThick(I,J) = 0.0 _d 0
d_HEFFbyGREASE(I,J) = 0.0 _d 0
#endif
DO IT=1,SEAICE_multDim
ticeInMult(I,J,IT) = 0.0 _d 0
ticeOutMult(I,J,IT) = 0.0 _d 0
a_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
a_QSWbyATMmult_cover(I,J,IT) = 0.0 _d 0
a_FWbySublimMult(I,J,IT) = 0.0 _d 0
#ifdef SEAICE_CAP_SUBLIM
latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
#endif
#ifdef SEAICE_ITD
d_HEFFbySublim_ITD(I,J,IT) = 0.0 _d 0
d_HSNWbySublim_ITD(I,J,IT) = 0.0 _d 0
d_HEFFbyOCNonICE_ITD(I,J,IT) = 0.0 _d 0
d_HSNWbyATMonSNW_ITD(I,J,IT) = 0.0 _d 0
d_HEFFbyATMonOCN_ITD(I,J,IT) = 0.0 _d 0
d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = 0.0 _d 0
d_HEFFbyATMonOCN_open_ITD(I,J,IT) = 0.0 _d 0
d_HSNWbyRAIN_ITD(I,J,IT) = 0.0 _d 0
d_HSNWbyOCNonSNW_ITD(I,J,IT) = 0.0 _d 0
d_HEFFbyFLOODING_ITD(I,J,IT) = 0.0 _d 0
r_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
r_FWbySublimMult(I,J,IT) = 0.0 _d 0
C for lateral melt parameterization:
latMeltFrac(I,J,IT) = 0.0 _d 0
latMeltRate(I,J,IT) = 0.0 _d 0
#endif
ENDDO
ENDDO
ENDDO
#if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
DO J=1-oLy,sNy+oLy
DO I=1-oLx,sNx+oLx
frWtrAtm(I,J,bi,bj) = 0.0 _d 0
ENDDO
ENDDO
#endif
C =====================================================================
C ===========PART 1: treat pathological cases (post advdiff)===========
C =====================================================================
C This part has been mostly moved to S/R seaice_reg_ridge, which is
C called before S/R seaice_growth
C store regularized values of heff, hsnow, area at the onset of thermo.
DO J=1,sNy
DO I=1,sNx
HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
AREApreTH(I,J)=AREA(I,J,bi,bj)
#ifdef ALLOW_DIAGNOSTICS
DIAGarrayB(I,J) = AREA(I,J,bi,bj)
DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
#endif
#ifdef ALLOW_SITRACER
SItrHEFF(I,J,bi,bj,1)=HEFF(I,J,bi,bj)
SItrAREA(I,J,bi,bj,2)=AREA(I,J,bi,bj)
#endif
ENDDO
ENDDO
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
C keep track of areal and volume fraction of each ITD category
IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
ELSE
C if there is no ice, potential growth starts in 1st category
IF (IT .EQ. 1) THEN
areaFracFactor(I,J,IT)=ONE
ELSE
areaFracFactor(I,J,IT)=ZERO
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
#ifdef ALLOW_SITRACER
C prepare SItrHEFF to be computed as cumulative sum
DO iTr=2,5
DO J=1,sNy
DO I=1,sNx
SItrHEFF(I,J,bi,bj,iTr)=ZERO
ENDDO
ENDDO
ENDDO
C prepare SItrAREA to be computed as cumulative sum
DO J=1,sNy
DO I=1,sNx
SItrAREA(I,J,bi,bj,3)=ZERO
ENDDO
ENDDO
#endif
#endif /* SEAICE_ITD */
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
#ifdef ALLOW_SITRACER
DO iTr = 1, SItrNumInUse
WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT'
IF (SItrMate(iTr).EQ.'HEFF') THEN
CALL DIAGNOSTICS_FRACT_FILL(
I SItracer(1-OLx,1-OLy,bi,bj,iTr),HEFF(1-OLx,1-OLy,bi,bj),
I ONE, 1, diagName,0,1,2,bi,bj,myThid )
ELSE
CALL DIAGNOSTICS_FRACT_FILL(
I SItracer(1-OLx,1-OLy,bi,bj,iTr),AREA(1-OLx,1-OLy,bi,bj),
I ONE, 1, diagName,0,1,2,bi,bj,myThid )
ENDIF
ENDDO
#endif /* ALLOW_SITRACER */
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
Cgf no additional dependency of air-sea fluxes to ice
IF ( SEAICEadjMODE.GE.1 ) THEN
DO J=1,sNy
DO I=1,sNx
HEFFpreTH(I,J) = 0. _d 0
HSNWpreTH(I,J) = 0. _d 0
AREApreTH(I,J) = 0. _d 0
ENDDO
ENDDO
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
HEFFITDpreTH(I,J,IT) = 0. _d 0
HSNWITDpreTH(I,J,IT) = 0. _d 0
AREAITDpreTH(I,J,IT) = 0. _d 0
ENDDO
ENDDO
ENDDO
#endif
ENDIF
#endif
#if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
DO J=1,sNy
DO I=1,sNx
AREAforAtmFW(I,J,bi,bj) = AREApreTH(I,J)
ENDDO
ENDDO
#endif
C COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
C regularize AREA with SEAICE_area_reg
tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
& + area_reg_sq)
C heffActual calculated with the regularized AREA
tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
C regularize heffActual with SEAICE_hice_reg (add lower bound)
heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
& + hice_reg_sq)
C hsnowActual calculated with the regularized AREA
hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
C regularize the inverse of heffActual by hice_reg
recip_heffActualMult(I,J,IT) = AREAITDpreTH(I,J,IT) /
& sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
& + hice_reg_sq)
C Do not regularize when HEFFpreTH = 0
ELSE
heffActualMult(I,J,IT) = ZERO
hsnowActualMult(I,J,IT) = ZERO
recip_heffActualMult(I,J,IT) = ZERO
ENDIF
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
IF (HEFFpreTH(I,J) .GT. ZERO) THEN
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
tmpscal2 = HEFFpreTH(I,J) / tmpscal1
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
hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
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
ELSE
heffActual(I,J) = ZERO
hsnowActual(I,J) = ZERO
recip_heffActual(I,J) = ZERO
ENDIF
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
#endif
#ifdef SEAICE_CAP_SUBLIM
C COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
C AND SNOW THICKNESS
C The latent heat flux over the sea ice which
C will sublimate all of the snow and ice over one time
C step (W/m^2)
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
& (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
& HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
& /AREAITDpreTH(I,J,IT)
ELSE
latentHeatFluxMaxMult(I,J,IT) = ZERO
ENDIF
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
IF (HEFFpreTH(I,J) .GT. ZERO) THEN
latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
& (HEFFpreTH(I,J) * SEAICE_rhoIce +
& HSNWpreTH(I,J) * SEAICE_rhoSnow)/AREApreTH(I,J)
ELSE
latentHeatFluxMax(I,J) = ZERO
ENDIF
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#endif /* SEAICE_CAP_SUBLIM */
C ===================================================================
C ================PART 2: determine heat fluxes/stocks===============
C ===================================================================
C determine available heat due to the atmosphere -- for open water
C ================================================================
DO j=1,sNy
DO i=1,sNx
C ocean surface/mixed layer temperature
TmixLoc(i,j) = theta(i,j,kSurface,bi,bj)+celsius2K
C wind speed from exf
UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
CALL SEAICE_BUDGET_OCEAN(
I UG,
I TmixLoc,
O a_QbyATM_open, a_QSWbyATM_open,
I bi, bj, myTime, myIter, myThid )
C determine available heat due to the atmosphere -- for ice covered water
C =======================================================================
IF (useRelativeWind.AND.useAtmWind) THEN
C Compute relative wind speed over sea ice.
DO J=1,sNy
DO I=1,sNx
SPEED_SQ =
& (uWind(I,J,bi,bj)
& +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
& +uVel(i+1,j,kSurface,bi,bj))
& -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
& +(vWind(I,J,bi,bj)
& +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
& +vVel(i,j+1,kSurface,bi,bj))
& -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
UG(I,J)=SEAICE_EPS
ELSE
UG(I,J)=SQRT(SPEED_SQ)
ENDIF
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE heffActual = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE tices(:,:,:,bi,bj)
CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C-- Start loop over multi-categories
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
ticeInMult(I,J,IT) = TICES(I,J,IT,bi,bj)
ticeOutMult(I,J,IT) = TICES(I,J,IT,bi,bj)
TICES(I,J,IT,bi,bj) = ZERO
ENDDO
ENDDO
#ifndef SEAICE_ITD
C for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult have been
C calculated above (instead of heffActual and latentHeatFluxMax)
C-- assume homogeneous distribution between 0 and 2 x heffActual
pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_denominator
pFacSnow = 1. _d 0
IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
DO J=1,sNy
DO I=1,sNx
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
ENDDO
ENDDO
#endif /* ndef SEAICE_ITD */
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 ticeInMult = comlev1_bibj, key = iicekey, byte = isbyte
# ifdef SEAICE_CAP_SUBLIM
CADJ STORE latentHeatFluxMaxMult
CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
# endif
CADJ STORE a_QbyATMmult_cover =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATMmult_cover =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_FWbySublimMult =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO IT=1,SEAICE_multDim
CALL SEAICE_SOLVE4TEMP(
I UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
#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),
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
CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
# endif
CADJ STORE a_QbyATMmult_cover =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATMmult_cover =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_FWbySublimMult =
CADJ & comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
C update TICES
CMLC and tIce, if required later on (currently not the case)
CML#ifdef SEAICE_ITD
CMLC calculate area weighted mean
CMLC (although the ice temperature relates to its energy content
CMLC and hence should be averaged weighted by ice volume,
CMLC the temperature here is a result of the fluxes through the ice surface
CMLC computed individually for each single category in SEAICE_SOLVE4TEMP
CMLC and hence is averaged area weighted [areaFracFactor])
CML tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
CML & + ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
CML#else
CML tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
CML & + ticeOutMult(I,J,IT)*SEAICE_PDF(IT)
CML#endif
TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
C average over categories
#ifdef SEAICE_ITD
C calculate area weighted mean
C (fluxes are per unit (ice surface) area and are thus area weighted)
a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
& + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
& + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
a_FWbySublim (I,J) = a_FWbySublim(I,J)
& + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
#else
a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
& + a_QbyATMmult_cover(I,J,IT)*SEAICE_PDF(IT)
a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
& + a_QSWbyATMmult_cover(I,J,IT)*SEAICE_PDF(IT)
a_FWbySublim (I,J) = a_FWbySublim(I,J)
& + a_FWbySublimMult(I,J,IT)*SEAICE_PDF(IT)
#endif
ENDDO
ENDDO
ENDDO
#ifdef SEAICE_CAP_SUBLIM
# ifdef ALLOW_DIAGNOSTICS
DO J=1,sNy
DO I=1,sNx
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
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(DIAGarrayA,
& 'SIactLHF',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(latentHeatFluxMax,
& 'SImaxLHF',0,1,3,bi,bj,myThid)
ENDIF
# endif /* ALLOW_DIAGNOSTICS */
#endif /* SEAICE_CAP_SUBLIM */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C switch heat fluxes from W/m2 to 'effective' ice meters
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
& * convertQ2HI * AREAITDpreTH(I,J,IT)
a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
& * convertQ2HI * AREAITDpreTH(I,J,IT)
C and initialize r_QbyATMmult_cover
r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
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
a_FWbySublimMult(I,J,IT) = ZERO
#endif
a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
& * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
ENDDO
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
& * convertQ2HI * ( ONE - AREApreTH(I,J) )
a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
& * convertQ2HI * ( ONE - AREApreTH(I,J) )
C and initialize r_QbyATM_open
r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
ENDDO
ENDDO
#else /* SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
& * convertQ2HI * AREApreTH(I,J)
a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
& * convertQ2HI * AREApreTH(I,J)
a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
& * convertQ2HI * ( ONE - AREApreTH(I,J) )
a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
& * convertQ2HI * ( ONE - AREApreTH(I,J) )
C and initialize r_QbyATM_cover/r_QbyATM_open
r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
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
a_FWbySublim(I,J) = ZERO
#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)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE r_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE r_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE r_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
Cgf no additional dependency through ice cover
IF ( SEAICEadjMODE.GE.3 ) THEN
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
a_QbyATMmult_cover(I,J,IT) = 0. _d 0
r_QbyATMmult_cover(I,J,IT) = 0. _d 0
a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
a_QbyATM_cover(I,J) = 0. _d 0
r_QbyATM_cover(I,J) = 0. _d 0
a_QSWbyATM_cover(I,J) = 0. _d 0
ENDDO
ENDDO
#endif /* SEAICE_ITD */
ENDIF
#endif
C determine available heat due to the ice pack tying the
C underlying surface water temperature to freezing point
C ======================================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
#endif
DO J=1,sNy
DO I=1,sNx
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
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)
IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
& (.NOT.SEAICE_mcPheeStepFunc) ) THEN
MixedLayerTurbulenceFactor = ONE -
& SEAICE_mcPheeTaper * AREApreTH(I,J)
ELSEIF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
& (SEAICE_mcPheeStepFunc) ) THEN
MixedLayerTurbulenceFactor = ONE - SEAICE_mcPheeTaper
ELSE
MixedLayerTurbulenceFactor = ONE
ENDIF
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
a_QbyOCN(i,j) =
& tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
r_QbyOCN(i,j) = a_QbyOCN(i,j)
ENDDO
ENDDO
#ifdef SEAICE_ITD
C determine lateral melt rate at floe edges based on an
C average floe diameter or a floe size distribution
C following Steele (1992, Tab. 2)
C ======================================================
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
tempFrz = SEAICE_tempFrz0 +
& SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
tmpscal1=(theta(I,J,kSurface,bi,bj)-tempFrz)
tmpscal2=sqrt(0.87 + 0.067*UG(i,j)) * UG(i,j)
C variable floe diameter following Luepkes et al. (2012, JGR, Equ. 26)
C with beta=1
CML tmpscal3=ONE/(ONE-(floeDiameterMin/floeDiameterMax))
CML floeDiameter = floeDiameterMin
CML & * (tmpscal3 / (tmpscal3-AREApreTH(I,J)))
C this form involves fewer divisions but gives the same result
floeDiameter = floeDiameterMin * floeDiameterMax
& / ( floeDiameterMax*( 1. _d 0 - AREApreTH(I,J) )
& + floeDiameterMin*AREApreTH(I,J) )
C following the idea of SEAICE_areaLossFormula == 1:
IF (a_QbyATMmult_cover(i,j,it).LT.ZERO .OR.
& a_QbyATM_open(i,j) .LT.ZERO .OR.
& a_QbyOCN(i,j) .LT.ZERO) THEN
C lateral melt rate as suggested by Perovich, 1983 (PhD thesis)
C latMeltRate(i,j,it) = 1.6 _d -6 * tmpscal1**1.36
C The following for does the same, but is faster
latMeltRate(i,j,it) = ZERO
IF (tmpscal1 .GT. ZERO)
& latMeltRate(i,j,it) = 1.6 _d -6 * exp(1.36*log(tmpscal1))
C lateral melt rate as suggested by Maykut and Perovich, 1987
C (JGR 92(C7)), Equ. 24
c latMeltRate(i,j,it) = 13.5 _d -6 * tmpscal2 * tmpscal1**1.3
C further suggestion by Maykut and Perovich to avoid
C latMeltRate -> 0 for UG -> 0
c latMeltRate(i,j,it) = (1.6 _d -6 + 13.5 _d -6 * tmpscal2)
c & * tmpscal1**1.3
C factor determining fraction of area and ice volume reduction
C due to lateral melt
latMeltFrac(i,j,it) =
& latMeltRate(i,j,it)*SEAICE_deltaTtherm*PI /
& (floeAlpha * floeDiameter)
latMeltFrac(i,j,it)=max(ZERO, min(latMeltFrac(i,j,it),ONE))
ELSE
latMeltRate(i,j,it)=0.0 _d 0
latMeltFrac(i,j,it)=0.0 _d 0
ENDIF
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
#endif
C ===================================================================
C =========PART 3: determine effective thicknesses increments========
C ===================================================================
#ifdef SEAICE_GREASE
C convert SItracer 'grease' from ratio to grease ice volume:
C ==========================================================
DO J=1,sNy
DO I=1,sNx
SItracer(I,J,bi,bj,iTrGrease) =
& SItracer(I,J,bi,bj,iTrGrease) * HEFF(I,J,bi,bj)
ENDDO
ENDDO
C compute actual grease ice layer thickness
C as a function of wind stress
C following Smedsrud [2011, Ann.Glac.]
C =========================================
DO J=1,sNy
DO I=1,sNx
C
C computing compaction force acting on grease
C (= air and water stress acting on ice)
C wind stress (using calculation of SPEED_SQ & UG as template)
C u_a**2 + v_a**2:
tmpscal1 = uWind(I,J,bi,bj)*uWind(I,J,bi,bj)
& + vWind(I,J,bi,bj)*vWind(I,J,bi,bj)
IF ( tmpscal1 .LE. SEAICE_EPS_SQ ) THEN
tmpscal1=SEAICE_EPS
ELSE
tmpscal1=SQRT(tmpscal2)
ENDIF
tmpscal1 = 1.4 _d 0 * 1.3 _d -3 * tmpscal1
C water stress
C u_w - u_i:
tmpscal2 =
& 0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
& +uVel(i+1,j,kSurface,bi,bj))
& -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj))
C v_w - v_i:
tmpscal3 =
& 0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
& +vVel(i,j+1,kSurface,bi,bj))
& -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj))
C (u_w - u_i)**2 + (v_w - v_i)**2:
tmpscal4 = (tmpscal2*tmpscal2 + tmpscal3*tmpscal3)
IF ( tmpscal4 .LE. SEAICE_EPS_SQ ) THEN
tmpscal4=SEAICE_EPS
ELSE
tmpscal4=SQRT(tmpscal4)
ENDIF
tmpscal4 = 1027.0 _d 0 * 6.0 _d -3 * tmpscal4
C magnitude of compined stresses:
tmpscal0 =
& ( tmpscal1 * uWind(I,J,bi,bj) + tmpscal4 * tmpscal2 )**2
& + ( tmpscal1 * vWind(I,J,bi,bj) + tmpscal4 * tmpscal3 )**2
IF ( tmpscal0 .LE. SEAICE_EPS_SQ ) THEN
tmpscal0=SEAICE_EPS
ELSE
tmpscal0=SQRT(tmpscal0)
ENDIF
C
C mean grid cell width between tracer points
tmpscal3 = 0.5 _d 0 * (dxC(I,J,bi,bj)+dyC(I,J,bi,bj))
C grease ice volume Vg [m^3/m] as in Smedsrud [2011]
C is SItracer * gridcell_area / lead_width
C where lead width is lead extent along lead,
C i.e. perpendicular to L_lead in Smedsrud (2011),
C which is in the absence of lead length statistics
C simply the grid cell length
tmpscal4 = 4. _d 0 * SItracer(I,J,bi,bj,iTrGrease) * tmpscal3
C
C mean grease ice layer thickness , Eq.10 in Smedsrud [2011] but incl. water stress
greaseLayerThick(I,J) = 0.763 _d 0
C grease ice volume
& * ( tmpscal4
C times magnitude of vector sum of air and water stresses
C (in fact, only the component perpendicular to the lead edge, i.e. along L_lead counts)
C ATTENTION: since we do not have lead orientation with respect to wind
C we use 80% of the total force assuming angles 45-90deg between wind and lead edge
& * 0.8 _d 0 * tmpscal0
C devided by K_r = 100 N/m^3 (resistance of grease against compression)
& * 0.01 _d 0 )**THIRD
C
C assure a minimum thickness of 4 cm (equals HO=0.01):
greaseLayerThick(I,J)=max(4. _d -2, greaseLayerThick(I,J))
C ... and a maximum thickness of 4 m (equals HO=1.0):
greaseLayerThick(I,J)=min(4. _d 0 , greaseLayerThick(I,J))
C
ENDDO
ENDDO
#endif /* SEAICE_GREASE */
C compute snow/ice tendency due to sublimation
C ============================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
#endif /* SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
C First sublimate/deposite snow
tmpscal2 =
#ifdef SEAICE_ITD
& MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
& *SNOW2ICE),ZERO)
d_HSNWbySublim_ITD(I,J,IT) = - tmpscal2 * ICE2SNOW
C accumulate change over ITD categories
d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2
& *ICE2SNOW
r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
#else /* ndef SEAICE_ITD */
& MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) - tmpscal2*ICE2SNOW
r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
#endif /* SEAICE_ITD */
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
C If anything is left, sublimate ice
tmpscal2 =
#ifdef SEAICE_ITD
& MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
d_HEFFbySublim_ITD(I,J,IT) = - tmpscal2
C accumulate change over ITD categories
d_HEFFbySublim(I,J) = d_HEFFbySublim(I,J) - tmpscal2
r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
#else /* ndef SEAICE_ITD */
& MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
d_HEFFbySublim(I,J) = - tmpscal2
HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal2
r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
#endif /* SEAICE_ITD */
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
C If anything is left, it will be evaporated from the ocean rather than sublimated.
C Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
C remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
#ifdef SEAICE_ITD
a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
& - r_FWbySublimMult(I,J,IT)
r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
& - r_FWbySublimMult(I,J,IT)
#else /* ndef SEAICE_ITD */
a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
#endif /* SEAICE_ITD */
ENDDO
ENDDO
#ifdef SEAICE_ITD
C end IT loop
ENDDO
#endif /* SEAICE_ITD */
C compute ice thickness tendency due to ice-ocean interaction
C ===========================================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF (.NOT.SEAICE_growMeltByConv) THEN
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
C ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
C equally distributed under the ice and hence weighted by
C fractional area of each thickness category
tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
& -HEFFITD(I,J,IT,bi,bj))
d_HEFFbyOCNonICE_ITD(I,J,IT)=tmpscal1
d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
ENDDO
ENDDO
ENDDO
#ifdef ALLOW_SITRACER
DO J=1,sNy
DO I=1,sNx
SItrHEFF(I,J,bi,bj,2) = HEFFpreTH(I,J)
& + d_HEFFbySublim(I,J)
& + d_HEFFbyOCNonICE(I,J)
ENDDO
ENDDO
#endif
DO J=1,sNy
DO I=1,sNx
r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
ENDDO
ENDDO
#else /* SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
#ifdef ALLOW_SITRACER
SItrHEFF(I,J,bi,bj,2)=HEFF(I,J,bi,bj)
#endif
ENDDO
ENDDO
#endif /* SEAICE_ITD */
ENDIF !SEAICE_growMeltByConv
C compute snow melt tendency due to snow-atmosphere interaction
C ==================================================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
C Convert to standard units (meters of ice) rather than to meters
C of snow. This appears to be more robust.
tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
& -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
tmpscal2=MIN(tmpscal1,0. _d 0)
#ifdef SEAICE_MODIFY_GROWTH_ADJ
Cgf no additional dependency through snow
IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
#endif
d_HSNWbyATMonSNW_ITD(I,J,IT) = tmpscal2*ICE2SNOW
d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
& + tmpscal2*ICE2SNOW
r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
& - tmpscal2
ENDDO
ENDDO
ENDDO
#else /* SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
C Convert to standard units (meters of ice) rather than to meters
C of snow. This appears to be more robust.
tmpscal1=MAX(r_QbyATM_cover(I,J),-HSNOW(I,J,bi,bj)*SNOW2ICE)
tmpscal2=MIN(tmpscal1,0. _d 0)
#ifdef SEAICE_MODIFY_GROWTH_ADJ
Cgf no additional dependency through snow
IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
#endif
d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
ENDDO
ENDDO
#endif /* SEAICE_ITD */
C compute ice thickness tendency due to the atmosphere
C ====================================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
Cgf note: this block is not actually tested by lab_sea
Cgf where all experiments start in January. So even though
Cgf the v1.81=>v1.82 revision would change results in
Cgf warming conditions, the lab_sea results were not changed.
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
tmpscal1 = HEFFITDpreTH(I,J,IT)
& + d_HEFFbySublim_ITD(I,J,IT)
& + d_HEFFbyOCNonICE_ITD(I,J,IT)
tmpscal2 = MAX(-tmpscal1,
& r_QbyATMmult_cover(I,J,IT)
C Limit ice growth by potential melt by ocean
& + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2
d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
& + tmpscal2
d_HEFFbyATMonOCN_ITD(I,J,IT) = d_HEFFbyATMonOCN_ITD(I,J,IT)
& + tmpscal2
d_HEFFbyATMonOCN(I,J) = d_HEFFbyATMonOCN(I,J)
& + tmpscal2
r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
& - tmpscal2
ENDDO
ENDDO
ENDDO
#ifdef ALLOW_SITRACER
DO J=1,sNy
DO I=1,sNx
SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,2)
& + d_HEFFbyATMonOCN_cover(I,J)
ENDDO
ENDDO
#endif
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
C Limit ice growth by potential melt by ocean
& AREApreTH(I,J) * r_QbyOCN(I,J))
d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
#ifdef ALLOW_SITRACER
SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
#endif
ENDDO
ENDDO
#endif /* SEAICE_ITD */
C add snow precipitation to HSNOW.
C =================================================
#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 */
IF ( snowPrecipFile .NE. ' ' ) THEN
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 ! would require computing tIce
c IF (atemp(I,J,bi,bj) .LT. celsius2K) 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
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
ENDDO
ENDDO
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
d_HSNWbyRAIN_ITD(I,J,IT)
& = d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
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
C compute snow melt due to heat available from ocean.
C =================================================================
Cgf do we need to keep this comment and cpp bracket?
Cph( very sensitive bit here by JZ
#ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF (.NOT.SEAICE_growMeltByConv) THEN
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
tmpscal4 = HSNWITDpreTH(I,J,IT)
& + d_HSNWbySublim_ITD(I,J,IT)
& + d_HSNWbyATMonSNW_ITD(I,J,IT)
& + d_HSNWbyRAIN_ITD(I,J,IT)
tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
& -tmpscal4)
tmpscal2=MIN(tmpscal1,0. _d 0)
#ifdef SEAICE_MODIFY_GROWTH_ADJ
Cgf no additional dependency through snow
if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
#endif
d_HSNWbyOCNonSNW_ITD(I,J,IT) = tmpscal2
d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
tmpscal2=MIN(tmpscal1,0. _d 0)
#ifdef SEAICE_MODIFY_GROWTH_ADJ
Cgf no additional dependency through snow
if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
#endif
d_HSNWbyOCNonSNW(I,J) = tmpscal2
r_QbyOCN(I,J)=r_QbyOCN(I,J)
& -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
ENDIF !SEAICE_growMeltByConv
#endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
Cph)
C gain of new ice over open water
C ===============================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
#ifdef SEAICE_ITD
C HEFF will be updated at the end of PART 3,
C hence sum of tendencies so far is needed
tmpscal4 = HEFFpreTH(I,J)
& + d_HEFFbySublim(I,J)
& + d_HEFFbyOCNonICE(I,J)
& + d_HEFFbyATMonOCN(I,J)
#else /* ndef SEAICE_ITD */
C HEFF is updated step by step throughout seaice_growth
tmpscal4 = HEFF(I,J,bi,bj)
#endif /* SEAICE_ITD */
C Initial ice growth is triggered by open water
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 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)
tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
& -tmpscal4*facOpenMelt)*HEFFM(I,J,bi,bj)
#ifdef SEAICE_GREASE
C Grease ice is a tracer or "bucket" for newly formed frazil ice
C that instead of becoming solid sea ice instantly has a half-time
C of 1 day (see greaseDecayTime) before solidifying.
C The most important effect is that for fluxes the grease ice area/volume
C acts like open water.
C
C store freezing/melting condition:
greaseNewFrazil = max(0.0 _d 0, tmpscal3)
C
C 1) mechanical removal of grease by ridging:
C if there is no open water left after advection, there cannot be any grease ice
if ((1.0 _d 0 - AREApreTH(I,J)).LE.siEps) then
tmpscal3 = tmpscal3 + SItracer(I,J,bi,bj,iTrGrease)
SItracer(I,J,bi,bj,iTrGrease) = 0. _d 0
elseif (greaseNewFrazil .GT. 0. _d 0) then
C new ice growth goes into grease tracer not HEFF:
tmpscal3=0. _d 0
C
C 2) solidification of "old" grease ice
C (only when cold enough for ice growth)
C
C time scale dependent solidification
C (50% of grease ice area become solid ice within 24h):
tmpscal1=exp(-SEAICE_deltaTtherm/greaseDecayTime)
C gain in solid sea ice volume due to solidified grease:
d_HEFFbyGREASE(I,J) =
& SItracer(I,J,bi,bj,iTrGrease)
& * (1.0 _d 0 - tmpscal1)
C ... and related loss of grease ice (tracer) volume
SItracer(I,J,bi,bj,iTrGrease) =
& SItracer(I,J,bi,bj,iTrGrease) * tmpscal1
C the solidified grease ice volume needs to be added to HEFF:
SItrBucket(I,J,bi,bj,iTrGrease) =
& SItrBucket(I,J,bi,bj,iTrGrease)
& + d_HEFFbyGREASE(I,J)
C
C 3) grease ice growth from new frazil ice:
C
SItracer(i,j,bi,bj,iTrGrease) =
& SItracer(i,j,bi,bj,iTrGrease) + greaseNewFrazil
endif
C 4) mapping SItrBucket to external variable, in this case HEFF, ...
tmpscal3=tmpscal3+SItrBucket(I,J,bi,bj,iTrGrease)
C ... and empty SItrBucket for tracer 'grease'
SItrBucket(I,J,bi,bj,iTrGrease)=0. _d 0
#endif /* SEAICE_GREASE */
#ifdef SEAICE_ITD
C ice growth in open water adds to first category
d_HEFFbyATMonOCN_open_ITD(I,J,1)=tmpscal3
d_HEFFbyATMonOCN_ITD(I,J,1) =d_HEFFbyATMonOCN_ITD(I,J,1)
& +tmpscal3
#endif /* SEAICE_ITD */
d_HEFFbyATMonOCN_open(I,J)=tmpscal3
d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
ENDDO
ENDDO
#ifdef ALLOW_SITRACER
DO J=1,sNy
DO I=1,sNx
C needs to be here to allow use also with LEGACY branch
#ifdef SEAICE_ITD
SItrHEFF(I,J,bi,bj,4)=SItrHEFF(I,J,bi,bj,3)
& +d_HEFFbyATMonOCN_open(I,J)
#else /* ndef SEAICE_ITD */
SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
#endif /* SEAICE_ITD */
ENDDO
ENDDO
#endif /* ALLOW_SITRACER */
C convert snow to ice if submerged.
C =================================
C note: in legacy, this process is done at the end
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( SEAICEuseFlooding ) THEN
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
tmpscal3 = HEFFITDpreTH(I,J,IT)
& + d_HEFFbySublim_ITD(I,J,IT)
& + d_HEFFbyOCNonICE_ITD(I,J,IT)
& + d_HEFFbyATMonOCN_ITD(I,J,IT)
tmpscal4 = HSNWITDpreTH(I,J,IT)
& + d_HSNWbySublim_ITD(I,J,IT)
& + d_HSNWbyATMonSNW_ITD(I,J,IT)
& + d_HSNWbyRAIN_ITD(I,J,IT)
tmpscal0 = (tmpscal4*SEAICE_rhoSnow
& + tmpscal3*SEAICE_rhoIce)
& * recip_rhoConst
tmpscal1 = MAX( 0. _d 0, tmpscal0 - tmpscal3)
d_HEFFbyFLOODING_ITD(I,J,IT) = tmpscal1
d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
& +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
d_HEFFbyFLOODING(I,J)=tmpscal1
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 /* SEAICE_ITD */
ENDIF
#ifdef SEAICE_ITD
C apply ice and snow thickness changes
C =================================================================
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)
& + d_HEFFbySublim_ITD(I,J,IT)
& + d_HEFFbyOCNonICE_ITD(I,J,IT)
& + d_HEFFbyATMonOCN_ITD(I,J,IT)
& + d_HEFFbyFLOODING_ITD(I,J,IT)
HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
& + d_HSNWbySublim_ITD(I,J,IT)
& + d_HSNWbyATMonSNW_ITD(I,J,IT)
& + d_HSNWbyRAIN_ITD(I,J,IT)
& + d_HSNWbyOCNonSNW_ITD(I,J,IT)
& - d_HEFFbyFLOODING_ITD(I,J,IT)
& * ICE2SNOW
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ITD */
C ===================================================================
C ==========PART 4: determine ice cover fraction increments=========-
C ===================================================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE recip_heffActual = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_hsnwbyatmonsnw = comlev1_bibj,key=iicekey,byte=isbyte
cph(
cphCADJ STORE d_AREAbyATM = comlev1_bibj,key=iicekey,byte=isbyte
cphCADJ STORE d_AREAbyICE = comlev1_bibj,key=iicekey,byte=isbyte
cphCADJ STORE d_AREAbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
cph)
CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef SEAICE_ITD
C-- in thinnest category account for lateral ice growth and melt the
C-- "non-ITD" way, so that the ITD simulation with SEAICE_multDim=1
C-- is identical with the non-ITD simulation;
C-- use HEFF, ARE, HSNOW, etc. as temporal storage for 1st category
DO J=1,sNy
DO I=1,sNx
HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
AREApreTH(I,J)=AREAITDpreTH(I,J,1)
recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
#ifdef SEAICE_GREASE
C grease ice layer thickness (Hg) includes
C 1 part frazil ice and 3 parts sea water
C i.e. HO = 0.25 * Hg
recip_HO=4. _d 0 / greaseLayerThick(I,J)
#else /* SEAICE_GREASE */
IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
recip_HO=1. _d 0 / HO_south
ELSE
recip_HO=1. _d 0 / HO
ENDIF
#endif /* SEAICE_GREASE */
recip_HH = recip_heffActual(I,J)
C gain of ice over open water : computed from
#ifdef SEAICE_GREASE
C from growth by ATM with grease ice time delay
tmpscal4 = MAX(ZERO,d_HEFFbyGREASE(I,J))
#else /* SEAICE_GREASE */
C (SEAICE_areaGainFormula.EQ.1) from growth by ATM
C (SEAICE_areaGainFormula.EQ.2) from predicted growth by ATM
IF (SEAICE_areaGainFormula.EQ.1) THEN
tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
ELSE
tmpscal4=MAX(ZERO,a_QbyATM_open(I,J))
ENDIF
#endif /* SEAICE_GREASE */
C loss of ice cover by melting : computed from
C (SEAICE_areaLossFormula.EQ.1) from all but only melt conributions by ATM and OCN
C (SEAICE_areaLossFormula.EQ.2) from net melt-growth>0 by ATM and OCN
C (SEAICE_areaLossFormula.EQ.3) from predicted melt by ATM
IF (SEAICE_areaLossFormula.EQ.1) THEN
tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J) )
& + MIN( 0. _d 0 , d_HEFFbyATMonOCN_open(I,J) )
& + MIN( 0. _d 0 , d_HEFFbyOCNonICE(I,J) )
ELSEIF (SEAICE_areaLossFormula.EQ.2) THEN
tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J)
& + d_HEFFbyATMonOCN_open(I,J) + d_HEFFbyOCNonICE(I,J) )
ELSE
C compute heff after ice melt by ocn:
tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
C compute available heat left after snow melt by atm:
tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
& - d_HSNWbyATMonSNW(I,J)*SNOW2ICE
C could not melt more than all the ice
tmpscal2 = MAX(-tmpscal0,tmpscal1)
tmpscal3 = MIN(ZERO,tmpscal2)
ENDIF
C apply tendency
IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
& (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
AREA(I,J,bi,bj)=MAX(0. _d 0,
& MIN( SEAICE_area_max, AREA(I,J,bi,bj)
& + recip_HO*tmpscal4+HALF*recip_HH*tmpscal3
& * areaPDFfac ))
ELSE
AREA(I,J,bi,bj)=0. _d 0
ENDIF
#ifdef ALLOW_SITRACER
SItrAREA(I,J,bi,bj,3)=AREA(I,J,bi,bj)
#endif /* ALLOW_SITRACER */
#ifdef ALLOW_DIAGNOSTICS
d_AREAbyATM(I,J)=
& recip_HO*MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
& +HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_open(I,J))
& *areaPDFfac
d_AREAbyICE(I,J)=
& HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(I,J))
& *areaPDFfac
d_AREAbyOCN(I,J)=
& HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
& *areaPDFfac
#endif /* ALLOW_DIAGNOSTICS */
ENDDO
ENDDO
#ifdef SEAICE_ITD
C transfer 1st category values back into ITD variables
DO J=1,sNy
DO I=1,sNx
HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
ENDDO
ENDDO
C now melt ice laterally in all other thickness categories
C (areal growth, i.e. new ice formation, only occurrs in 1st category)
IF (SEAICE_multDim .gt. 1) THEN
DO IT=2,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
C when thickness is zero, area should be zero, too:
AREAITD(I,J,IT,bi,bj)=ZERO
ELSE
C tmpscal1 is the minimal ice concentration after lateral melt that will
C not lead to an unphysical increase of ice thickness by lateral melt;
C estimated as the concentration before thermodynamics scaled by the
C ratio of new ice thickness and ice thickness before thermodynamics
IF ( HEFFITDpreTH(I,J,IT).LE.ZERO ) THEN
tmpscal1=0. _d 0
ELSE
tmpscal1=AREAITDpreTH(I,J,IT)*
& HEFFITD(I,J,IT,bi,bj)/HEFFITDpreTH(I,J,IT)
ENDIF
C melt ice laterally based on an average floe sice
C following Steele (1992)
AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
& * (ONE - latMeltFrac(I,J,IT))
CML not necessary:
CML AREAITD(I,J,IT,bi,bj) = max(ZERO,AREAITD(I,J,IT,bi,bj))
C limit area reduction so that actual ice thickness does not increase
AREAITD(I,J,IT,bi,bj) = max(AREAITD(I,J,IT,bi,bj),
& tmpscal1)
ENDIF
#ifdef ALLOW_SITRACER
SItrAREA(I,J,bi,bj,3)=SItrAREA(I,J,bi,bj,3)
& +AREAITD(I,J,IT,bi,bj)
#endif /* ALLOW_SITRACER */
ENDDO
ENDDO
ENDDO
ENDIF
#endif /* SEAICE_ITD */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
Cgf 'bulk' linearization of area=f(HEFF)
IF ( SEAICEadjMODE.GE.1 ) THEN
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
& ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
DO J=1,sNy
DO I=1,sNx
C AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
AREA(I,J,bi,bj) = AREApreTH(I,J) + 0.1 _d 0 *
& ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
ENDDO
ENDDO
#endif /* SEAICE_ITD */
ENDIF
#endif
#ifdef SEAICE_ITD
C check categories for consistency with limits after growth/melt ...
IF ( SEAICEuseLinRemapITD ) CALL SEAICE_ITD_REMAP(
I heffitdPreTH, areaitdPreTH,
I bi, bj, myTime, myIter, myThid )
CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
C ... and update total AREA, HEFF, HSNOW
C (the updated HEFF is used below for ice salinity increments)
CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
#endif /* SEAICE_ITD */
#ifdef SEAICE_GREASE
C convert SItracer 'grease' from grease ice volume back to ratio:
C ===============================================================
DO J=1,sNy
DO I=1,sNx
if (HEFF(I,J,bi,bj).GT.siEps) then
SItracer(I,J,bi,bj,iTrGrease) =
& SItracer(I,J,bi,bj,iTrGrease) / HEFF(I,J,bi,bj)
else
SItracer(I,J,bi,bj,iTrGrease) = 0. _d 0
endif
ENDDO
ENDDO
#endif /* SEAICE_GREASE */
C ===================================================================
C =============PART 5: determine ice salinity increments=============
C ===================================================================
#ifndef SEAICE_VARIABLE_SALINITY
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE d_HEFFbyNEG(:,:,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
# ifdef ALLOW_SALT_PLUME
CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbySublim = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
# endif /* ALLOW_SALT_PLUME */
# endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
tmpscal1 = d_HEFFbyNEG(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J) +
& d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
& + d_HEFFbySublim(I,J)
#ifdef EXF_SEAICE_FRACTION
& + d_HEFFbyRLX(I,J,bi,bj)
#endif
Catn: can not take out more that surface salinity when SSS
tmpscal3 = max( 0. _d 0,
& min(SEAICE_salt0,salt(I,J,kSurface,bi,bj)) )
tmpscal2 = tmpscal1 * tmpscal3 * HEFFM(I,J,bi,bj)
& * recip_deltaTtherm * SEAICE_rhoIce
saltFlux(I,J,bi,bj) = tmpscal2
#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_SPLIT_BASIN
catn attempt to split East/West basins in Arctic
localSPfrac(I,J) = SPsalFRAC(1)
IF ( SaltPlumeSplitBasin ) THEN
localSPfrac(I,J) = SPsalFRAC(2)
IF(YC(I,J,bi,bj).LT. 85.0 .AND. YC(I,J,bi,bj).GT. 71.0
& .AND. XC(I,J,bi,bj) .LT. -90.0) THEN
localSPfrac(I,J) = SPsalFRAC(1)
ENDIF
ENDIF
#else
localSPfrac(I,J) = SPsalFRAC
#endif /* SALT_PLUME_SPLIT_BASIN */
#ifdef SALT_PLUME_IN_LEADS
Catn: Only d_HEFFbyATMonOCN should contribute to plume.
C By redefining tmpscal1 here, saltPlumeFlux is smaller in case
C define inLeads than case undef inLeads. Physical interpretation
C is that when d_HEFF is formed from below via ocean freezing, it
C occurs more uniform over grid cell and not inLeads, thus not
C participating in pkg/salt_plume.
C Note: tmpscal1 is defined only after saltFlux is calculated.
IceGrowthRateInLeads(I,J)=max(0. _d 0,d_HEFFbyATMonOCN(I,J))
tmpscal1 = IceGrowthRateInLeads(I,J)
leadPlumeFraction(I,J) =
& (ONE + EXP( (SPinflectionPoint - AREApreTH(I,J))*5.0
& /(ONE - SPinflectionPoint) ))**(-ONE)
localSPfrac(I,J)=localSPfrac(I,J)*leadPlumeFraction(I,J)
#endif /* SALT_PLUME_IN_LEADS */
tmpscal3 = tmpscal1*salt(I,J,kSurface,bi,bj)*HEFFM(I,J,bi,bj)
& * recip_deltaTtherm * SEAICE_rhoIce
saltPlumeFlux(I,J,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
& *localSPfrac(I,J)
C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
IF ( .NOT. SaltPlumeSouthernOcean ) THEN
IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
& saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
ENDIF
#endif /* ALLOW_SALT_PLUME */
ENDDO
ENDDO
#endif /* ndef SEAICE_VARIABLE_SALINITY */
#ifdef SEAICE_VARIABLE_SALINITY
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
C sum up the terms that affect the salt content of the ice pack
tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
C recompute HEFF before thermodynamic updates (which is not AREApreTH in legacy code)
tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
C tmpscal1 > 0 : m of sea ice that is created
IF ( tmpscal1 .GE. 0.0 ) THEN
saltFlux(I,J,bi,bj) =
& HEFFM(I,J,bi,bj)*recip_deltaTtherm
& *SEAICE_saltFrac*salt(I,J,kSurface,bi,bj)
& *tmpscal1*SEAICE_rhoIce
#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_SPLIT_BASIN
catn attempt to split East/West basins in Arctic
localSPfrac(I,J) = SPsalFRAC(1)
IF ( SaltPlumeSplitBasin ) THEN
localSPfrac(I,J) = SPsalFRAC(2)
IF(YC(I,J,bi,bj).LT. 85.0 .AND. YC(I,J,bi,bj).GT. 71.0
& .AND. XC(I,J,bi,bj) .LT. -90.0) THEN
localSPfrac(I,J) = SPsalFRAC(1)
ENDIF
ENDIF
#else
localSPfrac(I,J) = SPsalFRAC
#endif /* SALT_PLUME_SPLIT_BASIN */
#ifndef SALT_PLUME_IN_LEADS
C saltPlumeFlux is defined only during freezing:
saltPlumeFlux(I,J,bi,bj)=
& HEFFM(I,J,bi,bj)*recip_deltaTtherm
& *(ONE-SEAICE_saltFrac)*salt(I,J,kSurface,bi,bj)
& *tmpscal1*SEAICE_rhoIce
& *localSPfrac(I,J)
#endif /* ndef SALT_PLUME_IN_LEADS */
#endif /* ALLOW_SALT_PLUME */
C tmpscal1 < 0 : m of sea ice that is melted
ELSE
saltFlux(I,J,bi,bj) =
& HEFFM(I,J,bi,bj)*recip_deltaTtherm
& *HSALT(I,J,bi,bj)
& *tmpscal1/tmpscal2
#ifdef ALLOW_SALT_PLUME
#ifndef SALT_PLUME_IN_LEADS
saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
#endif /* ndef SALT_PLUME_IN_LEADS */
#endif /* ALLOW_SALT_PLUME */
ENDIF
#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_IN_LEADS
Catn: only d_HEFFbyATMonOCN should contribute to plume
C By redefining tmpscal1 here, saltPlumeFlux is smaller in case
C define inLeads than case undef inLeads. Physical interpretation
C is that when d_HEFF is formed from below via ocean freezing, it
C occurs more uniform over grid cell and not inLeads, thus not
C participating in pkg/salt_plume.
C Note: tmpscal1 is defined only after saltFlux is calculated.
IceGrowthRateInLeads(I,J)=max(0. _d 0,d_HEFFbyATMonOCN(I,J))
tmpscal1 = IceGrowthRateInLeads(I,J)
leadPlumeFraction(I,J) =
& (ONE + EXP( (SPinflectionPoint - AREApreTH(I,J))*5.0
& /(ONE - SPinflectionPoint) ))**(-ONE)
localSPfrac(I,J)=localSPfrac(I,J)*leadPlumeFraction(I,J)
IF ( tmpscal1 .GE. 0.0) THEN
C saltPlumeFlux is defined only during freezing:
saltPlumeFlux(I,J,bi,bj)=
& HEFFM(I,J,bi,bj)*recip_deltaTtherm
& *(ONE-SEAICE_saltFrac)*salt(I,J,kSurface,bi,bj)
& *tmpscal1*SEAICE_rhoIce
& *localSPfrac(I,J)
ELSE
saltPlumeFlux(I,J,bi,bj) = 0. _d 0
ENDIF
#endif /* SALT_PLUME_IN_LEADS */
C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
IF ( .NOT. SaltPlumeSouthernOcean ) THEN
IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
& saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
ENDIF
#endif /* ALLOW_SALT_PLUME */
C update HSALT based on surface saltFlux
HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
& saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
saltFlux(I,J,bi,bj) =
& saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J,bi,bj)
ENDDO
ENDDO
#endif /* SEAICE_VARIABLE_SALINITY */
#ifdef ALLOW_SITRACER
DO J=1,sNy
DO I=1,sNx
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
#endif /* ALLOW_SITRACER */
C ===================================================================
C ==============PART 7: determine ocean model forcing================
C ===================================================================
C compute net heat flux leaving/entering the ocean,
C accounting for the part used in melt/freeze processes
C =====================================================
#ifdef SEAICE_ITD
C compute total of "mult" fluxes for ocean forcing
DO J=1,sNy
DO I=1,sNx
a_QbyATM_cover(I,J) = 0.0 _d 0
r_QbyATM_cover(I,J) = 0.0 _d 0
a_QSWbyATM_cover(I,J) = 0.0 _d 0
r_FWbySublim(I,J) = 0.0 _d 0
ENDDO
ENDDO
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
C if fluxes in W/m^2 then use:
c a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
c & + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
c r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
c & + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
c a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
c & + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
c r_FWbySublim(I,J)=r_FWbySublim(I,J)
c & + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
C if fluxes in effective ice meters, i.e. ice volume per area, then use:
a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
& + a_QbyATMmult_cover(I,J,IT)
r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
& + r_QbyATMmult_cover(I,J,IT)
a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
& + a_QSWbyATMmult_cover(I,J,IT)
r_FWbySublim(I,J)=r_FWbySublim(I,J)
& + r_FWbySublimMult(I,J,IT)
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE d_hsnwbyneg(:,:,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
& + a_QSWbyATM_cover(I,J)
& - ( d_HEFFbyOCNonICE(I,J)
& + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
& + d_HEFFbyNEG(I,J,bi,bj)
#ifdef EXF_SEAICE_FRACTION
& + d_HEFFbyRLX(I,J,bi,bj)
#endif
& + d_HSNWbyNEG(I,J,bi,bj)*SNOW2ICE
& - 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
C switch heat fluxes from 'effective' ice meters to W/m2
C ======================================================
DO J=1,sNy
DO I=1,sNx
QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
ENDDO
ENDDO
#ifndef SEAICE_DISABLE_HEATCONSFIX
C treat advective heat flux by ocean to ice water exchange (at 0decC)
C ===================================================================
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE d_HEFFbyNEG(:,:,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HSNWbyNEG(:,:,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
CADJ STORE d_HSNWbyOCNonSNW = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE d_HSNWbyATMonSNW = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ & key = iicekey, byte = isbyte
# endif /* ALLOW_AUTODIFF_TAMC */
Cgf 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 arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
C at a different temperature, it leads to a loss of conservation in the
C ocean+ice system. While this is mostly a serious issue in the
C real fresh water + non linear free surface framework, a mismatch
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
catn: initialize tmpscal1
tmpscal1 = ZERO
C 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,bi,bj)+ d_HSNWbyNEG(I,J,bi,bj)*SNOW2ICE )
& * convertHI2PRECIP
& - snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J)) )
C factor in the heat content as done in external_forcing_surf.F
IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
& useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
tmpscal1 = - tmpscal3*
& HeatCapacity_Cp * temp_EvPrRn
ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
& useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
tmpscal1 = - tmpscal3*
& HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
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 in all cases, diagnoze the boundary condition mismatch to SIaaflux
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
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
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*(
& a_QSWbyATM_cover(I,J) +
& 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 ==================================================
DO J=1,sNy
DO I=1,sNx
tmpscal1= d_HSNWbyATMonSNW(I,J)*SNOW2ICE
& +d_HFRWbyRAIN(I,J)
& +d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
& +d_HEFFbyOCNonICE(I,J)
& +d_HEFFbyATMonOCN(I,J)
& +d_HEFFbyNEG(I,J,bi,bj)
#ifdef EXF_SEAICE_FRACTION
& +d_HEFFbyRLX(I,J,bi,bj)
#endif
& +d_HSNWbyNEG(I,J,bi,bj)*SNOW2ICE
C If r_FWbySublim>0, then it is evaporated from ocean.
& +r_FWbySublim(I,J)
EmPmR(I,J,bi,bj) = 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 */
& + tmpscal1*convertHI2PRECIP
& )*rhoConstFresh
#ifdef SEAICE_ITD
C beware of the sign: fw2ObyRidge is snow mass moved into the ocean
C by ridging, so requires a minus sign
& - fw2ObyRidge(I,J,bi,bj)*recip_deltaTtherm
& * maskC(I,J,kSurface,bi,bj)
#endif /* SEAICE_ITD */
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--
DO J=1,sNy
DO I=1,sNx
frWtrAtm(I,J,bi,bj) = 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
# ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
& - a_FWbySublim(I,J)*AREApreTH(I,J)
# endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
ENDDO
ENDDO
C--
#else /* ndef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
C--
# ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
DO J=1,sNy
DO I=1,sNx
frWtrAtm(I,J,bi,bj) = 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
# endif
C--
#endif /* ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
#ifdef SEAICE_DEBUG
IF ( plotLevel.GE.debLevC ) THEN
CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
ENDIF
#endif /* SEAICE_DEBUG */
C Sea Ice Load on the sea surface.
C =================================
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( useRealFreshWaterFlux ) THEN
DO J=1,sNy
DO I=1,sNx
#ifdef SEAICE_CAP_ICELOAD
tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
& + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
#else
tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
& + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
#endif
sIceLoad(i,j,bi,bj) = tmpscal2
ENDDO
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 /* ALLOW_BALANCE_FLUXES */
C ===================================================================
C ======================PART 8: diagnostics==========================
C ===================================================================
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
tmpscal1=1. _d 0 * recip_deltaTtherm
CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_cover,
& tmpscal1,1,'SIaQbATC',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_open,
& tmpscal1,1,'SIaQbATO',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(a_QbyOCN,
& tmpscal1,1,'SIaQbOCN',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyOCNonICE,
& tmpscal1,1,'SIdHbOCN',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_cover,
& tmpscal1,1,'SIdHbATC',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_open,
& tmpscal1,1,'SIdHbATO',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
& tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyOCNonSNW,
& tmpscal1,1,'SIdSbOCN',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyATMonSNW,
& tmpscal1,1,'SIdSbATC',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyATM,
& tmpscal1,1,'SIdAbATO',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyICE,
& tmpscal1,1,'SIdAbATC',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyOCN,
& tmpscal1,1,'SIdAbOCN',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_open,
& convertHI2Q,1, 'SIqneto ',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_cover,
& convertHI2Q,1, 'SIqneti ',0,1,3,bi,bj,myThid)
C three that actually need intermediate storage
DO J=1,sNy
DO I=1,sNx
DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
& * d_HSNWbyRAIN(I,J)*SEAICE_rhoSnow*recip_deltaTtherm
DIAGarrayB(I,J) = AREA(I,J,bi,bj)-AREApreTH(I,J)
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(DIAGarrayA,
& 'SIsnPrcp',0,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_SCALE_FILL(DIAGarrayB,
& tmpscal1,1,'SIdA ',0,1,3,bi,bj,myThid)
DO J=1,sNy
DO I=1,sNx
DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
& a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(DIAGarrayB,
& 'SIfwSubl',0,1,3,bi,bj,myThid)
C
DO J=1,sNy
DO I=1,sNx
C the actual Freshwater flux of sublimated ice, >0 decreases ice
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
DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
& * r_FWbySublim(I,J)
& * SEAICE_rhoIce * recip_deltaTtherm
C the latent heat flux
tmpscal1= EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
& + r_FWbySublim(I,J)*convertHI2PRECIP
tmpscal2= ( a_FWbySublim(I,J)-r_FWbySublim(I,J) )
& * convertHI2PRECIP
tmpscal3= SEAICE_lhEvap+SEAICE_lhFusion
DIAGarrayB(I,J) = -maskC(I,J,kSurface,bi,bj)*rhoConstFresh
& * ( tmpscal1*SEAICE_lhEvap + tmpscal2*tmpscal3 )
ENDDO
ENDDO
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)
# ifdef SEAICE_GREASE
C actual grease ice layer thickness
CALL DIAGNOSTICS_FILL(greaseLayerThick,
& 'SIgrsLT ',0,1,3,bi,bj,myThid)
# endif /* SEAICE_GREASE */
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
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 /* ALLOW_BALANCE_FLUXES */
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
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
#endif /* ALLOW_DIAGNOSTICS */
#else /* ALLOW_EXF and ALLOW_ATM_TEMP */
STOP 'SEAICE_GROWTH not compiled without EXF and ALLOW_ATM_TEMP'
#endif /* ALLOW_EXF and ALLOW_ATM_TEMP */
RETURN
END