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 4634 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 – NEMO

Ignore:
Timestamp:
2014-05-12T22:46:18+02:00 (10 years ago)
Author:
clem
Message:

major changes in heat budget

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4332 r4634  
    3232   USE par_ice 
    3333   USE limitd_th 
     34   USE limitd_me 
    3435   USE limvar 
    3536   USE prtctl           ! Print control 
     
    4950 
    5051      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    51       REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    52       REAL(wp)  ::   rone   = 1._wp       !    -       - 
    5352          
    5453   !! * Substitutions 
     
    8079      INTEGER ::   jbnd1, jbnd2 
    8180      INTEGER ::   i_ice_switch 
    82       INTEGER ::   ind_im, layer      ! indices for internal melt 
    83       REAL(wp) ::   zweight, zesum, z_da_i, zhimax 
    8481      REAL(wp) ::   zinda, zindb, zindsn, zindic 
    85       REAL(wp) ::   zindg, zh, zdvres, zviold2 
    86       REAL(wp) ::   zbigvalue, zvsold2, z_da_ex 
    87       REAL(wp) ::   z_prescr_hi, zat_i_old, ztmelts, ze_s 
    88  
    89       REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for 
    90       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     82      REAL(wp) ::   zh, zdvres, zsal 
     83      REAL(wp) ::   zat_i_old, ztmelts 
     84 
     85      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    9186      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    92       ! mass and salt flux (clem) 
    93       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
    9487      !!------------------------------------------------------------------- 
    9588      IF( nn_timing == 1 )  CALL timing_start('limupdate1') 
    96  
    97       CALL wrk_alloc( jkmax, zthick0, zqm0 ) 
    98  
    99       CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
    10089 
    10190      !------------------------------------------------------------------------------ 
     
    10695      !  Trend terms 
    10796      !----------------- 
    108       d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
    109       d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
    110       d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
    111       d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
    112       d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
    113       d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
    114       d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
    115       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    116       d_smv_i_trp(:,:,:)   = 0._wp 
    117       IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    118  
    119       ! mass and salt flux init (clem) 
    120       zviold(:,:,:) = v_i(:,:,:) 
    121       zvsold(:,:,:) = v_s(:,:,:) 
    122       zsmvold(:,:,:) = smv_i(:,:,:) 
    12397 
    12498      ! ------------------------------- 
    12599      !- check conservation (C Rousset) 
    126100      IF (ln_limdiahsb) THEN 
    127          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     101         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    128102         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    129          zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
    130          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     103         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
     104         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
     105         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
     106         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    131107      ENDIF 
    132108      !- check conservation (C Rousset) 
    133109      ! ------------------------------- 
    134110 
     111      ! zap small values 
     112      !----------------- 
     113      CALL lim_itd_me_zapsmall 
     114 
    135115      CALL lim_var_glo2eqv 
    136  
    137       !-------------------------------------- 
    138       ! 2. Review of all pathological cases 
    139       !-------------------------------------- 
    140  
    141 ! clem: useless now 
    142       !------------------------------------------- 
    143       ! 2.1) Advection of ice in an ice-free cell 
    144       !------------------------------------------- 
    145       ! should be removed since it is treated after dynamics now 
    146 !      zhimax = 5._wp 
    147 !      ! first category 
    148 !      DO jj = 1, jpj 
    149 !         DO ji = 1, jpi 
    150 !            !--- the thickness of such an ice is often out of bounds 
    151 !            !--- thus we recompute a new area while conserving ice volume 
    152 !            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
    153 !            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) )  
    154 !            IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 
    155 !              &   .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
    156 !              &   .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
    157 !               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
    158 !               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
    159 !            ENDIF 
    160 !         END DO 
    161 !      END DO 
    162 ! 
    163 !      zhimax = 20._wp 
    164 !      ! other categories 
    165 !      DO jl = 2, jpl 
    166 !         jm = ice_types(jl) 
    167 !         DO jj = 1, jpj 
    168 !            DO ji = 1, jpi 
    169 !               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) )  
    170 !               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
    171 !               ! it makes problems when the advected volume and concentration do not seem to be  
    172 !               ! related with each other 
    173 !               ! the new thickness is sometimes very big! 
    174 !               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
    175 !               ! which of course is plausible 
    176 !               ! but fuck! it fucks everything up :) 
    177 !               IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
    178 !                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
    179 !                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
    180 !                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
    181 !               ENDIF 
    182 !            END DO ! ji 
    183 !         END DO !jj 
    184 !      END DO !jl 
    185116      
    186117      at_i(:,:) = 0._wp 
     
    203134      END DO 
    204135 
    205       zbigvalue      = 1.0e+20 
    206  
    207       DO jl = 1, jpl 
    208          DO jj = 1, jpj  
     136      !--- 2.13 ice concentration should not exceed amax  
     137      !----------------------------------------------------- 
     138      DO jl  = 1, jpl 
     139         DO jj = 1, jpj 
    209140            DO ji = 1, jpi 
    210  
    211                !switches 
    212                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    213                !switch = 1 if a_i > 1e-06 and 0 if not 
    214                zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 
    215                zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 
    216                ! bug fix 25 avril 2007 
    217                zindb         = zindb*zindic 
    218  
    219                !--- 2.3 Correction to ice age  
    220                !------------------------------ 
    221                !                IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN 
    222                !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday 
    223                !                ENDIF 
    224                IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 
    225                   oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl) 
     141               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     142                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 
     143                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    226144               ENDIF 
    227                oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 
    228  
    229                !--- 2.4 Correction to snow thickness 
    230                !------------------------------------- 
    231                !          ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0 
    232                !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)  
    233                ! snow thickness cannot be smaller than 1e-6 
    234                zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl) 
    235                v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 
    236  
    237                !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn 
    238   
    239                !--- 2.5 Correction to ice thickness 
    240                !------------------------------------- 
    241                zdvres = (zindb - 1._wp) * v_i(ji,jj,jl) 
    242                v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 
    243  
    244                !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic 
    245                !sfx_res(ji,jj)  = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice ) 
    246  
    247                !--- 2.6 Snow is transformed into ice if the original ice cover disappears 
    248                !---------------------------------------------------------------------------- 
    249                zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
    250                zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0 
    251                v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres 
    252  
    253                zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 
    254                v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres 
    255  
    256                !--- 2.7 Correction to ice concentrations  
    257                !-------------------------------------------- 
    258                ! if greater than 0, ice concentration cannot be smaller than 1e-10 
    259                a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 
    260  
    261                !------------------------- 
    262                ! 2.8) Snow heat content 
    263                !------------------------- 
    264                e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0._wp, e_s(ji,jj,1,jl) ), zbigvalue ) ) 
    265  
    266             END DO ! ji 
    267          END DO ! jj 
    268       END DO ! jl 
    269  
    270       !------------------------ 
    271       ! 2.9) Ice heat content  
    272       !------------------------ 
    273  
    274       DO jl = 1, jpl 
    275          DO jk = 1, nlay_i 
    276             DO jj = 1, jpj  
    277                DO ji = 1, jpi 
    278                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    279                   e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 
    280                END DO ! ji 
    281             END DO ! jj 
    282          END DO !jk 
    283       END DO !jl 
    284   
    285       at_i(:,:) = 0._wp 
     145            END DO 
     146         END DO 
     147      END DO 
     148 
     149      at_i(:,:) = 0.0 
    286150      DO jl = 1, jpl 
    287151         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    288152      END DO 
    289  
    290       !--- 2.13 ice concentration should not exceed amax  
    291       !         (it should not be the case)  
    292       !----------------------------------------------------- 
    293       DO jj = 1, jpj 
    294          DO ji = 1, jpi 
    295             z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
    296             zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) )  
    297             DO jl  = 1, jpl 
    298                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 
    299                a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
    300                ! 
    301                zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    302                ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 
    303                !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    304             END DO 
    305          END DO 
    306       END DO 
    307       at_i(:,:) = a_i(:,:,1) 
    308       DO jl = 2, jpl 
    309          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    310       END DO 
    311  
    312  
     153   
     154   
    313155      ! Final thickness distribution rebinning 
    314156      ! -------------------------------------- 
     
    322164 
    323165 
     166      ! zap small values 
     167      !----------------- 
     168      CALL lim_itd_me_zapsmall 
     169 
    324170      !--------------------- 
    325171      ! 2.11) Ice salinity 
    326172      !--------------------- 
    327       ! clem correct bug on smv_i 
    328       smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    329  
    330       IF (  num_sal == 2  ) THEN ! general case 
     173      IF (  num_sal == 2  ) THEN  
    331174         DO jl = 1, jpl 
    332             !DO jk = 1, nlay_i 
    333                DO jj = 1, jpj  
    334                   DO ji = 1, jpi 
    335                      ! salinity stays in bounds 
    336                      !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 
    337                      smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
    338                      i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
    339                      smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
    340                   END DO ! ji 
    341                END DO ! jj 
    342             !END DO !jk 
     175            DO jj = 1, jpj  
     176               DO ji = 1, jpi 
     177                  zsal            = smv_i(ji,jj,jl) 
     178                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 
     179                  ! salinity stays in bounds 
     180                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
     181                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
     182                  ! associated salt flux 
     183                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     184               END DO ! ji 
     185            END DO ! jj 
    343186         END DO !jl 
    344187      ENDIF 
    345188 
     189      ! ------------------- 
    346190      at_i(:,:) = a_i(:,:,1) 
    347191      DO jl = 2, jpl 
    348192         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    349193      END DO 
    350  
    351  
    352       !-------------------------------- 
    353       ! Update mass/salt fluxes (clem) 
    354       !-------------------------------- 
    355       DO jl = 1, jpl 
    356          DO jj = 1, jpj  
    357             DO ji = 1, jpi 
    358                diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice  
    359                rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic  
    360                rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn  
    361                sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice  
    362             END DO 
    363          END DO 
    364       END DO 
    365194   
     195 
     196      ! ------------------------------------------------- 
     197      ! Diagnostics 
     198      ! ------------------------------------------------- 
     199      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
     200      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
     201      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
     202      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
     203      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
     204      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
     205      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:) 
     206      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
     207      d_smv_i_trp(:,:,:)   = 0._wp 
     208      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
     209 
    366210      ! ------------------------------- 
    367211      !- check conservation (C Rousset) 
    368       IF (ln_limdiahsb) THEN 
    369  
    370          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    371          zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     212      IF( ln_limdiahsb ) THEN 
     213         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     214         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     215         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    372216  
    373          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     217         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    374218         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     219         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    375220 
    376221         zchk_vmin = glob_min(v_i) 
    377222         zchk_amax = glob_max(SUM(a_i,dim=3)) 
    378223         zchk_amin = glob_min(a_i) 
    379         
     224 
    380225         IF(lwp) THEN 
    381             IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * rday) 
     226            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate1) = ',(zchk_v_i * rday) 
    382227            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday) 
     228            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate1) = ',(zchk_e_i) 
    383229            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3) 
    384230            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax 
    385231            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin 
    386232         ENDIF 
    387       ENDIF 
     233       ENDIF 
    388234      !- check conservation (C Rousset) 
    389235      ! ------------------------------- 
     
    446292         CALL prt_ctl_info(' - Heat / FW fluxes : ') 
    447293         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ') 
    448          CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ') 
    449294         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ') 
    450          CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ') 
    451295 
    452296         CALL prt_ctl_info(' ') 
     
    458302      ENDIF 
    459303    
    460  
    461       CALL wrk_dealloc( jkmax, zthick0, zqm0 ) 
    462  
    463       CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
    464  
    465304      IF( nn_timing == 1 )  CALL timing_stop('limupdate1') 
    466305   END SUBROUTINE lim_update1 
Note: See TracChangeset for help on using the changeset viewer.