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.
limupdate.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

  • Property svn:keywords set to Id
File size: 47.8 KB
Line 
1MODULE limupdate
2   !!======================================================================
3   !!                     ***  MODULE  limupdate  ***
4   !!    Update of sea-ice global variables
5   !!    at the end of the time step
6   !!   
7   !!    Ice speed from ice dynamics
8   !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
9   !!      from advection and ice thermodynamics
10   !!======================================================================
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!    lim_update   : computes update of sea-ice global variables from trend terms
16   !!----------------------------------------------------------------------
17   USE limrhg          ! ice rheology
18   USE lbclnk
19
20   USE dom_oce
21   USE oce             ! dynamics and tracers variables
22   USE in_out_manager
23   USE sbc_oce         ! Surface boundary condition: ocean fields
24   USE sbc_ice         ! Surface boundary condition: ice fields
25   USE dom_ice
26   USE phycst          ! Define parameters for the routines
27   USE ice
28   USE lbclnk
29   USE limdyn
30   USE limtrp
31   USE limthd
32   USE limsbc
33   USE limdia
34   USE limwri
35   USE limrst
36   USE thd_ice         ! LIM thermodynamic sea-ice variables
37   USE par_ice
38   USE limitd_th
39   USE limvar
40   USE prtctl          ! Print control
41
42
43   IMPLICIT NONE
44   PRIVATE
45
46   !! * Accessibility
47   PUBLIC lim_update ! routine called by ice_step
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51
52   !!----------------------------------------------------------------------
53   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE lim_update
61      !!-------------------------------------------------------------------
62      !!               ***  ROUTINE lim_update  ***
63      !!               
64      !! ** Purpose :  Computes update of sea-ice global variables at
65      !!               the end of the time step.
66      !!               Address pathological cases
67      !!               This place is very important
68      !!               
69      !! ** Method  :  Mathematical
70      !!
71      !! ** Action  : -
72      !!
73      !! History : This routine was new for LIM 3.0
74      !!   3.0  !  04-06  (M. Vancoppenolle) Tendencies
75      !!---------------------------------------------------------------------
76      !! * Local variables
77      INTEGER ::      &
78         ji, jj,     & ! geographical indices
79         jk, jl, jm    ! layer, category and type indices
80      INTEGER ::      &
81         jbnd1, jbnd2
82      INTEGER ::      &
83         i_ice_switch
84
85      REAL(wp)  ::           &  ! constant values
86         epsi06 = 1.e-06  ,  &
87         epsi03 = 1.e-03  ,  &
88         epsi16 = 1.e-16  ,  &
89         epsi20 = 1.e-20  ,  &
90         epsi04 = 1.e-04  ,  &
91         epsi10 = 1.e-10  ,  &
92         rzero  = 0.e0    ,  &
93         rone   = 1.e0    ,  &
94         zhimax                   ! maximum thickness tolerated for advection of
95      ! in an ice-free cell
96      REAL(wp) ::            &  ! dummy switches and arguments
97         zindb, zindsn, zindic, zacrith,  &
98         zrtt, zindg, zh, zdvres, zviold,                       &
99         zbigvalue, zvsold, z_da_ex, zamax,                     &
100         z_prescr_hi, zat_i_old,                                &
101         ztmelts, ze_s
102
103      REAL(wp), DIMENSION(jpl) :: z_da_i, z_dv_i
104
105      LOGICAL, DIMENSION(jpi,jpj,jpl) ::  &
106         internal_melt
107
108      INTEGER ::      &
109         ind_im, layer      ! indices for internal melt
110      REAL(wp), DIMENSION(jkmax) :: &
111         zthick0, zqm0      ! thickness of the layers and heat contents for
112      ! internal melt
113      REAL(wp) ::                   &
114         zweight, zesum
115
116
117      !!-------------------------------------------------------------------
118
119      IF( ln_nicep ) THEN 
120         WRITE(numout,*) ' lim_update '
121         WRITE(numout,*) ' ~~~~~~~~~~ '
122
123         WRITE(numout,*) ' O) Initial values '
124         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
125         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
126         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
127         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
128         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
129         DO jk = 1, nlay_i
130            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
131         END DO
132      ENDIF
133
134      !------------------------------------------------------------------------------
135      ! 1. Update of Global variables                                               |
136      !------------------------------------------------------------------------------
137
138      !---------------------
139      ! Ice dynamics 
140      !---------------------
141
142      u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
143      v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
144
145      !-----------------------------
146      ! Update ice and snow volumes 
147      !-----------------------------
148
149      DO jl = 1, jpl
150         DO jj = 1, jpj
151            DO ji = 1, jpi
152
153               v_i(ji,jj,jl)  = v_i(ji,jj,jl) + d_v_i_trp(ji,jj,jl)  &
154                  + d_v_i_thd(ji,jj,jl) 
155               v_s(ji,jj,jl)  = v_s(ji,jj,jl) + d_v_s_trp(ji,jj,jl)  &
156                  + d_v_s_thd(ji,jj,jl)
157            END DO
158         END DO
159      END DO
160
161      !---------------------------------
162      ! Classify the pathological cases
163      !---------------------------------
164      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
165      ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
166      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
167      ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation
168      ! with negative advection, very pathological )
169      ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
170
171      DO jl = 1, jpl
172         DO jj = 1, jpj
173            DO ji = 1, jpi
174               patho_case(ji,jj,jl) = 1
175               IF ( v_i(ji,jj,jl) .GE. 0.0 ) THEN
176                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
177                     patho_case(ji,jj,jl) = 2
178                  ENDIF
179               ELSE
180                  patho_case(ji,jj,jl) = 3
181                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
182                     patho_case(ji,jj,jl) = 4
183                  ENDIF
184               ENDIF
185               IF ( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. &
186                  ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
187                  patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free
188                  ! cell
189                  IF( ln_nicep ) THEN 
190                     WRITE(numout,*) ' ALERTE patho_case still equal to 5 '
191                     WRITE(numout,*) ' ji , jj   : ', ji, jj
192                     WRITE(numout,*) ' old_v_i   : ', old_v_i(ji,jj,jl)
193                     WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl)
194                  ENDIF
195
196               ENDIF
197            END DO
198         END DO
199      END DO
200
201      !--------------------
202      ! Excessive ablation
203      !--------------------
204
205      DO jl = 1, jpl
206         DO jj = 1, jpj
207            DO ji = 1, jpi
208               IF (      ( patho_case(ji,jj,jl) .EQ. 3 ) &
209                  .OR. ( patho_case(ji,jj,jl) .EQ. 4 ) ) THEN
210                  zviold         = old_v_i(ji,jj,jl)
211                  zvsold         = old_v_s(ji,jj,jl)
212                  ! in cases 3 ( combined total ablation )
213                  !      and 4 ( total ablation with negative advection )
214                  ! there is excessive total ablation
215                  ! advection is chosen to be prioritary in order to conserve mass.
216                  ! dv_i_thd is computed as a residual
217                  ! negative energy has to be kept in memory and to be given to the ocean
218                  ! equivalent salt flux is given to the ocean
219                  !
220                  ! This was the best solution found. Otherwise, mass conservation in advection
221                  ! scheme should have been revised, which could have been a big problem
222                  ! Martin Vancoppenolle (2006, updated 2007)
223
224                  ! is there any ice left ?
225                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
226                  !=1 if hi > 1e-3 and 0 if not
227                  zdvres        = MAX(0.0,-v_i(ji,jj,jl)) !residual volume if too much ice was molten
228                  !this quantity is positive
229                  v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)    !ice volume cannot be negative
230                  !correct thermodynamic ablation
231                  d_v_i_thd(ji,jj,jl)  = zindic  *  d_v_i_thd(ji,jj,jl) + & 
232                     (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl)) 
233                  ! THIS IS NEW
234                  d_a_i_thd(ji,jj,jl)  = zindic  *  d_a_i_thd(ji,jj,jl) + & 
235                     (1.0-zindic) * (-old_a_i(ji,jj,jl)) 
236
237                  !residual salt flux if ice is over-molten
238                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * & 
239                     ( rhoic * zdvres / rdt_ice )
240                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres / rdt_ice
241
242                  ! is there any snow left ?
243                  zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 
244                  zvsold        = v_s(ji,jj,jl)
245                  zdvres        = MAX(0.0,-v_s(ji,jj,jl)) !residual volume if too much ice was molten
246                  !this quantity is positive
247                  v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl)    !snow volume cannot be negative
248                  !correct thermodynamic ablation
249                  d_v_s_thd(ji,jj,jl)  = zindsn  *  d_v_s_thd(ji,jj,jl) + & 
250                     (1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl)) 
251                  !unsure correction on salt flux.... maybe future will tell it was not that right
252
253                  !residual salt flux if snow is over-molten
254                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sss_m(ji,jj) * & 
255                     ( rhosn * zdvres / rdt_ice )
256                  !this flux will be positive if snow was over-molten
257                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres / rdt_ice
258               ENDIF
259            END DO !ji
260         END DO !jj
261      END DO !jl
262
263      IF( ln_nicep ) THEN 
264         DO jj = 1, jpj
265            DO ji = 1, jpi
266               IF ( ABS(fsalt_res(ji,jj)) .GT. 1.0 ) THEN
267                  WRITE(numout,*) ' ALERTE 1000 : residual salt flux of -> ', &
268                     fsalt_res(ji,jj)
269                  WRITE(numout,*) ' ji, jj : ', ji, jj, ' gphit, glamt : ', &
270                     gphit(ji,jj), glamt(ji,jj)
271               ENDIF
272            END DO
273         END DO
274
275         WRITE(numout,*) ' 1. Before update of Global variables '
276         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
277         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
278         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
279         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
280         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
281         DO jk = 1, nlay_i
282            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
283         END DO
284      ENDIF
285
286      !---------------------------------------------
287      ! Ice concentration and ice heat content
288      !---------------------------------------------
289
290      a_i (:,:,:) = a_i (:,:,:)   + d_a_i_trp(:,:,:)     &
291         + d_a_i_thd(:,:,:)
292      CALL lim_var_glo2eqv ! useless, just for debug
293      IF( ln_nicep ) THEN
294         DO jk = 1, nlay_i
295            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
296         END DO
297      ENDIF
298      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:) 
299      CALL lim_var_glo2eqv ! useless, just for debug
300      IF( ln_nicep) THEN
301      WRITE(numout,*) ' After transport update '
302         DO jk = 1, nlay_i
303            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
304         END DO
305      ENDIF
306      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:) 
307      CALL lim_var_glo2eqv ! useless, just for debug
308      IF( ln_nicep ) THEN
309         WRITE(numout,*) ' After thermodyn update '
310         DO jk = 1, nlay_i
311            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
312         END DO
313      ENDIF
314
315      at_i(:,:) = 0.0
316      DO jl = 1, jpl
317         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
318      END DO
319
320      IF( ln_nicep ) THEN 
321         WRITE(numout,*) ' 1. After update of Global variables (2) '
322         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
323         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
324         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
325         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
326         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
327         WRITE(numout,*) ' oa_i : ', oa_i(jiindx, jjindx, 1:jpl)
328         WRITE(numout,*) ' e_s : ', e_s(jiindx, jjindx, 1, 1:jpl)
329         DO jk = 1, nlay_i
330            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
331         END DO
332      ENDIF
333
334      !------------------------------
335      ! Snow temperature and ice age
336      !------------------------------
337
338      e_s(:,:,:,:) = e_s(:,:,:,:)        + &
339         d_e_s_trp(:,:,:,:)  + &
340         d_e_s_thd(:,:,:,:)
341
342      oa_i(:,:,:)  = oa_i(:,:,:)         + &
343         d_oa_i_trp(:,:,:)   + &
344         d_oa_i_thd(:,:,:)
345
346      !--------------
347      ! Ice salinity   
348      !--------------
349
350      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
351
352         IF( ln_nicep ) THEN 
353            WRITE(numout,*) ' Before everything '
354            WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
355            WRITE(numout,*) ' oa_i:  ', oa_i(jiindx, jjindx, 1:jpl)
356            DO jk = 1, nlay_i
357               WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
358            END DO
359            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
360         ENDIF
361
362         smv_i(:,:,:) = smv_i(:,:,:)       + &
363            d_smv_i_thd(:,:,:) + &
364            d_smv_i_trp(:,:,:)
365
366         IF( ln_nicep ) THEN 
367            WRITE(numout,*) ' After advection   '
368            WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
369            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
370         ENDIF
371
372      ENDIF ! num_sal .EQ. 2
373
374      CALL lim_var_glo2eqv
375
376      !--------------------------------------
377      ! 2. Review of all pathological cases
378      !--------------------------------------
379
380      zrtt          = 173.15 * rone
381      zacrith       = 1.0e-6
382
383      !-------------------------------------------
384      ! 2.1) Advection of ice in an ice-free cell
385      !-------------------------------------------
386      ! should be removed since it is treated after dynamics now
387
388      zhimax = 5.0
389      ! first category
390      DO jj = 1, jpj
391         DO ji = 1, jpi
392            !--- the thickness of such an ice is often out of bounds
393            !--- thus we recompute a new area while conserving ice volume
394            zat_i_old = SUM(old_a_i(ji,jj,:))
395            zindb          =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
396            IF (      ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb.GT.zhimax) &
397               .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb).GT.zhimax ) &
398               .AND.( zat_i_old.LT.zacrith ) )  THEN ! new line
399               z_prescr_hi      = hi_max(1) / 2.0
400               a_i(ji,jj,1)     = v_i(ji,jj,1) / z_prescr_hi
401            ENDIF
402         END DO
403      END DO
404
405      IF( ln_nicep ) THEN 
406         WRITE(numout,*) ' 2.1 '
407         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
408         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
409         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
410         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
411         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
412         DO jk = 1, nlay_i
413            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
414         END DO
415      ENDIF
416
417      !change this 14h44
418      zhimax = 20.0 ! line added up
419      ! change this also 17 aug
420      zhimax = 30.0 ! line added up
421
422      DO jl = 2, jpl
423         jm = ice_types(jl)
424         DO jj = 1, jpj
425            DO ji = 1, jpi
426               zindb          =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
427               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
428               ! it makes problems when the advected volume and concentration do not seem to be
429               ! related with each other
430               ! the new thickness is sometimes very big!
431               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
432               ! which of course is plausible
433               ! but fuck! it fucks everything up :)
434               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb.GT.zhimax) &
435                  .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb).GT.zhimax ) THEN
436                  z_prescr_hi  =  ( hi_max_typ(jl-ice_cat_bounds(jm,1)  ,jm) + &
437                     hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / &
438                     2.0
439                  a_i(ji,jj,jl) = v_i(ji,jj,jl) / z_prescr_hi
440                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
441               ENDIF
442               zat_i_old = SUM(old_a_i(ji,jj,:))
443
444            END DO ! ji
445         END DO !jj
446      END DO !jl
447
448      IF( ln_nicep ) THEN 
449         WRITE(numout,*) ' 2.1 initial '
450         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
451         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
452         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
453         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
454         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
455         DO jk = 1, nlay_i
456            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
457         END DO
458      ENDIF
459
460      at_i(:,:) = 0.0
461      DO jl = 1, jpl
462         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
463      END DO
464
465      !----------------------------------------------------
466      ! 2.2) Rebin categories with thickness out of bounds
467      !----------------------------------------------------
468      IF( ln_nicep ) THEN 
469         WRITE(numout,*) ' 2.1 before rebinning '
470         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
471         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
472         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
473         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
474         DO jk = 1, nlay_i
475            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
476         END DO
477         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
478      ENDIF
479
480      DO jm = 1, jpm
481         jbnd1 = ice_cat_bounds(jm,1)
482         jbnd2 = ice_cat_bounds(jm,2)
483         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
484      END DO
485
486
487      IF( ln_nicep ) THEN 
488         WRITE(numout,*) ' 2.1 after rebinning'
489         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
490         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
491         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
492         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
493         DO jk = 1, nlay_i
494            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
495            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
496         END DO
497         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
498      ENDIF
499
500      at_i(:,:) = 0.0
501      DO jl = 1, jpl
502         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
503      END DO
504
505      !---------------------------------
506      ! 2.3) Melt of an internal layer
507      !---------------------------------
508      internal_melt(:,:,:) = .false.
509
510      DO jl = 1, jpl
511         DO jk = 1, nlay_i
512            DO jj = 1, jpj 
513               DO ji = 1, jpi
514                  ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
515                  IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
516                     ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) .AND. &
517                     ( v_i(ji,jj,jl) .GT. 0.0 ) .AND. &
518                     ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
519                     !                    WRITE(numout,*) ' Internal layer melt : '
520                     !                    WRITE(numout,*) ' ji, jj, jk, jl : ', ji,jj,jk,jl
521                     !                    WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
522                     !                    WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
523                     internal_melt(ji,jj,jl) = .true.
524                  ENDIF
525               END DO ! ji
526            END DO ! jj
527         END DO !jk
528      END DO !jl
529
530      DO jl = 1, jpl
531         DO jj = 1, jpj 
532            DO ji = 1, jpi
533               IF ( internal_melt(ji,jj,jl) ) THEN
534                  ! initial ice thickness
535                  !-----------------------
536                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
537                  !             WRITE(numout,*) ' ji,jj,jl : ', ji,jj,jl
538                  !             WRITE(numout,*) ' old ht_i: ', ht_i(ji,jj,jl)
539                  !             WRITE(numout,*) ' Enthalpy at the beg : ', e_i(ji,jj,1:nlay_i,jl)
540                  !             WRITE(numout,*) ' smv_i   : ', smv_i(ji,jj,jl)
541
542                  ! reduce ice thickness
543                  !-----------------------
544                  ind_im = 0
545                  zesum = 0.0
546                  DO jk = 1, nlay_i
547                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
548                     IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR.  & 
549                        ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
550                        ind_im = ind_im + 1
551                     zesum = zesum + e_i(ji,jj,jk,jl)
552                  END DO
553                  IF (ind_im .LT.nlay_i ) smv_i(ji,jj,jl)= smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * & 
554                     ( ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i )
555                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i
556                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
557
558                  !             WRITE(numout,*) ' ind_im  : ', ind_im
559                  !             WRITE(numout,*) ' new ht_i: ', ht_i(ji,jj,jl)
560                  !             WRITE(numout,*) ' smv_i   : ', smv_i(ji,jj,jl)
561                  !             WRITE(numout,*) ' zesum   : ', zesum
562
563                  ! redistribute heat
564                  !-----------------------
565                  ! old thicknesses and enthalpies
566                  ind_im = 0
567                  DO jk = 1, nlay_i
568                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
569                     IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND.  & 
570                        ( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
571                        ind_im = ind_im + 1
572                        zthick0(ind_im) = ht_i(ji,jj,jl) * ind_im / nlay_i
573                        zqm0   (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
574                     ENDIF
575                  END DO
576
577                  !             WRITE(numout,*) ' Old thickness, enthalpy '
578                  !             WRITE(numout,*) ' Number of layer : ind_im ', ind_im
579                  !             WRITE(numout,*) ' zthick0 : ', zthick0(1:ind_im)
580                  !             WRITE(numout,*) ' zqm0    : ', zqm0(1:ind_im)
581
582                  ! Redistributing energy on the new grid
583                  IF ( ind_im .GT. 0 ) THEN
584
585                     DO jk = 1, nlay_i
586                        e_i(ji,jj,jk,jl) = 0.0
587                        DO layer = 1, ind_im
588                           zweight         = MAX (  &
589                              MIN( ht_i(ji,jj,jl) * layer / ind_im , ht_i(ji,jj,jl) * jk / nlay_i ) -       &
590                              MAX( ht_i(ji,jj,jl) * (layer-1) / ind_im , ht_i(ji,jj,jl) * (jk-1) / nlay_i ) , 0.0 ) &
591                              /  ( ht_i(ji,jj,jl) / ind_im )
592
593                           e_i(ji,jj,jk,jl) =  e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
594                        END DO !layer
595                     END DO ! jk
596
597                     zesum = 0.0
598                     DO jk = 1, nlay_i
599                        zesum = zesum + e_i(ji,jj,jk,jl)
600                     END DO
601
602                     !             WRITE(numout,*) ' Enthalpy at the end : ', e_i(ji,jj,1:nlay_i,jl)
603                     !             WRITE(numout,*) ' Volume   at the end : ', v_i(ji,jj,jl)
604                     !             WRITE(numout,*) ' zesum : ', zesum
605
606                  ELSE ! ind_im .EQ. 0, total melt
607                     e_i(ji,jj,jk,jl) = 0.0
608                  ENDIF
609
610               ENDIF ! internal_melt
611
612            END DO ! ji
613         END DO !jj
614      END DO !jl
615      IF( ln_nicep ) THEN 
616         WRITE(numout,*) ' 2.3 after melt of an internal ice layer '
617         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
618         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
619         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
620         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
621         DO jk = 1, nlay_i
622            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
623            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
624         END DO
625         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
626      ENDIF
627
628      internal_melt(:,:,:) = .false.
629
630      ! Melt of snow
631      !--------------
632      DO jl = 1, jpl
633         DO jj = 1, jpj 
634            DO ji = 1, jpi
635               ! snow energy of melting
636               ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) /              &
637                  MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting
638
639               ! If snow energy of melting smaller then Lf
640               ! Then all snow melts and meltwater, heat go to the ocean
641               IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = .true.
642
643               IF( ln_nicep ) THEN 
644                  IF ( (ji.eq.jiindx) .AND. (jj.eq.jjindx) ) THEN
645                     WRITE(numout,*) ' jl    : ', jl
646                     WRITE(numout,*) ' ze_s  : ', ze_s
647                     WRITE(numout,*) ' v_s   : ', v_s(ji,jj,jl)
648                     WRITE(numout,*) ' rhosn : ', rhosn
649                     WRITE(numout,*) ' rhosn : ', lfus 
650                     WRITE(numout,*) ' area  : ', area(ji,jj)
651                     WRITE(numout,*) ' rhosn * lfus : ', rhosn * lfus 
652                     WRITE(numout,*) ' internal_melt : ', internal_melt(ji,jj,jl)
653                  ENDIF
654               ENDIF
655
656            END DO
657         END DO
658      END DO
659
660      DO jl = 1, jpl
661         DO jj = 1, jpj 
662            DO ji = 1, jpi
663               IF ( internal_melt(ji,jj,jl) ) THEN
664                  v_s(ji,jj,jl)   = 0.0
665                  e_s(ji,jj,1,jl) = 0.0
666                  !   ! release heat
667                  fheat_res(ji,jj) = fheat_res(ji,jj)  &
668                     + ze_s * v_s(ji,jj,jl) / rdt_ice
669                  ! release mass
670                  rdmsnif(ji,jj) =  rdmsnif(ji,jj) + rhosn * v_s(ji,jj,jl)
671               ENDIF
672            END DO
673         END DO
674      END DO
675
676      zbigvalue      = 1.0d+20
677
678      DO jl = 1, jpl
679         DO jj = 1, jpj 
680            DO ji = 1, jpi
681
682               !switches
683               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
684               !switch = 1 if a_i > 1e-06 and 0 if not
685               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
686               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
687               ! bug fix 25 avril 2007
688               zindb         = zindb*zindic
689
690               !--- 2.3 Correction to ice age
691               !------------------------------
692               !                IF ((o_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*float(numit))) THEN
693               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
694               !                ENDIF
695               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
696                  oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl)
697               ENDIF
698               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
699
700               !--- 2.4 Correction to snow thickness
701               !-------------------------------------
702               !          ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
703               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)
704               ! snow thickness cannot be smaller than 1e-6
705               v_s(ji,jj,jl)  = zindsn*v_s(ji,jj,jl)*zindb
706               v_s(ji,jj,jl)  = v_s(ji,jj,jl) *  MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) )
707
708               !--- 2.5 Correction to ice thickness
709               !-------------------------------------
710               ! ice thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hi = 0
711               v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl), 0.0)
712               ! ice thickness cannot be smaller than 1e-3
713               v_i(ji,jj,jl)  = zindic*v_i(ji,jj,jl)
714
715               !--- 2.6 Snow is transformed into ice if the original ice cover disappears
716               !----------------------------------------------------------------------------
717               zindg          = tms(ji,jj) *  MAX( rzero , SIGN( rone , -v_i(ji,jj,jl) ) )
718               v_i(ji,jj,jl)  = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0
719               v_s(ji,jj,jl)  = ( rone - zindg ) * v_s(ji,jj,jl) + & 
720                  zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
721
722               !--- 2.7 Correction to ice concentrations
723               !--------------------------------------------
724               ! if greater than 0, ice concentration cannot be smaller than 1e-10
725               a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
726               ! then ice volume has to be corrected too...
727               ! instead, zap small areas
728
729               !-------------------------
730               ! 2.8) Snow heat content
731               !-------------------------
732
733               e_s(ji,jj,1,jl) = zindsn *                                &
734                  ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) ) + &
735                  ( 1.0 - zindsn ) * 0.0
736
737            END DO ! ji
738         END DO ! jj
739      END DO ! jl
740
741      IF( ln_nicep ) THEN 
742         WRITE(numout,*) ' 2.8 '
743         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
744         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
745         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
746         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
747         DO jk = 1, nlay_i
748            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
749         END DO
750         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
751      ENDIF
752
753      !------------------------
754      ! 2.9) Ice heat content
755      !------------------------
756
757      DO jl = 1, jpl
758         DO jk = 1, nlay_i
759            DO jj = 1, jpj 
760               DO ji = 1, jpi
761                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi06 ) ) 
762                  ! =1 if v_i > 1e-6 and 0 if not
763                  e_i(ji,jj,jk,jl)= zindic * & 
764                     ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) + &
765                     ( 1.0 - zindic ) * 0.0
766               END DO ! ji
767            END DO ! jj
768         END DO !jk
769      END DO !jl
770
771      IF( ln_nicep ) THEN 
772         WRITE(numout,*) ' 2.9 '
773         DO jk = 1, nlay_i
774            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
775         END DO
776         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
777
778         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
779      ENDIF
780
781      !---------------------
782      ! 2.11) Ice salinity
783      !---------------------
784
785      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
786
787         DO jl = 1, jpl
788            DO jk = 1, nlay_i
789               DO jj = 1, jpj 
790                  DO ji = 1, jpi
791                     ! salinity stays in bounds
792                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), &
793                        0.1 * v_i(ji,jj,jl) )
794                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
795                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + &
796                        0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
797                  END DO ! ji
798               END DO ! jj
799            END DO !jk
800         END DO !jl
801
802      ENDIF
803
804      IF( ln_nicep ) THEN 
805         WRITE(numout,*) ' 2.11 '
806         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
807         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
808         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
809         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
810         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
811      ENDIF
812
813      DO jm = 1, jpm
814         DO jj = 1, jpj 
815            DO ji = 1, jpi
816               jl = ice_cat_bounds(jm,1)
817               !--- 2.12 Constrain the thickness of the smallest category above 5 cm
818               !----------------------------------------------------------------------
819               ! the ice thickness of the smallest category should be higher than 5 cm
820               ! we changed hiclim to 10
821               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
822               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
823               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) )
824               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
825               !             v_s(ji,jj,jl)  = v_s(ji,jj,jl) * zh
826               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
827               !             v_i(ji,jj,jl)  = v_i(ji,jj,jl) * zh
828               a_i (ji,jj,jl) = a_i(ji,jj,jl) /zh
829            END DO !ji
830         END DO !jj
831      END DO !jm
832      IF( ln_nicep ) THEN 
833         WRITE(numout,*) ' 2.12 '
834         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
835         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
836         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
837         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
838         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
839      ENDIF
840
841      !--- 2.13 Total ice concentration should not exceed 1
842      !-----------------------------------------------------
843      zamax = amax
844      ! 2.13.1) individual concentrations cannot exceed zamax
845      !------------------------------------------------------
846
847      at_i(:,:) = 0.0
848      DO jl = 1, jpl
849         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
850      END DO
851
852      ! 2.13.2) Total ice concentration cannot exceed zamax
853      !----------------------------------------------------
854      at_i(:,:) = 0.0
855      DO jl = 1, jpl
856         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
857      END DO
858
859      DO jj = 1, jpj
860         DO ji = 1, jpi
861
862            ! 0) Excessive area ?
863            z_da_ex =  MAX( at_i(ji,jj) - zamax , 0.0 )       
864
865            ! 1) Count the number of existing categories
866            DO jl  = 1, jpl
867               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) ) 
868               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) ) 
869               z_da_i(jl) = a_i(ji,jj,jl)*zindb*z_da_ex/MAX(at_i(ji,jj),epsi06)
870               z_dv_i(jl) = v_i(ji,jj,jl)*z_da_i(jl)/MAX(at_i(ji,jj),epsi06)
871               a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i(jl)
872               v_i(ji,jj,jl) = v_i(ji,jj,jl) + z_dv_i(jl)
873
874            END DO
875
876         END DO !ji
877      END DO !jj
878
879      IF( ln_nicep ) THEN 
880         WRITE(numout,*) ' 2.13 '
881         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
882         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
883         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
884         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
885         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
886      ENDIF
887
888      at_i(:,:) = 0.0
889      DO jl = 1, jpl
890         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
891      END DO
892
893      IF( ln_nicep ) THEN 
894         DO jj = 1, jpj
895            DO ji = 1, jpi
896               IF (at_i(ji,jj).GT.1.0) THEN
897                  WRITE(numout,*) ' lim_update ! : at_i > 1 -> PAS BIEN -> ALERTE '
898                  WRITE(numout,*) ' ~~~~~~~~~~   at_i     ', at_i(ji,jj)
899                  WRITE(numout,*) ' Point ', ji, jj
900                  WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
901                  DO jl = 1, jpl
902                     WRITE(numout,*) ' a_i ***         ', a_i(ji,jj,jl), ' CAT no ', jl
903                     WRITE(numout,*) ' a_i_old ***     ', old_a_i(ji,jj,jl), ' CAT no ', jl
904                     WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
905                  END DO
906                  !             WRITE(numout,*) ' CORRECTION BARBARE '
907                  !             z_da_ex =  MAX( at_i(ji,jj) - zamax , 0.0 )       
908               ENDIF
909            END DO
910         END DO
911      ENDIF
912
913      ! Final thickness distribution rebinning
914      ! --------------------------------------
915      IF( ln_nicep ) THEN 
916         WRITE(numout,*) ' rebinning before'
917         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
918         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
919         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
920         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
921         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
922      ENDIF
923      !old version
924      !    CALL lim_itd_th_reb(1,jpl)
925
926      DO jm = 1, jpm
927         jbnd1 = ice_cat_bounds(jm,1)
928         jbnd2 = ice_cat_bounds(jm,2)
929         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
930         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
931         ENDIF
932      END DO
933
934      IF( ln_nicep ) THEN 
935         WRITE(numout,*) ' rebinning final'
936         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
937         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
938         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
939         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
940         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
941      ENDIF
942
943      at_i(:,:) = 0.0
944      DO jl = 1, jpl
945         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
946      END DO
947
948      !------------------------------------------------------------------------------
949      ! 2) Corrections to avoid wrong values                                        |
950      !------------------------------------------------------------------------------
951      ! Ice drift
952      !------------
953
954      DO jj = 2, jpjm1
955         DO ji = fs_2, fs_jpim1
956            IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
957               IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side
958               IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
959               IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side
960               IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side
961            ENDIF
962         END DO
963      END DO
964      !mask velocities
965      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
966      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
967      !lateral boundary conditions
968      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
969      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
970
971      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
972      ! ALERTES
973      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
974
975      IF( ln_nicep ) THEN 
976         DO jj = 1, jpj
977            DO ji = 1, jpi
978               DO jl = 1, jpl
979                  !                IF ((v_i(ji,jj,jl).NE.0.0).AND.(a_i(ji,jj,jl).EQ.0.0)) THEN
980                  !                   WRITE(numout,*) ' lim_update : incompatible volume and concentration '
981               END DO ! jl
982
983               DO jl = 1, jpl
984                  IF ( (a_i(ji,jj,jl).GT.1.0).OR.(at_i(ji,jj).GT.1.0) ) THEN
985                     zindb          =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
986                     WRITE(numout,*) ' lim_update : a_i > 1 '
987                     WRITE(numout,*) ' PAS BIEN ----> ALERTE !!! '
988                     WRITE(numout,*) ' ~~~~~~~~~~   at_i     ', at_i(ji,jj)
989                     WRITE(numout,*) ' Point - category', ji, jj, jl
990                     WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
991                     WRITE(numout,*) ' a_i *** a_i_old ', a_i(ji,jj,jl), old_a_i(ji,jj,jl)
992                     WRITE(numout,*) ' v_i *** v_i_old ', v_i(ji,jj,jl), old_v_i(ji,jj,jl)
993                     WRITE(numout,*) ' ht_i ***        ', v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi06)*zindb 
994                     WRITE(numout,*) ' hi_max(jl), hi_max(jl-1) ', hi_max(jl), hi_max(jl-1)
995                     WRITE(numout,*) ' d_v_i_thd / trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
996                     WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
997                  ENDIF
998               END DO
999
1000            END DO !jj
1001         END DO !ji
1002
1003         WRITE(numout,*) ' TESTOSC1 ', tio_u(jiindx,jjindx), tio_v(jiindx,jjindx)
1004         WRITE(numout,*) ' TESTOSC2 ', u_ice(jiindx,jjindx), v_ice(jiindx,jjindx)
1005         WRITE(numout,*) ' TESTOSC3 ', u_oce(jiindx,jjindx), v_oce(jiindx,jjindx)
1006         WRITE(numout,*) ' TESTOSC4 ', utau (jiindx,jjindx), vtau (jiindx,jjindx)
1007      ENDIF
1008
1009
1010      IF(ln_ctl) THEN   ! Control print
1011         CALL prt_ctl_info(' ')
1012         CALL prt_ctl_info(' - Cell values : ')
1013         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
1014         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update  : cell area   :')
1015         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update  : at_i        :')
1016         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update  : vt_i        :')
1017         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update  : vt_s        :')
1018         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update  : strength    :')
1019         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
1020         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
1021         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
1022
1023         DO jl = 1, jpl
1024            CALL prt_ctl_info(' ')
1025            CALL prt_ctl_info(' - Category : ', ivar1=jl)
1026            CALL prt_ctl_info('   ~~~~~~~~~~')
1027            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update  : ht_i        : ')
1028            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update  : ht_s        : ')
1029            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update  : t_su        : ')
1030            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update  : t_snow      : ')
1031            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update  : sm_i        : ')
1032            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update  : o_i         : ')
1033            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update  : a_i         : ')
1034            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update  : old_a_i     : ')
1035            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_a_i_trp   : ')
1036            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_a_i_thd   : ')
1037            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update  : v_i         : ')
1038            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update  : old_v_i     : ')
1039            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_v_i_trp   : ')
1040            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_v_i_thd   : ')
1041            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update  : v_s         : ')
1042            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update  : old_v_s     : ')
1043            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_v_s_trp   : ')
1044            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_v_s_thd   : ')
1045            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : e_i1        : ')
1046            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : old_e_i1    : ')
1047            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : de_i1_trp   : ')
1048            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : de_i1_thd   : ')
1049            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : e_i2        : ')
1050            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : old_e_i2    : ')
1051            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : de_i2_trp   : ')
1052            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : de_i2_thd   : ')
1053            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update  : e_snow      : ')
1054            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update  : old_e_snow  : ')
1055            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : d_e_s_trp   : ')
1056            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : d_e_s_thd   : ')
1057            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update  : smv_i       : ')
1058            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update  : old_smv_i   : ')
1059            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update  : d_smv_i_trp : ')
1060            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update  : d_smv_i_thd : ')
1061            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update  : oa_i        : ')
1062            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update  : old_oa_i    : ')
1063            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update  : d_oa_i_trp  : ')
1064            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update  : d_oa_i_thd  : ')
1065            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update  : Path. case  : ')
1066
1067            DO jk = 1, nlay_i
1068               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
1069               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update  : t_i       : ')
1070            END DO
1071         END DO
1072
1073         CALL prt_ctl_info(' ')
1074         CALL prt_ctl_info(' - Heat / FW fluxes : ')
1075         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
1076         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
1077         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
1078         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
1079
1080         CALL prt_ctl_info(' ')
1081         CALL prt_ctl_info(' - Stresses : ')
1082         CALL prt_ctl_info('   ~~~~~~~~~~ ')
1083         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
1084         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
1085         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
1086      ENDIF
1087
1088      !---------------------
1089
1090   END SUBROUTINE lim_update
1091#else
1092   !!----------------------------------------------------------------------
1093   !!   Default option         Empty Module               No sea-ice model
1094   !!----------------------------------------------------------------------
1095CONTAINS
1096   SUBROUTINE lim_update     ! Empty routine
1097   END SUBROUTINE lim_update
1098
1099#endif
1100
1101END MODULE limupdate
Note: See TracBrowser for help on using the repository browser.