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_2.F90 in trunk/NEMO/LIM_SRC_2 – NEMO

source: trunk/NEMO/LIM_SRC_2/limthd_2.F90 @ 1156

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

Update Id and licence information, see ticket #210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 25.0 KB
Line 
1MODULE limthd_2
2   !!======================================================================
3   !!                  ***  MODULE limthd_2   ***
4   !!              LIM thermo ice model : ice thermodynamic
5   !!======================================================================
6   !! History :  1.0  !  00-01 (LIM)
7   !!            2.0  !  02-07 (C. Ethe, G. Madec) F90
8   !!            2.0  !  03-08 (C. Ethe)  add lim_thd_init
9   !!---------------------------------------------------------------------
10#if defined key_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_thd_2      : thermodynamic of sea ice
15   !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE phycst          ! physical constants
19   USE dom_oce         ! ocean space and time domain variables
20   USE lbclnk
21   USE in_out_manager  ! I/O manager
22   USE ice_2           ! LIM sea-ice variables
23   USE ice_oce         ! sea-ice/ocean variables
24   USE sbc_oce         !
25   USE sbc_ice         !
26   USE thd_ice_2       ! LIM thermodynamic sea-ice variables
27   USE dom_ice_2       ! LIM sea-ice domain
28   USE iceini_2
29   USE limthd_zdf_2
30   USE limthd_lac_2
31   USE limtab_2
32   USE prtctl          ! Print control
33     
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_thd_2  ! called by lim_step
38
39   REAL(wp)  ::   epsi20 = 1.e-20   ,  &  ! constant values
40      &           epsi16 = 1.e-16   ,  &
41      &           epsi04 = 1.e-04   ,  &
42      &           zzero  = 0.e0     ,  &
43      &           zone   = 1.e0
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!-------- -------------------------------------------------------------
49   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
50   !! $Id$
51   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE lim_thd_2( kt )
57      !!-------------------------------------------------------------------
58      !!                ***  ROUTINE lim_thd_2  ***       
59      !! 
60      !! ** Purpose : This routine manages the ice thermodynamic.
61      !!         
62      !! ** Action : - Initialisation of some variables
63      !!             - Some preliminary computation (oceanic heat flux
64      !!               at the ice base, snow acc.,heat budget of the leads)
65      !!             - selection of the icy points and put them in an array
66      !!             - call lim_vert_ther for vert ice thermodynamic
67      !!             - back to the geographic grid
68      !!             - selection of points for lateral accretion
69      !!             - call lim_lat_acc  for the ice accretion
70      !!             - back to the geographic grid
71      !!
72      !! References :   Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
73      !!---------------------------------------------------------------------
74      INTEGER, INTENT(in) ::   kt     ! number of iteration
75      !!
76      INTEGER  ::   ji, jj,    &   ! dummy loop indices
77         nbpb  ,               &   ! nb of icy pts for thermo. cal.
78         nbpac                     ! nb of pts for lateral accretion
79      CHARACTER (len=22) :: charout
80      REAL(wp) ::  &
81         zfric_umin = 5e-03 ,  &   ! lower bound for the friction velocity
82         zfric_umax = 2e-02        ! upper bound for the friction velocity
83      REAL(wp) ::   &
84         zinda              ,  &   ! switch for test. the val. of concen.
85         zindb, zindg       ,  &   ! switches for test. the val of arg
86         za , zh, zthsnice  ,  &
87         zfric_u            ,  &   ! friction velocity
88         zfnsol             ,  &   ! total non solar heat
89         zfontn             ,  &   ! heat flux from snow thickness
90         zfntlat, zpareff          ! test. the val. of lead heat budget
91      REAL(wp), DIMENSION(jpi,jpj) ::   zhicifp,   &  ! ice thickness for outputs
92         &                              zqlbsbq       ! link with lead energy budget qldif
93      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmsk      ! working array
94      !!-------------------------------------------------------------------
95
96      IF( kt == nit000 )   CALL lim_thd_init_2  ! Initialization (first time-step only)
97   
98      !-------------------------------------------!
99      !   Initilization of diagnostic variables   !
100      !-------------------------------------------!
101     
102!i est-ce utile?  oui au moins en partie
103      rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
104      rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
105      fdvolif(:,:) = 0.e0   ! total variation of ice volume
106      rdvonif(:,:) = 0.e0   ! lateral variation of ice volume
107      fstric (:,:) = 0.e0   ! part of solar radiation absorbing inside the ice
108      fscmbq (:,:) = 0.e0   ! linked with fstric
109      ffltbif(:,:) = 0.e0   ! linked with fstric
110      qfvbq  (:,:) = 0.e0   ! linked with fstric
111      rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area
112      rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area
113      hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.
114      zmsk (:,:,:) = 0.e0
115
116      DO jj = 1, jpj
117         DO ji = 1, jpi
118            hsnif(ji,jj)  = hsnif(ji,jj) *  MAX( zzero, SIGN( zone , hsnif(ji,jj) - epsi04 ) )
119         END DO
120      END DO
121
122      IF(ln_ctl)   CALL prt_ctl(tab2d_1=hsnif     , clinfo1=' lim_thd: hsnif   : ')
123     
124      !-----------------------------------!
125      !   Treatment of particular cases   !
126      !-----------------------------------!
127     
128      DO jj = 1, jpj
129         DO ji = 1, jpi
130            !  snow is transformed into ice if the original ice cover disappears.
131            zindg         = tms(ji,jj) *  MAX( zzero , SIGN( zone , -hicif(ji,jj) ) )
132            hicif(ji,jj)  = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0
133            hsnif(ji,jj)  = ( zone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
134            dmgwi(ji,jj)  = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj)   ! snow/ice mass
135           
136            !  the lead fraction, frld, must be little than or equal to amax (ice ridging).
137            zthsnice      = hsnif(ji,jj) + hicif(ji,jj)
138            zindb         = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
139            za            = zindb * MIN( zone, ( 1.0 - frld(ji,jj) ) * uscomi )
140            hsnif (ji,jj) = hsnif(ji,jj)  * za
141            hicif (ji,jj) = hicif(ji,jj)  * za
142            qstoif(ji,jj) = qstoif(ji,jj) * za
143            frld  (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za , epsi20 )
144           
145            !  the in situ ice thickness, hicif, must be equal to or greater than hiclim.
146            zh            = MAX( zone , zindb * hiclim  / MAX( hicif(ji,jj) , epsi20 ) )
147            hsnif (ji,jj) = hsnif(ji,jj)  * zh
148            hicif (ji,jj) = hicif(ji,jj)  * zh
149            qstoif(ji,jj) = qstoif(ji,jj) * zh
150            frld  (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh
151         END DO
152      END DO
153
154      IF(ln_ctl) THEN
155         CALL prt_ctl(tab2d_1=hicif  , clinfo1=' lim_thd: hicif   : ')
156         CALL prt_ctl(tab2d_1=hsnif  , clinfo1=' lim_thd: hsnif   : ')
157         CALL prt_ctl(tab2d_1=dmgwi  , clinfo1=' lim_thd: dmgwi   : ')
158         CALL prt_ctl(tab2d_1=qstoif , clinfo1=' lim_thd: qstoif  : ')
159         CALL prt_ctl(tab2d_1=frld   , clinfo1=' lim_thd: frld    : ')
160      ENDIF
161
162     
163      !-------------------------------!
164      !   Thermodynamics of sea ice   !
165      !-------------------------------!
166     
167      !      Partial computation of forcing for the thermodynamic sea ice model.
168      !--------------------------------------------------------------------------
169
170      sst_m(:,:) = sst_m(:,:) + rt0
171
172      !CDIR NOVERRCHK
173      DO jj = 1, jpj
174         !CDIR NOVERRCHK
175         DO ji = 1, jpi
176            zthsnice       = hsnif(ji,jj) + hicif(ji,jj)
177            zindb          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
178            pfrld(ji,jj)   = frld(ji,jj)
179            zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
180           
181            !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
182            thcm(ji,jj)    = 0.e0 
183           
184            !  net downward heat flux from the ice to the ocean, expressed as a function of ocean
185            !  temperature and turbulent mixing (McPhee, 1992)
186            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity
187            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) 
188            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
189                       
190            !  partial computation of the lead energy budget (qldif)
191            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting
192            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean
193            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
194               &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
195               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
196               &                               * frld(ji,jj) * rdt_ice   
197            !  parlat : percentage of energy used for lateral ablation (0.0)
198            zfntlat        = 1.0 - MAX( zzero , SIGN( zone ,  - qldif(ji,jj) ) )
199            zpareff        = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat
200            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 )
201            qldif  (ji,jj) = zpareff *  qldif(ji,jj)
202            qdtcn  (ji,jj) = zpareff * qdtcn(ji,jj)
203           
204            !  energy needed to bring ocean surface layer until its freezing
205            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
206           
207            !  calculate oceanic heat flux.
208            fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) )
209           
210            ! computation of the daily thermodynamic ice production (only needed for output)
211            zhicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) )
212         END DO
213      END DO
214     
215      sst_m(:,:) = sst_m(:,:) - rt0
216     
217      !         Select icy points and fulfill arrays for the vectorial grid.
218      !----------------------------------------------------------------------
219      nbpb = 0
220      DO jj = 1, jpj
221         DO ji = 1, jpi
222            IF ( frld(ji,jj) < 1.0 ) THEN     
223               nbpb      = nbpb + 1
224               npb(nbpb) = (jj - 1) * jpi + ji
225            ENDIF
226         END DO
227      END DO
228
229      IF(ln_ctl) THEN
230         CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld   : ', tab2d_2=thcm   , clinfo2='  thcm    : ')
231         CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn  , clinfo2='  qdtcn   : ')
232         CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif   : ', tab2d_2=zqlbsbq, clinfo2='  zqlbsbq : ')
233         CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif   : ', tab2d_2=fbif   , clinfo2='  fbif    : ')
234         zmsk(:,:,1) = tms(:,:)
235         CALL prt_ctl(tab2d_1=qcmif  , clinfo1=' lim_thd: qcmif   : ', mask1=zmsk)
236         CALL prt_ctl(tab2d_1=zhicifp, clinfo1=' lim_thd: zhicifp : ')
237         WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb
238         CALL prt_ctl_info(charout)
239      ENDIF
240     
241     
242      ! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation
243      !------------------------------------------------------------------------------------
244
245      IF ( nbpb > 0) THEN
246         
247         !  put the variable in a 1-D array for thermodynamics process
248         CALL tab_2d_1d_2( nbpb, frld_1d    (1:nbpb)     , frld       , jpi, jpj, npb(1:nbpb) )
249         CALL tab_2d_1d_2( nbpb, h_ice_1d   (1:nbpb)     , hicif      , jpi, jpj, npb(1:nbpb) )
250         CALL tab_2d_1d_2( nbpb, h_snow_1d  (1:nbpb)     , hsnif      , jpi, jpj, npb(1:nbpb) )
251         CALL tab_2d_1d_2( nbpb, sist_1d    (1:nbpb)     , sist       , jpi, jpj, npb(1:nbpb) )
252         CALL tab_2d_1d_2( nbpb, tbif_1d    (1:nbpb , 1 ), tbif(:,:,1), jpi, jpj, npb(1:nbpb) )
253         CALL tab_2d_1d_2( nbpb, tbif_1d    (1:nbpb , 2 ), tbif(:,:,2), jpi, jpj, npb(1:nbpb) )
254         CALL tab_2d_1d_2( nbpb, tbif_1d    (1:nbpb , 3 ), tbif(:,:,3), jpi, jpj, npb(1:nbpb) )
255         CALL tab_2d_1d_2( nbpb, qsr_ice_1d (1:nbpb)     , qsr_ice    , jpi, jpj, npb(1:nbpb) )
256         CALL tab_2d_1d_2( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) )
257         CALL tab_2d_1d_2( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) )
258         CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) )
259#if ! defined key_coupled
260         CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) )
261         CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb)     , dqla_ice   , jpi, jpj, npb(1:nbpb) )
262#endif
263         CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb)     , dqns_ice   , jpi, jpj, npb(1:nbpb) )
264         CALL tab_2d_1d_2( nbpb, tfu_1d     (1:nbpb)     , tfu        , jpi, jpj, npb(1:nbpb) )
265         CALL tab_2d_1d_2( nbpb, sprecip_1d (1:nbpb)     , sprecip    , jpi, jpj, npb(1:nbpb) ) 
266         CALL tab_2d_1d_2( nbpb, fbif_1d    (1:nbpb)     , fbif       , jpi, jpj, npb(1:nbpb) )
267         CALL tab_2d_1d_2( nbpb, thcm_1d    (1:nbpb)     , thcm       , jpi, jpj, npb(1:nbpb) )
268         CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
269         CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) )
270         CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
271         CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
272         CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
273 
274         CALL lim_thd_zdf_2( 1, nbpb )       !  compute ice growth
275         
276         !  back to the geographic grid.
277         CALL tab_1d_2d_2( nbpb, frld       , npb, frld_1d   (1:nbpb)     , jpi, jpj )
278         CALL tab_1d_2d_2( nbpb, hicif      , npb, h_ice_1d  (1:nbpb)     , jpi, jpj )
279         CALL tab_1d_2d_2( nbpb, hsnif      , npb, h_snow_1d (1:nbpb)     , jpi, jpj )
280         CALL tab_1d_2d_2( nbpb, sist       , npb, sist_1d   (1:nbpb)     , jpi, jpj )
281         CALL tab_1d_2d_2( nbpb, tbif(:,:,1), npb, tbif_1d   (1:nbpb , 1 ), jpi, jpj )   
282         CALL tab_1d_2d_2( nbpb, tbif(:,:,2), npb, tbif_1d   (1:nbpb , 2 ), jpi, jpj )   
283         CALL tab_1d_2d_2( nbpb, tbif(:,:,3), npb, tbif_1d   (1:nbpb , 3 ), jpi, jpj )   
284         CALL tab_1d_2d_2( nbpb, fscmbq     , npb, fscbq_1d  (1:nbpb)     , jpi, jpj )
285         CALL tab_1d_2d_2( nbpb, ffltbif    , npb, fltbif_1d (1:nbpb)     , jpi, jpj )
286         CALL tab_1d_2d_2( nbpb, fstric     , npb, fstbif_1d (1:nbpb)     , jpi, jpj )
287         CALL tab_1d_2d_2( nbpb, qldif      , npb, qldif_1d  (1:nbpb)     , jpi, jpj )
288         CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
289         CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj )
290         CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
291         CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
292         CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
293         CALL tab_1d_2d_2( nbpb, rdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
294         CALL tab_1d_2d_2( nbpb, rdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
295         CALL tab_1d_2d_2( nbpb, fdvolif    , npb, dvlbq_1d  (1:nbpb)     , jpi, jpj )
296         CALL tab_1d_2d_2( nbpb, rdvonif    , npb, dvnbq_1d  (1:nbpb)     , jpi, jpj ) 
297
298 
299      ENDIF
300
301     
302      !      Up-date sea ice thickness.
303      !---------------------------------
304      DO jj = 1, jpj
305         DO ji = 1, jpi
306            phicif(ji,jj) = hicif(ji,jj) 
307            hicif(ji,jj)  = hicif(ji,jj) *  ( zone -  MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) )
308         END DO
309      END DO
310
311     
312      !      Tricky trick : add 2 to frld in the Southern Hemisphere.
313      !----------------------------------------------------------
314      IF( fcor(1,1) < 0.e0 ) THEN
315         DO jj = 1, njeqm1
316            DO ji = 1, jpi
317               frld(ji,jj) = frld(ji,jj) + 2.0
318            END DO
319         END DO
320      ENDIF
321     
322     
323      !     Select points for lateral accretion (this occurs when heat exchange
324      !     between ice and ocean is negative; ocean losing heat)
325      !-----------------------------------------------------------------
326      nbpac = 0
327      DO jj = 1, jpj
328         DO ji = 1, jpi
329!i yes!     IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
330            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
331               nbpac = nbpac + 1
332               npac( nbpac ) = (jj - 1) * jpi + ji
333            ENDIF
334         END DO
335      END DO
336     
337      IF(ln_ctl) THEN
338         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=hicif, clinfo2=' hicif : ')
339         WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac
340         CALL prt_ctl_info(charout)
341      ENDIF
342
343     
344      !
345      !     If ocean gains heat do nothing ; otherwise, one performs lateral accretion
346      !--------------------------------------------------------------------------------
347
348      IF( nbpac > 0 ) THEN
349         
350         !...Put the variable in a 1-D array for lateral accretion
351         CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) )
352         CALL tab_2d_1d_2( nbpac, h_snow_1d (1:nbpac)     , hsnif      , jpi, jpj, npac(1:nbpac) )
353         CALL tab_2d_1d_2( nbpac, h_ice_1d  (1:nbpac)     , hicif      , jpi, jpj, npac(1:nbpac) )
354         CALL tab_2d_1d_2( nbpac, tbif_1d   (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) )   
355         CALL tab_2d_1d_2( nbpac, tbif_1d   (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) )   
356         CALL tab_2d_1d_2( nbpac, tbif_1d   (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) )   
357         CALL tab_2d_1d_2( nbpac, qldif_1d  (1:nbpac)     , qldif      , jpi, jpj, npac(1:nbpac) )
358         CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) )
359         CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) )
360         CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) )
361         CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , fdvolif    , jpi, jpj, npac(1:nbpac) )
362         CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) )
363       
364         !  call lateral accretion routine.
365         CALL lim_thd_lac_2( 1 , nbpac )
366         
367         !   back to the geographic grid
368         CALL tab_1d_2d_2( nbpac, frld       , npac(1:nbpac), frld_1d   (1:nbpac)     , jpi, jpj )
369         CALL tab_1d_2d_2( nbpac, hsnif      , npac(1:nbpac), h_snow_1d (1:nbpac)     , jpi, jpj )
370         CALL tab_1d_2d_2( nbpac, hicif      , npac(1:nbpac), h_ice_1d  (1:nbpac)     , jpi, jpj )
371         CALL tab_1d_2d_2( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d   (1:nbpac , 1 ), jpi, jpj )
372         CALL tab_1d_2d_2( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d   (1:nbpac , 2 ), jpi, jpj )
373         CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj )
374         CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj )
375         CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj )
376         CALL tab_1d_2d_2( nbpac, fdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj )
377       
378      ENDIF
379       
380       
381      !      Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
382      !      Update daily thermodynamic ice production.   
383      !------------------------------------------------------------------------------
384       
385      DO jj = 1, jpj
386         DO ji = 1, jpi
387            frld  (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) )
388            hicifp(ji,jj) =  hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) - zhicifp(ji,jj) + hicifp(ji,jj)
389         END DO
390      END DO
391
392      IF(ln_ctl) THEN
393         CALL prt_ctl_info(' lim_thd  end  ')
394         CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ')
395         CALL prt_ctl(tab2d_1=frld  , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ')
396         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ')
397         CALL prt_ctl(tab2d_1=sist  , clinfo1=' lim_thd: sist    : ')
398         CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ')
399         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ')
400         CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ')
401         CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ')
402         CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ')
403      ENDIF
404       !
405    END SUBROUTINE lim_thd_2
406
407
408    SUBROUTINE lim_thd_init_2
409      !!-------------------------------------------------------------------
410      !!                   ***  ROUTINE lim_thd_init_2 ***
411      !!                 
412      !! ** Purpose :   Physical constants and parameters linked to the ice
413      !!      thermodynamics
414      !!
415      !! ** Method  :   Read the namicethd namelist and check the ice-thermo
416      !!       parameter values called at the first timestep (nit000)
417      !!
418      !! ** input   :   Namelist namicether
419      !!-------------------------------------------------------------------
420      NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax  ,        &
421         &                swiqst, sbeta  , parlat, hakspl, hibspl, exld,  &
422         &                hakdif, hnzst  , thth  , parsub, alphs
423      !!-------------------------------------------------------------------
424     
425
426      ! Define the initial parameters
427      ! -------------------------
428      REWIND( numnam_ice )
429      READ  ( numnam_ice , namicethd )
430      IF(lwp) THEN
431         WRITE(numout,*)
432         WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation '
433         WRITE(numout,*)'~~~~~~~~~~~~~~'
434         WRITE(numout,*)'       maximum melting at the bottom                           hmelt        = ', hmelt
435         WRITE(numout,*)'       ice thick. for lateral accretion in NH (SH)             hiccrit(1/2) = ', hiccrit
436         WRITE(numout,*)'       ice thick. corr. to max. energy stored in brine pocket  hicmin       = ', hicmin 
437         WRITE(numout,*)'       minimum ice thickness                                   hiclim       = ', hiclim 
438         WRITE(numout,*)'       maximum lead fraction                                   amax         = ', amax
439         WRITE(numout,*)'       energy stored in brine pocket (=1) or not (=0)     swiqst       = ', swiqst 
440         WRITE(numout,*)'       numerical carac. of the scheme for diffusion in ice '
441         WRITE(numout,*)'       Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta
442         WRITE(numout,*)'       percentage of energy used for lateral ablation          parlat       = ', parlat
443         WRITE(numout,*)'       slope of distr. for Hakkinen-Mellor lateral melting     hakspl       = ', hakspl 
444         WRITE(numout,*)'       slope of distribution for Hibler lateral melting        hibspl       = ', hibspl
445         WRITE(numout,*)'       exponent for leads-closure rate                         exld         = ', exld
446         WRITE(numout,*)'       coefficient for diffusions of ice and snow              hakdif       = ', hakdif
447         WRITE(numout,*)'       threshold thick. for comp. of eq. thermal conductivity  zhth         = ', thth 
448         WRITE(numout,*)'       thickness of the surf. layer in temp. computation       hnzst        = ', hnzst
449         WRITE(numout,*)'       switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub 
450         WRITE(numout,*)'       coefficient for snow density when snow ice formation    alphs        = ', alphs
451      ENDIF
452           
453      uscomi = 1.0 / ( 1.0 - amax )   ! inverse of minimum lead fraction
454      rcdsn = hakdif * rcdsn 
455      rcdic = hakdif * rcdic
456     
457      IF ( ( hsndif > 100.e0 ) .OR. ( hicdif > 100.e0 ) ) THEN
458         cnscg = 0.e0
459      ELSE
460         cnscg = rcpsn / rcpic   ! ratio  rcpsn/rcpic
461      ENDIF
462 
463   END SUBROUTINE lim_thd_init_2
464
465#else
466   !!----------------------------------------------------------------------
467   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
468   !!----------------------------------------------------------------------
469CONTAINS
470   SUBROUTINE lim_thd_2         ! Dummy routine
471   END SUBROUTINE lim_thd_2
472#endif
473
474   !!======================================================================
475END MODULE limthd_2
Note: See TracBrowser for help on using the repository browser.