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.
limupdate1.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/limupdate1.F90 @ 4072

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

few rewritings
File size: 25.4 KB
RevLine 
[4045]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 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 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   ! Check budget (Rousset)
41   USE in_out_manager   ! I/O manager
42   USE iom              ! I/O manager
43   USE lib_mpp          ! MPP library
[4072]44   USE timing          ! Timing
[4045]45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   lim_update1   ! routine called by ice_step
50
51      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
52      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
53      REAL(wp)  ::   epsi03 = 1.e-03_wp   !    -       -
54      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
55      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
56      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
57      REAL(wp)  ::   rzero  = 0._wp       !    -       -
58      REAL(wp)  ::   rone   = 1._wp       !    -       -
59         
60   !! * Substitutions
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
64   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE lim_update1
70      !!-------------------------------------------------------------------
71      !!               ***  ROUTINE lim_update1  ***
72      !!               
73      !! ** Purpose :  Computes update of sea-ice global variables at
74      !!               the end of the time step.
75      !!               Address pathological cases
76      !!               This place is very important
77      !!               
78      !! ** Method  : 
79      !!    Ice speed from ice dynamics
80      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
81      !!      from advection and ice thermodynamics
82      !!
83      !! ** Action  : -
84      !!---------------------------------------------------------------------
85      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
86      INTEGER ::   jbnd1, jbnd2
87      INTEGER ::   i_ice_switch
88      INTEGER ::   ind_im, layer      ! indices for internal melt
89      REAL(wp) ::   zweight, zesum, z_da_i, zhimax
90      REAL(wp) ::   zinda, zindb, zindsn, zindic
91      REAL(wp) ::   zindg, zh, zdvres, zviold2
92      REAL(wp) ::   zbigvalue, zvsold2, z_da_ex
93      REAL(wp) ::   z_prescr_hi, zat_i_old, ztmelts, ze_s
94
95      REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for
96      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)
97      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
98      ! mass and salt flux (clem)
99      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
100      !!-------------------------------------------------------------------
[4072]101      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
[4045]102
103      CALL wrk_alloc( jkmax, zthick0, zqm0 )
104
105      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
106
107      !------------------------------------------------------------------------------
108      ! 1. Update of Global variables                                               |
109      !------------------------------------------------------------------------------
110
111      !---------------------
112      ! Ice dynamics 
113      !---------------------
114      u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
115      v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
116 
117      !-----------------------------
118      ! Update ice and snow volumes 
119      !-----------------------------
120      DO jl = 1, jpl
121         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl) 
122         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl)
123      END DO
124
125
126      !---------------------------------------------
127      ! Ice concentration and ice heat content
128      !---------------------------------------------
129      a_i (:,:,:)  = a_i (:,:,:) + d_a_i_trp(:,:,:)
130      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:) 
131
132      !------------------------------
133      ! Snow temperature and ice age
134      !------------------------------
135      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:)
136      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:) 
137
138      !--------------
139      ! Ice salinity   
140      !--------------
141      IF(  num_sal == 2  ) THEN      ! general case
142         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:)
143      ENDIF
144
145      ! mass and salt flux init (clem)
146      zviold(:,:,:) = v_i(:,:,:)
147      zvsold(:,:,:) = v_s(:,:,:)
148      zsmvold(:,:,:) = smv_i(:,:,:)
149
150      ! -------------------------------
151      !- check conservation (C Rousset)
152      IF (ln_limdiahsb) THEN
153         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
154         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
155         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
156         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
157      ENDIF
158      !- check conservation (C Rousset)
159      ! -------------------------------
160
161      CALL lim_var_glo2eqv
162
163      !---------------------------------
164      ! Classify the pathological cases
165      !---------------------------------
166      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
167      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
168      ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
169      !
170      DO jl = 1, jpl
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173               patho_case(ji,jj,jl) = 1
174               IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN
175                  patho_case(ji,jj,jl) = 3
176               ENDIF
177               IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
178                  patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell
179               ENDIF
180            END DO
181         END DO
182      END DO
183
184      !--------------------------------------
185      ! 2. Review of all pathological cases
186      !--------------------------------------
187
188      !-------------------------------------------
189      ! 2.1) Advection of ice in an ice-free cell
190      !-------------------------------------------
191      ! should be removed since it is treated after dynamics now
192
193!      !IF( ln_nicep ) THEN 
194!         WRITE(numout,*) ' limupdate1 - before h correction '
195!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
196!         WRITE(numout,*) ' at_i : ', at_i(12,44)
197!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
198!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
199!      !ENDIF
200!
201      zhimax = .3_wp
202      ! first category
203      DO jj = 1, jpj
204         DO ji = 1, jpi
205            !--- the thickness of such an ice is often out of bounds
206            !--- thus we recompute a new area while conserving ice volume
207            zat_i_old = SUM( old_a_i(ji,jj,:) )
208            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
209            IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
210                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
211                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
212               ht_i(ji,jj,1) = hi_max(1) / 2.0
213               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
214            ENDIF
215         END DO
216      END DO
217
218
219!      !IF( ln_nicep ) THEN 
220!         at_i(:,:) = 0._wp
221!         DO jl = 1, jpl
222!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
223!         END DO
224!         WRITE(numout,*) ' limupdate1 - after h correction 1 '
225!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
226!         WRITE(numout,*) ' at_i : ', at_i(12,44)
227!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
228!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
229!      !ENDIF
230
231      zhimax = 1._wp
232      ! other categories
233      DO jl = 2, jpl
234         jm = ice_types(jl)
235         DO jj = 1, jpj
236            DO ji = 1, jpi
237               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
238               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
239               ! it makes problems when the advected volume and concentration do not seem to be
240               ! related with each other
241               ! the new thickness is sometimes very big!
242               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
243               ! which of course is plausible
244               ! but fuck! it fucks everything up :)
245               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
246                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
247                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
248                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
249               ENDIF
250            END DO ! ji
251         END DO !jj
252      END DO !jl
253     
254      at_i(:,:) = 0._wp
255      DO jl = 1, jpl
256         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
257      END DO
258
259!      !IF( ln_nicep ) THEN 
260!         WRITE(numout,*) ' limupdate1 - after h correction 2 '
261!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
262!         WRITE(numout,*) ' at_i : ', at_i(12,44)
263!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
264!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
265!      !ENDIF
266
267      !----------------------------------------------------
268      ! 2.2) Rebin categories with thickness out of bounds
269      !----------------------------------------------------
270      DO jm = 1, jpm
271         jbnd1 = ice_cat_bounds(jm,1)
272         jbnd2 = ice_cat_bounds(jm,2)
273         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
274      END DO
275
276      at_i(:,:) = 0._wp
277      DO jl = 1, jpl
278         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
279      END DO
280
281      zbigvalue      = 1.0e+20
282
283      DO jl = 1, jpl
284         DO jj = 1, jpj 
285            DO ji = 1, jpi
286
287               !switches
288               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
289               !switch = 1 if a_i > 1e-06 and 0 if not
290               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
291               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
292               ! bug fix 25 avril 2007
293               zindb         = zindb*zindic
294
295               !--- 2.3 Correction to ice age
296               !------------------------------
297               !                IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN
298               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday
299               !                ENDIF
300               IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
301                  oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl)
302               ENDIF
303               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
304
305               !--- 2.4 Correction to snow thickness
306               !-------------------------------------
307               !          ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
308               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)
309               ! snow thickness cannot be smaller than 1e-6
310               zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl)
311               v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres
312
313               !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn
314 
315               !--- 2.5 Correction to ice thickness
316               !-------------------------------------
317               zdvres = (zindb - 1._wp) * v_i(ji,jj,jl)
318               v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres
319
320               !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic
321               !sfx_res(ji,jj)  = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice )
322
323               !--- 2.6 Snow is transformed into ice if the original ice cover disappears
324               !----------------------------------------------------------------------------
325               zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
326               zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0
327               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres
328
329               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn )
330               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres
331
332               !--- 2.7 Correction to ice concentrations
333               !--------------------------------------------
334               ! if greater than 0, ice concentration cannot be smaller than 1e-10
335               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
336               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
337
338               !-------------------------
339               ! 2.8) Snow heat content
340               !-------------------------
341               e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0._wp, e_s(ji,jj,1,jl) ), zbigvalue ) )
342
343            END DO ! ji
344         END DO ! jj
345      END DO ! jl
346
347      !------------------------
348      ! 2.9) Ice heat content
349      !------------------------
350
351      DO jl = 1, jpl
352         DO jk = 1, nlay_i
353            DO jj = 1, jpj 
354               DO ji = 1, jpi
355                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
356                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
357               END DO ! ji
358            END DO ! jj
359         END DO !jk
360      END DO !jl
361 
362      at_i(:,:) = 0._wp
363      DO jl = 1, jpl
364         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
365      END DO
366
367      !--- 2.13 ice concentration should not exceed amax
368      !         (it should not be the case)
369      !-----------------------------------------------------
370      DO jj = 1, jpj
371         DO ji = 1, jpi
372            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )       
373            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
374            DO jl  = 1, jpl
375               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
376               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
377               !
378               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
379               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
380               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
381            END DO
382         END DO
383      END DO
384      at_i(:,:) = a_i(:,:,1)
385      DO jl = 2, jpl
386         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
387      END DO
388
389
390      ! Final thickness distribution rebinning
391      ! --------------------------------------
392      DO jm = 1, jpm
393         jbnd1 = ice_cat_bounds(jm,1)
394         jbnd2 = ice_cat_bounds(jm,2)
395         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
396         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
397         ENDIF
398      END DO
399
400
401      !---------------------
402      ! 2.11) Ice salinity
403      !---------------------
404      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
405      IF (  num_sal == 2  ) THEN ! general case
406
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                     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                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
414                     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)
415                  END DO ! ji
416               END DO ! jj
417            !END DO !jk
418         END DO !jl
419
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
[4072]542      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
[4045]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.