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 @ 4332

Last change on this file since 4332 was 4332, checked in by clem, 10 years ago

update LIM3 to fix remaining bugs. Now working in global and regional config.

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