C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/torge/itd/code/seaice_growth.F,v 1.2 2012/04/27 22:25:23 dimitri Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_EXF # include "EXF_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 C !FUNCTIONS: #ifdef ALLOW_DIAGNOSTICS LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON #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 constants _RL tempFrz, ICE2SNOW, SNOW2ICE _RL QI, QS, recip_QI C-- TmixLoc :: ocean surface/mixed-layer temperature (in K) _RL TmixLoc (1:sNx,1:sNy) 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 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 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 #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 #ifdef SEAICE_ITD c The change of mean ice area due to out-of-bounds values following c sea ice dynamics _RL d_AREAbyNEG (1:sNx,1:sNy) #endif 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 _RL d_HEFFbyOCNonICE (1:sNx,1:sNy) c The sum of mean ice thickness increments due to atmospheric fluxes over the open water c 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 water c 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_HSNWbyNEG (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 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,MULTDIM) #endif 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 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 heffFracFactor (1:sNx,1:sNy,1:nITD) #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 CToM<<< C INTEGER it INTEGER IT, K C>>>ToM _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 _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 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 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 */ 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 #ifdef SEAICE_ALLOW_AREA_RELAXATION d_AREAbyRLX(I,J) = 0.0 _d 0 d_HEFFbyRLX(I,J) = 0.0 _d 0 #endif #ifdef SEAICE_ITD d_AREAbyNEG(I,J) = 0.0 _d 0 #endif d_HEFFbyNEG(I,J) = 0.0 _d 0 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_HSNWbyNEG(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 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 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_ITD r_QbyATMmult_cover (I,J,IT) = 0.0 _d 0 r_FWbySublimMult(I,J,IT) = 0.0 _d 0 #endif #ifdef SEAICE_CAP_SUBLIM latentHeatFluxMaxMult(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 ===================================================================== #ifdef SEAICE_GROWTH_LEGACY DO J=1,sNy DO I=1,sNx HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj) HSNWpreTH(I,J)=HSNOW(I,J,bi,bj) AREApreTH(I,J)=AREANM1(I,J,bi,bj) d_HEFFbyNEG(I,J)=0. _d 0 d_HSNWbyNEG(I,J)=0. _d 0 #ifdef ALLOW_DIAGNOSTICS DIAGarrayA(I,J) = AREANM1(I,J,bi,bj) DIAGarrayB(I,J) = AREANM1(I,J,bi,bj) DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj) DIAGarrayD(I,J) = HSNOW(I,J,bi,bj) #endif ENDDO ENDDO #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj) HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj) AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj) IF (AREA(I,J,bi,bj).GT.0) THEN areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj) ELSE areaFracFactor(I,J,K)=ZERO ENDIF IF (HEFF(I,J,bi,bj).GT.0) THEN heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj) ELSE heffFracFactor(I,J,K)=ZERO ENDIF ENDDO ENDDO ENDDO #endif #else /* SEAICE_GROWTH_LEGACY */ #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ) Cgf no dependency through pathological cases treatment IF ( SEAICEadjMODE.EQ.0 ) THEN #endif #ifdef SEAICE_ALLOW_AREA_RELAXATION 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 IF (SEAICE_tauAreaObsRelax .GT. 0.) THEN DO J=1,sNy DO I=1,sNx C d_AREAbyRLX(i,j) = 0. _d 0 C d_HEFFbyRLX(i,j) = 0. _d 0 IF ( obsSIce(I,J,bi,bj).GT.AREA(I,J,bi,bj)) THEN d_AREAbyRLX(i,j) = & SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax & * (obsSIce(I,J,bi,bj) - AREA(I,J,bi,bj)) ENDIF IF ( obsSIce(I,J,bi,bj).GT.0. _d 0 .AND. & AREA(I,J,bi,bj).EQ.0. _d 0) THEN C d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j) d_HEFFbyRLX(i,j) = 1. _d 1 * siEps ENDIF #ifdef SEAICE_ITD AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj) & + d_AREAbyRLX(i,j) HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj) & + d_HEFFbyRLX(i,j) #endif AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j) HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j) ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_AREA_RELAXATION */ C 1) treat the case of negative values: #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 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx #ifdef SEAICE_ITD DO K=1,nITD tmpscal2=0. _d 0 tmpscal3=0. _d 0 tmpscal4=0. _d 0 tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0) HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2 tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0) HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3 tmpscal4=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0) AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4 d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal4 ENDDO AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyNEG(I,J) #else d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0) d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0) AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0) #endif HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J) HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J) ENDDO ENDDO C 1.25) treat the case of very thin ice: #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx #ifdef SEAICE_ITD DO K=1,nITD tmpscal2=0. _d 0 tmpscal3=0. _d 0 IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN tmpscal2=-HEFFITD(I,J,K,bi,bj) tmpscal3=-HSNOWITD(I,J,K,bi,bj) TICES(I,J,K,bi,bj)=celsius2K HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2 HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3 c TICE(I,J,bi,bj)=celsius2K c HEFF(I,J,bi,bj) =HEFF(I,J,bi,bj) +tmpscal2 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3 ENDIF ENDDO #else tmpscal2=0. _d 0 tmpscal3=0. _d 0 IF (HEFF(I,J,bi,bj).LE.siEps) THEN tmpscal2=-HEFF(I,J,bi,bj) tmpscal3=-HSNOW(I,J,bi,bj) TICE(I,J,bi,bj)=celsius2K DO IT=1,SEAICE_multDim TICES(I,J,IT,bi,bj)=celsius2K ENDDO ENDIF HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3 #endif ENDDO ENDDO C 1.5) treat the case of area but no ice/snow: #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 */ DO J=1,sNy DO I=1,sNx CToM<<< IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND. C & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0 & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) THEN AREA(I,J,bi,bj)=0. _d 0 #ifdef SEAICE_ITD DO K=1,nITD AREAITD(I,J,K,bi,bj)=0. _d 0 HEFFITD(I,J,K,bi,bj)=0. _d 0 HSNOWITD(I,J,K,bi,bj)=0. _d 0 ENDDO #endif ENDIF C>>>ToM ENDDO ENDDO C 2) treat the case of very small area: #ifndef DISABLE_AREA_FLOOR #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN #ifdef SEAICE_ITD tmpscal2=AREA(I,J,bi,bj) #endif AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor) #ifdef SEAICE_ITD c ice area added (tmpscal3 is .ge.0): tmpscal3=AREA(I,J,bi,bj)-tmpscal2 c distribute this gain proportionally over categories DO K=1,nITD tmpscal4=AREAITD(I,J,K,bi,bj)/tmpscal2*tmpscal3 AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4 ENDDO #endif ENDIF ENDDO ENDDO #endif /* DISABLE_AREA_FLOOR */ C 2.5) treat case of excessive ice cover, e.g., due to ridging: #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx #ifdef SEAICE_ITD tmpscal2=AREA(I,J,bi,bj) #endif #ifdef ALLOW_DIAGNOSTICS DIAGarrayA(I,J) = AREA(I,J,bi,bj) #endif #ifdef ALLOW_SITRACER SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj) #endif AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max) #ifdef SEAICE_ITD c ice area subtracted (tmpscal3 is .ge.0): tmpscal3=tmpscal2-AREA(I,J,bi,bj) c distribute this loss proportionally over categories DO K=1,nITD AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj) & -tmpscal3*areaFracFactor(I,J,K) ENDDO #endif ENDDO ENDDO #ifdef SEAICE_ITD C If AREAITD is changed due to regularization (but HEFFITD not) then the C actual ice thickness (HEFFITD/AREAITD) in a category can be changed so C that it does not fit its category limits anymore and redistribution is necessary CALL SEAICE_ITD_REDIST(myTime, myIter, myThid) C this should not affect the respective sums (AREA, HEFF, ...) C ... except a non-conserving redistribution scheme is used; then call: c CALL SEAICE_ITD_SUM(myTime, myIter, myThid) #endif #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ) ENDIF #endif C 3) 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 K=1,nITD DO J=1,sNy DO I=1,sNx HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj) HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj) AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj) ENDDO ENDDO ENDDO #endif C 4) treat sea ice salinity pathological cases #ifdef SEAICE_VARIABLE_SALINITY #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR. & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) * & HSALT(I,J,bi,bj) * recip_deltaTtherm HSALT(I,J,bi,bj) = 0.0 _d 0 ENDIF ENDDO ENDDO #endif /* SEAICE_VARIABLE_SALINITY */ #endif /* SEAICE_GROWTH_LEGACY */ #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid) 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 K=1,nITD DO J=1,sNy DO I=1,sNx HEFFITDpreTH(I,J,K) = 0. _d 0 HSNWITDpreTH(I,J,K) = 0. _d 0 AREAITDpreTH(I,J,K) = 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 4) 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 K=1,nITD #endif DO J=1,sNy DO I=1,sNx #ifdef SEAICE_ITD IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN #ifdef SEAICE_GROWTH_LEGACY tmpscal1 = MAX(SEAICE_area_reg,AREAITDpreTH(I,J,K)) hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1 tmpscal2 = HEFFITDpreTH(I,J,K)/tmpscal1 heffActualMult(I,J,K) = MAX(tmpscal2,SEAICE_hice_reg) #else /* SEAICE_GROWTH_LEGACY */ cif regularize AREA with SEAICE_area_reg tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K) & + area_reg_sq) cif heffActual calculated with the regularized AREA tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1 cif regularize heffActual with SEAICE_hice_reg (add lower bound) heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2 & + hice_reg_sq) cif hsnowActual calculated with the regularized AREA hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1 #endif /* SEAICE_GROWTH_LEGACY */ cif regularize the inverse of heffActual by hice_reg recip_heffActualMult(I,J,K) = AREAITDpreTH(I,J,K) / & sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K) & + hice_reg_sq) cif Do not regularize when HEFFpreTH = 0 ELSE heffActualMult(I,J,K) = ZERO hsnowActualMult(I,J,K) = ZERO recip_heffActualMult(I,J,K) = ZERO ENDIF #else IF (HEFFpreTH(I,J) .GT. ZERO) THEN #ifdef SEAICE_GROWTH_LEGACY tmpscal1 = MAX(SEAICE_area_reg,AREApreTH(I,J)) hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1 tmpscal2 = HEFFpreTH(I,J)/tmpscal1 heffActual(I,J) = MAX(tmpscal2,SEAICE_hice_reg) #else /* SEAICE_GROWTH_LEGACY */ 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 #endif /* SEAICE_GROWTH_LEGACY */ 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 #endif ENDDO ENDDO #ifdef SEAICE_ITD ENDDO #endif #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 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE C AND SNOW THICKNESS #ifdef SEAICE_ITD DO K=1,nITD #endif DO J=1,sNy DO I=1,sNx 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 IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm * & (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce + & HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K) ELSE latentHeatFluxMaxMult(I,J,K) = ZERO ENDIF #else 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 #endif ENDDO ENDDO #ifdef SEAICE_ITD ENDDO #endif #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 ======================================================================= #ifdef ALLOW_ATM_WIND IF (useRelativeWind) 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 #endif /* ALLOW_ATM_WIND */ #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE tice(:,:,bi,bj) CADJ & = comlev1_bibj, key = iicekey, byte = isbyte 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 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 #endif DO J=1,sNy DO I=1,sNx #ifndef SEAICE_ITD CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above C (instead of heffActual and latentHeatFluxMax) heffActualMult(I,J,IT)= heffActual(I,J)*pFac #ifdef SEAICE_CAP_SUBLIM latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac #endif #endif ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj) ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj) TICE(I,J,bi,bj) = ZERO TICES(I,J,IT,bi,bj) = ZERO ENDDO ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE heffActualMult = 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( #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_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 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 TICE & TICES #ifdef SEAICE_ITD C calculate area weighted mean TICE(I,J,bi,bj) = TICE(I,J,bi,bj) & + ticeOutMult(I,J,IT)*areaFracFactor(I,J,K) #else TICE(I,J,bi,bj) = TICE(I,J,bi,bj) & + ticeOutMult(I,J,IT)*recip_multDim #endif TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT) C average over categories #ifdef SEAICE_ITD C calculate area weighted mean a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J) & + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K) a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J) & + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K) a_FWbySublim (I,J) = a_FWbySublim(I,J) & + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K) #else a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J) & + a_QbyATMmult_cover(I,J,IT)*recip_multDim a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J) & + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim a_FWbySublim (I,J) = a_FWbySublim(I,J) & + a_FWbySublimMult(I,J,IT)*recip_multDim #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 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 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 #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K) & * convertQ2HI * AREAITDpreTH(I,J,K) a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K) & * convertQ2HI * AREAITDpreTH(I,J,K) r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K) #ifdef SEAICE_DISABLE_SUBLIM a_FWbySublimMult(I,J,K) = ZERO #endif a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce & * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K) r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K) ENDDO ENDDO ENDDO #endif #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 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 #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx a_QbyATMmult_cover(I,J,K) = 0. _d 0 r_QbyATMmult_cover(I,J,K) = 0. _d 0 a_QSWbyATMmult_cover(I,J,K) = 0. _d 0 ENDDO ENDDO ENDDO #endif 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 #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 =================================================================== 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 K=1,nITD #endif DO J=1,sNy DO I=1,sNx C First sublimate/deposite snow tmpscal2 = #ifdef SEAICE_ITD & MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj) & *SNOW2ICE),ZERO) C accumulate change over ITD categories d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2 & *ICE2SNOW HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal2 & *ICE2SNOW r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2 C keep total up to date, too r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2 #else & 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 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,K),HEFFITD(I,J,K,bi,bj)),ZERO) d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) - tmpscal2 r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2 C keep total up to date, too r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2 #else & 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 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,K) = a_QbyATMmult_cover(I,J,K) & - r_FWbySublimMult(I,J,K) r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K) & - r_FWbySublimMult(I,J,K) ENDDO ENDDO C end K loop ENDDO C then update totals DO J=1,sNy DO I=1,sNx #endif 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) ENDDO ENDDO 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 */ #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx C ice growth/melt due to ocean heat is equally distributed under the ice C and hence weighted by fractional area of each thickness category tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K), & -HEFFITD(I,J,K,bi,bj)) HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1 d_HEFFbyOCNonICE(I,J)=d_HEFFbyOCNonICE(I,J) + tmpscal1 ENDDO ENDDO ENDDO #endif DO J=1,sNy DO I=1,sNx #ifndef SEAICE_ITD d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj)) #endif 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 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 K=1,nITD 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,K),-HSNOWITD(I,J,K,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)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2 C keep the total up to date, too r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2 ENDDO ENDDO ENDDO #else 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 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2 ENDDO ENDDO #endif /* SEAICE_ITD */ DO J=1,sNy DO I=1,sNx HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyATMonSNW(I,J) ENDDO ENDDO 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 K=1,nITD DO J=1,sNy DO I=1,sNx #ifdef SEAICE_GROWTH_LEGACY tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)) #else tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K) c Limit ice growth by potential melt by ocean & + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J)*areaFracFactor(I,J,K)) #endif /* SEAICE_GROWTH_LEGACY */ d_HEFFbyATMonOCN_cover(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 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2 ENDDO ENDDO ENDDO #else DO J=1,sNy DO I=1,sNx #ifdef SEAICE_GROWTH_LEGACY 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 & AREApreTH(I,J) * r_QbyOCN(I,J)) #endif /* SEAICE_GROWTH_LEGACY */ 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 ENDDO ENDDO #endif /* SEAICE_ITD */ DO J=1,sNy DO I=1,sNx HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +d_HEFFbyATMonOCN_cover(I,J) #ifdef ALLOW_SITRACER SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj) #endif ENDDO ENDDO C attribute precip to fresh water or snow stock, C depending on atmospheric conditions. C ================================================= #ifdef ALLOW_ATM_TEMP #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 */ 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 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 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J) ENDDO ENDDO #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) & + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K) ENDDO ENDDO ENDDO #endif 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. #endif /* ALLOW_ATM_TEMP */ 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 */ #ifdef SEAICE_ITD DO K=1,nITD DO J=1,sNy DO I=1,sNx tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K), & -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 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2 d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2 ENDDO ENDDO ENDDO #endif DO J=1,sNy DO I=1,sNx #ifndef SEAICE_ITD 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 #endif 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_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 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, & -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj) 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 #ifdef SEAICE_ITD C determine thickness of new ice C considering the entire open water area to refreeze tmpscal4 = tmpscal3/(ONE-AREApreTH(I,J)) C then add new ice volume to appropriate thickness category DO K=1,nITD IF (tmpscal4.LT.Hlimit(K)) THEN HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3 AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj) & + ONE-AREApreTH(I,J) ENDIF ENDDO C in this case no open water is left after this step AREA(I,J,bi,bj) = ONE #endif 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 SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj) ENDDO ENDDO #endif /* ALLOW_SITRACER */ C convert snow to ice if submerged. C ================================= #ifndef SEAICE_GROWTH_LEGACY 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 K=1,nITD DO J=1,sNy DO I=1,sNx tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow & +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj)) d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1 & * ICE2SNOW ENDDO ENDDO ENDDO #else 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 ENDDO ENDDO #endif DO J=1,sNy DO I=1,sNx 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 /* SEAICE_GROWTH_LEGACY */ 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 replaces Hibler '79 scheme and lead closing parameter C because ITD accounts explicitly for lead openings and C different melt rates due to varying ice thickness C C only consider ice area loss due to total ice thickness loss C ice area gain due to freezing of open water as handled above C under "gain of new ice over open water" C C does not account for lateral melt of ice floes C DO K=1,nITD DO J=1,sNy DO I=1,sNx IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN AREAITD(I,J,K,bi,bj)=ZERO ENDIF ENDDO ENDDO ENDDO C update total AREA, HEFF, HSNOW CALL SEAICE_ITD_SUM(myTime,myIter,myThid) #else DO J=1,sNy DO I=1,sNx IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN recip_HO=1. _d 0 / HO_south ELSE recip_HO=1. _d 0 / HO ENDIF #ifdef SEAICE_GROWTH_LEGACY tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J) recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0) #else recip_HH = recip_heffActual(I,J) #endif C gain of ice over open water : computed from 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 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 )) 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)) d_AREAbyICE(I,J)= & HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(I,J)) d_AREAbyOCN(I,J)= & HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) ) #endif /* ALLOW_DIAGNOSTICS */ ENDDO ENDDO #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 K=1,nITD DO J=1,sNy DO I=1,sNx AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 * & ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) ) ENDDO ENDDO ENDDO C update total AREA, HEFF, HSNOW CALL SEAICE_ITD_SUM(myTime,myIter,myThid) #else 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 ENDIF #endif C =================================================================== C =============PART 5: determine ice salinity increments============= C =================================================================== #ifndef SEAICE_VARIABLE_SALINITY # if (defined ALLOW_AUTODIFF_TAMC && defined ALLOW_SALT_PLUME) CADJ STORE d_HEFFbyNEG = comlev1_bibj,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_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_AUTODIFF_TAMC and ALLOW_SALT_PLUME */ DO J=1,sNy DO I=1,sNx 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 + d_HEFFbyRLX(I,J) #endif tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj) & * recip_deltaTtherm * SEAICE_rhoIce saltFlux(I,J,bi,bj) = tmpscal2 #ifdef ALLOW_SALT_PLUME 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) & *SPsalFRAC #endif /* ALLOW_SALT_PLUME */ ENDDO ENDDO #endif /* ndef SEAICE_VARIABLE_SALINITY */ #ifdef SEAICE_VARIABLE_SALINITY #ifdef SEAICE_GROWTH_LEGACY # 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 set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) * & HSALT(I,J,bi,bj) * recip_deltaTtherm HSALT(I,J,bi,bj) = 0.0 _d 0 ENDIF ENDDO ENDDO #endif /* SEAICE_GROWTH_LEGACY */ #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 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 & *SPsalFRAC 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 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 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0 #endif /* ALLOW_SALT_PLUME */ ENDIF 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) #ifdef SEAICE_GROWTH_LEGACY C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) - & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm HSALT(I,J,bi,bj) = 0.0 _d 0 #ifdef ALLOW_SALT_PLUME saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0 #endif /* ALLOW_SALT_PLUME */ ENDIF #endif /* SEAICE_GROWTH_LEGACY */ ENDDO ENDDO #endif /* SEAICE_VARIABLE_SALINITY */ C ======================================================================= C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding === C ======================================================================= #ifdef SEAICE_GROWTH_LEGACY C treat values of ice cover fraction oustide C the [0 1] range, and other such issues. C =========================================== Cgf note: this part cannot be heat and water conserving #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE area(:,:,bi,bj) = comlev1_bibj, CADJ & key = iicekey, byte = isbyte CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, CADJ & key = iicekey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx C NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE CML replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably CML meant to be something like a minimum thickness AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4) ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE area(:,:,bi,bj) = comlev1_bibj, CADJ & key = iicekey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx C NOW TRUNCATE AREA AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj)) ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE area(:,:,bi,bj) = comlev1_bibj, CADJ & key = iicekey, byte = isbyte CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, CADJ & key = iicekey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO J=1,sNy DO I=1,sNx AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj)) HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj)) AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj) HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj) #ifdef SEAICE_CAP_HEFF C This is not energy conserving, but at least it conserves fresh water tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0) d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0 #endif /* SEAICE_CAP_HEFF */ HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj) ENDDO ENDDO C convert snow to ice if submerged. C ================================= IF ( SEAICEuseFlooding ) THEN 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 #endif /* SEAICE_GROWTH_LEGACY */ #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 ALLOW_AUTODIFF_TAMC CADJ STORE d_hsnwbyneg = comlev1_bibj,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) #ifndef SEAICE_GROWTH_LEGACY C in principle a_QSWbyATM_cover should always be included here, however 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_HSNWbyOCNonSNW(I,J)*SNOW2ICE + & d_HEFFbyNEG(I,J) + #ifdef SEAICE_ALLOW_AREA_RELAXATION & d_HEFFbyRLX(I,J) + #endif & d_HSNWbyNEG(I,J)*SNOW2ICE ) & * maskC(I,J,kSurface,bi,bj) 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 = comlev1_bibj,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 = comlev1_bibj,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 */ IF ( SEAICEheatConsFix ) THEN 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 to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates c on external_forcing_surf.F and applies the correction to QNET. IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and c the small error of the synchronous time stepping case (see external_forcing_surf.F). DO J=1,sNy DO I=1,sNx #ifdef ALLOW_DIAGNOSTICS c store unaltered QNET for diagnostic purposes 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 c will associate with EMPMR, and remove it from QNET, so that c melt/freez water is in effect consistently gained/lost at 0degC IF (temp_EvPrRn.NE.UNSET_RL) THEN QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3* & HeatCapacity_Cp * temp_EvPrRn ELSE QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3* & HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) ENDIF #ifdef ALLOW_DIAGNOSTICS c back out the eventual TFLUX adjustement and fill diag DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J) #endif 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 net fresh water flux leaving/entering C the ocean, accounting for fresh/salt water stocks. C ================================================== #ifdef ALLOW_ATM_TEMP 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) #ifdef SEAICE_ALLOW_AREA_RELAXATION & +d_HEFFbyRLX(I,J) #endif & +d_HSNWbyNEG(I,J)*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 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 */ #endif /* ALLOW_ATM_TEMP */ #ifdef SEAICE_DEBUG 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 /* 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 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) #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 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) 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 #endif /* ALLOW_DIAGNOSTICS */ C close bi,bj loops ENDDO ENDDO RETURN END