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

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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