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

Last change on this file since 834 was 834, checked in by ctlod, 16 years ago

Clean comments and useless lines, see ticket:#72

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