--- MITgcm/pkg/thsice/thsice_step_fwd.F 2004/04/07 23:40:34 1.3 +++ MITgcm/pkg/thsice/thsice_step_fwd.F 2010/12/17 04:00:14 1.30 @@ -1,19 +1,21 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.3 2004/04/07 23:40:34 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.30 2010/12/17 04:00:14 gforget Exp $ C $Name: $ #include "THSICE_OPTIONS.h" +#ifdef ALLOW_ATM2D +# include "ctrparam.h" +#endif CBOP C !ROUTINE: THSICE_STEP_FWD C !INTERFACE: - SUBROUTINE THSICE_STEP_FWD( + SUBROUTINE THSICE_STEP_FWD( I bi, bj, iMin, iMax, jMin, jMax, I prcAtm, - U evpAtm, flxSW, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | S/R THSICE_STEP_FWD +C | S/R THSICE_STEP_FWD C | o Step Forward Therm-SeaIce model. C *==========================================================* C \ev @@ -26,33 +28,45 @@ #include "EEPARAMS.h" #include "PARAMS.h" #include "FFIELDS.h" +#ifdef ALLOW_ATM2D +# include "ATMSIZE.h" +# include "ATM2D_VARS.h" +#endif #include "THSICE_SIZE.h" #include "THSICE_PARAMS.h" #include "THSICE_VARS.h" #include "THSICE_TAVE.h" - +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +# include "tamc_keys.h" +#endif + + INTEGER siLo, siHi, sjLo, sjHi + PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx ) + PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy ) + C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === +C- input: C bi,bj :: tile indices C iMin,iMax :: computation domain: 1rst index range C jMin,jMax :: computation domain: 2nd index range +C prcAtm :: total precip from the atmosphere [kg/m2/s] +C myTime :: current Time of simulation [s] +C myIter :: current Iteration number in simulation +C myThid :: my Thread Id number +C-- Use fluxes hold in commom blocks C- input: -C prcAtm :: total precip from the atmosphere [kg/m2/s] -C evpAtm :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate) -C flxSW :: (Inp) short-wave heat flux (+=down): downward comp. only -C (part.1), becomes net SW flux into ocean (part.2). +C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean +C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down) +C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up) C- output -C evpAtm :: (Out) net fresh-water flux (E-P) from the atmosphere [m/s] (+=up) -C flxSW :: (Out) net surf. heat flux from the atmosphere [W/m2], (+=down) -C myTime :: time counter for this thread -C myIter :: iteration counter for this thread -C myThid :: thread number for this instance of the routine. +C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down) +C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up) INTEGER bi,bj INTEGER iMin, iMax INTEGER jMin, jMax _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL myTime INTEGER myIter INTEGER myThid @@ -61,269 +75,330 @@ #ifdef ALLOW_THSICE C !LOCAL VARIABLES: C === Local variables === -C snowPr :: snow precipitation [kg/m2/s] -C agingTime :: aging time scale (s) -C ageFac :: snow aging factor [1] -C albedo :: surface albedo [0-1] -C flxAtm :: net heat flux from the atmosphere (+=down) [W/m2] -C frwAtm :: net fresh-water flux (E-P) to the atmosphere [kg/m2/s] -C Fbot :: the oceanic heat flux already incorporated (ice_therm) +C iceFrac :: fraction of grid area covered in ice C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2] C frw2oc :: fresh-water flux from the ice to the ocean C fsalt :: mass salt flux to the ocean C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2] -C TFrzOce :: sea-water freezing temperature [oC] (function of S) - INTEGER i,j - _RL snowPr - _RL agingTime, ageFac - _RL albedo - _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL frwAtm - _RL flx2oc - _RL frw2oc - _RL fsalt - _RL TFrzOce, cphm, frzmltMxL - _RL Fbot, esurp +C tFrzOce :: sea-water freezing temperature [oC] (function of S) +C isIceFree :: true for ice-free grid-cell that remains ice-free +C ageFac :: snow aging factor [1] +C snowFac :: snowing refreshing-age factor [units of 1/snowPr] + LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL ageFac + _RL snowFac + _RL cphm _RL opFrac, icFrac - _RL oceV2s, oceTs - _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr) - _RL tmpflx(0:2), tmpdTs - - LOGICAL dBug +#ifdef ALLOW_DIAGNOSTICS + _RL tmpFac +#endif + INTEGER i,j + LOGICAL dBugFlag - dBug = .FALSE. - 1010 FORMAT(A,1P4E11.3) +C- define grid-point location where to print debugging values +#include "THSICE_DEBUG.h" - IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN - DO j = jMin, jMax - DO i = iMin, iMax -c dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.15 ) + 1010 FORMAT(A,1P4E14.6) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C part.1 : ice-covered fraction ; -C Solve for surface and ice temperature (implicitly) ; compute surf. fluxes -C------- - IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN - icFrac = iceMask(i,j,bi,bj) - TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj) - hIce = iceHeight(i,j,bi,bj) - hSnow = snowHeight(i,j,bi,bj) - Tsf = Tsrf(i,j,bi,bj) - qicen(1)= Qice1(i,j,bi,bj) - qicen(2)= Qice2(i,j,bi,bj) - - CALL THSICE_ALBEDO( - I hIce, hSnow, Tsf, snowAge(i,j,bi,bj), - O albedo, - I myThid ) - flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo) - - CALL THSICE_SOLVE4TEMP( - I useBulkforce, tmpflx, TFrzOce, hIce, hSnow, - U flxSW(i,j), Tsf, qicen, - O Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), - O tmpdTs, flxAtm(i,j), evpAtm(i,j), - I i,j, bi,bj, myThid) -#ifdef SHORTWAVE_HEATING -C-- Update Fluxes : - opFrac= 1. _d 0-icFrac - Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj) +#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 + ticekey = (act1 + 1) + act2*max1 + & + act3*max1*max2 + & + act4*max1*max2*max3 +#endif /* ALLOW_AUTODIFF_TAMC */ + +C- Initialise + dBugFlag = debugLevel.GE.debLevB + DO j = 1-OLy, sNy+OLy + DO i = 1-OLx, sNx+OLx + isIceFree(i,j) = .FALSE. +#ifdef ALLOW_ATM2D + sFluxFromIce(i,j) = 0. _d 0 +#else + saltFlux(i,j,bi,bj) = 0. _d 0 +#endif +#ifdef ALLOW_AUTODIFF_TAMC + iceFrac(i,j) = 0. #endif -C-- Update Sea-Ice state : - Tsrf(i,j,bi,bj) =Tsf - Tice1(i,j,bi,bj)=Tice(1) - Tice2(i,j,bi,bj)=Tice(2) - Qice1(i,j,bi,bj)=qicen(1) - Qice2(i,j,bi,bj)=qicen(2) -#ifdef ALLOW_TIMEAVE - ice_albedo_Ave(i,j,bi,bj) = ice_albedo_Ave(i,j,bi,bj) - & + icFrac*albedo*thSIce_deltaT -#endif /*ALLOW_TIMEAVE*/ - ENDIF ENDDO + ENDDO + + ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime + snowFac = thSIce_deltaT/(rhos*hNewSnowAge) + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +#endif + DO j = jMin, jMax + DO i = iMin, iMax + IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN +C-- Snow aging : + snowAge(i,j,bi,bj) = thSIce_deltaT + & + snowAge(i,j,bi,bj)*ageFac + IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 ) + & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj) + & * EXP( - snowFac*snowPrc(i,j,bi,bj) ) +c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos) +c & /hNewSnowAge ) +C------- +C note: Any flux of mass (here fresh water) that enter or leave the system +C with a non zero energy HAS TO be counted: add snow precip. + icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj) + & - Lfresh*snowPrc(i,j,bi,bj) +C-- + ENDIF ENDDO + ENDDO + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + tmpFac = 1. _d 0 + CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid) + CALL DIAGNOSTICS_FRACT_FILL( + I snowPrc, iceMask,tmpFac,1,'SIsnwPrc', + I 0,1,1,bi,bj,myThid) + CALL DIAGNOSTICS_FRACT_FILL( + I siceAlb, iceMask,tmpFac,1,'SIalbedo', + I 0,1,1,bi,bj,myThid) ENDIF - dBug = .FALSE. +#endif /* ALLOW_DIAGNOSTICS */ + DO j = jMin, jMax + DO i = iMin, iMax + siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj) + ENDDO + ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C part.2 : ice-covered fraction ; -C change in ice/snow thickness and ice-fraction +C part.2 : ice-covered fraction ; +C change in ice/snow thickness and ice-fraction C note: can only reduce the ice-fraction but not increase it. C------- - agingTime = 50. _d 0 * 86400. _d 0 - ageFac = 1. _d 0 - thSIce_deltaT/agingTime DO j = jMin, jMax DO i = iMin, iMax -c dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.15 ) - TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj) - oceTs = tOceMxL(i,j,bi,bj) + tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj) cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj) - frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT - - Fbot = 0. _d 0 - saltFlux(i,j,bi,bj) = 0. _d 0 - compact= iceMask(i,j,bi,bj) + frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) ) + & * cphm/ocean_deltaT + iceFrac(i,j) = iceMask(i,j,bi,bj) + flx2oc(i,j) = icFlxSW(i,j,bi,bj) C------- - IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN - WRITE(6,1010) 'ThSI_FWD:-1- iceMask,hIc,hSn,Qnet=', - & compact, hIce, hSnow, Qnet(i,j,bi,bj) - WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL=', - & oceTs,TFrzOce,frzmltMxL +#ifdef ALLOW_DBUG_THSICE + IF ( dBug(i,j,bi,bj) ) THEN + IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN + WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj + WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =', + & iceFrac(i,j), iceHeight(i,j,bi,bj), + & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj) + WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=', + & tOceMxL(i,j,bi,bj), tFrzOce(i,j), + & frzmltMxL(i,j), Qnet(i,j,bi,bj) + ENDIF + IF (iceFrac(i,j).GT.0.) + & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=', + & iceFrac(i,j), icFlxAtm(i,j,bi,bj), + & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj) ENDIF -C------- - IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN - - oceV2s = v2ocMxL(i,j,bi,bj) - snowPr = snowPrc(i,j,bi,bj) - hIce = iceHeight(i,j,bi,bj) - hSnow = snowHeight(i,j,bi,bj) - Tsf = Tsrf(i,j,bi,bj) - qicen(1)= Qice1(i,j,bi,bj) - qicen(2)= Qice2(i,j,bi,bj) - flx2oc = flxSW(i,j) - - CALL THSICE_CALC_THICKN( - I frzmltMxL, TFrzOce, oceTs, oceV2s, snowPr, - I sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), evpAtm(i,j), - U compact, hIce, hSnow, Tsf, qicen, flx2oc, - O frw2oc, fsalt, Fbot, - I dBug, myThid) - -C- note : snowPr was not supposed to be modified in THSICE_THERM ; -C but to reproduce old results, is reset to zero if Tsf >= 0 - snowPrc(i,j,bi,bj) = snowPr - -C-- Snow aging : - snowAge(i,j,bi,bj) = thSIce_deltaT - & + snowAge(i,j,bi,bj)*ageFac - IF ( snowPr.GT.0. _d 0 ) - & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj) - & * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge ) -C-- - -C-- Diagnostic of Atmospheric Fluxes over sea-ice : - frwAtm = evpAtm(i,j) - prcAtm(i,j) -C note: Any flux of mass (here fresh water) that enter or leave the system -C with a non zero energy HAS TO be counted: add snow precip. - flxAtm(i,j) = flxAtm(i,j) - Lfresh*snowPrc(i,j,bi,bj) - -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=', - & iceMask(i,j,bi,bj),flxAtm(i,j),evpAtm(i,j),-Lfresh*snowPr - IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=', - & compact,flx2oc,fsalt,frw2oc -#ifdef CHECK_ENERGY_CONSERV - icFrac = iceMask(i,j,bi,bj) - CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0, - I icFrac, compact, hIce, hSnow, qicen, - I flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm, - I myTime, myIter, myThid ) -#endif /* CHECK_ENERGY_CONSERV */ -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +#endif + ENDDO + ENDDO -C-- Update Sea-Ice state : -c iceMask(i,j,bi,bj)=compact - iceheight(i,j,bi,bj) = hIce - snowheight(i,j,bi,bj)= hSnow - Tsrf(i,j,bi,bj) =Tsf - Qice1(i,j,bi,bj)=qicen(1) - Qice2(i,j,bi,bj)=qicen(2) +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +#endif + CALL THSICE_CALC_THICKN( + I bi, bj, + I iMin,iMax, jMin,jMax, dBugFlag, + I iceMask(siLo,sjLo,bi,bj), tFrzOce, + I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj), + I snowPrc(siLo,sjLo,bi,bj), prcAtm, + I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj), + U iceFrac, iceHeight(siLo,sjLo,bi,bj), + U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj), + U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj), + U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc, + O frw2oc, fsalt, + I myTime, myIter, myThid ) + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE fsalt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE flx2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte +CADJ STORE frw2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte +#endif C-- Net fluxes : - frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj)) + DO j = jMin, jMax + DO i = iMin, iMax +c#ifdef ALLOW_AUTODIFF_TAMC +c ikey_1 = i +c & + sNx*(j-1) +c & + sNx*sNy*act1 +c & + sNx*sNy*max1*act2 +c & + sNx*sNy*max1*max2*act3 +c & + sNx*sNy*max1*max2*max3*act4 +c#endif /* ALLOW_AUTODIFF_TAMC */ +c#ifdef ALLOW_AUTODIFF_TAMC +cCADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1 +c#endif + IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN C- weighted average net fluxes: +c#ifdef ALLOW_AUTODIFF_TAMC +cCADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1 +cCADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1 +cCADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1 +cCADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1 +c#endif icFrac = iceMask(i,j,bi,bj) opFrac= 1. _d 0-icFrac - flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj) - frwAtm = icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj) - Qnet(i,j,bi,bj)=-icFrac*flx2oc +opFrac*Qnet(i,j,bi,bj) - EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj) - saltFlux(i,j,bi,bj)=-icFrac*fsalt - - IF (dBug) WRITE(6,1010)'ThSI_FWD:-3- compact,hIc,hSn,Qnet=', - & compact,hIce,hSnow,Qnet(i,j,bi,bj) - - ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN - flxAtm(i,j) = -Qnet(i,j,bi,bj) - frwAtm = rhofw*EmPmR(i,j,bi,bj) +#ifdef ALLOW_ATM2D + pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j) + pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw + sFluxFromIce(i,j) = -icFrac*fsalt(i,j) +#else + icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj) + & - opFrac*Qnet(i,j,bi,bj) + icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj) + & + opFrac*EmPmR(i,j,bi,bj) + Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj) + EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j) + & + opFrac*EmPmR(i,j,bi,bj) + saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j) +#endif + +#ifdef ALLOW_DBUG_THSICE + IF (dBug(i,j,bi,bj)) WRITE(6,1010) + & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =', + & iceFrac(i,j), iceHeight(i,j,bi,bj), + & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj) +#endif + + ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN + icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj) + icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj) ELSE - flxAtm(i,j) = 0. _d 0 - frwAtm = 0. _d 0 + icFlxAtm(i,j,bi,bj) = 0. _d 0 + icFrwAtm(i,j,bi,bj) = 0. _d 0 ENDIF + ENDDO + ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C part.3 : freezing of sea-water C over ice-free fraction and what is left from ice-covered fraction C------- -c compact= iceMask(i,j,bi,bj) - hIce = iceHeight(i,j,bi,bj) - hSnow = snowHeight(i,j,bi,bj) - - esurp = frzmltMxL - Fbot*iceMask(i,j,bi,bj) - IF (esurp.GT.0. _d 0) THEN - icFrac = compact - qicen(1)= Qice1(i,j,bi,bj) - qicen(2)= Qice2(i,j,bi,bj) - CALL THSICE_EXTEND( - I esurp, TFrzOce, - U oceTs, compact, hIce, hSnow, qicen, - O flx2oc, frw2oc, fsalt, - I dBug, myThid ) -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=' - & ,compact,flx2oc,fsalt,frw2oc -#ifdef CHECK_ENERGY_CONSERV - tmpflx(1) = 0. - tmpflx(2) = 0. - CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1, - I icFrac, compact, hIce, hSnow, qicen, - I flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2), - I myTime, myIter, myThid ) -#endif /* CHECK_ENERGY_CONSERV */ -C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C-- Update Sea-Ice state : - IF ( compact.GT.0. _d 0 .AND. icFrac.EQ.0. _d 0) THEN - Tsrf(i,j,bi,bj) = TFrzOce - Tice1(i,j,bi,bj) = TFrzOce - Tice2(i,j,bi,bj) = TFrzOce - Qice1(i,j,bi,bj) = qicen(1) - Qice2(i,j,bi,bj) = qicen(2) - ENDIF - iceheight(i,j,bi,bj) = hIce - snowheight(i,j,bi,bj)= hSnow + CALL THSICE_EXTEND( + I bi, bj, + I iMin,iMax, jMin,jMax, dBugFlag, + I frzmltMxL, tFrzOce, + I tOceMxL(siLo,sjLo,bi,bj), + U iceFrac, iceHeight(siLo,sjLo,bi,bj), + U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj), + U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj), + U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj), + O flx2oc, frw2oc, fsalt, + I myTime, myIter, myThid ) + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE snowHeight(:,:,bi,bj) = +CADJ & comlev1_bibj,key=ticekey,byte=isbyte +#endif + DO j = jMin, jMax + DO i = iMin, iMax + IF (frzmltMxL(i,j).GT.0. _d 0) THEN C-- Net fluxes : - Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc - EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw - saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt - - IF (dBug) WRITE(6,1010)'ThSI_FWD:-4- compact,hIc,hSn,Qnet=', - & compact,hIce,hSnow,Qnet(i,j,bi,bj) -C-- - if esurp > 0 : end +#ifdef ALLOW_ATM2D + pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j) + pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw + sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j) +#else + Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j) + EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j) + saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j) +#endif + +#ifdef ALLOW_DBUG_THSICE + IF (dBug(i,j,bi,bj)) WRITE(6,1010) + & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =', + & iceFrac(i,j), iceHeight(i,j,bi,bj), + & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj) +#endif ENDIF - IF ( compact .GT. 0. _d 0 ) THEN - iceMask(i,j,bi,bj)=compact - IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0 + IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) + & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0 + & .AND. iceFrac(i,j) .LE.0. _d 0 + IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN + iceMask(i,j,bi,bj)=iceFrac(i,j) + IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 ) + & snowAge(i,j,bi,bj) = 0. _d 0 ELSE iceMask(i,j,bi,bj) = 0. _d 0 iceHeight(i,j,bi,bj)= 0. _d 0 snowHeight(i,j,bi,bj)=0. _d 0 snowAge(i,j,bi,bj) = 0. _d 0 - Tsrf(i,j,bi,bj) = oceTs + Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj) Tice1(i,j,bi,bj) = 0. _d 0 Tice2(i,j,bi,bj) = 0. _d 0 - Qice1(i,j,bi,bj) = 0. _d 0 - Qice2(i,j,bi,bj) = 0. _d 0 + Qice1(i,j,bi,bj) = Lfresh + Qice2(i,j,bi,bj) = Lfresh ENDIF + ENDDO + ENDDO -C-- Return atmospheric fluxes in evpAtm & flxSW (same sign and units): - evpAtm(i,j) = frwAtm - flxSW (i,j) = flxAtm(i,j) +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE snowHeight(:,:,bi,bj) = +CADJ & comlev1_bibj,key=ticekey,byte=isbyte +# endif + DO j = jMin, jMax + DO i = iMin, iMax +C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit) + sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos + & + iceHeight(i,j,bi,bj)*rhoi + & )*iceMask(i,j,bi,bj) ENDDO ENDDO + IF ( thSIceAdvScheme.GT.0 ) THEN +C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux + DO j = jMin, jMax + DO i = iMin, iMax + IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN + Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj) + EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj) + saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj) + ENDIF + ENDDO + ENDDO + ENDIF + +#ifdef ALLOW_BULK_FORCE + IF ( useBulkForce ) THEN + CALL BULKF_FLUX_ADJUST( + I bi, bj, iMin, iMax, jMin, jMax, + I isIceFree, myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_BULK_FORCE */ + C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #endif /* ALLOW_THSICE */