C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/thsice/thsice_solve4temp.F,v 1.2 2004/04/18 22:06:14 jmc Exp $ C $Name: $ #include "THSICE_OPTIONS.h" CBOP C !ROUTINE: THSICE_SOLVE4TEMP C !INTERFACE: SUBROUTINE THSICE_SOLVE4TEMP( I useBlkFlx, flxExcSw, Tf, hi, hs, U flxSW, Tsf, qicen, O Tice, sHeating, flxCnB, O dTsf, flxAtm, evpAtm, I i,j,bi,bj, myThid) C !DESCRIPTION: \bv C *==========================================================* C | S/R THSICE_SOLVE4TEMP C *==========================================================* C | Solve (implicitly) for sea-ice and surface temperature C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables === #include "THSICE_SIZE.h" #include "THSICE_PARAMS.h" C !INPUT/OUTPUT PARAMETERS: C == Routine Arguments == C useBlkFlx :: use surf. fluxes from bulk-forcing external S/R C flxExcSw :: 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 Tf :: freezing temperature (oC) of local sea-water C hi :: ice height C hs :: snow height C flxSW :: net Short-Wave flux (+=down) [W/m2]: input= at surface C :: output= at the sea-ice base to the ocean C Tsf :: surface (ice or snow) temperature C qicen :: ice enthalpy (J m-3) C Tice :: internal ice temperatures C sHeating :: surf heating left to melt snow or ice (= Atmos-conduction) C flxCnB :: heat flux conducted through the ice to bottom surface C dTsf :: surf. temp adjusment: Ts^n+1 - Ts^n 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 i,j,bi,bj :: indices of current grid point C myThid :: Thread no. that called this routine. LOGICAL useBlkFlx _RL flxExcSw(0:2) _RL Tf _RL hi _RL hs _RL flxSW _RL Tsf _RL qicen(nlyr) _RL Tice (nlyr) _RL sHeating _RL flxCnB _RL dTsf _RL flxAtm _RL evpAtm INTEGER i,j, bi,bj INTEGER myThid CEOP #ifdef ALLOW_THSICE 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 == Local Variables == INTEGER k _RL frsnow ! fractional snow cover _RL fswpen ! SW penetrating beneath surface (W m-2) _RL fswdn ! SW absorbed at surface (W m-2) _RL fswint ! SW absorbed in ice (W m-2) _RL fswocn ! SW passed through ice to ocean (W m-2) _RL 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] _RL evap, dEvdT _RL flx0 ! net surf heat flux, from Atmos. to sea-ice (W m-2) _RL fct ! heat conducted to top surface _RL df0dT ! deriv of flx0 wrt Tsf (W m-2 deg-1) _RL k12, k32 ! thermal conductivity terms _RL a10, b10 ! coefficients in quadratic eqn for T1 _RL a1, b1, c1 ! coefficients in quadratic eqn for T1 c _RL Tsf_start ! old value of Tsf _RL dt ! timestep INTEGER iceornot LOGICAL dBug 1010 FORMAT(A,I3,3F8.3) 1020 FORMAT(A,1P4E11.3) dt = thSIce_deltaT dBug = .FALSE. c dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 ) c dBug = ( bi.EQ.6 .AND. i.EQ.18 .AND. j.EQ.10 ) IF ( dBug ) WRITE(6,'(A,2I4,2I2)') 'ThSI_THERM: i,j=',i,j,bi,bj if ( hi.lt.himin ) then C If hi < himin, melt the ice. STOP 'THSICE_THERM: should not enter if hi 0 = ',Tice(1) write (6,*) 'BBerr Tice(2) > 0 = ',Tice(2) endif IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,Tice=',0,Tsf,Tice C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 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 Compute new surface and internal temperatures; iterate until C Tsfc converges. C ----- begin iteration ----- do 100 k = 1,nitMaxTsf C Save temperatures at start of iteration. c Tsf_start = Tsf IF ( useBlkFlx ) THEN C Compute top surface flux. if (hs.gt.3. _d -1) then iceornot=2 else iceornot=1 endif CALL THSICE_GET_BULKF( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) flx0 = fswdn + flxExceptSw ELSE flx0 = fswdn + flxExcSw(1) df0dT = flxExcSw(2) 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. 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 IF(dBug) WRITE(6,1010) 'ThSI_THERM: k,tice=', & k,Tsf,Tice(1),dTsf a1 = a10 + k12 b1 = b10 ! note b1 = b10 - k12*Tf0 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 CALL THSICE_GET_BULKF( I iceornot, Tsf, O flxExceptSw, df0dT, evap, dEvdT, I i,j,bi,bj,myThid ) flx0 = fswdn + flxExceptSw dTsf = 0. _d 0 ELSE flx0 = fswdn + flxExcSw(0) dTsf = 1000. df0dT = 0. ENDIF Tsf = 0. _d 0 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. IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,tice=', & k,Tsf,Tice(1),dTsf IF ( useBlkFlx ) THEN if (abs(dTsf).lt.Terrmax) then goto 150 elseif (k.eq.nitMaxTsf) then write (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf write (6,*) 'BB: thermw conv err, iceheight ', hi write (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0 if (Tsf.lt.-70. _d 0) stop !QQQQ fails endif ELSE goto 150 ENDIF 100 continue ! surface temperature iteration 150 continue C ------ end iteration ------------ c if (Tsf.lt.-70. _d 0) then cQQ print*,'QQQQQ Tsf =',Tsf cQQ stop !QQQQ fails c endif C Compute new bottom layer temperature. 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) IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,Tice=',k,Tsf,Tice C Compute final flux values at surfaces. fct = k12*(Tsf-Tice(1)) flxCnB = 4. _d 0*kice *(Tice(2)-Tf)/hi flx0 = flx0 + df0dT*dTsf IF ( useBlkFlx ) THEN evpAtm = evap C-- needs to update also Evap (Tsf changes) since Latent heat has been updated evpAtm = evap + dEvdT*dTsf 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 = flxSW + flxExceptSw + df0dT*dTsf & + evpAtm*Lfresh ELSE flxAtm = 0. evpAtm = 0. ENDIF sHeating = flx0 - fct C- SW flux at sea-ice base left to the ocean flxSW = fswocn IF (dBug) WRITE(6,1020) 'ThSI_THERM: flx0,fct,Dif,flxCnB=', & flx0,fct,flx0-fct,flxCnB 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 i0.) 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 #endif /* ALLOW_THSICE */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| RETURN END