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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 12 years ago

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