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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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