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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5075 for branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 – NEMO

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4333 r5075  
    1212   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation 
    1313   !!             -   ! 2012-10  (C. Rousset)  add lim_diahsb 
     14   !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface 
    1415   !!---------------------------------------------------------------------- 
    1516#if defined key_lim3 
     
    6869 
    6970   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
     71   PUBLIC lim_prt_state 
    7072    
    7173   !! * Substitutions 
     
    7880   !!---------------------------------------------------------------------- 
    7981CONTAINS 
    80  
    81    FUNCTION fice_cell_ave ( ptab) 
    82       !!-------------------------------------------------------------------------- 
    83       !! * Compute average over categories, for grid cell (ice covered and free ocean) 
    84       !!-------------------------------------------------------------------------- 
    85       REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 
    86       REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 
    87       INTEGER :: jl ! Dummy loop index 
    88        
    89       fice_cell_ave (:,:) = 0.0_wp 
    90        
    91       DO jl = 1, jpl 
    92          fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
    93             &                  + a_i (:,:,jl) * ptab (:,:,jl) 
    94       END DO 
    95        
    96    END FUNCTION fice_cell_ave 
    97     
    98    FUNCTION fice_ice_ave ( ptab) 
    99       !!-------------------------------------------------------------------------- 
    100       !! * Compute average over categories, for ice covered part of grid cell 
    101       !!-------------------------------------------------------------------------- 
    102       REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 
    103       REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 
    104  
    105       fice_ice_ave (:,:) = 0.0_wp 
    106       WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
    107  
    108    END FUNCTION fice_ice_ave 
    10982 
    11083   !!====================================================================== 
     
    131104      !!--------------------------------------------------------------------- 
    132105      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    133       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE) 
     106      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED) 
    134107      !! 
    135108      INTEGER  ::   jl      ! dummy loop index 
    136109      REAL(wp) ::   zcoef   ! local scalar 
    137       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
    138       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled) 
    139  
    140       REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories 
    141       REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories 
    142        
    143       REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories 
    144       REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories 
    145       REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories 
    146       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories 
    147       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories 
     110      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
     111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
    148112      !!---------------------------------------------------------------------- 
    149113 
    150       !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ????? 
    151  
    152114      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    153  
    154       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    155  
    156 #if defined key_coupled 
    157       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
    158       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    159          &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    160 #endif 
    161115 
    162116      IF( kt == nit000 ) THEN 
     
    168122         ! 
    169123         IF( ln_nicep ) THEN      ! control print at a given point 
    170             jiindx = 177   ;   jjindx = 112 
     124            jiindx = 15    ;   jjindx =  44 
    171125            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    172126         ENDIF 
     
    176130      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    177131         !                                     !----------------------! 
    178          !                                           !  Bulk Formulea ! 
     132         !                                           !  Bulk Formulae ! 
    179133         !                                           !----------------! 
    180134         ! 
    181          u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point 
    182          v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    183          ! 
    184          t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin] 
    185          !                                           ! (set to rt0 over land) 
    186          CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
    187  
     135         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point 
     136         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean) 
     137         ! 
     138         t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )  ! masked sea surface freezing temperature [Kelvin] 
     139         !                                                                                  ! (set to rt0 over land) 
     140         !                                           ! Ice albedo 
     141         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )       
     142 
     143         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     144 
     145         SELECT CASE( kblk ) 
     146         CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations 
     147 
     148            ! albedo depends on cloud fraction because of non-linear spectral effects 
     149            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     150            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     151            ! (zalb_ice) is computed within the bulk routine 
     152             
     153         END SELECT 
     154          
     155         !                                           ! Mask sea ice surface temperature 
    188156         DO jl = 1, jpl 
    189157            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
    190158         END DO 
    191  
    192          IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
    193           
    194 #if defined key_coupled 
    195          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    196             ! 
    197             ! Compute mean albedo and temperature 
    198             zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
    199             ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
    200             ! 
    201          ENDIF 
    202 #endif 
    203                                                ! Bulk formulea - provides the following fields: 
     159      
     160         ! Bulk formulae  - provides the following fields: 
    204161         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
    205162         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2] 
     
    210167         ! 
    211168         SELECT CASE( kblk ) 
    212          CASE( 3 )                                       ! CLIO bulk formulation 
    213             CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           & 
     169         CASE( jp_clio )                                       ! CLIO bulk formulation 
     170            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               & 
    214171               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   & 
    215172               &                      qla_ice    , dqns_ice   , dqla_ice  ,               & 
     
    217174               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  ) 
    218175            !          
    219          CASE( 4 )                                       ! CORE bulk formulation 
    220             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               & 
     176            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     177               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     178 
     179         CASE( jp_core )                                       ! CORE bulk formulation 
     180            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    221181               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
    222182               &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
    223183               &                      tprecip   , sprecip   ,                            & 
    224184               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  ) 
     185               ! 
     186            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     187               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
    225188            ! 
    226          CASE ( 5 ) 
    227             zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
     189         CASE ( jp_cpl ) 
    228190             
    229191            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 
    230192 
    231             CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    ) 
    232  
    233             ! Latent heat flux is forced to 0 in coupled : 
    234             !  it is included in qns (non-solar heat flux) 
    235             qla_ice  (:,:,:) = 0.0e0_wp 
    236             dqla_ice (:,:,:) = 0.0e0_wp 
     193            ! MV -> seb  
     194!           CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
     195 
     196!           IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     197!              &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     198!           ! Latent heat flux is forced to 0 in coupled : 
     199!           !  it is included in qns (non-solar heat flux) 
     200!           qla_ice  (:,:,:) = 0._wp 
     201!           dqla_ice (:,:,:) = 0._wp 
     202            ! END MV -> seb 
    237203            ! 
    238204         END SELECT 
    239  
    240          ! Average over all categories 
    241 #if defined key_coupled 
    242          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    243  
    244             z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) ) 
    245             z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) ) 
    246             z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) ) 
    247             z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) ) 
    248             z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) ) 
    249  
    250             DO jl = 1, jpl 
    251                dqns_ice (:,:,jl) = z_dqns_ice_all (:,:) 
    252                dqla_ice (:,:,jl) = z_dqla_ice_all (:,:) 
    253             END DO 
    254             ! 
    255             IF ( ln_iceflx_ave ) THEN 
    256                DO jl = 1, jpl 
    257                   qns_ice  (:,:,jl) = z_qns_ice_all  (:,:) 
    258                   qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:) 
    259                   qla_ice  (:,:,jl) = z_qla_ice_all  (:,:) 
    260                END DO 
    261             END IF 
    262             ! 
    263             IF ( ln_iceflx_linear ) THEN 
    264                DO jl = 1, jpl 
    265                   qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
    266                   qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
    267                   qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:) 
    268                END DO 
    269             END IF 
    270          END IF 
    271 #endif 
     205          
    272206         !                                           !----------------------! 
    273207         !                                           ! LIM-3  time-stepping ! 
     
    277211         ! 
    278212         !                                           ! Store previous ice values 
    279 !!gm : remark   old_...   should becomes ...b  as tn versus tb   
    280          old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area 
    281          old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
    282          old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
    283          old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume  
    284          old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
    285          old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content 
    286          old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content 
    287          ! 
    288          old_u_ice(:,:) = u_ice(:,:) 
    289          old_v_ice(:,:) = v_ice(:,:) 
    290          !                                           ! intialisation to zero    !!gm is it truly necessary ??? 
    291          d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
    292          d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp 
    293          d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp 
    294          d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp 
    295          d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp 
    296          d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
    297          d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
    298          ! 
    299          d_u_ice_dyn(:,:) = 0._wp 
    300          d_v_ice_dyn(:,:) = 0._wp 
    301          ! 
    302          sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp 
    303          sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp 
    304          fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
    305          fhmec  (:,:) = 0._wp   ;    
    306          fmmec  (:,:) = 0._wp 
    307          fmmflx (:,:) = 0._wp      
    308          focea2D(:,:) = 0._wp 
    309          fsup2D (:,:) = 0._wp 
    310  
    311          ! used in limthd.F90 
    312          rdvosif(:,:) = 0._wp   ! variation of ice volume at surface 
    313          rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom 
    314          fdvolif(:,:) = 0._wp   ! total variation of ice volume 
    315          rdvonif(:,:) = 0._wp   ! lateral variation of ice volume 
    316          fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice 
    317          ffltbif(:,:) = 0._wp   ! linked with fstric 
    318          qfvbq  (:,:) = 0._wp   ! linked with fstric 
    319          rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area 
    320          rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area 
    321          hicifp (:,:) = 0._wp   ! daily thermodynamic ice production.  
    322          ! 
    323          diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp 
    324          diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp 
    325          diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp 
    326          diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp 
    327          ! dynamical invariants 
    328          delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp 
     213         a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area 
     214         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
     215         v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
     216         v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume  
     217         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
     218         smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content 
     219         oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content 
     220         u_ice_b(:,:)     = u_ice(:,:) 
     221         v_ice_b(:,:)     = v_ice(:,:) 
     222 
     223         ! salt, heat and mass fluxes 
     224         sfx    (:,:) = 0._wp   ; 
     225         sfx_bri(:,:) = 0._wp   ;  
     226         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     227         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     228         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     229         sfx_res(:,:) = 0._wp 
     230 
     231         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     232         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     233         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     234         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     235         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     236         wfx_spr(:,:) = 0._wp   ;    
     237 
     238         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     239         hfx_thd(:,:) = 0._wp   ;    
     240         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     241         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     242         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     243         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     244         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     245         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
    329246 
    330247                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
     
    352269         ENDIF 
    353270!                         !- Change old values for new values 
    354                           old_u_ice(:,:)   = u_ice (:,:) 
    355                           old_v_ice(:,:)   = v_ice (:,:) 
    356                           old_a_i(:,:,:)   = a_i (:,:,:) 
    357                           old_v_s(:,:,:)   = v_s (:,:,:) 
    358                           old_v_i(:,:,:)   = v_i (:,:,:) 
    359                           old_e_s(:,:,:,:) = e_s (:,:,:,:) 
    360                           old_e_i(:,:,:,:) = e_i (:,:,:,:) 
    361                           old_oa_i(:,:,:)  = oa_i(:,:,:) 
    362                           old_smv_i(:,:,:) = smv_i (:,:,:) 
     271                          u_ice_b(:,:)     = u_ice(:,:) 
     272                          v_ice_b(:,:)     = v_ice(:,:) 
     273                          a_i_b  (:,:,:)   = a_i (:,:,:) 
     274                          v_s_b  (:,:,:)   = v_s (:,:,:) 
     275                          v_i_b  (:,:,:)   = v_i (:,:,:) 
     276                          e_s_b  (:,:,:,:) = e_s (:,:,:,:) 
     277                          e_i_b  (:,:,:,:) = e_i (:,:,:,:) 
     278                          oa_i_b (:,:,:)   = oa_i (:,:,:) 
     279                          smv_i_b(:,:,:)   = smv_i(:,:,:) 
    363280  
    364281         ! ---------------------------------------------- 
     
    370287                          pfrld(:,:)   = 1._wp - at_i(:,:) 
    371288                          phicif(:,:)  = vt_i(:,:) 
     289 
     290                          ! MV -> seb 
     291                          SELECT CASE( kblk ) 
     292                             CASE ( jp_cpl ) 
     293                             CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
     294                             IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     295                          &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     296                           ! Latent heat flux is forced to 0 in coupled : 
     297                           !  it is included in qns (non-solar heat flux) 
     298                             qla_ice  (:,:,:) = 0._wp 
     299                             dqla_ice (:,:,:) = 0._wp 
     300                          END SELECT 
     301                          ! END MV -> seb 
    372302                          ! 
    373303                          CALL lim_var_bv                 ! bulk brine volume (diag) 
     
    375305                          zcoef = rdt_ice /rday           !  Ice natural aging 
    376306                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    377                           CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    378307         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    379308                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
     
    391320         !                                           ! Diagnostics and outputs  
    392321         IF (ln_limdiaout) CALL lim_diahsb 
    393 !clem # if ! defined key_iomput 
     322 
    394323                          CALL lim_wri( 1  )              ! Ice outputs  
    395 !clem # endif 
     324 
    396325         IF( kt == nit000 .AND. ln_rstart )   & 
    397326            &             CALL iom_close( numrir )        ! clem: close input ice restart file 
     
    401330         ! 
    402331         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash 
     332         ! 
     333         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
    403334         ! 
    404335      ENDIF                                    ! End sea-ice time step only 
     
    411342      !                                                   ! otherwise the atm.-ocean stresses are used everywhere 
    412343      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
    413        
    414344!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    415       ! 
    416       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    417  
    418 #if defined key_coupled 
    419       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
    420       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    421          &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    422 #endif 
     345 
    423346      ! 
    424347      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
    425348      ! 
    426349   END SUBROUTINE sbc_ice_lim 
    427  
    428  
     350    
     351    
     352      SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
     353         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
     354      !!--------------------------------------------------------------------- 
     355      !!                  ***  ROUTINE sbc_ice_lim  *** 
     356      !!                    
     357      !! ** Purpose :   update the ice surface boundary condition by averaging and / or 
     358      !!                redistributing fluxes on ice categories                    
     359      !! 
     360      !! ** Method  :   average then redistribute  
     361      !! 
     362      !! ** Action  :    
     363      !!--------------------------------------------------------------------- 
     364      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;  
     365                                                                ! =1 average and redistribute ; =2 redistribute 
     366      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature  
     367      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo 
     368      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux 
     369      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux 
     370      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity 
     371      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux 
     372      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity 
     373      ! 
     374      INTEGER  ::   jl      ! dummy loop index 
     375      ! 
     376      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m    ! Mean albedo over all categories 
     377      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m    ! Mean temperature over all categories 
     378      ! 
     379      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories 
     380      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories 
     381      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories 
     382      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories 
     383      REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories 
     384      !!---------------------------------------------------------------------- 
     385 
     386      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx') 
     387      ! 
     388      ! 
     389      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==! 
     390      CASE( 0 , 1 ) 
     391         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     392         ! 
     393         z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 
     394         z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 
     395         z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 
     396         z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) ) 
     397         z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) ) 
     398         DO jl = 1, jpl 
     399            pdqn_ice(:,:,jl) = z_dqn_m(:,:) 
     400            pdql_ice(:,:,jl) = z_dql_m(:,:) 
     401         END DO 
     402         ! 
     403         DO jl = 1, jpl 
     404            pqns_ice(:,:,jl) = z_qns_m(:,:) 
     405            pqsr_ice(:,:,jl) = z_qsr_m(:,:) 
     406            pqla_ice(:,:,jl) = z_qla_m(:,:) 
     407         END DO 
     408         ! 
     409         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     410      END SELECT 
     411 
     412      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==! 
     413      CASE( 1 , 2 ) 
     414         CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m ) 
     415         ! 
     416         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )  
     417         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )  
     418         DO jl = 1, jpl 
     419            pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
     420            pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
     421            pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )  
     422         END DO 
     423         ! 
     424         CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m ) 
     425      END SELECT 
     426      ! 
     427      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx') 
     428      ! 
     429   END SUBROUTINE ice_lim_flx 
     430    
     431    
    429432   SUBROUTINE lim_ctl( kt ) 
    430433      !!----------------------------------------------------------------------- 
     
    456459                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    457460                  !WRITE(numout,*) ' Point - category', ji, jj, jl 
    458                   !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
    459                   !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
     461                  !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
     462                  !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
    460463                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    461464                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
     
    534537!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    535538!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    536 !                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl) 
    537539!                 WRITE(numout,*)  
    538540                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    568570               !DO jl = 1, jpl 
    569571                  !WRITE(numout,*) ' Category no: ', jl 
    570                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
     572                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    571573                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    572                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
     574                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    573575                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    574576                  !WRITE(numout,*) ' ' 
     
    591593               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    592594               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    593                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
    594                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
    595                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
    596                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
    597                !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
    598                !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
    599                !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
    600                !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
    601                !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
    602                !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
    603                !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
    604                !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
    605                !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
    606595               ! 
    607596               !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
     
    672661      !!                n : number of the option 
    673662      !!------------------------------------------------------------------- 
    674       INTEGER         , INTENT(in) ::   kt      ! ocean time step 
     663      INTEGER         , INTENT(in) ::   kt            ! ocean time step 
    675664      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    676665      CHARACTER(len=*), INTENT(in) ::   cd1           ! 
     
    759748               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    760749               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
    761                WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)   
     750               WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    762751               WRITE(numout,*) 
    763752                
     
    769758                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
    770759                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
    771                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)    
     760                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)    
    772761                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    773                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)    
     762                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    774763                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    775                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)   
     764                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    776765                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
    777                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9  
     766                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    778767                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9 
    779                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9   
     768                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
    780769                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9 
    781                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl)  
     770                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    782771                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl) 
    783                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)    
     772                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
    784773                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
    785                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl) 
     774                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    786775                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
    787776               END DO !jl 
     
    790779               WRITE(numout,*) ' - Heat / FW fluxes ' 
    791780               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    792                WRITE(numout,*) ' emp        : ', emp      (ji,jj) 
    793                WRITE(numout,*) ' sfx        : ', sfx      (ji,jj) 
    794                WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj) 
    795                WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj) 
    796                WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj) 
    797                WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj) 
    798                WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj) 
    799                WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
    800                WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
    801                WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
     781               WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
     782               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) ) 
     783               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 
     784               WRITE(numout,*) 
    802785               WRITE(numout,*)  
    803786               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
     
    829812               WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    830813               WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    831                WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
    832                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice 
    833                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice 
     814               WRITE(numout,*) 
     815               WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj) 
     816               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
     817               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
     818               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
     819               WRITE(numout,*) 
     820               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
     821               WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj) 
     822               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
     823               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
     824               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
    834825               WRITE(numout,*) 
    835826               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    836827               WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    837                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    838828               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    839829               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    840                WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj) 
    841                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    842                WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
     830               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
     831               WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    843832               WRITE(numout,*) 
    844833               WRITE(numout,*) ' - Momentum fluxes ' 
    845834               WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    846835               WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    847             ENDIF 
     836            ENDIF  
    848837            WRITE(numout,*) ' ' 
    849838            ! 
    850839         END DO 
    851840      END DO 
    852  
     841      ! 
    853842   END SUBROUTINE lim_prt_state 
     843    
     844      
     845   FUNCTION fice_cell_ave ( ptab ) 
     846      !!-------------------------------------------------------------------------- 
     847      !! * Compute average over categories, for grid cell (ice covered and free ocean) 
     848      !!-------------------------------------------------------------------------- 
     849      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 
     850      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 
     851      INTEGER :: jl ! Dummy loop index 
     852       
     853      fice_cell_ave (:,:) = 0.0_wp 
     854       
     855      DO jl = 1, jpl 
     856         fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
     857            &                  + a_i (:,:,jl) * ptab (:,:,jl) 
     858      END DO 
     859       
     860   END FUNCTION fice_cell_ave 
     861    
     862    
     863   FUNCTION fice_ice_ave ( ptab ) 
     864      !!-------------------------------------------------------------------------- 
     865      !! * Compute average over categories, for ice covered part of grid cell 
     866      !!-------------------------------------------------------------------------- 
     867      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 
     868      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 
     869 
     870      fice_ice_ave (:,:) = 0.0_wp 
     871      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
     872 
     873   END FUNCTION fice_ice_ave 
     874 
    854875 
    855876#else 
Note: See TracChangeset for help on using the changeset viewer.