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

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]
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 : ')
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.