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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5575 for branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 – NEMO

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

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

File:
1 edited

Legend:

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

    r5473 r5575  
    2222   USE phycst         ! physical constants 
    2323   USE dom_oce        ! ocean space and time domain variables 
    24    USE oce     , ONLY : fraqsr_1lev  
    2524   USE ice            ! LIM: sea-ice variables 
    26    USE par_ice        ! LIM: sea-ice parameters 
    2725   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2826   USE sbc_ice        ! Surface boundary condition: ice fields 
    2927   USE thd_ice        ! LIM thermodynamic sea-ice variables 
    3028   USE dom_ice        ! LIM sea-ice domain 
    31    USE domvvl         ! domain: variable volume level 
    3229   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion 
    3330   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation 
    3431   USE limthd_sal     ! LIM: thermodynamics, ice salinity 
    3532   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution 
     33   USE limthd_lac     ! LIM-3 lateral accretion 
     34   USE limitd_th      ! remapping thickness distribution 
    3635   USE limtab         ! LIM: 1D <==> 2D transformation 
    3736   USE limvar         ! LIM: sea-ice variables 
     
    4443   USE timing         ! Timing 
    4544   USE limcons        ! conservation tests 
     45   USE limctl 
    4646 
    4747   IMPLICIT NONE 
    4848   PRIVATE 
    4949 
    50    PUBLIC   lim_thd        ! called by limstp module 
    51    PUBLIC   lim_thd_init   ! called by iceini module 
     50   PUBLIC   lim_thd         ! called by limstp module 
     51   PUBLIC   lim_thd_init    ! called by sbc_lim_init 
    5252 
    5353   !! * Substitutions 
     
    8080      !! ** References :  
    8181      !!--------------------------------------------------------------------- 
    82       INTEGER, INTENT(in) ::   kt    ! number of iteration 
     82      INTEGER, INTENT(in) :: kt    ! number of iteration 
    8383      !! 
    8484      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    85       INTEGER  :: nbpb             ! nb of icy pts for thermo. cal. 
     85      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations 
    8686      INTEGER  :: ii, ij           ! temporary dummy loop index 
    87       REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
    88       REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient 
    89       REAL(wp) :: zareamin  
    9087      REAL(wp) :: zfric_u, zqld, zqfr 
    91       ! 
    9288      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    93       ! 
    94       REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
     89      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
     90      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
     91      ! 
    9592      !!------------------------------------------------------------------- 
    96       CALL wrk_alloc( jpi, jpj, zqsr, zqns ) 
    9793 
    9894      IF( nn_timing == 1 )  CALL timing_start('limthd') 
     
    10197      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    10298 
     99      CALL lim_var_glo2eqv 
    103100      !------------------------------------------------------------------------! 
    104101      ! 1) Initialization of some variables                                    ! 
     
    106103      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice 
    107104 
    108  
    109105      !-------------------- 
    110106      ! 1.2) Heat content     
    111107      !-------------------- 
    112       ! Change the units of heat content; from global units to J.m3 
     108      ! Change the units of heat content; from J/m2 to J/m3 
    113109      DO jl = 1, jpl 
    114110         DO jk = 1, nlay_i 
     
    116112               DO ji = 1, jpi 
    117113                  !0 if no ice and 1 if yes 
    118                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  ) 
     114                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  ) 
    119115                  !Energy of melting q(S,T) [J.m-3] 
    120                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 
    121                   !convert units ! very important that this line is here         
    122                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac  
     116                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 
    123117               END DO 
    124118            END DO 
     
    128122               DO ji = 1, jpi 
    129123                  !0 if no ice and 1 if yes 
    130                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  ) 
     124                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  ) 
    131125                  !Energy of melting q(S,T) [J.m-3] 
    132                   e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 
    133                   !convert units ! very important that this line is here 
    134                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac  
     126                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 
    135127               END DO 
    136128            END DO 
     
    140132      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      ! 
    141133      !-----------------------------------------------------------------------------! 
    142  
    143       !--- Ocean solar and non solar fluxes to be used in zqld 
    144       IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean 
    145          ! 
    146          zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:) 
    147          ! 
    148       ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted 
    149          ! 
    150          zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:) 
    151          ! 
    152          DO jl = 1, jpl 
    153             DO jj = 1, jpj 
    154                DO ji = 1, jpi 
    155                   zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    156                   zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    157                END DO 
    158             END DO 
    159          END DO 
    160          ! 
    161       ENDIF 
    162  
    163 !CDIR NOVERRCHK 
    164134      DO jj = 1, jpj 
    165 !CDIR NOVERRCHK 
    166135         DO ji = 1, jpi 
    167             rswitch          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
     136            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    168137            ! 
    169138            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
     
    173142            !           !  temperature and turbulent mixing (McPhee, 1992) 
    174143            ! 
    175  
    176144            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 
    177             ! REMARK valid at least in forced mode from clem 
    178             ! precip is included in qns but not in qns_ice 
    179             IF ( lk_cpl ) THEN 
    180                zqld =  tms(ji,jj) * rdt_ice *  & 
    181                   &    (   zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj)               &   ! pfrld already included in coupled mode 
    182                   &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *     &   ! heat content of precip 
    183                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )   & 
    184                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 
    185             ELSE 
    186                zqld =  tms(ji,jj) * rdt_ice *  & 
    187                   &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) )    & 
    188                   &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *             &  ! heat content of precip 
    189                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )           & 
    190                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 
    191             ENDIF 
    192  
    193             !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
    194             zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
     145            zqld =  tmask(ji,jj,1) * rdt_ice *  & 
     146               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
     147 
     148            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
     149            zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
     150 
     151            ! --- Energy from the turbulent oceanic heat flux (W/m2) --- ! 
     152            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin )  
     153            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 
     154            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
     155            ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
     156            !                        the freezing point, so that we do not have SST < T_freeze 
     157            !                        This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 
    195158 
    196159            !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice 
    197             qlead(ji,jj) = MIN( 0._wp , zqld - zqfr )  
     160            qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    198161 
    199162            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    200             IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN 
    201                fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 
     163            IF( zqld > 0._wp ) THEN 
     164               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 
    202165               qlead(ji,jj) = 0._wp 
    203166            ELSE 
     
    205168            ENDIF 
    206169            ! 
    207             !-- Energy from the turbulent oceanic heat flux --- ! 
    208             !clem zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 
    209             zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin )  
    210             fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2  
    211             ! upper bound for fhtur: we do not want SST to drop below Tfreeze.  
    212             ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr)    
    213             ! This is not a clean budget, so that should be corrected at some point 
    214             fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
    215  
    216170            ! ----------------------------------------- 
    217171            ! Net heat flux on top of ice-ocean [W.m-2] 
    218172            ! ----------------------------------------- 
    219             !     First  step here      : heat flux at the ocean surface + precip 
    220             !     Second step below     : heat flux at the ice   surface (after limthd_dif)  
    221             hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         &  
    222                ! heat flux above the ocean 
    223                &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  & 
    224                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    225                &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  & 
    226                &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) 
     173            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
    227174 
    228175            ! ----------------------------------------------------------------------------- 
    229             ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2] 
     176            ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    230177            ! ----------------------------------------------------------------------------- 
    231178            !     First  step here              :  non solar + precip - qlead - qturb 
    232179            !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema)  
    233180            !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar 
    234             hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       &  
    235                ! Non solar heat flux received by the ocean 
    236                &    +        pfrld(ji,jj) * qns(ji,jj)                                                                            & 
    237                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    238                &    +      ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj)       & 
    239                &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  & 
    240                &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt )       & 
    241                ! heat flux taken from the ocean where there is open water ice formation 
    242                &    -      qlead(ji,jj) * r1_rdtice                                                                               & 
    243                ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 
    244                &    -      at_i(ji,jj) * fhtur(ji,jj)                                                                             & 
    245                &    -      at_i(ji,jj) *  fhld(ji,jj) 
    246  
     181            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean                
     182               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation 
     183               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence 
     184               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt  
     185                                                                                   !    (fhld should be 0 while bott growth) 
    247186         END DO 
    248187      END DO 
     
    259198         ENDIF 
    260199 
    261          zareamin = epsi10 
    262200         nbpb = 0 
    263201         DO jj = 1, jpj 
    264202            DO ji = 1, jpi 
    265                IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN      
     203               IF ( a_i(ji,jj,jl) > epsi10 ) THEN      
    266204                  nbpb      = nbpb  + 1 
    267205                  npb(nbpb) = (jj - 1) * jpi + ji 
     
    272210         ! debug point to follow 
    273211         jiindex_1d = 0 
    274          IF( ln_nicep ) THEN 
    275             DO ji = mi0(jiindx), mi1(jiindx) 
    276                DO jj = mj0(jjindx), mj1(jjindx) 
     212         IF( ln_icectl ) THEN 
     213            DO ji = mi0(iiceprt), mi1(iiceprt) 
     214               DO jj = mj0(jiceprt), mj1(jiceprt) 
    277215                  jiindex_1d = (jj - 1) * jpi + ji 
    278216                  WRITE(numout,*) ' lim_thd : Category no : ', jl  
     
    289227         IF( nbpb > 0 ) THEN  ! If there is no ice, do nothing. 
    290228 
    291             !------------------------- 
    292             ! 4.1 Move to 1D arrays 
    293             !------------------------- 
    294  
    295             CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) ) 
    296             CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) ) 
    297             CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    298             CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    299  
    300             CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    301             CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    302             DO jk = 1, nlay_s 
    303                CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    304                CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    305             END DO 
    306             DO jk = 1, nlay_i 
    307                CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    308                CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    309                CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
    310             END DO 
    311  
    312             CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) ) 
    313             CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    314             CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
    315             CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) ) 
    316             CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    317             CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    318             IF( .NOT. lk_cpl ) THEN 
    319                CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    320                CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    321             ENDIF 
    322             CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    323             CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
    324             CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) )  
    325             CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) ) 
    326             CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) ) 
    327             CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) ) 
    328  
    329             CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) ) 
    330             CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) ) 
    331  
    332             CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    333             CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    334             CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    335             CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) ) 
    336             CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) ) 
    337             CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    338  
    339             CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    340             CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    341             CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    342             CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) ) 
    343             CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
    344             CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
    345  
    346             CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
    347             CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    348             CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    349             CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) ) 
    350             CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) ) 
    351             CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) ) 
    352             CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) ) 
    353             CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) ) 
    354             CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) ) 
    355             CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) ) 
    356             CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) ) 
    357             CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 
    358  
    359             !-------------------------------- 
    360             ! 4.3) Thermodynamic processes 
    361             !-------------------------------- 
     229            !-------------------------! 
     230            ! --- Move to 1D arrays --- 
     231            !-------------------------! 
     232            CALL lim_thd_1d2d( nbpb, jl, 1 ) 
     233 
     234            !--------------------------------------! 
     235            ! --- Ice/Snow Temperature profile --- ! 
     236            !--------------------------------------! 
     237            CALL lim_thd_dif( 1, nbpb ) 
    362238 
    363239            !---------------------------------! 
    364             ! Ice/Snow Temperature profile    ! 
    365             !---------------------------------! 
    366             CALL lim_thd_dif( 1, nbpb ) 
    367  
    368             !---------------------------------! 
    369             ! Ice/Snow thicnkess              ! 
     240            ! --- Ice/Snow thickness ---      ! 
    370241            !---------------------------------! 
    371242            CALL lim_thd_dh( 1, nbpb )     
     
    375246                                             
    376247            !---------------------------------! 
    377             ! --- Ice salinity --- ! 
     248            ! --- Ice salinity ---            ! 
    378249            !---------------------------------! 
    379250            CALL lim_thd_sal( 1, nbpb )     
    380251 
    381252            !---------------------------------! 
    382             ! --- temperature update --- ! 
     253            ! --- temperature update ---      ! 
    383254            !---------------------------------! 
    384255            CALL lim_thd_temp( 1, nbpb ) 
    385256 
    386             !-------------------------------- 
    387             ! 4.4) Move 1D to 2D vectors 
    388             !-------------------------------- 
    389  
    390                CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj ) 
    391                CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj ) 
    392                CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj ) 
    393                CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj ) 
    394                CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj ) 
    395                CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj ) 
    396             DO jk = 1, nlay_s 
    397                CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj) 
    398                CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj) 
    399             END DO 
    400             DO jk = 1, nlay_i 
    401                CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj) 
    402                CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj) 
    403                CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj) 
    404             END DO 
    405                CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj ) 
    406  
    407                CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj ) 
    408                CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj ) 
    409  
    410                CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    411                CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    412                CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    413                CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj ) 
    414                CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj ) 
    415                CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    416  
    417                CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    418                CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    419                CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    420                CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj ) 
    421                CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
    422                CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
    423  
    424               CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
    425               CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    426               CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    427               CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj ) 
    428               CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj ) 
    429               CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj ) 
    430               CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj ) 
    431               CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj ) 
    432               CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj ) 
    433               CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj ) 
    434               CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    435               CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb)   , jpi, jpj ) 
    436             ! 
    437               CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    438               CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
     257            !------------------------------------! 
     258            ! --- lateral melting if monocat --- ! 
     259            !------------------------------------! 
     260            IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 
     261               CALL lim_thd_lam( 1, nbpb ) 
     262            END IF 
     263 
     264            !-------------------------! 
     265            ! --- Move to 2D arrays --- 
     266            !-------------------------! 
     267            CALL lim_thd_1d2d( nbpb, jl, 2 ) 
     268 
    439269            ! 
    440270            IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 
    441271         ENDIF 
    442272         ! 
    443       END DO 
     273      END DO !jl 
    444274 
    445275      !------------------------------------------------------------------------------! 
     
    448278 
    449279      !------------------------ 
    450       ! 5.1) Ice heat content               
     280      ! Ice heat content               
    451281      !------------------------ 
    452       ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     282      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 
    453283      DO jl = 1, jpl 
    454284         DO jk = 1, nlay_i 
    455             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) ) 
     285            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i 
    456286         END DO 
    457287      END DO 
    458288 
    459289      !------------------------ 
    460       ! 5.2) Snow heat content               
     290      ! Snow heat content               
    461291      !------------------------ 
    462       ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     292      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 
    463293      DO jl = 1, jpl 
    464294         DO jk = 1, nlay_s 
    465             e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) ) 
     295            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s 
    466296         END DO 
    467297      END DO 
    468  
     298  
    469299      !---------------------------------- 
    470       ! 5.3) Change thickness to volume 
     300      ! Change thickness to volume 
    471301      !---------------------------------- 
    472       CALL lim_var_eqv2glo 
     302      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:) 
     303      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
     304      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
     305 
     306      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat) 
     307      DO jl  = 1, jpl 
     308         DO jj = 1, jpj 
     309            DO ji = 1, jpi 
     310               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) ) 
     311               oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) 
     312            END DO 
     313         END DO 
     314      END DO 
     315 
     316      CALL lim_var_zapsmall 
    473317 
    474318      !-------------------------------------------- 
    475       ! 5.4) Diagnostic thermodynamic growth rates 
     319      ! Diagnostic thermodynamic growth rates 
    476320      !-------------------------------------------- 
     321      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print 
     322 
    477323      IF(ln_ctl) THEN            ! Control print 
    478324         CALL prt_ctl_info(' ') 
    479325         CALL prt_ctl_info(' - Cell values : ') 
    480326         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    481          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd  : cell area :') 
     327         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd  : cell area :') 
    482328         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :') 
    483329         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :') 
     
    508354      ! 
    509355      ! 
    510       CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    511  
    512       ! 
    513       ! conservation test 
    514356      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     357 
     358      !------------------------------------------------------------------------------| 
     359      !  6) Transport of ice between thickness categories.                           | 
     360      !------------------------------------------------------------------------------| 
     361      ! Given thermodynamic growth rates, transport ice between thickness categories. 
     362      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     363 
     364      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt ) 
     365 
     366      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     367 
     368      !------------------------------------------------------------------------------| 
     369      !  7) Add frazil ice growing in leads. 
     370      !------------------------------------------------------------------------------| 
     371      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     372 
     373      CALL lim_thd_lac 
     374       
     375      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     376 
     377      ! Control print 
     378      IF(ln_ctl) THEN 
     379         CALL lim_var_glo2eqv 
     380 
     381         CALL prt_ctl_info(' ') 
     382         CALL prt_ctl_info(' - Cell values : ') 
     383         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
     384         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th  : cell area :') 
     385         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
     386         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
     387         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :') 
     388         DO jl = 1, jpl 
     389            CALL prt_ctl_info(' ') 
     390            CALL prt_ctl_info(' - Category : ', ivar1=jl) 
     391            CALL prt_ctl_info('   ~~~~~~~~~~') 
     392            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ') 
     393            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ') 
     394            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ') 
     395            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ') 
     396            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ') 
     397            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ') 
     398            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ') 
     399            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ') 
     400            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
     401            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
     402            DO jk = 1, nlay_i 
     403               CALL prt_ctl_info(' ') 
     404               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
     405               CALL prt_ctl_info('   ~~~~~~~') 
     406               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
     407               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
     408            END DO 
     409         END DO 
     410      ENDIF 
    515411      ! 
    516412      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
     
    518414   END SUBROUTINE lim_thd  
    519415 
     416  
    520417   SUBROUTINE lim_thd_temp( kideb, kiut ) 
    521418      !!----------------------------------------------------------------------- 
     
    534431      DO jk = 1, nlay_i 
    535432         DO ji = kideb, kiut 
    536             ztmelts       =  -tmut * s_i_1d(ji,jk) + rtt 
     433            ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0 
    537434            ! Conversion q(S,T) -> T (second order equation) 
    538435            zaaa          =  cpic 
    539             zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 
    540             zccc          =  lfus * ( ztmelts - rtt ) 
     436            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus 
     437            zccc          =  lfus * ( ztmelts - rt0 ) 
    541438            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 
    542             t_i_1d(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 
     439            t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 
    543440             
    544441            ! mask temperature 
    545442            rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    546             t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rtt 
     443            t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 
    547444         END DO  
    548445      END DO  
    549446 
    550447   END SUBROUTINE lim_thd_temp 
     448 
     449   SUBROUTINE lim_thd_lam( kideb, kiut ) 
     450      !!----------------------------------------------------------------------- 
     451      !!                   ***  ROUTINE lim_thd_lam ***  
     452      !!                  
     453      !! ** Purpose :   Lateral melting in case monocategory 
     454      !!                          ( dA = A/2h dh ) 
     455      !!----------------------------------------------------------------------- 
     456      INTEGER, INTENT(in) ::   kideb, kiut        ! bounds for the spatial loop 
     457      INTEGER             ::   ji                 ! dummy loop indices 
     458      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo 
     459      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting 
     460      REAL(wp)            ::   zvi, zvs           ! ice/snow volumes  
     461 
     462      DO ji = kideb, kiut 
     463         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     464         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
     465            zvi          = a_i_1d(ji) * ht_i_1d(ji) 
     466            zvs          = a_i_1d(ji) * ht_s_1d(ji) 
     467            ! lateral melting = concentration change 
     468            zhi_bef     = ht_i_1d(ji) - zdh_mel 
     469            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 
     470            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
     471            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel )  
     472             ! adjust thickness 
     473            ht_i_1d(ji) = zvi / a_i_1d(ji)             
     474            ht_s_1d(ji) = zvs / a_i_1d(ji)             
     475            ! retrieve total concentration 
     476            at_i_1d(ji) = a_i_1d(ji) 
     477         END IF 
     478      END DO 
     479       
     480   END SUBROUTINE lim_thd_lam 
     481 
     482   SUBROUTINE lim_thd_1d2d( nbpb, jl, kn ) 
     483      !!----------------------------------------------------------------------- 
     484      !!                   ***  ROUTINE lim_thd_1d2d ***  
     485      !!                  
     486      !! ** Purpose :   move arrays from 1d to 2d and the reverse 
     487      !!----------------------------------------------------------------------- 
     488      INTEGER, INTENT(in) ::   kn       ! 1= from 2D to 1D 
     489                                        ! 2= from 1D to 2D 
     490      INTEGER, INTENT(in) ::   nbpb     ! size of 1D arrays 
     491      INTEGER, INTENT(in) ::   jl       ! ice cat 
     492      INTEGER             ::   jk       ! dummy loop indices 
     493 
     494      SELECT CASE( kn ) 
     495 
     496      CASE( 1 ) 
     497 
     498         CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) ) 
     499         CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) ) 
     500         CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     501         CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     502          
     503         CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     504         CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
     505         DO jk = 1, nlay_s 
     506            CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     507            CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     508         END DO 
     509         DO jk = 1, nlay_i 
     510            CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     511            CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     512            CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) ) 
     513         END DO 
     514          
     515         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 
     516         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     517         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
     518         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) ) 
     519         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     520         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     521         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
     522         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
     523         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
     524         CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) )  
     525         CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) ) 
     526         CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) ) 
     527         CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) ) 
     528          
     529         CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) ) 
     530         CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) ) 
     531          
     532         CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     533         CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     534         CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     535         CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) ) 
     536         CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) ) 
     537         CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) ) 
     538          
     539         CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     540         CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     541         CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     542         CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) ) 
     543         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
     544         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
     545          
     546         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
     547         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
     548         CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     549         CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     550         CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     551         CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) ) 
     552         CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) ) 
     553         CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) ) 
     554         CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) ) 
     555         CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) ) 
     556         CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) ) 
     557         CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) ) 
     558         CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) ) 
     559 
     560      CASE( 2 ) 
     561 
     562         CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj ) 
     563         CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj ) 
     564         CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj ) 
     565         CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj ) 
     566         CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj ) 
     567         CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj ) 
     568         DO jk = 1, nlay_s 
     569            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj) 
     570            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj) 
     571         END DO 
     572         DO jk = 1, nlay_i 
     573            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj) 
     574            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj) 
     575            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj) 
     576         END DO 
     577         CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj ) 
     578          
     579         CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj ) 
     580         CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj ) 
     581          
     582         CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     583         CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     584         CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     585         CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     586         CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj ) 
     587         CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj ) 
     588          
     589         CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     590         CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     591         CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     592         CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     593         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
     594         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
     595          
     596         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
     597         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 
     598         CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     599         CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     600         CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     601         CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj ) 
     602         CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj ) 
     603         CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj ) 
     604         CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj ) 
     605         CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj ) 
     606         CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
     607         CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
     608         CALL tab_1d_2d( nbpb, hfx_err_dif   , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj ) 
     609         ! 
     610         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
     611         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
     612         !          
     613      END SELECT 
     614 
     615   END SUBROUTINE lim_thd_1d2d 
     616 
    551617 
    552618   SUBROUTINE lim_thd_init 
     
    563629      !!------------------------------------------------------------------- 
    564630      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    565       NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   & 
    566          &                hiclim, hnzst, parsub, betas,                          &  
    567          &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 
     631      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
     632         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
     633         &                nn_monocat, ln_it_qnsice 
    568634      !!------------------------------------------------------------------- 
    569635      ! 
     
    582648902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 
    583649      IF(lwm) WRITE ( numoni, namicethd ) 
    584  
    585       IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 
     650      ! 
     651      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN  
     652         nn_monocat = 0 
     653         IF(lwp) WRITE(numout, *) '   nn_monocat must be 0 in multi-category case ' 
     654      ENDIF 
     655 
    586656      ! 
    587657      IF(lwp) THEN                          ! control print 
    588658         WRITE(numout,*) 
    589659         WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation ' 
    590          WRITE(numout,*)'      maximum melting at the bottom                           hmelt        = ', hmelt 
    591          WRITE(numout,*)'      ice thick. for lateral accretion                        hiccrit      = ', hiccrit 
    592          WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       fraz_swi     = ', fraz_swi 
    593          WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   maxfrazb     = ', maxfrazb 
    594          WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  vfrazb       = ', vfrazb 
    595          WRITE(numout,*)'      Squeezing coefficient for collection of frazil          Cfrazb       = ', Cfrazb 
    596          WRITE(numout,*)'      minimum ice thickness                                   hiclim       = ', hiclim  
     660         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice 
     661         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil 
     662         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb 
     663         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb    = ', rn_vfrazb 
     664         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb 
     665         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin  
    597666         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice ' 
    598          WRITE(numout,*)'      thickness of the surf. layer in temp. computation       hnzst        = ', hnzst 
    599          WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub   
    600          WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          betas        = ', betas 
    601          WRITE(numout,*)'      extinction radiation parameter in sea ice (1.0)         kappa_i      = ', kappa_i 
    602          WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nconv_i_thd  = ', nconv_i_thd 
    603          WRITE(numout,*)'      maximal err. on T for heat diffusion computation        maxer_i_thd  = ', maxer_i_thd 
    604          WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     thcon_i_swi  = ', thcon_i_swi 
     667         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas 
     668         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
     669         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nn_conv_dif  = ', nn_conv_dif 
     670         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif 
     671         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon 
    605672         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i 
     673         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat 
     674         WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice 
    606675      ENDIF 
    607676      ! 
Note: See TracChangeset for help on using the changeset viewer.