C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.7 2004/07/22 22:52:59 jmc Exp $ C $Name: $ #include "THSICE_OPTIONS.h" CBOP C !ROUTINE: THSICE_STEP_FWD C !INTERFACE: 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 | o Step Forward Therm-SeaIce model. C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "FFIELDS.h" #include "THSICE_SIZE.h" #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 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- 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. 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 CEOP #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 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 _RL opFrac, icFrac _RL oceV2s, oceTs _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr) _RL tmpflx(0:2), tmpdTs LOGICAL dBug 1010 FORMAT(A,1P4E11.3) dBug = .FALSE. C- Initialise flxAtm DO j = 1-Oly, sNy+Oly DO i = 1-Olx, sNx+Olx flxAtm(i,j) = 0. ENDDO ENDDO 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.11 ) 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) IF ( dBug ) THEN WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj WRITE(6,1010) 'ThSI_FWD:-0- iceMask, hIc, hSn, Tsf =', & icFrac, hIce,hSnow,Tsf ENDIF 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) #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*/ IF ( dBug ) THEN WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), frzmltMxL =', & Tsf, Tice, frzmltMxL 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) ENDIF ENDIF ENDDO ENDDO ENDIF dBug = .FALSE. C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 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.11 ) TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj) oceTs = tOceMxL(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) C------- IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.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), 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 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-|--+----| 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) C-- Net fluxes : frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj)) C- weighted average net fluxes: 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) ELSE flxAtm(i,j) = 0. _d 0 frwAtm = 0. _d 0 ENDIF 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 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 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 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 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 ENDIF C-- Return atmospheric fluxes in evpAtm & flxSW (same sign and units): evpAtm(i,j) = frwAtm flxSW (i,j) = flxAtm(i,j) #ifdef ATMOSPHERIC_LOADING 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) #endif ENDDO ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #endif /* ALLOW_THSICE */ RETURN END