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 @ 3319

Last change on this file since 3319 was 3294, checked in by rblod, 12 years ago

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