New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limupdate.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limupdate.F90 @ 867

Last change on this file since 867 was 867, checked in by ctlod, 16 years ago

Add control prints for sea-ice

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