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

Last change on this file since 1649 was 1565, checked in by rblod, 15 years ago

Use appropiate e3 for divergence computation and for sea-ice, see ticket #507

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