New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limupdate2.F90 in branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 3938 was 3938, checked in by flavoni, 11 years ago

dev_r3406_CNRS_LIM3: update LIM3, see ticket #1116
  • Property svn:executable set to *
File size: 27.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 limdia
29   USE limdiahsb
30   USE limwri
31   USE limrst
32   USE thd_ice         ! LIM thermodynamic sea-ice variables
33   USE par_ice
34   USE limitd_th
35   USE limitd_me
36   USE limvar
37   USE prtctl           ! Print control
38   USE lbclnk           ! lateral boundary condition - MPP exchanges
39   USE wrk_nemo         ! work arrays
40   USE lib_fortran     ! glob_sum
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   lim_update2   ! routine called by ice_step
46
47      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
48      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
49      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
50      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
51      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
52      REAL(wp)  ::   rzero  = 0._wp       !    -       -
53      REAL(wp)  ::   rone   = 1._wp       !    -       -
54         
55   !! * Substitutions
56#  include "vectopt_loop_substitute.h90"
57   !!----------------------------------------------------------------------
58   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
59   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE lim_update2
65      !!-------------------------------------------------------------------
66      !!               ***  ROUTINE lim_update2  ***
67      !!               
68      !! ** Purpose :  Computes update of sea-ice global variables at
69      !!               the end of the time step.
70      !!               Address pathological cases
71      !!               This place is very important
72      !!               
73      !! ** Method  : 
74      !!    Ice speed from ice dynamics
75      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
76      !!      from advection and ice thermodynamics
77      !!
78      !! ** Action  : -
79      !!---------------------------------------------------------------------
80      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
81      INTEGER ::   jbnd1, jbnd2
82      INTEGER ::   i_ice_switch
83      INTEGER ::   ind_im, layer      ! indices for internal melt
84      REAL(wp) ::   zweight, zesum, zhimax, z_da_i
85      REAL(wp) ::   zindb, zindsn, zindic
86      REAL(wp) ::   zindg, zh, zdvres, zviold2
87      REAL(wp) ::   zbigvalue, zvsold2, z_da_ex
88      REAL(wp) ::   z_prescr_hi, zat_i_old, ztmelts, ze_s
89
90      INTEGER , POINTER, DIMENSION(:,:,:) ::  internal_melt
91      REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for
92      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
93
94      ! mass and salt flux (clem)
95      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
96      !!-------------------------------------------------------------------
97
98      CALL wrk_alloc( jpi,jpj,jpl, internal_melt )   ! integer
99      CALL wrk_alloc( jkmax, zthick0, zqm0 )
100
101      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
102
103      !------------------------------------------------------------------------------
104      ! 1. Update of Global variables                                               |
105      !------------------------------------------------------------------------------
106
107      !-----------------------------
108      ! Update ice and snow volumes 
109      !-----------------------------
110      DO jl = 1, jpl
111         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_thd(:,:,jl) 
112         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_thd(:,:,jl)
113      END DO
114 
115      !---------------------------------------------
116      ! Ice concentration and ice heat content
117      !---------------------------------------------
118      a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:)
119      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:) 
120 
121      !------------------------------
122      ! Snow temperature and ice age
123      !------------------------------
124      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:)
125      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_thd(:,:,:)
126
127      !--------------
128      ! Ice salinity   
129      !--------------
130      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case
131         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:)
132      ENDIF
133
134      ! mass and salt flux init (clem)
135      zviold(:,:,:) = v_i(:,:,:)
136      zvsold(:,:,:) = v_s(:,:,:)
137      zsmvold(:,:,:) = smv_i(:,:,:)
138
139      ! -------------------------------
140      !- check conservation (C Rousset)
141      IF (ln_limdiahsb) THEN
142         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
143         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
144         zchk_fw_b  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) )
145         zchk_fs_b  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) )
146      ENDIF
147      !- check conservation (C Rousset)
148      ! -------------------------------
149
150      CALL lim_var_glo2eqv
151
152      !---------------------------------
153      ! Classify the pathological cases
154      !---------------------------------
155      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
156      ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
157      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
158      ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation
159      ! with negative advection, very pathological )
160      !
161      DO jl = 1, jpl
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               patho_case(ji,jj,jl) = 1
165               IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN
166                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
167                     patho_case(ji,jj,jl) = 2
168                  ENDIF
169               ELSE
170                  patho_case(ji,jj,jl) = 3
171                  IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
172                     patho_case(ji,jj,jl) = 4
173                  ENDIF
174               ENDIF
175             END DO
176         END DO
177      END DO
178
179
180      !--------------------------------------
181      ! 2. Review of all pathological cases
182      !--------------------------------------
183
184
185      !----------------------------------------------------
186      ! 2.2) Rebin categories with thickness out of bounds
187      !----------------------------------------------------
188      DO jm = 1, jpm
189         jbnd1 = ice_cat_bounds(jm,1)
190         jbnd2 = ice_cat_bounds(jm,2)
191         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
192      END DO
193
194      !---------------------------------
195      ! 2.3) Melt of an internal layer
196      !---------------------------------
197      internal_melt(:,:,:) = 0
198
199      DO jl = 1, jpl
200         DO jk = 1, nlay_i
201            DO jj = 1, jpj 
202               DO ji = 1, jpi
203                  ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
204                  IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
205                    & .AND. ( v_i(ji,jj,jl) .GT. 0.0 ) .AND. ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
206                     internal_melt(ji,jj,jl) = 1
207                  ENDIF
208               END DO ! ji
209            END DO ! jj
210         END DO !jk
211      END DO !jl
212
213      DO jl = 1, jpl
214         DO jj = 1, jpj 
215            DO ji = 1, jpi
216               IF( internal_melt(ji,jj,jl) == 1 ) THEN
217                  ! initial ice thickness
218                  !-----------------------
219                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
220
221                  ! reduce ice thickness
222                  !-----------------------
223                  ind_im = 0
224                  zesum = 0.0
225                  DO jk = 1, nlay_i
226                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
227                     IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) ind_im = ind_im + 1
228                     zesum = zesum + e_i(ji,jj,jk,jl)
229                  END DO
230                  IF (ind_im .LT. nlay_i ) THEN
231                     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) )
232                  ENDIF
233                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i)
234                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
235
236                  !CLEM
237                  zdvres = REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) * a_i(ji,jj,jl)
238                  !rdmicif(ji,jj) = rdmicif(ji,jj) - zdvres * rhoic
239                  !fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice )
240
241                  ! redistribute heat
242                  !-----------------------
243                  ! old thicknesses and enthalpies
244                  ind_im = 0
245                  DO jk = 1, nlay_i
246                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
247                     IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND.  & 
248                        ( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
249                        ind_im = ind_im + 1
250                        zthick0(ind_im) = ht_i(ji,jj,jl) * REAL(ind_im / nlay_i)
251                        zqm0   (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
252                     ENDIF
253                  END DO
254
255                  ! Redistributing energy on the new grid
256                  IF ( ind_im .GT. 0 ) THEN
257
258                     DO jk = 1, nlay_i
259                        e_i(ji,jj,jk,jl) = 0.0
260                        DO layer = 1, ind_im
261                           zweight         = MAX (  &
262                              MIN( ht_i(ji,jj,jl) * REAL(layer/ind_im) , ht_i(ji,jj,jl) * REAL(jk / nlay_i) ) -       &
263                              MAX( ht_i(ji,jj,jl) * REAL((layer-1)/ind_im) , ht_i(ji,jj,jl) * REAL((jk-1) / nlay_i) ) , 0.0 ) &
264                              /  ( ht_i(ji,jj,jl) / REAL(ind_im) )
265
266                           e_i(ji,jj,jk,jl) =  e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
267                        END DO !layer
268                     END DO ! jk
269
270                     zesum = 0.0
271                     DO jk = 1, nlay_i
272                        zesum = zesum + e_i(ji,jj,jk,jl)
273                     END DO
274
275                  ELSE ! ind_im .EQ. 0, total melt
276                     e_i(ji,jj,jk,jl) = 0.0
277                  ENDIF
278
279               ENDIF ! internal_melt
280
281            END DO ! ji
282         END DO !jj
283      END DO !jl
284
285      internal_melt(:,:,:) = 0
286
287
288      ! Melt of snow
289      !--------------
290      DO jl = 1, jpl
291         DO jj = 1, jpj 
292            DO ji = 1, jpi
293               ! snow energy of melting
294               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
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                  !rdmsnif(ji,jj) =  rdmsnif(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) - epsi06 ) ) 
328               !switch = 1 if a_i > 1e-06 and 0 if not
329               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
330               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 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)*86400.0.gt.(rdt_ice*float(numit))) THEN
337               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
338               !                ENDIF
339               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
340                  oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*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               !rdmsnif(ji,jj) = rdmsnif(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               !rdmicif(ji,jj) = rdmicif(ji,jj) + zdvres * rhoic
357               !fsalt_res(ji,jj)  = fsalt_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         = - 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               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
371               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
372
373               !-------------------------
374               ! 2.8) Snow heat content
375               !-------------------------
376               e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) )
377
378            END DO ! ji
379         END DO ! jj
380      END DO ! jl
381
382      !------------------------
383      ! 2.9) Ice heat content
384      !------------------------
385
386      DO jl = 1, jpl
387         DO jk = 1, nlay_i
388            DO jj = 1, jpj 
389               DO ji = 1, jpi
390                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
391                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
392               END DO ! ji
393            END DO ! jj
394         END DO !jk
395      END DO !jl
396
397
398      DO jm = 1, jpm
399         DO jj = 1, jpj 
400            DO ji = 1, jpi
401               jl = ice_cat_bounds(jm,1)
402               !--- 2.12 Constrain the thickness of the smallest category above 5 cm
403               !----------------------------------------------------------------------
404               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
405               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
406               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) )
407               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
408               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
409               a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh
410               !CLEM
411               v_i (ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl)
412               v_s (ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl)
413            END DO !ji
414         END DO !jj
415      END DO !jm
416
417      at_i(:,:) = 0.0
418      DO jl = 1, jpl
419         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
420      END DO
421     
422      !--- 2.13 ice concentration should not exceed amax
423      !-----------------------------------------------------
424      DO jj = 1, jpj
425         DO ji = 1, jpi
426            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )       
427            DO jl  = 1, jpl
428               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
429               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
430               a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i
431            END DO
432         END DO
433      END DO
434      at_i(:,:) = 0.0
435      DO jl = 1, jpl
436         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
437      END DO
438
439      ! Final thickness distribution rebinning
440      ! --------------------------------------
441      DO jm = 1, jpm
442         jbnd1 = ice_cat_bounds(jm,1)
443         jbnd2 = ice_cat_bounds(jm,2)
444         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
445         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
446         ENDIF
447      END DO
448
449      !---------------------
450      ! 2.11) Ice salinity
451      !---------------------
452      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
453      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
454         DO jl = 1, jpl
455            !DO jk = 1, nlay_i
456               DO jj = 1, jpj 
457                  DO ji = 1, jpi
458                     ! salinity stays in bounds
459                     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) )
460                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
461                     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)
462                  END DO ! ji
463               END DO ! jj
464            !END DO !jk
465         END DO !jl
466      ENDIF
467
468      ! -------------------
469      at_i(:,:) = a_i(:,:,1)
470      DO jl = 2, jpl
471         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
472      END DO
473
474      !------------------------------------------------------------------------------
475      ! 2) Corrections to avoid wrong values                                        |
476      !------------------------------------------------------------------------------
477      ! Ice drift
478      !------------
479      DO jj = 2, jpjm1
480         DO ji = fs_2, fs_jpim1
481            IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
482               IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side
483               IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
484               IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side
485               IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side
486            ENDIF
487         END DO
488      END DO
489      !mask velocities
490      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
491      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
492      !lateral boundary conditions
493      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
494      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
495 
496      !--------------------------------
497      ! Update mass/salt fluxes (clem)
498      !--------------------------------
499      DO jl = 1, jpl
500         DO jj = 1, jpj 
501            DO ji = 1, jpi
502               diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 
503               rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
504               rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
505               fsalt_res(ji,jj) = fsalt_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 
506            END DO
507         END DO
508      END DO
509
510      ! -------------------------------
511      !- check conservation (C Rousset)
512      IF (ln_limdiahsb) THEN
513         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
514         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
515 
516         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
517         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
518
519         IF(lwp) THEN
520            IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.)
521            IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.)
522            IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(MINVAL(v_i) * 1.e-3)
523            IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limupdate2) = ',MAXVAL(SUM(a_i,dim=3))
524            IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',MINVAL(a_i)
525         ENDIF
526      ENDIF
527      !- check conservation (C Rousset)
528      ! -------------------------------
529
530      IF(ln_ctl) THEN   ! Control print
531         CALL prt_ctl_info(' ')
532         CALL prt_ctl_info(' - Cell values : ')
533         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
534         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :')
535         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
536         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
537         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
538         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
539         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
540         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update2  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
541
542         DO jl = 1, jpl
543            CALL prt_ctl_info(' ')
544            CALL prt_ctl_info(' - Category : ', ivar1=jl)
545            CALL prt_ctl_info('   ~~~~~~~~~~')
546            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
547            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
548            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
549            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
550            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
551            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
552            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
553            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_a_i     : ')
554            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
555            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
556            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_i     : ')
557            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
558            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
559            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_s     : ')
560            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
561            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
562            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i1    : ')
563            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
564            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
565            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i2    : ')
566            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
567            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
568            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update2  : old_e_snow  : ')
569            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
570            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
571            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update2  : old_smv_i   : ')
572            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
573            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
574            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ')
575            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
576            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update2  : Path. case  : ')
577
578            DO jk = 1, nlay_i
579               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
580               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
581            END DO
582         END DO
583
584         CALL prt_ctl_info(' ')
585         CALL prt_ctl_info(' - Heat / FW fluxes : ')
586         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
587         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
588         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
589         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
590
591         CALL prt_ctl_info(' ')
592         CALL prt_ctl_info(' - Stresses : ')
593         CALL prt_ctl_info('   ~~~~~~~~~~ ')
594         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
595         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
596         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
597      ENDIF
598   
599      CALL wrk_dealloc( jpi,jpj,jpl, internal_melt )   ! integer
600      CALL wrk_dealloc( jkmax, zthick0, zqm0 )
601
602      CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
603
604   END SUBROUTINE lim_update2
605#else
606   !!----------------------------------------------------------------------
607   !!   Default option         Empty Module               No sea-ice model
608   !!----------------------------------------------------------------------
609CONTAINS
610   SUBROUTINE lim_update2     ! Empty routine
611   END SUBROUTINE lim_update2
612
613#endif
614
615END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.