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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 5134

Last change on this file since 5134 was 5134, checked in by clem, 9 years ago
LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)
Property svn:keywords set to `Id`
File size: 39.7 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
[5123]	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	!!-------------------------------------------------------------------
[4990]	96	CALL wrk_alloc( jpi, jpj, zqsr, zqns )
	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
[4990]	103	!------------------------------------------------------------------------!
	104	! 1) Initialization of some variables !
	105	!------------------------------------------------------------------------!
	106	ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice
[825]	107
	108	!--------------------
	109	! 1.2) Heat content
	110	!--------------------
[5123]	111	! Change the units of heat content; from J/m2 to J/m3
[825]	112	DO jl = 1, jpl
[921]	113	DO jk = 1, nlay_i
	114	DO jj = 1, jpj
	115	DO ji = 1, jpi
	116	!0 if no ice and 1 if yes
[5134]	117	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) )
[4688]	118	!Energy of melting q(S,T) [J.m-3]
[5123]	119	e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
[921]	120	END DO
[825]	121	END DO
[921]	122	END DO
	123	DO jk = 1, nlay_s
	124	DO jj = 1, jpj
	125	DO ji = 1, jpi
	126	!0 if no ice and 1 if yes
[5134]	127	rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )
[4688]	128	!Energy of melting q(S,T) [J.m-3]
[5123]	129	e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
[921]	130	END DO
[825]	131	END DO
[921]	132	END DO
[825]	133	END DO
	134
[921]	135	! 2) Partial computation of forcing for the thermodynamic sea ice model. !
	136	!-----------------------------------------------------------------------------!
[825]	137
[4990]	138	!--- Ocean solar and non solar fluxes to be used in zqld
	139	IF ( .NOT. lk_cpl ) THEN ! --- forced case, fluxes to the lead are the same as over the ocean
	140	!
	141	zqsr(:,:) = qsr(:,:) ; zqns(:,:) = qns(:,:)
	142	!
	143	ELSE ! --- coupled case, fluxes to the lead are total - intercepted
	144	!
	145	zqsr(:,:) = qsr_tot(:,:) ; zqns(:,:) = qns_tot(:,:)
	146	!
	147	DO jl = 1, jpl
	148	DO jj = 1, jpj
	149	DO ji = 1, jpi
	150	zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
	151	zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
	152	END DO
	153	END DO
	154	END DO
	155	!
	156	ENDIF
	157
[921]	158	DO jj = 1, jpj
	159	DO ji = 1, jpi
[5134]	160	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
[2528]	161	!
[921]	162	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	163	! ! practically no "direct lateral ablation"
	164	!
	165	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	166	! ! temperature and turbulent mixing (McPhee, 1992)
[4688]	167	!
[4990]	168
[4688]	169	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
[4990]	170	! REMARK valid at least in forced mode from clem
	171	! precip is included in qns but not in qns_ice
	172	IF ( lk_cpl ) THEN
[5123]	173	zqld = tmask(ji,jj,1) * rdt_ice * &
[4990]	174	& ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode
[5123]	175	& + ( pfrld(ji,jj)*rn_betas - pfrld(ji,jj) ) sprecip(ji,jj) * & ! heat content of precip
	176	& ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) &
	177	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
[4990]	178	ELSE
[5123]	179	zqld = tmask(ji,jj,1) * rdt_ice * &
[4990]	180	& ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) &
[5123]	181	& + ( pfrld(ji,jj)*rn_betas - pfrld(ji,jj) ) sprecip(ji,jj) * & ! heat content of precip
	182	& ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) &
	183	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
[4990]	184	ENDIF
[825]	185
[5123]	186	! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
	187	zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
[4688]	188
[5123]	189	! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
	190	zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin )
	191	fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
	192	fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
	193	! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
	194	! the freezing point, so that we do not have SST < T_freeze
	195	! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
	196
[4688]	197	!-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
[5123]	198	qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
[4688]	199
	200	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
[5123]	201	IF( zqld > 0._wp ) THEN
	202	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]	203	qlead(ji,jj) = 0._wp
[4990]	204	ELSE
	205	fhld (ji,jj) = 0._wp
[4688]	206	ENDIF
[2528]	207	!
[4688]	208	! -----------------------------------------
	209	! Net heat flux on top of ice-ocean [W.m-2]
	210	! -----------------------------------------
	211	! First step here : heat flux at the ocean surface + precip
	212	! Second step below : heat flux at the ice surface (after limthd_dif)
	213	hfx_in(ji,jj) = hfx_in(ji,jj) &
	214	! heat flux above the ocean
[4990]	215	& + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) &
[4688]	216	! latent heat of precip (note that precip is included in qns but not in qns_ice)
[5123]	217	& + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) &
	218	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )
[4688]	219
	220	! -----------------------------------------------------------------------------
	221	! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2]
	222	! -----------------------------------------------------------------------------
	223	! First step here : non solar + precip - qlead - qturb
	224	! Second step in limthd_dh : heat remaining if total melt (zq_rema)
	225	! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar
[4990]	226	hfx_out(ji,jj) = hfx_out(ji,jj) &
[4688]	227	! Non solar heat flux received by the ocean
[4990]	228	& + pfrld(ji,jj) * qns(ji,jj) &
[4688]	229	! latent heat of precip (note that precip is included in qns but not in qns_ice)
[5123]	230	& + ( pfrld(ji,jj)*rn_betas - pfrld(ji,jj) ) sprecip(ji,jj) &
	231	& * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) &
	232	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) &
[4688]	233	! heat flux taken from the ocean where there is open water ice formation
[4990]	234	& - qlead(ji,jj) * r1_rdtice &
[4688]	235	! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth)
[4990]	236	& - at_i(ji,jj) * fhtur(ji,jj) &
[4688]	237	& - at_i(ji,jj) * fhld(ji,jj)
	238
[825]	239	END DO
	240	END DO
	241
[921]	242	!------------------------------------------------------------------------------!
	243	! 3) Select icy points and fulfill arrays for the vectorial grid.
	244	!------------------------------------------------------------------------------!
[825]	245
	246	DO jl = 1, jpl !loop over ice categories
	247
[921]	248	IF( kt == nit000 .AND. lwp ) THEN
	249	WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl
	250	WRITE(numout,*) ' ~~~~~~~~'
	251	ENDIF
[825]	252
	253	nbpb = 0
	254	DO jj = 1, jpj
	255	DO ji = 1, jpi
[5123]	256	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
[825]	257	nbpb = nbpb + 1
	258	npb(nbpb) = (jj - 1) * jpi + ji
	259	ENDIF
	260	END DO
	261	END DO
	262
[4333]	263	! debug point to follow
	264	jiindex_1d = 0
[5128]	265	IF( ln_icectl ) THEN
	266	DO ji = mi0(iiceprt), mi1(iiceprt)
	267	DO jj = mj0(jiceprt), mj1(jiceprt)
[4333]	268	jiindex_1d = (jj - 1) * jpi + ji
[4688]	269	WRITE(numout,*) ' lim_thd : Category no : ', jl
[4333]	270	END DO
	271	END DO
	272	ENDIF
	273
[921]	274	!------------------------------------------------------------------------------!
	275	! 4) Thermodynamic computation
	276	!------------------------------------------------------------------------------!
[825]	277
[2715]	278	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
[869]	279
[1572]	280	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
[825]	281
[5123]	282	!-------------------------!
	283	! --- Move to 1D arrays ---
	284	!-------------------------!
	285	CALL lim_thd_1d2d( nbpb, jl, 1 )
[825]	286
[5123]	287	!--------------------------------------!
	288	! --- Ice/Snow Temperature profile --- !
	289	!--------------------------------------!
[4688]	290	CALL lim_thd_dif( 1, nbpb )
[921]	291
[4688]	292	!---------------------------------!
[5123]	293	! --- Ice/Snow thickness --- !
[4688]	294	!---------------------------------!
	295	CALL lim_thd_dh( 1, nbpb )
[825]	296
[4688]	297	! --- Ice enthalpy remapping --- !
[4872]	298	CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) )
[4688]	299
	300	!---------------------------------!
[5123]	301	! --- Ice salinity --- !
[4688]	302	!---------------------------------!
	303	CALL lim_thd_sal( 1, nbpb )
[825]	304
[4688]	305	!---------------------------------!
[5123]	306	! --- temperature update --- !
[4688]	307	!---------------------------------!
	308	CALL lim_thd_temp( 1, nbpb )
[825]	309
[5123]	310	!------------------------------------!
	311	! --- lateral melting if monocat --- !
	312	!------------------------------------!
	313	IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN
	314	CALL lim_thd_lam( 1, nbpb )
	315	END IF
[825]	316
[5123]	317	!-------------------------!
	318	! --- Move to 2D arrays ---
	319	!-------------------------!
	320	CALL lim_thd_1d2d( nbpb, jl, 2 )
[4688]	321
[2528]	322	!
[1572]	323	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	324	ENDIF
	325	!
	326	END DO
[825]	327
[921]	328	!------------------------------------------------------------------------------!
	329	! 5) Global variables, diagnostics
	330	!------------------------------------------------------------------------------!
[825]	331
	332	!------------------------
[5123]	333	! Ice heat content
[825]	334	!------------------------
[5123]	335	! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
[825]	336	DO jl = 1, jpl
[921]	337	DO jk = 1, nlay_i
[5123]	338	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
[1572]	339	END DO
	340	END DO
[825]	341
	342	!------------------------
[5123]	343	! Snow heat content
[825]	344	!------------------------
[5123]	345	! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
[825]	346	DO jl = 1, jpl
	347	DO jk = 1, nlay_s
[5123]	348	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
[1572]	349	END DO
	350	END DO
[5123]	351
	352	!------------------------
	353	! Ice natural aging
	354	!------------------------
	355	oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday
[825]	356
	357	!----------------------------------
[5123]	358	! Change thickness to volume
[825]	359	!----------------------------------
	360	CALL lim_var_eqv2glo
	361
[5134]	362	CALL lim_var_zapsmall
[825]	363	!--------------------------------------------
[5123]	364	! Diagnostic thermodynamic growth rates
[825]	365	!--------------------------------------------
[5128]	366	IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' ) ! control print
[5123]	367
[2528]	368	IF(ln_ctl) THEN ! Control print
[867]	369	CALL prt_ctl_info(' ')
	370	CALL prt_ctl_info(' - Cell values : ')
	371	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[5123]	372	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd : cell area :')
[863]	373	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :')
	374	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :')
	375	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd : vt_s :')
	376	DO jl = 1, jpl
[867]	377	CALL prt_ctl_info(' ')
[863]	378	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	379	CALL prt_ctl_info(' ~~~~~~~~~~')
	380	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_thd : a_i : ')
	381	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_thd : ht_i : ')
	382	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_thd : ht_s : ')
	383	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_thd : v_i : ')
	384	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_thd : v_s : ')
	385	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_thd : e_s : ')
	386	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_thd : t_su : ')
	387	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_thd : t_snow : ')
	388	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_thd : sm_i : ')
	389	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_thd : smv_i : ')
	390	DO jk = 1, nlay_i
[867]	391	CALL prt_ctl_info(' ')
[863]	392	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	393	CALL prt_ctl_info(' ~~~~~~~')
	394	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd : t_i : ')
	395	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd : e_i : ')
	396	END DO
	397	END DO
	398	ENDIF
[2528]	399	!
[4990]	400	!
	401	CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
	402
[5123]	403	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	404	!------------------------------------------------------------------------------\|
	405	! 6) Transport of ice between thickness categories. \|
	406	!------------------------------------------------------------------------------\|
	407	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	408
	409	! Given thermodynamic growth rates, transport ice between thickness categories.
	410	IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )
[4990]	411	!
[5123]	412	CALL lim_var_glo2eqv ! only for info
	413	CALL lim_var_agg(1)
	414
	415	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	416	!------------------------------------------------------------------------------\|
	417	! 7) Add frazil ice growing in leads.
	418	!------------------------------------------------------------------------------\|
[5134]	419	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[5123]	420	CALL lim_thd_lac
	421	CALL lim_var_glo2eqv ! only for info
	422
[4688]	423	! conservation test
[5123]	424	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	425
	426	IF(ln_ctl) THEN ! Control print
	427	CALL prt_ctl_info(' ')
	428	CALL prt_ctl_info(' - Cell values : ')
	429	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
	430	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th : cell area :')
	431	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :')
	432	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :')
	433	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :')
	434	DO jl = 1, jpl
	435	CALL prt_ctl_info(' ')
	436	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	437	CALL prt_ctl_info(' ~~~~~~~~~~')
	438	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ')
	439	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ')
	440	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ')
	441	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ')
	442	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ')
	443	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ')
	444	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ')
	445	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ')
	446	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ')
	447	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ')
	448	DO jk = 1, nlay_i
	449	CALL prt_ctl_info(' ')
	450	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	451	CALL prt_ctl_info(' ~~~~~~~')
	452	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ')
	453	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ')
	454	END DO
	455	END DO
	456	ENDIF
[4161]	457	!
	458	IF( nn_timing == 1 ) CALL timing_stop('limthd')
[4990]	459
[4688]	460	END SUBROUTINE lim_thd
[825]	461
[4688]	462	SUBROUTINE lim_thd_temp( kideb, kiut )
[825]	463	!!-----------------------------------------------------------------------
[4688]	464	!! * ROUTINE lim_thd_temp *
[825]	465	!!
[4688]	466	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
[825]	467	!!
	468	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	469	!!-------------------------------------------------------------------
[1572]	470	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	471	!!
[2715]	472	INTEGER :: ji, jk ! dummy loop indices
[4990]	473	REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar
[825]	474	!!-------------------------------------------------------------------
[4688]	475	! Recover ice temperature
	476	DO jk = 1, nlay_i
[825]	477	DO ji = kideb, kiut
[5123]	478	ztmelts = -tmut * s_i_1d(ji,jk) + rt0
[4688]	479	! Conversion q(S,T) -> T (second order equation)
	480	zaaa = cpic
[5123]	481	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
	482	zccc = lfus * ( ztmelts - rt0 )
[4688]	483	zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
[5123]	484	t_i_1d(ji,jk) = rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
[4688]	485
	486	! mask temperature
[4990]	487	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )
[5123]	488	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
[4688]	489	END DO
	490	END DO
[825]	491
[4688]	492	END SUBROUTINE lim_thd_temp
[825]	493
[5123]	494	SUBROUTINE lim_thd_lam( kideb, kiut )
	495	!!-----------------------------------------------------------------------
	496	!! * ROUTINE lim_thd_lam *
	497	!!
	498	!! ** Purpose : Lateral melting in case monocategory
	499	!! ( dA = A/2h dh )
	500	!!-----------------------------------------------------------------------
	501	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	502	INTEGER :: ji ! dummy loop indices
	503	REAL(wp) :: zhi_bef ! ice thickness before thermo
	504	REAL(wp) :: zdh_mel, zda_mel ! net melting
	505	REAL(wp) :: zv ! ice volume
	506
	507	DO ji = kideb, kiut
	508	zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
	509	IF( zdh_mel < 0._wp ) THEN
	510	zv = a_i_1d(ji) * ht_i_1d(ji)
	511	! lateral melting = concentration change
	512	zhi_bef = ht_i_1d(ji) - zdh_mel
	513	zda_mel = a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) )
	514	a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel )
	515	! adjust thickness
[5134]	516	rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) )
[5123]	517	ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
	518	! retrieve total concentration
	519	at_i_1d(ji) = a_i_1d(ji)
	520	END IF
	521	END DO
	522
	523	END SUBROUTINE lim_thd_lam
	524
	525	SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
	526	!!-----------------------------------------------------------------------
	527	!! * ROUTINE lim_thd_1d2d *
	528	!!
	529	!! ** Purpose : move arrays from 1d to 2d and the reverse
	530	!!-----------------------------------------------------------------------
	531	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D
	532	! 2= from 1D to 2D
	533	INTEGER, INTENT(in) :: nbpb ! size of 1D arrays
	534	INTEGER, INTENT(in) :: jl ! ice cat
	535	INTEGER :: jk ! dummy loop indices
	536
	537	SELECT CASE( kn )
	538
	539	CASE( 1 )
	540
	541	CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
	542	CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	543	CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	544	CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	545
	546	CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	547	CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	548	DO jk = 1, nlay_s
	549	CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	550	CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	551	END DO
	552	DO jk = 1, nlay_i
	553	CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	554	CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	555	CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	556	END DO
	557
	558	CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) )
	559	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	560	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
	561	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
	562	CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	563	CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	564	IF( .NOT. lk_cpl ) THEN
	565	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	566	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	567	ENDIF
	568	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	569	CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
	570	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
	571	CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) )
	572	CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) )
	573	CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) )
	574
	575	CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) )
	576	CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) )
	577
	578	CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) )
	579	CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) )
	580	CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) )
	581	CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) )
	582	CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) )
	583	CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) )
	584
	585	CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) )
	586	CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) )
	587	CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) )
	588	CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) )
	589	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
	590	CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) )
	591
	592	CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) )
	593	CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) )
	594	CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) )
	595	CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) )
	596	CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) )
	597	CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) )
	598	CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) )
	599	CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) )
	600	CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) )
	601	CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) )
	602	CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) )
	603	CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
	604
	605	CASE( 2 )
	606
	607	CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj )
	608	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj )
	609	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj )
	610	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj )
	611	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj )
	612	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj )
	613	DO jk = 1, nlay_s
	614	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj)
	615	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj)
	616	END DO
	617	DO jk = 1, nlay_i
	618	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj)
	619	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj)
	620	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj)
	621	END DO
	622	CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj )
	623
	624	CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj )
	625	CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj )
	626
	627	CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj )
	628	CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj )
	629	CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj )
	630	CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj )
	631	CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj )
	632	CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj )
	633
	634	CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj )
	635	CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj )
	636	CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj )
	637	CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj )
	638	CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj )
	639	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
	640
	641	CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj )
	642	CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj )
	643	CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj )
	644	CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj )
	645	CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj )
	646	CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj )
	647	CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj )
	648	CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj )
	649	CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj )
	650	CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj )
	651	CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj )
	652	CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
	653	!
	654	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
	655	CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
	656
	657	END SELECT
	658
	659	END SUBROUTINE lim_thd_1d2d
	660
	661
[825]	662	SUBROUTINE lim_thd_init
	663	!!-----------------------------------------------------------------------
	664	!! * ROUTINE lim_thd_init *
	665	!!
	666	!! ** Purpose : Physical constants and parameters linked to the ice
[1572]	667	!! thermodynamics
[825]	668	!!
	669	!! ** Method : Read the namicethd namelist and check the ice-thermo
[1572]	670	!! parameter values called at the first timestep (nit000)
[825]	671	!!
	672	!! ** input : Namelist namicether
[2528]	673	!!-------------------------------------------------------------------
[4147]	674	INTEGER :: ios ! Local integer output status for namelist read
[5123]	675	NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, &
	676	& rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
	677	& nn_monocat
[825]	678	!!-------------------------------------------------------------------
[2528]	679	!
[1572]	680	IF(lwp) THEN
	681	WRITE(numout,*)
	682	WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
	683	WRITE(numout,*) '~~~~~~~'
	684	ENDIF
[2528]	685	!
[4147]	686	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
	687	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
	688	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
	689
	690	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
	691	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
	692	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
[4624]	693	IF(lwm) WRITE ( numoni, namicethd )
[5123]	694	!
	695	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
	696	nn_monocat = 0
	697	IF(lwp) WRITE(numout, *) ' nn_monocat must be 0 in multi-category case '
	698	ENDIF
[4990]	699
	700	IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
[2528]	701	!
[1572]	702	IF(lwp) THEN ! control print
[825]	703	WRITE(numout,*)
[1572]	704	WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation '
[5123]	705	WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice
	706	WRITE(numout,*)' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
	707	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb
	708	WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb
	709	WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb
	710	WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin
[1572]	711	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
	712	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
[5123]	713	WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas
	714	WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i
	715	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif
	716	WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif
	717	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon
[4688]	718	WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i
[5123]	719	WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat
[825]	720	ENDIF
[1572]	721	!
[825]	722	END SUBROUTINE lim_thd_init
	723
	724	#else
[1572]	725	!!----------------------------------------------------------------------
[2528]	726	!! Default option Dummy module NO LIM3 sea-ice model
[1572]	727	!!----------------------------------------------------------------------
[825]	728	#endif
	729
	730	!!======================================================================
	731	END MODULE limthd

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

Download in other formats: