New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 8984 – NEMO

Changeset 8984


Ignore:
Timestamp:
2017-12-12T00:32:40+01:00 (5 years ago)
Author:
clem
Message:

clarify vertical heat diffusion. There are still issues with heat conservation in case nice_jules=2 but it will be solved in due time. Sette tests are all OK except ISOMIP

Location:
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8933 r8984  
    154154   !!---------------------------------------------------------------------- 
    155155   !                                     !!** ice-generic parameters namelist (nampar) ** 
    156    INTEGER           , PUBLIC ::   jpl             !: number of ice  categories  
    157    INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers  
    158    INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers  
    159    INTEGER           , PUBLIC ::   nn_monocat      !: virtual ITD mono-category parameterizations (1-4) or not (0) 
    160    LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F) 
    161    LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F) 
    162    REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere 
    163    REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere 
    164    CHARACTER(len=256), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
    165    CHARACTER(len=256), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
    166    CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory 
    167    CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory 
     156   INTEGER           , PUBLIC ::   jpl              !: number of ice  categories  
     157   INTEGER           , PUBLIC ::   nlay_i           !: number of ice  layers  
     158   INTEGER           , PUBLIC ::   nlay_s           !: number of snow layers  
     159   INTEGER           , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1-4) or not (0) 
     160   LOGICAL           , PUBLIC ::   ln_icedyn        !: flag for ice dynamics (T) or not (F) 
     161   LOGICAL           , PUBLIC ::   ln_icethd        !: flag for ice thermo   (T) or not (F) 
     162   REAL(wp)          , PUBLIC ::   rn_amax_n        !: maximum ice concentration Northern hemisphere 
     163   REAL(wp)          , PUBLIC ::   rn_amax_s        !: maximum ice concentration Southern hemisphere 
     164   CHARACTER(len=256), PUBLIC ::   cn_icerst_in     !: suffix of ice restart name (input) 
     165   CHARACTER(len=256), PUBLIC ::   cn_icerst_out    !: suffix of ice restart name (output) 
     166   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir  !: ice restart input directory 
     167   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir !: ice restart output directory 
    168168 
    169169   !                                     !!** ice-itd namelist (namitd) ** 
    … …  
    192192   !                                      !   =-1  Do nothing (needs N(cat) fluxes) 
    193193   !                                      !   = 0  Average N(cat) fluxes then apply the average over the N(cat) ice  
    194    !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice 
    195    !                                      !                                   using T-ice and albedo sensitivity 
     194   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity 
    196195   !                                      !   = 2  Redistribute a single flux over categories 
    197  
    198    INTEGER , PUBLIC            ::   nice_jules            !: Choice of jules coupling  
    199    !                                                      !  Associated indices 
    200    INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   !: no Jules coupling (ice thermodynamics forced via qsr and qns) 
    201    INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 
    202    INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 
     196   INTEGER , PUBLIC            ::   nice_jules           !: Choice of jules coupling  
     197   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0  !: no Jules coupling (ice thermodynamics forced via qsr and qns) 
     198   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1  !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 
     199   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2  !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 
     200 
     201   !                                     !!** ice-vertical diffusion namelist (namthd_zdf) ** 
     202   LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964) 
     203   LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007) 
     204   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
     205   REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]    
    203206 
    204207   !                                     !!** ice-salinity namelist (namthd_sal) ** 
    … …  
    211214   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU] 
    212215 
    213    !                                     !!** namelist namthd_pnd 
     216   !                                     !!** ice-ponds namelist (namthd_pnd) 
    214217   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
    215218   LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater 
    … …  
    255258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange   [kg.m-2.s-1] 
    256259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sni !: snow ice growth component of wfx_snw [kg.m-2.s-1] 
    257    ! MV MP 2016 
    258260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1] 
    259261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1] 
    260    ! END MV MP 2016 
    261262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice  [kg.m-2.s-1] 
    262263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice    [kg.m-2.s-1] 
    … …  
    353354   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sz_i     !: ice salinity          [PSU] 
    354355 
    355    ! MV MP 2016 
    356356   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond fraction per grid cell area 
    357357   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area [m] 
    … …  
    361361   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction 
    362362   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area [m] 
    363    ! END MV MP 2016 
    364363 
    365364   !!---------------------------------------------------------------------- 
    … …  
    376375   !! * Ice thickness distribution variables 
    377376   !!---------------------------------------------------------------------- 
    378    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max         !: Boundary of ice thickness categories in thickness space 
    379    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean        !: Mean ice thickness in catgories  
     377   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_max         !: Boundary of ice thickness categories in thickness space 
     378   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_mean        !: Mean ice thickness in catgories  
    380379   ! 
    381380   !!---------------------------------------------------------------------- 
    … …  
    385384   ! trp   ''         ''     ''       advection (transport of ice) 
    386385   ! 
    387    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vi   !: transport of ice volume 
    388    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vs   !: transport of snw volume 
    389    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2) 
    390    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2) 
    391    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_sv   !: transport of salt content 
    392    ! 
    393    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat     !: snw/ice heat content variation   [W/m2]  
    394    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sice     !: ice salt content variation   []  
    395    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s]  
    396    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s]  
     386   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume 
     387   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume 
     388   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_ei   !: transport of ice enthalpy (W/m2) 
     389   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_es   !: transport of snw enthalpy (W/m2) 
     390   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_sv   !: transport of salt content 
     391   ! 
     392   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_heat     !: snw/ice heat content variation   [W/m2]  
     393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_sice     !: ice salt content variation   []  
     394   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice     !: ice volume variation   [m/s]  
     395   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s]  
    397396 
    398397   ! 
    … …  
    401400   !!---------------------------------------------------------------------- 
    402401   ! Extra sea ice diagnostics to address the data request 
    403    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_si          !: Temperature at Snow-ice interface (K)  
    404    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K)  
    405    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2) 
    406    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2) 
     402   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_si          !: Temperature at Snow-ice interface (K)  
     403   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tm_si         !: mean temperature at the snow-ice interface (K)  
     404   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_fc_bo    !: Bottom conduction flux (W/m2) 
     405   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_fc_su    !: Surface conduction flux (W/m2) 
    407406 
    408407   ! 

  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90

    r8882 r8984  
    251251            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 
    252252         END DO 
    253          ! MV MP 2016 
    254253         IF ( ln_pnd_H12 ) THEN 
    255             pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:) 
    256             pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:) 
     254            pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) 
     255            pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) 
    257256         ENDIF 
    258          ! END MV MP 2016 
    259257      END DO 
    260258      ! 

  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8939 r8984  
    1717   !!  ice_thd_zdf      : select the appropriate routine for vertical heat diffusion calculation 
    1818   !!  ice_thd_zdf_BL99 :  
    19    !!  ice_thd_enmelt   : 
    2019   !!  ice_thd_zdf_init : 
    2120   !!---------------------------------------------------------------------- 
    22    USE dom_oce        ! ocean space and time domain 
    23    USE phycst         ! physical constants (ocean directory)  
    24    USE ice            ! sea-ice: variables 
    25    USE ice1D          ! sea-ice: thermodynamics variables 
     21   USE dom_oce         ! ocean space and time domain 
     22   USE phycst          ! physical constants (ocean directory)  
     23   USE ice             ! sea-ice: variables 
     24   USE icethd_zdf_BL99 ! sea-ice: vertical diffusion (Bitz and Lipscomb, 1999)  
    2625   ! 
    27    USE in_out_manager ! I/O manager 
    28    USE lib_mpp        ! MPP library 
    29    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
     26   USE in_out_manager  ! I/O manager 
     27   USE lib_mpp         ! MPP library 
     28   USE lib_fortran     ! fortran utilities (glob_sum + no signed zero) 
    3029 
    3130   IMPLICIT NONE 
    … …  
    3534   PUBLIC   ice_thd_zdf_init   ! called by icestp 
    3635 
     36   INTEGER ::   nice_zdf       ! Choice of the type of vertical heat diffusion formulation 
     37   !                                 ! associated indices: 
     38   INTEGER, PARAMETER ::   np_BL99 = 1   ! Bitz and Lipscomb (1999) 
     39!! INTEGER, PARAMETER ::   np_XXXX = 2 
     40 
    3741   !!** namelist (namthd_zdf) ** 
    38    LOGICAL  ::   ln_zdf_BL99      ! Heat diffusion follows Bitz and Lipscomb (1999) 
    39    LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964) 
    40    LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007) 
    41    REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
    42    REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K] 
     42   LOGICAL ::   ln_zdf_BL99    ! Heat diffusion follows Bitz and Lipscomb (1999) 
    4343 
    44    INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation 
    45    !                                           ! associated indices: 
    46    INTEGER, PARAMETER ::   np_BL99 = 1         ! Bitz and Lipscomb (1999) 
    47  
    48    INTEGER, PARAMETER ::   np_zdf_jules_OFF = 0   !   compute all temperatures from qsr and qns 
    49    INTEGER, PARAMETER ::   np_zdf_jules_SND = 1   !   compute conductive heat flux and surface temperature from qsr and qns 
    50    INTEGER, PARAMETER ::   np_zdf_jules_RCV = 2   !   compute snow and inner ice temperatures from qcnd 
    51     
    5244   !!---------------------------------------------------------------------- 
    5345   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    … …  
    6759      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver 
    6860      ! 
    69       !                             !-------------       
    70       CASE( np_BL99 )               ! BL99 solver 
    71          !                          !------------- 
     61      !                             !-------------!       
     62      CASE( np_BL99 )               ! BL99 solver ! 
     63         !                          !-------------! 
    7264         SELECT CASE( nice_jules ) 
    73          CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )   ! No Jules coupler      
    74          CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )   ! Jules coupler is emulated (send/recieve) 
    75                                        CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
    76          CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )   ! Jules coupler is active (receive only) 
     65         !                         ! No Jules coupler ==> default option 
     66         CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_jules_OFF    ) 
     67         ! 
     68         !                         ! Jules coupler is emulated => 1st call to get the needed fields (conduction...) 
     69         !                                                        2nd call to use these fields to calculate heat diffusion    
     70         CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_jules_EMULE  ) 
     71                                       CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE ) 
     72         ! 
     73         !                         ! Jules coupler is active ==> Met Office default option 
     74         CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE ) 
     75         ! 
    7776         END SELECT 
     77         ! 
    7878      END SELECT 
    7979       
    … …  
    8181    
    8282    
    83    SUBROUTINE ice_thd_zdf_BL99( k_jules ) 
    84       !!------------------------------------------------------------------- 
    85       !!                ***  ROUTINE ice_thd_zdf_BL99  *** 
    86       !! 
    87       !! ** Purpose :   computes the time evolution of snow and sea-ice temperature 
    88       !!              profiles, using the original Bitz and Lipscomb (1999) algorithm 
    89       !! 
    90       !! ** Method  :   solves the heat equation diffusion with a Neumann boundary 
    91       !!              condition at the surface and a Dirichlet one at the bottom.  
    92       !!              Solar radiation is partially absorbed into the ice. 
    93       !!              The specific heat and thermal conductivities depend on ice  
    94       !!              salinity and temperature to take into account brine pocket    
    95       !!              melting. The numerical scheme is an iterative Crank-Nicolson 
    96       !!              on a non-uniform multilayer grid in the ice and snow system. 
    97       !! 
    98       !!           The successive steps of this routine are 
    99       !!           1.  initialization of ice-snow layers thicknesses 
    100       !!           2.  Internal absorbed and transmitted radiation 
    101       !!           Then iterative procedure begins 
    102       !!           3.  Thermal conductivity 
    103       !!           4.  Kappa factors 
    104       !!           5.  specific heat in the ice 
    105       !!           6.  eta factors 
    106       !!           7.  surface flux computation 
    107       !!           8.  tridiagonal system terms 
    108       !!           9.  solving the tridiagonal system with Gauss elimination 
    109       !!           Iterative procedure ends according to a criterion on evolution 
    110       !!           of temperature 
    111       !!           10. Fluxes at the interfaces 
    112       !! 
    113       !! ** Inputs / Ouputs : (global commons) 
    114       !!           surface temperature : t_su_1d 
    115       !!           ice/snow temperatures   : t_i_1d, t_s_1d 
    116       !!           ice salinities          : sz_i_1d 
    117       !!           number of layers in the ice/snow: nlay_i, nlay_s 
    118       !!           total ice/snow thickness : h_i_1d, h_s_1d 
    119       !!------------------------------------------------------------------- 
    120       INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 
    121       ! 
    122       INTEGER ::   ji, jk         ! spatial loop index 
    123       INTEGER ::   jm             ! current reference number of equation 
    124       INTEGER ::   jm_mint, jm_maxt 
    125       INTEGER ::   iconv          ! number of iterations in iterative procedure 
    126       INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure 
    127       ! 
    128       INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation 
    129       INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation 
    130       ! 
    131       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    132       REAL(wp) ::   zg1       =  2._wp        ! 
    133       REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat 
    134       REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13) 
    135       REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow 
    136       REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    137       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C  
    138       REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    139       REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    140       REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    141       REAL(wp) ::   zcpi                      ! Ice specific heat 
    142       REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
    143       REAL(wp) ::   zfac                      ! dummy factor 
    144       ! 
    145       REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow 
    146       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration 
    147       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness 
    148       REAL(wp), DIMENSION(jpij) ::   zh_s, z1_h_s ! snow layer thickness 
    149       REAL(wp), DIMENSION(jpij) ::   zqns_ice_b   ! solar radiation absorbed at the surface 
    150       REAL(wp), DIMENSION(jpij) ::   zfnet        ! surface flux function 
    151       REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function 
    152       ! 
    153       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    154       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    155       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    156       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
    157       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    158       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    159       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    160       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    161       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    162       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    163       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    164       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    165       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    166       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    167       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    168       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    169       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    170       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    171       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    172       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    173       ! 
    174       REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor 
    175       ! 
    176       ! Mono-category 
    177       REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done 
    178       REAL(wp) :: zhe           ! dummy factor 
    179       REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity 
    180       !!------------------------------------------------------------------      
    181  
    182       ! --- diag error on heat diffusion - PART 1 --- ! 
    183       DO ji = 1, npti 
    184          zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
    185             &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )  
    186       END DO 
    187  
    188       !------------------ 
    189       ! 1) Initialization 
    190       !------------------ 
    191       DO ji = 1, npti 
    192          isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not 
    193          ! layer thickness 
    194          zh_i(ji) = h_i_1d(ji) * r1_nlay_i 
    195          zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
    196       END DO 
    197       ! 
    198       WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 
    199       ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp 
    200       END WHERE 
    201       ! 
    202       WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
    203       ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp 
    204       END WHERE 
    205       ! 
    206       ! Store initial temperatures and non solar heat fluxes 
    207       IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
    208          ! 
    209          ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1 
    210          ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value 
    211          zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux  
    212          zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value 
    213          ! 
    214          t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not 
    215          ! 
    216       ENDIF    
    217       ! 
    218       ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature 
    219       ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature 
    220  
    221       !------------- 
    222       ! 2) Radiation 
    223       !------------- 
    224       ! --- Transmission/absorption of solar radiation in the ice --- ! 
    225 !     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value 
    226 !     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1 
    227 ! 
    228 !     DO ji = 1, npti 
    229 ! 
    230 !        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )      
    231 ! 
    232 !              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
    233 !              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque 
    234 ! 
    235 !              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation  
    236 ! 
    237 !              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 
    238 !              zftrice(ji) = qsr_ice_tr_1d(ji) 
    239 !              i0(ji) = zfrqsr_tr_i 
    240 ! 
    241 !     END DO 
    242  
    243       zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 
    244       DO jk = 1, nlay_s 
    245          DO ji = 1, npti 
    246             !                             ! radiation transmitted below the layer-th snow layer 
    247             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) 
    248             !                             ! radiation absorbed by the layer-th snow layer 
    249             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
    250          END DO 
    251       END DO 
    252       ! 
    253       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) ) 
    254       DO jk = 1, nlay_i  
    255          DO ji = 1, npti 
    256             !                             ! radiation transmitted below the layer-th ice layer 
    257             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 
    258             !                             ! radiation absorbed by the layer-th ice layer 
    259             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
    260          END DO 
    261       END DO 
    262       ! 
    263       ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice 
    264       ! 
    265       iconv    =  0          ! number of iterations 
    266       zdti_max =  1000._wp   ! maximal value of error on all points 
    267       ! 
    268       !                                                          !============================! 
    269       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    270          !                                                       !============================! 
    271          iconv = iconv + 1 
    272          ! 
    273          ztib(1:npti,:) = t_i_1d(1:npti,:) 
    274          ztsb(1:npti,:) = t_s_1d(1:npti,:) 
    275          ! 
    276          !-------------------------------- 
    277          ! 3) Sea ice thermal conductivity 
    278          !-------------------------------- 
    279          IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
    280             ! 
    281             DO ji = 1, npti 
    282                ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    283                ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    284             END DO 
    285             DO jk = 1, nlay_i-1 
    286                DO ji = 1, npti 
    287                   ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    288                      &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    289                END DO 
    290             END DO 
    291             ! 
    292          ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
    293             ! 
    294             DO ji = 1, npti 
    295                ztcond_i(ji,0)      = rcdic + 0.09_wp * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    296                   &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    297                ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    298                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    299             END DO 
    300             DO jk = 1, nlay_i-1 
    301                DO ji = 1, npti 
    302                   ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    303                      &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   & 
    304                      &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    305                END DO 
    306             END DO 
    307             ! 
    308          ENDIF 
    309          ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )         
    310          ! 
    311          !--- G(he) : enhancement of thermal conductivity in mono-category case 
    312          ! Computation of effective thermal conductivity G(h) 
    313          ! Used in mono-category case only to simulate an ITD implicitly 
    314          ! Fichefet and Morales Maqueda, JGR 1997 
    315          zghe(1:npti) = 1._wp 
    316          ! 
    317          SELECT CASE ( nn_monocat ) 
    318          ! 
    319          CASE ( 1 , 3 ) 
    320             ! 
    321             zepsilon = 0.1_wp 
    322             DO ji = 1, npti 
    323                zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity 
    324                zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )     ! Effective thickness he (zhe) 
    325                IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN 
    326                   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he) 
    327                ENDIF 
    328             END DO 
    329             ! 
    330          END SELECT 
    331          ! 
    332          !----------------- 
    333          ! 4) kappa factors 
    334          !----------------- 
    335          !--- Snow 
    336          DO jk = 0, nlay_s-1 
    337             DO ji = 1, npti 
    338                zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
    339             END DO 
    340          END DO 
    341          DO ji = 1, npti   ! Snow-ice interface 
    342             zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
    343             IF( zfac > epsi10 ) THEN 
    344                zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
    345             ELSE 
    346                zkappa_s(ji,nlay_s) = 0._wp 
    347             ENDIF 
    348          END DO 
    349  
    350          !--- Ice 
    351          DO jk = 0, nlay_i 
    352             DO ji = 1, npti 
    353                zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
    354             END DO 
    355          END DO 
    356          DO ji = 1, npti   ! Snow-ice interface 
    357             zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    358          END DO 
    359          ! 
    360          !-------------------------------------- 
    361          ! 5) Sea ice specific heat, eta factors 
    362          !-------------------------------------- 
    363          DO jk = 1, nlay_i 
    364             DO ji = 1, npti 
    365                zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
    366                zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi )  
    367             END DO 
    368          END DO 
    369  
    370          DO jk = 1, nlay_s 
    371             DO ji = 1, npti 
    372                zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
    373             END DO 
    374          END DO 
    375           
    376          ! 
    377          !----------------------------------------! 
    378          !                                        ! 
    379          !       JULES IF (OFF or SND MODE)       ! 
    380          !                                        ! 
    381          !----------------------------------------! 
    382          ! 
    383           
    384          IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
    385           
    386             ! 
    387             ! In OFF mode the original BL99 temperature computation is used 
    388             ! (with qsr_ice, qns_ice and dqns_ice as inputs) 
    389             ! 
    390             ! In SND mode, the computation is required to compute the conduction fluxes 
    391             ! 
    392           
    393             !---------------------------- 
    394             ! 6) surface flux computation 
    395             !---------------------------- 
    396  
    397             ! update of the non solar flux according to the update in T_su 
    398             DO ji = 1, npti 
    399                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    400             END DO 
    401  
    402             DO ji = 1, npti 
    403                zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 
    404             END DO 
    405             ! 
    406             !---------------------------- 
    407             ! 7) tridiagonal system terms 
    408             !---------------------------- 
    409             !!layer denotes the number of the layer in the snow or in the ice 
    410             !!jm denotes the reference number of the equation in the tridiagonal 
    411             !!system, terms of tridiagonal system are indexed as following : 
    412             !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
    413  
    414             !!ice interior terms (top equation has the same form as the others) 
    415             ztrid   (1:npti,:,:) = 0._wp 
    416             zindterm(1:npti,:)   = 0._wp 
    417             zindtbis(1:npti,:)   = 0._wp 
    418             zdiagbis(1:npti,:)   = 0._wp 
    419  
    420             DO jm = nlay_s + 2, nlay_s + nlay_i  
    421             DO ji = 1, npti 
    422                jk = jm - nlay_s - 1 
    423                ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
    424                ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    425                ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
    426                zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    427             END DO 
    428             END DO 
    429  
    430             jm =  nlay_s + nlay_i + 1 
    431             DO ji = 1, npti 
    432             !!ice bottom term 
    433             ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    434             ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
    435             ztrid(ji,jm,3)  = 0.0 
    436             zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    437                &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    438             END DO 
    439  
    440  
    441             DO ji = 1, npti 
    442             !                               !---------------------! 
    443             IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells ! 
    444                !                            !---------------------! 
    445                ! snow interior terms (bottom equation has the same form as the others) 
    446                DO jm = 3, nlay_s + 1 
    447                   jk = jm - 1 
    448                   ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
    449                   ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    450                   ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    451                   zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    452                END DO 
    453  
    454                ! case of only one layer in the ice (ice equation is altered) 
    455                IF ( nlay_i == 1 ) THEN 
    456                   ztrid(ji,nlay_s+2,3)  = 0.0 
    457                   zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    458                ENDIF 
    459  
    460                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    461  
    462                   jm_min(ji) = 1 
    463                   jm_max(ji) = nlay_i + nlay_s + 1 
    464  
    465                   ! surface equation 
    466                   ztrid(ji,1,1)  = 0.0 
    467                   ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    468                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    469                   zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
    470  
    471                   ! first layer of snow equation 
    472                   ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
    473                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    474                   ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
    475                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    476  
    477                ELSE                            !--  case 2 : surface is melting 
    478                   ! 
    479                   jm_min(ji) = 2 
    480                   jm_max(ji) = nlay_i + nlay_s + 1 
    481  
    482                   ! first layer of snow equation 
    483                   ztrid(ji,2,1)  = 0.0 
    484                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    485                   ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
    486                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    487                ENDIF 
    488                !                            !---------------------! 
    489             ELSE                            ! cells without snow  ! 
    490                !                            !---------------------! 
    491                ! 
    492                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    493                   ! 
    494                   jm_min(ji) = nlay_s + 1 
    495                   jm_max(ji) = nlay_i + nlay_s + 1 
    496  
    497                   ! surface equation    
    498                   ztrid(ji,jm_min(ji),1)  = 0.0 
    499                   ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    500                   ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
    501                   zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
    502  
    503                   ! first layer of ice equation 
    504                   ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    505                   ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    506                   ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
    507                   zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    508  
    509                   ! case of only one layer in the ice (surface & ice equations are altered) 
    510                   IF ( nlay_i == 1 ) THEN 
    511                      ztrid(ji,jm_min(ji),1)    = 0.0 
    512                      ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    513                      ztrid(ji,jm_min(ji),3)    =  zkappa_i(ji,0) * 2.0 
    514                      ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
    515                      ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    516                      ztrid(ji,jm_min(ji)+1,3)  = 0.0 
    517                      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) ) 
    518                   ENDIF 
    519  
    520                ELSE                            !--  case 2 : surface is melting 
    521  
    522                   jm_min(ji)    =  nlay_s + 2 
    523                   jm_max(ji)    =  nlay_i + nlay_s + 1 
    524  
    525                   ! first layer of ice equation 
    526                   ztrid(ji,jm_min(ji),1)  = 0.0 
    527                   ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    528                   ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
    529                   zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
    530                      &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    531  
    532                   ! case of only one layer in the ice (surface & ice equations are altered) 
    533                   IF ( nlay_i == 1 ) THEN 
    534                      ztrid(ji,jm_min(ji),1)  = 0.0 
    535                      ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    536                      ztrid(ji,jm_min(ji),3)  = 0.0 
    537                      zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
    538                         &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    539                   ENDIF 
    540  
    541                ENDIF 
    542             ENDIF 
    543             ! 
    544             zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
    545             zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
    546             ! 
    547             END DO 
    548             ! 
    549             !------------------------------ 
    550             ! 8) tridiagonal system solving 
    551             !------------------------------ 
    552             ! Solve the tridiagonal system with Gauss elimination method. 
    553             ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    554             jm_maxt = 0 
    555             jm_mint = nlay_i+5 
    556             DO ji = 1, npti 
    557                jm_mint = MIN(jm_min(ji),jm_mint) 
    558                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    559             END DO 
    560  
    561             DO jk = jm_mint+1, jm_maxt 
    562                DO ji = 1, npti 
    563                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
    564                   zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
    565                   zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    566                END DO 
    567             END DO 
    568  
    569             ! ice temperatures 
    570             DO ji = 1, npti 
    571                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    572             END DO 
    573  
    574             DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    575                DO ji = 1, npti 
    576                   jk = jm - nlay_s - 1 
    577                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    578                END DO 
    579             END DO 
    580  
    581             DO ji = 1, npti 
    582                ! snow temperatures       
    583                IF( h_s_1d(ji) > 0._wp ) THEN 
    584                   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) 
    585                ENDIF 
    586                ! surface temperature 
    587                ztsub(ji) = t_su_1d(ji) 
    588                IF( t_su_1d(ji) < rt0 ) THEN 
    589                   t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    590                      &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    591                ENDIF 
    592             END DO 
    593             ! 
    594             !-------------------------------------------------------------- 
    595             ! 9) Has the scheme converged ?, end of the iterative procedure 
    596             !-------------------------------------------------------------- 
    597             ! check that nowhere it has started to melt 
    598             ! zdti_max is a measure of error, it has to be under zdti_bnd 
    599             zdti_max = 0._wp 
    600             DO ji = 1, npti 
    601                t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    602                zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    603             END DO 
    604  
    605             DO jk = 1, nlay_s 
    606                DO ji = 1, npti 
    607                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    608                   zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    609                END DO 
    610             END DO 
    611  
    612             DO jk = 1, nlay_i 
    613                DO ji = 1, npti 
    614                   ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
    615                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
    616                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    617                END DO 
    618             END DO 
    619  
    620             ! Compute spatial maximum over all errors 
    621             ! note that this could be optimized substantially by iterating only the non-converging points 
    622             IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    623          ! 
    624          !----------------------------------------! 
    625          !                                        ! 
    626          !       JULES IF (RCV MODE)              ! 
    627          !                                        ! 
    628          !----------------------------------------! 
    629          ! 
    630  
    631          ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode 
    632           
    633             ! 
    634             ! In RCV mode, we use a modified BL99 solver  
    635             ! with conduction flux (qcn_ice) as forcing term 
    636             ! 
    637             !---------------------------- 
    638             ! 7) tridiagonal system terms 
    639             !---------------------------- 
    640             !!layer denotes the number of the layer in the snow or in the ice 
    641             !!jm denotes the reference number of the equation in the tridiagonal 
    642             !!system, terms of tridiagonal system are indexed as following : 
    643             !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
    644  
    645             !!ice interior terms (top equation has the same form as the others) 
    646             ztrid   (1:npti,:,:) = 0._wp 
    647             zindterm(1:npti,:)   = 0._wp 
    648             zindtbis(1:npti,:)   = 0._wp 
    649             zdiagbis(1:npti,:)   = 0._wp 
    650  
    651             DO jm = nlay_s + 2, nlay_s + nlay_i  
    652                DO ji = 1, npti 
    653                   jk = jm - nlay_s - 1 
    654                   ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
    655                   ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    656                   ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
    657                   zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    658                END DO 
    659             ENDDO 
    660  
    661             jm =  nlay_s + nlay_i + 1 
    662             DO ji = 1, npti 
    663                !!ice bottom term 
    664                ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    665                ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
    666                ztrid(ji,jm,3)  = 0.0 
    667                zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    668                   &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    669             ENDDO 
    670  
    671  
    672             DO ji = 1, npti 
    673                !                              !---------------------! 
    674                IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
    675                   !                           !---------------------! 
    676                   ! snow interior terms (bottom equation has the same form as the others) 
    677                   DO jm = 3, nlay_s + 1 
    678                      jk = jm - 1 
    679                      ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
    680                      ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    681                      ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    682                      zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    683                   END DO 
    684                    
    685                   ! case of only one layer in the ice (ice equation is altered) 
    686                   IF ( nlay_i == 1 ) THEN 
    687                      ztrid(ji,nlay_s+2,3)  = 0.0 
    688                      zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    689                   ENDIF 
    690                    
    691                   jm_min(ji) = 2 
    692                   jm_max(ji) = nlay_i + nlay_s + 1 
    693                    
    694                   ! first layer of snow equation 
    695                   ztrid(ji,2,1)  = 0.0 
    696                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 
    697                   ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
    698                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )  
    699                    
    700                   !                            !---------------------! 
    701                ELSE                            ! cells without snow  ! 
    702                   !                            !---------------------! 
    703                   jm_min(ji)    =  nlay_s + 2 
    704                   jm_max(ji)    =  nlay_i + nlay_s + 1 
    705                    
    706                   ! first layer of ice equation 
    707                   ztrid(ji,jm_min(ji),1)  = 0.0 
    708                   ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)   
    709                   ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
    710                   zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 
    711                    
    712                   ! case of only one layer in the ice (surface & ice equations are altered) 
    713                   IF ( nlay_i == 1 ) THEN 
    714                      ztrid(ji,jm_min(ji),1)  = 0.0 
    715                      ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 
    716                      ztrid(ji,jm_min(ji),3)  = 0.0 
    717                      zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)  & 
    718                         &                    + qcn_ice_1d(ji) ) 
    719                   ENDIF 
    720                    
    721                ENDIF 
    722                ! 
    723                zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
    724                zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
    725                ! 
    726             END DO 
    727             ! 
    728             !------------------------------ 
    729             ! 8) tridiagonal system solving 
    730             !------------------------------ 
    731             ! Solve the tridiagonal system with Gauss elimination method. 
    732             ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    733             jm_maxt = 0 
    734             jm_mint = nlay_i+5 
    735             DO ji = 1, npti 
    736                jm_mint = MIN(jm_min(ji),jm_mint) 
    737                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    738             END DO 
    739              
    740             DO jk = jm_mint+1, jm_maxt 
    741                DO ji = 1, npti 
    742                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
    743                   zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
    744                   zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    745                END DO 
    746             END DO 
    747              
    748             ! ice temperatures 
    749            DO ji = 1, npti 
    750                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    751             END DO 
    752  
    753             DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    754                DO ji = 1, npti 
    755                   jk = jm - nlay_s - 1 
    756                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    757                END DO 
    758             END DO 
    759              
    760             ! snow temperatures       
    761             DO ji = 1, npti 
    762                IF( h_s_1d(ji) > 0._wp ) THEN 
    763                   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) 
    764                ENDIF 
    765             END DO 
    766             ! 
    767             !-------------------------------------------------------------- 
    768             ! 9) Has the scheme converged ?, end of the iterative procedure 
    769             !-------------------------------------------------------------- 
    770             ! check that nowhere it has started to melt 
    771             ! zdti_max is a measure of error, it has to be under zdti_bnd 
    772             zdti_max = 0._wp 
    773  
    774             DO jk = 1, nlay_s 
    775                DO ji = 1, npti 
    776                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    777                   zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    778                END DO 
    779             END DO 
    780              
    781             DO jk = 1, nlay_i 
    782                DO ji = 1, npti 
    783                   ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
    784                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
    785                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    786                END DO 
    787             END DO 
    788  
    789             ! Compute spatial maximum over all errors 
    790             ! note that this could be optimized substantially by iterating only the non-converging points 
    791             IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    792  
    793          ENDIF ! k_jules 
    794           
    795       END DO  ! End of the do while iterative procedure 
    796        
    797       IF( ln_icectl .AND. lwp ) THEN 
    798          WRITE(numout,*) ' zdti_max : ', zdti_max 
    799          WRITE(numout,*) ' iconv    : ', iconv 
    800       ENDIF 
    801        
    802       ! 
    803       !----------------------------- 
    804       ! 10) Fluxes at the interfaces 
    805       !----------------------------- 
    806       ! 
    807       ! --- update conduction fluxes 
    808       ! 
    809       DO ji = 1, npti 
    810          !                                ! surface ice conduction flux 
    811          fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    812             &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    813          !                                ! bottom ice conduction flux 
    814          fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    815       END DO 
    816        
    817       ! 
    818       ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
    819       ! 
    820       IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  ! OFF or SND mode 
    821          DO ji = 1, npti 
    822             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)  
    823          END DO 
    824       ELSE        ! RCV mode 
    825          DO ji = 1, npti 
    826             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji)  
    827          END DO 
    828       ENDIF 
    829        
    830       ! 
    831       ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 
    832       ! 
    833  
    834       IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 
    835           
    836          CALL ice_thd_enmelt        
    837           
    838          !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    839          DO ji = 1, npti 
    840             zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
    841                &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
    842              
    843             IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 
    844                 
    845                IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    846                   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) 
    847                ELSE                          ! case T_su = 0degC 
    848                   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) 
    849                ENDIF 
    850                 
    851             ELSE ! RCV CASE 
    852              
    853                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) 
    854              
    855             ENDIF 
    856             ! 
    857             ! total heat sink to be sent to the ocean 
    858             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
    859             ! 
    860             ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2    
    861             hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
    862             ! 
    863          END DO 
    864          ! 
    865          ! --- SIMIP diagnostics 
    866          ! 
    867          DO ji = 1, npti 
    868             !--- Conduction fluxes (positive downwards) 
    869             diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
    870             diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
    871     
    872             !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    873             zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
    874             IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
    875                t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    876                   &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
    877             ELSE 
    878                t_si_1d(ji) = t_su_1d(ji) 
    879             ENDIF 
    880          END DO 
    881          ! 
    882       ENDIF 
    883       ! 
    884       !--------------------------------------------------------------------------------------- 
    885       ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 
    886       !--------------------------------------------------------------------------------------- 
    887       ! effective conductivity and 1st layer temperature (needed by Met Office) 
    888       DO ji = 1, npti 
    889          IF( h_s_1d(ji) > 0.1_wp ) THEN  
    890             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 
    891          ELSE 
    892             IF( h_i_1d(ji) > 0.1_wp ) THEN 
    893                cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    894             ELSE 
    895                cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp 
    896             ENDIF 
    897          ENDIF 
    898          t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 
    899       END DO 
    900       ! 
    901       IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode 
    902          ! Restore temperatures to their initial values 
    903          t_s_1d    (1:npti,:) = ztsold(1:npti,:) 
    904          t_i_1d    (1:npti,:) = ztiold(1:npti,:) 
    905          qcn_ice_1d(1:npti)   = fc_su (1:npti) 
    906       ENDIF 
    907       ! 
    908    END SUBROUTINE ice_thd_zdf_BL99 
    909  
    910  
    911    SUBROUTINE ice_thd_enmelt 
    912       !!------------------------------------------------------------------- 
    913       !!                   ***  ROUTINE ice_thd_enmelt ***  
    914       !!                  
    915       !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
    916       !! 
    917       !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
    918       !!------------------------------------------------------------------- 
    919       INTEGER  ::   ji, jk   ! dummy loop indices 
    920       REAL(wp) ::   ztmelts  ! local scalar  
    921       !!------------------------------------------------------------------- 
    922       ! 
    923       DO jk = 1, nlay_i             ! Sea ice energy of melting 
    924          DO ji = 1, npti 
    925             ztmelts      = - tmut  * sz_i_1d(ji,jk) 
    926             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 
    927                                                                 !   (sometimes dif scheme produces abnormally high temperatures)    
    928             e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
    929                &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
    930                &                    - rcp  *   ztmelts ) 
    931          END DO 
    932       END DO 
    933       DO jk = 1, nlay_s             ! Snow energy of melting 
    934          DO ji = 1, npti 
    935             e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 
    936          END DO 
    937       END DO 
    938       ! 
    939    END SUBROUTINE ice_thd_enmelt 
    940  
    941  
    94283   SUBROUTINE ice_thd_zdf_init 
    94384      !!----------------------------------------------------------------------- 
    … …  
    976117         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s 
    977118         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
    978      ENDIF 
    979  
     119      ENDIF 
    980120      ! 
    981121      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 
    982122         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 
    983123      ENDIF 
    984  
    985124      !                             !== set the choice of ice vertical thermodynamic formulation ==! 
    986125      ioptio = 0  
    987       !      !--- BL99 thermo dynamics                              (linear liquidus + constant thermal properties) 
    988       IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF 
    989       IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 
     126      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF   ! BL99 thermodynamics (linear liquidus + constant thermal properties) 
     127!!    IF( ln_zdf_XXXX ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_XXXX   ;   ENDIF 
     128      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 
    990129      ! 
    991130   END SUBROUTINE ice_thd_zdf_init 

  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8906 r8984  
    4747   !!   ice_var_itd2      : convert N-cat to jpl-cat 
    4848   !!   ice_var_bv        : brine volume 
     49   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature 
    4950   !!---------------------------------------------------------------------- 
    5051   USE dom_oce        ! ocean space and time domain 
    … …  
    7071   PUBLIC   ice_var_itd2 
    7172   PUBLIC   ice_var_bv            
     73   PUBLIC   ice_var_enthalpy            
    7274 
    7375   !!---------------------------------------------------------------------- 
    … …  
    396398         !                                      ! Slope of the linear profile  
    397399         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti) 
    398          ELSEWHERE                           ;   z_slope_s(1:npti) = 0._wp 
     400         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp 
    399401         END WHERE 
    400402          
    … …  
    822824 
    823825 
     826   SUBROUTINE ice_var_enthalpy 
     827      !!------------------------------------------------------------------- 
     828      !!                   ***  ROUTINE ice_var_enthalpy ***  
     829      !!                  
     830      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
     831      !! 
     832      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     833      !!------------------------------------------------------------------- 
     834      INTEGER  ::   ji, jk   ! dummy loop indices 
     835      REAL(wp) ::   ztmelts  ! local scalar  
     836      !!------------------------------------------------------------------- 
     837      ! 
     838      DO jk = 1, nlay_i             ! Sea ice energy of melting 
     839         DO ji = 1, npti 
     840            ztmelts      = - tmut  * sz_i_1d(ji,jk) 
     841            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 
     842                                                                !   (sometimes zdf scheme produces abnormally high temperatures)    
     843            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
     844               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
     845               &                    - rcp  *   ztmelts ) 
     846         END DO 
     847      END DO 
     848      DO jk = 1, nlay_s             ! Snow energy of melting 
     849         DO ji = 1, npti 
     850            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 
     851         END DO 
     852      END DO 
     853      ! 
     854   END SUBROUTINE ice_var_enthalpy 
     855 
    824856#else 
    825857   !!---------------------------------------------------------------------- 

  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8971 r8984  
    4646   USE lib_fortran    ! to use key_nosignedzero 
    4747#if defined key_lim3 
    48    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
     48   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s 
    4949   USE icethd_dh      ! for CALL ice_thd_snwblow 
    50    USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    5150#endif 
    5251   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
Note: See TracChangeset for help on using the changeset viewer.