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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 3938

Last change on this file since 3938 was 3938, checked in by flavoni, 11 years ago

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