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.
limthd_lac.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 3558

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

Fix issues when using key_nosignedzeo, see ticket #996
  • Property svn:keywords set to Id
File size: 33.0 KB
Line 
1MODULE limthd_lac
2   !!======================================================================
3   !!                       ***  MODULE limthd_lac   ***
4   !!                lateral thermodynamic growth of the ice
5   !!======================================================================
6   !! History :  LIM  ! 2005-12 (M. Vancoppenolle)  Original code
7   !!             -   ! 2006-01 (M. Vancoppenolle)  add ITD
8   !!            3.0  ! 2007-07 (M. Vancoppenolle)  Mass and energy conservation tested
9   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_lat_acr    : lateral accretion of ice
16   !!----------------------------------------------------------------------
17   USE par_oce          ! ocean parameters
18   USE dom_oce          ! domain variables
19   USE phycst           ! physical constants
20   USE sbc_oce          ! Surface boundary condition: ocean fields
21   USE sbc_ice          ! Surface boundary condition: ice fields
22   USE thd_ice          ! LIM thermodynamics
23   USE dom_ice          ! LIM domain
24   USE par_ice          ! LIM parameters
25   USE ice              ! LIM variables
26   USE limtab           ! LIM 2D <==> 1D
27   USE limcons          ! LIM conservation
28   USE in_out_manager   ! I/O manager
29   USE lib_mpp          ! MPP library
30   USE wrk_nemo         ! work arrays
31   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC lim_thd_lac     ! called by lim_thd
37
38   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
39   REAL(wp) ::   epsi13 = 1e-13_wp   !
40   REAL(wp) ::   epsi11 = 1e-11_wp   !
41   REAL(wp) ::   epsi10 = 1e-10_wp   !
42   REAL(wp) ::   epsi06 = 1e-06_wp   !
43   REAL(wp) ::   epsi03 = 1e-03_wp   !
44   REAL(wp) ::   zzero  = 0._wp      !
45   REAL(wp) ::   zone   = 1._wp      !
46
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE lim_thd_lac
55      !!-------------------------------------------------------------------
56      !!               ***   ROUTINE lim_thd_lac  ***
57      !! 
58      !! ** Purpose : Computation of the evolution of the ice thickness and
59      !!      concentration as a function of the heat balance in the leads.
60      !!      It is only used for lateral accretion
61      !!       
62      !! ** Method  : Ice is formed in the open water when ocean lose heat
63      !!      (heat budget of open water Bl is negative) .
64      !!      Computation of the increase of 1-A (ice concentration) fol-
65      !!      lowing the law :
66      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
67      !!       where - h0 is the thickness of ice created in the lead
68      !!             - a is a minimum fraction for leads
69      !!             - F is a monotonic non-increasing function defined as:
70      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
71      !!             - exld is the exponent closure rate (=2 default val.)
72      !!
73      !! ** Action : - Adjustment of snow and ice thicknesses and heat
74      !!                content in brine pockets
75      !!             - Updating ice internal temperature
76      !!             - Computation of variation of ice volume and mass
77      !!             - Computation of frldb after lateral accretion and
78      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)     
79      !!------------------------------------------------------------------------
80      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
81      INTEGER ::   layer, nbpac     ! local integers
82      INTEGER ::   zji, zjj, iter   !   -       -
83      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde  ! local scalars
84      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
85      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
86      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
87      CHARACTER (len = 15) :: fieldid
88      !
89      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
90      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
91
92      REAL(wp), POINTER, DIMENSION(:) ::   zv_newice   ! volume of accreted ice
93      REAL(wp), POINTER, DIMENSION(:) ::   za_newice   ! fractional area of accreted ice
94      REAL(wp), POINTER, DIMENSION(:) ::   zh_newice   ! thickness of accreted ice
95      REAL(wp), POINTER, DIMENSION(:) ::   ze_newice   ! heat content of accreted ice
96      REAL(wp), POINTER, DIMENSION(:) ::   zs_newice   ! salinity of accreted ice
97      REAL(wp), POINTER, DIMENSION(:) ::   zo_newice   ! age of accreted ice
98      REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget
99      REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
100      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction   
101      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
102      REAL(wp), POINTER, DIMENSION(:) ::   zdh_frazb   ! accretion of frazil ice at the ice bottom
103      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector)
104
105      REAL(wp), POINTER, DIMENSION(:,:) ::   zhice_old   ! previous ice thickness
106      REAL(wp), POINTER, DIMENSION(:,:) ::   zdummy      ! dummy thickness of new ice
107      REAL(wp), POINTER, DIMENSION(:,:) ::   zdhicbot    ! thickness of new ice which is accreted vertically
108      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
109      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
110      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i
111      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i
112      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i
113      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i
114
115      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i
116
117      REAL(wp), POINTER, DIMENSION(:) ::   zqbgow    ! heat budget of the open water (negative)
118      REAL(wp), POINTER, DIMENSION(:) ::   zdhex     ! excessively thick accreted sea ice (hlead-hice)
119
120      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content
121      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness
122
123      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
124      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories
125      REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories
126      REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories
127      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
128      !!-----------------------------------------------------------------------!
129
130      CALL wrk_alloc( jpij, zcatac )   ! integer
131      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
132      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
133      CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
134      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
135      CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
136      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
137
138      et_i_init(:,:) = 0._wp
139      et_s_init(:,:) = 0._wp
140      vt_i_init(:,:) = 0._wp
141      vt_s_init(:,:) = 0._wp
142
143      !------------------------------------------------------------------------------!
144      ! 1) Conservation check and changes in each ice category
145      !------------------------------------------------------------------------------!
146      IF ( con_i ) THEN
147         CALL lim_column_sum (jpl, v_i, vt_i_init)
148         CALL lim_column_sum (jpl, v_s, vt_s_init)
149         CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
150         CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
151      ENDIF
152
153      !------------------------------------------------------------------------------|
154      ! 2) Convert units for ice internal energy
155      !------------------------------------------------------------------------------|
156      DO jl = 1, jpl
157         DO jk = 1, nlay_i
158            DO jj = 1, jpj
159               DO ji = 1, jpi
160                  !Energy of melting q(S,T) [J.m-3]
161                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / &
162                     MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * &
163                     nlay_i
164                  zindb      = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes
165                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb
166               END DO
167            END DO
168         END DO
169      END DO
170
171      !------------------------------------------------------------------------------!
172      ! 3) Collection thickness of ice formed in leads and polynyas
173      !------------------------------------------------------------------------------!   
174      ! hicol is the thickness of new ice formed in open water
175      ! hicol can be either prescribed (frazswi = 0)
176      ! or computed (frazswi = 1)
177      ! Frazil ice forms in open water, is transported by wind
178      ! accumulates at the edge of the consolidated ice edge
179      ! where it forms aggregates of a specific thickness called
180      ! collection thickness.
181
182      ! Note : the following algorithm currently breaks vectorization
183      !
184
185      zvrel(:,:) = 0.0
186
187      ! Default new ice thickness
188      DO jj = 1, jpj
189         DO ji = 1, jpi
190            hicol(ji,jj) = hiccrit(1)
191         END DO
192      END DO
193
194      IF (fraz_swi.eq.1.0) THEN
195
196         !--------------------
197         ! Physical constants
198         !--------------------
199         hicol(:,:) = 0.0
200
201         zhicrit = 0.04 ! frazil ice thickness
202         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
203         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
204         zgamafr = 0.03
205
206         DO jj = 1, jpj
207            DO ji = 1, jpi
208
209               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
210                  !-------------
211                  ! Wind stress
212                  !-------------
213                  ! C-grid wind stress components
214                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  ) &
215                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) / 2.0
216                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1) &
217                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) / 2.0
218                  ! Square root of wind stress
219                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
220
221                  !---------------------
222                  ! Frazil ice velocity
223                  !---------------------
224                  zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
225                  zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
226
227                  !-------------------
228                  ! Pack ice velocity
229                  !-------------------
230                  ! C-grid ice velocity
231                  zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) ))
232                  zvgx  = zindb * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  ) &
233                     + u_ice(ji,jj    ) * tmu(ji  ,jj  ) ) / 2.0
234                  zvgy  = zindb * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1) &
235                     + v_ice(ji,jj    ) * tmv(ji  ,jj  ) ) / 2.0
236
237                  !-----------------------------------
238                  ! Relative frazil/pack ice velocity
239                  !-----------------------------------
240                  ! absolute relative velocity
241                  zvrel2        = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + &
242                     ( zvfry - zvgy ) * ( zvfry - zvgy )   &
243                     , 0.15 * 0.15 )
244                  zvrel(ji,jj)  = SQRT(zvrel2)
245
246                  !---------------------
247                  ! Iterative procedure
248                  !---------------------
249                  hicol(ji,jj) = zhicrit + 0.1 
250                  hicol(ji,jj) = zhicrit + hicol(ji,jj) /      & 
251                     ( hicol(ji,jj) * hicol(ji,jj) - &
252                     zhicrit * zhicrit ) * ztwogp * zvrel2
253
254                  iter = 1
255                  iterate_frazil = .true.
256
257                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
258                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
259                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
260                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
261                        - zhicrit * ztwogp * zvrel2
262                     zhicol_new = hicol(ji,jj) - zf/zfp
263                     hicol(ji,jj)   = zhicol_new
264
265                     iter = iter + 1
266
267                  END DO ! do while
268
269               ENDIF ! end of selection of pixels where ice forms
270
271            END DO ! loop on ji ends
272         END DO ! loop on jj ends
273
274      ENDIF ! End of computation of frazil ice collection thickness
275
276      !------------------------------------------------------------------------------!
277      ! 4) Identify grid points where new ice forms
278      !------------------------------------------------------------------------------!
279
280      !-------------------------------------
281      ! Select points for new ice formation
282      !-------------------------------------
283      ! This occurs if open water energy budget is negative
284      nbpac = 0
285      DO jj = 1, jpj
286         DO ji = 1, jpi
287            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
288               nbpac = nbpac + 1
289               npac( nbpac ) = (jj - 1) * jpi + ji
290               IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN
291                  jiindex_1d = nbpac
292               ENDIF
293            ENDIF
294         END DO
295      END DO
296
297      IF( ln_nicep ) THEN
298         WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
299      ENDIF
300
301      !------------------------------
302      ! Move from 2-D to 1-D vectors
303      !------------------------------
304      ! If ocean gains heat do nothing
305      ! 0therwise compute new ice formation
306
307      IF ( nbpac > 0 ) THEN
308
309         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         ,       &
310            jpi, jpj, npac(1:nbpac) )
311         DO jl = 1, jpl
312            CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  ,       &
313               jpi, jpj, npac(1:nbpac) )
314            CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  ,       &
315               jpi, jpj, npac(1:nbpac) )
316            CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) ,       &
317               jpi, jpj, npac(1:nbpac) )
318            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl),       &
319               jpi, jpj, npac(1:nbpac) )
320            DO jk = 1, nlay_i
321               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , &
322                  jpi, jpj, npac(1:nbpac) )
323            END DO ! jk
324         END DO ! jl
325
326         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif ,              &
327            jpi, jpj, npac(1:nbpac) )
328         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif ,              &
329            jpi, jpj, npac(1:nbpac) )
330         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  ,              &
331            jpi, jpj, npac(1:nbpac) )
332         CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv ,              &
333            jpi, jpj, npac(1:nbpac) )
334         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol ,              &
335            jpi, jpj, npac(1:nbpac) )
336         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel ,              &
337            jpi, jpj, npac(1:nbpac) )
338
339         !------------------------------------------------------------------------------!
340         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
341         !------------------------------------------------------------------------------!
342
343         !----------------------
344         ! Thickness of new ice
345         !----------------------
346         DO ji = 1, nbpac
347            zh_newice(ji)     = hiccrit(1)
348         END DO
349         IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:)
350
351         !----------------------
352         ! Salinity of new ice
353         !----------------------
354
355         IF ( num_sal .EQ. 1 ) THEN
356            zs_newice(:)      =   bulk_sal
357         ENDIF ! num_sal
358
359         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
360
361            DO ji = 1, nbpac
362               zs_newice(ji)  =   MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max )
363               zji            =   MOD( npac(ji) - 1, jpi ) + 1
364               zjj            =   ( npac(ji) - 1 ) / jpi + 1
365               zs_newice(ji)  =   MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) )
366            END DO ! jl
367
368         ENDIF ! num_sal
369
370         IF ( num_sal .EQ. 3 ) THEN
371            zs_newice(:)      =   2.3
372         ENDIF ! num_sal
373
374         !-------------------------
375         ! Heat content of new ice
376         !-------------------------
377         ! We assume that new ice is formed at the seawater freezing point
378         DO ji = 1, nbpac
379            ztmelts           = - tmut * zs_newice(ji) + rtt ! Melting point (K)
380            ze_newice(ji)     =   rhoic * ( cpic * ( ztmelts - t_bo_b(ji) )    &
381               + lfus * ( 1.0 - ( ztmelts - rtt )   &
382               / ( t_bo_b(ji) - rtt ) )           &
383               - rcp * ( ztmelts-rtt ) )
384            ze_newice(ji)     =   MAX( ze_newice(ji) , 0.0 ) +                 &
385               MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) )   & 
386               * rhoic * lfus
387         END DO ! ji
388         !----------------
389         ! Age of new ice
390         !----------------
391         DO ji = 1, nbpac
392            zo_newice(ji)     = 0.0
393         END DO ! ji
394
395         !--------------------------
396         ! Open water energy budget
397         !--------------------------
398         DO ji = 1, nbpac
399            zqbgow(ji)        = qldif_1d(ji) - qcmif_1d(ji) !<0
400         END DO ! ji
401
402         !-------------------
403         ! Volume of new ice
404         !-------------------
405         DO ji = 1, nbpac
406            zv_newice(ji)     = - zqbgow(ji) / ze_newice(ji)
407
408            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
409            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) )     & 
410               + 1.0 ) / 2.0 * maxfrazb
411            zdh_frazb(ji) = zfrazb*zv_newice(ji)
412            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
413         END DO
414
415         !---------------------------------
416         ! Salt flux due to new ice growth
417         !---------------------------------
418         IF ( ( num_sal .EQ. 4 ) ) THEN
419            DO ji = 1, nbpac
420               zji            = MOD( npac(ji) - 1, jpi ) + 1
421               zjj            = ( npac(ji) - 1 ) / jpi + 1
422               fseqv_1d(ji)   = fseqv_1d(ji) +                                     &
423                  ( sss_m(zji,zjj) - bulk_sal      ) * rhoic *       &
424                  zv_newice(ji) / rdt_ice
425            END DO
426         ELSE
427            DO ji = 1, nbpac
428               zji            = MOD( npac(ji) - 1, jpi ) + 1
429               zjj            = ( npac(ji) - 1 ) / jpi + 1
430               fseqv_1d(ji)   = fseqv_1d(ji) +                                     &
431                  ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic *       &
432                  zv_newice(ji) / rdt_ice
433            END DO ! ji
434         ENDIF
435
436         !------------------------------------
437         ! Diags for energy conservation test
438         !------------------------------------
439         DO ji = 1, nbpac
440            ! Volume
441            zji                  = MOD( npac(ji) - 1, jpi ) + 1
442            zjj                  = ( npac(ji) - 1 ) / jpi + 1
443            vt_i_init(zji,zjj)   = vt_i_init(zji,zjj) + zv_newice(ji)
444            ! Energy
445            zde                  = ze_newice(ji) / unit_fac
446            zde                  = zde * area(zji,zjj) * zv_newice(ji)
447            et_i_init(zji,zjj)   = et_i_init(zji,zjj) + zde
448         END DO
449
450         ! keep new ice volume in memory
451         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , &
452            jpi, jpj )
453
454         !-----------------
455         ! Area of new ice
456         !-----------------
457         DO ji = 1, nbpac
458            za_newice(ji)     = zv_newice(ji) / zh_newice(ji)
459            ! diagnostic
460            zji                  = MOD( npac(ji) - 1, jpi ) + 1
461            zjj                  = ( npac(ji) - 1 ) / jpi + 1
462            diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice
463         END DO !ji
464
465         !------------------------------------------------------------------------------!
466         ! 6) Redistribute new ice area and volume into ice categories                  !
467         !------------------------------------------------------------------------------!
468
469         !-----------------------------------------
470         ! Keep old ice areas and volume in memory
471         !-----------------------------------------
472         zv_old(:,:) = zv_i_ac(:,:) 
473         za_old(:,:) = za_i_ac(:,:)
474
475         !-------------------------------------------
476         ! Compute excessive new ice area and volume
477         !-------------------------------------------
478         ! If lateral ice growth gives an ice concentration gt 1, then
479         ! we keep the excessive volume in memory and attribute it later
480         ! to bottom accretion
481         DO ji = 1, nbpac
482            ! vectorize
483            IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN
484               zda_res(ji)    = za_newice(ji) - (1.0 - zat_i_ac(ji) )
485               zdv_res(ji)    = zda_res(ji) * zh_newice(ji) 
486               za_newice(ji)  = za_newice(ji) - zda_res(ji)
487               zv_newice(ji)  = zv_newice(ji) - zdv_res(ji)
488            ELSE
489               zda_res(ji) = 0.0
490               zdv_res(ji) = 0.0
491            ENDIF
492         END DO ! ji
493
494         !------------------------------------------------
495         ! Laterally redistribute new ice volume and area
496         !------------------------------------------------
497         zat_i_ac(:) = 0._wp
498         DO jl = 1, jpl
499            DO ji = 1, nbpac
500               IF(  hi_max   (jl-1)  <  zh_newice(ji)   .AND.   &
501                  & zh_newice(ji)    <= hi_max   (jl)         ) THEN
502                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
503                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
504                  zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
505                  zcatac  (ji)    = jl
506               ENDIF
507            END DO ! ji
508         END DO ! jl
509
510         !----------------------------------
511         ! Heat content - lateral accretion
512         !----------------------------------
513         DO ji = 1, nbpac
514            jl = zcatac(ji)                                                           ! categroy in which new ice is put
515            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) )             ! zindb=1 if ice =0 otherwise
516            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
517            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
518            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi11 ) )   ! ice totally new in jl category
519         END DO
520
521         DO jk = 1, nlay_i
522            DO ji = 1, nbpac
523               jl = zcatac(ji)
524               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
525               zalphai = MIN( zhice_old(ji,jl) *   jk       / nlay_i , zh_newice(ji) )   &
526                  &    - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) )
527               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
528                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / nlay_i   &
529                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
530                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i )
531            END DO
532         END DO
533
534         !-----------------------------------------------
535         ! Add excessive volume of new ice at the bottom
536         !-----------------------------------------------
537         ! If the ice concentration exceeds 1, the remaining volume of new ice
538         ! is equally redistributed among all ice categories in which there is
539         ! ice
540
541         ! Fraction of level ice
542         jm = 1
543         zat_i_lev(:) = 0._wp
544
545         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
546            DO ji = 1, nbpac
547               zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 
548            END DO
549         END DO
550
551         IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
552         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
553            DO ji = 1, nbpac
554               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
555               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
556            END DO
557         END DO
558         IF( ln_nicep )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
559
560         !---------------------------------
561         ! Heat content - bottom accretion
562         !---------------------------------
563         jm = 1
564         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
565            DO ji = 1, nbpac
566               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) )       ! zindb=1 if ice =0 otherwise
567               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
568               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
569                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce
570               zdummy(ji,jl)    = zv_i_ac(ji,jl)/MAX(za_i_ac(ji,jl),epsi10)*zindb      ! thickness of residual ice
571            END DO
572         END DO
573
574         ! old layers thicknesses and enthalpies
575         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
576            DO jk = 1, nlay_i
577               DO ji = 1, nbpac
578                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / nlay_i
579                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
580               END DO
581            END DO
582         END DO
583!!gm ???  why the previous do loop  if ocerwriten by the following one ?
584         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
585            DO ji = 1, nbpac
586               zthick0(ji,nlay_i+1,jl) =  zdhicbot(ji,jl)
587               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zdhicbot(ji,jl)
588            END DO ! ji
589         END DO ! jl
590
591         ! Redistributing energy on the new grid
592         ze_i_ac(:,:,:) = 0._wp
593         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
594            DO jk = 1, nlay_i
595               DO layer = 1, nlay_i + 1
596                  DO ji = 1, nbpac
597                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) 
598                     ! Redistributing energy on the new grid
599                     zweight = MAX (  MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk )   &
600                        &    - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) * ( jk - 1 ) ) , 0._wp )   &
601                        &    /( MAX(nlay_i * zthick0(ji,layer,jl),epsi10) ) * zindb
602                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 
603                  END DO ! ji
604               END DO ! layer
605            END DO ! jk
606         END DO ! jl
607
608         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
609            DO jk = 1, nlay_i
610               DO ji = 1, nbpac
611                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) 
612                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
613                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i * zindb
614               END DO
615            END DO
616         END DO
617
618         !------------
619         ! Update age
620         !------------
621         DO jl = 1, jpl
622            DO ji = 1, nbpac
623               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes
624               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
625            END DO
626         END DO   
627
628         !-----------------
629         ! Update salinity
630         !-----------------
631         IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
632            DO jl = 1, jpl
633               DO ji = 1, nbpac
634                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes
635                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
636                  zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb
637               END DO
638            END DO   
639         ENDIF
640
641         !------------------------------------------------------------------------------!
642         ! 8) Change 2D vectors to 1D vectors
643         !------------------------------------------------------------------------------!
644         DO jl = 1, jpl
645            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
646            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
647            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
648            IF (  num_sal == 2  .OR.  num_sal == 4  )   &
649               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
650            DO jk = 1, nlay_i
651               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
652            END DO
653         END DO
654         CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d  (1:nbpac) , jpi, jpj )
655         !
656      ENDIF ! nbpac > 0
657
658      !------------------------------------------------------------------------------!
659      ! 9) Change units for e_i
660      !------------------------------------------------------------------------------!   
661      DO jl = 1, jpl
662         DO jk = 1, nlay_i          ! heat content in 10^9 Joules
663            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i  / unit_fac 
664         END DO
665      END DO
666
667      !------------------------------------------------------------------------------|
668      ! 10) Conservation check and changes in each ice category
669      !------------------------------------------------------------------------------|
670      IF( con_i ) THEN
671         CALL lim_column_sum (jpl,   v_i, vt_i_final)
672         fieldid = 'v_i, limthd_lac'
673         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
674         !
675         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
676         fieldid = 'e_i, limthd_lac'
677         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
678         !
679         CALL lim_column_sum (jpl,   v_s, vt_s_final)
680         fieldid = 'v_s, limthd_lac'
681         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
682         !
683         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
684         !     fieldid = 'e_s, limthd_lac'
685         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
686         IF( ln_nicep ) THEN
687            WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx)
688            WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx)
689            WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx)
690            WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx)
691         ENDIF
692         !
693      ENDIF
694      !
695      CALL wrk_dealloc( jpij, zcatac )   ! integer
696      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
697      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
698      CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
699      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
700      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
701      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
702      !
703   END SUBROUTINE lim_thd_lac
704
705#else
706   !!----------------------------------------------------------------------
707   !!   Default option                               NO  LIM3 sea-ice model
708   !!----------------------------------------------------------------------
709CONTAINS
710   SUBROUTINE lim_thd_lac           ! Empty routine
711   END SUBROUTINE lim_thd_lac
712#endif
713
714   !!======================================================================
715END MODULE limthd_lac
Note: See TracBrowser for help on using the repository browser.