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

source: trunk/NEMO/LIM_SRC_3/limthd_lac.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210

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