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

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 8239

Last change on this file since 8239 was 8239, checked in by clem, 7 years ago
merge with v3_6_CMIP6_ice_diagnostics@r8238
Property svn:keywords set to `Id`
File size: 35.9 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
[7646]	24	USE ice ! sea-ice variables
[3625]	25	USE sbc_oce ! Surface boundary condition: ocean fields
	26	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	27	USE thd_ice ! thermodynamic sea-ice variables
	28	USE limthd_dif ! vertical diffusion
	29	USE limthd_dh ! ice-snow growth and melt
	30	USE limthd_da ! lateral melting
	31	USE limthd_sal ! ice salinity
	32	USE limthd_ent ! ice enthalpy redistribution
	33	USE limthd_lac ! lateral accretion
	34	USE limitd_th ! remapping thickness distribution
	35	USE limtab ! 1D <==> 2D transformation
	36	USE limvar !
	37	USE limcons ! conservation tests
	38	USE limctl ! control print
[6140]	39	!
	40	USE in_out_manager ! I/O manager
[3625]	41	USE lbclnk ! lateral boundary condition - MPP links
	42	USE lib_mpp ! MPP library
	43	USE wrk_nemo ! work arrays
	44	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4161]	45	USE timing ! Timing
[825]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
[5407]	50	PUBLIC lim_thd ! called by limstp module
	51	PUBLIC lim_thd_init ! called by sbc_lim_init
[825]	52
	53	!! * Substitutions
	54	# include "vectopt_loop_substitute.h90"
	55	!!----------------------------------------------------------------------
[2528]	56	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
[1156]	57	!! $Id$
[2528]	58	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	59	!!----------------------------------------------------------------------
	60	CONTAINS
	61
[921]	62	SUBROUTINE lim_thd( kt )
[825]	63	!!-------------------------------------------------------------------
	64	!! * ROUTINE lim_thd *
	65	!!
[4990]	66	!! ** Purpose : This routine manages ice thermodynamics
[825]	67	!!
	68	!! ** Action : - Initialisation of some variables
	69	!! - Some preliminary computation (oceanic heat flux
	70	!! at the ice base, snow acc.,heat budget of the leads)
	71	!! - selection of the icy points and put them in an array
[4990]	72	!! - call lim_thd_dif for vertical heat diffusion
	73	!! - call lim_thd_dh for vertical ice growth and melt
	74	!! - call lim_thd_ent for enthalpy remapping
	75	!! - call lim_thd_sal for ice desalination
	76	!! - call lim_thd_temp to retrieve temperature from ice enthalpy
[825]	77	!! - back to the geographic grid
	78	!!
[4990]	79	!! ** References :
[1572]	80	!!---------------------------------------------------------------------
[5123]	81	INTEGER, INTENT(in) :: kt ! number of iteration
[6140]	82	!
[4688]	83	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[5123]	84	INTEGER :: nbpb ! nb of icy pts for vertical thermo calculations
[8239]	85	REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
[4688]	86	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[5202]	87	REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
	88	REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
[7646]	89	REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
	90	!
[825]	91	!!-------------------------------------------------------------------
[4990]	92
[6140]	93	IF( nn_timing == 1 ) CALL timing_start('limthd')
[2715]	94
[7646]	95	CALL wrk_alloc( jpi,jpj, zu_io, zv_io, zfric )
	96
	97	IF( kt == nit000 .AND. lwp ) THEN
	98	WRITE(numout,*)''
	99	WRITE(numout,*)' lim_thd '
	100	WRITE(numout,*)' ~~~~~~~~'
	101	ENDIF
	102
[4688]	103	! conservation test
[7646]	104	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	105
[5167]	106	CALL lim_var_glo2eqv
[7646]	107
	108	!---------------------------------------------!
	109	! computation of friction velocity at T points
	110	!---------------------------------------------!
	111	IF( ln_limdyn ) THEN
[7753]	112	zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
	113	zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
[7646]	114	DO jj = 2, jpjm1
	115	DO ji = fs_2, fs_jpim1
	116	zfric(ji,jj) = rn_cio * ( 0.5_wp * &
	117	& ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
	118	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
	119	END DO
	120	END DO
	121	ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean
	122	DO jj = 2, jpjm1
	123	DO ji = fs_2, fs_jpim1
	124	zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp * &
	125	& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
	126	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
	127	END DO
	128	END DO
	129	ENDIF
	130	CALL lbc_lnk( zfric, 'T', 1. )
	131	!
	132	!----------------------------------!
	133	! Initialization and units change
	134	!----------------------------------!
[7753]	135	ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice
[825]	136
[5123]	137	! Change the units of heat content; from J/m2 to J/m3
[825]	138	DO jl = 1, jpl
[921]	139	DO jk = 1, nlay_i
	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
[5134]	142	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) )
[4688]	143	!Energy of melting q(S,T) [J.m-3]
[5123]	144	e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
[921]	145	END DO
[825]	146	END DO
[921]	147	END DO
	148	DO jk = 1, nlay_s
	149	DO jj = 1, jpj
	150	DO ji = 1, jpi
[5134]	151	rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )
[4688]	152	!Energy of melting q(S,T) [J.m-3]
[5123]	153	e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
[921]	154	END DO
[825]	155	END DO
[921]	156	END DO
[825]	157	END DO
	158
[7646]	159	!--------------------------------------------------------------------!
	160	! Partial computation of forcing for the thermodynamic sea ice model
	161	!--------------------------------------------------------------------!
[921]	162	DO jj = 1, jpj
	163	DO ji = 1, jpi
[5134]	164	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
[2528]	165	!
[921]	166	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	167	! ! practically no "direct lateral ablation"
	168	!
	169	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	170	! ! temperature and turbulent mixing (McPhee, 1992)
[4688]	171	!
	172	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
[5407]	173	zqld = tmask(ji,jj,1) * rdt_ice * &
	174	& ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
[825]	175
[5123]	176	! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
[8239]	177	! includes supercooling potential energy (>0) or "above-freezing" energy (<0)
[6140]	178	zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
[4688]	179
[8239]	180	! --- Above-freezing sensible heat content (J/m2 grid)
	181	zqfr_neg = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp )
	182
	183	! --- Sensible ocean-to-ice heat flux (W/m2)
[7646]	184	zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
[8239]	185	fhtur(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
	186
	187	fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
[5123]	188	! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
	189	! the freezing point, so that we do not have SST < T_freeze
	190	! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
	191
[8239]	192	!-- Energy Budget of the leads (J.m-2), source of lateral accretion. Must be < 0 to form ice
[5123]	193	qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
[4688]	194
	195	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
[5123]	196	IF( zqld > 0._wp ) THEN
	197	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]	198	qlead(ji,jj) = 0._wp
[4990]	199	ELSE
	200	fhld (ji,jj) = 0._wp
[4688]	201	ENDIF
[2528]	202	!
[7646]	203	! Net heat flux on top of the ice-ocean [W.m-2]
	204	! ---------------------------------------------
[5407]	205	hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
[7646]	206	END DO
	207	END DO
	208
	209	! In case we bypass open-water ice formation
[7753]	210	IF( .NOT. ln_limdO ) qlead(:,:) = 0._wp
[7646]	211	! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
[7753]	212	IF( .NOT. ln_limdH ) hfx_in(:,:) = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
	213	IF( .NOT. ln_limdH ) fhtur (:,:) = 0._wp ; fhld (:,:) = 0._wp
[4688]	214
[7646]	215	! ---------------------------------------------------------------------
	216	! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
	217	! ---------------------------------------------------------------------
	218	! First step here : non solar + precip - qlead - qturb
	219	! Second step in limthd_dh : heat remaining if total melt (zq_rema)
	220	! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar
	221	DO jj = 1, jpj
	222	DO ji = 1, jpi
[5407]	223	hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean
	224	& - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation
	225	& - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence
	226	& - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt
	227	! (fhld should be 0 while bott growth)
[825]	228	END DO
	229	END DO
	230
[921]	231	!------------------------------------------------------------------------------!
[7646]	232	! Thermodynamic computation (only on grid points covered by ice)
[921]	233	!------------------------------------------------------------------------------!
[825]	234
	235	DO jl = 1, jpl !loop over ice categories
	236
[7646]	237	! select ice covered grid points
[825]	238	nbpb = 0
	239	DO jj = 1, jpj
	240	DO ji = 1, jpi
[5123]	241	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
[825]	242	nbpb = nbpb + 1
	243	npb(nbpb) = (jj - 1) * jpi + ji
	244	ENDIF
	245	END DO
	246	END DO
	247
[4333]	248	! debug point to follow
	249	jiindex_1d = 0
[7646]	250	IF( ln_limctl ) THEN
[5128]	251	DO ji = mi0(iiceprt), mi1(iiceprt)
	252	DO jj = mj0(jiceprt), mj1(jiceprt)
[4333]	253	jiindex_1d = (jj - 1) * jpi + ji
[4688]	254	WRITE(numout,*) ' lim_thd : Category no : ', jl
[4333]	255	END DO
	256	END DO
	257	ENDIF
	258
[7646]	259	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
[825]	260
[7646]	261	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
	262	!
	263	s_i_new (:) = 0._wp ; dh_s_tot (:) = 0._wp ! --- some init --- !
	264	dh_i_surf (:) = 0._wp ; dh_i_bott(:) = 0._wp
	265	dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp
[869]	266
[7646]	267	CALL lim_thd_1d2d( nbpb, jl, 1 ) ! --- Move to 1D arrays --- !
[6140]	268	!
[7646]	269	IF( ln_limdH ) CALL lim_thd_dif( 1, nbpb ) ! --- Ice/Snow Temperature profile --- !
[6140]	270	!
[7646]	271	IF( ln_limdH ) CALL lim_thd_dh( 1, nbpb ) ! --- Ice/Snow thickness --- !
[6140]	272	!
[7646]	273	IF( ln_limdH ) CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) ! --- Ice enthalpy remapping --- !
[6140]	274	!
[7646]	275	CALL lim_thd_sal( 1, nbpb ) ! --- Ice salinity --- !
[6140]	276	!
[7646]	277	CALL lim_thd_temp( 1, nbpb ) ! --- temperature update --- !
[6140]	278	!
[7646]	279	IF( ln_limdH ) THEN
	280	IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
	281	CALL lim_thd_lam( 1, nbpb ) ! --- extra lateral melting if monocat --- !
	282	END IF
[5123]	283	END IF
[2528]	284	!
[7646]	285	CALL lim_thd_1d2d( nbpb, jl, 2 ) ! --- Move to 2D arrays --- !
[6140]	286	!
[7646]	287	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
[1572]	288	ENDIF
	289	!
[5167]	290	END DO !jl
[825]	291
[7646]	292	IF( ln_limdA) CALL lim_thd_da ! --- lateral melting --- !
[825]	293
[5123]	294	! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
[825]	295	DO jl = 1, jpl
[921]	296	DO jk = 1, nlay_i
[7753]	297	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
[1572]	298	END DO
[825]	299	DO jk = 1, nlay_s
[7753]	300	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
[1572]	301	END DO
	302	END DO
[7753]	303
	304	! Change thickness to volume
	305	v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
	306	v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
	307	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
[825]	308
[5202]	309	! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat)
	310	DO jl = 1, jpl
	311	DO jj = 1, jpj
	312	DO ji = 1, jpi
	313	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) )
	314	oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 )
	315	END DO
	316	END DO
	317	END DO
	318
[5134]	319	CALL lim_var_zapsmall
[5202]	320
[7646]	321	! control checks
	322	IF( ln_limctl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' ) ! control print
[2528]	323	!
[7646]	324	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4990]	325
[7646]	326	!------------------------------------------------!
	327	! Transport ice between thickness categories
	328	!------------------------------------------------!
[5167]	329	! Given thermodynamic growth rates, transport ice between thickness categories.
[7646]	330	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[5123]	331
[5167]	332	IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )
[5123]	333
[7646]	334	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[5167]	335
[7646]	336	!------------------------------------------------!
	337	! Add frazil ice growing in leads
	338	!------------------------------------------------!
	339	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[5167]	340
[7646]	341	IF( ln_limdO ) CALL lim_thd_lac
[5123]	342
[7646]	343	! conservation test
	344	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[5123]	345
[5167]	346	! Control print
[7646]	347	IF( ln_ctl ) CALL lim_prt3D( 'limthd' )
[4161]	348	!
[7646]	349	CALL wrk_dealloc( jpi,jpj, zu_io, zv_io, zfric )
[6140]	350	!
[7646]	351	IF( nn_timing == 1 ) CALL timing_stop('limthd')
	352
[4688]	353	END SUBROUTINE lim_thd
[825]	354
[5407]	355
[4688]	356	SUBROUTINE lim_thd_temp( kideb, kiut )
[825]	357	!!-----------------------------------------------------------------------
[4688]	358	!! * ROUTINE lim_thd_temp *
[825]	359	!!
[4688]	360	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
[825]	361	!!
	362	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	363	!!-------------------------------------------------------------------
[1572]	364	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
[6140]	365	!
[2715]	366	INTEGER :: ji, jk ! dummy loop indices
[4990]	367	REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar
[825]	368	!!-------------------------------------------------------------------
[4688]	369	! Recover ice temperature
	370	DO jk = 1, nlay_i
[825]	371	DO ji = kideb, kiut
[5123]	372	ztmelts = -tmut * s_i_1d(ji,jk) + rt0
[4688]	373	! Conversion q(S,T) -> T (second order equation)
	374	zaaa = cpic
[5123]	375	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
	376	zccc = lfus * ( ztmelts - rt0 )
[4688]	377	zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
[5123]	378	t_i_1d(ji,jk) = rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
[4688]	379
	380	! mask temperature
[4990]	381	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )
[5123]	382	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
[4688]	383	END DO
	384	END DO
[6140]	385	!
[4688]	386	END SUBROUTINE lim_thd_temp
[825]	387
[6140]	388
[5123]	389	SUBROUTINE lim_thd_lam( kideb, kiut )
	390	!!-----------------------------------------------------------------------
	391	!! * ROUTINE lim_thd_lam *
	392	!!
	393	!! ** Purpose : Lateral melting in case monocategory
	394	!! ( dA = A/2h dh )
	395	!!-----------------------------------------------------------------------
	396	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
[6140]	397	!
	398	INTEGER :: ji ! dummy loop indices
	399	REAL(wp) :: zhi_bef ! ice thickness before thermo
	400	REAL(wp) :: zdh_mel, zda_mel ! net melting
	401	REAL(wp) :: zvi, zvs ! ice/snow volumes
	402	!!-----------------------------------------------------------------------
	403	!
[5123]	404	DO ji = kideb, kiut
[6416]	405	zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
[5167]	406	IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN
	407	zvi = a_i_1d(ji) * ht_i_1d(ji)
	408	zvs = a_i_1d(ji) * ht_s_1d(ji)
[5123]	409	! lateral melting = concentration change
	410	zhi_bef = ht_i_1d(ji) - zdh_mel
[5167]	411	rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
	412	zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
	413	a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel )
[7646]	414	! adjust thickness
[5167]	415	ht_i_1d(ji) = zvi / a_i_1d(ji)
	416	ht_s_1d(ji) = zvs / a_i_1d(ji)
[5123]	417	! retrieve total concentration
	418	at_i_1d(ji) = a_i_1d(ji)
	419	END IF
	420	END DO
[6140]	421	!
[5123]	422	END SUBROUTINE lim_thd_lam
	423
[6140]	424
[5123]	425	SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
	426	!!-----------------------------------------------------------------------
	427	!! * ROUTINE lim_thd_1d2d *
	428	!!
	429	!! ** Purpose : move arrays from 1d to 2d and the reverse
	430	!!-----------------------------------------------------------------------
[6140]	431	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
[5123]	432	INTEGER, INTENT(in) :: nbpb ! size of 1D arrays
	433	INTEGER, INTENT(in) :: jl ! ice cat
[6140]	434	!
[5123]	435	INTEGER :: jk ! dummy loop indices
[6140]	436	!!-----------------------------------------------------------------------
	437	!
[5123]	438	SELECT CASE( kn )
[6140]	439	!
	440	CASE( 1 ) ! from 2D to 1D
	441	!
[5123]	442	CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
	443	CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	444	CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	445	CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[6140]	446	!
[5123]	447	CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	448	CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	449	DO jk = 1, nlay_s
	450	CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	451	CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	452	END DO
	453	DO jk = 1, nlay_i
	454	CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	455	CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	456	CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	457	END DO
[6140]	458	!
[5407]	459	CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
[6416]	460	CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[5123]	461	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	462	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
	463	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
	464	CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	465	CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[5407]	466	CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
[5123]	467	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	468	CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
	469	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
	470	CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) )
	471	CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) )
	472	CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) )
[6140]	473	!
[5123]	474	CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) )
[8233]	475	CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum , jpi, jpj, npb(1:nbpb) )
[5123]	476	CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) )
[8239]	477	CALL tab_2d_1d( nbpb, wfx_snw_sub_1d(1:nbpb), wfx_snw_sub , jpi, jpj, npb(1:nbpb) )
	478	CALL tab_2d_1d( nbpb, wfx_ice_sub_1d(1:nbpb), wfx_ice_sub , jpi, jpj, npb(1:nbpb) )
[6140]	479	!
[5123]	480	CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) )
	481	CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) )
	482	CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) )
	483	CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) )
	484	CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) )
	485	CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) )
[6140]	486	!
[5123]	487	CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) )
	488	CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) )
	489	CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) )
	490	CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) )
	491	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
	492	CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) )
[6416]	493	CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub , jpi, jpj,npb(1:nbpb) )
[6140]	494	!
[5123]	495	CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) )
	496	CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) )
	497	CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) )
	498	CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) )
	499	CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) )
	500	CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) )
	501	CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) )
	502	CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) )
	503	CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) )
	504	CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) )
	505	CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) )
[5146]	506	CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
[5123]	507	CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
[6140]	508	!
[8239]	509	! SIMIP diagnostics
	510	CALL tab_2d_1d( nbpb, diag_fc_bo_1d (1:nbpb), diag_fc_bo , jpi, jpj, npb(1:nbpb) )
	511	CALL tab_2d_1d( nbpb, diag_fc_su_1d (1:nbpb), diag_fc_su , jpi, jpj, npb(1:nbpb) )
	512	!
[6140]	513	CASE( 2 ) ! from 1D to 2D
	514	!
[5123]	515	CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj )
	516	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj )
	517	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj )
	518	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj )
	519	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj )
	520	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj )
	521	DO jk = 1, nlay_s
	522	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj)
	523	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj)
	524	END DO
	525	DO jk = 1, nlay_i
	526	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj)
	527	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj)
	528	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj)
	529	END DO
	530	CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj )
[6140]	531	!
[5123]	532	CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj )
[8233]	533	CALL tab_1d_2d( nbpb, wfx_snw_sum , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
[5123]	534	CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj )
[8239]	535	CALL tab_1d_2d( nbpb, wfx_snw_sub , npb, wfx_snw_sub_1d(1:nbpb), jpi, jpj )
	536	CALL tab_1d_2d( nbpb, wfx_ice_sub , npb, wfx_ice_sub_1d(1:nbpb), jpi, jpj )
[6140]	537	!
[5123]	538	CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj )
	539	CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj )
	540	CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj )
	541	CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj )
	542	CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj )
	543	CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj )
[6140]	544	!
[5123]	545	CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj )
	546	CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj )
	547	CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj )
	548	CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj )
	549	CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj )
	550	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
[6416]	551	CALL tab_1d_2d( nbpb, sfx_sub , npb, sfx_sub_1d(1:nbpb) , jpi, jpj )
[6140]	552	!
[5123]	553	CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj )
	554	CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj )
	555	CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj )
	556	CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj )
	557	CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj )
	558	CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj )
	559	CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj )
	560	CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj )
	561	CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj )
	562	CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj )
	563	CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj )
	564	CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
[5146]	565	CALL tab_1d_2d( nbpb, hfx_err_dif , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
[5123]	566	!
	567	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
	568	CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
[8239]	569	!
	570	! SIMIP diagnostics
	571	CALL tab_1d_2d( nbpb, t_si(:,:,jl) , npb, t_si_1d (1:nbpb) , jpi, jpj )
	572	CALL tab_1d_2d( nbpb, diag_fc_bo , npb, diag_fc_bo_1d(1:nbpb) , jpi, jpj )
	573	CALL tab_1d_2d( nbpb, diag_fc_su , npb, diag_fc_su_1d(1:nbpb) , jpi, jpj )
[5123]	574	END SELECT
[6140]	575	!
[5123]	576	END SUBROUTINE lim_thd_1d2d
	577
	578
[825]	579	SUBROUTINE lim_thd_init
	580	!!-----------------------------------------------------------------------
	581	!! * ROUTINE lim_thd_init *
	582	!!
	583	!! ** Purpose : Physical constants and parameters linked to the ice
[1572]	584	!! thermodynamics
[825]	585	!!
	586	!! ** Method : Read the namicethd namelist and check the ice-thermo
[1572]	587	!! parameter values called at the first timestep (nit000)
[825]	588	!!
	589	!! ** input : Namelist namicether
[2528]	590	!!-------------------------------------------------------------------
[4147]	591	INTEGER :: ios ! Local integer output status for namelist read
[7813]	592	NAMELIST/namicethd/ rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, ln_it_qnsice, nn_monocat, rn_cdsn, &
	593	& ln_limdH, rn_betas, &
	594	& ln_limdA, rn_beta, rn_dmin, &
[7646]	595	& ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin
[825]	596	!!-------------------------------------------------------------------
[2528]	597	!
[4147]	598	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
	599	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
	600	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
	601
	602	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
	603	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
	604	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
[4624]	605	IF(lwm) WRITE ( numoni, namicethd )
[5123]	606	!
[7813]	607	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
	608	nn_monocat = 0
	609	IF(lwp) WRITE(numout,*)
	610	IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
	611	ENDIF
	612	!
[1572]	613	IF(lwp) THEN ! control print
[7646]	614	WRITE(numout,*) 'lim_thd_init : Ice Thermodynamics'
	615	WRITE(numout,*) '~~~~~~~~~~~~~'
	616	WRITE(numout,*)' -- limthd_dif --'
	617	WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i
	618	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif
	619	WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif
	620	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon
	621	WRITE(numout,*)' iterate the surface non-solar flux (T) or not (F) ln_it_qnsice = ', ln_it_qnsice
	622	WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat
[7813]	623	WRITE(numout,*)' thermal conductivity of the snow rn_cdsn = ', rn_cdsn
[7646]	624	WRITE(numout,*)' -- limthd_dh --'
	625	WRITE(numout,*)' activate ice thick change from top/bot (T) or not (F) ln_limdH = ', ln_limdH
	626	WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas
	627	WRITE(numout,*)' -- limthd_da --'
	628	WRITE(numout,*)' activate lateral melting (T) or not (F) ln_limdA = ', ln_limdA
	629	WRITE(numout,*)' Coef. beta for lateral melting param. rn_beta = ', rn_beta
	630	WRITE(numout,*)' Minimum floe diameter for lateral melting param. rn_dmin = ', rn_dmin
	631	WRITE(numout,*)' -- limthd_lac --'
	632	WRITE(numout,*)' activate ice growth in open-water (T) or not (F) ln_limdO = ', ln_limdO
[5123]	633	WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice
	634	WRITE(numout,*)' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
	635	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb
	636	WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb
	637	WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb
[7646]	638	WRITE(numout,*)' -- limitd_th --'
[5123]	639	WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin
[4688]	640	WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i
[825]	641	ENDIF
[1572]	642	!
[825]	643	END SUBROUTINE lim_thd_init
	644
	645	#else
[1572]	646	!!----------------------------------------------------------------------
[2528]	647	!! Default option Dummy module NO LIM3 sea-ice model
[1572]	648	!!----------------------------------------------------------------------
[825]	649	#endif
	650
	651	!!======================================================================
	652	END MODULE limthd

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

Download in other formats: