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 5086 for branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 – NEMO

Ignore:
Timestamp:
2015-02-17T10:06:39+01:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4333 r5086  
    55   !!====================================================================== 
    66   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code 
     7   !!            3.6  !  2014-06  (C. Rousset)       Complete rewriting/cleaning 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim3 
     
    3233   USE par_ice 
    3334   USE limitd_th 
     35   USE limitd_me 
    3436   USE limvar 
    3537   USE prtctl           ! Print control 
     
    3739   USE wrk_nemo         ! work arrays 
    3840   USE lib_fortran     ! glob_sum 
    39    ! Check budget (Rousset) 
    4041   USE in_out_manager   ! I/O manager 
    4142   USE iom              ! I/O manager 
    4243   USE lib_mpp          ! MPP library 
    4344   USE timing          ! Timing 
     45   USE limcons        ! conservation tests 
    4446 
    4547   IMPLICIT NONE 
     
    4850   PUBLIC   lim_update1   ! routine called by ice_step 
    4951 
    50       REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    51       REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    52       REAL(wp)  ::   rone   = 1._wp       !    -       - 
    53           
    5452   !! * Substitutions 
    5553#  include "vectopt_loop_substitute.h90" 
     
    6664      !!                
    6765      !! ** Purpose :  Computes update of sea-ice global variables at  
    68       !!               the end of the time step. 
    69       !!               Address pathological cases 
    70       !!               This place is very important 
     66      !!               the end of the dynamics. 
    7167      !!                 
    72       !! ** Method  :   
    73       !!    Ice speed from ice dynamics 
    74       !!    Ice thickness, Snow thickness, Temperatures, Lead fraction 
    75       !!      from advection and ice thermodynamics  
    76       !! 
    77       !! ** Action  : -  
    7868      !!--------------------------------------------------------------------- 
    79       INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices 
    80       INTEGER ::   jbnd1, jbnd2 
    81       INTEGER ::   i_ice_switch 
    82       INTEGER ::   ind_im, layer      ! indices for internal melt 
    83       REAL(wp) ::   zweight, zesum, z_da_i, zhimax 
    84       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) 
    91       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... 
     69      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     70      INTEGER  ::   i_ice_switch 
     71      REAL(wp) ::   zsal 
     72      ! 
     73      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    9474      !!------------------------------------------------------------------- 
    9575      IF( nn_timing == 1 )  CALL timing_start('limupdate1') 
    9676 
    97       CALL wrk_alloc( jkmax, zthick0, zqm0 ) 
    98  
    99       CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
    100  
    101       !------------------------------------------------------------------------------ 
    102       ! 1. Update of Global variables                                               | 
    103       !------------------------------------------------------------------------------ 
    104  
    105       !----------------- 
    106       !  Trend terms 
    107       !----------------- 
    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(:,:,:) 
    123  
    124       ! ------------------------------- 
    125       !- check conservation (C Rousset) 
    126       IF (ln_limdiahsb) THEN 
    127          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    128          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(:,:) ) 
    131       ENDIF 
    132       !- check conservation (C Rousset) 
    133       ! ------------------------------- 
     77      IF( ln_limdyn ) THEN  
     78 
     79      ! conservation test 
     80      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     81 
     82      !----------------- 
     83      ! zap small values 
     84      !----------------- 
     85      CALL lim_itd_me_zapsmall 
    13486 
    13587      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 
    18588      
     89      !---------------------------------------------------- 
     90      ! Rebin categories with thickness out of bounds 
     91      !---------------------------------------------------- 
     92      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl) 
     93 
    18694      at_i(:,:) = 0._wp 
    18795      DO jl = 1, jpl 
     
    19098 
    19199      !---------------------------------------------------- 
    192       ! 2.2) Rebin categories with thickness out of bounds 
    193       !---------------------------------------------------- 
    194       DO jm = 1, jpm 
    195          jbnd1 = ice_cat_bounds(jm,1) 
    196          jbnd2 = ice_cat_bounds(jm,2) 
    197          IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 
     100      ! ice concentration should not exceed amax  
     101      !----------------------------------------------------- 
     102      DO jl  = 1, jpl 
     103         DO jj = 1, jpj 
     104            DO ji = 1, jpi 
     105               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     106                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 
     107                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     108               ENDIF 
     109            END DO 
     110         END DO 
    198111      END DO 
    199112 
     
    202115         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    203116      END DO 
    204  
    205       zbigvalue      = 1.0e+20 
    206  
    207       DO jl = 1, jpl 
    208          DO jj = 1, jpj  
    209             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) 
    226                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 
     117     
     118      ! -------------------------------------- 
     119      ! Final thickness distribution rebinning 
     120      ! -------------------------------------- 
     121      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
     122 
     123      !----------------- 
     124      ! zap small values 
     125      !----------------- 
     126      CALL lim_itd_me_zapsmall 
     127 
     128      !--------------------- 
     129      ! Ice salinity bounds 
     130      !--------------------- 
     131      IF (  num_sal == 2  ) THEN  
     132         DO jl = 1, jpl 
    276133            DO jj = 1, jpj  
    277134               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 
    286       DO jl = 1, jpl 
    287          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    288       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 
     135                  zsal            = smv_i(ji,jj,jl) 
     136                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 
     137                  ! salinity stays in bounds 
     138                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
     139                  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) ) 
     140                  ! associated salt flux 
     141                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     142               END DO 
    304143            END DO 
    305144         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  
    313       ! Final thickness distribution rebinning 
    314       ! -------------------------------------- 
    315       DO jm = 1, jpm 
    316          jbnd1 = ice_cat_bounds(jm,1) 
    317          jbnd2 = ice_cat_bounds(jm,2) 
    318          IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 
    319          IF (ice_ncat_types(jm) .EQ. 1 ) THEN 
    320          ENDIF 
    321       END DO 
    322  
    323  
    324       !--------------------- 
    325       ! 2.11) Ice salinity 
    326       !--------------------- 
    327       ! clem correct bug on smv_i 
    328       smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    329  
    330       IF (  num_sal == 2  ) THEN ! general case 
    331          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 
    343          END DO !jl 
    344145      ENDIF 
    345146 
    346       at_i(:,:) = a_i(:,:,1) 
    347       DO jl = 2, jpl 
    348          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    349       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 
    365    
    366       ! ------------------------------- 
    367       !- 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 
    372   
    373          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
    374          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    375  
    376          zchk_vmin = glob_min(v_i) 
    377          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    378          zchk_amin = glob_min(a_i) 
    379         
    380          IF(lwp) THEN 
    381             IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * rday) 
    382             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday) 
    383             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3) 
    384             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax 
    385             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin 
    386          ENDIF 
    387       ENDIF 
    388       !- check conservation (C Rousset) 
    389       ! ------------------------------- 
     147      ! ------------------------------------------------- 
     148      ! Diagnostics 
     149      ! ------------------------------------------------- 
     150      d_u_ice_dyn(:,:)     = u_ice(:,:)     - u_ice_b(:,:) 
     151      d_v_ice_dyn(:,:)     = v_ice(:,:)     - v_ice_b(:,:) 
     152      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - a_i_b  (:,:,:) 
     153      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - v_s_b  (:,:,:)   
     154      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - v_i_b  (:,:,:)    
     155      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - e_s_b  (:,:,:,:)   
     156      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:) 
     157      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - oa_i_b (:,:,:) 
     158      d_smv_i_trp(:,:,:)   = 0._wp 
     159      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:) 
     160 
     161      ! conservation test 
     162      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    390163 
    391164      IF(ln_ctl) THEN   ! Control print 
     
    400173         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :') 
    401174         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :') 
    402          CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :') 
     175         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :') 
    403176 
    404177         DO jl = 1, jpl 
     
    413186            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ') 
    414187            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ') 
    415             CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ') 
     188            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ') 
    416189            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ') 
    417190            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ') 
    418             CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ') 
     191            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ') 
    419192            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ') 
    420193            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ') 
    421             CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ') 
     194            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ') 
    422195            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ') 
    423196            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ') 
    424             CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ') 
     197            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ') 
    425198            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ') 
    426199            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ') 
    427             CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ') 
     200            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ') 
    428201            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ') 
    429202            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ') 
    430             CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ') 
     203            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ') 
    431204            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ') 
    432205            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ') 
    433             CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ') 
     206            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ') 
    434207            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ') 
    435208            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ') 
    436             CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ') 
     209            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ') 
    437210            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ') 
    438211 
     
    446219         CALL prt_ctl_info(' - Heat / FW fluxes : ') 
    447220         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ') 
    448          CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ') 
    449221         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 : ') 
    451222 
    452223         CALL prt_ctl_info(' ') 
     
    458229      ENDIF 
    459230    
    460  
    461       CALL wrk_dealloc( jkmax, zthick0, zqm0 ) 
    462  
    463       CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
     231      ENDIF ! ln_limdyn 
    464232 
    465233      IF( nn_timing == 1 )  CALL timing_stop('limupdate1') 
Note: See TracChangeset for help on using the changeset viewer.