MODULE icethd_zdf !!====================================================================== !! *** MODULE icethd_zdf *** !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) !!====================================================================== !! History : LIM ! 02-2003 (M. Vancoppenolle) original 1D code !! ! 06-2005 (M. Vancoppenolle) 3d version !! ! 11-2006 (X Fettweis) Vectorization by Xavier !! ! 04-2007 (M. Vancoppenolle) Energy conservation !! 4.0 ! 2011-02 (G. Madec) dynamical allocation !! - ! 2012-05 (C. Rousset) add penetration solar flux !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' ESIM sea-ice model !!---------------------------------------------------------------------- USE dom_oce ! ocean space and time domain USE phycst ! physical constants (ocean directory) USE ice ! sea-ice: variables USE ice1D ! sea-ice: thermodynamics variables ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) IMPLICIT NONE PRIVATE PUBLIC ice_thd_zdf ! called by icethd PUBLIC ice_thd_zdf_init ! called by icestp !!** namelist (namthd_zdf) ** LOGICAL :: ln_zdf_BL99 ! Heat diffusion follows Bitz and Lipscomb (1999) LOGICAL :: ln_cndi_U64 ! thermal conductivity: Untersteiner (1964) LOGICAL :: ln_cndi_P07 ! thermal conductivity: Pringle et al (2007) REAL(wp) :: rn_kappa_i ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] REAL(wp), PUBLIC :: rn_cnd_s ! thermal conductivity of the snow [W/m/K] INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation ! ! associated indices: INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) INTEGER , PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns INTEGER , PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns INTEGER , PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_thd_zdf !! !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_zdf *** !! ** Purpose : !! This chooses between the appropriate routine for the !! computation of vertical diffusion !! !!------------------------------------------------------------------- !! SELECT CASE ( nice_zdf ) ! Choose the vertical heat diffusion solver !------------- CASE( np_BL99 ) ! BL99 solver !------------- IF ( nice_jules == np_jules_OFF ) THEN ! No Jules coupler CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) ENDIF IF ( nice_jules == np_jules_EMULE ) THEN ! Jules coupler is emulated CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) ENDIF IF ( nice_jules == np_jules_ACTIVE ) THEN ! Jules coupler is emulated CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) ENDIF END SELECT END SUBROUTINE ice_thd_zdf SUBROUTINE ice_thd_zdf_BL99(k_jules) !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_zdf_BL99 *** !! ** Purpose : !! This routine determines the time evolution of snow and sea-ice !! temperature profiles, using the original Bitz and Lipscomb (1999) algorithm !! !! ** 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. initialization of ice-snow layers thicknesses !! 2. Internal absorbed and transmitted radiation !! Then iterative procedure begins !! 3. Thermal conductivity !! 4. Kappa factors !! 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 !! 10. Fluxes at the interfaces !! !! ** Inputs / Ouputs : (global commons) !! surface temperature : t_su_1d !! ice/snow temperatures : t_i_1d, t_s_1d !! ice salinities : sz_i_1d !! number of layers in the ice/snow: nlay_i, nlay_s !! total ice/snow thickness : h_i_1d, h_s_1d !!------------------------------------------------------------------- INTEGER, INTENT(in) :: k_jules ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) INTEGER :: ji, jk ! spatial loop index INTEGER :: jm ! current reference number of equation INTEGER :: jm_mint, jm_maxt INTEGER :: iconv ! number of iterations in iterative procedure INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system REAL(wp) :: zg1 = 2._wp ! REAL(wp) :: zgamma = 18009._wp ! for specific heat REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature REAL(wp) :: ztmelt_i ! ice melting temperature REAL(wp) :: zdti_max ! current maximal error on temperature REAL(wp) :: zcpi ! Ice specific heat REAL(wp) :: zhfx_err, zdq ! diag errors on heat REAL(wp) :: zfac ! dummy factor REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat REAL(wp) :: zfr1, zfr2, zfrqsr_tr_i ! dummy factor ! Mono-category REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done REAL(wp) :: zhe ! dummy factor REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity !!------------------------------------------------------------------ ! --- diag error on heat diffusion - PART 1 --- ! DO ji = 1, npti zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) END DO !------------------ ! 1) Initialization !------------------ DO ji = 1, npti isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not ! layer thickness zh_i(ji) = h_i_1d(ji) * r1_nlay_i zh_s(ji) = h_s_1d(ji) * r1_nlay_s END DO ! WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) ELSEWHERE ; z1_h_i(1:npti) = 0._wp END WHERE WHERE( zh_s(1:npti) >= epsi10 ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) ELSEWHERE ; z1_h_s(1:npti) = 0._wp END WHERE ! Store initial temperatures and non solar heat fluxes IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 ztsuold(1:npti) = t_su_1d(1:npti) ! surface temperature initial value zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not ENDIF ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature !------------- ! 2) Radiation !------------- ! --- Transmission/absorption of solar radiation in the ice --- ! ! zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 ) ! standard value ! zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 ) ! zfr2 such that zfr1 + zfr2 to equal 1 ! DO ji = 1, npti ! zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) ) ! zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness ! IF ( h_s_1d(ji) >= 0.0_wp ) zfrqsr_tr_i = 0._wp ! snow fully opaque ! qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji) ! transmitted solar radiation ! zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) ! zftrice(ji) = qsr_ice_tr_1d(ji) ! i0(ji) = zfrqsr_tr_i ! END DO zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) DO jk = 1, nlay_s DO ji = 1, npti ! ! radiation transmitted below the layer-th snow layer zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) ! ! radiation absorbed by the layer-th snow layer zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) END DO END DO zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) ) DO jk = 1, nlay_i DO ji = 1, npti ! ! radiation transmitted below the layer-th ice layer zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) ! ! radiation absorbed by the layer-th ice layer zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) END DO END DO ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice ! iconv = 0 ! number of iterations zdti_max = 1000._wp ! maximal value of error on all points ! !============================! DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! ! !============================! iconv = iconv + 1 ! ztib(1:npti,:) = t_i_1d(1:npti,:) ztsb(1:npti,:) = t_s_1d(1:npti,:) ! !-------------------------------- ! 3) Sea ice thermal conductivity !-------------------------------- IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T ! DO ji = 1, npti ztcond_i(ji,0) = rcdic + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) END DO DO jk = 1, nlay_i-1 DO ji = 1, npti ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) END DO END DO ! ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T ! DO ji = 1, npti ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) END DO DO jk = 1, nlay_i-1 DO ji = 1, npti ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) & & - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) END DO END DO ! ENDIF ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) ! !--- G(he) : enhancement of thermal conductivity in mono-category case ! Computation of effective thermal conductivity G(h) ! Used in mono-category case only to simulate an ITD implicitly ! Fichefet and Morales Maqueda, JGR 1997 zghe(1:npti) = 1._wp ! SELECT CASE ( nn_monocat ) CASE ( 1 , 3 ) zepsilon = 0.1_wp DO ji = 1, npti zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) THEN zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he) ENDIF END DO END SELECT ! !----------------- ! 4) kappa factors !----------------- !--- Snow DO jk = 0, nlay_s-1 DO ji = 1, npti zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) END DO END DO DO ji = 1, npti ! Snow-ice interface zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) IF( zfac > epsi10 ) THEN zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac ELSE zkappa_s(ji,nlay_s) = 0._wp ENDIF END DO !--- Ice DO jk = 0, nlay_i DO ji = 1, npti zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) END DO END DO DO ji = 1, npti ! Snow-ice interface zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) END DO ! !-------------------------------------- ! 5) Sea ice specific heat, eta factors !-------------------------------------- DO jk = 1, nlay_i DO ji = 1, npti zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) END DO END DO DO jk = 1, nlay_s DO ji = 1, npti zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) END DO END DO ! !----------------------------------------! ! ! ! JULES IF (OFF or SND MODE) ! ! ! !----------------------------------------! ! IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode ! ! In OFF mode the original BL99 temperature computation is used ! (with qsr_ice, qns_ice and dqns_ice as inputs) ! ! In SND mode, the computation is required to compute the conduction fluxes ! !---------------------------- ! 6) surface flux computation !---------------------------- DO ji = 1, npti ! update of the non solar flux according to the update in T_su qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) END DO DO ji = 1, npti zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar END DO ! !---------------------------- ! 7) tridiagonal system terms !---------------------------- !!layer denotes the number of the layer in the snow or in the ice !!jm 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) ztrid (1:npti,:,:) = 0._wp zindterm(1:npti,:) = 0._wp zindtbis(1:npti,:) = 0._wp zdiagbis(1:npti,:) = 0._wp DO jm = nlay_s + 2, nlay_s + nlay_i DO ji = 1, npti jk = jm - nlay_s - 1 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) END DO ENDDO jm = nlay_s + nlay_i + 1 DO ji = 1, npti !!ice bottom term ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) ztrid(ji,jm,3) = 0.0 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) ENDDO DO ji = 1, npti ! !---------------------! IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! ! !---------------------! ! snow interior terms (bottom equation has the same form as the others) DO jm = 3, nlay_s + 1 jk = jm - 1 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) END DO ! case of only one layer in the ice (ice equation is altered) IF ( nlay_i == 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_1d(ji) ENDIF IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting jm_min(ji) = 1 jm_max(ji) = nlay_i + nlay_s + 1 ! surface equation ztrid(ji,1,1) = 0.0 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(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 ! jm_min(ji) = 2 jm_max(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_1d(ji) ) ENDIF ! !---------------------! ELSE ! cells without snow ! ! !---------------------! ! IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting ! jm_min(ji) = nlay_s + 1 jm_max(ji) = nlay_i + nlay_s + 1 ! surface equation ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0)*zg1 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) ! first layer of ice equation ztrid(ji,jm_min(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) ztrid(ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) zindterm(ji,jm_min(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 == 1 ) THEN ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) ztrid(ji,jm_min(ji)+1,3) = 0.0 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) ENDIF ELSE !-- case 2 : surface is melting jm_min(ji) = nlay_s + 2 jm_max(ji) = nlay_i + nlay_s + 1 ! first layer of ice equation ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) ! case of only one layer in the ice (surface & ice equations are altered) IF ( nlay_i == 1 ) THEN ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) ztrid(ji,jm_min(ji),3) = 0.0 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 ENDIF ENDIF ENDIF ! zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) ! END DO ! !------------------------------ ! 8) tridiagonal system solving !------------------------------ ! Solve the tridiagonal system with Gauss elimination method. ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 jm_maxt = 0 jm_mint = nlay_i+5 DO ji = 1, npti jm_mint = MIN(jm_min(ji),jm_mint) jm_maxt = MAX(jm_max(ji),jm_maxt) END DO DO jk = jm_mint+1, jm_maxt DO ji = 1, npti jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) END DO END DO DO ji = 1, npti ! ice temperatures t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) END DO DO jm = nlay_i + nlay_s, nlay_s + 2, -1 DO ji = 1, npti jk = jm - nlay_s - 1 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) END DO END DO DO ji = 1, npti ! snow temperatures IF( h_s_1d(ji) > 0._wp ) THEN t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & & / zdiagbis(ji,nlay_s+1) ENDIF ! surface temperature ztsub(ji) = t_su_1d(ji) IF( t_su_1d(ji) < rt0 ) THEN t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) ENDIF END DO ! !-------------------------------------------------------------- ! 9) Has the scheme converged ?, end of the iterative procedure !-------------------------------------------------------------- ! check that nowhere it has started to melt ! zdti_max is a measure of error, it has to be under zdti_bnd zdti_max = 0._wp DO ji = 1, npti t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) END DO DO jk = 1, nlay_s DO ji = 1, npti t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) END DO END DO DO jk = 1, nlay_i DO ji = 1, npti ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) END DO END DO ! Compute spatial maximum over all errors ! note that this could be optimized substantially by iterating only the non-converging points IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) ! !----------------------------------------! ! ! ! JULES IF (RCV MODE) ! ! ! !----------------------------------------! ! ELSE IF ( k_jules == np_zdf_jules_RCV ) THEN ! RCV mode ! ! In RCV mode, we use a modified BL99 solver ! with conduction flux (qcn_ice) as forcing term ! !---------------------------- ! 7) tridiagonal system terms !---------------------------- !!layer denotes the number of the layer in the snow or in the ice !!jm 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) ztrid (1:npti,:,:) = 0._wp zindterm(1:npti,:) = 0._wp zindtbis(1:npti,:) = 0._wp zdiagbis(1:npti,:) = 0._wp DO jm = nlay_s + 2, nlay_s + nlay_i DO ji = 1, npti jk = jm - nlay_s - 1 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) END DO ENDDO jm = nlay_s + nlay_i + 1 DO ji = 1, npti !!ice bottom term ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) ztrid(ji,jm,3) = 0.0 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) ENDDO DO ji = 1, npti ! !---------------------! IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! ! !---------------------! ! snow interior terms (bottom equation has the same form as the others) DO jm = 3, nlay_s + 1 jk = jm - 1 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) END DO ! case of only one layer in the ice (ice equation is altered) IF ( nlay_i == 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_1d(ji) ENDIF jm_min(ji) = 2 jm_max(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) 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) + qcn_ice_1d(ji) ) ! !---------------------! ELSE ! cells without snow ! ! !---------------------! jm_min(ji) = nlay_s + 2 jm_max(ji) = nlay_i + nlay_s + 1 ! first layer of ice equation ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & & ( zradab_i(ji,1) + qcn_ice_1d(ji) ) ! case of only one layer in the ice (surface & ice equations are altered) IF ( nlay_i == 1 ) THEN ztrid(ji,jm_min(ji),1) = 0.0 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) ztrid(ji,jm_min(ji),3) = 0.0 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & & + qcn_ice_1d(ji) ) ENDIF ENDIF ! zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) ! END DO ! !------------------------------ ! 8) tridiagonal system solving !------------------------------ ! Solve the tridiagonal system with Gauss elimination method. ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 jm_maxt = 0 jm_mint = nlay_i+5 DO ji = 1, npti jm_mint = MIN(jm_min(ji),jm_mint) jm_maxt = MAX(jm_max(ji),jm_maxt) END DO DO jk = jm_mint+1, jm_maxt DO ji = 1, npti jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) END DO END DO DO ji = 1, npti ! ice temperatures t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) END DO DO jm = nlay_i + nlay_s, nlay_s + 2, -1 DO ji = 1, npti jk = jm - nlay_s - 1 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) END DO END DO DO ji = 1, npti ! snow temperatures IF( h_s_1d(ji) > 0._wp ) THEN t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & & / zdiagbis(ji,nlay_s+1) ENDIF END DO ! !-------------------------------------------------------------- ! 9) Has the scheme converged ?, end of the iterative procedure !-------------------------------------------------------------- ! check that nowhere it has started to melt ! zdti_max is a measure of error, it has to be under zdti_bnd zdti_max = 0._wp DO jk = 1, nlay_s DO ji = 1, npti t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) END DO END DO DO jk = 1, nlay_i DO ji = 1, npti ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) END DO END DO ! Compute spatial maximum over all errors ! note that this could be optimized substantially by iterating only the non-converging points IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) ENDIF ! k_jules END DO ! End of the do while iterative procedure IF( ln_icectl .AND. lwp ) THEN WRITE(numout,*) ' zdti_max : ', zdti_max WRITE(numout,*) ' iconv : ', iconv ENDIF ! !----------------------------- ! 10) Fluxes at the interfaces !----------------------------- ! ! --- update conduction fluxes ! DO ji = 1, npti ! ! surface ice conduction flux fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) ! ! bottom ice conduction flux fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) END DO ! ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! ! DO ji = 1, npti IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) ELSE ! RCV mode hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) ENDIF END DO ! ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! ! IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF CALL ice_thd_enmelt ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation DO ji = 1, npti zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) ELSE ! case T_su = 0degC zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) ENDIF ELSE ! RCV CASE zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) ENDIF ! total heat sink to be sent to the ocean hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) END DO ! ! --- SIMIP diagnostics ! DO ji = 1, npti !--- Conduction fluxes (positive downwards) diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) !--- Snow-ice interfacial temperature (diagnostic SIMIP) zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac ELSE t_si_1d(ji) = t_su_1d(ji) ENDIF END DO ENDIF ! !--------------------------------------------------------------------------------------- ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes !--------------------------------------------------------------------------------------- ! effective conductivity and 1st layer temperature (Jules coupling) DO ji = 1, npti cnd_ice_1d(ji) = 2._wp * ( isnow(ji) * zkappa_s(ji,0) + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) ) t1_ice_1d (ji) = ( isnow(ji) * t_s_1d (ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d (ji,1) ) END DO ! IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode ! Restore temperatures to their initial values t_s_1d(1:npti,:) = ztsold (1:npti,:) t_i_1d(1:npti,:) = ztiold (1:npti,:) qcn_ice_1d(1:npti) = fc_su(1:npti) ENDIF END SUBROUTINE ice_thd_zdf_BL99 SUBROUTINE ice_thd_enmelt !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_enmelt *** !! !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature !! !! ** Method : Formula (Bitz and Lipscomb, 1999) !!------------------------------------------------------------------- INTEGER :: ji, jk ! dummy loop indices REAL(wp) :: ztmelts ! local scalar !!------------------------------------------------------------------- ! DO jk = 1, nlay_i ! Sea ice energy of melting DO ji = 1, npti ztmelts = - tmut * sz_i_1d(ji,jk) t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point ! (sometimes dif scheme produces abnormally high temperatures) e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & & + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & & - rcp * ztmelts ) END DO END DO DO jk = 1, nlay_s ! Snow energy of melting DO ji = 1, npti e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) END DO END DO ! END SUBROUTINE ice_thd_enmelt SUBROUTINE ice_thd_zdf_init !!----------------------------------------------------------------------- !! *** ROUTINE ice_thd_zdf_init *** !! !! ** Purpose : Physical constants and parameters associated with !! ice thermodynamics !! !! ** Method : Read the namthd_zdf namelist and check the parameters !! called at the first timestep (nit000) !! !! ** input : Namelist namthd_zdf !!------------------------------------------------------------------- INTEGER :: ios, ioptio ! Local integer output status for namelist read !! NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i !!------------------------------------------------------------------- ! REWIND( numnam_ice_ref ) ! Namelist namthd_zdf in reference namelist : Ice thermodynamics READ ( numnam_ice_ref, namthd_zdf, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in reference namelist', lwp ) REWIND( numnam_ice_cfg ) ! Namelist namthd_zdf in configuration namelist : Ice thermodynamics READ ( numnam_ice_cfg, namthd_zdf, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in configuration namelist', lwp ) IF(lwm) WRITE ( numoni, namthd_zdf ) ! ! IF(lwp) THEN ! control print WRITE(numout,*) 'ice_thd_zdf_init: Ice vertical heat diffusion' WRITE(numout,*) '~~~~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namthd_zdf:' WRITE(numout,*) ' Bitz and Lipscomb (1999) formulation ln_zdf_BL99 = ', ln_zdf_BL99 WRITE(numout,*) ' thermal conductivity in the ice (Untersteiner 1964) ln_cndi_U64 = ', ln_cndi_U64 WRITE(numout,*) ' thermal conductivity in the ice (Pringle et al 2007) ln_cndi_P07 = ', ln_cndi_P07 WRITE(numout,*) ' thermal conductivity in the snow rn_cnd_s = ', rn_cnd_s WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i ENDIF ! IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) ENDIF ! !== set the choice of ice vertical thermodynamic formulation ==! ioptio = 0 ! !--- BL99 thermo dynamics (linear liquidus + constant thermal properties) IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99 ; ENDIF IF( ioptio /= 1 ) CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) ! END SUBROUTINE ice_thd_zdf_init #else !!---------------------------------------------------------------------- !! Default option Dummy Module No ESIM sea-ice model !!--------------------------------------------------------------------- #endif !!====================================================================== END MODULE icethd_zdf