Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5034 r5600  
    1919   !!---------------------------------------------------------------------- 
    2020   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area 
    21    !!   lim_ctl       : alerts in case of ice model crash 
    22    !!   lim_prt_state : ice control print at a given grid point 
    2321   !!---------------------------------------------------------------------- 
    2422   USE oce             ! ocean dynamics and tracers 
    2523   USE dom_oce         ! ocean space and time domain 
    26    USE par_ice         ! sea-ice parameters 
    2724   USE ice             ! LIM-3: ice variables 
    28    USE iceini          ! LIM-3: ice initialisation 
     25   USE thd_ice         ! LIM-3: thermodynamical variables 
    2926   USE dom_ice         ! LIM-3: ice domain 
    3027 
     
    4037   USE limdyn          ! Ice dynamics 
    4138   USE limtrp          ! Ice transport 
     39   USE limhdf          ! Ice horizontal diffusion 
    4240   USE limthd          ! Ice thermodynamics 
    43    USE limitd_th       ! Thermodynamics on ice thickness distribution  
    4441   USE limitd_me       ! Mechanics on ice thickness distribution 
    4542   USE limsbc          ! sea surface boundary condition 
     
    4744   USE limwri          ! Ice outputs 
    4845   USE limrst          ! Ice restarts 
    49    USE limupdate1       ! update of global variables 
    50    USE limupdate2       ! update of global variables 
     46   USE limupdate1      ! update of global variables 
     47   USE limupdate2      ! update of global variables 
    5148   USE limvar          ! Ice variables switch 
     49 
     50   USE limmsh          ! LIM mesh 
     51   USE limistate       ! LIM initial state 
     52   USE limthd_sal      ! LIM ice thermodynamics: salinity 
    5253 
    5354   USE c1d             ! 1D vertical configuration 
     
    6061   USE prtctl          ! Print control 
    6162   USE lib_fortran     !  
     63   USE limctl 
    6264 
    6365#if defined key_bdy  
     
    6971 
    7072   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
    71    PUBLIC lim_prt_state 
     73   PUBLIC sbc_lim_init ! routine called by sbcmod.F90 
    7274    
    7375   !! * Substitutions 
     
    106108      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED) 
    107109      !! 
    108       INTEGER  ::   jl      ! dummy loop index 
    109       REAL(wp) ::   zcoef   ! local scalar 
     110      INTEGER  ::   jl                 ! dummy loop index 
    110111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    111112      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
     113      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice  
    112114      !!---------------------------------------------------------------------- 
    113115 
    114116      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    115117 
    116       IF( kt == nit000 ) THEN 
    117          IF(lwp) WRITE(numout,*) 
    118          IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
    119          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    120          ! 
    121          CALL ice_init 
    122          ! 
    123          IF( ln_nicep ) THEN      ! control print at a given point 
    124             jiindx = 15    ;   jjindx =  44 
    125             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    126          ENDIF 
    127       ENDIF 
    128  
    129       !                                        !----------------------! 
    130       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    131          !                                     !----------------------! 
    132          !                                           !  Bulk Formulae ! 
    133          !                                           !----------------! 
    134          ! 
    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  
     118      !-----------------------! 
     119      ! --- Ice time step --- ! 
     120      !-----------------------! 
     121      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     122 
     123         ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean) 
     124         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 
     125         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 
     126          
     127         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
     128         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 
     129         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
     130           
     131         ! Mask sea ice surface temperature (set to rt0 over land) 
     132         DO jl = 1, jpl 
     133            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
     134         END DO      
     135         ! 
     136         !------------------------------------------------!                                            
     137         ! --- Dynamical coupling with the atmosphere --- !                                            
     138         !------------------------------------------------! 
     139         ! It provides the following fields: 
     140         ! utau_ice, vtau_ice : surface ice stress (U- & V-points)   [N/m2] 
     141         !----------------------------------------------------------------- 
    145142         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              
     143         CASE( jp_clio    )   ;   CALL blk_ice_clio_tau                         ! CLIO bulk formulation             
     144         CASE( jp_core    )   ;   CALL blk_ice_core_tau                         ! CORE bulk formulation 
     145         CASE( jp_purecpl )   ;   CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation 
    153146         END SELECT 
    154147          
    155          !                                           ! Mask sea ice surface temperature 
    156          DO jl = 1, jpl 
    157             t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
    158          END DO 
    159       
    160          ! Bulk formulae  - provides the following fields: 
    161          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
     148         IF( ln_mixcpl) THEN   ! Case of a mixed Bulk/Coupled formulation 
     149            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice) 
     150            CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
     151            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
     152            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
     153            CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice) 
     154         ENDIF 
     155 
     156         !-------------------------------------------------------! 
     157         ! --- ice dynamics and transport (except in 1D case) ---! 
     158         !-------------------------------------------------------! 
     159         numit = numit + nn_fsbc                  ! Ice model time step 
     160         !                                                    
     161         CALL sbc_lim_bef                         ! Store previous ice values 
     162         CALL sbc_lim_diag0                       ! set diag of mass, heat and salt fluxes to 0 
     163         CALL lim_rst_opn( kt )                   ! Open Ice restart file 
     164         ! 
     165         IF( .NOT. lk_c1d ) THEN 
     166            ! 
     167            CALL lim_dyn( kt )                    ! Ice dynamics    ( rheology/dynamics )    
     168            ! 
     169            CALL lim_trp( kt )                    ! Ice transport   ( Advection/diffusion ) 
     170            ! 
     171            IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 
     172            ! 
     173#if defined key_bdy 
     174            CALL bdy_ice_lim( kt )                ! bdy ice thermo  
     175            IF( ln_icectl )       CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
     176#endif 
     177            ! 
     178            CALL lim_update1( kt )                ! Corrections 
     179            ! 
     180         ENDIF 
     181          
     182         ! previous lead fraction and ice volume for flux calculations 
     183         CALL sbc_lim_bef                         
     184         CALL lim_var_glo2eqv                     ! ht_i and ht_s for ice albedo calculation 
     185         CALL lim_var_agg(1)                      ! at_i for coupling (via pfrld)  
     186         pfrld(:,:)   = 1._wp - at_i(:,:) 
     187         phicif(:,:)  = vt_i(:,:) 
     188          
     189         !------------------------------------------------------!                                            
     190         ! --- Thermodynamical coupling with the atmosphere --- !                                            
     191         !------------------------------------------------------! 
     192         ! It provides the following fields: 
    162193         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2] 
    163194         ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2] 
     
    165196         ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s] 
    166197         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%] 
    167          ! 
     198         !---------------------------------------------------------------------------------------- 
     199         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     200         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
     201 
    168202         SELECT CASE( kblk ) 
    169203         CASE( jp_clio )                                       ! CLIO bulk formulation 
    170             CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               & 
    171                &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   & 
    172                &                      qla_ice    , dqns_ice   , dqla_ice  ,               & 
    173                &                      tprecip    , sprecip    ,                           & 
    174                &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  ) 
    175             !          
    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  
     204            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     205            ! (zalb_ice) is computed within the bulk routine 
     206            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice ) 
     207            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
     208            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    179209         CASE( jp_core )                                       ! CORE bulk formulation 
    180             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    181                &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
    182                &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
    183                &                      tprecip   , sprecip   ,                            & 
    184                &                      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 ) 
    188             ! 
    189          CASE ( jp_cpl ) 
    190              
    191             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 
    192  
    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 
    203             ! 
     210            ! albedo depends on cloud fraction because of non-linear spectral effects 
     211            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     212            CALL blk_ice_core_flx( t_su, zalb_ice ) 
     213            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
     214            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     215         CASE ( jp_purecpl ) 
     216            ! albedo depends on cloud fraction because of non-linear spectral effects 
     217            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     218                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
     219            ! clem: evap_ice is forced to 0 in coupled mode for now  
     220            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 
     221            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp 
     222            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    204223         END SELECT 
    205           
    206          !                                           !----------------------! 
    207          !                                           ! LIM-3  time-stepping ! 
    208          !                                           !----------------------! 
    209          !  
    210          numit = numit + nn_fsbc                     ! Ice model time step 
    211          ! 
    212          !                                           ! Store previous ice values 
    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 
    246  
    247                           CALL lim_rst_opn( kt )     ! Open Ice restart file 
    248          ! 
    249          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    250          ! ---------------------------------------------- 
    251          ! ice dynamics and transport (except in 1D case) 
    252          ! ---------------------------------------------- 
    253          IF( .NOT. lk_c1d ) THEN 
    254                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    255                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    256                           CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    257          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    258                           CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    259                           CALL lim_var_agg( 1 )  
    260 #if defined key_bdy 
    261                           ! bdy ice thermo  
    262                           CALL lim_var_glo2eqv            ! equivalent variables 
    263                           CALL bdy_ice_lim( kt ) 
    264                           CALL lim_itd_me_zapsmall 
    265                           CALL lim_var_agg(1) 
    266          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
    267 #endif 
    268                           CALL lim_update1 
    269          ENDIF 
    270 !                         !- Change old values for new values 
    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(:,:,:) 
    280   
    281          ! ---------------------------------------------- 
    282          ! ice thermodynamic 
    283          ! ---------------------------------------------- 
    284                           CALL lim_var_glo2eqv            ! equivalent variables 
    285                           CALL lim_var_agg(1)             ! aggregate ice categories 
    286                           ! previous lead fraction and ice volume for flux calculations 
    287                           pfrld(:,:)   = 1._wp - at_i(:,:) 
    288                           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 
    302                           ! 
    303                           CALL lim_var_bv                 ! bulk brine volume (diag) 
    304                           CALL lim_thd( kt )              ! Ice thermodynamics  
    305                           zcoef = rdt_ice /rday           !  Ice natural aging 
    306                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    307          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    308                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
    309                           CALL lim_var_agg( 1 )           ! requested by limupdate 
    310                           CALL lim_update2                ! Global variables update 
    311  
    312                           CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    313                           CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    314          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    315          ! 
    316                           CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    317          ! 
    318          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    319          ! 
    320          !                                           ! Diagnostics and outputs  
    321          IF (ln_limdiaout) CALL lim_diahsb 
    322  
    323                           CALL lim_wri( 1  )              ! Ice outputs  
    324  
     224         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     225 
     226         !----------------------------! 
     227         ! --- ice thermodynamics --- ! 
     228         !----------------------------! 
     229         CALL lim_thd( kt )                         ! Ice thermodynamics       
     230         ! 
     231         CALL lim_update2( kt )                     ! Corrections 
     232         ! 
     233         CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes 
     234         ! 
     235         IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs  
     236         ! 
     237         CALL lim_wri( 1 )                          ! Ice outputs  
     238         ! 
    325239         IF( kt == nit000 .AND. ln_rstart )   & 
    326             &             CALL iom_close( numrir )        ! clem: close input ice restart file 
    327          ! 
    328          IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file  
    329                           CALL lim_var_glo2eqv            ! ??? 
    330          ! 
    331          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 ) 
    334          ! 
    335       ENDIF                                    ! End sea-ice time step only 
    336  
    337       !                                        !--------------------------! 
    338       !                                        !  at all ocean time step  ! 
    339       !                                        !--------------------------! 
    340       !                                                
    341       !                                              ! Update surface ocean stresses (only in ice-dynamic case) 
    342       !                                                   ! otherwise the atm.-ocean stresses are used everywhere 
     240            &             CALL iom_close( numrir )  ! close input ice restart file 
     241         ! 
     242         IF( lrst_ice )   CALL lim_rst_write( kt )  ! Ice restart file  
     243         ! 
     244         IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash 
     245         ! 
     246      ENDIF   ! End sea-ice time step only 
     247 
     248      !-------------------------! 
     249      ! --- Ocean time step --- ! 
     250      !-------------------------! 
     251      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 
    343252      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
    344253!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    345  
    346       ! 
    347       IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
     254      ! 
     255      IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim') 
    348256      ! 
    349257   END SUBROUTINE sbc_ice_lim 
    350258    
     259 
     260   SUBROUTINE sbc_lim_init 
     261      !!---------------------------------------------------------------------- 
     262      !!                  ***  ROUTINE sbc_lim_init  *** 
     263      !! 
     264      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules 
     265      !!---------------------------------------------------------------------- 
     266      INTEGER :: ierr 
     267      !!---------------------------------------------------------------------- 
     268      IF(lwp) WRITE(numout,*) 
     269      IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
     270      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
     271      ! 
     272                                       ! Open the reference and configuration namelist files and namelist output file  
     273      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )  
     274      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     275      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) 
     276 
     277      CALL ice_run                     ! set some ice run parameters 
     278      ! 
     279      !                                ! Allocate the ice arrays 
     280      ierr =        ice_alloc        ()      ! ice variables 
     281      ierr = ierr + dom_ice_alloc    ()      ! domain 
     282      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
     283      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
     284      ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
     285      ! 
     286      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     287      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays') 
     288      ! 
     289      !                                ! adequation jpk versus ice/snow layers/categories 
     290      IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   & 
     291         &      CALL ctl_stop( 'STOP',                          & 
     292         &     'sbc_lim_init: the 3rd dimension of workspace arrays is too small.',   & 
     293         &     'use more ocean levels or less ice/snow layers/categories.' ) 
     294      ! 
     295      CALL lim_itd_init                ! ice thickness distribution initialization 
     296      ! 
     297      CALL lim_hdf_init                ! set ice horizontal diffusion computation parameters 
     298      ! 
     299      CALL lim_thd_init                ! set ice thermodynics parameters 
     300      ! 
     301      CALL lim_thd_sal_init            ! set ice salinity parameters 
     302      ! 
     303      CALL lim_msh                     ! ice mesh initialization 
     304      ! 
     305      CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
     306      !                                ! Initial sea-ice state 
     307      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     308         numit = 0 
     309         numit = nit000 - 1 
     310         CALL lim_istate 
     311      ELSE                                    ! start from a restart file 
     312         CALL lim_rst_read 
     313         numit = nit000 - 1 
     314      ENDIF 
     315      CALL lim_var_agg(1) 
     316      CALL lim_var_glo2eqv 
     317      ! 
     318      CALL lim_sbc_init                 ! ice surface boundary condition    
     319      ! 
     320      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     321      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
     322      ! 
     323      nstart = numit  + nn_fsbc       
     324      nitrun = nitend - nit000 + 1  
     325      nlast  = numit  + nitrun  
     326      ! 
     327      IF( nstock == 0 )   nstock = nlast + 1 
     328      ! 
     329   END SUBROUTINE sbc_lim_init 
     330 
     331 
     332   SUBROUTINE ice_run 
     333      !!------------------------------------------------------------------- 
     334      !!                  ***  ROUTINE ice_run *** 
     335      !!                  
     336      !! ** Purpose :   Definition some run parameter for ice model 
     337      !! 
     338      !! ** Method  :   Read the namicerun namelist and check the parameter  
     339      !!              values called at the first timestep (nit000) 
     340      !! 
     341      !! ** input   :   Namelist namicerun 
     342      !!------------------------------------------------------------------- 
     343      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     344      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  & 
     345         &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
     346      !!------------------------------------------------------------------- 
     347      !                     
     348      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice 
     349      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 
     350901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
     351 
     352      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice 
     353      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 
     354902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
     355      IF(lwm) WRITE ( numoni, namicerun ) 
     356      ! 
     357      ! 
     358      IF(lwp) THEN                        ! control print 
     359         WRITE(numout,*) 
     360         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     361         WRITE(numout,*) ' ~~~~~~' 
     362         WRITE(numout,*) '   number of ice  categories                               = ', jpl 
     363         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i 
     364         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
     365         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
     366         WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax  
     367         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
     368         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
     369         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_icectl 
     370         WRITE(numout,*) '   i-index for control prints (ln_icectl=true)             = ', iiceprt 
     371         WRITE(numout,*) '   j-index for control prints (ln_icectl=true)             = ', jiceprt 
     372      ENDIF 
     373      ! 
     374      ! sea-ice timestep and inverse 
     375      rdt_ice   = nn_fsbc * rdttra(1)   
     376      r1_rdtice = 1._wp / rdt_ice  
     377 
     378      ! inverse of nlay_i and nlay_s 
     379      r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 
     380      r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 
     381      ! 
     382#if defined key_bdy 
     383      IF( lwp .AND. ln_limdiahsb )  CALL ctl_warn('online conservation check activated but it does not work with BDY') 
     384#endif 
     385      ! 
     386   END SUBROUTINE ice_run 
     387 
     388 
     389   SUBROUTINE lim_itd_init 
     390      !!------------------------------------------------------------------ 
     391      !!                ***  ROUTINE lim_itd_init *** 
     392      !! 
     393      !! ** Purpose :   Initializes the ice thickness distribution 
     394      !! ** Method  :   ... 
     395      !! ** input   :   Namelist namiceitd 
     396      !!------------------------------------------------------------------- 
     397      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     398      NAMELIST/namiceitd/ nn_catbnd, rn_himean 
     399      ! 
     400      INTEGER  ::   jl                   ! dummy loop index 
     401      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
     402      REAL(wp) ::   zhmax, znum, zden, zalpha ! 
     403      !!------------------------------------------------------------------ 
     404      ! 
     405      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice 
     406      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 903) 
     407903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
     408 
     409      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice 
     410      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 904 ) 
     411904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
     412      IF(lwm) WRITE ( numoni, namiceitd ) 
     413      ! 
     414      ! 
     415      IF(lwp) THEN                        ! control print 
     416         WRITE(numout,*) 
     417         WRITE(numout,*) 'ice_itd : ice cat distribution' 
     418         WRITE(numout,*) ' ~~~~~~' 
     419         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd 
     420         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean 
     421      ENDIF 
     422 
     423      !---------------------------------- 
     424      !- Thickness categories boundaries  
     425      !---------------------------------- 
     426      IF(lwp) WRITE(numout,*) 
     427      IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 
     428      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     429 
     430      hi_max(:) = 0._wp 
     431 
     432      SELECT CASE ( nn_catbnd  )        
     433                                   !---------------------- 
     434         CASE (1)                  ! tanh function (CICE) 
     435                                   !---------------------- 
     436         zc1 =  3._wp / REAL( jpl, wp ) 
     437         zc2 = 10._wp * zc1 
     438         zc3 =  3._wp 
     439 
     440         DO jl = 1, jpl 
     441            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 
     442            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 
     443         END DO 
     444 
     445                                   !---------------------- 
     446         CASE (2)                  ! h^(-alpha) function 
     447                                   !---------------------- 
     448         zalpha = 0.05             ! exponent of the transform function 
     449 
     450         zhmax  = 3.*rn_himean 
     451 
     452         DO jl = 1, jpl  
     453            znum = jpl * ( zhmax+1 )**zalpha 
     454            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl 
     455            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     456         END DO 
     457 
     458      END SELECT 
     459 
     460      DO jl = 1, jpl 
     461         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     462      END DO 
     463 
     464      ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 
     465      hi_max(jpl) = 99._wp 
     466 
     467      IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 
     468      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 
     469      ! 
     470   END SUBROUTINE lim_itd_init 
     471 
    351472    
    352       SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
    353          &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
     473   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 
    354474      !!--------------------------------------------------------------------- 
    355       !!                  ***  ROUTINE sbc_ice_lim  *** 
     475      !!                  ***  ROUTINE ice_lim_flx  *** 
    356476      !!                    
    357477      !! ** Purpose :   update the ice surface boundary condition by averaging and / or 
     
    369489      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux 
    370490      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 
     491      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation 
     492      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity 
    373493      ! 
    374494      INTEGER  ::   jl      ! dummy loop index 
     
    379499      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories 
    380500      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 
     501      REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m  ! Mean sublimation over all categories 
    382502      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 
     503      REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 
    384504      !!---------------------------------------------------------------------- 
    385505 
     
    389509      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==! 
    390510      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 (:,:,:) ) 
     511         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 
     512         ! 
     513         z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 
     514         z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 
     515         z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 
     516         z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) ) 
     517         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) ) 
    398518         DO jl = 1, jpl 
    399             pdqn_ice(:,:,jl) = z_dqn_m(:,:) 
    400             pdql_ice(:,:,jl) = z_dql_m(:,:) 
     519            pdqn_ice  (:,:,jl) = z_dqn_m(:,:) 
     520            pdevap_ice(:,:,jl) = z_devap_m(:,:) 
    401521         END DO 
    402522         ! 
    403523         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(:,:) 
     524            pqns_ice (:,:,jl) = z_qns_m(:,:) 
     525            pqsr_ice (:,:,jl) = z_qsr_m(:,:) 
     526            pevap_ice(:,:,jl) = z_evap_m(:,:) 
    407527         END DO 
    408528         ! 
    409          CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     529         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 
    410530      END SELECT 
    411531 
     
    417537         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )  
    418538         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(:,:) )  
     539            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 
     540            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 
     541            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )  
    422542         END DO 
    423543         ! 
     
    428548      ! 
    429549   END SUBROUTINE ice_lim_flx 
    430     
    431     
    432    SUBROUTINE lim_ctl( kt ) 
    433       !!----------------------------------------------------------------------- 
    434       !!                   ***  ROUTINE lim_ctl ***  
    435       !!                  
    436       !! ** Purpose :   Alerts in case of model crash 
    437       !!------------------------------------------------------------------- 
    438       INTEGER, INTENT(in) ::   kt      ! ocean time step 
    439       INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    440       INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    441       INTEGER  ::   ialert_id         ! number of the current alert 
    442       REAL(wp) ::   ztmelts           ! ice layer melting point 
    443       CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert 
    444       INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive 
    445       !!------------------------------------------------------------------- 
    446  
    447       inb_altests = 10 
    448       inb_alp(:)  =  0 
    449  
    450       ! Alert if incompatible volume and concentration 
    451       ialert_id = 2 ! reference number of this alert 
    452       cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    453  
    454       DO jl = 1, jpl 
    455          DO jj = 1, jpj 
    456             DO ji = 1, jpi 
    457                IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    458                   !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    459                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    460                   !WRITE(numout,*) ' Point - category', 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) 
    463                   !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    464                   !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
    465                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    466                ENDIF 
    467             END DO 
    468          END DO 
    469       END DO 
    470  
    471       ! Alerte if very thick ice 
    472       ialert_id = 3 ! reference number of this alert 
    473       cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert 
    474       jl = jpl  
    475       DO jj = 1, jpj 
    476          DO ji = 1, jpi 
    477             IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    478                !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    479                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    480             ENDIF 
    481          END DO 
    482       END DO 
    483  
    484       ! Alert if very fast ice 
    485       ialert_id = 4 ! reference number of this alert 
    486       cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert 
    487       DO jj = 1, jpj 
    488          DO ji = 1, jpi 
    489             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
    490                &  at_i(ji,jj) > 0._wp   ) THEN 
    491                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    492                !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    493                !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    494                !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    495                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    496                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    497                !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    498                !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    499                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    500                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    501                !WRITE(numout,*)  
    502                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    503             ENDIF 
    504          END DO 
    505       END DO 
    506  
    507       ! Alert if there is ice on continents 
    508       ialert_id = 6 ! reference number of this alert 
    509       cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert 
    510       DO jj = 1, jpj 
    511          DO ji = 1, jpi 
    512             IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    513                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    514                !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    515                !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    516                !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    517                !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    518                !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    519                !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    520                !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    521                !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    522                ! 
    523                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    524             ENDIF 
    525          END DO 
    526       END DO 
    527  
    528 ! 
    529 !     ! Alert if very fresh ice 
    530       ialert_id = 7 ! reference number of this alert 
    531       cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert 
    532       DO jl = 1, jpl 
    533          DO jj = 1, jpj 
    534             DO ji = 1, jpi 
    535                IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    536 !                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    537 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    538 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    539 !                 WRITE(numout,*)  
    540                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    541                ENDIF 
    542             END DO 
    543          END DO 
    544       END DO 
    545 ! 
    546  
    547 !     ! Alert if too old ice 
    548       ialert_id = 9 ! reference number of this alert 
    549       cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert 
    550       DO jl = 1, jpl 
    551          DO jj = 1, jpj 
    552             DO ji = 1, jpi 
    553                IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. & 
    554                       ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    555                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    556                   !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    557                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    558                ENDIF 
    559             END DO 
    560          END DO 
    561       END DO 
    562   
    563       ! Alert on salt flux 
    564       ialert_id = 5 ! reference number of this alert 
    565       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    566       DO jj = 1, jpj 
    567          DO ji = 1, jpi 
    568             IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    569                !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    570                !DO jl = 1, jpl 
    571                   !WRITE(numout,*) ' Category no: ', jl 
    572                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    573                   !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    574                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    575                   !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    576                   !WRITE(numout,*) ' ' 
    577                !END DO 
    578                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    579             ENDIF 
    580          END DO 
    581       END DO 
    582  
    583       ! Alert if qns very big 
    584       ialert_id = 8 ! reference number of this alert 
    585       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    586       DO jj = 1, jpj 
    587          DO ji = 1, jpi 
    588             IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    589                ! 
    590                !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    591                !WRITE(numout,*) ' ji, jj    : ', ji, jj 
    592                !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    593                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    594                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    595                ! 
    596                !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
    597                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    598                ! 
    599             ENDIF 
    600          END DO 
    601       END DO 
    602       !+++++ 
    603   
    604       ! Alert if very warm ice 
    605       ialert_id = 10 ! reference number of this alert 
    606       cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert 
    607       inb_alp(ialert_id) = 0 
    608       DO jl = 1, jpl 
    609          DO jk = 1, nlay_i 
    610             DO jj = 1, jpj 
    611                DO ji = 1, jpi 
    612                   ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt 
    613                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    614                      &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    615                      !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    616                      !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    617                      !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    618                      !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    619                      !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
    620                      !WRITE(numout,*) ' ztmelts : ', ztmelts 
    621                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    622                   ENDIF 
    623                END DO 
    624             END DO 
    625          END DO 
    626       END DO 
    627  
    628       ! sum of the alerts on all processors 
    629       IF( lk_mpp ) THEN 
    630          DO ialert_id = 1, inb_altests 
    631             CALL mpp_sum(inb_alp(ialert_id)) 
    632          END DO 
    633       ENDIF 
    634  
    635       ! print alerts 
    636       IF( lwp ) THEN 
    637          ialert_id = 1                                 ! reference number of this alert 
    638          cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    639          WRITE(numout,*) ' time step ',kt 
    640          WRITE(numout,*) ' All alerts at the end of ice model ' 
    641          DO ialert_id = 1, inb_altests 
    642             WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
    643          END DO 
    644       ENDIF 
    645      ! 
    646    END SUBROUTINE lim_ctl 
    647   
    648     
    649    SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 
    650       !!----------------------------------------------------------------------- 
    651       !!                   ***  ROUTINE lim_prt_state ***  
    652       !!                  
    653       !! ** Purpose :   Writes global ice state on the (i,j) point  
    654       !!                in ocean.ouput  
    655       !!                3 possibilities exist  
    656       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
    657       !!                n = 2    -> exhaustive state 
    658       !!                n = 3    -> ice/ocean salt fluxes 
    659       !! 
    660       !! ** input   :   point coordinates (i,j)  
    661       !!                n : number of the option 
    662       !!------------------------------------------------------------------- 
    663       INTEGER         , INTENT(in) ::   kt            ! ocean time step 
    664       INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    665       CHARACTER(len=*), INTENT(in) ::   cd1           ! 
    666       !! 
    667       INTEGER :: jl, ji, jj 
    668       !!------------------------------------------------------------------- 
    669  
    670       DO ji = mi0(ki), mi1(ki) 
    671          DO jj = mj0(kj), mj1(kj) 
    672  
    673             WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title 
    674  
    675             !---------------- 
    676             !  Simple state 
    677             !---------------- 
    678              
    679             IF ( kn == 1 .OR. kn == -1 ) THEN 
    680                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    681                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    682                WRITE(numout,*) ' Simple state ' 
    683                WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
    684                WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj) 
    685                WRITE(numout,*) ' Time step     : ', numit 
    686                WRITE(numout,*) ' - Ice drift   ' 
    687                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    688                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    689                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    690                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    691                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    692                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    693                WRITE(numout,*) 
    694                WRITE(numout,*) ' - Cell values ' 
    695                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    696                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    697                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    698                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    699                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    700                DO jl = 1, jpl 
    701                   WRITE(numout,*) ' - Category (', jl,')' 
    702                   WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
    703                   WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl) 
    704                   WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl) 
    705                   WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
    706                   WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
    707                   WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
    708                   WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
    709                   WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
    710                   WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
    711                   WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl) 
    712                   WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl) 
    713                   WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl) 
    714                   WRITE(numout,*) 
    715                END DO 
    716             ENDIF 
    717             IF( kn == -1 ) THEN 
    718                WRITE(numout,*) ' Mechanical Check ************** ' 
    719                WRITE(numout,*) ' Check what means ice divergence ' 
    720                WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
    721                WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
    722                WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
    723                WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
    724             ENDIF 
    725              
    726  
    727             !-------------------- 
    728             !  Exhaustive state 
    729             !-------------------- 
    730              
    731             IF ( kn .EQ. 2 ) THEN 
    732                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    733                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    734                WRITE(numout,*) ' Exhaustive state ' 
    735                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    736                WRITE(numout,*) ' Time step ', numit 
    737                WRITE(numout,*)  
    738                WRITE(numout,*) ' - Cell values ' 
    739                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    740                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    741                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    742                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    743                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    744                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    745                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    746                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    747                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    748                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    749                WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
    750                WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    751                WRITE(numout,*) 
    752                 
    753                DO jl = 1, jpl 
    754                   WRITE(numout,*) ' - Category (',jl,')' 
    755                   WRITE(numout,*) '   ~~~~~~~~         '  
    756                   WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl) 
    757                   WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl) 
    758                   WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
    759                   WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
    760                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)    
    761                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    762                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    763                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    764                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    765                   WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
    766                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    767                   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 
    768                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
    769                   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 
    770                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    771                   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) 
    772                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
    773                   WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
    774                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    775                   WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
    776                END DO !jl 
    777                 
    778                WRITE(numout,*) 
    779                WRITE(numout,*) ' - Heat / FW fluxes ' 
    780                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    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,*) 
    785                WRITE(numout,*)  
    786                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
    787                WRITE(numout,*) ' sss        : ', sss_m(ji,jj)   
    788                WRITE(numout,*)  
    789                WRITE(numout,*) ' - Stresses ' 
    790                WRITE(numout,*) '   ~~~~~~~~ ' 
    791                WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj)  
    792                WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj) 
    793                WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
    794                WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
    795                WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
    796                WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
    797             ENDIF 
    798              
    799             !--------------------- 
    800             ! Salt / heat fluxes 
    801             !--------------------- 
    802              
    803             IF ( kn .EQ. 3 ) THEN 
    804                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    805                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    806                WRITE(numout,*) ' - Salt / Heat Fluxes ' 
    807                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    808                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    809                WRITE(numout,*) ' Time step ', numit 
    810                WRITE(numout,*) 
    811                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    812                WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    813                WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    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 
    825                WRITE(numout,*) 
    826                WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    827                WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    828                WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    829                WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    830                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    831                WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    832                WRITE(numout,*) 
    833                WRITE(numout,*) ' - Momentum fluxes ' 
    834                WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    835                WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    836             ENDIF  
    837             WRITE(numout,*) ' ' 
    838             ! 
    839          END DO 
    840       END DO 
    841       ! 
    842    END SUBROUTINE lim_prt_state 
    843     
     550 
     551   SUBROUTINE sbc_lim_bef 
     552      !!---------------------------------------------------------------------- 
     553      !!                  ***  ROUTINE sbc_lim_bef  *** 
     554      !! 
     555      !! ** purpose :  store ice variables at "before" time step  
     556      !!---------------------------------------------------------------------- 
     557      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area 
     558      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
     559      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
     560      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume  
     561      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
     562      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content 
     563      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content 
     564      u_ice_b(:,:)     = u_ice(:,:) 
     565      v_ice_b(:,:)     = v_ice(:,:) 
     566       
     567   END SUBROUTINE sbc_lim_bef 
     568 
     569   SUBROUTINE sbc_lim_diag0 
     570      !!---------------------------------------------------------------------- 
     571      !!                  ***  ROUTINE sbc_lim_diag0  *** 
     572      !! 
     573      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining 
     574      !!               of the time step 
     575      !!---------------------------------------------------------------------- 
     576      sfx    (:,:) = 0._wp   ; 
     577      sfx_bri(:,:) = 0._wp   ;  
     578      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     579      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     580      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     581      sfx_res(:,:) = 0._wp 
     582       
     583      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     584      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     585      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     586      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     587      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     588      wfx_spr(:,:) = 0._wp   ;    
     589       
     590      hfx_thd(:,:) = 0._wp   ;    
     591      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     592      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     593      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     594      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     595      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     596      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     597      hfx_err_dif(:,:) = 0._wp   ; 
     598 
     599      afx_tot(:,:) = 0._wp   ; 
     600      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp 
     601 
     602      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ; 
     603      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ; 
     604       
     605   END SUBROUTINE sbc_lim_diag0 
     606 
    844607      
    845608   FUNCTION fice_cell_ave ( ptab ) 
     
    852615       
    853616      fice_cell_ave (:,:) = 0.0_wp 
    854        
    855617      DO jl = 1, jpl 
    856          fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
    857             &                  + a_i (:,:,jl) * ptab (:,:,jl) 
     618         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl) 
    858619      END DO 
    859620       
     
    869630 
    870631      fice_ice_ave (:,:) = 0.0_wp 
    871       WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
     632      WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
    872633 
    873634   END FUNCTION fice_ice_ave 
     
    882643      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk 
    883644   END SUBROUTINE sbc_ice_lim 
     645   SUBROUTINE sbc_lim_init                 ! Dummy routine 
     646   END SUBROUTINE sbc_lim_init 
    884647#endif 
    885648 
Note: See TracChangeset for help on using the changeset viewer.