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

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 4166

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

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

Download in other formats: