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

Context Navigation

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago
[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]
Property svn:keywords set to `Id`
File size: 46.1 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
[3625]	10	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_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
[4045]	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
[4220]	24	USE oce , ONLY : iatte, oatte
[3625]	25	USE ice ! LIM: sea-ice variables
	26	USE par_ice ! LIM: sea-ice parameters
	27	USE sbc_oce ! Surface boundary condition: ocean fields
	28	USE sbc_ice ! Surface boundary condition: ice fields
	29	USE thd_ice ! LIM thermodynamic sea-ice variables
	30	USE dom_ice ! LIM sea-ice domain
	31	USE domvvl ! domain: variable volume level
	32	USE limthd_dif ! LIM: thermodynamics, vertical diffusion
	33	USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation
	34	USE limthd_sal ! LIM: thermodynamics, ice salinity
	35	USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution
	36	USE limtab ! LIM: 1D <==> 2D transformation
	37	USE limvar ! LIM: sea-ice variables
	38	USE lbclnk ! lateral boundary condition - MPP links
	39	USE lib_mpp ! MPP library
	40	USE wrk_nemo ! work arrays
	41	USE in_out_manager ! I/O manager
	42	USE prtctl ! Print control
	43	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4072]	44	USE timing ! Timing
[825]	45
	46	IMPLICIT NONE
	47	PRIVATE
	48
[2528]	49	PUBLIC lim_thd ! called by limstp module
	50	PUBLIC lim_thd_init ! called by iceini module
[825]	51
[4332]	52	REAL(wp) :: epsi10 = 1.e-10_wp !
[2528]	53	REAL(wp) :: zzero = 0._wp !
	54	REAL(wp) :: zone = 1._wp !
[825]	55
	56	!! * Substitutions
	57	# include "domzgr_substitute.h90"
	58	# include "vectopt_loop_substitute.h90"
	59	!!----------------------------------------------------------------------
[2528]	60	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
[1156]	61	!! $Id$
[2528]	62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	63	!!----------------------------------------------------------------------
	64	CONTAINS
	65
[921]	66	SUBROUTINE lim_thd( kt )
[825]	67	!!-------------------------------------------------------------------
	68	!! * ROUTINE lim_thd *
	69	!!
	70	!! ** Purpose : This routine manages the ice thermodynamic.
	71	!!
	72	!! ** Action : - Initialisation of some variables
	73	!! - Some preliminary computation (oceanic heat flux
	74	!! at the ice base, snow acc.,heat budget of the leads)
	75	!! - selection of the icy points and put them in an array
	76	!! - call lim_vert_ther for vert ice thermodynamic
	77	!! - back to the geographic grid
	78	!! - selection of points for lateral accretion
	79	!! - call lim_lat_acc for the ice accretion
	80	!! - back to the geographic grid
	81	!!
[1572]	82	!! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
	83	!!---------------------------------------------------------------------
	84	INTEGER, INTENT(in) :: kt ! number of iteration
[825]	85	!!
[1572]	86	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	87	INTEGER :: nbpb ! nb of icy pts for thermo. cal.
[2715]	88	REAL(wp) :: zfric_umin = 5e-03_wp ! lower bound for the friction velocity
	89	REAL(wp) :: zfric_umax = 2e-02_wp ! upper bound for the friction velocity
	90	REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! local scalar
	91	REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - -
[3294]	92	REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif
[4045]	93	REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
	94	REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
[825]	95	!!-------------------------------------------------------------------
[4072]	96	IF( nn_timing == 1 ) CALL timing_start('limthd')
[2715]	97
[3294]	98	CALL wrk_alloc( jpi, jpj, zqlbsbq )
[2715]	99
[4045]	100	! -------------------------------
	101	!- check conservation (C Rousset)
	102	IF (ln_limdiahsb) THEN
	103	zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
	104	zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
	105	zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
	106	zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
	107	ENDIF
	108	!- check conservation (C Rousset)
	109	! -------------------------------
	110
[921]	111	!------------------------------------------------------------------------------!
	112	! 1) Initialization of diagnostic variables !
	113	!------------------------------------------------------------------------------!
[825]	114
	115	!--------------------
	116	! 1.2) Heat content
	117	!--------------------
[1572]	118	! Change the units of heat content; from global units to J.m3
[825]	119	DO jl = 1, jpl
[921]	120	DO jk = 1, nlay_i
	121	DO jj = 1, jpj
	122	DO ji = 1, jpi
	123	!Energy of melting q(S,T) [J.m-3]
[4332]	124	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )
[921]	125	!0 if no ice and 1 if yes
[4332]	126	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )
[921]	127	!convert units ! very important that this line is here
	128	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
	129	END DO
[825]	130	END DO
[921]	131	END DO
	132	DO jk = 1, nlay_s
	133	DO jj = 1, jpj
	134	DO ji = 1, jpi
	135	!Energy of melting q(S,T) [J.m-3]
[4332]	136	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )
[921]	137	!0 if no ice and 1 if yes
[4332]	138	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) )
[921]	139	!convert units ! very important that this line is here
	140	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb
	141	END DO
[825]	142	END DO
[921]	143	END DO
[825]	144	END DO
	145
	146	!-----------------------------------
	147	! 1.4) Compute global heat content
	148	!-----------------------------------
[1572]	149	qt_i_in (:,:) = 0.e0
	150	qt_s_in (:,:) = 0.e0
	151	qt_i_fin (:,:) = 0.e0
	152	qt_s_fin (:,:) = 0.e0
[869]	153	sum_fluxq(:,:) = 0.e0
[1572]	154	fatm (:,:) = 0.e0
[825]	155
[921]	156	! 2) Partial computation of forcing for the thermodynamic sea ice model. !
	157	!-----------------------------------------------------------------------------!
[825]	158
[921]	159	!CDIR NOVERRCHK
	160	DO jj = 1, jpj
	161	!CDIR NOVERRCHK
	162	DO ji = 1, jpi
[4332]	163	zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) )
[2528]	164	!
[921]	165	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	166	! ! practically no "direct lateral ablation"
	167	!
	168	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	169	! ! temperature and turbulent mixing (McPhee, 1992)
[825]	170	! friction velocity
	171	zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )
	172
	173	! here the drag will depend on ice thickness and type (0.006)
[4332]	174	fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
[825]	175	! also category dependent
[921]	176	! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead
[4332]	177	qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice
[921]	178	!
[2528]	179	! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)
	180	! ! caution: exponent betas used as more snow can fallinto leads
	181	qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( &
[4045]	182	& pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif
[2528]	183	& + qns(ji,jj) & ! non solar heat
	184	& + fdtcn(ji,jj) & ! turbulent ice-ocean heat
[4045]	185	& + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step
[4506]	186	! MV HC 2014 check that this is good.
	187	! We should remove the heat content of precip that has fallen on sea ice
[2387]	188	& - pfrld(ji,jj)*betas sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting
[4506]	189	! MV HC 2014 partie heat content manque
[2528]	190	!
[825]	191	! Positive heat budget is used for bottom ablation
	192	zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) )
	193	!= 1 if positive heat budget
	194	zpareff = 1.0 - zinda * zfntlat
[1572]	195	!= 0 if ice and positive heat budget and 1 if one of those two is false
[4332]	196	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) )
[2528]	197	!
[825]	198	! Heat budget of the lead, energy transferred from ice to ocean
	199	qldif (ji,jj) = zpareff * qldif(ji,jj)
	200	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
[2528]	201	!
[1572]	202	! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
[4332]	203	qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
[2528]	204	!
[1572]	205	! oceanic heat flux (limthd_dh)
[4332]	206	fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) )
[1571]	207	!
[825]	208	END DO
	209	END DO
	210
[921]	211	!------------------------------------------------------------------------------!
	212	! 3) Select icy points and fulfill arrays for the vectorial grid.
	213	!------------------------------------------------------------------------------!
[825]	214
	215	DO jl = 1, jpl !loop over ice categories
	216
[921]	217	IF( kt == nit000 .AND. lwp ) THEN
	218	WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl
	219	WRITE(numout,*) ' ~~~~~~~~'
	220	ENDIF
[825]	221
[4332]	222	zareamin = epsi10
[825]	223	nbpb = 0
	224	DO jj = 1, jpj
	225	DO ji = 1, jpi
	226	IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN
	227	nbpb = nbpb + 1
	228	npb(nbpb) = (jj - 1) * jpi + ji
	229	ENDIF
	230	END DO
	231	END DO
	232
[4332]	233	! debug point to follow
	234	jiindex_1d = 0
	235	IF( ln_nicep ) THEN
	236	DO ji = mi0(jiindx), mi1(jiindx)
	237	DO jj = mj0(jjindx), mj1(jjindx)
	238	jiindex_1d = (jj - 1) * jpi + ji
	239	END DO
	240	END DO
	241	ENDIF
	242
[921]	243	!------------------------------------------------------------------------------!
	244	! 4) Thermodynamic computation
	245	!------------------------------------------------------------------------------!
[825]	246
[2715]	247	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
[869]	248
[1572]	249	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
[825]	250
[921]	251	!-------------------------
	252	! 4.1 Move to 1D arrays
	253	!-------------------------
[825]	254
[1572]	255	CALL tab_2d_1d( nbpb, at_i_b (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
	256	CALL tab_2d_1d( nbpb, a_i_b (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	257	CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	258	CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	259
[1572]	260	CALL tab_2d_1d( nbpb, t_su_b (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	261	CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	262	DO jk = 1, nlay_s
[1572]	263	CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	264	CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	265	END DO
	266	DO jk = 1, nlay_i
[1572]	267	CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	268	CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	269	CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	270	END DO
	271
[1572]	272	CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) )
	273	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	274	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
	275	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
	276	CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	277	#if ! defined key_coupled
[3625]	278	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	279	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
[825]	280	#endif
[3625]	281	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	282	CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
	283	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
	284	CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb), fbif , jpi, jpj, npb(1:nbpb) )
	285	CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb), qldif , jpi, jpj, npb(1:nbpb) )
	286	CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice , jpi, jpj, npb(1:nbpb) )
[4506]	287	CALL tab_2d_1d( nbpb, rdq_ice_1d (1:nbpb), rdq_ice , jpi, jpj, npb(1:nbpb) )
[3625]	288	CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw , jpi, jpj, npb(1:nbpb) )
[4506]	289	CALL tab_2d_1d( nbpb, rdq_snw_1d (1:nbpb), rdq_snw , jpi, jpj, npb(1:nbpb) )
[3625]	290	CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb), dmgwi , jpi, jpj, npb(1:nbpb) )
	291	CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb), zqlbsbq , jpi, jpj, npb(1:nbpb) )
[825]	292
[3625]	293	CALL tab_2d_1d( nbpb, sfx_thd_1d (1:nbpb), sfx_thd , jpi, jpj, npb(1:nbpb) )
	294	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
	295	CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb), fhbri , jpi, jpj, npb(1:nbpb) )
	296	CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) )
	297	CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) )
[825]	298
[4045]	299	CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif
	300	CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif
[921]	301	!--------------------------------
	302	! 4.3) Thermodynamic processes
	303	!--------------------------------
	304
[4332]	305	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting
	306	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl )
[921]	307
[1572]	308	! !---------------------------------!
	309	CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile !
	310	! !---------------------------------!
[825]	311
[1572]	312	CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh
[825]	313
[4332]	314	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )
	315	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dif( 1 , nbpb , jl )
[825]	316
[1572]	317	! !---------------------------------!
	318	CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness !
	319	! !---------------------------------!
[825]	320
[1572]	321	! !---------------------------------!
	322	CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping !
	323	! !---------------------------------!
[825]	324
[1572]	325	! !---------------------------------!
	326	CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation !
	327	! !---------------------------------!
[825]	328
[921]	329	! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting
[4332]	330	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )
	331	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dh ( 1 , nbpb , jl )
[825]	332
[921]	333	!--------------------------------
	334	! 4.4) Move 1D to 2D vectors
	335	!--------------------------------
[825]	336
[3625]	337	CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb) , jpi, jpj )
	338	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_b (1:nbpb) , jpi, jpj )
	339	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_b (1:nbpb) , jpi, jpj )
	340	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_b (1:nbpb) , jpi, jpj )
	341	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_b (1:nbpb) , jpi, jpj )
	342	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_b (1:nbpb) , jpi, jpj )
[825]	343	DO jk = 1, nlay_s
[3625]	344	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b (1:nbpb,jk), jpi, jpj)
	345	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b (1:nbpb,jk), jpi, jpj)
[825]	346	END DO
	347	DO jk = 1, nlay_i
[3625]	348	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b (1:nbpb,jk), jpi, jpj)
	349	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b (1:nbpb,jk), jpi, jpj)
	350	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b (1:nbpb,jk), jpi, jpj)
[825]	351	END DO
[3625]	352	CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
	353	CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
	354	CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
	355	CALL tab_1d_2d( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj )
[4506]	356	CALL tab_1d_2d( nbpb, rdq_ice , npb, rdq_ice_1d(1:nbpb) , jpi, jpj )
[3625]	357	CALL tab_1d_2d( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj )
[4506]	358	CALL tab_1d_2d( nbpb, rdq_snw , npb, rdq_snw_1d(1:nbpb) , jpi, jpj )
[3625]	359	CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
	360	CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
	361	CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
	362	CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
	363	CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
	364	CALL tab_1d_2d( nbpb, sfx_thd , npb, sfx_thd_1d(1:nbpb) , jpi, jpj )
[2528]	365	!
[1572]	366	IF( num_sal == 2 ) THEN
[3625]	367	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
	368	CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb) , jpi, jpj )
[825]	369	ENDIF
[2528]	370	!
[3625]	371	!+++++ temporary stuff for a dummy version
[825]	372	CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )
	373	CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )
	374	CALL tab_1d_2d( nbpb, fsup2D , npb, fsup (1:nbpb) , jpi, jpj )
	375	CALL tab_1d_2d( nbpb, focea2D , npb, focea (1:nbpb) , jpi, jpj )
	376	CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )
	377	CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )
[888]	378	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
[825]	379	!+++++
[2528]	380	!
[1572]	381	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	382	ENDIF
	383	!
	384	END DO
[825]	385
[921]	386	!------------------------------------------------------------------------------!
	387	! 5) Global variables, diagnostics
	388	!------------------------------------------------------------------------------!
[825]	389
	390	!------------------------
	391	! 5.1) Ice heat content
	392	!------------------------
[3625]	393	! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules)
[2715]	394	zcoef = 1._wp / ( unit_fac * REAL( nlay_i ) )
[825]	395	DO jl = 1, jpl
[921]	396	DO jk = 1, nlay_i
[1572]	397	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef
	398	END DO
	399	END DO
[825]	400
	401	!------------------------
	402	! 5.2) Snow heat content
	403	!------------------------
[3625]	404	! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules)
[2715]	405	zcoef = 1._wp / ( unit_fac * REAL( nlay_s ) )
[825]	406	DO jl = 1, jpl
	407	DO jk = 1, nlay_s
[1572]	408	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef
	409	END DO
	410	END DO
[825]	411
	412	!----------------------------------
	413	! 5.3) Change thickness to volume
	414	!----------------------------------
	415	CALL lim_var_eqv2glo
	416
	417	!--------------------------------------------
	418	! 5.4) Diagnostic thermodynamic growth rates
	419	!--------------------------------------------
[4045]	420	!clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes
	421	!clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
[825]	422
[4332]	423	IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
[825]	424
[2528]	425	IF(ln_ctl) THEN ! Control print
[867]	426	CALL prt_ctl_info(' ')
	427	CALL prt_ctl_info(' - Cell values : ')
	428	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	429	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd : cell area :')
	430	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :')
	431	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :')
	432	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd : vt_s :')
	433	DO jl = 1, jpl
[867]	434	CALL prt_ctl_info(' ')
[863]	435	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	436	CALL prt_ctl_info(' ~~~~~~~~~~')
	437	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_thd : a_i : ')
	438	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_thd : ht_i : ')
	439	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_thd : ht_s : ')
	440	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_thd : v_i : ')
	441	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_thd : v_s : ')
	442	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_thd : e_s : ')
	443	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_thd : t_su : ')
	444	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_thd : t_snow : ')
	445	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_thd : sm_i : ')
	446	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_thd : smv_i : ')
	447	DO jk = 1, nlay_i
[867]	448	CALL prt_ctl_info(' ')
[863]	449	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	450	CALL prt_ctl_info(' ~~~~~~~')
	451	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd : t_i : ')
	452	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd : e_i : ')
	453	END DO
	454	END DO
	455	ENDIF
[2528]	456	!
[4045]	457	! -------------------------------
	458	!- check conservation (C Rousset)
	459	IF (ln_limdiahsb) THEN
	460	zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
	461	zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
	462
	463	zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
	464	zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
	465
	466	zchk_vmin = glob_min(v_i)
	467	zchk_amax = glob_max(SUM(a_i,dim=3))
	468	zchk_amin = glob_min(a_i)
	469
	470	IF(lwp) THEN
	471	IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,) 'violation volume [m3/day] (limthd) = ',(zchk_v_i rday)
	472	IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,) 'violation saline [psum3/day] (limthd) = ',(zchk_smv * rday)
	473	IF ( zchk_vmin < 0. ) WRITE(numout,) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin 1.e-3)
	474	IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax
	475	IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin
	476	ENDIF
	477	ENDIF
	478	!- check conservation (C Rousset)
	479	! -------------------------------
	480	!
[3294]	481	CALL wrk_dealloc( jpi, jpj, zqlbsbq )
[2715]	482	!
[4072]	483	IF( nn_timing == 1 ) CALL timing_stop('limthd')
[825]	484	END SUBROUTINE lim_thd
	485
	486
[1572]	487	SUBROUTINE lim_thd_glohec( eti, ets, etilayer, kideb, kiut, jl )
[825]	488	!!-----------------------------------------------------------------------
	489	!! * ROUTINE lim_thd_glohec *
	490	!!
	491	!! ** Purpose : Compute total heat content for each category
	492	!! Works with 1d vectors only
[1572]	493	!!-----------------------------------------------------------------------
	494	INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop
	495	INTEGER , INTENT(in ) :: jl ! category number
	496	REAL(wp), INTENT( out), DIMENSION (jpij,jpl ) :: eti, ets ! vertically-summed heat content for ice & snow
	497	REAL(wp), INTENT( out), DIMENSION (jpij,jkmax) :: etilayer ! heat content for ice layers
[825]	498	!!
[1572]	499	INTEGER :: ji,jk ! loop indices
[825]	500	!!-----------------------------------------------------------------------
[2715]	501	eti(:,:) = 0._wp
	502	ets(:,:) = 0._wp
	503	!
[1572]	504	DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2]
[825]	505	DO ji = kideb, kiut
[4045]	506	etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
[1572]	507	eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk)
[825]	508	END DO
	509	END DO
[1572]	510	DO ji = kideb, kiut ! total q over all layers, snow [J.m-2]
[4045]	511	ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s )
[825]	512	END DO
[2715]	513	!
[4332]	514	WRITE(numout,*) ' lim_thd_glohec '
	515	WRITE(numout,) ' qt_i_in : ', eti(jiindex_1d,jl) r1_rdtice
	516	WRITE(numout,) ' qt_s_in : ', ets(jiindex_1d,jl) r1_rdtice
	517	WRITE(numout,) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) r1_rdtice
[1572]	518	!
[825]	519	END SUBROUTINE lim_thd_glohec
	520
	521
[1572]	522	SUBROUTINE lim_thd_con_dif( kideb, kiut, jl )
[825]	523	!!-----------------------------------------------------------------------
	524	!! * ROUTINE lim_thd_con_dif *
	525	!!
	526	!! ** Purpose : Test energy conservation after heat diffusion
	527	!!-------------------------------------------------------------------
[1572]	528	INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop
	529	INTEGER , INTENT(in ) :: jl ! category number
[825]	530
[1572]	531	INTEGER :: ji, jk ! loop indices
[4045]	532	INTEGER :: ii, ij
[1572]	533	INTEGER :: numce ! number of points for which conservation is violated
	534	REAL(wp) :: meance ! mean conservation error
	535	REAL(wp) :: max_cons_err, max_surf_err
[825]	536	!!---------------------------------------------------------------------
	537
[2715]	538	max_cons_err = 1.0_wp ! maximum tolerated conservation error
	539	max_surf_err = 0.001_wp ! maximum tolerated surface error
[921]	540
[825]	541	!--------------------------
	542	! Increment of energy
	543	!--------------------------
	544	! global
	545	DO ji = kideb, kiut
[1572]	546	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
[825]	547	END DO
	548	! layer by layer
[2528]	549	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
[825]	550
	551	!----------------------------------------
	552	! Atmospheric heat flux, ice heat budget
	553	!----------------------------------------
	554	DO ji = kideb, kiut
[4045]	555	ii = MOD( npb(ji) - 1 , jpi ) + 1
	556	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	557	fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji)
[4045]	558	sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl)
[825]	559	END DO
	560
	561	!--------------------
	562	! Conservation error
	563	!--------------------
	564	DO ji = kideb, kiut
[3625]	565	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	566	END DO
	567
[2528]	568	numce = 0
[2715]	569	meance = 0._wp
[825]	570	DO ji = kideb, kiut
[921]	571	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
	572	numce = numce + 1
	573	meance = meance + cons_error(ji,jl)
	574	ENDIF
[2528]	575	END DO
[2715]	576	IF( numce > 0 ) meance = meance / numce
[825]	577
	578	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
	579	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
[1572]	580	WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit
[825]	581	WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit
	582
	583	!-------------------------------------------------------
	584	! Surface error due to imbalance between Fatm and Fcsu
	585	!-------------------------------------------------------
[2528]	586	numce = 0
[2715]	587	meance = 0._wp
[825]	588
	589	DO ji = kideb, kiut
	590	surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
[2528]	591	IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN
[825]	592	numce = numce + 1
	593	meance = meance + surf_error(ji,jl)
	594	ENDIF
	595	ENDDO
[2715]	596	IF( numce > 0 ) meance = meance / numce
[825]	597
	598	WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
	599	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
	600	WRITE(numout,*) ' Mean surface error on big error points ', meance, numit
	601	WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
	602
[4332]	603	WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)
	604	WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl)
	605	WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d)
[825]	606
	607	!---------------------------------------
	608	! Write ice state in case of big errors
	609	!---------------------------------------
	610	DO ji = kideb, kiut
	611	IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
[921]	612	( cons_error(ji,jl) .GT. max_cons_err ) ) THEN
[4045]	613	ii = MOD( npb(ji) - 1, jpi ) + 1
	614	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	615	!
[921]	616	WRITE(numout,*) ' alerte 1 '
	617	WRITE(numout,*) ' Untolerated conservation / surface error after '
	618	WRITE(numout,*) ' heat diffusion in the ice '
	619	WRITE(numout,*) ' Category : ', jl
[4045]	620	WRITE(numout,*) ' ii , ij : ', ii, ij
	621	WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij)
[921]	622	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
	623	WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
[3625]	624	WRITE(numout,) ' dq_i : ', - dq_i(ji,jl) r1_rdtice
[921]	625	WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl)
	626	WRITE(numout,*)
	627	! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl)
	628	! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl)
	629	! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl)
	630	! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl)
[2528]	631	! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl)
[921]	632	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
	633	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
	634	WRITE(numout,*) ' t_su : ', t_su_b(ji)
	635	WRITE(numout,*) ' t_s : ', t_s_b(ji,1)
	636	WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i)
	637	WRITE(numout,*) ' t_bo : ', t_bo_b(ji)
	638	WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i)
	639	WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i)
	640	WRITE(numout,) ' tmelts : ', rtt - tmuts_i_b(ji,1:nlay_i)
	641	WRITE(numout,*)
	642	WRITE(numout,*) ' Fluxes '
	643	WRITE(numout,*) ' ~~~~~~ '
	644	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
	645	WRITE(numout,*) ' fc_su : ', fc_su (ji)
	646	WRITE(numout,) ' fstr_inice : ', qsr_ice_1d(ji)i0(ji)
	647	WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji)
	648	WRITE(numout,*) ' foc : ', fbif_1d(ji)
[4045]	649	WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl)
[921]	650	WRITE(numout,*) ' i0 : ', i0(ji)
	651	WRITE(numout,) ' qsr_ice : ', (1.0-i0(ji))qsr_ice_1d(ji)
	652	WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji)
	653	WRITE(numout,*) ' Conduction fluxes : '
	654	WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s)
	655	WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i)
	656	WRITE(numout,*)
	657	WRITE(numout,*) ' Layer by layer ... '
[3625]	658	WRITE(numout,) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) r1_rdtice
[2528]	659	WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0)
[921]	660	DO jk = 1, nlay_i
	661	WRITE(numout,*) ' layer : ', jk
[3625]	662	WRITE(numout,) ' dq_ice : ', dq_i_layer(ji,jk) r1_rdtice
[921]	663	WRITE(numout,*) ' radab : ', radab(ji,jk)
[2528]	664	WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1)
	665	WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk)
[921]	666	END DO
[825]	667
	668	ENDIF
[2715]	669	!
[825]	670	END DO
[1572]	671	!
[825]	672	END SUBROUTINE lim_thd_con_dif
	673
	674
[2528]	675	SUBROUTINE lim_thd_con_dh( kideb, kiut, jl )
[825]	676	!!-----------------------------------------------------------------------
	677	!! * ROUTINE lim_thd_con_dh *
	678	!!
	679	!! ** Purpose : Test energy conservation after enthalpy redistr.
	680	!!-----------------------------------------------------------------------
[2715]	681	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	682	INTEGER, INTENT(in) :: jl ! category number
	683	!
	684	INTEGER :: ji ! loop indices
[4045]	685	INTEGER :: ii, ij, numce ! local integers
[2715]	686	REAL(wp) :: meance, max_cons_err !local scalar
[825]	687	!!---------------------------------------------------------------------
	688
[2715]	689	max_cons_err = 1._wp
[921]	690
[825]	691	!--------------------------
	692	! Increment of energy
	693	!--------------------------
	694	DO ji = kideb, kiut
[2715]	695	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) + qt_s_fin(ji,jl) - qt_s_in(ji,jl) ! global
[825]	696	END DO
[2715]	697	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) ! layer by layer
[825]	698
	699	!----------------------------------------
	700	! Atmospheric heat flux, ice heat budget
	701	!----------------------------------------
	702	DO ji = kideb, kiut
[4045]	703	ii = MOD( npb(ji) - 1 , jpi ) + 1
	704	ij = ( npb(ji) - 1 ) / jpi + 1
[825]	705
[2715]	706	fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux
[4045]	707	sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl)
[3625]	708	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	709	END DO
	710
	711	!--------------------
	712	! Conservation error
	713	!--------------------
	714	DO ji = kideb, kiut
[3625]	715	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	716	END DO
	717
	718	numce = 0
[2715]	719	meance = 0._wp
[825]	720	DO ji = kideb, kiut
[2715]	721	IF( cons_error(ji,jl) .GT. max_cons_err ) THEN
[921]	722	numce = numce + 1
	723	meance = meance + cons_error(ji,jl)
	724	ENDIF
[825]	725	ENDDO
[2715]	726	IF(numce > 0 ) meance = meance / numce
[825]	727
	728	WRITE(numout,*) ' Error report - Category : ', jl
	729	WRITE(numout,*) ' ~~~~~~~~~~~~ '
	730	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
	731	WRITE(numout,*) ' After lim_thd_ent, category : ', jl
[2715]	732	WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit
[825]	733	WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit
	734
	735	!---------------------------------------
	736	! Write ice state in case of big errors
	737	!---------------------------------------
	738	DO ji = kideb, kiut
	739	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
[4045]	740	ii = MOD( npb(ji) - 1, jpi ) + 1
	741	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	742	!
[921]	743	WRITE(numout,*) ' alerte 1 - category : ', jl
	744	WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
[4045]	745	WRITE(numout,*) ' ii , ij : ', ii, ij
	746	WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij)
[921]	747	WRITE(numout,) ' '
	748	WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl)
[3625]	749	WRITE(numout,) ' dq_t : ', - dq_i(ji,jl) r1_rdtice
	750	WRITE(numout,) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) r1_rdtice
	751	WRITE(numout,) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) r1_rdtice
[921]	752	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
	753	WRITE(numout,) ' '
	754	WRITE(numout,*) ' Fluxes --- : '
	755	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
	756	WRITE(numout,*) ' foce : ', fbif_1d(ji)
	757	WRITE(numout,*) ' fres : ', ftotal_fin(ji)
[4045]	758	WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl)
[921]	759	WRITE(numout,) ' '
	760	WRITE(numout,*) ' Heat contents --- : '
[3625]	761	WRITE(numout,) ' qt_s_in : ', qt_s_in(ji,jl) r1_rdtice
	762	WRITE(numout,) ' qt_i_in : ', qt_i_in(ji,jl) r1_rdtice
	763	WRITE(numout,) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) r1_rdtice
	764	WRITE(numout,) ' qt_s_fin : ', qt_s_fin(ji,jl) r1_rdtice
	765	WRITE(numout,) ' qt_i_fin : ', qt_i_fin(ji,jl) r1_rdtice
	766	WRITE(numout,) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) r1_rdtice
[921]	767	WRITE(numout,) ' '
	768	WRITE(numout,*) ' Ice variables --- : '
	769	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
	770	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
	771	WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
	772	WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)
	773	WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
	774	WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
[825]	775	ENDIF
[1572]	776	!
[825]	777	END DO
[1572]	778	!
[825]	779	END SUBROUTINE lim_thd_con_dh
	780
[1572]	781
	782	SUBROUTINE lim_thd_enmelt( kideb, kiut )
[825]	783	!!-----------------------------------------------------------------------
	784	!! * ROUTINE lim_thd_enmelt *
	785	!!
	786	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3)
	787	!!
	788	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	789	!!-------------------------------------------------------------------
[1572]	790	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	791	!!
[2715]	792	INTEGER :: ji, jk ! dummy loop indices
	793	REAL(wp) :: ztmelts ! local scalar
[825]	794	!!-------------------------------------------------------------------
[1572]	795	!
	796	DO jk = 1, nlay_i ! Sea ice energy of melting
[825]	797	DO ji = kideb, kiut
[1572]	798	ztmelts = - tmut * s_i_b(ji,jk) + rtt
	799	q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &
[2715]	800	& + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) &
[1572]	801	& - rcp * ( ztmelts-rtt ) )
	802	END DO
	803	END DO
	804	DO jk = 1, nlay_s ! Snow energy of melting
[2715]	805	DO ji = kideb, kiut
[825]	806	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
[1572]	807	END DO
	808	END DO
	809	!
[825]	810	END SUBROUTINE lim_thd_enmelt
	811
	812
	813	SUBROUTINE lim_thd_init
	814	!!-----------------------------------------------------------------------
	815	!! * ROUTINE lim_thd_init *
	816	!!
	817	!! ** Purpose : Physical constants and parameters linked to the ice
[1572]	818	!! thermodynamics
[825]	819	!!
	820	!! ** Method : Read the namicethd namelist and check the ice-thermo
[1572]	821	!! parameter values called at the first timestep (nit000)
[825]	822	!!
	823	!! ** input : Namelist namicether
[2528]	824	!!-------------------------------------------------------------------
[1572]	825	NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
[4045]	826	& hicmin, hiclim, &
[1572]	827	& sbeta , parlat, hakspl, hibspl, exld, &
	828	& hakdif, hnzst , thth , parsub, alphs, betas, &
[825]	829	& kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
	830	!!-------------------------------------------------------------------
[2528]	831	!
[1572]	832	IF(lwp) THEN
	833	WRITE(numout,*)
	834	WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
	835	WRITE(numout,*) '~~~~~~~'
	836	ENDIF
[2528]	837	!
[1572]	838	REWIND( numnam_ice ) ! read Namelist numnam_ice
[825]	839	READ ( numnam_ice , namicethd )
[2528]	840	!
[1572]	841	IF(lwp) THEN ! control print
[825]	842	WRITE(numout,*)
[1572]	843	WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation '
	844	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
	845	WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit
	846	WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi
	847	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb
	848	WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb
	849	WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb
	850	WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin
	851	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
	852	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
	853	WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
	854	WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat
	855	WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl
	856	WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl
	857	WRITE(numout,*)' exponent for leads-closure rate exld = ', exld
	858	WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif
	859	WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth
	860	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
	861	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
	862	WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs
	863	WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas
	864	WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i
	865	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd
	866	WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd
	867	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi
[825]	868	ENDIF
[1572]	869	!
[825]	870	rcdsn = hakdif * rcdsn
	871	rcdic = hakdif * rcdic
[1572]	872	!
[825]	873	END SUBROUTINE lim_thd_init
	874
	875	#else
[1572]	876	!!----------------------------------------------------------------------
[2528]	877	!! Default option Dummy module NO LIM3 sea-ice model
[1572]	878	!!----------------------------------------------------------------------
[825]	879	#endif
	880
	881	!!======================================================================
	882	END MODULE limthd

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

Download in other formats: