C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_solve4temp.F,v 1.18 2009/09/14 20:54:33 jmc Exp $ C $Name: $ #include "THSICE_OPTIONS.h" CBOP C !ROUTINE: THSICE_SOLVE4TEMP C !INTERFACE: SUBROUTINE THSICE_SOLVE4TEMP( I bi, bj, siLo, siHi, sjLo, sjHi, I iMin,iMax, jMin,jMax, dBugFlag, I useBulkForce, useEXF, I iceMask, hIce, hSnow, tFrz, flxExSW, U flxSW, tSrf, qIc1, qIc2, O tIc1, tIc2, dTsrf, O sHeat, flxCnB, flxAtm, evpAtm, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R THSICE_SOLVE4TEMP C *==========================================================* C | Solve (implicitly) for sea-ice and surface temperature C *==========================================================* C \ev C ADAPTED FROM: C LANL CICE.v2.0.2 C----------------------------------------------------------------------- C.. thermodynamics (vertical physics) based on M. Winton 3-layer model C.. See Bitz, C. M. and W. H. Lipscomb, 1999: "An energy-conserving C.. thermodynamic sea ice model for climate study." J. Geophys. C.. Res., 104, 15669 - 15677. C.. Winton, M., 1999: "A reformulated three-layer sea ice model." C.. Submitted to J. Atmos. Ocean. Technol. C.. authors Elizabeth C. Hunke and William Lipscomb C.. Fluid Dynamics Group, Los Alamos National Laboratory C----------------------------------------------------------------------- Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc) C.. Compute temperature change using Winton model with 2 ice layers, of C.. which only the top layer has a variable heat capacity. C !USES: IMPLICIT NONE C == Global variables === #include "EEPARAMS.h" #include "THSICE_SIZE.h" #include "THSICE_PARAMS.h" #ifdef ALLOW_AUTODIFF_TAMC # include "SIZE.h" # include "tamc.h" # include "tamc_keys.h" #endif C !INPUT/OUTPUT PARAMETERS: C == Routine Arguments == C siLo,siHi :: size of input/output array: 1rst dim. lower,higher bounds C sjLo,sjHi :: size of input/output array: 2nd dim. lower,higher bounds C bi,bj :: tile indices C iMin,iMax :: computation domain: 1rst index range C jMin,jMax :: computation domain: 2nd index range C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point). C useBulkForce:: use surf. fluxes from bulk-forcing external S/R C useEXF :: use surf. fluxes from exf external S/R C--- Input: C iceMask :: sea-ice fractional mask [0-1] C hIce (hi) :: ice height [m] C hSnow (hs) :: snow height [m] C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S) C flxExSW (=) :: surf. heat flux (+=down) except SW, function of surf. temp Ts: C 0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n) C--- Modified (input&output): C flxSW (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface C (=) :: output= below sea-ice, into the ocean C tSrf (Tsf) :: surface (ice or snow) temperature C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level C--- Output C tIc1 (Tice) :: temperature of ice layer 1 [oC] C tIc2 (Tice) :: temperature of ice layer 2 [oC] C dTsrf (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n C sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction) C flxCnB (=) :: heat flux conducted through the ice to bottom surface C flxAtm (=) :: net flux of energy from the atmosphere [W/m2] (+=down) C without snow precip. (energy=0 for liquid water at 0.oC) C evpAtm (=) :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate) C--- Input: C myTime :: current Time of simulation [s] C myIter :: current Iteration number in simulation C myThid :: my Thread Id number INTEGER siLo, siHi, sjLo, sjHi INTEGER bi,bj INTEGER iMin, iMax INTEGER jMin, jMax LOGICAL dBugFlag LOGICAL useBulkForce LOGICAL useEXF _RL iceMask(siLo:siHi,sjLo:sjHi) _RL hIce (siLo:siHi,sjLo:sjHi) _RL hSnow (siLo:siHi,sjLo:sjHi) _RL tFrz (siLo:siHi,sjLo:sjHi) _RL flxExSW(iMin:iMax,jMin:jMax,0:2) _RL flxSW (siLo:siHi,sjLo:sjHi) _RL tSrf (siLo:siHi,sjLo:sjHi) _RL qIc1 (siLo:siHi,sjLo:sjHi) _RL qIc2 (siLo:siHi,sjLo:sjHi) _RL tIc1 (siLo:siHi,sjLo:sjHi) _RL tIc2 (siLo:siHi,sjLo:sjHi) c _RL dTsrf (siLo:siHi,sjLo:sjHi) _RL dTsrf (iMin:iMax,jMin:jMax) _RL sHeat (siLo:siHi,sjLo:sjHi) _RL flxCnB (siLo:siHi,sjLo:sjHi) _RL flxAtm (siLo:siHi,sjLo:sjHi) _RL evpAtm (siLo:siHi,sjLo:sjHi) _RL myTime INTEGER myIter INTEGER myThid CEOP #ifdef ALLOW_THSICE C !LOCAL VARIABLES: C--- local copy of input/output argument list variables (see description above) c _RL flxExcSw(0:2) _RL Tf _RL hi _RL hs _RL netSW _RL Tsf _RL qicen(nlyr) _RL Tice (nlyr) c _RL sHeating c _RL flxCnB _RL dTsf c _RL flxAtm c _RL evpAtm C == Local Variables == C frsnow :: fractional snow cover C fswpen :: SW penetrating beneath surface (W m-2) C fswdn :: SW absorbed at surface (W m-2) C fswint :: SW absorbed in ice (W m-2) C fswocn :: SW passed through ice to ocean (W m-2) C flxExceptSw :: net surface heat flux, except short-wave (W/m2) C evap :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate) C dEvdT :: derivative of evap. with respect to Tsf [kg/m2/s/K] C flx0 :: net surf heat flux, from Atmos. to sea-ice (W m-2) C fct :: heat conducted to top surface C df0dT :: deriv of flx0 wrt Tsf (W m-2 deg-1) C k12, k32 :: thermal conductivity terms C a10, b10 :: coefficients in quadratic eqn for T1 C a1, b1, c1 :: coefficients in quadratic eqn for T1 C Tsf_start :: old value of Tsf C dt :: timestep INTEGER i,j INTEGER k, iterMax _RL frsnow _RL fswpen _RL fswdn _RL fswint _RL fswocn _RL flxExceptSw _RL evap, dEvdT _RL flx0 _RL fct _RL df0dT _RL k12, k32 _RL a10, b10 _RL a1, b1, c1 c _RL Tsf_start _RL dt _RL recip_dhSnowLin INTEGER iceornot LOGICAL useBlkFlx C- define grid-point location where to print debugging values #include "THSICE_DEBUG.h" 1010 FORMAT(A,I3,3F11.6) 1020 FORMAT(A,1P4E14.6) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #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 #endif /* ALLOW_AUTODIFF_TAMC */ useBlkFlx = useEXF .OR. useBulkForce IF ( dhSnowLin.GT.0. _d 0 ) THEN recip_dhSnowLin = 1. _d 0 / dhSnowLin ELSE recip_dhSnowLin = 0. _d 0 ENDIF dt = thSIce_dtTemp DO j = jMin, jMax DO i = iMin, iMax #ifdef ALLOW_AUTODIFF_TAMC ikey_1 = i & + sNx*(j-1) & + sNx*sNy*act1 & + sNx*sNy*max1*act2 & + sNx*sNy*max1*max2*act3 & + sNx*sNy*max1*max2*max3*act4 #endif /* ALLOW_AUTODIFF_TAMC */ C-- #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE devdt = comlev1_thsice_1, key=ikey_1 CADJ STORE df0dt = comlev1_thsice_1, key=ikey_1 CADJ STORE flxexceptsw = comlev1_thsice_1, key=ikey_1 CADJ STORE flxsw(i,j) = comlev1_thsice_1, key=ikey_1 CADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1 CADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1 CADJ STORE tsrf(i,j) = comlev1_thsice_1, key=ikey_1 #endif IF ( iceMask(i,j).GT.0. _d 0) THEN hi = hIce(i,j) hs = hSnow(i,j) Tf = tFrz(i,j) netSW = flxSW(i,j) Tsf = tSrf(i,j) qicen(1)= qIc1(i,j) qicen(2)= qIc2(i,j) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)') & 'ThSI_SOLVE4T: i,j=',i,j,bi,bj #endif IF ( hi.LT.hIceMin ) THEN C If hi < hIceMin, melt the ice. STOP 'THSICE_SOLVE4TEMP: should not enter if hi 0 ; it=', myIter, qicen(1) WRITE (standardMessageUnit,'(A,4I5,2F11.4)') & ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice ENDIF IF ( Tice(2).GT.0. _d 0) THEN WRITE (standardMessageUnit,'(A,I12,1PE14.6)') & ' BBerr: Tice(2) > 0 ; it=', myIter, qicen(2) WRITE (standardMessageUnit,'(A,4I5,2F11.4)') & ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice ENDIF #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1010) & 'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice #endif C Compute coefficients used in quadratic formula. a10 = rhoi*cpIce *hi/(2. _d 0*dt) + & k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi) & / (6. _d 0*dt*k32 + rhoi*cpIce *hi) b10 = -hi* & (rhoi*cpIce*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1)) & /(2. _d 0*dt) & - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2)) & / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C Compute new surface and internal temperatures; iterate until C Tsfc converges. IF ( useBlkFlx ) THEN iterMax = nitMaxTsf ELSE iterMax = 1 ENDIF dTsf = Terrmax C ----- begin iteration ----- DO k = 1,iterMax #ifdef ALLOW_AUTODIFF_TAMC ikey_3 = (ikey_1-1)*MaxTsf + k #endif #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE tsf = comlev1_thsice_3, key=ikey_3 CADJ STORE dtsf = comlev1_thsice_3, key=ikey_3 CADJ STORE df0dt = comlev1_thsice_3, key=ikey_3 CADJ STORE flxexceptsw = comlev1_thsice_3, key=ikey_3 #endif IF ( ABS(dTsf).GE.Terrmax ) THEN C Save temperatures at start of iteration. c Tsf_start = Tsf IF ( useBlkFlx ) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE tsf = comlev1_thsice_3, key=ikey_3 #endif C Compute top surface flux. IF (hs.GT.3. _d -1) THEN iceornot=2 ELSE iceornot=1 ENDIF IF ( useBulkForce ) THEN CALL THSICE_GET_BULKF( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) ELSEIF ( useEXF ) THEN CALL THSICE_GET_EXF ( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) ENDIF ELSE flxExceptSw = flxExSW(i,j,1) df0dT = flxExSW(i,j,2) ENDIF flx0 = fswdn + flxExceptSw #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1020) & 'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT #endif C Compute new top layer and surface temperatures. C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1 C with different coefficients. #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE tsf = comlev1_thsice_3, key=ikey_3 #endif a1 = a10 - k12*df0dT / (k12-df0dT) b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT) Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1) dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT) Tsf = Tsf + dTsf IF (Tsf .GT. 0. _d 0) THEN #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1010) & 'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf #endif a1 = a10 + k12 C note: b1 = b10 - k12*Tf0 b1 = b10 Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1) Tsf = 0. _d 0 IF ( useBlkFlx ) THEN IF (hs.GT.3. _d -1) THEN iceornot=2 ELSE iceornot=1 ENDIF IF ( useBulkForce ) THEN CALL THSICE_GET_BULKF( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) ELSEIF ( useEXF ) THEN CALL THSICE_GET_EXF ( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) ENDIF dTsf = 0. _d 0 ELSE flxExceptSw = flxExSW(i,j,0) dTsf = 1000. df0dT = 0. ENDIF flx0 = fswdn + flxExceptSw ENDIF C Check for convergence. If no convergence, then repeat. C C Convergence test: Make sure Tsfc has converged, within prescribed error. C (Energy conservation is guaranteed within machine roundoff, even C if Tsfc has not converged.) C If no convergence, then repeat. #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1010) & 'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf #endif IF ( useBlkFlx .AND. k.EQ.nitMaxTsf & .AND. ABS(dTsf).GE.Terrmax ) THEN WRITE (6,'(A,4I4,I12,F15.9)') & ' BB: not converge: i,j,it,hi=',i,j,bi,bj, & myIter,hi WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT IF (Tsf.LT.-70. _d 0) STOP ENDIF ENDIF ENDDO C ------ end iteration ------------ C Compute new bottom layer temperature. #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE Tice(:) = comlev1_thsice_1, key=ikey_1 CADJ STORE df0dt = comlev1_thsice_1, key=ikey_1 #endif Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf) & + rhoi*cpIce *hi*Tice(2)) & /(6. _d 0*dt*k32 + rhoi*cpIce *hi) #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1010) & 'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice #endif C Compute final flux values at surfaces. fct = k12*(Tsf-Tice(1)) flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi flx0 = flx0 + df0dT*dTsf IF ( useBlkFlx ) THEN C-- needs to update also Evap (Tsf changes) since Latent heat has been updated evpAtm(i,j) = evap + dEvdT*dTsf ELSE C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics) evpAtm(i,j) = 0. ENDIF C- energy flux to Atmos: use net short-wave flux at surf. and C use latent heat = Lvap (energy=0 for liq. water at 0.oC) flxAtm(i,j) = netSW + flxExceptSw & + df0dT*dTsf + evpAtm(i,j)*Lfresh C- excess of energy @ surface (used for surface melting): sHeat(i,j) = flx0 - fct C- SW flux at sea-ice base left to the ocean flxSW(i,j) = fswocn #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) WRITE(6,1020) & 'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=', & flx0,fct,flx0-fct,flxCnB(i,j) #endif C Compute new enthalpy for each layer. qicen(1) = -cpWater*Tmlt1 + cpIce *(Tmlt1-Tice(1)) & + Lfresh*(1. _d 0-Tmlt1/Tice(1)) qicen(2) = -cpIce *Tice(2) + Lfresh C Make sure internal ice temperatures do not exceed Tmlt. C (This should not happen for reasonable values of i0swFrac) IF (Tice(1) .GE. Tmlt1) THEN WRITE (6,'(A,2I4,2I3,1P2E14.6)') & ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1 ENDIF IF (Tice(2) .GE. 0. _d 0) THEN WRITE (6,'(A,2I4,2I3,1P2E14.6)') & ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2) ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Update Sea-Ice state : tSrf(i,j) = Tsf tIc1(i,j) = Tice(1) tic2(i,j) = Tice(2) qIc1(i,j) = qicen(1) qIc2(i,j) = qicen(2) c dTsrf(i,j) = dTsf IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf c sHeat(i,j) = sHeating c flxCnB(i,j)= flxCnB c flxAtm(i,j)= flxAtm c evpAtm(i,j)= evpAtm #ifdef ALLOW_DBUG_THSICE IF ( dBug(i,j,bi,bj) ) THEN WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=', & Tsf, Tice, dTsf WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =', & sHeat(i,j), flxCnB(i,j), qicen WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=', & flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j) ENDIF #endif ELSE IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0 ENDIF ENDDO ENDDO #endif /* ALLOW_THSICE */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| RETURN END