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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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