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.
icethd_do.F90 in NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icethd_do.F90 @ 11243

Last change on this file since 11243 was 11229, checked in by clem, 5 years ago

debug BDY-ice by allowing any number of categories as input. Ticket #2300 may also be solved

  • Property svn:keywords set to Id
File size: 23.6 KB
Line 
1MODULE icethd_do
2   !!======================================================================
3   !!                       ***  MODULE icethd_do   ***
4   !!   sea-ice: sea ice growth in the leads (open water) 
5   !!======================================================================
6   !! History :       !  2005-12  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_thd_do        : ice growth in open water (=lateral accretion of ice)
14   !!   ice_thd_do_init   : initialization
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE sbc_oce , ONLY : sss_m
19   USE sbc_ice , ONLY : utau_ice, vtau_ice
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE ice            ! sea-ice: variables
22   USE icetab         ! sea-ice: 2D <==> 1D
23   USE icectl         ! sea-ice: conservation
24   USE icethd_ent     ! sea-ice: thermodynamics, enthalpy
25   USE icevar         ! sea-ice: operations
26   USE icethd_sal     ! sea-ice: salinity profiles
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
31   USE lbclnk         ! lateral boundary conditions (or mpp links)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_thd_do        ! called by ice_thd
37   PUBLIC   ice_thd_do_init   ! called by ice_stp
38
39   !                          !!** namelist (namthd_do) **
40   REAL(wp) ::   rn_hinew      ! thickness for new ice formation (m)
41   LOGICAL  ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F)
42   REAL(wp) ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
43   REAL(wp) ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
44   REAL(wp) ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
45
46   !!----------------------------------------------------------------------
47   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE ice_thd_do
54      !!-------------------------------------------------------------------
55      !!               ***   ROUTINE ice_thd_do  ***
56      !! 
57      !! ** Purpose : Computation of the evolution of the ice thickness and
58      !!              concentration as a function of the heat balance in the leads
59      !!       
60      !! ** Method  : Ice is formed in the open water when ocean looses heat
61      !!              (heat budget of open water is negative) following
62      !!
63      !!       (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
64      !!          where - h0 is the thickness of ice created in the lead
65      !!                - a is a minimum fraction for leads
66      !!                - F is a monotonic non-increasing function defined as:
67      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
68      !!                - exld is the exponent closure rate (=2 default val.)
69      !!
70      !! ** Action : - Adjustment of snow and ice thicknesses and heat
71      !!                content in brine pockets
72      !!             - Updating ice internal temperature
73      !!             - Computation of variation of ice volume and mass
74      !!             - Computation of a_i after lateral accretion and
75      !!               update h_s_1d, h_i_1d     
76      !!------------------------------------------------------------------------
77      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
78      INTEGER  ::   iter             !   -       -
79      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                               ! local scalars
80      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
81      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
82      !
83      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
84      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
85      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
86      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
87      !
88      REAL(wp) ::   zv_newfra
89      !
90      INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
91      REAL(wp), DIMENSION(jpij) ::   zswinew     ! switch for new ice or not
92      !
93      REAL(wp), DIMENSION(jpij) ::   zv_newice   ! volume of accreted ice
94      REAL(wp), DIMENSION(jpij) ::   za_newice   ! fractional area of accreted ice
95      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! thickness of accreted ice
96      REAL(wp), DIMENSION(jpij) ::   ze_newice   ! heat content of accreted ice
97      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of accreted ice
98      REAL(wp), DIMENSION(jpij) ::   zo_newice   ! age of accreted ice
99      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
100      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
101      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
102      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
103      !
104      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b    ! old volume of ice in category jl
105      REAL(wp), DIMENSION(jpij,jpl) ::   za_b    ! old area of ice in category jl
106      !
107      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
108      !
109      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel    ! relative ice / frazil velocity
110      !
111      REAL(wp) :: zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used)
112      !!-----------------------------------------------------------------------!
113
114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
115
116      at_i(:,:) = SUM( a_i, dim=3 )
117      !------------------------------------------------------------------------------!
118      ! 1) Collection thickness of ice formed in leads and polynyas
119      !------------------------------------------------------------------------------!   
120      ! ht_i_new is the thickness of new ice formed in open water
121      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
122      ! Frazil ice forms in open water, is transported by wind
123      ! accumulates at the edge of the consolidated ice edge
124      ! where it forms aggregates of a specific thickness called
125      ! collection thickness.
126
127      zvrel(:,:) = 0._wp
128
129      ! Default new ice thickness
130      WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast
131      ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp
132      END WHERE
133
134      IF( ln_frazil ) THEN
135         !
136         ht_i_new(:,:) = 0._wp
137         !
138         ! Physical constants
139         zhicrit = 0.04                                          ! frazil ice thickness
140         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoi ) )  ! reduced grav
141         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
142         zgamafr = 0.03
143         !
144         DO jj = 2, jpjm1
145            DO ji = 2, jpim1
146               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
147                  ! -- Wind stress -- !
148                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
149                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
150                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
151                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
152                  ! Square root of wind stress
153                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
154
155                  ! -- Frazil ice velocity -- !
156                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
157                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
158                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
159
160                  ! -- Pack ice velocity -- !
161                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
162                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
163
164                  ! -- Relative frazil/pack ice velocity -- !
165                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
166                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
167                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
168                  zvrel(ji,jj) = SQRT( zvrel2 )
169
170                  ! -- new ice thickness (iterative loop) -- !
171                  ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
172                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
173
174                  iter = 1
175                  DO WHILE ( iter < 20 ) 
176                     zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
177                        &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
178                     zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
179
180                     ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
181                     iter = iter + 1
182                  END DO
183                  !
184                  ! bound ht_i_new (though I don't see why it should be necessary)
185                  ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
186                  !
187               ENDIF
188               !
189            END DO
190         END DO 
191         !
192         CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1.  )
193
194      ENDIF
195
196      !------------------------------------------------------------------------------!
197      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
198      !------------------------------------------------------------------------------!
199      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
200
201      ! Identify grid points where new ice forms
202      npti = 0   ;   nptidx(:) = 0
203      DO jj = 1, jpj
204         DO ji = 1, jpi
205            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
206               npti = npti + 1
207               nptidx( npti ) = (jj - 1) * jpi + ji
208            ENDIF
209         END DO
210      END DO
211
212      ! Move from 2-D to 1-D vectors
213      IF ( npti > 0 ) THEN
214
215         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)      , at_i        )
216         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
217         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
218         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
219         DO jl = 1, jpl
220            DO jk = 1, nlay_i
221               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
222            END DO
223         END DO
224         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      )
225         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       )
226         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    )
227         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    )
228         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   )
229         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      )
230
231         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
232         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw    )
233         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
234         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
235
236         ! Convert units for ice internal energy
237         DO jl = 1, jpl
238            DO jk = 1, nlay_i               
239               WHERE( v_i_2d(1:npti,jl) > 0._wp )
240                  ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
241               ELSEWHERE
242                  ze_i_2d(1:npti,jk,jl) = 0._wp
243               END WHERE
244            END DO
245         END DO
246
247         ! Keep old ice areas and volume in memory
248         zv_b(1:npti,:) = v_i_2d(1:npti,:) 
249         za_b(1:npti,:) = a_i_2d(1:npti,:)
250
251         ! --- Salinity of new ice --- !
252         SELECT CASE ( nn_icesal )
253         CASE ( 1 )                    ! Sice = constant
254            zs_newice(1:npti) = rn_icesal
255         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
256            DO ji = 1, npti
257               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
258            END DO
259         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
260            zs_newice(1:npti) =   2.3
261         END SELECT
262
263         ! --- Heat content of new ice --- !
264         ! We assume that new ice is formed at the seawater freezing point
265         DO ji = 1, npti
266            ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
267            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
268               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
269               &                      - rcp   *         ztmelts )
270         END DO
271
272         ! --- Age of new ice --- !
273         zo_newice(1:npti) = 0._wp
274
275         ! --- Volume of new ice --- !
276         DO ji = 1, npti
277
278            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg]
279
280            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
281                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
282                                                                   
283            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
284                                             
285            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)
286                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
287            zv_newice(ji) = - zfmdt * r1_rhoi
288
289            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
290
291            ! Contribution to heat flux to the ocean [W.m-2], >0 
292            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
293            ! Total heat flux used in this process [W.m-2] 
294            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
295            ! mass flux
296            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_rdtice
297            ! salt flux
298            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_rdtice
299         END DO
300         
301         zv_frazb(1:npti) = 0._wp
302         IF( ln_frazil ) THEN
303            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
304            DO ji = 1, npti
305               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
306               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
307               zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
308               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
309            END DO
310         END IF
311         
312         ! --- Area of new ice --- !
313         DO ji = 1, npti
314            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
315         END DO
316
317         !------------------------------------------------------------------------------!
318         ! 3) Redistribute new ice area and volume into ice categories                  !
319         !------------------------------------------------------------------------------!
320
321         ! --- lateral ice growth --- !
322         ! If lateral ice growth gives an ice concentration > amax, then
323         ! we keep the excessive volume in memory and attribute it later to bottom accretion
324         DO ji = 1, npti
325            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
326               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
327               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
328               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
329               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
330            ELSE
331               zda_res(ji) = 0._wp
332               zdv_res(ji) = 0._wp
333            ENDIF
334         END DO
335
336         ! find which category to fill
337         DO jl = 1, jpl
338            DO ji = 1, npti
339               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
340                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
341                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
342                  jcat(ji) = jl
343               ENDIF
344            END DO
345         END DO
346         at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
347
348         ! Heat content
349         DO ji = 1, npti
350            jl = jcat(ji)                                                    ! categroy in which new ice is put
351            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
352         END DO
353
354         DO jk = 1, nlay_i
355            DO ji = 1, npti
356               jl = jcat(ji)
357               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
358               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
359                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
360                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
361            END DO
362         END DO
363
364         ! --- bottom ice growth + ice enthalpy remapping --- !
365         DO jl = 1, jpl
366
367            ! for remapping
368            h_i_old (1:npti,0:nlay_i+1) = 0._wp
369            eh_i_old(1:npti,0:nlay_i+1) = 0._wp
370            DO jk = 1, nlay_i
371               DO ji = 1, npti
372                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
373                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
374               END DO
375            END DO
376
377            ! new volumes including lateral/bottom accretion + residual
378            DO ji = 1, npti
379               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
380               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
381               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
382               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
383               ! for remapping
384               h_i_old (ji,nlay_i+1) = zv_newfra
385               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
386            END DO
387            ! --- Ice enthalpy remapping --- !
388            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
389         END DO
390
391         ! --- Update salinity --- !
392         DO jl = 1, jpl
393            DO ji = 1, npti
394               sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
395            END DO
396         END DO
397
398         ! Change units for e_i
399         DO jl = 1, jpl
400            DO jk = 1, nlay_i
401               ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 
402            END DO
403         END DO
404
405         ! Move 2D vectors to 1D vectors
406         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
407         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
408         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
409          DO jl = 1, jpl
410            DO jk = 1, nlay_i
411               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
412            END DO
413         END DO
414         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
415         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
416         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
417         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
418         !
419      ENDIF ! npti > 0
420      !
421      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
422      !
423   END SUBROUTINE ice_thd_do
424
425
426   SUBROUTINE ice_thd_do_init
427      !!-----------------------------------------------------------------------
428      !!                   ***  ROUTINE ice_thd_do_init ***
429      !!                 
430      !! ** Purpose :   Physical constants and parameters associated with
431      !!                ice growth in the leads
432      !!
433      !! ** Method  :   Read the namthd_do namelist and check the parameters
434      !!                called at the first timestep (nit000)
435      !!
436      !! ** input   :   Namelist namthd_do
437      !!-------------------------------------------------------------------
438      INTEGER  ::   ios   ! Local integer
439      !!
440      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
441      !!-------------------------------------------------------------------
442      !
443      REWIND( numnam_ice_ref )              ! Namelist namthd_do in reference namelist : Ice thermodynamics
444      READ  ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
445901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_do in reference namelist', lwp )
446      REWIND( numnam_ice_cfg )              ! Namelist namthd_do in configuration namelist : Ice thermodynamics
447      READ  ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
448902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_do in configuration namelist', lwp )
449      IF(lwm) WRITE( numoni, namthd_do )
450      !
451      IF(lwp) THEN                          ! control print
452         WRITE(numout,*)
453         WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
454         WRITE(numout,*) '~~~~~~~~~~~~~~~'
455         WRITE(numout,*) '   Namelist namthd_do:'
456         WRITE(numout,*) '      ice thickness for lateral accretion                       rn_hinew   = ', rn_hinew
457         WRITE(numout,*) '      Frazil ice thickness as a function of wind or not         ln_frazil  = ', ln_frazil
458         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom     rn_maxfraz = ', rn_maxfraz
459         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil   rn_vfraz   = ', rn_vfraz
460         WRITE(numout,*) '      Squeezing coefficient for collection of frazil            rn_Cfraz   = ', rn_Cfraz
461      ENDIF
462      !
463      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
464      !
465   END SUBROUTINE ice_thd_do_init
466   
467#else
468   !!----------------------------------------------------------------------
469   !!   Default option                                NO SI3 sea-ice model
470   !!----------------------------------------------------------------------
471#endif
472
473   !!======================================================================
474END MODULE icethd_do
Note: See TracBrowser for help on using the repository browser.