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/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 5350

Last change on this file since 5350 was 5350, checked in by hadcv, 9 years ago
Update to head of the trunk (r5344).
Property svn:keywords set to `Id`
File size: 40.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: