--- MITgcm/pkg/thsice/thsice_step_fwd.F 2004/07/12 20:20:03 1.6 +++ MITgcm/pkg/thsice/thsice_step_fwd.F 2006/04/09 17:35:30 1.16 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.6 2004/07/12 20:20:03 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.16 2006/04/09 17:35:30 heimbach Exp $ C $Name: $ #include "THSICE_OPTIONS.h" @@ -6,14 +6,14 @@ 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 @@ -30,14 +30,14 @@ #include "THSICE_PARAMS.h" #include "THSICE_VARS.h" #include "THSICE_TAVE.h" - + C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C bi,bj :: tile indices C iMin,iMax :: computation domain: 1rst index range C jMin,jMax :: computation domain: 2nd index range C- input: -C prcAtm :: total precip from the atmosphere [kg/m2/s] +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). @@ -73,11 +73,11 @@ 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) +C isIceFree :: true for ice-free grid-cell that remains ice-free 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 @@ -86,33 +86,48 @@ _RL Fbot, esurp _RL opFrac, icFrac _RL oceV2s, oceTs - _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr) + _RL Tsf, Tice(nlyr), qicen(nlyr) _RL tmpflx(0:2), tmpdTs + _RL hSnow +c + _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL compact(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL hIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy) +c + LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy) +#ifdef ALLOW_DIAGNOSTICS + _RL tmpFac +#endif LOGICAL dBug - 1010 FORMAT(A,1P4E11.3) + 1010 FORMAT(A,1P4E14.6) dBug = .FALSE. C- Initialise flxAtm DO j = 1-Oly, sNy+Oly DO i = 1-Olx, sNx+Olx flxAtm(i,j) = 0. + isIceFree(i,j) = .FALSE. +#ifdef ALLOW_AUTODIFF_TAMC + compact(i,j) = 0. + hIce(i,j) = 0. +#endif ENDDO ENDDO - IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN + IF ( fluidIsWater ) THEN DO j = jMin, jMax DO i = iMin, iMax c dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 ) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C part.1 : ice-covered fraction ; +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 + 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) + hIce(i,j) = iceHeight(i,j,bi,bj) hSnow = snowHeight(i,j,bi,bj) Tsf = Tsrf(i,j,bi,bj) qicen(1)= Qice1(i,j,bi,bj) @@ -124,13 +139,14 @@ ENDIF CALL THSICE_ALBEDO( - I hIce, hSnow, Tsf, snowAge(i,j,bi,bj), + I hIce(i,j), hSnow, Tsf, snowAge(i,j,bi,bj), O albedo, I myThid ) flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo) + siceAlb(i,j,bi,bj) = albedo CALL THSICE_SOLVE4TEMP( - I useBulkforce, tmpflx, TFrzOce, hIce, hSnow, + I useBulkForce, tmpflx, TFrzOce, hIce(i,j), 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), @@ -147,13 +163,9 @@ 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*/ IF ( dBug ) THEN - WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =', - & Tsf, Tice, frzmltMxL + WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), TFrzOce =', + & Tsf, Tice, TFrzOce WRITE(6,1010) 'ThSI_FWD: sHeat,fxCndBt, fxAtm,evAtm=', & sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), & flxAtm(i,j), evpAtm(i,j) @@ -164,9 +176,26 @@ ENDIF dBug = .FALSE. +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + tmpFac = 1. _d 0 + 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 +#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 @@ -182,22 +211,22 @@ Fbot = 0. _d 0 saltFlux(i,j,bi,bj) = 0. _d 0 - compact= iceMask(i,j,bi,bj) + compact(i,j)= iceMask(i,j,bi,bj) C------- - IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN + IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact(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 =', - & compact, iceHeight(i,j,bi,bj), - & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj) - WRITE(6,1010) 'ThSI_FWD: ocTs, TFrzOce, frzmltMxL =', - & oceTs, TFrzOce, frzmltMxL + & compact(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=', + & oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj) ENDIF C------- - IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN + 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) + hIce(i,j) = iceHeight(i,j,bi,bj) hSnow = snowHeight(i,j,bi,bj) Tsf = Tsrf(i,j,bi,bj) qicen(1)= Qice1(i,j,bi,bj) @@ -207,7 +236,7 @@ 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, + U compact(i,j), hIce(i,j), hSnow, Tsf, qicen, flx2oc, O frw2oc, fsalt, Fbot, I dBug, myThid) @@ -216,10 +245,10 @@ snowPrc(i,j,bi,bj) = snowPr C-- Snow aging : - snowAge(i,j,bi,bj) = thSIce_deltaT + 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) + IF ( snowPr.GT.0. _d 0 ) + & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj) & * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge ) C-- @@ -233,11 +262,11 @@ 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 + & compact(i,j),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 icFrac, compact(i,j), hIce(i,j), hSnow, qicen, I flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm, I myTime, myIter, myThid ) #endif /* CHECK_ENERGY_CONSERV */ @@ -245,7 +274,7 @@ C-- Update Sea-Ice state : c iceMask(i,j,bi,bj)=compact - iceHeight(i,j,bi,bj) = hIce + iceHeight(i,j,bi,bj) = hIce(i,j) snowHeight(i,j,bi,bj)= hSnow Tsrf(i,j,bi,bj) =Tsf Qice1(i,j,bi,bj)=qicen(1) @@ -264,7 +293,8 @@ IF (dBug) WRITE(6,1010) & 'ThSI_FWD:-3- compact, hIc, hSn, Qnet =', - & compact,hIce,hSnow,Qnet(i,j,bi,bj) + & compact(i,j),hIce(i,j),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) @@ -278,55 +308,59 @@ 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) +c compact(i,j) = iceMask(i,j,bi,bj) + hIce(i,j) = 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 + icFrac = compact(i,j) 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, + U oceTs, compact(i,j), hIce(i,j), 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. + & ,compact(i,j),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 icFrac, compact(i,j), hIce(i,j), 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 + IF ( compact(i,j).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 + iceHeight(i,j,bi,bj) = hIce(i,j) snowHeight(i,j,bi,bj)= hSnow 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 + 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) + & 'ThSI_FWD:-4- compact, hIc, hSn, Qnet =', + & compact(i,j),hIce(i,j),hSnow, + & Qnet(i,j,bi,bj) C-- - if esurp > 0 : end ENDIF - IF ( compact .GT. 0. _d 0 ) THEN - iceMask(i,j,bi,bj)=compact + IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) + & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0 + & .AND. compact(i,j) .LE.0. _d 0 + IF ( compact(i,j) .GT. 0. _d 0 ) THEN + iceMask(i,j,bi,bj)=compact(i,j) IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0 ELSE iceMask(i,j,bi,bj) = 0. _d 0 @@ -354,6 +388,14 @@ ENDDO ENDDO +#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 */