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/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 3838

Last change on this file since 3838 was 3838, checked in by gm, 11 years ago

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