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 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

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