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 5575 for branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 – NEMO

Ignore:
Timestamp:
2015-07-09T12:44:22+02:00 (9 years ago)
Author:
davestorkey
Message:

Update UKMO/dev_r5107_hadgem3_cplfld branch to trunk revision 5518
(= branching point of NEMO 3.6_stable).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5473 r5575  
    1919   USE phycst         ! physical constants (ocean directory)  
    2020   USE ice            ! LIM-3 variables 
    21    USE par_ice        ! LIM-3 parameters 
    2221   USE thd_ice        ! LIM-3: thermodynamics 
    2322   USE in_out_manager ! I/O manager 
     
    2524   USE wrk_nemo       ! work arrays 
    2625   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    27    USE sbc_oce, ONLY : lk_cpl 
    2826 
    2927   IMPLICIT NONE 
     
    10098      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    10199      INTEGER ::   minnumeqmin, maxnumeqmax 
     100       
    102101      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
    103102      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
    104       INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     103       
    105104      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    106105      REAL(wp) ::   zg1       =  2._wp        ! 
     
    112111      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    113112      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    114       REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point 
    115       REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    116       REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration 
    117       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
    118       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    119       REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface 
    120       REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function 
    121       REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function 
    122       REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature 
    123       REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4) 
    124       REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice 
    125       REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu 
    126       REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity 
    127       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice 
    128       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice 
    130       REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zswiterm    ! Independent term 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitbis    ! temporary independent term 
    144       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
    145       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms 
     113      REAL(wp) ::   zhsu 
     114       
     115      REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
     116      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
     117      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration 
     118      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness 
     119      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness 
     120      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface 
     121      REAL(wp), POINTER, DIMENSION(:)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
     122      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function 
     123      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function 
     124      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature 
     125      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4) 
     126      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice 
     127      REAL(wp), POINTER, DIMENSION(:)     ::   zihic 
     128       
     129      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity 
     130      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice 
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice 
     133      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice 
     134      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice 
     135      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     136      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat 
     137      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice 
     138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow 
     139      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow 
     140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow 
     141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow 
     142      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
     143      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow 
     144      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow 
     145      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term 
     146      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term 
     147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     148      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms 
     149       
    146150      ! diag errors on heat 
    147       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 
     151      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err 
     152       
     153      ! Mono-category 
     154      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done 
     155      REAL(wp)                            :: zratio_s      ! dummy factor 
     156      REAL(wp)                            :: zratio_i      ! dummy factor 
     157      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation 
     158      REAL(wp)                            :: zhe           ! dummy factor 
     159      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity 
     160      REAL(wp)                            :: zfac          ! dummy factor 
     161      REAL(wp)                            :: zihe          ! dummy factor 
     162      REAL(wp)                            :: zheshth       ! dummy factor 
     163       
     164      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat 
     165       
    148166      !!------------------------------------------------------------------      
    149167      !  
    150       CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
    151       CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    152       CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
    153       CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
    154       CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 
    155       CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis  ) 
    156       CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 
     168      CALL wrk_alloc( jpij, numeqmin, numeqmax ) 
     169      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     170      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 
     171      CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 
     172      CALL wrk_alloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 
     173      CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  ) 
     174      CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 
    157175 
    158176      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 
     
    161179      zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
    162180      DO ji = kideb, kiut 
    163          zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
    164             &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) )  
     181         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     182            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )  
    165183      END DO 
    166184 
     
    168186      ! 1) Initialization                                                            ! 
    169187      !------------------------------------------------------------------------------! 
    170       ! clem clean: replace just ztfs by rtt 
    171188      DO ji = kideb , kiut 
    172          ! is there snow or not 
    173          isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ) 
    174          ! surface temperature of fusion 
    175          ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 
     189         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not 
    176190         ! layer thickness 
    177          zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 
    178          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     191         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i 
     192         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    179193      END DO 
    180194 
     
    188202      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    189203         DO ji = kideb , kiut 
    190             z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 
     204            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 
    191205         END DO 
    192206      END DO 
     
    194208      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    195209         DO ji = kideb , kiut 
    196             z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 
     210            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 
    197211         END DO 
    198212      END DO 
    199213      ! 
    200214      !------------------------------------------------------------------------------| 
    201       ! 2) Radiations                                                                | 
     215      ! 2) Radiation                                                       | 
    202216      !------------------------------------------------------------------------------| 
    203217      ! 
     
    212226      ! zftrice = io.qsr_ice       is below the surface  
    213227      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    214  
     228      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
     229      zhsu = 0.1_wp ! threshold for the computation of i0 
    215230      DO ji = kideb , kiut 
    216231         ! switches 
    217          isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) )  
     232         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  
    218233         ! hs > 0, isnow = 1 
    219          zhsu (ji) = hnzst  ! threshold for the computation of i0 
    220          zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) )      
    221  
    222          i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    223          !fr1_i0_1d = i0 for a thin ice surface 
    224          !fr1_i0_2d = i0 for a thick ice surface 
    225          !            a function of the cloud cover 
    226          ! 
    227          !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 
    228          !formula used in Cice 
     234         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )      
     235 
     236         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    229237      END DO 
    230238 
     
    234242      !------------------------------------------------------- 
    235243      DO ji = kideb , kiut 
    236          zfsw   (ji) =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
    237          zftrice(ji) =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
    238          dzf    (ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
     244         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
     245         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
     246         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
     247         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value 
    239248      END DO 
    240249 
     
    257266 
    258267      DO ji = kideb, kiut           ! ice initialization 
    259          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 
     268         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 
    260269      END DO 
    261270 
     
    263272         DO ji = kideb, kiut 
    264273            !                             ! radiation transmitted below the layer-th ice layer 
    265             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 
     274            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 
    266275            !                             ! radiation absorbed by the layer-th ice layer 
    267276            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    273282      END DO 
    274283 
    275       ! 
    276284      !------------------------------------------------------------------------------| 
    277285      !  3) Iterative procedure begins                                               | 
     
    281289         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
    282290         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter 
    283          t_su_1d   (ji) =  MIN( t_su_1d(ji), ztfs(ji) - ztsu_err )  ! necessary 
    284          zerrit   (ji) =  1000._wp                                ! initial value of error 
     291         t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary 
     292         zerrit (ji) =  1000._wp                                 ! initial value of error 
    285293      END DO 
    286294 
     
    300308      zerritmax =  1000._wp    ! maximal value of error on all points 
    301309 
    302       DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd ) 
     310      DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif ) 
    303311         ! 
    304312         nconv = nconv + 1 
     
    308316         !------------------------------------------------------------------------------| 
    309317         ! 
    310          IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula 
    311             DO ji = kideb , kiut 
    312                ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 
    313                ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
     318         IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula 
     319            DO ji = kideb , kiut 
     320               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     321               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    314322            END DO 
    315323            DO jk = 1, nlay_i-1 
    316324               DO ji = kideb , kiut 
    317                   ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
    318                      MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 
    319                   ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 
     325                  ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
     326                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 
     327                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
    320328               END DO 
    321329            END DO 
    322330         ENDIF 
    323331 
    324          IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    325             DO ji = kideb , kiut 
    326                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt )   & 
    327                   &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt 
     332         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
     333            DO ji = kideb , kiut 
     334               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
     335                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 
    328336               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    329337            END DO 
    330338            DO jk = 1, nlay_i-1 
    331339               DO ji = kideb , kiut 
    332                   ztcond_i(ji,jk) = rcdic +                                                                     &  
    333                      &                 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                          & 
    334                      &                 / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)   & 
    335                      &               - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt 
     340                  ztcond_i(ji,jk) = rcdic +                                                                       &  
     341                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              & 
     342                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   & 
     343                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 
    336344                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
    337345               END DO 
    338346            END DO 
    339347            DO ji = kideb , kiut 
    340                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt)   & 
    341                   &                        - 0.011_wp * ( t_bo_1d(ji) - rtt 
     348               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
     349                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 
    342350               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    343351            END DO 
    344352         ENDIF 
    345          ! 
    346          !------------------------------------------------------------------------------| 
    347          !  5) kappa factors                                                            | 
    348          !------------------------------------------------------------------------------| 
    349          ! 
     353          
     354         ! 
     355         !------------------------------------------------------------------------------| 
     356         !  5) G(he) - enhancement of thermal conductivity in mono-category case        | 
     357         !------------------------------------------------------------------------------| 
     358         ! 
     359         ! Computation of effective thermal conductivity G(h) 
     360         ! Used in mono-category case only to simulate an ITD implicitly 
     361         ! Fichefet and Morales Maqueda, JGR 1997 
     362 
     363         zghe(:) = 1._wp 
     364 
     365         SELECT CASE ( nn_monocat ) 
     366 
     367         CASE (1,3) ! LIM3 
     368 
     369            zepsilon = 0.1_wp 
     370            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 
     371 
     372            DO ji = kideb, kiut 
     373    
     374               ! Mean sea ice thermal conductivity 
     375               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp ) 
     376 
     377               ! Effective thickness he (zhe) 
     378               zfac     = 1._wp / ( rcdsn + zkimean ) 
     379               zratio_s = rcdsn   * zfac 
     380               zratio_i = zkimean * zfac 
     381               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     382 
     383               ! G(he) 
     384               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
     385               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 
     386    
     387               ! Impose G(he) < 2. 
     388               zghe(ji) = MIN( zghe(ji), 2._wp ) 
     389 
     390            END DO 
     391 
     392         END SELECT 
     393 
     394         ! 
     395         !------------------------------------------------------------------------------| 
     396         !  6) kappa factors                                                            | 
     397         !------------------------------------------------------------------------------| 
     398         ! 
     399         !--- Snow 
    350400         DO ji = kideb, kiut 
    351  
    352             !-- Snow kappa factors 
    353             zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
    354             zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
     401            zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
     402            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac 
     403            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac 
    355404         END DO 
    356405 
    357406         DO jk = 1, nlay_s-1 
    358407            DO ji = kideb , kiut 
    359                zkappa_s(ji,jk)  = 2.0 * rcdsn / & 
    360                   MAX(epsi10,2.0*zh_s(ji)) 
    361             END DO 
    362          END DO 
    363  
     408               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
     409            END DO 
     410         END DO 
     411 
     412         !--- Ice 
    364413         DO jk = 1, nlay_i-1 
    365414            DO ji = kideb , kiut 
    366                !-- Ice kappa factors 
    367                zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ & 
    368                   MAX(epsi10,2.0*zh_i(ji))  
    369             END DO 
    370          END DO 
    371  
    372          DO ji = kideb , kiut 
    373             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
    374             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    375             !-- Interface 
    376             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    377                (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    378             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
    379                + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 
    380          END DO 
    381          ! 
    382          !------------------------------------------------------------------------------| 
    383          ! 6) Sea ice specific heat, eta factors                                        | 
     415               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 
     416            END DO 
     417         END DO 
     418 
     419         !--- Snow-ice interface 
     420         DO ji = kideb , kiut 
     421            zfac                  = 1./ MAX( epsi10 , zh_i(ji) ) 
     422            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
     423            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
     424            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
     425           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
     426            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     427         END DO 
     428 
     429         ! 
     430         !------------------------------------------------------------------------------| 
     431         ! 7) Sea ice specific heat, eta factors                                        | 
    384432         !------------------------------------------------------------------------------| 
    385433         ! 
     
    387435            DO ji = kideb , kiut 
    388436               ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    389                zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 
    390                   MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 
    391                zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 
    392                   epsi10) 
     437               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
     438               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 
    393439            END DO 
    394440         END DO 
     
    397443            DO ji = kideb , kiut 
    398444               ztstemp(ji,jk) = t_s_1d(ji,jk) 
    399                zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
    400             END DO 
    401          END DO 
    402          ! 
    403          !------------------------------------------------------------------------------| 
    404          ! 7) surface flux computation                                                  | 
    405          !------------------------------------------------------------------------------| 
    406          ! 
    407          IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case 
     445               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 
     446            END DO 
     447         END DO 
     448 
     449         ! 
     450         !------------------------------------------------------------------------------| 
     451         ! 8) surface flux computation                                                  | 
     452         !------------------------------------------------------------------------------| 
     453         ! 
     454         IF ( ln_it_qnsice ) THEN  
    408455            DO ji = kideb , kiut 
    409456               ! update of the non solar flux according to the update in T_su 
     
    415462         DO ji = kideb , kiut 
    416463            ! update incoming flux 
    417             zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation 
    418                + qns_ice_1d(ji)                   ! non solar total flux  
    419             ! (LWup, LWdw, SH, LH) 
    420          END DO 
    421  
    422          ! 
    423          !------------------------------------------------------------------------------| 
    424          ! 8) tridiagonal system terms                                                  | 
     464            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation 
     465               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH) 
     466         END DO 
     467 
     468         ! 
     469         !------------------------------------------------------------------------------| 
     470         ! 9) tridiagonal system terms                                                  | 
    425471         !------------------------------------------------------------------------------| 
    426472         ! 
     
    437483               ztrid(ji,numeq,2) = 0. 
    438484               ztrid(ji,numeq,3) = 0. 
    439                zswiterm(ji,numeq)= 0. 
    440                zswitbis(ji,numeq)= 0. 
     485               zindterm(ji,numeq)= 0. 
     486               zindtbis(ji,numeq)= 0. 
    441487               zdiagbis(ji,numeq)= 0. 
    442488            ENDDO 
     
    445491         DO numeq = nlay_s + 2, nlay_s + nlay_i  
    446492            DO ji = kideb , kiut 
    447                jk              = numeq - nlay_s - 1 
    448                ztrid(ji,numeq,1)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 
    449                ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 
    450                   zkappa_i(ji,jk)) 
    451                ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk) 
    452                zswiterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
    453                   zradab_i(ji,jk) 
     493               jk                 = numeq - nlay_s - 1 
     494               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
     495               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
     496               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk) 
     497               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    454498            END DO 
    455499         ENDDO 
     
    459503            !!ice bottom term 
    460504            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    461             ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & 
    462                +  zkappa_i(ji,nlay_i-1) ) 
     505            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
    463506            ztrid(ji,numeq,3)  =  0.0 
    464             zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
    465                ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 
    466                *  t_bo_1d(ji) )  
     507            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     508               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    467509         ENDDO 
    468510 
    469511 
    470512         DO ji = kideb , kiut 
    471             IF ( ht_s_1d(ji).gt.0.0 ) THEN 
     513            IF ( ht_s_1d(ji) > 0.0 ) THEN 
    472514               ! 
    473515               !------------------------------------------------------------------------------| 
     
    477519               !!snow interior terms (bottom equation has the same form as the others) 
    478520               DO numeq = 3, nlay_s + 1 
    479                   jk =  numeq - 1 
    480                   ztrid(ji,numeq,1)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 
    481                   ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 
    482                      zkappa_s(ji,jk) ) 
     521                  jk                  =  numeq - 1 
     522                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     523                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    483524                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    484                   zswiterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
    485                      zradab_s(ji,jk) 
     525                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    486526               END DO 
    487527 
     
    489529               IF ( nlay_i.eq.1 ) THEN 
    490530                  ztrid(ji,nlay_s+2,3)    =  0.0 
    491                   zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
    492                      t_bo_1d(ji)  
     531                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    493532               ENDIF 
    494533 
    495                IF ( t_su_1d(ji) .LT. rtt ) THEN 
     534               IF ( t_su_1d(ji) < rt0 ) THEN 
    496535 
    497536                  !------------------------------------------------------------------------------| 
     
    503542 
    504543                  !!surface equation 
    505                   ztrid(ji,1,1) = 0.0 
    506                   ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 
    507                   ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 
    508                   zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)  - zf(ji) 
     544                  ztrid(ji,1,1)  = 0.0 
     545                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0) 
     546                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
     547                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji) 
    509548 
    510549                  !!first layer of snow equation 
    511                   ztrid(ji,2,1)  =  - zkappa_s(ji,0)*zg1s*zeta_s(ji,1) 
    512                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 
     550                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
     551                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    513552                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    514                   zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
     553                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    515554 
    516555               ELSE  
     
    526565                  !!first layer of snow equation 
    527566                  ztrid(ji,2,1)  =  0.0 
    528                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & 
    529                      zkappa_s(ji,0) * zg1s ) 
     567                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    530568                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    531                   zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
    532                      ( zradab_s(ji,1) +                         & 
    533                      zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     569                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  & 
     570                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    534571               ENDIF 
    535572            ELSE 
     
    539576               !------------------------------------------------------------------------------| 
    540577               ! 
    541                IF (t_su_1d(ji) .LT. rtt) THEN 
     578               IF ( t_su_1d(ji) < rt0 ) THEN 
    542579                  ! 
    543580                  !------------------------------------------------------------------------------| 
     
    553590                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
    554591                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    555                   zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
     592                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
    556593 
    557594                  !!first layer of ice equation 
    558595                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    559                   ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &  
    560                      + zkappa_i(ji,0) * zg1 ) 
    561                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)   
    562                   zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
     596                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
     597                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)   
     598                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    563599 
    564600                  !!case of only one layer in the ice (surface & ice equations are altered) 
    565601 
    566                   IF (nlay_i.eq.1) THEN 
     602                  IF ( nlay_i == 1 ) THEN 
    567603                     ztrid(ji,numeqmin(ji),1)    =  0.0 
    568                      ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0)*2.0 
    569                      ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0)*2.0 
    570                      ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0)*2.0*zeta_i(ji,1) 
    571                      ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 
    572                         zkappa_i(ji,1)) 
     604                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0 
     605                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0 
     606                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     607                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    573608                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    574609 
    575                      zswiterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
    576                         ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 
     610                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) * & 
     611                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
    577612                  ENDIF 
    578613 
     
    590625                  !!first layer of ice equation 
    591626                  ztrid(ji,numeqmin(ji),1)      =  0.0 
    592                   ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & 
    593                      zg1)   
     627                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    594628                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    595                   zswiterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
    596                      zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
     629                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) * & 
     630                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    597631 
    598632                  !!case of only one layer in the ice (surface & ice equations are altered) 
    599                   IF (nlay_i.eq.1) THEN 
     633                  IF ( nlay_i == 1 ) THEN 
    600634                     ztrid(ji,numeqmin(ji),1)  =  0.0 
    601                      ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 
    602                         zkappa_i(ji,1)) 
     635                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    603636                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    604                      zswiterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
    605                         (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 
    606                         + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 
     637                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
     638                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    607639                  ENDIF 
    608640 
     
    614646         ! 
    615647         !------------------------------------------------------------------------------| 
    616          ! 9) tridiagonal system solving                                                | 
     648         ! 10) tridiagonal system solving                                               | 
    617649         !------------------------------------------------------------------------------| 
    618650         ! 
     
    626658 
    627659         DO ji = kideb , kiut 
    628             zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji)) 
     660            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
    629661            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
    630662            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin) 
     
    635667            DO ji = kideb , kiut 
    636668               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 
    637                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 
    638                   ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 
    639                zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 
    640                   zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
     669               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1) 
     670               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1) 
    641671            END DO 
    642672         END DO 
     
    644674         DO ji = kideb , kiut 
    645675            ! ice temperatures 
    646             t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
    647          END DO 
    648  
    649          DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 
     676            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 
     677         END DO 
     678 
     679         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 
    650680            DO ji = kideb , kiut 
    651681               jk    =  numeq - nlay_s - 1 
    652                t_i_1d(ji,jk)  =  (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 
    653                   t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 
     682               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) 
    654683            END DO 
    655684         END DO 
     
    657686         DO ji = kideb , kiut 
    658687            ! snow temperatures       
    659             IF (ht_s_1d(ji).GT.0._wp) & 
    660                t_s_1d(ji,nlay_s)     =  (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    661                *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 
    662                *        MAX(0.0,SIGN(1.0,ht_s_1d(ji)))  
     688            IF (ht_s_1d(ji) > 0._wp) & 
     689               t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     690               &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) )  
    663691 
    664692            ! surface temperature 
    665             isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  ) ) 
     693            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 
    666694            ztsubit(ji) = t_su_1d(ji) 
    667             IF( t_su_1d(ji) < ztfs(ji) ) & 
    668                t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
    669                &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
     695            IF( t_su_1d(ji) < rt0 ) & 
     696               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & 
     697               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    670698         END DO 
    671699         ! 
    672700         !-------------------------------------------------------------------------- 
    673          !  10) Has the scheme converged ?, end of the iterative procedure         | 
     701         !  11) Has the scheme converged ?, end of the iterative procedure         | 
    674702         !-------------------------------------------------------------------------- 
    675703         ! 
    676704         ! check that nowhere it has started to melt 
    677          ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 
    678          DO ji = kideb , kiut 
    679             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp  ) 
    680             zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )      
     705         ! zerrit(ji) is a measure of error, it has to be under terr_dif 
     706         DO ji = kideb , kiut 
     707            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  ) 
     708            zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) )      
    681709         END DO 
    682710 
    683711         DO jk  =  1, nlay_s 
    684712            DO ji = kideb , kiut 
    685                t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  ) 
    686                zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 
     713               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  ) 
     714               zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 
    687715            END DO 
    688716         END DO 
     
    690718         DO jk  =  1, nlay_i 
    691719            DO ji = kideb , kiut 
    692                ztmelt_i        = -tmut * s_i_1d(ji,jk) + rtt  
    693                t_i_1d(ji,jk) =  MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 
    694                zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 
     720               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0  
     721               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) 
     722               zerrit(ji)    =  MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) 
    695723            END DO 
    696724         END DO 
     
    706734      END DO  ! End of the do while iterative procedure 
    707735 
    708       IF( ln_nicep .AND. lwp ) THEN 
     736      IF( ln_icectl .AND. lwp ) THEN 
    709737         WRITE(numout,*) ' zerritmax : ', zerritmax 
    710738         WRITE(numout,*) ' nconv     : ', nconv 
     
    713741      ! 
    714742      !-------------------------------------------------------------------------! 
    715       !   11) Fluxes at the interfaces                                          ! 
     743      !   12) Fluxes at the interfaces                                          ! 
    716744      !-------------------------------------------------------------------------! 
    717745      DO ji = kideb, kiut 
    718          ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)  
    719          IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 
    720746         !                                ! surface ice conduction flux 
    721          isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 
    722          fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    723             &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
     747         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 
     748         fc_su(ji)       =  -           isnow(ji)  * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
     749            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    724750         !                                ! bottom ice conduction flux 
    725751         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    726752      END DO 
     753 
     754      ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
     755      CALL lim_thd_enmelt( kideb, kiut ) 
     756 
     757      ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 
     758      IF ( ln_it_qnsice ) THEN 
     759         DO ji = kideb, kiut 
     760            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
     761         END DO 
     762      END IF 
     763 
     764      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
     765      DO ji = kideb, kiut 
     766         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     767            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
     768         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
     769            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     770         ELSE                          ! case T_su = 0degC 
     771            zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     772         ENDIF 
     773         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     774 
     775         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
     776         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     777      END DO  
    727778 
    728779      !----------------------------------------- 
     
    730781      !----------------------------------------- 
    731782      DO ji = kideb, kiut 
    732          IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC 
     783         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    733784            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & 
    734785               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    735          ELSE                         ! case T_su = 0degC 
     786         ELSE                          ! case T_su = 0degC 
    736787            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    & 
    737788               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    738789         ENDIF 
    739       END DO 
    740  
    741       ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
    742       CALL lim_thd_enmelt( kideb, kiut ) 
    743  
    744       ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    745       DO ji = kideb, kiut 
    746          zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
    747             &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
    748          zhfx_err(ji)   = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice )  
    749          hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
    750       END DO  
    751  
    752       ! diagnose external surface (forced case) or bottom (forced case) from heat conservation 
    753       IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed 
    754          ! 
    755          DO ji = kideb, kiut 
    756             qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji) 
    757             fc_su     (ji) = fc_su(ji)      - zhfx_err(ji) 
    758          END DO 
    759          ! 
    760       ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed 
    761          ! 
    762          DO ji = kideb, kiut 
    763             fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji) 
    764          END DO 
    765          ! 
    766       ENDIF 
    767  
    768       ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2) 
    769       DO ji = kideb, kiut 
    770          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    771          hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 
    772       END DO 
    773     
     790         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
     791         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 
     792      END DO    
    774793      ! 
    775       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    776       CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    777       CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
    778       CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 
    779          &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
    780       CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
    781       CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 
    782       CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 
     794      CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 
     795      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     796      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 
     797      CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
     798      CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
     799      CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 
     800      CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 
    783801      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 
    784802 
     
    801819      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    802820         DO ji = kideb, kiut 
    803             ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt  
    804             rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 
    805             q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                             & 
    806                &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) )   & 
    807                &                   - rcp  *                 ( ztmelts-rtt )  )  
     821            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0 
     822            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 
     823                                                          !   (sometimes dif scheme produces abnormally high temperatures)    
     824            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           & 
     825               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   & 
     826               &                    - rcp  *         ( ztmelts-rt0 )  )  
    808827         END DO 
    809828      END DO 
    810829      DO jk = 1, nlay_s             ! Snow energy of melting 
    811830         DO ji = kideb, kiut 
    812             q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 
     831            q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 
    813832         END DO 
    814833      END DO 
Note: See TracChangeset for help on using the changeset viewer.