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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 3963

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

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116

  • Property svn:executable set to *
File size: 25.3 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 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
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
101      CALL wrk_alloc( jkmax, zthick0, zqm0 )
102
103      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
104
105      !------------------------------------------------------------------------------
106      ! 1. Update of Global variables                                               |
107      !------------------------------------------------------------------------------
108
109      !---------------------
110      ! Ice dynamics 
111      !---------------------
112      u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
113      v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
114 
115      !-----------------------------
116      ! Update ice and snow volumes 
117      !-----------------------------
118      DO jl = 1, jpl
119         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl) 
120         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl)
121      END DO
122
123
124      !---------------------------------------------
125      ! Ice concentration and ice heat content
126      !---------------------------------------------
127      a_i (:,:,:)  = a_i (:,:,:) + d_a_i_trp(:,:,:)
128      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:) 
129
130      !------------------------------
131      ! Snow temperature and ice age
132      !------------------------------
133      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:)
134      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:) 
135
136      !--------------
137      ! Ice salinity   
138      !--------------
139      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case
140         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:)
141      ENDIF
142
143      ! mass and salt flux init (clem)
144      zviold(:,:,:) = v_i(:,:,:)
145      zvsold(:,:,:) = v_s(:,:,:)
146      zsmvold(:,:,:) = smv_i(:,:,:)
147
148      ! -------------------------------
149      !- check conservation (C Rousset)
150      IF (ln_limdiahsb) THEN
151         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
152         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
153         zchk_fw_b  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) )
154         zchk_fs_b  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) )
155      ENDIF
156      !- check conservation (C Rousset)
157      ! -------------------------------
158
159      CALL lim_var_glo2eqv
160
161      !---------------------------------
162      ! Classify the pathological cases
163      !---------------------------------
164      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
165      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
166      ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
167      !
168      DO jl = 1, jpl
169         DO jj = 1, jpj
170            DO ji = 1, jpi
171               patho_case(ji,jj,jl) = 1
172               IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN
173                  patho_case(ji,jj,jl) = 3
174               ENDIF
175               IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
176                  patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell
177               ENDIF
178            END DO
179         END DO
180      END DO
181
182      !--------------------------------------
183      ! 2. Review of all pathological cases
184      !--------------------------------------
185
186      !-------------------------------------------
187      ! 2.1) Advection of ice in an ice-free cell
188      !-------------------------------------------
189      ! should be removed since it is treated after dynamics now
190
191!      !IF( ln_nicep ) THEN 
192!         WRITE(numout,*) ' limupdate1 - before h correction '
193!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
194!         WRITE(numout,*) ' at_i : ', at_i(12,44)
195!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
196!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
197!      !ENDIF
198!
199      zhimax = 1._wp
200      ! first category
201      DO jj = 1, jpj
202         DO ji = 1, jpi
203            !--- the thickness of such an ice is often out of bounds
204            !--- thus we recompute a new area while conserving ice volume
205            zat_i_old = SUM( old_a_i(ji,jj,:) )
206            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
207            IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
208                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
209                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
210               ht_i(ji,jj,1) = hi_max(1) / 2.0
211               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
212            ENDIF
213         END DO
214      END DO
215
216
217!      !IF( ln_nicep ) THEN 
218!         at_i(:,:) = 0._wp
219!         DO jl = 1, jpl
220!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
221!         END DO
222!         WRITE(numout,*) ' limupdate1 - after h correction 1 '
223!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
224!         WRITE(numout,*) ' at_i : ', at_i(12,44)
225!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
226!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
227!      !ENDIF
228
229      zhimax = 10._wp
230      ! other categories
231      DO jl = 2, jpl
232         jm = ice_types(jl)
233         DO jj = 1, jpj
234            DO ji = 1, jpi
235               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
236               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
237               ! it makes problems when the advected volume and concentration do not seem to be
238               ! related with each other
239               ! the new thickness is sometimes very big!
240               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
241               ! which of course is plausible
242               ! but fuck! it fucks everything up :)
243               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
244                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
245                  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
246                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
247               ENDIF
248            END DO ! ji
249         END DO !jj
250      END DO !jl
251     
252      at_i(:,:) = 0._wp
253      DO jl = 1, jpl
254         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
255      END DO
256
257!      !IF( ln_nicep ) THEN 
258!         WRITE(numout,*) ' limupdate1 - after h correction 2 '
259!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
260!         WRITE(numout,*) ' at_i : ', at_i(12,44)
261!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
262!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
263!      !ENDIF
264
265      !----------------------------------------------------
266      ! 2.2) Rebin categories with thickness out of bounds
267      !----------------------------------------------------
268      DO jm = 1, jpm
269         jbnd1 = ice_cat_bounds(jm,1)
270         jbnd2 = ice_cat_bounds(jm,2)
271         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
272      END DO
273
274      at_i(:,:) = 0._wp
275      DO jl = 1, jpl
276         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
277      END DO
278
279      zbigvalue      = 1.0e+20
280
281      DO jl = 1, jpl
282         DO jj = 1, jpj 
283            DO ji = 1, jpi
284
285               !switches
286               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
287               !switch = 1 if a_i > 1e-06 and 0 if not
288               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
289               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
290               ! bug fix 25 avril 2007
291               zindb         = zindb*zindic
292
293               !--- 2.3 Correction to ice age
294               !------------------------------
295               !                IF ((o_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*float(numit))) THEN
296               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
297               !                ENDIF
298               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
299                  oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl)
300               ENDIF
301               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
302
303               !--- 2.4 Correction to snow thickness
304               !-------------------------------------
305               !          ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
306               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)
307               ! snow thickness cannot be smaller than 1e-6
308               zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl)
309               v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres
310
311               !rdmsnif(ji,jj) = rdmsnif(ji,jj) + zdvres * rhosn
312 
313               !--- 2.5 Correction to ice thickness
314               !-------------------------------------
315               zdvres = (zindb - 1._wp) * v_i(ji,jj,jl)
316               v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres
317
318               !rdmicif(ji,jj) = rdmicif(ji,jj) + zdvres * rhoic
319               !fsalt_res(ji,jj)  = fsalt_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice )
320
321               !--- 2.6 Snow is transformed into ice if the original ice cover disappears
322               !----------------------------------------------------------------------------
323               zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
324               zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0
325               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres
326
327               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn )
328               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres
329
330               !--- 2.7 Correction to ice concentrations
331               !--------------------------------------------
332               ! if greater than 0, ice concentration cannot be smaller than 1e-10
333               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
334               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
335
336               !-------------------------
337               ! 2.8) Snow heat content
338               !-------------------------
339               e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0._wp, e_s(ji,jj,1,jl) ), zbigvalue ) )
340
341            END DO ! ji
342         END DO ! jj
343      END DO ! jl
344
345      !------------------------
346      ! 2.9) Ice heat content
347      !------------------------
348
349      DO jl = 1, jpl
350         DO jk = 1, nlay_i
351            DO jj = 1, jpj 
352               DO ji = 1, jpi
353                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
354                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
355               END DO ! ji
356            END DO ! jj
357         END DO !jk
358      END DO !jl
359 
360      at_i(:,:) = 0._wp
361      DO jl = 1, jpl
362         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
363      END DO
364
365      !--- 2.13 ice concentration should not exceed amax
366      !         (it should not be the case)
367      !-----------------------------------------------------
368      DO jj = 1, jpj
369         DO ji = 1, jpi
370            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )       
371            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
372            DO jl  = 1, jpl
373               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
374               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
375               !
376               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
377               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
378               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
379            END DO
380         END DO
381      END DO
382      at_i(:,:) = a_i(:,:,1)
383      DO jl = 2, jpl
384         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
385      END DO
386
387
388      ! Final thickness distribution rebinning
389      ! --------------------------------------
390      DO jm = 1, jpm
391         jbnd1 = ice_cat_bounds(jm,1)
392         jbnd2 = ice_cat_bounds(jm,2)
393         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
394         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
395         ENDIF
396      END DO
397
398
399      !---------------------
400      ! 2.11) Ice salinity
401      !---------------------
402      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
403      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
404
405         DO jl = 1, jpl
406            !DO jk = 1, nlay_i
407               DO jj = 1, jpj 
408                  DO ji = 1, jpi
409                     ! salinity stays in bounds
410                     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) )
411                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
412                     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)
413                  END DO ! ji
414               END DO ! jj
415            !END DO !jk
416         END DO !jl
417
418      ENDIF
419
420      at_i(:,:) = a_i(:,:,1)
421      DO jl = 2, jpl
422         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
423      END DO
424
425
426      !--------------------------------
427      ! Update mass/salt fluxes (clem)
428      !--------------------------------
429      DO jl = 1, jpl
430         DO jj = 1, jpj 
431            DO ji = 1, jpi
432               diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 
433               rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
434               rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
435               fsalt_res(ji,jj) = fsalt_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 
436            END DO
437         END DO
438      END DO
439 
440      ! -------------------------------
441      !- check conservation (C Rousset)
442      IF (ln_limdiahsb) THEN
443
444         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
445         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
446 
447         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
448         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
449
450         zchk_vmin = glob_min(v_i)
451         zchk_amax = glob_max(SUM(a_i,dim=3))
452         zchk_amin = glob_min(a_i)
453
454         IF(lwp) THEN
455            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * 86400.)
456            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.)
457            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3)
458            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax
459            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin
460         ENDIF
461      ENDIF
462      !- check conservation (C Rousset)
463      ! -------------------------------
464
465      IF(ln_ctl) THEN   ! Control print
466         CALL prt_ctl_info(' ')
467         CALL prt_ctl_info(' - Cell values : ')
468         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
469         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :')
470         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
471         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
472         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
473         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
474         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
475         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 :')
476         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
477
478         DO jl = 1, jpl
479            CALL prt_ctl_info(' ')
480            CALL prt_ctl_info(' - Category : ', ivar1=jl)
481            CALL prt_ctl_info('   ~~~~~~~~~~')
482            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
483            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
484            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
485            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
486            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
487            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
488            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
489            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ')
490            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
491            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
492            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ')
493            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
494            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
495            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ')
496            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
497            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
498            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ')
499            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
500            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
501            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ')
502            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
503            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
504            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ')
505            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
506            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
507            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ')
508            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
509            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
510            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ')
511            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
512            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update1  : Path. case  : ')
513
514            DO jk = 1, nlay_i
515               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
516               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
517            END DO
518         END DO
519
520         CALL prt_ctl_info(' ')
521         CALL prt_ctl_info(' - Heat / FW fluxes : ')
522         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
523         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
524         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
525         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
526
527         CALL prt_ctl_info(' ')
528         CALL prt_ctl_info(' - Stresses : ')
529         CALL prt_ctl_info('   ~~~~~~~~~~ ')
530         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
531         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
532         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
533      ENDIF
534   
535
536      CALL wrk_dealloc( jkmax, zthick0, zqm0 )
537
538      CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
539
540   END SUBROUTINE lim_update1
541#else
542   !!----------------------------------------------------------------------
543   !!   Default option         Empty Module               No sea-ice model
544   !!----------------------------------------------------------------------
545CONTAINS
546   SUBROUTINE lim_update1     ! Empty routine
547   END SUBROUTINE lim_update1
548
549#endif
550
551END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.