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

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

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