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 branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2 – NEMO

source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90 @ 1858

Last change on this file since 1858 was 1858, checked in by gm, 14 years ago

ticket:#665 : step 1 - heat content of freezing-melting ice

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