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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago
remove remaining bugs in LIM3, so that it can run in both regional and global config
Property svn:keywords set to `Id`
File size: 46.0 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
[4205]	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)
[4161]	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
[4333]	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
[4161]	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	!!-------------------------------------------------------------------
[4161]	96	IF( nn_timing == 1 ) CALL timing_start('limthd')
[2715]	97
[3294]	98	CALL wrk_alloc( jpi, jpj, zqlbsbq )
[2715]	99
[4161]	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]
[4333]	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
[4333]	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]
[4333]	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
[4333]	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
[4333]	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)
[4333]	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
[4333]	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 * ( &
[4161]	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
[4161]	185	& + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step
[2387]	186	& - pfrld(ji,jj)*betas sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting
[2528]	187	!
[825]	188	! Positive heat budget is used for bottom ablation
	189	zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) )
	190	!= 1 if positive heat budget
	191	zpareff = 1.0 - zinda * zfntlat
[1572]	192	!= 0 if ice and positive heat budget and 1 if one of those two is false
[4333]	193	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) )
[2528]	194	!
[825]	195	! Heat budget of the lead, energy transferred from ice to ocean
	196	qldif (ji,jj) = zpareff * qldif(ji,jj)
	197	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
[2528]	198	!
[1572]	199	! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
[4333]	200	qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
[2528]	201	!
[1572]	202	! oceanic heat flux (limthd_dh)
[4333]	203	fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) )
[1571]	204	!
[825]	205	END DO
	206	END DO
	207
[921]	208	!------------------------------------------------------------------------------!
	209	! 3) Select icy points and fulfill arrays for the vectorial grid.
	210	!------------------------------------------------------------------------------!
[825]	211
	212	DO jl = 1, jpl !loop over ice categories
	213
[921]	214	IF( kt == nit000 .AND. lwp ) THEN
	215	WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl
	216	WRITE(numout,*) ' ~~~~~~~~'
	217	ENDIF
[825]	218
[4333]	219	zareamin = epsi10
[825]	220	nbpb = 0
	221	DO jj = 1, jpj
	222	DO ji = 1, jpi
	223	IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN
	224	nbpb = nbpb + 1
	225	npb(nbpb) = (jj - 1) * jpi + ji
	226	ENDIF
	227	END DO
	228	END DO
	229
[4333]	230	! debug point to follow
	231	jiindex_1d = 0
	232	IF( ln_nicep ) THEN
	233	DO ji = mi0(jiindx), mi1(jiindx)
	234	DO jj = mj0(jjindx), mj1(jjindx)
	235	jiindex_1d = (jj - 1) * jpi + ji
	236	END DO
	237	END DO
	238	ENDIF
	239
[921]	240	!------------------------------------------------------------------------------!
	241	! 4) Thermodynamic computation
	242	!------------------------------------------------------------------------------!
[825]	243
[2715]	244	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
[869]	245
[1572]	246	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
[825]	247
[921]	248	!-------------------------
	249	! 4.1 Move to 1D arrays
	250	!-------------------------
[825]	251
[1572]	252	CALL tab_2d_1d( nbpb, at_i_b (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
	253	CALL tab_2d_1d( nbpb, a_i_b (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	254	CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	255	CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	256
[1572]	257	CALL tab_2d_1d( nbpb, t_su_b (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	258	CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	259	DO jk = 1, nlay_s
[1572]	260	CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	261	CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	262	END DO
	263	DO jk = 1, nlay_i
[1572]	264	CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	265	CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	266	CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	267	END DO
	268
[1572]	269	CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) )
	270	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	271	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
	272	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
	273	CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	274	#if ! defined key_coupled
[3625]	275	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	276	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
[825]	277	#endif
[3625]	278	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	279	CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
	280	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
	281	CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb), fbif , jpi, jpj, npb(1:nbpb) )
	282	CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb), qldif , jpi, jpj, npb(1:nbpb) )
	283	CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice , jpi, jpj, npb(1:nbpb) )
	284	CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw , jpi, jpj, npb(1:nbpb) )
	285	CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb), dmgwi , jpi, jpj, npb(1:nbpb) )
	286	CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb), zqlbsbq , jpi, jpj, npb(1:nbpb) )
[825]	287
[3625]	288	CALL tab_2d_1d( nbpb, sfx_thd_1d (1:nbpb), sfx_thd , jpi, jpj, npb(1:nbpb) )
	289	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
	290	CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb), fhbri , jpi, jpj, npb(1:nbpb) )
	291	CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb), fstric , jpi, jpj, npb(1:nbpb) )
	292	CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) )
[825]	293
[4161]	294	CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif
	295	CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif
[921]	296	!--------------------------------
	297	! 4.3) Thermodynamic processes
	298	!--------------------------------
	299
[4333]	300	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting
	301	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]	302
[1572]	303	! !---------------------------------!
	304	CALL lim_thd_dif( 1, nbpb, jl ) ! Ice/Snow Temperature profile !
	305	! !---------------------------------!
[825]	306
[1572]	307	CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh
[825]	308
[4333]	309	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )
	310	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dif( 1 , nbpb , jl )
[825]	311
[1572]	312	! !---------------------------------!
	313	CALL lim_thd_dh( 1, nbpb, jl ) ! Ice/Snow thickness !
	314	! !---------------------------------!
[825]	315
[1572]	316	! !---------------------------------!
	317	CALL lim_thd_ent( 1, nbpb, jl ) ! Ice/Snow enthalpy remapping !
	318	! !---------------------------------!
[825]	319
[1572]	320	! !---------------------------------!
	321	CALL lim_thd_sal( 1, nbpb ) ! Ice salinity computation !
	322	! !---------------------------------!
[825]	323
[921]	324	! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting
[4333]	325	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )
	326	IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dh ( 1 , nbpb , jl )
[825]	327
[921]	328	!--------------------------------
	329	! 4.4) Move 1D to 2D vectors
	330	!--------------------------------
[825]	331
[3625]	332	CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb) , jpi, jpj )
	333	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_b (1:nbpb) , jpi, jpj )
	334	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_b (1:nbpb) , jpi, jpj )
	335	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_b (1:nbpb) , jpi, jpj )
	336	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_b (1:nbpb) , jpi, jpj )
	337	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_b (1:nbpb) , jpi, jpj )
[825]	338	DO jk = 1, nlay_s
[3625]	339	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b (1:nbpb,jk), jpi, jpj)
	340	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b (1:nbpb,jk), jpi, jpj)
[825]	341	END DO
	342	DO jk = 1, nlay_i
[3625]	343	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b (1:nbpb,jk), jpi, jpj)
	344	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b (1:nbpb,jk), jpi, jpj)
	345	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b (1:nbpb,jk), jpi, jpj)
[825]	346	END DO
[3625]	347	CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
	348	CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
	349	CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
	350	CALL tab_1d_2d( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj )
	351	CALL tab_1d_2d( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj )
	352	CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
	353	CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
	354	CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
	355	CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
	356	CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
	357	CALL tab_1d_2d( nbpb, sfx_thd , npb, sfx_thd_1d(1:nbpb) , jpi, jpj )
[2528]	358	!
[1572]	359	IF( num_sal == 2 ) THEN
[3625]	360	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
	361	CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb) , jpi, jpj )
[825]	362	ENDIF
[2528]	363	!
[3625]	364	!+++++ temporary stuff for a dummy version
[825]	365	CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )
	366	CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )
	367	CALL tab_1d_2d( nbpb, fsup2D , npb, fsup (1:nbpb) , jpi, jpj )
	368	CALL tab_1d_2d( nbpb, focea2D , npb, focea (1:nbpb) , jpi, jpj )
	369	CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )
	370	CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )
[888]	371	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
[825]	372	!+++++
[2528]	373	!
[1572]	374	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	375	ENDIF
	376	!
	377	END DO
[825]	378
[921]	379	!------------------------------------------------------------------------------!
	380	! 5) Global variables, diagnostics
	381	!------------------------------------------------------------------------------!
[825]	382
	383	!------------------------
	384	! 5.1) Ice heat content
	385	!------------------------
[3625]	386	! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules)
[2715]	387	zcoef = 1._wp / ( unit_fac * REAL( nlay_i ) )
[825]	388	DO jl = 1, jpl
[921]	389	DO jk = 1, nlay_i
[1572]	390	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef
	391	END DO
	392	END DO
[825]	393
	394	!------------------------
	395	! 5.2) Snow heat content
	396	!------------------------
[3625]	397	! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules)
[2715]	398	zcoef = 1._wp / ( unit_fac * REAL( nlay_s ) )
[825]	399	DO jl = 1, jpl
	400	DO jk = 1, nlay_s
[1572]	401	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef
	402	END DO
	403	END DO
[825]	404
	405	!----------------------------------
	406	! 5.3) Change thickness to volume
	407	!----------------------------------
	408	CALL lim_var_eqv2glo
	409
	410	!--------------------------------------------
	411	! 5.4) Diagnostic thermodynamic growth rates
	412	!--------------------------------------------
[4161]	413	!clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes
	414	!clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
[825]	415
[4333]	416	IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
[825]	417
[2528]	418	IF(ln_ctl) THEN ! Control print
[867]	419	CALL prt_ctl_info(' ')
	420	CALL prt_ctl_info(' - Cell values : ')
	421	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	422	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd : cell area :')
	423	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :')
	424	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :')
	425	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd : vt_s :')
	426	DO jl = 1, jpl
[867]	427	CALL prt_ctl_info(' ')
[863]	428	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	429	CALL prt_ctl_info(' ~~~~~~~~~~')
	430	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_thd : a_i : ')
	431	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_thd : ht_i : ')
	432	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_thd : ht_s : ')
	433	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_thd : v_i : ')
	434	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_thd : v_s : ')
	435	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_thd : e_s : ')
	436	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_thd : t_su : ')
	437	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_thd : t_snow : ')
	438	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_thd : sm_i : ')
	439	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_thd : smv_i : ')
	440	DO jk = 1, nlay_i
[867]	441	CALL prt_ctl_info(' ')
[863]	442	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	443	CALL prt_ctl_info(' ~~~~~~~')
	444	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd : t_i : ')
	445	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd : e_i : ')
	446	END DO
	447	END DO
	448	ENDIF
[2528]	449	!
[4161]	450	! -------------------------------
	451	!- check conservation (C Rousset)
	452	IF (ln_limdiahsb) THEN
	453	zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
	454	zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
	455
	456	zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
	457	zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
	458
	459	zchk_vmin = glob_min(v_i)
	460	zchk_amax = glob_max(SUM(a_i,dim=3))
	461	zchk_amin = glob_min(a_i)
	462
	463	IF(lwp) THEN
	464	IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,) 'violation volume [m3/day] (limthd) = ',(zchk_v_i rday)
	465	IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,) 'violation saline [psum3/day] (limthd) = ',(zchk_smv * rday)
	466	IF ( zchk_vmin < 0. ) WRITE(numout,) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin 1.e-3)
	467	IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax
	468	IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin
	469	ENDIF
	470	ENDIF
	471	!- check conservation (C Rousset)
	472	! -------------------------------
	473	!
[3294]	474	CALL wrk_dealloc( jpi, jpj, zqlbsbq )
[2715]	475	!
[4161]	476	IF( nn_timing == 1 ) CALL timing_stop('limthd')
[825]	477	END SUBROUTINE lim_thd
	478
	479
[1572]	480	SUBROUTINE lim_thd_glohec( eti, ets, etilayer, kideb, kiut, jl )
[825]	481	!!-----------------------------------------------------------------------
	482	!! * ROUTINE lim_thd_glohec *
	483	!!
	484	!! ** Purpose : Compute total heat content for each category
	485	!! Works with 1d vectors only
[1572]	486	!!-----------------------------------------------------------------------
	487	INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop
	488	INTEGER , INTENT(in ) :: jl ! category number
	489	REAL(wp), INTENT( out), DIMENSION (jpij,jpl ) :: eti, ets ! vertically-summed heat content for ice & snow
	490	REAL(wp), INTENT( out), DIMENSION (jpij,jkmax) :: etilayer ! heat content for ice layers
[825]	491	!!
[1572]	492	INTEGER :: ji,jk ! loop indices
[825]	493	!!-----------------------------------------------------------------------
[2715]	494	eti(:,:) = 0._wp
	495	ets(:,:) = 0._wp
	496	!
[1572]	497	DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2]
[825]	498	DO ji = kideb, kiut
[4161]	499	etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
[1572]	500	eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk)
[825]	501	END DO
	502	END DO
[1572]	503	DO ji = kideb, kiut ! total q over all layers, snow [J.m-2]
[4161]	504	ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s )
[825]	505	END DO
[2715]	506	!
[4333]	507	WRITE(numout,*) ' lim_thd_glohec '
	508	WRITE(numout,) ' qt_i_in : ', eti(jiindex_1d,jl) r1_rdtice
	509	WRITE(numout,) ' qt_s_in : ', ets(jiindex_1d,jl) r1_rdtice
	510	WRITE(numout,) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) r1_rdtice
[1572]	511	!
[825]	512	END SUBROUTINE lim_thd_glohec
	513
	514
[1572]	515	SUBROUTINE lim_thd_con_dif( kideb, kiut, jl )
[825]	516	!!-----------------------------------------------------------------------
	517	!! * ROUTINE lim_thd_con_dif *
	518	!!
	519	!! ** Purpose : Test energy conservation after heat diffusion
	520	!!-------------------------------------------------------------------
[1572]	521	INTEGER , INTENT(in ) :: kideb, kiut ! bounds for the spatial loop
	522	INTEGER , INTENT(in ) :: jl ! category number
[825]	523
[1572]	524	INTEGER :: ji, jk ! loop indices
[4161]	525	INTEGER :: ii, ij
[1572]	526	INTEGER :: numce ! number of points for which conservation is violated
	527	REAL(wp) :: meance ! mean conservation error
	528	REAL(wp) :: max_cons_err, max_surf_err
[825]	529	!!---------------------------------------------------------------------
	530
[2715]	531	max_cons_err = 1.0_wp ! maximum tolerated conservation error
	532	max_surf_err = 0.001_wp ! maximum tolerated surface error
[921]	533
[825]	534	!--------------------------
	535	! Increment of energy
	536	!--------------------------
	537	! global
	538	DO ji = kideb, kiut
[1572]	539	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]	540	END DO
	541	! layer by layer
[2528]	542	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
[825]	543
	544	!----------------------------------------
	545	! Atmospheric heat flux, ice heat budget
	546	!----------------------------------------
	547	DO ji = kideb, kiut
[4161]	548	ii = MOD( npb(ji) - 1 , jpi ) + 1
	549	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	550	fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji)
[4161]	551	sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl)
[825]	552	END DO
	553
	554	!--------------------
	555	! Conservation error
	556	!--------------------
	557	DO ji = kideb, kiut
[3625]	558	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	559	END DO
	560
[2528]	561	numce = 0
[2715]	562	meance = 0._wp
[825]	563	DO ji = kideb, kiut
[921]	564	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
	565	numce = numce + 1
	566	meance = meance + cons_error(ji,jl)
	567	ENDIF
[2528]	568	END DO
[2715]	569	IF( numce > 0 ) meance = meance / numce
[825]	570
	571	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
	572	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
[1572]	573	WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit
[825]	574	WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit
	575
	576	!-------------------------------------------------------
	577	! Surface error due to imbalance between Fatm and Fcsu
	578	!-------------------------------------------------------
[2528]	579	numce = 0
[2715]	580	meance = 0._wp
[825]	581
	582	DO ji = kideb, kiut
	583	surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
[2528]	584	IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN
[825]	585	numce = numce + 1
	586	meance = meance + surf_error(ji,jl)
	587	ENDIF
	588	ENDDO
[2715]	589	IF( numce > 0 ) meance = meance / numce
[825]	590
	591	WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
	592	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
	593	WRITE(numout,*) ' Mean surface error on big error points ', meance, numit
	594	WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
	595
[4333]	596	WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)
	597	WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl)
	598	WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d)
[825]	599
	600	!---------------------------------------
	601	! Write ice state in case of big errors
	602	!---------------------------------------
	603	DO ji = kideb, kiut
	604	IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
[921]	605	( cons_error(ji,jl) .GT. max_cons_err ) ) THEN
[4161]	606	ii = MOD( npb(ji) - 1, jpi ) + 1
	607	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	608	!
[921]	609	WRITE(numout,*) ' alerte 1 '
	610	WRITE(numout,*) ' Untolerated conservation / surface error after '
	611	WRITE(numout,*) ' heat diffusion in the ice '
	612	WRITE(numout,*) ' Category : ', jl
[4161]	613	WRITE(numout,*) ' ii , ij : ', ii, ij
	614	WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij)
[921]	615	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
	616	WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
[3625]	617	WRITE(numout,) ' dq_i : ', - dq_i(ji,jl) r1_rdtice
[921]	618	WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl)
	619	WRITE(numout,*)
	620	! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl)
	621	! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl)
	622	! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl)
	623	! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl)
[2528]	624	! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl)
[921]	625	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
	626	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
	627	WRITE(numout,*) ' t_su : ', t_su_b(ji)
	628	WRITE(numout,*) ' t_s : ', t_s_b(ji,1)
	629	WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i)
	630	WRITE(numout,*) ' t_bo : ', t_bo_b(ji)
	631	WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i)
	632	WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i)
	633	WRITE(numout,) ' tmelts : ', rtt - tmuts_i_b(ji,1:nlay_i)
	634	WRITE(numout,*)
	635	WRITE(numout,*) ' Fluxes '
	636	WRITE(numout,*) ' ~~~~~~ '
	637	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
	638	WRITE(numout,*) ' fc_su : ', fc_su (ji)
	639	WRITE(numout,) ' fstr_inice : ', qsr_ice_1d(ji)i0(ji)
	640	WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji)
	641	WRITE(numout,*) ' foc : ', fbif_1d(ji)
[4161]	642	WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl)
[921]	643	WRITE(numout,*) ' i0 : ', i0(ji)
	644	WRITE(numout,) ' qsr_ice : ', (1.0-i0(ji))qsr_ice_1d(ji)
	645	WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji)
	646	WRITE(numout,*) ' Conduction fluxes : '
	647	WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s)
	648	WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i)
	649	WRITE(numout,*)
	650	WRITE(numout,*) ' Layer by layer ... '
[3625]	651	WRITE(numout,) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) r1_rdtice
[2528]	652	WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0)
[921]	653	DO jk = 1, nlay_i
	654	WRITE(numout,*) ' layer : ', jk
[3625]	655	WRITE(numout,) ' dq_ice : ', dq_i_layer(ji,jk) r1_rdtice
[921]	656	WRITE(numout,*) ' radab : ', radab(ji,jk)
[2528]	657	WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1)
	658	WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk)
[921]	659	END DO
[825]	660
	661	ENDIF
[2715]	662	!
[825]	663	END DO
[1572]	664	!
[825]	665	END SUBROUTINE lim_thd_con_dif
	666
	667
[2528]	668	SUBROUTINE lim_thd_con_dh( kideb, kiut, jl )
[825]	669	!!-----------------------------------------------------------------------
	670	!! * ROUTINE lim_thd_con_dh *
	671	!!
	672	!! ** Purpose : Test energy conservation after enthalpy redistr.
	673	!!-----------------------------------------------------------------------
[2715]	674	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	675	INTEGER, INTENT(in) :: jl ! category number
	676	!
	677	INTEGER :: ji ! loop indices
[4161]	678	INTEGER :: ii, ij, numce ! local integers
[2715]	679	REAL(wp) :: meance, max_cons_err !local scalar
[825]	680	!!---------------------------------------------------------------------
	681
[2715]	682	max_cons_err = 1._wp
[921]	683
[825]	684	!--------------------------
	685	! Increment of energy
	686	!--------------------------
	687	DO ji = kideb, kiut
[2715]	688	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]	689	END DO
[2715]	690	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) ! layer by layer
[825]	691
	692	!----------------------------------------
	693	! Atmospheric heat flux, ice heat budget
	694	!----------------------------------------
	695	DO ji = kideb, kiut
[4161]	696	ii = MOD( npb(ji) - 1 , jpi ) + 1
	697	ij = ( npb(ji) - 1 ) / jpi + 1
[825]	698
[2715]	699	fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux
[4161]	700	sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl)
[3625]	701	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	702	END DO
	703
	704	!--------------------
	705	! Conservation error
	706	!--------------------
	707	DO ji = kideb, kiut
[3625]	708	cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
[825]	709	END DO
	710
	711	numce = 0
[2715]	712	meance = 0._wp
[825]	713	DO ji = kideb, kiut
[2715]	714	IF( cons_error(ji,jl) .GT. max_cons_err ) THEN
[921]	715	numce = numce + 1
	716	meance = meance + cons_error(ji,jl)
	717	ENDIF
[825]	718	ENDDO
[2715]	719	IF(numce > 0 ) meance = meance / numce
[825]	720
	721	WRITE(numout,*) ' Error report - Category : ', jl
	722	WRITE(numout,*) ' ~~~~~~~~~~~~ '
	723	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
	724	WRITE(numout,*) ' After lim_thd_ent, category : ', jl
[2715]	725	WRITE(numout,*) ' Mean conservation error on big error points ', meance, numit
[825]	726	WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit
	727
	728	!---------------------------------------
	729	! Write ice state in case of big errors
	730	!---------------------------------------
	731	DO ji = kideb, kiut
	732	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
[4161]	733	ii = MOD( npb(ji) - 1, jpi ) + 1
	734	ij = ( npb(ji) - 1 ) / jpi + 1
[2528]	735	!
[921]	736	WRITE(numout,*) ' alerte 1 - category : ', jl
	737	WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
[4161]	738	WRITE(numout,*) ' ii , ij : ', ii, ij
	739	WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij)
[921]	740	WRITE(numout,) ' '
	741	WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl)
[3625]	742	WRITE(numout,) ' dq_t : ', - dq_i(ji,jl) r1_rdtice
	743	WRITE(numout,) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) r1_rdtice
	744	WRITE(numout,) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) r1_rdtice
[921]	745	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
	746	WRITE(numout,) ' '
	747	WRITE(numout,*) ' Fluxes --- : '
	748	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
	749	WRITE(numout,*) ' foce : ', fbif_1d(ji)
	750	WRITE(numout,*) ' fres : ', ftotal_fin(ji)
[4161]	751	WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl)
[921]	752	WRITE(numout,) ' '
	753	WRITE(numout,*) ' Heat contents --- : '
[3625]	754	WRITE(numout,) ' qt_s_in : ', qt_s_in(ji,jl) r1_rdtice
	755	WRITE(numout,) ' qt_i_in : ', qt_i_in(ji,jl) r1_rdtice
	756	WRITE(numout,) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) r1_rdtice
	757	WRITE(numout,) ' qt_s_fin : ', qt_s_fin(ji,jl) r1_rdtice
	758	WRITE(numout,) ' qt_i_fin : ', qt_i_fin(ji,jl) r1_rdtice
	759	WRITE(numout,) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) r1_rdtice
[921]	760	WRITE(numout,) ' '
	761	WRITE(numout,*) ' Ice variables --- : '
	762	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
	763	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
	764	WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
	765	WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)
	766	WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
	767	WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
[825]	768	ENDIF
[1572]	769	!
[825]	770	END DO
[1572]	771	!
[825]	772	END SUBROUTINE lim_thd_con_dh
	773
[1572]	774
	775	SUBROUTINE lim_thd_enmelt( kideb, kiut )
[825]	776	!!-----------------------------------------------------------------------
	777	!! * ROUTINE lim_thd_enmelt *
	778	!!
	779	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3)
	780	!!
	781	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	782	!!-------------------------------------------------------------------
[1572]	783	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	784	!!
[2715]	785	INTEGER :: ji, jk ! dummy loop indices
	786	REAL(wp) :: ztmelts ! local scalar
[825]	787	!!-------------------------------------------------------------------
[1572]	788	!
	789	DO jk = 1, nlay_i ! Sea ice energy of melting
[825]	790	DO ji = kideb, kiut
[1572]	791	ztmelts = - tmut * s_i_b(ji,jk) + rtt
	792	q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &
[2715]	793	& + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) &
[1572]	794	& - rcp * ( ztmelts-rtt ) )
	795	END DO
	796	END DO
	797	DO jk = 1, nlay_s ! Snow energy of melting
[2715]	798	DO ji = kideb, kiut
[825]	799	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
[1572]	800	END DO
	801	END DO
	802	!
[825]	803	END SUBROUTINE lim_thd_enmelt
	804
	805
	806	SUBROUTINE lim_thd_init
	807	!!-----------------------------------------------------------------------
	808	!! * ROUTINE lim_thd_init *
	809	!!
	810	!! ** Purpose : Physical constants and parameters linked to the ice
[1572]	811	!! thermodynamics
[825]	812	!!
	813	!! ** Method : Read the namicethd namelist and check the ice-thermo
[1572]	814	!! parameter values called at the first timestep (nit000)
[825]	815	!!
	816	!! ** input : Namelist namicether
[2528]	817	!!-------------------------------------------------------------------
[4147]	818	INTEGER :: ios ! Local integer output status for namelist read
[1572]	819	NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
[4161]	820	& hicmin, hiclim, &
[1572]	821	& sbeta , parlat, hakspl, hibspl, exld, &
	822	& hakdif, hnzst , thth , parsub, alphs, betas, &
[825]	823	& kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
	824	!!-------------------------------------------------------------------
[2528]	825	!
[1572]	826	IF(lwp) THEN
	827	WRITE(numout,*)
	828	WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
	829	WRITE(numout,*) '~~~~~~~'
	830	ENDIF
[2528]	831	!
[4147]	832	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
	833	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
	834	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
	835
	836	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
	837	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
	838	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
	839	WRITE ( numoni, 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: