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.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 36.7 KB
Line 
1MODULE limthd
2   !!======================================================================
3   !!                  ***  MODULE limthd   ***
4   !!  LIM-3 :   ice thermodynamic
5   !!======================================================================
6   !! History :  LIM  ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
7   !!            2.0  ! 2002-07 (C. Ethe, G. Madec)  LIM-2 (F90 rewriting)
8   !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations
9   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif
10   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
11   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
12   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
14   !!----------------------------------------------------------------------
15#if defined key_lim3
16   !!----------------------------------------------------------------------
17   !!   'key_lim3'                                      LIM3 sea-ice model
18   !!----------------------------------------------------------------------
19   !!   lim_thd       : thermodynamic of sea ice
20   !!   lim_thd_init  : initialisation of sea-ice thermodynamic
21   !!----------------------------------------------------------------------
22   USE phycst         ! physical constants
23   USE dom_oce        ! ocean space and time domain variables
24   USE ice            ! LIM: sea-ice variables
25   USE sbc_oce        ! Surface boundary condition: ocean fields
26   USE sbc_ice        ! Surface boundary condition: ice fields
27   USE thd_ice        ! LIM thermodynamic sea-ice variables
28   USE dom_ice        ! LIM sea-ice domain
29   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion
30   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation
31   USE limthd_sal     ! LIM: thermodynamics, ice salinity
32   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution
33   USE limthd_lac     ! LIM-3 lateral accretion
34   USE limitd_th      ! remapping thickness distribution
35   USE limtab         ! LIM: 1D <==> 2D transformation
36   USE limvar         ! LIM: sea-ice variables
37   USE lbclnk         ! lateral boundary condition - MPP links
38   USE lib_mpp        ! MPP library
39   USE wrk_nemo       ! work arrays
40   USE in_out_manager ! I/O manager
41   USE prtctl         ! Print control
42   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
43   USE timing         ! Timing
44   USE limcons        ! conservation tests
45   USE limctl
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   lim_thd         ! called by limstp module
51   PUBLIC   lim_thd_init    ! called by sbc_lim_init
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE lim_thd( kt )
64      !!-------------------------------------------------------------------
65      !!                ***  ROUTINE lim_thd  ***       
66      !! 
67      !! ** Purpose : This routine manages ice thermodynamics
68      !!         
69      !! ** Action : - Initialisation of some variables
70      !!             - Some preliminary computation (oceanic heat flux
71      !!               at the ice base, snow acc.,heat budget of the leads)
72      !!             - selection of the icy points and put them in an array
73      !!             - call lim_thd_dif  for vertical heat diffusion
74      !!             - call lim_thd_dh   for vertical ice growth and melt
75      !!             - call lim_thd_ent  for enthalpy remapping
76      !!             - call lim_thd_sal  for ice desalination
77      !!             - call lim_thd_temp to  retrieve temperature from ice enthalpy
78      !!             - back to the geographic grid
79      !!     
80      !! ** References :
81      !!---------------------------------------------------------------------
82      INTEGER, INTENT(in) :: kt    ! number of iteration
83      !!
84      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
85      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
86      INTEGER  :: ii, ij           ! temporary dummy loop index
87      REAL(wp) :: zfric_u, zqld, zqfr
88      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
89      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
90      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
91      !
92      !!-------------------------------------------------------------------
93
94      IF( nn_timing == 1 )  CALL timing_start('limthd')
95
96      ! conservation test
97      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
98
99      CALL lim_var_glo2eqv
100      !------------------------------------------------------------------------!
101      ! 1) Initialization of some variables                                    !
102      !------------------------------------------------------------------------!
103      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice
104
105      !--------------------
106      ! 1.2) Heat content   
107      !--------------------
108      ! Change the units of heat content; from J/m2 to J/m3
109      DO jl = 1, jpl
110         DO jk = 1, nlay_i
111            DO jj = 1, jpj
112               DO ji = 1, jpi
113                  !0 if no ice and 1 if yes
114                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )
115                  !Energy of melting q(S,T) [J.m-3]
116                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
117               END DO
118            END DO
119         END DO
120         DO jk = 1, nlay_s
121            DO jj = 1, jpj
122               DO ji = 1, jpi
123                  !0 if no ice and 1 if yes
124                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  )
125                  !Energy of melting q(S,T) [J.m-3]
126                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
127               END DO
128            END DO
129         END DO
130      END DO
131
132      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      !
133      !-----------------------------------------------------------------------------!
134      DO jj = 1, jpj
135         DO ji = 1, jpi
136            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
137            !
138            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
139            !           !  practically no "direct lateral ablation"
140            !           
141            !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean
142            !           !  temperature and turbulent mixing (McPhee, 1992)
143            !
144            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
145            zqld =  tmask(ji,jj,1) * rdt_ice *  &
146               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
147
148            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
149            zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
150
151            ! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
152            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
153            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
154            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
155            ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
156            !                        the freezing point, so that we do not have SST < T_freeze
157            !                        This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
158
159            !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
160            qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
161
162            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
163            IF( zqld > 0._wp ) THEN
164               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90
165               qlead(ji,jj) = 0._wp
166            ELSE
167               fhld (ji,jj) = 0._wp
168            ENDIF
169            !
170            ! -----------------------------------------
171            ! Net heat flux on top of ice-ocean [W.m-2]
172            ! -----------------------------------------
173            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
174
175            ! -----------------------------------------------------------------------------
176            ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
177            ! -----------------------------------------------------------------------------
178            !     First  step here              :  non solar + precip - qlead - qturb
179            !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema)
180            !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar
181            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean               
182               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation
183               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence
184               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt
185                                                                                   !    (fhld should be 0 while bott growth)
186         END DO
187      END DO
188
189      !------------------------------------------------------------------------------!
190      ! 3) Select icy points and fulfill arrays for the vectorial grid.           
191      !------------------------------------------------------------------------------!
192
193      DO jl = 1, jpl      !loop over ice categories
194
195         IF( kt == nit000 .AND. lwp ) THEN
196            WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl 
197            WRITE(numout,*) ' ~~~~~~~~'
198         ENDIF
199
200         nbpb = 0
201         DO jj = 1, jpj
202            DO ji = 1, jpi
203               IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
204                  nbpb      = nbpb  + 1
205                  npb(nbpb) = (jj - 1) * jpi + ji
206               ENDIF
207            END DO
208         END DO
209
210         ! debug point to follow
211         jiindex_1d = 0
212         IF( ln_icectl ) THEN
213            DO ji = mi0(iiceprt), mi1(iiceprt)
214               DO jj = mj0(jiceprt), mj1(jiceprt)
215                  jiindex_1d = (jj - 1) * jpi + ji
216                  WRITE(numout,*) ' lim_thd : Category no : ', jl 
217               END DO
218            END DO
219         ENDIF
220
221         !------------------------------------------------------------------------------!
222         ! 4) Thermodynamic computation
223         !------------------------------------------------------------------------------!
224
225         IF( lk_mpp )   CALL mpp_ini_ice( nbpb , numout )
226
227         IF( nbpb > 0 ) THEN  ! If there is no ice, do nothing.
228
229            !-------------------------!
230            ! --- Move to 1D arrays ---
231            !-------------------------!
232            CALL lim_thd_1d2d( nbpb, jl, 1 )
233
234            !--------------------------------------!
235            ! --- Ice/Snow Temperature profile --- !
236            !--------------------------------------!
237            CALL lim_thd_dif( 1, nbpb )
238
239            !---------------------------------!
240            ! --- Ice/Snow thickness ---      !
241            !---------------------------------!
242            CALL lim_thd_dh( 1, nbpb )   
243
244            ! --- Ice enthalpy remapping --- !
245            CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 
246                                           
247            !---------------------------------!
248            ! --- Ice salinity ---            !
249            !---------------------------------!
250            CALL lim_thd_sal( 1, nbpb )   
251
252            !---------------------------------!
253            ! --- temperature update ---      !
254            !---------------------------------!
255            CALL lim_thd_temp( 1, nbpb )
256
257            !------------------------------------!
258            ! --- lateral melting if monocat --- !
259            !------------------------------------!
260            IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
261               CALL lim_thd_lam( 1, nbpb )
262            END IF
263
264            !-------------------------!
265            ! --- Move to 2D arrays ---
266            !-------------------------!
267            CALL lim_thd_1d2d( nbpb, jl, 2 )
268
269            !
270            IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
271         ENDIF
272         !
273      END DO !jl
274
275      !------------------------------------------------------------------------------!
276      ! 5) Global variables, diagnostics
277      !------------------------------------------------------------------------------!
278
279      !------------------------
280      ! Ice heat content             
281      !------------------------
282      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
283      DO jl = 1, jpl
284         DO jk = 1, nlay_i
285            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
286         END DO
287      END DO
288
289      !------------------------
290      ! Snow heat content             
291      !------------------------
292      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
293      DO jl = 1, jpl
294         DO jk = 1, nlay_s
295            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
296         END DO
297      END DO
298 
299      !----------------------------------
300      ! Change thickness to volume
301      !----------------------------------
302      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
303      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
304      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
305
306      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat)
307      DO jl  = 1, jpl
308         DO jj = 1, jpj
309            DO ji = 1, jpi
310               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) )
311               oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 )
312            END DO
313         END DO
314      END DO
315
316      CALL lim_var_zapsmall
317
318      !--------------------------------------------
319      ! Diagnostic thermodynamic growth rates
320      !--------------------------------------------
321      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
322
323      IF(ln_ctl) THEN            ! Control print
324         CALL prt_ctl_info(' ')
325         CALL prt_ctl_info(' - Cell values : ')
326         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
327         CALL prt_ctl(tab2d_1=e1e2t, clinfo1=' lim_thd  : cell area :')
328         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :')
329         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :')
330         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd  : vt_s      :')
331         DO jl = 1, jpl
332            CALL prt_ctl_info(' ')
333            CALL prt_ctl_info(' - Category : ', ivar1=jl)
334            CALL prt_ctl_info('   ~~~~~~~~~~')
335            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_thd  : a_i      : ')
336            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_thd  : ht_i     : ')
337            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_thd  : ht_s     : ')
338            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_thd  : v_i      : ')
339            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_thd  : v_s      : ')
340            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_thd  : e_s      : ')
341            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_thd  : t_su     : ')
342            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_thd  : t_snow   : ')
343            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_thd  : sm_i     : ')
344            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_thd  : smv_i    : ')
345            DO jk = 1, nlay_i
346               CALL prt_ctl_info(' ')
347               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
348               CALL prt_ctl_info('   ~~~~~~~')
349               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd  : t_i      : ')
350               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd  : e_i      : ')
351            END DO
352         END DO
353      ENDIF
354      !
355      !
356      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
357
358      !------------------------------------------------------------------------------|
359      !  6) Transport of ice between thickness categories.                           |
360      !------------------------------------------------------------------------------|
361      ! Given thermodynamic growth rates, transport ice between thickness categories.
362      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
363
364      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt )
365
366      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
367
368      !------------------------------------------------------------------------------|
369      !  7) Add frazil ice growing in leads.
370      !------------------------------------------------------------------------------|
371      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
372
373      CALL lim_thd_lac
374     
375      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
376
377      ! Control print
378      IF(ln_ctl) THEN
379         CALL lim_var_glo2eqv
380
381         CALL prt_ctl_info(' ')
382         CALL prt_ctl_info(' - Cell values : ')
383         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
384         CALL prt_ctl(tab2d_1=e1e2t, clinfo1=' lim_itd_th  : cell area :')
385         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
386         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
387         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :')
388         DO jl = 1, jpl
389            CALL prt_ctl_info(' ')
390            CALL prt_ctl_info(' - Category : ', ivar1=jl)
391            CALL prt_ctl_info('   ~~~~~~~~~~')
392            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ')
393            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ')
394            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ')
395            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ')
396            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ')
397            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ')
398            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ')
399            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ')
400            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
401            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
402            DO jk = 1, nlay_i
403               CALL prt_ctl_info(' ')
404               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
405               CALL prt_ctl_info('   ~~~~~~~')
406               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
407               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
408            END DO
409         END DO
410      ENDIF
411      !
412      IF( nn_timing == 1 )  CALL timing_stop('limthd')
413
414   END SUBROUTINE lim_thd 
415
416 
417   SUBROUTINE lim_thd_temp( kideb, kiut )
418      !!-----------------------------------------------------------------------
419      !!                   ***  ROUTINE lim_thd_temp ***
420      !!                 
421      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy
422      !!
423      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
424      !!-------------------------------------------------------------------
425      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
426      !!
427      INTEGER  ::   ji, jk   ! dummy loop indices
428      REAL(wp) ::   ztmelts, zaaa, zbbb, zccc, zdiscrim  ! local scalar
429      !!-------------------------------------------------------------------
430      ! Recover ice temperature
431      DO jk = 1, nlay_i
432         DO ji = kideb, kiut
433            ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0
434            ! Conversion q(S,T) -> T (second order equation)
435            zaaa          =  cpic
436            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
437            zccc          =  lfus * ( ztmelts - rt0 )
438            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
439            t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
440           
441            ! mask temperature
442            rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
443            t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
444         END DO
445      END DO
446
447   END SUBROUTINE lim_thd_temp
448
449   SUBROUTINE lim_thd_lam( kideb, kiut )
450      !!-----------------------------------------------------------------------
451      !!                   ***  ROUTINE lim_thd_lam ***
452      !!                 
453      !! ** Purpose :   Lateral melting in case monocategory
454      !!                          ( dA = A/2h dh )
455      !!-----------------------------------------------------------------------
456      INTEGER, INTENT(in) ::   kideb, kiut        ! bounds for the spatial loop
457      INTEGER             ::   ji                 ! dummy loop indices
458      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo
459      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting
460      REAL(wp)            ::   zvi, zvs           ! ice/snow volumes
461
462      DO ji = kideb, kiut
463         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
464         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
465            zvi          = a_i_1d(ji) * ht_i_1d(ji)
466            zvs          = a_i_1d(ji) * ht_s_1d(ji)
467            ! lateral melting = concentration change
468            zhi_bef     = ht_i_1d(ji) - zdh_mel
469            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
470            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
471            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
472             ! adjust thickness
473            ht_i_1d(ji) = zvi / a_i_1d(ji)           
474            ht_s_1d(ji) = zvs / a_i_1d(ji)           
475            ! retrieve total concentration
476            at_i_1d(ji) = a_i_1d(ji)
477         END IF
478      END DO
479     
480   END SUBROUTINE lim_thd_lam
481
482   SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
483      !!-----------------------------------------------------------------------
484      !!                   ***  ROUTINE lim_thd_1d2d ***
485      !!                 
486      !! ** Purpose :   move arrays from 1d to 2d and the reverse
487      !!-----------------------------------------------------------------------
488      INTEGER, INTENT(in) ::   kn       ! 1= from 2D to 1D
489                                        ! 2= from 1D to 2D
490      INTEGER, INTENT(in) ::   nbpb     ! size of 1D arrays
491      INTEGER, INTENT(in) ::   jl       ! ice cat
492      INTEGER             ::   jk       ! dummy loop indices
493
494      SELECT CASE( kn )
495
496      CASE( 1 )
497
498         CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) )
499         CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
500         CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
501         CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
502         
503         CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
504         CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
505         DO jk = 1, nlay_s
506            CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
507            CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
508         END DO
509         DO jk = 1, nlay_i
510            CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
511            CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
512            CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
513         END DO
514         
515         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
516         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
517         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
518         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) )
519         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
520         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
521         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
522         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
523         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
524         CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) ) 
525         CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) )
526         CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) )
527         CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) )
528         
529         CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) )
530         CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
531         
532         CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
533         CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) )
534         CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) )
535         CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) )
536         CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) )
537         CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) )
538         
539         CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) )
540         CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) )
541         CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) )
542         CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) )
543         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
544         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
545         
546         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
547         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
548         CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) )
549         CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) )
550         CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) )
551         CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) )
552         CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) )
553         CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) )
554         CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) )
555         CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) )
556         CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) )
557         CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
558         CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
559
560      CASE( 2 )
561
562         CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj )
563         CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj )
564         CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj )
565         CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj )
566         CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj )
567         CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj )
568         DO jk = 1, nlay_s
569            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj)
570            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj)
571         END DO
572         DO jk = 1, nlay_i
573            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj)
574            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj)
575            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj)
576         END DO
577         CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj )
578         
579         CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj )
580         CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
581         
582         CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
583         CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj )
584         CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj )
585         CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj )
586         CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj )
587         CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj )
588         
589         CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj )
590         CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj )
591         CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj )
592         CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj )
593         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
594         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj )
595         
596         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
597         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
598         CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj )
599         CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj )
600         CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj )
601         CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj )
602         CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj )
603         CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj )
604         CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj )
605         CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj )
606         CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
607         CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
608         CALL tab_1d_2d( nbpb, hfx_err_dif   , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
609         !
610         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
611         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
612         !         
613      END SELECT
614
615   END SUBROUTINE lim_thd_1d2d
616
617
618   SUBROUTINE lim_thd_init
619      !!-----------------------------------------------------------------------
620      !!                   ***  ROUTINE lim_thd_init ***
621      !!                 
622      !! ** Purpose :   Physical constants and parameters linked to the ice
623      !!              thermodynamics
624      !!
625      !! ** Method  :   Read the namicethd namelist and check the ice-thermo
626      !!              parameter values called at the first timestep (nit000)
627      !!
628      !! ** input   :   Namelist namicether
629      !!-------------------------------------------------------------------
630      INTEGER  ::   ios                 ! Local integer output status for namelist read
631      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
632         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
633         &                nn_monocat, ln_it_qnsice
634      !!-------------------------------------------------------------------
635      !
636      IF(lwp) THEN
637         WRITE(numout,*)
638         WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
639         WRITE(numout,*) '~~~~~~~'
640      ENDIF
641      !
642      REWIND( numnam_ice_ref )              ! Namelist namicethd in reference namelist : Ice thermodynamics
643      READ  ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
644901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
645
646      REWIND( numnam_ice_cfg )              ! Namelist namicethd in configuration namelist : Ice thermodynamics
647      READ  ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
648902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
649      IF(lwm) WRITE ( numoni, namicethd )
650      !
651      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
652         nn_monocat = 0
653         IF(lwp) WRITE(numout, *) '   nn_monocat must be 0 in multi-category case '
654      ENDIF
655
656      !
657      IF(lwp) THEN                          ! control print
658         WRITE(numout,*)
659         WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation '
660         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice
661         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil
662         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb
663         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb    = ', rn_vfrazb
664         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb
665         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin 
666         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice '
667         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas
668         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
669         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nn_conv_dif  = ', nn_conv_dif
670         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif
671         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
672         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
673         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
674         WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice
675      ENDIF
676      !
677   END SUBROUTINE lim_thd_init
678
679#else
680   !!----------------------------------------------------------------------
681   !!   Default option         Dummy module          NO  LIM3 sea-ice model
682   !!----------------------------------------------------------------------
683#endif
684
685   !!======================================================================
686END MODULE limthd
Note: See TracBrowser for help on using the repository browser.