Changeset 5076


Ignore:
Timestamp:
2015-02-11T14:13:25+01:00 (6 years ago)
Author:
clem
Message:

LIM3: important bug fix on heat conservation

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5070 r5076  
    386386   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2) 
    387387   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2) 
     388   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_smv  !: transport of salt content 
    388389   ! 
    389390   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat_dhc !: snw/ice heat content variation   [W/m2]  
     
    482483      ! * Ice diagnostics 
    483484      ii = ii + 1 
    484       ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs  (jpi,jpj), diag_trp_ei(jpi,jpj),   &  
    485          &      diag_trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj),  STAT=ierr(ii) ) 
     485      ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei  (jpi,jpj),   &  
     486         &      diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat_dhc(jpi,jpj),  STAT=ierr(ii) ) 
    486487 
    487488      ice_alloc = MAXVAL( ierr(:) ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5070 r5076  
    223223            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
    224224            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
    225             IF ( ABS( zei    ) >  1    ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
     225            IF ( ABS( zei    ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
    226226            IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
    227227            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5070 r5076  
    281281         IF( nbpb > 0 ) THEN  ! If there is no ice, do nothing. 
    282282 
    283             !------------------------- 
    284             ! 4.1 Move to 1D arrays 
    285             !------------------------- 
    286  
    287             CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) ) 
    288             CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) ) 
    289             CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    290             CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    291  
    292             CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    293             CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    294             DO jk = 1, nlay_s 
    295                CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    296                CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    297             END DO 
    298             DO jk = 1, nlay_i 
    299                CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    300                CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    301                CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    302             END DO 
    303  
    304             CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) ) 
    305             CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    306             CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
    307             CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) ) 
    308             CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    309             CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    310             IF( .NOT. lk_cpl ) THEN 
    311                CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    312                CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    313             ENDIF 
    314             CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    315             CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
    316             CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) )  
    317             CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) ) 
    318             CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) ) 
    319             CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) ) 
    320  
    321             CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) ) 
    322             CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) ) 
    323  
    324             CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    325             CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    326             CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    327             CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) ) 
    328             CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) ) 
    329             CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    330  
    331             CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    332             CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    333             CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    334             CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) ) 
    335             CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
    336             CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
    337  
    338             CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
    339             CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    340             CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    341             CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    342             CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    343             CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) ) 
    344             CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) ) 
    345             CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) ) 
    346             CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) ) 
    347             CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) ) 
    348             CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) ) 
    349             CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 
    350  
    351             !-------------------------------- 
    352             ! 4.3) Thermodynamic processes 
    353             !-------------------------------- 
     283            !-------------------------! 
     284            ! --- Move to 1D arrays --- 
     285            !-------------------------! 
     286            CALL lim_thd_1d2d( nbpb, jl, 1 ) 
    354287 
    355288            !--------------------------------------! 
     
    383316            END IF 
    384317 
    385             !-------------------------------- 
    386             ! 4.4) Move 1D to 2D vectors 
    387             !-------------------------------- 
    388  
    389                CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj ) 
    390                CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj ) 
    391                CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj ) 
    392                CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj ) 
    393                CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj ) 
    394                CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj ) 
    395             DO jk = 1, nlay_s 
    396                CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj) 
    397                CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj) 
    398             END DO 
    399             DO jk = 1, nlay_i 
    400                CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj) 
    401                CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj) 
    402                CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj) 
    403             END DO 
    404                CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj ) 
    405  
    406                CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj ) 
    407                CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj ) 
    408  
    409                CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    410                CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    411                CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    412                CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj ) 
    413                CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj ) 
    414                CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    415  
    416                CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    417                CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    418                CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    419                CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj ) 
    420                CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
    421                CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
    422  
    423               CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
    424               CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    425               CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    426               CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    427               CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    428               CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj ) 
    429               CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj ) 
    430               CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj ) 
    431               CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj ) 
    432               CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj ) 
    433               CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    434               CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
    435             ! 
    436               CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    437               CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
     318            !-------------------------! 
     319            ! --- Move to 2D arrays --- 
     320            !-------------------------! 
     321            CALL lim_thd_1d2d( nbpb, jl, 2 ) 
     322 
    438323            ! 
    439324            IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 
     
    638523   END SUBROUTINE lim_thd_lam 
    639524 
     525   SUBROUTINE lim_thd_1d2d( nbpb, jl, kn ) 
     526      !!----------------------------------------------------------------------- 
     527      !!                   ***  ROUTINE lim_thd_1d2d ***  
     528      !!                  
     529      !! ** Purpose :   move arrays from 1d to 2d and the reverse 
     530      !!----------------------------------------------------------------------- 
     531      INTEGER, INTENT(in) ::   kn       ! 1= from 2D to 1D 
     532                                        ! 2= from 1D to 2D 
     533      INTEGER, INTENT(in) ::   nbpb     ! size of 1D arrays 
     534      INTEGER, INTENT(in) ::   jl       ! ice cat 
     535      INTEGER             ::   jk       ! dummy loop indices 
     536 
     537      SELECT CASE( kn ) 
     538 
     539      CASE( 1 ) 
     540 
     541         CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) ) 
     542         CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) ) 
     543         CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     544         CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     545          
     546         CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     547         CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     548         DO jk = 1, nlay_s 
     549            CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     550            CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     551         END DO 
     552         DO jk = 1, nlay_i 
     553            CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     554            CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     555            CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     556         END DO 
     557          
     558         CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) ) 
     559         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     560         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
     561         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) ) 
     562         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     563         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     564         IF( .NOT. lk_cpl ) THEN 
     565            CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     566            CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
     567         ENDIF 
     568         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
     569         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
     570         CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) )  
     571         CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) ) 
     572         CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) ) 
     573         CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) ) 
     574          
     575         CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) ) 
     576         CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) ) 
     577          
     578         CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     579         CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     580         CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     581         CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) ) 
     582         CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) ) 
     583         CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) ) 
     584          
     585         CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     586         CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     587         CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     588         CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) ) 
     589         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
     590         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
     591          
     592         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
     593         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
     594         CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     595         CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     596         CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     597         CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) ) 
     598         CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) ) 
     599         CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) ) 
     600         CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) ) 
     601         CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) ) 
     602         CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) ) 
     603         CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 
     604 
     605      CASE( 2 ) 
     606 
     607         CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj ) 
     608         CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj ) 
     609         CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj ) 
     610         CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj ) 
     611         CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj ) 
     612         CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj ) 
     613         DO jk = 1, nlay_s 
     614            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj) 
     615            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj) 
     616         END DO 
     617         DO jk = 1, nlay_i 
     618            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj) 
     619            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj) 
     620            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj) 
     621         END DO 
     622         CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj ) 
     623          
     624         CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj ) 
     625         CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj ) 
     626          
     627         CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     628         CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     629         CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     630         CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     631         CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj ) 
     632         CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj ) 
     633          
     634         CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     635         CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     636         CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     637         CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     638         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
     639         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
     640          
     641         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
     642         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 
     643         CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     644         CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     645         CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     646         CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj ) 
     647         CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj ) 
     648         CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj ) 
     649         CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj ) 
     650         CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj ) 
     651         CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
     652         CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
     653         ! 
     654         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
     655         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
     656                   
     657      END SELECT 
     658 
     659   END SUBROUTINE lim_thd_1d2d 
     660 
     661 
    640662   SUBROUTINE lim_thd_init 
    641663      !!----------------------------------------------------------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5070 r5076  
    102102      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
    103103      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
    104       INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    105104       
    106105      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
     
    115114      REAL(wp) ::   zhsu 
    116115       
     116      REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
    117117      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    118118      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration 
     
    166166      !!------------------------------------------------------------------      
    167167      !  
    168       CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
    169       CALL wrk_alloc( jpij, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     168      CALL wrk_alloc( jpij, numeqmin, numeqmax ) 
     169      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    170170      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    171171      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) 
     
    187187      !------------------------------------------------------------------------------! 
    188188      DO ji = kideb , kiut 
    189          ! is there snow or not 
    190          isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ) 
     189         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not 
    191190         ! layer thickness 
    192191         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 
     
    214213      ! 
    215214      !------------------------------------------------------------------------------| 
    216       ! 2) Radiation                                              | 
     215      ! 2) Radiation                          | 
    217216      !------------------------------------------------------------------------------| 
    218217      ! 
     
    231230      DO ji = kideb , kiut 
    232231         ! switches 
    233          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) ) )  
    234233         ! hs > 0, isnow = 1 
    235234         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )      
    236235 
    237          i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
     236         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    238237      END DO 
    239238 
     
    266265 
    267266      DO ji = kideb, kiut           ! ice initialization 
    268          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 
     267         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 
    269268      END DO 
    270269 
     
    290289         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
    291290         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter 
    292          t_su_1d   (ji) =  MIN( t_su_1d(ji), rtt - ztsu_err )    ! necessary 
    293          zerrit   (ji) =  1000._wp                               ! initial value of error 
     291         t_su_1d(ji) =  MIN( t_su_1d(ji), rtt - ztsu_err )       ! necessary 
     292         zerrit (ji) =  1000._wp                                 ! initial value of error 
    294293      END DO 
    295294 
     
    319318         IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula 
    320319            DO ji = kideb , kiut 
    321                ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 
    322                ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
     320               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rtt ) 
     321               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    323322            END DO 
    324323            DO jk = 1, nlay_i-1 
    325324               DO ji = kideb , kiut 
    326                   ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
    327                      MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 
    328                   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 * rtt) 
     327                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
    329328               END DO 
    330329            END DO 
     
    333332         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    334333            DO ji = kideb , kiut 
    335                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt )   & 
     334               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rtt )   & 
    336335                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt )   
    337336               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     
    339338            DO jk = 1, nlay_i-1 
    340339               DO ji = kideb , kiut 
    341                   ztcond_i(ji,jk) = rcdic +                                                                     &  
    342                      &                 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                          & 
    343                      &                 / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)   & 
    344                      &               - 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 * rtt )   & 
     343                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rtt )   
    345344                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 
    346345               END DO 
    347346            END DO 
    348347            DO ji = kideb , kiut 
    349                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,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) - rtt )   & 
    350349                  &                        - 0.011_wp * ( t_bo_1d(ji) - rtt )   
    351350               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
     
    368367         CASE (1,3) ! LIM3 
    369368 
    370             zepsilon = 0.1 
    371             zh_thres = EXP( 1._wp ) * zepsilon / 2. 
     369            zepsilon = 0.1_wp 
     370            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 
    372371 
    373372            DO ji = kideb, kiut 
    374373    
    375374               ! Mean sea ice thermal conductivity 
    376                zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp) 
     375               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp ) 
    377376 
    378377               ! Effective thickness he (zhe) 
     
    384383               ! G(he) 
    385384               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
    386                zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2.*zhe / zepsilon ) ) 
     385               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 
    387386    
    388387               ! Impose G(he) < 2. 
     
    407406         DO jk = 1, nlay_s-1 
    408407            DO ji = kideb , kiut 
    409                zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) ) 
     408               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
    410409            END DO 
    411410         END DO 
     
    414413         DO jk = 1, nlay_i-1 
    415414            DO ji = kideb , kiut 
    416                zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) ) 
     415               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 
    417416            END DO 
    418417         END DO 
     
    425424            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
    426425           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
    427             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp ) 
     426            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    428427         END DO 
    429428 
     
    436435            DO ji = kideb , kiut 
    437436               ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    438                zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 ) 
    439                zeta_i(ji,jk)    = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 ) 
     437               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rtt ) * ( ztib(ji,jk) - rtt ), epsi10 ) 
     438               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 
    440439            END DO 
    441440         END DO 
     
    444443            DO ji = kideb , kiut 
    445444               ztstemp(ji,jk) = t_s_1d(ji,jk) 
    446                zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
     445               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 
    447446            END DO 
    448447         END DO 
     
    463462         DO ji = kideb , kiut 
    464463            ! update incoming flux 
    465             zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation 
    466                + qns_ice_1d(ji)                   ! non solar total flux  
    467             ! (LWup, LWdw, SH, LH) 
     464            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation 
     465               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH) 
    468466         END DO 
    469467 
     
    493491         DO numeq = nlay_s + 2, nlay_s + nlay_i  
    494492            DO ji = kideb , kiut 
    495                jk              = numeq - nlay_s - 1 
    496                ztrid(ji,numeq,1)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 
    497                ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 
    498                   zkappa_i(ji,jk)) 
    499                ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk) 
    500                zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
    501                   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) 
    502498            END DO 
    503499         ENDDO 
     
    507503            !!ice bottom term 
    508504            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    509             ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & 
    510                +  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) ) 
    511506            ztrid(ji,numeq,3)  =  0.0 
    512             zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
    513                ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 
    514                *  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) )  
    515509         ENDDO 
    516510 
     
    525519               !!snow interior terms (bottom equation has the same form as the others) 
    526520               DO numeq = 3, nlay_s + 1 
    527                   jk =  numeq - 1 
    528                   ztrid(ji,numeq,1)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 
    529                   ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 
    530                      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) ) 
    531524                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    532                   zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
    533                      zradab_s(ji,jk) 
     525                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    534526               END DO 
    535527 
     
    537529               IF ( nlay_i.eq.1 ) THEN 
    538530                  ztrid(ji,nlay_s+2,3)    =  0.0 
    539                   zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
    540                      t_bo_1d(ji)  
     531                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    541532               ENDIF 
    542533 
     
    551542 
    552543                  !!surface equation 
    553                   ztrid(ji,1,1) = 0.0 
    554                   ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 
    555                   ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 
    556                   zindterm(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) 
    557548 
    558549                  !!first layer of snow equation 
    559                   ztrid(ji,2,1)  =  - zkappa_s(ji,0)*zg1s*zeta_s(ji,1) 
    560                   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 ) 
    561552                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    562                   zindterm(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) 
    563554 
    564555               ELSE  
     
    574565                  !!first layer of snow equation 
    575566                  ztrid(ji,2,1)  =  0.0 
    576                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & 
    577                      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 ) 
    578568                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    579                   zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
    580                      ( zradab_s(ji,1) +                         & 
    581                      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) )  
    582571               ENDIF 
    583572            ELSE 
     
    605594                  !!first layer of ice equation 
    606595                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    607                   ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &  
    608                      + zkappa_i(ji,0) * zg1 ) 
    609                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)   
    610                   zindterm(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)   
    611599 
    612600                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    614602                  IF (nlay_i.eq.1) THEN 
    615603                     ztrid(ji,numeqmin(ji),1)    =  0.0 
    616                      ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0)*2.0 
    617                      ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0)*2.0 
    618                      ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0)*2.0*zeta_i(ji,1) 
    619                      ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 
    620                         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) ) 
    621608                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    622609 
    623                      zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
    624                         ( 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) ) 
    625612                  ENDIF 
    626613 
     
    638625                  !!first layer of ice equation 
    639626                  ztrid(ji,numeqmin(ji),1)      =  0.0 
    640                   ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & 
    641                      zg1)   
     627                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    642628                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    643                   zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
    644                      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) )  
    645631 
    646632                  !!case of only one layer in the ice (surface & ice equations are altered) 
    647633                  IF (nlay_i.eq.1) THEN 
    648634                     ztrid(ji,numeqmin(ji),1)  =  0.0 
    649                      ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 
    650                         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) ) 
    651636                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    652                      zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
    653                         (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 
    654                         + 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 
    655639                  ENDIF 
    656640 
     
    683667            DO ji = kideb , kiut 
    684668               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 
    685                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 
    686                   ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 
    687                zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 
    688                   zindtbis(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) 
    689671            END DO 
    690672         END DO 
     
    692674         DO ji = kideb , kiut 
    693675            ! ice temperatures 
    694             t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
     676            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 
    695677         END DO 
    696678 
     
    698680            DO ji = kideb , kiut 
    699681               jk    =  numeq - nlay_s - 1 
    700                t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 
    701                   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) 
    702683            END DO 
    703684         END DO 
     
    706687            ! snow temperatures       
    707688            IF (ht_s_1d(ji) > 0._wp) & 
    708                t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    709                *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 
    710                *        MAX(0.0,SIGN(1.0,ht_s_1d(ji)))  
     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) ) )  
    711691 
    712692            ! surface temperature 
    713             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) ) ) 
    714694            ztsubit(ji) = t_su_1d(ji) 
    715695            IF( t_su_1d(ji) < rtt ) & 
    716                t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
    717                &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
     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))   
    718698         END DO 
    719699         ! 
     
    726706         DO ji = kideb , kiut 
    727707            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rtt ) , 190._wp  ) 
    728             zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )      
     708            zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) )      
    729709         END DO 
    730710 
     
    732712            DO ji = kideb , kiut 
    733713               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  ) 
    734                zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 
     714               zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 
    735715            END DO 
    736716         END DO 
     
    738718         DO jk  =  1, nlay_i 
    739719            DO ji = kideb , kiut 
    740                ztmelt_i        = -tmut * s_i_1d(ji,jk) + rtt  
    741                t_i_1d(ji,jk) =  MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 
    742                zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 
     720               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rtt  
     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) ) ) 
    743723            END DO 
    744724         END DO 
     
    767747         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) ) ) 
    768748         !                                ! surface ice conduction flux 
    769          isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 
    770          fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    771             &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
     749         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 
     750         fc_su(ji)       =  -           isnow(ji)  * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
     751            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    772752         !                                ! bottom ice conduction flux 
    773753         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    774754      END DO 
     755 
     756      ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
     757      CALL lim_thd_enmelt( kideb, kiut ) 
     758 
     759 
     760      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
     761      DO ji = kideb, kiut 
     762         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
     763            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
     764         IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC 
     765            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
     766         ELSE                          ! case T_su = 0degC 
     767            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  
     768         ENDIF 
     769         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     770      END DO  
     771 
     772      ! diagnose external surface (forced case) or bottom (forced case) from heat conservation 
     773      IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed 
     774         ! 
     775         DO ji = kideb, kiut 
     776            qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji) 
     777            fc_su     (ji) = fc_su(ji)      - zhfx_err(ji) 
     778         END DO 
     779         ! 
     780      ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed 
     781         ! 
     782         DO ji = kideb, kiut 
     783            fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji) 
     784         END DO 
     785         ! 
     786      ENDIF 
    775787 
    776788      !----------------------------------------- 
     
    781793            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & 
    782794               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    783          ELSE                         ! case T_su = 0degC 
     795         ELSE                          ! case T_su = 0degC 
    784796            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    & 
    785797               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
     
    787799      END DO 
    788800 
    789       ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
    790       CALL lim_thd_enmelt( kideb, kiut ) 
    791  
    792       ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    793       DO ji = kideb, kiut 
    794          zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
    795             &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
    796          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 )  
    797          hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
    798       END DO  
    799  
    800       ! diagnose external surface (forced case) or bottom (forced case) from heat conservation 
    801       IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed 
    802          ! 
    803          DO ji = kideb, kiut 
    804             qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji) 
    805             fc_su     (ji) = fc_su(ji)      - zhfx_err(ji) 
    806          END DO 
    807          ! 
    808       ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed 
    809          ! 
    810          DO ji = kideb, kiut 
    811             fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji) 
    812          END DO 
    813          ! 
    814       ENDIF 
    815  
    816       ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2) 
     801      ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m-2) 
    817802      DO ji = kideb, kiut 
    818803         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     
    821806    
    822807      ! 
    823       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    824       CALL wrk_dealloc( jpij, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     808      CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 
     809      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    825810      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    826811      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 
     
    850835         DO ji = kideb, kiut 
    851836            ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt  
    852             rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 
    853             q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                             & 
    854                &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) )   & 
    855                &                   - rcp  *                 ( ztmelts-rtt )  )  
     837            rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi20 ) ) 
     838            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                                & 
     839               &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk) - rtt, -epsi20 ) )   & 
     840               &                   - rcp  *                   ( ztmelts-rtt )  )  
    856841         END DO 
    857842      END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5070 r5076  
    7272      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ow 
    7373      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
    74       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold           ! old ice volume... 
     74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume... 
    7575      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness 
    7676      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies 
     
    8484      CALL wrk_alloc( jpi,jpj,1,         zs0ow ) 
    8585      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
    86       CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold ) 
     86      CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold, zsmvold ) 
    8787 
    8888      IF( numit == nstart .AND. lwp ) THEN 
     
    105105 
    106106         ! mass and salt flux init 
    107          zviold(:,:,:) = v_i(:,:,:) 
    108          zvsold(:,:,:) = v_s(:,:,:) 
    109          zeiold(:,:)   = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
    110          zesold(:,:)   = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
     107         zviold(:,:,:)  = v_i(:,:,:) 
     108         zvsold(:,:,:)  = v_s(:,:,:) 
     109         zsmvold(:,:,:) = smv_i(:,:,:) 
     110         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
     111         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
    111112 
    112113         !--- Thickness correction init. ------------------------------- 
     
    412413               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 
    413414 
    414                diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 
    415                diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 
     415               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice 
     416               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice 
     417               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 
    416418            END DO 
    417419         END DO 
     
    459461      CALL wrk_dealloc( jpi,jpj,1,         zs0ow ) 
    460462      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
    461       CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax ) 
     463      CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax, zsmvold ) 
    462464      ! 
    463465      IF( nn_timing == 1 )  CALL timing_stop('limtrp') 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5067 r5076  
    147147      DO jj = 1, jpj 
    148148         DO ji = 1, jpi             
    149             diag_heat_dhc(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    150                &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    151                &                   ) * r1_rdtice    
     149            diag_heat_dhc(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
     150               &                       SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
     151               &                     ) * r1_rdtice    
    152152         END DO 
    153153      END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5067 r5076  
    188188      DO jj = 1, jpj 
    189189         DO ji = 1, jpi             
    190             diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj)+   & 
     190            diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj) -  & 
    191191               &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    192192               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5067 r5076  
    185185      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport 
    186186      CALL iom_put( "snwtrp"      , diag_trp_vs * rday  )        ! snw volume transport 
     187      CALL iom_put( "saltrp"      , diag_trp_smv * rday * rhoic ) ! salt content transport 
    187188      CALL iom_put( "deitrp"      , diag_trp_ei         )        ! advected ice enthalpy (W/m2) 
    188189      CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2) 
Note: See TracChangeset for help on using the changeset viewer.