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/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90 @ 5422

Last change on this file since 5422 was 5422, checked in by deazer, 9 years ago

Pre removal of svn keywords

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