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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 10 years ago

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