source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 4293

Last change on this file since 4293 was 4293, checked in by clem, 7 years ago

small corrections for ice categories hi_max and for shifting ice between categories

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