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 branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 2833

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

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