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

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

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

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

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