MODULE limthd_dif #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !!====================================================================== !! *** MODULE limthd_dif *** !! heat diffusion in sea ice !! computation of surface and inner T !!====================================================================== !!---------------------------------------------------------------------- !! * Modules used USE par_oce ! ocean parameters USE phycst ! physical constants (ocean directory) USE ice_oce ! ice variables USE thd_ice USE iceini USE limistate USE in_out_manager USE ice USE par_ice USE lib_mpp IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_thd_dif ! called by lim_thd !! * Module variables REAL(wp) :: & ! constant values epsi20 = 1e-20 , & epsi13 = 1e-13 , & zzero = 0.e0 , & zone = 1.e0 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) !! (c) UCL-ASTR and Martin Vancoppenolle !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_dif( kideb , kiut , jl ) !!------------------------------------------------------------------ !! *** ROUTINE lim_thd_dif *** !! ** 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 !! 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 !! profile of the ice/snow layers : z_i, z_s !! total ice/snow thickness : ht_i_b, ht_s_b !! !! ** External : !! !! ** References : !! !! ** History : !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium !! (06-2005) Martin Vancoppenolle, 3d version !! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR) !! (04-2007) Energy conservation tested by M. Vancoppenolle !! !!------------------------------------------------------------------ !! * Arguments INTEGER , INTENT (in) :: & kideb , & ! Start point on which the the computation is applied kiut , & ! End point on which the the computation is applied jl ! Category number !! * Local variables INTEGER :: ji, & ! spatial loop index zji, zjj, & ! temporary dummy loop index numeq, & ! current reference number of equation layer, & ! vertical dummy loop index nconv, & ! number of iterations in iterative procedure minnumeqmin, & ! maxnumeqmax INTEGER , DIMENSION(jpij) :: & numeqmin, & ! reference number of top equation numeqmax, & ! reference number of bottom equation isnow ! switch for presence (1) or absence (0) of snow !! * New local variables REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & ztcond_i, & !Ice thermal conductivity zradtr_i, & !Radiation transmitted through the ice zradab_i, & !Radiation absorbed in the ice zkappa_i !Kappa factor in the ice REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & zradtr_s, & !Radiation transmited through the snow zradab_s, & !Radiation absorbed in the snow zkappa_s !Kappa factor in the snow REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & ztiold, & !Old temperature in the ice zeta_i, & !Eta factor in the ice ztitemp, & !Temporary temperature in the ice to check the convergence zspeche_i, & !Ice specific heat z_i !Vertical cotes of the layers in the ice REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & zeta_s, & !Eta factor in the snow ztstemp, & !Temporary temperature in the snow to check the convergence ztsold, & !Temporary temperature in the snow z_s !Vertical cotes of the layers in the snow REAL(wp) , DIMENSION(jpij,jkmax+2) :: & zindterm, & ! Independent term zindtbis, & ! temporary independent term zdiagbis REAL(wp) , DIMENSION(jpij,jkmax+2,3) :: & ztrid ! tridiagonal system terms REAL(wp), DIMENSION(jpij) :: & ztfs , & ! ice melting point ztsuold , & ! old surface temperature (before the iterative ! procedure ) ztsuoldit, & ! surface temperature at previous iteration zh_i , & !ice layer thickness zh_s , & !snow layer thickness zfsw , & !solar radiation absorbed at the surface zf , & ! surface flux function dzf ! derivative of the surface flux function REAL(wp) :: & ! constant values zeps = 1.0e-10, & ! zg1s = 2.0, & !: for the tridiagonal system zg1 = 2.0, & zgamma = 18009.0, & !: for specific heat zbeta = 0.117, & !: for thermal conductivity (could be 0.13) zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow zkimin = 0.10 , & !: minimum ice thermal conductivity zht_smin = 1.0e-4 !: minimum snow depth REAL(wp) :: & ! local variables ztmelt_i, & ! ice melting temperature zerritmax ! current maximal error on temperature REAL(wp), DIMENSION(jpij) :: & zerrit, & ! current error on temperature zdifcase, & ! case of the equation resolution (1->4) zftrice, & ! solar radiation transmitted through the ice zihic, zhsu !!-- End of declarations !!---------------------------------------------------------------------------------------------- IF(lwp) WRITE(numout,*)'lim_thd_dif : Heat diffusion in sea ice for cat :', jl ! !------------------------------------------------------------------------------! ! 1) Initialization ! !------------------------------------------------------------------------------! ! DO ji = kideb , kiut ! is there snow or not isnow(ji)= INT ( 1.0 - MAX( 0.0 , SIGN (1.0, - ht_s_b(ji) ) ) ) ! surface temperature of fusion ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt ! layer thickness zh_i(ji) = ht_i_b(ji) / nlay_i zh_s(ji) = ht_s_b(ji) / nlay_s END DO !-------------------- ! Ice / snow layers !-------------------- z_s(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st snow layer z_i(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st ice layer DO layer = 1, nlay_s DO ji = kideb , kiut ! vert. coord of the up. lim. of the layer-th snow layer z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s END DO END DO DO layer = 1, nlay_i DO ji = kideb , kiut ! vert. coord of the up. lim. of the layer-th ice layer z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i END DO END DO ! !------------------------------------------------------------------------------| ! 2) Radiations | !------------------------------------------------------------------------------| ! !------------------- ! Computation of i0 !------------------- ! i0 describes the fraction of solar radiation which does not contribute ! to the surface energy budget but rather penetrates inside the ice. ! We assume that no radiation is transmitted through the snow ! If there is no no snow ! zfsw = (1-i0).qsr_ice is absorbed at the surface ! zftrice = io.qsr_ice is below the surface ! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice DO ji = kideb , kiut ! switches isnow(ji) = INT ( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) ! hs > 0, isnow = 1 zhsu(ji) = hnzst !threshold for the computation of i0 zihic(ji) = MAX( zzero , 1.0 - ( ht_i_b(ji) / zhsu(ji) ) ) i0(ji) = ( 1.0 - isnow(ji) ) * & ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) !fr1_i0_1d = i0 for a thin ice surface !fr1_i0_2d = i0 for a thick ice surface ! a function of the cloud cover ! !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) !formula used in Cice END DO !------------------------------------------------------- ! Solar radiation absorbed / transmitted at the surface ! Derivative of the non solar flux !------------------------------------------------------- DO ji = kideb , kiut ! Shortwave radiation absorbed at surface zfsw(ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Solar radiation transmitted below the surface layer zftrice(ji)= qsr_ice_1d(ji) * i0(ji) ! derivative of incoming nonsolar flux dzf(ji) = dqns_ice_1d(ji) END DO !--------------------------------------------------------- ! Transmission - absorption of solar radiation in the ice !--------------------------------------------------------- DO ji = kideb , kiut ! Initialization zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow END DO ! Radiation through snow DO layer = 1, nlay_s DO ji = kideb , kiut ! radiation transmitted below the layer-th snow layer zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP ( - zraext_s * ( MAX ( 0.0 , & z_s(ji,layer) ) ) ) ! radiation absorbed by the layer-th snow layer zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer) END DO END DO ! Radiation through ice DO ji = kideb , kiut zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + & zftrice(ji) * ( 1 - isnow(ji) ) END DO DO layer = 1, nlay_i DO ji = kideb , kiut ! radiation transmitted below the layer-th ice layer zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP ( - kappa_i * ( MAX ( 0.0 , & z_i(ji,layer) ) ) ) ! radiation absorbed by the layer-th ice layer zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer) END DO END DO ! Radiation transmitted below the ice DO ji = kideb , kiut fstbif_1d(ji) = fstbif_1d(ji) + & zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) END DO ! +++++ ! just to check energy conservation DO ji = kideb , kiut zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 fstroc(zji,zjj,jl) = & zradtr_i(ji,nlay_i) END DO ! +++++ DO layer = 1, nlay_i DO ji = kideb , kiut radab(ji,layer) = zradab_i(ji,layer) END DO END DO ! !------------------------------------------------------------------------------| ! 3) Iterative procedure begins | !------------------------------------------------------------------------------| ! ! Old surface temperature DO ji = kideb, kiut ztsuold(ji) = t_su_b(ji) ! temperature at the beg of iter pr. ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter t_su_b(ji) = MIN(t_su_b(ji),ztfs(ji)-0.00001) !necessary zerrit(ji) = 1000.0 ! initial value of error END DO !RB Min global ?? ! Old snow temperature DO layer = 1, nlay_s DO ji = kideb , kiut ztsold(ji,layer) = t_s_b(ji,layer) END DO END DO ! Old ice temperature DO layer = 1, nlay_i DO ji = kideb , kiut ztiold(ji,layer) = t_i_b(ji,layer) END DO END DO nconv = 0 ! number of iterations zerritmax = 1000.0 ! maximal value of error on all points DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd)) nconv = nconv+1 ! !------------------------------------------------------------------------------| ! 4) Sea ice thermal conductivity | !------------------------------------------------------------------------------| ! IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner (1964) formula DO ji = kideb , kiut ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / & MIN(-zeps,t_i_b(ji,1)-rtt) ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) END DO ENDIF IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle et al formula included, ! 2.11 + 0.09 S/T - 0.011.T DO ji = kideb , kiut ztcond_i(ji,0) = rcdic + 0.09*s_i_b(ji,1) / & MIN(-zeps,t_i_b(ji,1)-rtt) - & 0.011* ( t_i_b(ji,1) - rtt ) ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) END DO ENDIF IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner DO layer = 1, nlay_i-1 DO ji = kideb , kiut ztcond_i(ji,layer) = rcdic + zbeta*( 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*rtt) ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) END DO END DO ENDIF IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle DO layer = 1, nlay_i-1 DO ji = kideb , kiut ztcond_i(ji,layer) = rcdic + 0.09*( 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*rtt) - & 0.011* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) END DO END DO ENDIF IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner DO ji = kideb , kiut ztcond_i(ji,nlay_i) = rcdic + zbeta*s_i_b(ji,nlay_i) / & MIN(-zeps,t_bo_b(ji)-rtt) ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) END DO ENDIF IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle DO ji = kideb , kiut ztcond_i(ji,nlay_i) = rcdic + 0.09*s_i_b(ji,nlay_i) / & MIN(-zeps,t_bo_b(ji)-rtt) - & 0.011* ( t_bo_b(ji) - rtt ) ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) END DO ENDIF ! !------------------------------------------------------------------------------| ! 5) kappa factors | !------------------------------------------------------------------------------| ! DO ji = kideb, kiut !-- Snow kappa factors zkappa_s(ji,0) = rcdsn / MAX(zeps,zh_s(ji)) zkappa_s(ji,nlay_s) = rcdsn / MAX(zeps,zh_s(ji)) END DO DO layer = 1, nlay_s-1 DO ji = kideb , kiut zkappa_s(ji,layer) = 2.0 * rcdsn / & MAX(zeps,2.0*zh_s(ji)) END DO END DO DO layer = 1, nlay_i-1 DO ji = kideb , kiut !-- Ice kappa factors zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ & MAX(zeps,2.0*zh_i(ji)) END DO END DO DO ji = kideb , kiut zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) !-- Interface zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*isnow(ji) & + zkappa_i(ji,0)*(1.0-isnow(ji)) END DO ! !------------------------------------------------------------------------------| ! 6) Sea ice specific heat, eta factors | !------------------------------------------------------------------------------| ! DO layer = 1, nlay_i DO ji = kideb , kiut ztitemp(ji,layer) = t_i_b(ji,layer) zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) zeta_i(ji,layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & zeps) END DO END DO DO layer = 1, nlay_s DO ji = kideb , kiut ztstemp(ji,layer) = t_s_b(ji,layer) zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) END DO END DO ! !------------------------------------------------------------------------------| ! 7) surface flux computation | !------------------------------------------------------------------------------| ! DO ji = kideb , kiut ! update of the non solar flux according to the update in T_su qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & ( t_su_b(ji) - ztsuoldit(ji) ) ! update incoming flux zf(ji) = zfsw(ji) & ! net absorbed solar radiation + qnsr_ice_1d(ji) ! non solar total flux ! (LWup, LWdw, SH, LH) END DO ! !------------------------------------------------------------------------------| ! 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,jkmax+2 DO ji = kideb , kiut ztrid(ji,numeq,1) = 0. ztrid(ji,numeq,2) = 0. ztrid(ji,numeq,3) = 0. zindterm(ji,numeq)= 0. zindtbis(ji,numeq)= 0. zdiagbis(ji,numeq)= 0. ENDDO ENDDO DO numeq = nlay_s + 2, nlay_s + nlay_i DO ji = kideb , kiut layer = numeq - nlay_s - 1 ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1) ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + & zkappa_i(ji,layer)) ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer) zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* & zradab_i(ji,layer) END DO ENDDO numeq = nlay_s + nlay_i + 1 DO ji = kideb , kiut !!ice bottom term ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & + zkappa_i(ji,nlay_i-1) ) ztrid(ji,numeq,3) = 0.0 zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & * t_bo_b(ji) ) ENDDO DO ji = kideb , kiut 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(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1) ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + & zkappa_s(ji,layer) ) ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer) zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* & zradab_s(ji,layer) END DO !!case of only one layer in the ice (ice equation is altered) IF ( nlay_i.eq.1 ) THEN ztrid(ji,nlay_s+2,3) = 0.0 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & t_bo_b(ji) ENDIF IF ( t_su_b(ji) .LT. rtt ) THEN !------------------------------------------------------------------------------| ! case 1 : no surface melting - snow present | !------------------------------------------------------------------------------| zdifcase(ji) = 1.0 numeqmin(ji) = 1 numeqmax(ji) = nlay_i + nlay_s + 1 !!surface equation ztrid(ji,1,1) = 0.0 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji) !!first layer of snow equation ztrid(ji,2,1) = - zkappa_s(ji,0)*zg1s*zeta_s(ji,1) ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) ELSE ! !------------------------------------------------------------------------------| ! case 2 : surface is melting - snow present | !------------------------------------------------------------------------------| ! zdifcase(ji) = 2.0 numeqmin(ji) = 2 numeqmax(ji) = nlay_i + nlay_s + 1 !!first layer of snow equation ztrid(ji,2,1) = 0.0 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & zkappa_s(ji,0) * zg1s ) ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & ( zradab_s(ji,1) + & zkappa_s(ji,0) * zg1s * t_su_b(ji) ) ENDIF ELSE ! !------------------------------------------------------------------------------| ! cells without snow | !------------------------------------------------------------------------------| ! IF (t_su_b(ji) .LT. rtt) THEN ! !------------------------------------------------------------------------------| ! case 3 : no surface melting - no snow | !------------------------------------------------------------------------------| ! zdifcase(ji) = 3.0 numeqmin(ji) = nlay_s + 1 numeqmax(ji) = nlay_i + nlay_s + 1 !!surface equation ztrid(ji,numeqmin(ji),1) = 0.0 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji) !!first layer of ice equation ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & + zkappa_i(ji,0) * zg1 ) ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) !!case of only one layer in the ice (surface & ice equations are altered) IF (nlay_i.eq.1) THEN ztrid(ji,numeqmin(ji),1) = 0.0 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)*2.0*zeta_i(ji,1) ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & zkappa_i(ji,1)) ztrid(ji,numeqmin(ji)+1,3) = 0.0 zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* & ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) ) ENDIF ELSE ! !------------------------------------------------------------------------------| ! case 4 : surface is melting - no snow | !------------------------------------------------------------------------------| ! zdifcase(ji) = 4.0 numeqmin(ji) = nlay_s + 2 numeqmax(ji) = nlay_i + nlay_s + 1 !!first layer of ice equation ztrid(ji,numeqmin(ji),1) = 0.0 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & zg1) ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & zkappa_i(ji,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(ji,numeqmin(ji),1) = 0.0 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & zkappa_i(ji,1)) ztrid(ji,numeqmin(ji),3) = 0.0 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* & (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) & + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 ENDIF ENDIF ENDIF END DO ! !------------------------------------------------------------------------------| ! 9) tridiagonal system solving | !------------------------------------------------------------------------------| ! ! Solve the tridiagonal system with Gauss elimination method. ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, ! McGraw-Hill 1984. maxnumeqmax = 0 minnumeqmin = jkmax+4 DO ji = kideb , kiut zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) minnumeqmin = MIN(numeqmin(ji),minnumeqmin) maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax) END DO DO layer = minnumeqmin+1, maxnumeqmax DO ji = kideb , kiut numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji)) zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) END DO END DO DO ji = kideb , kiut ! ice temperatures t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) END DO DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 DO ji = kideb , kiut layer = numeq - nlay_s - 1 t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & t_i_b(ji,layer+1))/zdiagbis(ji,numeq) END DO END DO DO ji = kideb , kiut ! snow temperatures IF (ht_s_b(ji).GT.0) & t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & * MAX(0.0,SIGN(1.0,ht_s_b(ji)-zeps)) ! surface temperature isnow(ji) = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) ztsuoldit(ji) = t_su_b(ji) IF (t_su_b(ji) .LT. ztfs(ji)) & t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* & ( isnow(ji)*t_s_b(ji,1) + & (1.0-isnow(ji))*t_i_b(ji,1) ) ) / & zdiagbis(ji,numeqmin(ji)) END DO ! !-------------------------------------------------------------------------- ! 10) Has the scheme converged ?, end of the iterative procedure | !-------------------------------------------------------------------------- ! ! check that nowhere it has started to melt ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd DO ji = kideb , kiut t_su_b(ji) = MAX(MIN(t_su_b(ji),ztfs(ji)),190.0) zerrit(ji) = ABS(t_su_b(ji)-ztsuoldit(ji)) END DO DO layer = 1, nlay_s DO ji = kideb , kiut zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 t_s_b(ji,layer) = MAX(MIN(t_s_b(ji,layer),rtt),190.0) zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) & - ztstemp(ji,layer))) END DO END DO DO layer = 1, nlay_i DO ji = kideb , kiut ztmelt_i = -tmut*s_i_b(ji,layer) +rtt t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0) zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) END DO END DO ! Compute spatial maximum over all errors ! note that this could be optimized substantially by iterating only ! the non-converging points zerritmax = 0.0 DO ji = kideb , kiut zerritmax = MAX(zerritmax,zerrit(ji)) END DO IF( lk_mpp ) CALL mpp_max(zerritmax, kcom=ncomm_ice) END DO ! End of the do while iterative procedure WRITE(numout,*) ' zerritmax : ', zerritmax WRITE(numout,*) ' nconv : ', nconv ! !-------------------------------------------------------------------------- ! 11) Fluxes at the interfaces | !-------------------------------------------------------------------------- ! DO ji = kideb, kiut ! update of latent heat fluxes qla_ice_1d (ji) = qla_ice_1d (ji) + & dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) ! surface ice conduction flux isnow(ji) = int(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) fc_su(ji) = - isnow(ji)*zkappa_s(ji,0)*zg1s*(t_s_b(ji,1) - & t_su_b(ji)) & - (1.0-isnow(ji))*zkappa_i(ji,0)*zg1* & (t_i_b(ji,1) - t_su_b(ji)) ! bottom ice conduction flux fc_bo_i(ji) = - zkappa_i(ji,nlay_i)* & ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) END DO !-------------------------! ! Heat conservation ! !-------------------------! IF ( con_i ) THEN DO ji = kideb, kiut ! Upper snow value fc_s(ji,0) = - isnow(ji)* & zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - & t_su_b(ji) ) ! Bott. snow value fc_s(ji,1) = - isnow(ji)* & zkappa_s(ji,1) * ( t_i_b(ji,1) - & t_s_b(ji,1) ) END DO ! Upper ice layer DO ji = kideb, kiut fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & - ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * & zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not END DO ! Internal ice layers DO layer = 1, nlay_i - 1 DO ji = kideb, kiut fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - & t_i_b(ji,layer) ) zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 END DO END DO ! Bottom ice layers DO ji = kideb, kiut fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i)* & ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) END DO ENDIF END SUBROUTINE lim_thd_dif #else !!====================================================================== !! *** MODULE limthd_dif *** !! no sea ice model !!====================================================================== CONTAINS SUBROUTINE lim_thd_dif ! Empty routine END SUBROUTINE lim_thd_dif #endif END MODULE limthd_dif