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 5075 for branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 – NEMO

Ignore:
Timestamp:
2015-02-11T11:50:34+01:00 (9 years ago)
Author:
timgraham
Message:

Upgraded branch to current head of trunk (r5072) so it can be used with the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4333 r5075  
    2525   USE wrk_nemo       ! work arrays 
    2626   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     27   USE sbc_oce, ONLY : lk_cpl 
    2728 
    2829   IMPLICIT NONE 
     
    3132   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3233 
    33    REAL(wp) ::   epsi10      =  1.e-10_wp    ! 
    3434   !!---------------------------------------------------------------------- 
    3535   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     
    3939CONTAINS 
    4040 
    41    SUBROUTINE lim_thd_dif( kideb , kiut , jl ) 
     41   SUBROUTINE lim_thd_dif( kideb , kiut ) 
    4242      !!------------------------------------------------------------------ 
    4343      !!                ***  ROUTINE lim_thd_dif  *** 
     
    7474      !! 
    7575      !! ** Inputs / Ouputs : (global commons) 
    76       !!           surface temperature : t_su_b 
    77       !!           ice/snow temperatures   : t_i_b, t_s_b 
    78       !!           ice salinities          : s_i_b 
     76      !!           surface temperature : t_su_1d 
     77      !!           ice/snow temperatures   : t_i_1d, t_s_1d 
     78      !!           ice salinities          : s_i_1d 
    7979      !!           number of layers in the ice/snow: nlay_i, nlay_s 
    8080      !!           profile of the ice/snow layers : z_i, z_s 
    81       !!           total ice/snow thickness : ht_i_b, ht_s_b 
     81      !!           total ice/snow thickness : ht_i_1d, ht_s_1d 
    8282      !! 
    8383      !! ** External :  
     
    9191      !!           (04-2007) Energy conservation tested by M. Vancoppenolle 
    9292      !!------------------------------------------------------------------ 
    93       INTEGER , INTENT (in) ::   kideb   ! Start point on which the  the computation is applied 
    94       INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied 
    95       INTEGER , INTENT (in) ::   jl      ! Category number 
     93      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    9694 
    9795      !! * Local variables 
     
    9997      INTEGER ::   ii, ij      ! temporary dummy loop index 
    10098      INTEGER ::   numeq       ! current reference number of equation 
    101       INTEGER ::   layer       ! vertical dummy loop index  
     99      INTEGER ::   jk       ! vertical dummy loop index  
    102100      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    103101      INTEGER ::   minnumeqmin, maxnumeqmax 
    104       INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation 
    105       INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    106       INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     102      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
     103      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
     104      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    107105      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    108106      REAL(wp) ::   zg1       =  2._wp        ! 
    109107      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat 
    110108      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13) 
    111       REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow 
     109      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow 
    112110      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
     111      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C  
    113112      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    114113      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    115       REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point 
    116       REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure ) 
    117       REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration 
    118       REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness 
    119       REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness 
    120       REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface 
    121       REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function 
    122       REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function 
    123       REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature 
    124       REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4) 
    125       REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice 
    126       REAL(wp), DIMENSION(kiut) ::   zihic, zhsu 
    127       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
    128       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
    129       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
    130       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
    131       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice 
    132       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
    133       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    134       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat 
    135       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice 
    136       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
    137       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
    138       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
    139       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow 
    140       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence 
    141       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow 
    142       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow 
    143       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term 
    144       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term 
    145       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis 
    146       REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms 
     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 
     146      ! diag errors on heat 
     147      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 
    147148      !!------------------------------------------------------------------      
    148149      !  
     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 ) 
     157 
     158      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 
     159 
     160      ! --- diag error on heat diffusion - PART 1 --- ! 
     161      zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
     162      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 ) )  
     165      END DO 
     166 
    149167      !------------------------------------------------------------------------------! 
    150168      ! 1) Initialization                                                            ! 
    151169      !------------------------------------------------------------------------------! 
    152       ! 
     170      ! clem clean: replace just ztfs by rtt 
    153171      DO ji = kideb , kiut 
    154172         ! is there snow or not 
    155          isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
     173         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ) 
    156174         ! surface temperature of fusion 
    157 !!gm ???  ztfs(ji) = rtt !!!???? 
    158175         ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 
    159176         ! layer thickness 
    160          zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
    161          zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 
     177         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 
     178         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
    162179      END DO 
    163180 
     
    169186      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer 
    170187 
    171       DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    172          DO ji = kideb , kiut 
    173             z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s ) 
    174          END DO 
    175       END DO 
    176  
    177       DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    178          DO ji = kideb , kiut 
    179             z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i ) 
     188      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
     189         DO ji = kideb , kiut 
     190            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 
     191         END DO 
     192      END DO 
     193 
     194      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
     195         DO ji = kideb , kiut 
     196            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 
    180197         END DO 
    181198      END DO 
     
    194211      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
    195212      ! zftrice = io.qsr_ice       is below the surface  
    196       ! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
     213      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    197214 
    198215      DO ji = kideb , kiut 
    199216         ! switches 
    200          isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
     217         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  )  
    201218         ! hs > 0, isnow = 1 
    202219         zhsu (ji) = hnzst  ! threshold for the computation of i0 
    203          zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )      
     220         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) )      
    204221 
    205222         i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
     
    208225         !            a function of the cloud cover 
    209226         ! 
    210          !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) 
     227         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 
    211228         !formula used in Cice 
    212229      END DO 
     
    230247      END DO 
    231248 
    232       DO layer = 1, nlay_s          ! Radiation through snow 
     249      DO jk = 1, nlay_s          ! Radiation through snow 
    233250         DO ji = kideb, kiut 
    234251            !                             ! radiation transmitted below the layer-th snow layer 
    235             zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) ) 
     252            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 
    236253            !                             ! radiation absorbed by the layer-th snow layer 
    237             zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer) 
     254            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
    238255         END DO 
    239256      END DO 
     
    243260      END DO 
    244261 
    245       DO layer = 1, nlay_i          ! Radiation through ice 
     262      DO jk = 1, nlay_i          ! Radiation through ice 
    246263         DO ji = kideb, kiut 
    247264            !                             ! radiation transmitted below the layer-th ice layer 
    248             zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) ) 
     265            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 
    249266            !                             ! radiation absorbed by the layer-th ice layer 
    250             zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer) 
     267            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
    251268         END DO 
    252269      END DO 
    253270 
    254271      DO ji = kideb, kiut           ! Radiation transmitted below the ice 
    255          fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 
    256       END DO 
    257  
    258       ! +++++ 
    259       ! just to check energy conservation 
    260       DO ji = kideb, kiut 
    261          ii = MOD( npb(ji) - 1 , jpi ) + 1 
    262          ij =    ( npb(ji) - 1 ) / jpi + 1 
    263          fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 
    264       END DO 
    265       ! +++++ 
    266  
    267       DO layer = 1, nlay_i 
    268          DO ji = kideb, kiut 
    269             radab(ji,layer) = zradab_i(ji,layer) 
    270          END DO 
     272         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i)  
    271273      END DO 
    272274 
     
    277279      ! 
    278280      DO ji = kideb, kiut        ! Old surface temperature 
    279          ztsuold  (ji) =  t_su_b(ji)                              ! temperature at the beg of iter pr. 
    280          ztsuoldit(ji) =  t_su_b(ji)                              ! temperature at the previous iter 
    281          t_su_b   (ji) =  MIN( t_su_b(ji), ztfs(ji)-0.00001 )     ! necessary 
     281         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
     282         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 
    282284         zerrit   (ji) =  1000._wp                                ! initial value of error 
    283285      END DO 
    284286 
    285       DO layer = 1, nlay_s       ! Old snow temperature 
    286          DO ji = kideb , kiut 
    287             ztsold(ji,layer) =  t_s_b(ji,layer) 
    288          END DO 
    289       END DO 
    290  
    291       DO layer = 1, nlay_i       ! Old ice temperature 
    292          DO ji = kideb , kiut 
    293             ztiold(ji,layer) =  t_i_b(ji,layer) 
     287      DO jk = 1, nlay_s       ! Old snow temperature 
     288         DO ji = kideb , kiut 
     289            ztsb(ji,jk) =  t_s_1d(ji,jk) 
     290         END DO 
     291      END DO 
     292 
     293      DO jk = 1, nlay_i       ! Old ice temperature 
     294         DO ji = kideb , kiut 
     295            ztib(ji,jk) =  t_i_1d(ji,jk) 
    294296         END DO 
    295297      END DO 
     
    308310         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula 
    309311            DO ji = kideb , kiut 
    310                ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt) 
     312               ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 
    311313               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
    312314            END DO 
    313             DO layer = 1, nlay_i-1 
     315            DO jk = 1, nlay_i-1 
    314316               DO ji = kideb , kiut 
    315                   ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  & 
    316                      MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
    317                   ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 
     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) 
    318320               END DO 
    319321            END DO 
     
    322324         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    323325            DO ji = kideb , kiut 
    324                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   & 
    325                   &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )   
     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 )   
    326328               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    327329            END DO 
    328             DO layer = 1, nlay_i-1 
     330            DO jk = 1, nlay_i-1 
    329331               DO ji = kideb , kiut 
    330                   ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    331                      &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
    332                      &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    333                   ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
     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 )   
     336                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
    334337               END DO 
    335338            END DO 
    336339            DO ji = kideb , kiut 
    337                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
    338                   &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
     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 )   
    339342               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    340343            END DO 
     
    352355         END DO 
    353356 
    354          DO layer = 1, nlay_s-1 
    355             DO ji = kideb , kiut 
    356                zkappa_s(ji,layer)  = 2.0 * rcdsn / & 
     357         DO jk = 1, nlay_s-1 
     358            DO ji = kideb , kiut 
     359               zkappa_s(ji,jk)  = 2.0 * rcdsn / & 
    357360                  MAX(epsi10,2.0*zh_s(ji)) 
    358361            END DO 
    359362         END DO 
    360363 
    361          DO layer = 1, nlay_i-1 
     364         DO jk = 1, nlay_i-1 
    362365            DO ji = kideb , kiut 
    363366               !-- Ice kappa factors 
    364                zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ & 
     367               zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ & 
    365368                  MAX(epsi10,2.0*zh_i(ji))  
    366369            END DO 
     
    381384         !------------------------------------------------------------------------------| 
    382385         ! 
    383          DO layer = 1, nlay_i 
    384             DO ji = kideb , kiut 
    385                ztitemp(ji,layer)   = t_i_b(ji,layer) 
    386                zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 
    387                   MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 
    388                zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 
     386         DO jk = 1, nlay_i 
     387            DO ji = kideb , kiut 
     388               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), & 
    389392                  epsi10) 
    390393            END DO 
    391394         END DO 
    392395 
    393          DO layer = 1, nlay_s 
    394             DO ji = kideb , kiut 
    395                ztstemp(ji,layer) = t_s_b(ji,layer) 
    396                zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
     396         DO jk = 1, nlay_s 
     397            DO ji = kideb , kiut 
     398               ztstemp(ji,jk) = t_s_1d(ji,jk) 
     399               zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
    397400            END DO 
    398401         END DO 
     
    402405         !------------------------------------------------------------------------------| 
    403406         ! 
    404          DO ji = kideb , kiut 
    405  
    406             ! update of the non solar flux according to the update in T_su 
    407             qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * &  
    408                ( t_su_b(ji) - ztsuoldit(ji) ) 
    409  
     407         IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case 
     408            DO ji = kideb , kiut 
     409               ! update of the non solar flux according to the update in T_su 
     410               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 
     411            END DO 
     412         ENDIF 
     413 
     414         ! Update incoming flux 
     415         DO ji = kideb , kiut 
    410416            ! update incoming flux 
    411417            zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation 
    412                + qnsr_ice_1d(ji)           ! non solar total flux  
     418               + qns_ice_1d(ji)                   ! non solar total flux  
    413419            ! (LWup, LWdw, SH, LH) 
    414  
    415420         END DO 
    416421 
     
    427432         !!ice interior terms (top equation has the same form as the others) 
    428433 
    429          DO numeq=1,jkmax+2 
     434         DO numeq=1,nlay_i+3 
    430435            DO ji = kideb , kiut 
    431436               ztrid(ji,numeq,1) = 0. 
    432437               ztrid(ji,numeq,2) = 0. 
    433438               ztrid(ji,numeq,3) = 0. 
    434                zindterm(ji,numeq)= 0. 
    435                zindtbis(ji,numeq)= 0. 
     439               zswiterm(ji,numeq)= 0. 
     440               zswitbis(ji,numeq)= 0. 
    436441               zdiagbis(ji,numeq)= 0. 
    437442            ENDDO 
     
    440445         DO numeq = nlay_s + 2, nlay_s + nlay_i  
    441446            DO ji = kideb , kiut 
    442                layer              = numeq - nlay_s - 1 
    443                ztrid(ji,numeq,1)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer-1) 
    444                ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + & 
    445                   zkappa_i(ji,layer)) 
    446                ztrid(ji,numeq,3)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer) 
    447                zindterm(ji,numeq) =  ztiold(ji,layer) + zeta_i(ji,layer)* & 
    448                   zradab_i(ji,layer) 
     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) 
    449454            END DO 
    450455         ENDDO 
     
    457462               +  zkappa_i(ji,nlay_i-1) ) 
    458463            ztrid(ji,numeq,3)  =  0.0 
    459             zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
     464            zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
    460465               ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 
    461                *  t_bo_b(ji) )  
     466               *  t_bo_1d(ji) )  
    462467         ENDDO 
    463468 
    464469 
    465470         DO ji = kideb , kiut 
    466             IF ( ht_s_b(ji).gt.0.0 ) THEN 
     471            IF ( ht_s_1d(ji).gt.0.0 ) THEN 
    467472               ! 
    468473               !------------------------------------------------------------------------------| 
     
    472477               !!snow interior terms (bottom equation has the same form as the others) 
    473478               DO numeq = 3, nlay_s + 1 
    474                   layer =  numeq - 1 
    475                   ztrid(ji,numeq,1)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer-1) 
    476                   ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + & 
    477                      zkappa_s(ji,layer) ) 
    478                   ztrid(ji,numeq,3)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer) 
    479                   zindterm(ji,numeq)  =  ztsold(ji,layer) + zeta_s(ji,layer)* & 
    480                      zradab_s(ji,layer) 
     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) ) 
     483                  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) 
    481486               END DO 
    482487 
     
    484489               IF ( nlay_i.eq.1 ) THEN 
    485490                  ztrid(ji,nlay_s+2,3)    =  0.0 
    486                   zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
    487                      t_bo_b(ji)  
     491                  zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
     492                     t_bo_1d(ji)  
    488493               ENDIF 
    489494 
    490                IF ( t_su_b(ji) .LT. rtt ) THEN 
     495               IF ( t_su_1d(ji) .LT. rtt ) THEN 
    491496 
    492497                  !------------------------------------------------------------------------------| 
     
    501506                  ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 
    502507                  ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 
    503                   zindterm(ji,1) = dzf(ji)*t_su_b(ji)   - zf(ji) 
     508                  zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji) 
    504509 
    505510                  !!first layer of snow equation 
     
    507512                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 
    508513                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    509                   zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
     514                  zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
    510515 
    511516               ELSE  
     
    524529                     zkappa_s(ji,0) * zg1s ) 
    525530                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    526                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *            & 
     531                  zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
    527532                     ( zradab_s(ji,1) +                         & 
    528                      zkappa_s(ji,0) * zg1s * t_su_b(ji) )  
     533                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    529534               ENDIF 
    530535            ELSE 
     
    534539               !------------------------------------------------------------------------------| 
    535540               ! 
    536                IF (t_su_b(ji) .LT. rtt) THEN 
     541               IF (t_su_1d(ji) .LT. rtt) THEN 
    537542                  ! 
    538543                  !------------------------------------------------------------------------------| 
     
    548553                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
    549554                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    550                   zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_b(ji) - zf(ji) 
     555                  zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
    551556 
    552557                  !!first layer of ice equation 
     
    555560                     + zkappa_i(ji,0) * zg1 ) 
    556561                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)   
    557                   zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
     562                  zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
    558563 
    559564                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    568573                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    569574 
    570                      zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1)* & 
    571                         ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) ) 
     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) ) 
    572577                  ENDIF 
    573578 
     
    588593                     zg1)   
    589594                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    590                   zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
    591                      zkappa_i(ji,0) * zg1 * t_su_b(ji) )  
     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) )  
    592597 
    593598                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    597602                        zkappa_i(ji,1)) 
    598603                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    599                      zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1)* & 
    600                         (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) & 
    601                         + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.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 
    602607                  ENDIF 
    603608 
     
    618623 
    619624         maxnumeqmax = 0 
    620          minnumeqmin = jkmax+4 
    621  
    622          DO ji = kideb , kiut 
    623             zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
     625         minnumeqmin = nlay_i+5 
     626 
     627         DO ji = kideb , kiut 
     628            zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji)) 
    624629            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
    625630            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin) 
     
    627632         END DO 
    628633 
    629          DO layer = minnumeqmin+1, maxnumeqmax 
    630             DO ji = kideb , kiut 
    631                numeq               =  min(max(numeqmin(ji)+1,layer),numeqmax(ji)) 
     634         DO jk = minnumeqmin+1, maxnumeqmax 
     635            DO ji = kideb , kiut 
     636               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 
    632637               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 
    633638                  ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 
    634                zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 
    635                   zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
     639               zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 
     640                  zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
    636641            END DO 
    637642         END DO 
     
    639644         DO ji = kideb , kiut 
    640645            ! ice temperatures 
    641             t_i_b(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
     646            t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
    642647         END DO 
    643648 
    644649         DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 
    645650            DO ji = kideb , kiut 
    646                layer    =  numeq - nlay_s - 1 
    647                t_i_b(ji,layer)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 
    648                   t_i_b(ji,layer+1))/zdiagbis(ji,numeq) 
     651               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) 
    649654            END DO 
    650655         END DO 
     
    652657         DO ji = kideb , kiut 
    653658            ! snow temperatures       
    654             IF (ht_s_b(ji).GT.0._wp) & 
    655                t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    656                *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 
    657                *        MAX(0.0,SIGN(1.0,ht_s_b(ji)))  
     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)))  
    658663 
    659664            ! surface temperature 
    660             isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  ) 
    661             ztsuoldit(ji) = t_su_b(ji) 
    662             IF( t_su_b(ji) < ztfs(ji) ) & 
    663                t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)   & 
    664                &          + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
     665            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  )  ) 
     666            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))   
    665670         END DO 
    666671         ! 
     
    672677         ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 
    673678         DO ji = kideb , kiut 
    674             t_su_b(ji) =  MAX(  MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp  ) 
    675             zerrit(ji) =  ABS( t_su_b(ji) - ztsuoldit(ji) )      
    676          END DO 
    677  
    678          DO layer  =  1, nlay_s 
    679             DO ji = kideb , kiut 
    680                ii = MOD( npb(ji) - 1, jpi ) + 1 
    681                ij = ( npb(ji) - 1 ) / jpi + 1 
    682                t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  ) 
    683                zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer))) 
    684             END DO 
    685          END DO 
    686  
    687          DO layer  =  1, nlay_i 
    688             DO ji = kideb , kiut 
    689                ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt  
    690                t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp) 
    691                zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) 
     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) )      
     681         END DO 
     682 
     683         DO jk  =  1, nlay_s 
     684            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))) 
     687            END DO 
     688         END DO 
     689 
     690         DO jk  =  1, nlay_i 
     691            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))) 
    692695            END DO 
    693696         END DO 
     
    713716      !-------------------------------------------------------------------------! 
    714717      DO ji = kideb, kiut 
    715 #if ! defined key_coupled 
    716718         ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)  
    717          qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) ) 
    718 #endif 
     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) ) ) 
    719720         !                                ! surface ice conduction flux 
    720          isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  ) 
    721          fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   & 
    722             &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji)) 
     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)) 
    723724         !                                ! bottom ice conduction flux 
    724          fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 
    725       END DO 
    726  
    727       !-------------------------! 
    728       ! Heat conservation       ! 
    729       !-------------------------! 
    730       IF( con_i .AND. jiindex_1d > 0 ) THEN 
     725         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
     726      END DO 
     727 
     728      !----------------------------------------- 
     729      ! Heat flux used to warm/cool ice in W.m-2 
     730      !----------------------------------------- 
     731      DO ji = kideb, kiut 
     732         IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC 
     733            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & 
     734               &            ( 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 
     736            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    & 
     737               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
     738         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         ! 
    731755         DO ji = kideb, kiut 
    732             ! Upper snow value 
    733             fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )  
    734             ! Bott. snow value 
    735             fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )  
    736          END DO 
    737          DO ji = kideb, kiut         ! Upper ice layer 
    738             fc_i(ji,0) = - REAL( isnow(ji) ) * &  ! interface flux if there is snow 
    739                ( zkappa_i(ji,0)  * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 
    740                - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * &  
    741                zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 
    742          END DO 
    743          DO layer = 1, nlay_i - 1         ! Internal ice layers 
    744             DO ji = kideb, kiut 
    745                fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) 
    746                ii = MOD( npb(ji) - 1, jpi ) + 1 
    747                ij = ( npb(ji) - 1 ) / jpi + 1 
    748             END DO 
    749          END DO 
    750          DO ji = kideb, kiut         ! Bottom ice layers 
    751             fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 
    752          END DO 
     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         ! 
    753766      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    
    754774      ! 
     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 ) 
     783      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 
     784 
    755785   END SUBROUTINE lim_thd_dif 
     786 
     787   SUBROUTINE lim_thd_enmelt( kideb, kiut ) 
     788      !!----------------------------------------------------------------------- 
     789      !!                   ***  ROUTINE lim_thd_enmelt ***  
     790      !!                  
     791      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
     792      !! 
     793      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     794      !!------------------------------------------------------------------- 
     795      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
     796      ! 
     797      INTEGER  ::   ji, jk   ! dummy loop indices 
     798      REAL(wp) ::   ztmelts  ! local scalar  
     799      !!------------------------------------------------------------------- 
     800      ! 
     801      DO jk = 1, nlay_i             ! Sea ice energy of melting 
     802         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 )  )  
     808         END DO 
     809      END DO 
     810      DO jk = 1, nlay_s             ! Snow energy of melting 
     811         DO ji = kideb, kiut 
     812            q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 
     813         END DO 
     814      END DO 
     815      ! 
     816   END SUBROUTINE lim_thd_enmelt 
    756817 
    757818#else 
Note: See TracChangeset for help on using the changeset viewer.