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.
limupdate2.F90 in branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116
  • Property svn:executable set to *
File size: 31.5 KB
Line 
1MODULE limupdate2
2   !!======================================================================
3   !!                     ***  MODULE  limupdate2  ***
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_update2   : 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 limdiahsb
30   USE limwri
31   USE limrst
32   USE thd_ice         ! LIM thermodynamic sea-ice variables
33   USE par_ice
34   USE limitd_th
35   USE limitd_me
36   USE limvar
37   USE prtctl           ! Print control
38   USE lbclnk           ! lateral boundary condition - MPP exchanges
39   USE wrk_nemo         ! work arrays
40   USE lib_fortran     ! glob_sum
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   lim_update2   ! routine called by ice_step
46
47      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
48      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
49      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
50      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
51      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
52      REAL(wp)  ::   rzero  = 0._wp       !    -       -
53      REAL(wp)  ::   rone   = 1._wp       !    -       -
54         
55   !! * Substitutions
56#  include "vectopt_loop_substitute.h90"
57   !!----------------------------------------------------------------------
58   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
59   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE lim_update2
65      !!-------------------------------------------------------------------
66      !!               ***  ROUTINE lim_update2  ***
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  : 
74      !!    Ice speed from ice dynamics
75      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
76      !!      from advection and ice thermodynamics
77      !!
78      !! ** Action  : -
79      !!---------------------------------------------------------------------
80      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
81      INTEGER ::   jbnd1, jbnd2
82      INTEGER ::   i_ice_switch
83      INTEGER ::   ind_im, layer      ! indices for internal melt
84      REAL(wp) ::   zweight, zesum, zhimax, z_da_i
85      REAL(wp) ::   zinda, zindb, zindsn, zindic
86      REAL(wp) ::   zindg, zh, zdvres, zviold2
87      REAL(wp) ::   zbigvalue, zvsold2, z_da_ex
88      REAL(wp) ::   z_prescr_hi, zat_i_old, ztmelts, ze_s
89
90      INTEGER , POINTER, DIMENSION(:,:,:) ::  internal_melt
91      REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for
92      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
93      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
94      ! mass and salt flux (clem)
95      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
96      !!-------------------------------------------------------------------
97
98      CALL wrk_alloc( jpi,jpj,jpl, internal_melt )   ! integer
99      CALL wrk_alloc( jkmax, zthick0, zqm0 )
100
101      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
102
103      !------------------------------------------------------------------------------
104      ! 1. Update of Global variables                                               |
105      !------------------------------------------------------------------------------
106
107      !-----------------------------
108      ! Update ice and snow volumes 
109      !-----------------------------
110      DO jl = 1, jpl
111         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_thd(:,:,jl) 
112         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_thd(:,:,jl)
113      END DO
114 
115      !---------------------------------------------
116      ! Ice concentration and ice heat content
117      !---------------------------------------------
118      a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:)
119      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:) 
120 
121      !------------------------------
122      ! Snow temperature and ice age
123      !------------------------------
124      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:)
125      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_thd(:,:,:)
126
127      !--------------
128      ! Ice salinity   
129      !--------------
130      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case
131         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:)
132      ENDIF
133
134      ! mass and salt flux init (clem)
135      zviold(:,:,:) = v_i(:,:,:)
136      zvsold(:,:,:) = v_s(:,:,:)
137      zsmvold(:,:,:) = smv_i(:,:,:)
138
139      ! -------------------------------
140      !- check conservation (C Rousset)
141      IF (ln_limdiahsb) THEN
142         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
143         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
144         zchk_fw_b  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) )
145         zchk_fs_b  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) )
146      ENDIF
147      !- check conservation (C Rousset)
148      ! -------------------------------
149
150      CALL lim_var_glo2eqv
151
152      !---------------------------------
153      ! Classify the pathological cases
154      !---------------------------------
155      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
156      ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
157      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
158      ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation
159      ! with negative advection, very pathological )
160      !
161      DO jl = 1, jpl
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               patho_case(ji,jj,jl) = 1
165               IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN
166                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
167                     patho_case(ji,jj,jl) = 2
168                  ENDIF
169               ELSE
170                  patho_case(ji,jj,jl) = 3
171                  IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
172                     patho_case(ji,jj,jl) = 4
173                  ENDIF
174               ENDIF
175             END DO
176         END DO
177      END DO
178
179      !--------------------------------------
180      ! 2. Review of all pathological cases
181      !--------------------------------------
182      !-------------------------------------------
183      ! 2.1) Advection of ice in an ice-free cell
184      !-------------------------------------------
185      ! should be removed since it is treated after dynamics now
186
187!      !IF( ln_nicep ) THEN 
188!         WRITE(numout,*) ' limupdate2 - before h correction '
189!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
190!         WRITE(numout,*) ' at_i : ', at_i(12,44)
191!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
192!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
193!      !ENDIF
194
195      zhimax = 1._wp
196      ! first category
197      DO jj = 1, jpj
198         DO ji = 1, jpi
199            !--- the thickness of such an ice is often out of bounds
200            !--- thus we recompute a new area while conserving ice volume
201            zat_i_old = SUM( old_a_i(ji,jj,:) )
202            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) ) 
203            IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
204                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
205                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
206               ht_i(ji,jj,1) = hi_max(1) / 2.0
207               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
208            ENDIF
209         END DO
210      END DO
211
212!      !IF( ln_nicep ) THEN 
213!         at_i(:,:) = 0._wp
214!         DO jl = 1, jpl
215!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
216!         END DO
217!         WRITE(numout,*) ' limupdate2 - after h correction 1 '
218!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
219!         WRITE(numout,*) ' at_i : ', at_i(12,44)
220!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
221!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
222!      !ENDIF
223!
224      zhimax = 10._wp
225      ! other categories
226      DO jl = 2, jpl
227         jm = ice_types(jl)
228         DO jj = 1, jpj
229            DO ji = 1, jpi
230               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 
231              ! this correction is very tricky... sometimes, advection gets wrong i don't know why
232               ! it makes problems when the advected volume and concentration do not seem to be
233               ! related with each other
234               ! the new thickness is sometimes very big!
235               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
236               ! which of course is plausible
237               ! but fuck! it fucks everything up :)
238               IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
239                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
240                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
241                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
242               ENDIF
243            END DO ! ji
244         END DO !jj
245      END DO !jl
246     
247      at_i(:,:) = 0._wp
248      DO jl = 1, jpl
249         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
250      END DO
251
252!      !IF( ln_nicep ) THEN 
253!         WRITE(numout,*) ' limupdate2 - after h correction 2 '
254!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
255!         WRITE(numout,*) ' at_i : ', at_i(12,44)
256!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
257!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
258!      !ENDIF
259
260      !----------------------------------------------------
261      ! 2.2) Rebin categories with thickness out of bounds
262      !----------------------------------------------------
263      DO jm = 1, jpm
264         jbnd1 = ice_cat_bounds(jm,1)
265         jbnd2 = ice_cat_bounds(jm,2)
266         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
267      END DO
268
269      !---------------------------------
270      ! 2.3) Melt of an internal layer
271      !---------------------------------
272      internal_melt(:,:,:) = 0
273
274      DO jl = 1, jpl
275         DO jk = 1, nlay_i
276            DO jj = 1, jpj 
277               DO ji = 1, jpi
278                  ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
279                  IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
280                    & .AND. ( v_i(ji,jj,jl) .GT. 0.0 ) .AND. ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
281                     internal_melt(ji,jj,jl) = 1
282                  ENDIF
283               END DO ! ji
284            END DO ! jj
285         END DO !jk
286      END DO !jl
287
288      DO jl = 1, jpl
289         DO jj = 1, jpj 
290            DO ji = 1, jpi
291               IF( internal_melt(ji,jj,jl) == 1 ) THEN
292                  ! initial ice thickness
293                  !-----------------------
294                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
295
296                  ! reduce ice thickness
297                  !-----------------------
298                  ind_im = 0
299                  zesum = 0.0
300                  DO jk = 1, nlay_i
301                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
302                     IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) ind_im = ind_im + 1
303                     zesum = zesum + e_i(ji,jj,jk,jl)
304                  END DO
305                  IF (ind_im .LT. nlay_i ) THEN
306                     smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )
307                  ENDIF
308                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i)
309                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
310
311                  !CLEM
312                  zdvres = REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) * a_i(ji,jj,jl)
313                  !rdmicif(ji,jj) = rdmicif(ji,jj) - zdvres * rhoic
314                  !fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice )
315
316                  ! redistribute heat
317                  !-----------------------
318                  ! old thicknesses and enthalpies
319                  ind_im = 0
320                  DO jk = 1, nlay_i
321                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
322                     IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND.  & 
323                        ( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
324                        ind_im = ind_im + 1
325                        zthick0(ind_im) = ht_i(ji,jj,jl) * REAL(ind_im / nlay_i)
326                        zqm0   (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
327                     ENDIF
328                  END DO
329
330                  ! Redistributing energy on the new grid
331                  IF ( ind_im .GT. 0 ) THEN
332
333                     DO jk = 1, nlay_i
334                        e_i(ji,jj,jk,jl) = 0.0
335                        DO layer = 1, ind_im
336                           zweight         = MAX (  &
337                              MIN( ht_i(ji,jj,jl) * REAL(layer/ind_im) , ht_i(ji,jj,jl) * REAL(jk / nlay_i) ) -       &
338                              MAX( ht_i(ji,jj,jl) * REAL((layer-1)/ind_im) , ht_i(ji,jj,jl) * REAL((jk-1) / nlay_i) ) , 0.0 ) &
339                              /  ( ht_i(ji,jj,jl) / REAL(ind_im) )
340
341                           e_i(ji,jj,jk,jl) =  e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
342                        END DO !layer
343                     END DO ! jk
344
345                     zesum = 0.0
346                     DO jk = 1, nlay_i
347                        zesum = zesum + e_i(ji,jj,jk,jl)
348                     END DO
349
350                  ELSE ! ind_im .EQ. 0, total melt
351                     e_i(ji,jj,jk,jl) = 0.0
352                  ENDIF
353
354               ENDIF ! internal_melt
355
356            END DO ! ji
357         END DO !jj
358      END DO !jl
359
360      internal_melt(:,:,:) = 0
361
362
363      ! Melt of snow
364      !--------------
365      DO jl = 1, jpl
366         DO jj = 1, jpj 
367            DO ji = 1, jpi
368               ! snow energy of melting
369               ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting
370
371               ! If snow energy of melting smaller then Lf
372               ! Then all snow melts and meltwater, heat go to the ocean
373               IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = 1
374
375            END DO
376         END DO
377      END DO
378
379      DO jl = 1, jpl
380         DO jj = 1, jpj 
381            DO ji = 1, jpi
382               IF ( internal_melt(ji,jj,jl) == 1 ) THEN
383                  zdvres = v_s(ji,jj,jl)
384                  ! release heat
385                  fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s * zdvres / rdt_ice
386                  ! release mass
387                  !rdmsnif(ji,jj) =  rdmsnif(ji,jj) - zdvres * rhosn
388                  !
389                  v_s(ji,jj,jl)   = 0.0
390                  e_s(ji,jj,1,jl) = 0.0
391                 ENDIF
392            END DO
393         END DO
394      END DO
395
396      zbigvalue      = 1.0e+20
397      DO jl = 1, jpl
398         DO jj = 1, jpj 
399            DO ji = 1, jpi
400
401               !switches
402               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
403               !switch = 1 if a_i > 1e-06 and 0 if not
404               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
405               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
406               ! bug fix 25 avril 2007
407               zindb         = zindb*zindic
408
409               !--- 2.3 Correction to ice age
410               !------------------------------
411               !                IF ((o_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*float(numit))) THEN
412               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
413               !                ENDIF
414               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
415                  oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl)
416               ENDIF
417               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
418
419               !--- 2.4 Correction to snow thickness
420               !-------------------------------------
421               zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl)
422               v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres
423
424               !rdmsnif(ji,jj) = rdmsnif(ji,jj) + zdvres * rhosn
425 
426               !--- 2.5 Correction to ice thickness
427               !-------------------------------------
428               zdvres = (zindb - 1._wp) * v_i(ji,jj,jl)
429               v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres
430
431               !rdmicif(ji,jj) = rdmicif(ji,jj) + zdvres * rhoic
432               !fsalt_res(ji,jj)  = fsalt_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice )
433
434               !--- 2.6 Snow is transformed into ice if the original ice cover disappears
435               !----------------------------------------------------------------------------
436               zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
437               zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0
438               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres
439
440               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn )
441               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres
442
443               !--- 2.7 Correction to ice concentrations
444               !--------------------------------------------
445               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
446               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
447
448               !-------------------------
449               ! 2.8) Snow heat content
450               !-------------------------
451               e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) )
452
453            END DO ! ji
454         END DO ! jj
455      END DO ! jl
456
457      !------------------------
458      ! 2.9) Ice heat content
459      !------------------------
460
461      DO jl = 1, jpl
462         DO jk = 1, nlay_i
463            DO jj = 1, jpj 
464               DO ji = 1, jpi
465                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
466                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
467               END DO ! ji
468            END DO ! jj
469         END DO !jk
470      END DO !jl
471
472
473      DO jm = 1, jpm
474         DO jj = 1, jpj 
475            DO ji = 1, jpi
476               jl = ice_cat_bounds(jm,1)
477               !--- 2.12 Constrain the thickness of the smallest category above 5 cm
478               !----------------------------------------------------------------------
479               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
480               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
481               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) )
482               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
483               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
484               a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh
485               !CLEM
486               v_i (ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl)
487               v_s (ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl)
488            END DO !ji
489         END DO !jj
490      END DO !jm
491
492      at_i(:,:) = 0.0
493      DO jl = 1, jpl
494         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
495      END DO
496     
497      !--- 2.13 ice concentration should not exceed amax
498      !         (it should not be the case)
499      !-----------------------------------------------------
500      DO jj = 1, jpj
501         DO ji = 1, jpi
502            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )       
503            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
504            DO jl  = 1, jpl
505               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
506               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
507               !
508               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
509               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
510               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
511            END DO
512         END DO
513      END DO
514      at_i(:,:) = 0.0
515      DO jl = 1, jpl
516         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
517      END DO
518
519      ! Final thickness distribution rebinning
520      ! --------------------------------------
521      DO jm = 1, jpm
522         jbnd1 = ice_cat_bounds(jm,1)
523         jbnd2 = ice_cat_bounds(jm,2)
524         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
525         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
526         ENDIF
527      END DO
528
529      !---------------------
530      ! 2.11) Ice salinity
531      !---------------------
532      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
533      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
534         DO jl = 1, jpl
535            !DO jk = 1, nlay_i
536               DO jj = 1, jpj 
537                  DO ji = 1, jpi
538                     ! salinity stays in bounds
539                     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) )
540                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
541                     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)
542                  END DO ! ji
543               END DO ! jj
544            !END DO !jk
545         END DO !jl
546      ENDIF
547
548      ! -------------------
549      at_i(:,:) = a_i(:,:,1)
550      DO jl = 2, jpl
551         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
552      END DO
553
554      !------------------------------------------------------------------------------
555      ! 2) Corrections to avoid wrong values                                        |
556      !------------------------------------------------------------------------------
557      ! Ice drift
558      !------------
559      DO jj = 2, jpjm1
560         DO ji = fs_2, fs_jpim1
561            IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
562               IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side
563               IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
564               IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side
565               IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side
566            ENDIF
567         END DO
568      END DO
569      !mask velocities
570      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
571      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
572      !lateral boundary conditions
573      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
574      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
575 
576      !--------------------------------
577      ! Update mass/salt fluxes (clem)
578      !--------------------------------
579      DO jl = 1, jpl
580         DO jj = 1, jpj 
581            DO ji = 1, jpi
582               diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 
583               rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
584               rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
585               fsalt_res(ji,jj) = fsalt_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 
586            END DO
587         END DO
588      END DO
589
590      ! -------------------------------
591      !- check conservation (C Rousset)
592      IF (ln_limdiahsb) THEN
593         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
594         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
595 
596         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
597         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
598
599         zchk_vmin = glob_min(v_i)
600         zchk_amax = glob_max(SUM(a_i,dim=3))
601         zchk_amin = glob_min(a_i)
602
603         IF(lwp) THEN
604            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.)
605            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.)
606            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(zchk_vmin * 1.e-3)
607            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate2) = ',zchk_amax
608            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',zchk_amin
609         ENDIF
610      ENDIF
611      !- check conservation (C Rousset)
612      ! -------------------------------
613
614      IF(ln_ctl) THEN   ! Control print
615         CALL prt_ctl_info(' ')
616         CALL prt_ctl_info(' - Cell values : ')
617         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
618         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :')
619         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
620         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
621         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
622         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
623         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
624         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update2  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
625
626         DO jl = 1, jpl
627            CALL prt_ctl_info(' ')
628            CALL prt_ctl_info(' - Category : ', ivar1=jl)
629            CALL prt_ctl_info('   ~~~~~~~~~~')
630            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
631            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
632            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
633            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
634            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
635            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
636            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
637            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_a_i     : ')
638            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
639            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
640            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_i     : ')
641            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
642            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
643            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_s     : ')
644            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
645            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
646            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i1    : ')
647            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
648            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
649            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i2    : ')
650            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
651            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
652            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update2  : old_e_snow  : ')
653            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
654            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
655            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update2  : old_smv_i   : ')
656            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
657            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
658            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ')
659            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
660            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update2  : Path. case  : ')
661
662            DO jk = 1, nlay_i
663               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
664               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
665            END DO
666         END DO
667
668         CALL prt_ctl_info(' ')
669         CALL prt_ctl_info(' - Heat / FW fluxes : ')
670         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
671         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
672         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
673         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
674
675         CALL prt_ctl_info(' ')
676         CALL prt_ctl_info(' - Stresses : ')
677         CALL prt_ctl_info('   ~~~~~~~~~~ ')
678         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
679         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
680         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
681      ENDIF
682   
683      CALL wrk_dealloc( jpi,jpj,jpl, internal_melt )   ! integer
684      CALL wrk_dealloc( jkmax, zthick0, zqm0 )
685
686      CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
687
688   END SUBROUTINE lim_update2
689#else
690   !!----------------------------------------------------------------------
691   !!   Default option         Empty Module               No sea-ice model
692   !!----------------------------------------------------------------------
693CONTAINS
694   SUBROUTINE lim_update2     ! Empty routine
695   END SUBROUTINE lim_update2
696
697#endif
698
699END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.