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_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE/icethd_do.F90 @ 11954

Last change on this file since 11954 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • 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      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft )
116
117      at_i(:,:) = SUM( a_i, dim=3 )
118      !------------------------------------------------------------------------------!
119      ! 1) Collection thickness of ice formed in leads and polynyas
120      !------------------------------------------------------------------------------!   
121      ! ht_i_new is the thickness of new ice formed in open water
122      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
123      ! Frazil ice forms in open water, is transported by wind
124      ! accumulates at the edge of the consolidated ice edge
125      ! where it forms aggregates of a specific thickness called
126      ! collection thickness.
127
128      zvrel(:,:) = 0._wp
129
130      ! Default new ice thickness
131      WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast
132      ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp
133      END WHERE
134
135      IF( ln_frazil ) THEN
136         !
137         ht_i_new(:,:) = 0._wp
138         !
139         ! Physical constants
140         zhicrit = 0.04                                          ! frazil ice thickness
141         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoi ) )  ! reduced grav
142         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
143         zgamafr = 0.03
144         !
145         DO jj = 2, jpjm1
146            DO ji = 2, jpim1
147               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
148                  ! -- Wind stress -- !
149                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
150                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
151                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
152                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
153                  ! Square root of wind stress
154                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
155
156                  ! -- Frazil ice velocity -- !
157                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
158                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
159                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
160
161                  ! -- Pack ice velocity -- !
162                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
163                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
164
165                  ! -- Relative frazil/pack ice velocity -- !
166                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
167                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
168                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
169                  zvrel(ji,jj) = SQRT( zvrel2 )
170
171                  ! -- new ice thickness (iterative loop) -- !
172                  ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
173                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
174
175                  iter = 1
176                  DO WHILE ( iter < 20 ) 
177                     zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
178                        &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
179                     zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
180
181                     ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
182                     iter = iter + 1
183                  END DO
184                  !
185                  ! bound ht_i_new (though I don't see why it should be necessary)
186                  ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
187                  !
188               ENDIF
189               !
190            END DO
191         END DO 
192         !
193         CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1.  )
194
195      ENDIF
196
197      !------------------------------------------------------------------------------!
198      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
199      !------------------------------------------------------------------------------!
200      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
201
202      ! Identify grid points where new ice forms
203      npti = 0   ;   nptidx(:) = 0
204      DO jj = 1, jpj
205         DO ji = 1, jpi
206            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
207               npti = npti + 1
208               nptidx( npti ) = (jj - 1) * jpi + ji
209            ENDIF
210         END DO
211      END DO
212
213      ! Move from 2-D to 1-D vectors
214      IF ( npti > 0 ) THEN
215
216         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)      , at_i        )
217         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
218         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
219         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
220         DO jl = 1, jpl
221            DO jk = 1, nlay_i
222               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
223            END DO
224         END DO
225         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      )
226         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       )
227         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    )
228         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    )
229         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   )
230         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      )
231
232         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
233         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw    )
234         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
235         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
236
237         ! Convert units for ice internal energy
238         DO jl = 1, jpl
239            DO jk = 1, nlay_i               
240               WHERE( v_i_2d(1:npti,jl) > 0._wp )
241                  ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
242               ELSEWHERE
243                  ze_i_2d(1:npti,jk,jl) = 0._wp
244               END WHERE
245            END DO
246         END DO
247
248         ! Keep old ice areas and volume in memory
249         zv_b(1:npti,:) = v_i_2d(1:npti,:) 
250         za_b(1:npti,:) = a_i_2d(1:npti,:)
251
252         ! --- Salinity of new ice --- !
253         SELECT CASE ( nn_icesal )
254         CASE ( 1 )                    ! Sice = constant
255            zs_newice(1:npti) = rn_icesal
256         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
257            DO ji = 1, npti
258               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
259            END DO
260         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
261            zs_newice(1:npti) =   2.3
262         END SELECT
263
264         ! --- Heat content of new ice --- !
265         ! We assume that new ice is formed at the seawater freezing point
266         DO ji = 1, npti
267            ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
268            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
269               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
270               &                      - rcp   *         ztmelts )
271         END DO
272
273         ! --- Age of new ice --- !
274         zo_newice(1:npti) = 0._wp
275
276         ! --- Volume of new ice --- !
277         DO ji = 1, npti
278
279            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg]
280
281            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
282                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
283                                                                   
284            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
285                                             
286            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)
287                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
288            zv_newice(ji) = - zfmdt * r1_rhoi
289
290            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
291
292            ! Contribution to heat flux to the ocean [W.m-2], >0 
293            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
294            ! Total heat flux used in this process [W.m-2] 
295            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
296            ! mass flux
297            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_rdtice
298            ! salt flux
299            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_rdtice
300         END DO
301         
302         zv_frazb(1:npti) = 0._wp
303         IF( ln_frazil ) THEN
304            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
305            DO ji = 1, npti
306               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
307               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
308               zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
309               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
310            END DO
311         END IF
312         
313         ! --- Area of new ice --- !
314         DO ji = 1, npti
315            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
316         END DO
317
318         !------------------------------------------------------------------------------!
319         ! 3) Redistribute new ice area and volume into ice categories                  !
320         !------------------------------------------------------------------------------!
321
322         ! --- lateral ice growth --- !
323         ! If lateral ice growth gives an ice concentration > amax, then
324         ! we keep the excessive volume in memory and attribute it later to bottom accretion
325         DO ji = 1, npti
326            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
327               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
328               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
329               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
330               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
331            ELSE
332               zda_res(ji) = 0._wp
333               zdv_res(ji) = 0._wp
334            ENDIF
335         END DO
336
337         ! find which category to fill
338         DO jl = 1, jpl
339            DO ji = 1, npti
340               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
341                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
342                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
343                  jcat(ji) = jl
344               ENDIF
345            END DO
346         END DO
347         at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
348
349         ! Heat content
350         DO ji = 1, npti
351            jl = jcat(ji)                                                    ! categroy in which new ice is put
352            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
353         END DO
354
355         DO jk = 1, nlay_i
356            DO ji = 1, npti
357               jl = jcat(ji)
358               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
359               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
360                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
361                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
362            END DO
363         END DO
364
365         ! --- bottom ice growth + ice enthalpy remapping --- !
366         DO jl = 1, jpl
367
368            ! for remapping
369            h_i_old (1:npti,0:nlay_i+1) = 0._wp
370            eh_i_old(1:npti,0:nlay_i+1) = 0._wp
371            DO jk = 1, nlay_i
372               DO ji = 1, npti
373                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
374                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
375               END DO
376            END DO
377
378            ! new volumes including lateral/bottom accretion + residual
379            DO ji = 1, npti
380               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
381               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
382               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
383               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
384               ! for remapping
385               h_i_old (ji,nlay_i+1) = zv_newfra
386               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
387            END DO
388            ! --- Ice enthalpy remapping --- !
389            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
390         END DO
391
392         ! --- Update salinity --- !
393         DO jl = 1, jpl
394            DO ji = 1, npti
395               sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
396            END DO
397         END DO
398
399         ! Change units for e_i
400         DO jl = 1, jpl
401            DO jk = 1, nlay_i
402               ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 
403            END DO
404         END DO
405
406         ! Move 2D vectors to 1D vectors
407         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
408         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
409         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
410          DO jl = 1, jpl
411            DO jk = 1, nlay_i
412               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
413            END DO
414         END DO
415         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
416         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
417         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
418         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
419         !
420      ENDIF ! npti > 0
421      !
422      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
423      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
424      !
425   END SUBROUTINE ice_thd_do
426
427
428   SUBROUTINE ice_thd_do_init
429      !!-----------------------------------------------------------------------
430      !!                   ***  ROUTINE ice_thd_do_init ***
431      !!                 
432      !! ** Purpose :   Physical constants and parameters associated with
433      !!                ice growth in the leads
434      !!
435      !! ** Method  :   Read the namthd_do namelist and check the parameters
436      !!                called at the first timestep (nit000)
437      !!
438      !! ** input   :   Namelist namthd_do
439      !!-------------------------------------------------------------------
440      INTEGER  ::   ios   ! Local integer
441      !!
442      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
443      !!-------------------------------------------------------------------
444      !
445      READ  ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
446901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_do in reference namelist' )
447      READ  ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
448902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
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.