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

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 4072

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

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