SUBROUTINE ice_th_diff(nlay_s,nlay_i,kideb,kiut) !!------------------------------------------------------------------ !! *** ROUTINE ice_th_diff *** !! ** Purpose : !! This routine determines the time evolution of snow and sea-ice !! temperature profiles. !! ** Method : !! This is done by solving the heat equation diffusion with !! a Neumann boundary condition at the surface and a Dirichlet one !! at the bottom. Solar radiation is partially absorbed into the ice. !! The specific heat and thermal conductivities depend on ice salinity !! and temperature to take into account brine pocket melting. The !! numerical !! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid !! in the ice and snow system. !! The successive steps of this routine are !! Vertical grid !! 1. Thermal conductivity at the interfaces of the ice layers !! 2. Internal absorbed radiation !! 3. Scale factors due to non-uniform grid !! 4. Kappa factors !! Then iterative procedure begins !! 5. specific heat in the ice !! 6. eta factors !! 7. surface flux computation !! 8. tridiagonal system terms !! 9. solving the tridiagonal system with Gauss elimination !! Iterative procedure ends according to a criterion on evolution !! of temperature !! !! ** Arguments : !! kideb , kiut : Starting and ending points on which the !! the computation is applied !! !! ** Inputs / Ouputs : (global commons) !! surface temperature : t_su_b !! ice/snow temperatures : t_i_b, t_s_b !! ice salinities : s_i_b !! number of layers in the ice/snow: nlay_i, nlay_s !! total ice/snow thickness : ht_i_b, ht_s_b !! !! ** External : !! !! ** References : !! !! ** History : !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium !! USE lib_fortran INCLUDE 'type.com' INCLUDE 'para.com' INCLUDE 'const.com' INCLUDE 'ice.com' INCLUDE 'thermo.com' ! Local variables INTEGER numeqmin, numeqmax, numeq DIMENSION ztcond_i(0:nlay_i), & zkappa_s(0:nlay_s), & zkappa_i(0:nlay_i),ztstemp(0:nlay_s),ztitemp(0:nlay_i), & zspeche_i(0:nlay_i),ztsold(0:nlay_s),ztiold(0:nlay_i), & zeta_s(0:nlay_s),zeta_i(0:nlay_i),ztrid(2*maxnlay+1,3), & zindterm(2*maxnlay+1),zindtbis(2*maxnlay+1), & zdiagbis(2*maxnlay+1),ykn(nbpt),ykg(nbpt), & zrchu1(nbpt),zrchu2(nbpt),zqsat(nbpt) LOGICAL :: ln_write ! Local constants zeps = 1.0d-20 zg1s = 2.0 zg1 = 2.0 zbeta = 0.117 ! factor for thermal conductivity zerrmax = 1.0d-11 ! max error at the surface ! new lines zerrmax = 1.0e-4 nconv_max = 50 ln_write = .TRUE. DO 5 ji = kideb, kiut IF ( ln_write ) THEN WRITE(numout,*) ' ** ice_th_diff : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' WRITE(numout,*) ' nlay_i: ', nlay_i WRITE(numout,*) ' nlay_s: ', nlay_s WRITE(numout,*) ' t_su_b: ', t_su_b(ji) WRITE(numout,*) ' t_s_b : ', ( t_s_b(ji,layer), & layer = 1, nlay_s ) WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), & layer = 1, nlay_i ) WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji) WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) ENDIF ! Switches ! isnow equals 1 if snow is present and 0 if absent isnow = int(1.0-max(0.0,sign(1.0d0,-ht_s_b(ji)))) ! imelt equals 1 if surface is melting and 0 if not imelt = int(max(0.0,sign(1.0d0,t_su_b(ji)-tpw))) tfs(ji) = real(isnow)*tfsn+ (1.0-real(isnow))*tfsg ! Oceanic heat flux and precipitations fbbqb(ji) = oce_flx IF ( ln_write ) THEN WRITE(numout,*) ' isnow : ', isnow WRITE(numout,*) ' imelt : ', imelt WRITE(numout,*) ' tfs : ', tfs(ji) ENDIF 5 CONTINUE ! !------------------------------------------------------------------------------ ! 1) Thermal conductivity at the ice interfaces !------------------------------------------------------------------------------ ! ! Pringle et al., JGR 2007 formula ! 2.11 + 0.09 S/T - 0.011.T DO 10 ji = kideb, kiut ! thermal conductivity in the snow ykn(ji) = xkn zkimin = 0.1 ztcond_i(0) = xkg + betak1*s_i_b(ji,1) & / MIN( -zeps , t_i_b(ji,1) - tpw ) & - betak2* ( t_i_b(ji,1) - tpw ) ztcond_i(0) = MAX( ztcond_i(0) , zkimin ) DO layer = 1, nlay_i-1 ztcond_i(layer) = xkg + betak1*( s_i_b(ji,layer) & + s_i_b(ji,layer+1) ) / MIN(-zeps, & t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*tpw) & - betak2 * 0.5 * ( t_i_b(ji,layer) + ! bugfix fred dupont add 0.5 & t_i_b(ji,layer+1) - 2.0*tpw ) ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin ) END DO ztcond_i(nlay_i) = xkg + betak1*s_i_b(ji,nlay_i) / & MIN( -zeps , t_bo_b(ji) - tpw ) & - betak2 * ( t_bo_b(ji) - tpw ) ztcond_i(nlay_i) = MAX( ztcond_i(nlay_i) , zkimin ) IF ( ln_write ) THEN WRITE(numout,*) ' ztcond_i : ', ztcond_i(0:nlay_i) ENDIF 10 CONTINUE ! !------------------------------------------------------------------------------ ! 3) kappa factors !------------------------------------------------------------------------------ ! DO 30 ji = kideb, kiut ! snow zkappa_s(0) = ykn(ji)/max(zeps,deltaz_s_phy(1)) do layer = 1, nlay_s-1 zkappa_s(layer) = 2.0*ykn(ji)/ & max(zeps,deltaz_s_phy(layer)+deltaz_s_phy(layer+1)) end do zkappa_s(nlay_s) = ykn(ji)/max(zeps,deltaz_s_phy(nlay_s)) ! ice zkappa_i(0) = ztcond_i(0)/max(zeps,deltaz_i_phy(1)) do layer = 1, nlay_i-1 zkappa_i(layer) = 2.0*ztcond_i(layer)/ & max(zeps,deltaz_i_phy(layer)+deltaz_i_phy(layer+1)) end do zkappa_i(nlay_i) = ztcond_i(nlay_i) / & MAX(zeps,deltaz_i_phy(nlay_i)) ! interface zkappa_s(nlay_s) = 2.0*ykn(ji)*ztcond_i(0)/max(zeps, & (ztcond_i(0)*deltaz_s_phy(nlay_s) + & ykn(ji)*deltaz_i_phy(1))) zkappa_i(0) = zkappa_s(nlay_s)*real(isnow) & + zkappa_i(0)*(1.0-real(isnow)) IF ( ln_write ) THEN WRITE(numout,*) ' nlay_s : ', nlay_s WRITE(numout,*) ' zkappa_s : ', zkappa_s(0:nlay_s) WRITE(numout,*) ' zkappa_i : ', zkappa_i(0:nlay_i) ENDIF 30 CONTINUE ! !------------------------------------------------------------------------------| ! 4) iterative procedure begins | !------------------------------------------------------------------------------| ! DO 40 ji = kideb, kiut !------------------------------ ! keeping old values in memory !------------------------------ ztsuold = t_su_b(ji) t_su_b(ji) = min(t_su_b(ji),tfs(ji)-0.00001) DO layer = 1, nlay_s ztsold(layer) = t_s_b(ji,layer) END DO DO layer = 1, nlay_i ztiold(layer) = t_i_b(ji,layer) ti_old(layer) = ztiold(layer) END DO nconv = 0 40 CONTINUE !------------------------------ ! Beginning of the loop !------------------------------ zerrit = 10000.0 ! DO WHILE ( zerrit .GT. zerrmax ) DO WHILE ( ( zerrit .GT. zerrmax ) .AND. ( nconv < nconv_max ) ) do 42 ji = kideb, kiut !45 CONTINUE nconv = nconv+1 ztsutemp = t_su_b(ji) DO layer = 1, nlay_s ztstemp(layer) = t_s_b(ji,layer) END DO DO layer = 1, nlay_i ztitemp(layer) = t_i_b(ji,layer) END DO IF ( ln_write ) THEN WRITE(numout,*) ' zerrit : ', zerrit WRITE(numout,*) ' nconv : ', nconv WRITE(numout,*) ' t_s_b : ', t_s_b(ji,1) WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), layer = 1, & nlay_i ) ENDIF ! !------------------------------------------------------------------------------| ! 5) specific heat in the ice | !------------------------------------------------------------------------------| ! DO layer = 1, nlay_i zspeche_i(layer) = cpg + lfus*tmut*s_i_b(ji,layer)/ & MAX((t_i_b(ji,layer)-tpw)*(ztiold(layer)-tpw),zeps) END DO ! !------------------------------------------------------------------------------| ! 6) eta factors | !------------------------------------------------------------------------------| ! DO layer = 1, nlay_s zeta_s(layer) = ddtb / max(rhon*cpg*deltaz_s_phy(layer),zeps) END DO DO layer = 1, nlay_i zeta_i(layer) = ddtb / max(rhog*deltaz_i_phy(layer) * & zspeche_i(layer) & ,zeps) END DO ! !------------------------------------------------------------------------------| ! 7) surface flux computation | !------------------------------------------------------------------------------| ! t_su_b(ji) = min(t_su_b(ji),tfs(ji)) imelt = int(max(0.0,sign(1.0d0,t_su_b(ji)-tpw))) ! pressure of water vapor saturation (Pa) es = 611.0*10.0**(9.5*(t_su_b(ji)-273.16)/ & (t_su_b(ji)-7.66)) ! net longwave radiative flux ! MV BUG ! fratsb(ji) = emig*(ratbqb(ji)- ! & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)) fratsb(ji) = ratbqb(ji) - emig * & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji) ! sensible and latent heat flux CALL flx(ht_i_b(ji),t_su_b(ji),tabqb(ji),qabqb(ji),fcsb(ji), & fleb(ji),qsfcb(ji),zrchu1(ji),zrchu2(ji),vabqb(ji),zref) fcsb(ji) = -fcsb(ji) fleb(ji) = MIN( -fleb(ji) , 0.0 ) ! always negative, as precip ! energy already added ! qsfcb(ji) = zqsat(ji) ! intermediate variable zssdqw = qsfcb(ji)*qsfcb(ji)*psbqb(ji)/ & (0.622*es)*alog(10.0)*9.5* & ((273.16-7.66)/(t_su_b(ji)-7.66)**2.0) ! derivative of the surface atmospheric net flux dzf = -4.0*emig*stefan*t_su_b(ji) & *t_su_b(ji)*t_su_b(ji) & -(zrchu1(ji)+zrchu2(ji)*zssdqw) ! surface atmospheric net flux zf = ab(ji)*fsolgb(ji)+fratsb(ji) & +fcsb(ji)+fleb(ji) ! !------------------------------------------------------------------------------| ! 8) tridiagonal system terms | !------------------------------------------------------------------------------| ! ! layer denotes the number of the layer in the snow or in the ice ! numeq denotes the reference number of the equation in the tridiagonal ! system, terms of tridiagonal system are indexed as following : ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one ! ice interior terms (top equation has the same form as the others) DO numeq = 1, maxnlay ztrid(numeq,1) = 0.0 ztrid(numeq,2) = 0.0 ztrid(numeq,3) = 0.0 zindterm(numeq) = 0.0 zindtbis(numeq) = 0.0 zdiagbis(numeq) = 0.0 END DO DO numeq = nlay_s + 2, nlay_s + nlay_i layer = numeq - nlay_s - 1 ztrid(numeq,1) = - zeta_i(layer)*zkappa_i(layer-1) ztrid(numeq,2) = 1.0 + zeta_i(layer)*(zkappa_i(layer-1) + & zkappa_i(layer)) ztrid(numeq,3) = - zeta_i(layer)*zkappa_i(layer) zindterm(numeq) = ztiold(layer) + zeta_i(layer)* ! & ( radab_phy_i(layer) + radab_alg_i(layer) ) & radab_i(layer) END DO ! ice bottom terms numeq = nlay_s + nlay_i + 1 ztrid(numeq,1) = - zeta_i(nlay_i)*zkappa_i(nlay_i-1) ztrid(numeq,2) = 1.0 + zeta_i(nlay_i)*( zkappa_i(nlay_i)*zg1 & + zkappa_i(nlay_i-1) ) ztrid(numeq,3) = 0.0 zindterm(numeq) = ztiold(nlay_i) + zeta_i(nlay_i)* ! & ( radab_phy_i(nlay_i) + radab_alg_i(nlay_i) & ( radab_i(nlay_i) & + zkappa_i(nlay_i)*zg1 & *t_bo_b(ji) ) IF (ht_s_b(ji).GT.0.0) THEN ! !------------------------------------------------------------------------------| ! snow-covered cells | !------------------------------------------------------------------------------| ! ! snow interior terms (bottom equation has the same form as the others) do numeq = 3, nlay_s + 1 layer = numeq - 1 ztrid(numeq,1) = - zeta_s(layer)*zkappa_s(layer-1) ztrid(numeq,2) = 1.0 + zeta_s(layer)*( zkappa_s(layer-1) + & zkappa_s(layer) ) ztrid(numeq,3) = - zeta_s(layer)*zkappa_s(layer) zindterm(numeq) = ztsold(layer) + zeta_s(layer)* & radab_s(layer) end do ! case of only one layer in the ice (ice equation is altered) if (nlay_i.eq.1) then ztrid(nlay_s+2,3) = 0.0 zindterm(nlay_s+2) = zindterm(nlay_s+2) + zkappa_i(1)* & t_bo_b(ji) endif IF (t_su_b(ji).LT.tpw) THEN ! !------------------------------------------------------------------------------| ! case 1 : no surface melting - snow present | !------------------------------------------------------------------------------| ! zdifcase = 1.0 numeqmin = 1 numeqmax = nlay_i + nlay_s + 1 ! surface equation ztrid(1,1) = 0.0 ztrid(1,2) = dzf - zg1s*zkappa_s(0) ztrid(1,3) = zg1s*zkappa_s(0) zindterm(1) = dzf*t_su_b(ji) - zf ! first layer of snow equation ztrid(2,1) = - zkappa_s(0)*zg1s*zeta_s(1) ztrid(2,2) = 1.0 + zeta_s(1)*(zkappa_s(1) + zkappa_s(0)*zg1s) ztrid(2,3) = - zeta_s(1)* zkappa_s(1) zindterm(2) = ztsold(1) + zeta_s(1)*radab_s(1) else ! !------------------------------------------------------------------------------| ! case 2 : surface is melting - snow present | !------------------------------------------------------------------------------| ! zdifcase = 2.0 numeqmin = 2 numeqmax = nlay_i + nlay_s + 1 ! first layer of snow equation ztrid(2,1) = 0.0 ztrid(2,2) = 1.0 + zeta_s(1)*(zkappa_s(1) + zkappa_s(0)*zg1s) ztrid(2,3) = - zeta_s(1)*zkappa_s(1) zindterm(2) = ztsold(1) + zeta_s(1)*(radab_s(1) + zkappa_s(0)* & zg1s*t_su_b(ji)) ENDIF ELSE ! !------------------------------------------------------------------------------| ! cells without snow | !------------------------------------------------------------------------------| ! IF (t_su_b(ji).lt.tpw) THEN ! !------------------------------------------------------------------------------| ! case 3 : no surface melting - no snow | !------------------------------------------------------------------------------| ! zdifcase = 3.0 numeqmin = nlay_s + 1 numeqmax = nlay_i + nlay_s + 1 ! surface equation ztrid(numeqmin,1) = 0.0 ztrid(numeqmin,2) = dzf - zkappa_i(0)*zg1 ztrid(numeqmin,3) = zkappa_i(0)*zg1 zindterm(numeqmin) = dzf*t_su_b(ji) - zf ! first layer of ice equation ztrid(numeqmin+1,1) = - zkappa_i(0)*zg1*zeta_i(1) ztrid(numeqmin+1,2) = 1.0 + zeta_i(1)*(zkappa_i(1) + zkappa_i(0)* & zg1) ztrid(numeqmin+1,3) = - zeta_i(1)*zkappa_i(1) zindterm(numeqmin+1)= ztiold(1) + zeta_i(1)*radab_i(1) ! ( radab_phy_i(1) + ! & radab_alg_i(1) ) ! case of only one layer in the ice (surface & ice equations are altered) if (nlay_i.eq.1) then ztrid(numeqmin,1) = 0.0 ztrid(numeqmin,2) = dzf - zkappa_i(0)*2.0 ztrid(numeqmin,3) = zkappa_i(0)*2.0 ztrid(numeqmin+1,1) = -zkappa_i(0)*2.0*zeta_i(1) ztrid(numeqmin+1,2) = 1.0 + zeta_i(1)*(zkappa_i(0)*2.0 + & zkappa_i(1)) ztrid(numeqmin+1,3) = 0.0 zindterm(numeqmin+1) = ztiold(1) + zeta_i(1)* ! & ( radab_phy_i(1) + radab_alg_i(1) + & ( radab_i(1) + & zkappa_i(1)*t_bo_b(ji) ) endif else ! !------------------------------------------------------------------------------| ! case 4 : surface is melting - no snow | !------------------------------------------------------------------------------| ! zdifcase = 4.0 numeqmin = nlay_s + 2 numeqmax = nlay_i + nlay_s + 1 ! first layer of ice equation ztrid(numeqmin,1) = 0.0 ztrid(numeqmin,2) = 1.0 + zeta_i(1)*(zkappa_i(1) + zkappa_i(0)* & zg1) ztrid(numeqmin,3) = - zeta_i(1)* zkappa_i(1) zindterm(numeqmin) = ztiold(1) + zeta_i(1)* ( radab_i(1) + !(radab_phy_i(1) + ! & radab_alg_i(1) + & zkappa_i(0)*zg1*t_su_b(ji) ) ! case of only one layer in the ice (surface & ice equations are altered) if (nlay_i.eq.1) then ztrid(numeqmin,1) = 0.0 ztrid(numeqmin,2) = 1.0 + zeta_i(1)*(zkappa_i(0)*2.0 + & zkappa_i(1)) ztrid(numeqmin,3) = 0.0 zindterm(numeqmin) = ztiold(1) + zeta_i(1)* ! & (radab_phy_i(1) + radab_alg_i(1) & ( radab_i(1) + & + zkappa_i(1)*t_bo_b(ji) ) & + t_su_b(ji)*zeta_i(1)*zkappa_i(0)*2.0 endif endif endif ! !------------------------------------------------------------------------------| ! 9) tridiagonal system solving | !------------------------------------------------------------------------------| ! ! solving the tridiagonal system with Gauss elimination ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, ! McGraw-Hill 1984. ! computing temporary results zindtbis(numeqmin) = zindterm(numeqmin) zdiagbis(numeqmin) = ztrid(numeqmin,2) do numeq = numeqmin+1, numeqmax zdiagbis(numeq) = ztrid(numeq,2) - ztrid(numeq,1)* & ztrid(numeq-1,3)/zdiagbis(numeq-1) zindtbis(numeq) = zindterm(numeq) - ztrid(numeq,1)* & zindtbis(numeq-1)/zdiagbis(numeq-1) end do ! ice temperatures t_i_b(ji,nlay_i) = zindtbis(numeqmax)/zdiagbis(numeqmax) ! do numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 do numeq = nlay_i + nlay_s, nlay_s + 2, -1 layer = numeq - nlay_s - 1 t_i_b(ji,layer) = (zindtbis(numeq) - ztrid(numeq,3)* & t_i_b(ji,layer+1))/zdiagbis(numeq) end do ! snow temperatures if (ht_s_b(ji).gt.0.0) then t_s_b(ji,nlay_s) = (zindtbis(nlay_s+1) - ztrid(nlay_s+1,3)* & t_i_b(ji,1))/zdiagbis(nlay_s+1) if (nlay_s.gt.1) then do numeq = nlay_s, 2, -1 layer = numeq - 1 t_s_b(ji,layer) = (zindtbis(numeq) - ztrid(numeq,3)* & t_s_b(ji,layer+1))/zdiagbis(numeq) end do endif endif ! surface temperature if (t_su_b(ji).lt.tfs(ji)) then t_su_b(ji) = ( zindtbis(numeqmin) - ztrid(numeqmin,3)* & ( real(isnow)*t_s_b(ji,1) + & (1.0-real(isnow))*t_i_b(ji,1) ) ) / & zdiagbis(numeqmin) endif ! !-------------------------------------------------------------------------- ! 10) Has the scheme converged ?, end of the iterative procedure | !-------------------------------------------------------------------------- ! ! we verify that nothing has started to melt t_su_b(ji) = MIN( t_su_b(ji) , tfs(ji) ) DO layer = 1, nlay_i ztmelt_i = MIN( -tmut*s_i_b(ji,layer) + tpw , & 273.149999999d0) t_i_b(ji,layer) = MIN( t_i_b(ji,layer) ,ztmelt_i ) END DO ! zerrit is a residual which has to be under zerrmax zerrit = ABS(t_su_b(ji)-ztsutemp) DO layer = 1, nlay_s zerrit = MAX(zerrit,ABS(t_s_b(ji,layer) & - ztstemp(layer))) END DO DO layer = 1, nlay_i zerrit = max(zerrit,abs(t_i_b(ji,layer) - ztitemp(layer))) END DO 42 CONTINUE END DO DO 45 ji = kideb, kiut ! compute available energy for internal melt f_s_im(ji) = 0. DO layer = 1, nlay_s IF ( t_s_b(ji,layer) .GE. tfs(ji) ) THEN f_s_im(ji) = - rhon * ( cpg*( tfs(ji) - t_s_b(ji,layer))) & * ht_s_b(ji) / ddtb t_s_b(ji,layer) = MIN( t_s_b(ji,layer) ,tfs(ji) ) ENDIF END DO 45 CONTINUE ! !-------------------------------------------------------------------------- ! 11) Heat conduction fluxes | !-------------------------------------------------------------------------- ! DO 50 ji = kideb, kiut ! surface conduction flux fc_su(ji) = - real(isnow)*zkappa_s(0)*zg1s*(t_s_b(ji,1) - & t_su_b(ji)) & - (1.0-real(isnow))*zkappa_i(0)*zg1* & (t_i_b(ji,1) - t_su_b(ji)) ! bottom conduction flux fc_bo_i(ji) = - zkappa_i(nlay_i)* & ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) ! internal conduction fluxes : snow !--upper snow value fc_s(ji,0) = - isnow* & zkappa_s(0) * zg1s * ( t_s_b(ji,1) - & t_su_b(ji) ) !--basal snow value fc_s(ji,1) = - isnow* & zkappa_s(1) * ( t_i_b(ji,1) - & t_s_b(ji,1) ) ! internal conduction fluxes : ice !--upper layer fc_i(ji,0) = - isnow * ! interface flux if there is snow & ( zkappa_i(0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & - ( 1.0 - isnow ) * ( zkappa_i(0) * & zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if no !--internal ice layers DO layer = 1, nlay_i - 1 fc_i(ji,layer) = - zkappa_i(layer) * ( t_i_b(ji,layer+1) - & t_i_b(ji,layer) ) END DO !--under the basal ice layer fc_i(ji,nlay_i) = fc_bo_i(ji) ! case of only one layer in the ice IF (nlay_i.EQ.1) THEN fc_su(ji) = -real(isnow)*(zkappa_s(0)*(zg1s*(t_s_b(ji,1)- & t_su_b(ji)))) & -(1.0-real(isnow))*zkappa_i(0)*zg1*(t_i_b(ji,1)- & t_su_b(ji)) fc_bo_i(ji) = -zkappa_i(nlay_i)*zg1* & (t_bo_b(ji)-t_i_b(ji,nlay_i)) ENDIF ! !-------------------------------------------------------------------------- ! 12) Update atmospheric heat fluxes and energy of melting | !-------------------------------------------------------------------------- ! ! pressure of water vapor saturation (Pa) es = 611.0*10.0**(9.5*(t_su_b(ji)-273.16)/ & (t_su_b(ji)-7.66)) ! net longwave radiative flux ! MV BUG ! fratsb(ji) = emig*(ratbqb(ji)- ! & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)) fratsb(ji) = ratbqb(ji) - emig * & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji) ! sensible and latent heat flux CALL flx(ht_i_b(ji),t_su_b(ji),tabqb(ji),qabqb(ji),fcsb(ji), & fleb(ji),qsfcb(ji),zrchu1(ji),zrchu2(ji),vabqb(ji),zref) fcsb(ji) = -fcsb(ji) fleb(ji) = MIN( -fleb(ji) , 0.0 ) ! always negative, as precip ! energy already added ! ice energy of melting CALL ice_th_enmelt(kideb, kiut, nlay_s, nlay_i) IF ( ln_write ) THEN WRITE(numout,*) ' nconv : ', nconv WRITE(numout,*) ' zerrit : ', zerrit WRITE(numout,*) ' t_su_b: ', t_su_b(ji) WRITE(numout,*) ' t_s_b : ', ( t_s_b(ji,layer), & layer = 1, nlay_s ) WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), & layer = 1, nlay_i ) WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji) WRITE(numout,*) ENDIF 50 CONTINUE ! !------------------------------------------------------------------------------ ! End of ice_th_diff END SUBROUTINE