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 8518 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90 – NEMO

Ignore:
Timestamp:
2017-09-13T18:46:56+02:00 (7 years ago)
Author:
clem
Message:

changes in style - part6 - commits of the day

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90

    r8514 r8518  
    9494      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C  
    9595      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    96       REAL(wp) ::   zhsu 
     96      REAL(wp) ::   z1_hsu 
    9797      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    9898      REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature  
     99      REAL(wp) ::   zcpi                      ! Ice specific heat 
     100      REAL(wp) ::   zhi                       !  
     101      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
    99102       
    100103      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
    101       REAL(wp), DIMENSION(jpij)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    102       REAL(wp), DIMENSION(jpij)     ::   ztsubit     ! surface temperature at previous iteration 
     104      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration 
    103105      REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness 
    104106      REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness 
     
    106108      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    107109      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function 
    108       REAL(wp), DIMENSION(jpij)     ::   dzf        ! derivative of the surface flux function 
     110      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    109111      REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature 
    110       REAL(wp), DIMENSION(jpij)     ::   zdifcase    ! case of the equation resolution (1->4) 
    111112      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
    112       REAL(wp), DIMENSION(jpij)     ::   zihic 
    113113       
     114      REAL(wp), DIMENSION(jpij,nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice 
     115      REAL(wp), DIMENSION(jpij,nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow 
     116      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Old temperature in the ice 
     117      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow 
     118      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     119      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    114120      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    115121      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    116122      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    117123      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztib        ! Old temperature in the ice 
    119124      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    120       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    121       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zspeche_i   ! Ice specific heat 
    122       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice 
    123125      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    124126      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    125127      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    126128      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    127       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    128       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztsb        ! Temporary temperature in the snow 
    129       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow 
    130129      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    131130      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    132131      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    133132      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    134       REAL(wp), DIMENSION(jpij)            ::   zdq, zq_ini, zhfx_err ! diag errors on heat 
     133      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    135134      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    136135       
     
    148147 
    149148      ! --- diag error on heat diffusion - PART 1 --- ! 
    150       zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
    151149      DO ji = 1, nidx 
    152150         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     
    167165      ! Ice / snow layers 
    168166      !-------------------- 
    169  
    170       z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer 
    171       z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer 
    172  
    173       DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    174          DO ji = 1 , nidx 
    175             z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 
    176          END DO 
    177       END DO 
    178  
    179       DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    180          DO ji = 1 , nidx 
    181             z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 
     167      z_s(1:nidx,1) = zh_s(1:nidx) 
     168      z_i(1:nidx,1) = zh_i(1:nidx) 
     169      DO jk = 2, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
     170         DO ji = 1 , nidx 
     171            z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 
     172         END DO 
     173      END DO 
     174 
     175      DO jk = 2, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
     176         DO ji = 1 , nidx 
     177            z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 
    182178         END DO 
    183179      END DO 
     
    187183      !------------------------------------------------------------------------------| 
    188184      ! 
    189       !------------------- 
    190       ! Computation of i0 
    191       !------------------- 
    192       ! i0 describes the fraction of solar radiation which does not contribute 
    193       ! to the surface energy budget but rather penetrates inside the ice. 
    194       ! We assume that no radiation is transmitted through the snow 
    195       ! If there is no no snow 
    196       ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
    197       ! zftrice = io.qsr_ice       is below the surface  
    198       ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    199       ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    200       zhsu = 0.1_wp ! threshold for the computation of i0 
     185      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    201186      DO ji = 1 , nidx 
    202          ! switches 
    203          isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  
    204          ! hs > 0, isnow = 1 
    205          zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )      
    206  
    207          i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    208       END DO 
    209  
    210       !------------------------------------------------------- 
    211       ! Solar radiation absorbed / transmitted at the surface 
    212       ! Derivative of the non solar flux 
    213       !------------------------------------------------------- 
    214       DO ji = 1 , nidx 
    215          zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
    216          zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
    217          dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
    218          zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value 
     187         !------------------- 
     188         ! Computation of i0 
     189         !------------------- 
     190         ! i0 describes the fraction of solar radiation which does not contribute 
     191         ! to the surface energy budget but rather penetrates inside the ice. 
     192         ! We assume that no radiation is transmitted through the snow 
     193         ! If there is no no snow 
     194         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
     195         ! zftrice = io.qsr_ice       is below the surface  
     196         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
     197         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
     198         zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )      
     199         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) ) 
     200 
     201         !------------------------------------------------------- 
     202         ! Solar radiation absorbed / transmitted at the surface 
     203         ! Derivative of the non solar flux 
     204         !------------------------------------------------------- 
     205         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
     206         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
     207         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
     208         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value 
    219209      END DO 
    220210 
     
    222212      ! Transmission - absorption of solar radiation in the ice 
    223213      !--------------------------------------------------------- 
    224  
    225       DO ji = 1, nidx           ! snow initialization 
    226          zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow 
    227       END DO 
    228  
    229       DO jk = 1, nlay_s          ! Radiation through snow 
     214      zradtr_s(1:nidx,0) = zftrice(1:nidx) 
     215      DO jk = 1, nlay_s 
    230216         DO ji = 1, nidx 
    231217            !                             ! radiation transmitted below the layer-th snow layer 
    232             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 
     218            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) ) 
    233219            !                             ! radiation absorbed by the layer-th snow layer 
    234220            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
    235221         END DO 
    236222      END DO 
    237  
    238       DO ji = 1, nidx           ! ice initialization 
    239          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 
    240       END DO 
    241  
    242       DO jk = 1, nlay_i          ! Radiation through ice 
     223          
     224      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) ) 
     225      DO jk = 1, nlay_i  
    243226         DO ji = 1, nidx 
    244227            !                             ! radiation transmitted below the layer-th ice layer 
    245             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 
     228            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) ) 
    246229            !                             ! radiation absorbed by the layer-th ice layer 
    247230            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    249232      END DO 
    250233 
    251       DO ji = 1, nidx           ! Radiation transmitted below the ice 
    252          ftr_ice_1d(ji) = zradtr_i(ji,nlay_i)  
    253       END DO 
     234      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice 
    254235 
    255236      !------------------------------------------------------------------------------| 
     
    257238      !------------------------------------------------------------------------------| 
    258239      ! 
    259       DO ji = 1, nidx        ! Old surface temperature 
    260          ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
    261          ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter 
    262          t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary 
    263          zdti   (ji) =  1000._wp                                 ! initial value of error 
    264       END DO 
    265  
    266       DO jk = 1, nlay_s       ! Old snow temperature 
    267          DO ji = 1 , nidx 
    268             ztsb(ji,jk) =  t_s_1d(ji,jk) 
    269          END DO 
    270       END DO 
    271  
    272       DO jk = 1, nlay_i       ! Old ice temperature 
    273          DO ji = 1 , nidx 
    274             ztib(ji,jk) =  t_i_1d(ji,jk) 
    275          END DO 
    276       END DO 
     240      ztsub  (1:nidx)   =  t_su_1d(1:nidx)                          ! temperature at the previous iter 
     241      t_su_1d(1:nidx)   =  MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary 
     242      zdti   (1:nidx)   =  1000._wp                                 ! initial value of error 
     243      ztsb   (1:nidx,:) =  t_s_1d(1:nidx,:)                         ! Old snow temperature 
     244      ztib   (1:nidx,:) =  t_i_1d(1:nidx,:)                         ! Old ice temperature 
    277245 
    278246      iconv    =  0           ! number of iterations 
     
    289257         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
    290258            DO ji = 1 , nidx 
    291                ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    292                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     259               ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 
    293260            END DO 
    294261            DO jk = 1, nlay_i-1 
    295262               DO ji = 1 , nidx 
    296                   ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
    297                      MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 
    298                   ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
     263                  ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
     264                     &                           MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ) 
    299265               END DO 
    300266            END DO 
     
    302268         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
    303269            DO ji = 1 , nidx 
    304                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
    305                   &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 )   
    306                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     270               ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
     271                  &                               - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 
    307272            END DO 
    308273            DO jk = 1, nlay_i-1 
    309274               DO ji = 1 , nidx 
    310                   ztcond_i(ji,jk) = rcdic +                                                                       &  
     275                  ztcond_i(ji,jk) = MAX( zkimin, rcdic +                                                           &  
    311276                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              & 
    312277                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   & 
    313                      &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )   
    314                   ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
     278                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  ) 
    315279               END DO 
    316280            END DO 
    317281            DO ji = 1 , nidx 
    318                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
    319                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )   
    320                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
     282               ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
     283                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) )   
    321284            END DO 
    322285         ENDIF 
     
    335298         SELECT CASE ( nn_monocat ) 
    336299 
    337          CASE (1,3) ! LIM3 
     300         CASE (1,3) 
    338301 
    339302            zepsilon = 0.1_wp 
     
    404367         DO jk = 1, nlay_i 
    405368            DO ji = 1 , nidx 
    406                ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    407                zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
    408                zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 
     369               ztitemp(ji,jk) = t_i_1d(ji,jk) 
     370               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
     371               zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 
    409372            END DO 
    410373         END DO 
     
    425388            DO ji = 1 , nidx 
    426389               ! update of the non solar flux according to the update in T_su 
    427                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 
     390               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    428391            END DO 
    429392         ENDIF 
     
    507470                  !  case 1 : no surface melting - snow present                                  | 
    508471                  !------------------------------------------------------------------------------| 
    509                   zdifcase(ji)    =  1.0 
    510472                  numeqmin(ji)    =  1 
    511473                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    513475                  !!surface equation 
    514476                  ztrid(ji,1,1)  = 0.0 
    515                   ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0) 
     477                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    516478                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    517                   zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji) 
     479                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji) 
    518480 
    519481                  !!first layer of snow equation 
     
    529491                  !------------------------------------------------------------------------------| 
    530492                  ! 
    531                   zdifcase(ji)    =  2.0 
    532493                  numeqmin(ji)    =  2 
    533494                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    552513                  !------------------------------------------------------------------------------| 
    553514                  ! 
    554                   zdifcase(ji)      =  3.0 
    555515                  numeqmin(ji)      =  nlay_s + 1 
    556516                  numeqmax(ji)      =  nlay_i + nlay_s + 1 
     
    558518                  !!surface equation    
    559519                  ztrid(ji,numeqmin(ji),1)   =  0.0 
    560                   ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
     520                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    561521                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    562                   zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
     522                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 
    563523 
    564524                  !!first layer of ice equation 
     
    572532                  IF ( nlay_i == 1 ) THEN 
    573533                     ztrid(ji,numeqmin(ji),1)    =  0.0 
    574                      ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0 
     534                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    575535                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0 
    576536                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     
    589549                  !------------------------------------------------------------------------------| 
    590550                  ! 
    591                   zdifcase(ji)    =  4.0 
    592551                  numeqmin(ji)    =  nlay_s + 2 
    593552                  numeqmax(ji)    =  nlay_i + nlay_s + 1 
     
    662621            ! surface temperature 
    663622            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 
    664             ztsubit(ji) = t_su_1d(ji) 
     623            ztsub(ji) = t_su_1d(ji) 
    665624            IF( t_su_1d(ji) < rt0 ) & 
    666625               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & 
     
    676635         DO ji = 1 , nidx 
    677636            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  ) 
    678             zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )      
     637            zdti   (ji) =  ABS( t_su_1d(ji) - ztsub(ji) )      
    679638         END DO 
    680639 
     
    754713 
    755714      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
     715      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2 
     716      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    756717      DO ji = 1, nidx 
    757          zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
    758             &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
     718         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     719            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
     720 
    759721         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    760             zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     722            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)  
    761723         ELSE                          ! case T_su = 0degC 
    762             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  
     724            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
    763725         ENDIF 
    764          hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     726         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
    765727 
    766728         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
    767          hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
    768       END DO  
    769  
    770       !----------------------------------------- 
    771       ! Heat flux used to warm/cool ice in W.m-2 
    772       !----------------------------------------- 
    773       DO ji = 1, nidx 
    774          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    775             hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & 
    776                &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    777          ELSE                          ! case T_su = 0degC 
    778             hfx_dif_1d(ji) = hfx_dif_1d(ji) +    & 
    779                &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    780          ENDIF 
    781          ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    782          hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 
     729         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
     730 
    783731      END DO    
    784732      ! 
Note: See TracChangeset for help on using the changeset viewer.