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.

limistate.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6772

Last change on this file since 6772 was 6772, checked in by cbricaud, 8 years ago
clean in coarsening branch
Property svn:keywords set to `Id`
File size: 44.2 KB

Rev	Line
[825]	1	MODULE limistate
	2	!!======================================================================
	3	!! * MODULE limistate *
	4	!! Initialisation of diagnostics ice variables
	5	!!======================================================================
[2528]	6	!! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code
[4161]	7	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[4990]	8	!! - ! 2014 (C. Rousset) add N/S initializations
[2528]	9	!!----------------------------------------------------------------------
[825]	10	#if defined key_lim3
	11	!!----------------------------------------------------------------------
[834]	12	!! 'key_lim3' : LIM3 sea-ice model
[825]	13	!!----------------------------------------------------------------------
	14	!! lim_istate : Initialisation of diagnostics ice variables
	15	!! lim_istate_init : initialization of ice state and namelist read
	16	!!----------------------------------------------------------------------
[2528]	17	USE phycst ! physical constant
	18	USE oce ! dynamics and tracers variables
	19	USE dom_oce ! ocean domain
	20	USE sbc_oce ! Surface boundary condition: ocean fields
[4161]	21	USE sbc_ice ! Surface boundary condition: ice fields
[2528]	22	USE eosbn2 ! equation of state
	23	USE ice ! sea-ice variables
[4161]	24	USE par_oce ! ocean parameters
[2528]	25	USE dom_ice ! sea-ice domain
	26	USE in_out_manager ! I/O manager
[2715]	27	USE lib_mpp ! MPP library
[4161]	28	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3294]	29	USE wrk_nemo ! work arrays
[6772]	30	USE fldread ! read input fields
	31	USE iom
[825]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
[2528]	36	PUBLIC lim_istate ! routine called by lim_init.F90
[825]	37
[4166]	38	! !! init namelist (namiceini)
[5602]	39	REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice
	40	REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north
	41	REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south
	42	REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north
	43	REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south
	44	REAL(wp) :: rn_ati_ini_n ! initial leads area in the north
	45	REAL(wp) :: rn_ati_ini_s ! initial leads area in the south
	46	REAL(wp) :: rn_smi_ini_n ! initial salinity
	47	REAL(wp) :: rn_smi_ini_s ! initial salinity
	48	REAL(wp) :: rn_tmi_ini_n ! initial temperature
	49	REAL(wp) :: rn_tmi_ini_s ! initial temperature
[4166]	50
[6772]	51	INTEGER , PARAMETER :: jpfldi = 7 ! maximum number of files to read
	52	INTEGER , PARAMETER :: jp_hicif = 1 ! index of thick (m) at T-point
	53	INTEGER , PARAMETER :: jp_hsnif = 2 ! index of thick (m) at T-point
	54	INTEGER , PARAMETER :: jp_frld = 3 ! index of ice fraction (%) at T-point
	55	INTEGER , PARAMETER :: jp_sist = 4 ! index of ice surface temp (K) at T-point
	56	INTEGER , PARAMETER :: jp_tbif1 = 5 ! index of ice temp lev1 (K) at T-point
	57	INTEGER , PARAMETER :: jp_tbif2 = 6 ! index of ice temp lev2 (K) at T-point
	58	INTEGER , PARAMETER :: jp_tbif3 = 7 ! index of ice temp lev3 (K) at T-point
	59	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	60
	61	REAL(wp),DIMENSION(:,:) ,ALLOCATABLE :: hicif_ini,hsnif_ini,frld_ini,sist_ini, zswitch
	62	REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: tbif_ini
	63
[5602]	64	LOGICAL :: ln_iceini ! initialization or not
[6772]	65	LOGICAL :: ln_limini_file ! Ice initialization state from 2D netcdf file
[825]	66	!!----------------------------------------------------------------------
[4161]	67	!! LIM 3.0, UCL-LOCEAN-IPSL (2008)
[1156]	68	!! $Id$
[4161]	69	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]	70	!!----------------------------------------------------------------------
	71	CONTAINS
	72
	73	SUBROUTINE lim_istate
	74	!!-------------------------------------------------------------------
	75	!! * ROUTINE lim_istate *
	76	!!
	77	!! ** Purpose : defined the sea-ice initial state
	78	!!
[4161]	79	!! ** Method :
	80	!! This routine will put some ice where ocean
	81	!! is at the freezing point, then fill in ice
	82	!! state variables using prescribed initial
	83	!! values in the namelist
	84	!!
	85	!! ** Steps :
	86	!! 1) Read namelist
	87	!! 2) Basal temperature; ice and hemisphere masks
	88	!! 3) Fill in the ice thickness distribution using gaussian
	89	!! 4) Fill in space-dependent arrays for state variables
	90	!! 5) Diagnostic arrays
	91	!! 6) Lateral boundary conditions
	92	!!
[4333]	93	!! ** Notes : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
[4990]	94	!! where there is no ice (clem: I do not know why, is it mandatory?)
[4333]	95	!!
[4161]	96	!! History :
	97	!! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code
	98	!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
	99	!! 4.0 ! 09-11 (M. Vancoppenolle) Enhanced version for ice cats
	100	!!--------------------------------------------------------------------
	101
	102	!! * Local variables
	103	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[5602]	104	REAL(wp) :: ztmelts, zdh
[4161]	105	INTEGER :: i_hemis, i_fill, jl0
	106	REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
[4688]	107	REAL(wp), POINTER, DIMENSION(:) :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
	108	REAL(wp), POINTER, DIMENSION(:,:) :: zh_i_ini, za_i_ini, zv_i_ini
[4161]	109	INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index
[825]	110	!--------------------------------------------------------------------
[921]	111
[4161]	112	CALL wrk_alloc( jpi, jpj, zhemis )
[4688]	113	CALL wrk_alloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini )
	114	CALL wrk_alloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[2715]	115
[4161]	116	IF(lwp) WRITE(numout,*)
	117	IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
	118	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
	119
[825]	120	!--------------------------------------------------------------------
[4161]	121	! 1) Read namelist
[825]	122	!--------------------------------------------------------------------
	123
	124	CALL lim_istate_init ! reading the initials parameters of the ice
	125
[4688]	126	! surface temperature
	127	DO jl = 1, jpl ! loop over categories
[5602]	128	t_su (:,:,jl) = rt0 * tmask(:,:,1)
	129	tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
[4688]	130	END DO
	131
[4990]	132	! basal temperature (considered at freezing point)
[5602]	133	t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1)
[4990]	134
[5602]	135	IF( ln_iceini ) THEN
[4688]	136
[825]	137	!--------------------------------------------------------------------
[4161]	138	! 2) Basal temperature, ice mask and hemispheric index
[825]	139	!--------------------------------------------------------------------
[4990]	140
	141	DO jj = 1, jpj ! ice if sst <= t-freez + ttest
[1037]	142	DO ji = 1, jpi
[5602]	143	IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN
	144	zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice
[4765]	145	ELSE
[5602]	146	zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice
[1037]	147	ENDIF
	148	END DO
	149	END DO
	150
	151
[4161]	152	! Hemispheric index
[825]	153	DO jj = 1, jpj
	154	DO ji = 1, jpi
[2528]	155	IF( fcor(ji,jj) >= 0._wp ) THEN
[4161]	156	zhemis(ji,jj) = 1 ! Northern hemisphere
	157	ELSE
	158	zhemis(ji,jj) = 2 ! Southern hemisphere
	159	ENDIF
	160	END DO
	161	END DO
[825]	162
[4161]	163	!--------------------------------------------------------------------
	164	! 3) Initialization of sea ice state variables
	165	!--------------------------------------------------------------------
[6772]	166	IF( ln_limini_file )THEN
[825]	167
[6772]	168	CALL limini_file
	169
	170	ELSE
	171
[4161]	172	!-----------------------------
	173	! 3.1) Hemisphere-dependent arrays
	174	!-----------------------------
[4688]	175	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
[5602]	176	zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s ! ice thickness
	177	zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s ! snow depth
	178	zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s ! ice concentration
	179	zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s ! bulk ice salinity
	180	ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s ! temperature (ice and snow)
[825]	181
[4688]	182	zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:) ! ice volume
	183
[4161]	184	!---------------------------------------------------------------------
	185	! 3.2) Distribute ice concentration and thickness into the categories
	186	!---------------------------------------------------------------------
	187	! a gaussian distribution for ice concentration is used
	188	! then we check whether the distribution fullfills
	189	! volume and area conservation, positivity and ice categories bounds
	190	DO i_hemis = 1, 2
[825]	191
[4161]	192	ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
[825]	193
[4161]	194	! note for the great nemo engineers:
	195	! only very few of the WRITE statements are necessary for the reference version
	196	! they were one day useful, but now i personally doubt of their
	197	! potential for bringing anything useful
[825]	198
[4161]	199	DO i_fill = jpl, 1, -1
	200	IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
	201	!----------------------------
	202	! fill the i_fill categories
	203	!----------------------------
	204	! *** 1 category to fill
	205	IF ( i_fill .EQ. 1 ) THEN
[4688]	206	zh_i_ini(1,i_hemis) = zht_i_ini(i_hemis)
	207	za_i_ini(1,i_hemis) = zat_i_ini(i_hemis)
	208	zh_i_ini(2:jpl,i_hemis) = 0._wp
	209	za_i_ini(2:jpl,i_hemis) = 0._wp
[4161]	210	ELSE
[825]	211
[4688]	212	! *** >1 categores to fill
	213	!--- Ice thicknesses in the i_fill - 1 first categories
[4161]	214	DO jl = 1, i_fill - 1
[5602]	215	zh_i_ini(jl,i_hemis) = hi_mean(jl)
[4161]	216	END DO
[4688]	217
	218	!--- jl0: most likely index where cc will be maximum
[825]	219	DO jl = 1, jpl
[5602]	220	IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. &
	221	& ( zht_i_ini(i_hemis) <= hi_max(jl) ) ) THEN
[4161]	222	jl0 = jl
	223	ENDIF
	224	END DO
	225	jl0 = MIN(jl0, i_fill)
[4688]	226
	227	!--- Concentrations
[4161]	228	za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
	229	DO jl = 1, i_fill - 1
	230	IF ( jl .NE. jl0 ) THEN
[4688]	231	zsigma = 0.5 * zht_i_ini(i_hemis)
	232	zarg = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
[4161]	233	za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
	234	ENDIF
[4688]	235	END DO
	236
[4161]	237	zA = 0. ! sum of the areas in the jpl categories
	238	DO jl = 1, i_fill - 1
	239	zA = zA + za_i_ini(jl,i_hemis)
	240	END DO
	241	za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category
	242	IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
	243
[4688]	244	!--- Ice thickness in the last category
[4161]	245	zV = 0. ! sum of the volumes of the N-1 categories
	246	DO jl = 1, i_fill - 1
[4688]	247	zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
[4161]	248	END DO
[4688]	249	zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)
	250	IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[825]	251
[4688]	252	!--- volumes
	253	zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
[4161]	254	IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[921]	255
[4161]	256	ENDIF ! i_fill
[825]	257
[4161]	258	!---------------------
	259	! Compatibility tests
	260	!---------------------
	261	! Test 1: area conservation
	262	zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
	263	IF ( zconv .LT. 1.0e-6 ) THEN
	264	ztest_1 = 1
	265	ELSE
	266	! this write is useful
	267	IF(lwp) WRITE(numout,) ' TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis)
	268	ztest_1 = 0
	269	ENDIF
[825]	270
[4161]	271	! Test 2: volume conservation
	272	zV_cons = SUM(zv_i_ini(:,i_hemis))
	273	zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
[825]	274
[4161]	275	IF ( zconv .LT. 1.0e-6 ) THEN
	276	ztest_2 = 1
	277	ELSE
	278	! this write is useful
	279	IF(lwp) WRITE(numout,) ' TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
	280	' zvt_i_ini = ', zvt_i_ini(i_hemis)
	281	ztest_2 = 0
	282	ENDIF
[825]	283
[4161]	284	! Test 3: thickness of the last category is in-bounds ?
[5602]	285	IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
[4161]	286	ztest_3 = 1
	287	ELSE
	288	! this write is useful
[4688]	289	IF(lwp) WRITE(numout,) ' TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
	290	zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
[4161]	291	ztest_3 = 0
	292	ENDIF
[825]	293
[4161]	294	! Test 4: positivity of ice concentrations
	295	ztest_4 = 1
	296	DO jl = 1, jpl
	297	IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN
	298	! this write is useful
	299	IF(lwp) WRITE(numout,) ' TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
	300	ztest_4 = 0
[825]	301	ENDIF
[4161]	302	END DO
[825]	303
[4161]	304	ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
	305
	306	ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
[825]	307
[4161]	308	END DO ! i_fill
[825]	309
[4161]	310	IF(lwp) THEN
[4908]	311	WRITE(numout,*) ' ztests : ', ztests
[4161]	312	IF ( ztests .NE. 4 ) THEN
	313	WRITE(numout,*)
[4908]	314	WRITE(numout,*) ' !!!! ALERT !!! '
	315	WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
[4161]	316	WRITE(numout,*)
[4908]	317	WRITE(numout,) ' ** ztests is not equal to 4 '
	318	WRITE(numout,) ' ** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
	319	WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
	320	WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
[4161]	321	ENDIF ! ztests .NE. 4
	322	ENDIF
	323
	324	END DO ! i_hemis
[825]	325
[4161]	326	!---------------------------------------------------------------------
	327	! 3.3) Space-dependent arrays for ice state variables
	328	!---------------------------------------------------------------------
[825]	329
[4333]	330	! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
[4161]	331	DO jl = 1, jpl ! loop over categories
	332	DO jj = 1, jpj
	333	DO ji = 1, jpi
[4688]	334	a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration
[5602]	335	ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness
[4688]	336	ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) ) ! snow depth
[5602]	337	sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity
	338	o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day)
	339	t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
[4333]	340
	341	! This case below should not be used if (ht_s/ht_i) is ok in namelist
	342	! In case snow load is in excess that would lead to transformation from snow to ice
	343	! Then, transfer the snow excess into the ice (different from limthd_dh)
	344	zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
	345	! recompute ht_i, ht_s avoiding out of bounds values
	346	ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
[5602]	347	ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
[4333]	348
	349	! ice volume, salt content, age content
[4161]	350	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume
	351	v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume
	352	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
	353	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content
[5602]	354	END DO
	355	END DO
	356	END DO
[825]	357
[4161]	358	! Snow temperature and heat content
	359	DO jk = 1, nlay_s
	360	DO jl = 1, jpl ! loop over categories
	361	DO jj = 1, jpj
	362	DO ji = 1, jpi
[5602]	363	t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
[4161]	364	! Snow energy of melting
[5602]	365	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
[825]	366
[5602]	367	! Mutliply by volume, and divide by number of layers to get heat content in J/m2
	368	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
	369	END DO
	370	END DO
	371	END DO
	372	END DO
	373
[4161]	374	! Ice salinity, temperature and heat content
	375	DO jk = 1, nlay_i
	376	DO jl = 1, jpl ! loop over categories
	377	DO jj = 1, jpj
	378	DO ji = 1, jpi
[5602]	379	t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
	380	s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
	381	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
[825]	382
[4161]	383	! heat content per unit volume
[4688]	384	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
[5602]	385	+ lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
	386	- rcp * ( ztmelts - rt0 ) )
[825]	387
[5602]	388	! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
	389	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
	390	END DO
	391	END DO
	392	END DO
	393	END DO
[825]	394
[4688]	395	tn_ice (:,:,:) = t_su (:,:,:)
	396
[6772]	397	ENDIF !ln_limini_file
	398
[4688]	399	ELSE
[5602]	400	! if ln_iceini=false
[4688]	401	a_i (:,:,:) = 0._wp
	402	v_i (:,:,:) = 0._wp
	403	v_s (:,:,:) = 0._wp
	404	smv_i(:,:,:) = 0._wp
	405	oa_i (:,:,:) = 0._wp
	406	ht_i (:,:,:) = 0._wp
	407	ht_s (:,:,:) = 0._wp
	408	sm_i (:,:,:) = 0._wp
	409	o_i (:,:,:) = 0._wp
	410
	411	e_i(:,:,:,:) = 0._wp
	412	e_s(:,:,:,:) = 0._wp
	413
	414	DO jl = 1, jpl
	415	DO jk = 1, nlay_i
[5602]	416	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]	417	END DO
	418	DO jk = 1, nlay_s
[5602]	419	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]	420	END DO
	421	END DO
[6772]	422
[5602]	423	ENDIF ! ln_iceini
[4688]	424
	425	at_i (:,:) = 0.0_wp
	426	DO jl = 1, jpl
	427	at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
	428	END DO
	429	!
[825]	430	!--------------------------------------------------------------------
[4161]	431	! 4) Global ice variables for output diagnostics \|
[825]	432	!--------------------------------------------------------------------
[4161]	433	u_ice (:,:) = 0._wp
	434	v_ice (:,:) = 0._wp
	435	stress1_i(:,:) = 0._wp
	436	stress2_i(:,:) = 0._wp
	437	stress12_i(:,:) = 0._wp
[825]	438
	439	!--------------------------------------------------------------------
[4161]	440	! 5) Moments for advection
[825]	441	!--------------------------------------------------------------------
	442
[4161]	443	sxopw (:,:) = 0._wp
	444	syopw (:,:) = 0._wp
	445	sxxopw(:,:) = 0._wp
	446	syyopw(:,:) = 0._wp
	447	sxyopw(:,:) = 0._wp
[3610]	448
[4161]	449	sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp
	450	syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp
	451	sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp
	452	syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp
	453	sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp
[825]	454
[4161]	455	sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp
	456	syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp
	457	sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp
	458	syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp
	459	sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp
[825]	460
[4161]	461	sxsal (:,:,:) = 0._wp
	462	sysal (:,:,:) = 0._wp
	463	sxxsal (:,:,:) = 0._wp
	464	syysal (:,:,:) = 0._wp
	465	sxysal (:,:,:) = 0._wp
[825]	466
[4335]	467	sxage (:,:,:) = 0._wp
	468	syage (:,:,:) = 0._wp
	469	sxxage (:,:,:) = 0._wp
	470	syyage (:,:,:) = 0._wp
	471	sxyage (:,:,:) = 0._wp
	472
[825]	473
[4161]	474	CALL wrk_dealloc( jpi, jpj, zhemis )
[4688]	475	CALL wrk_dealloc( jpl, 2, zh_i_ini, za_i_ini, zv_i_ini )
	476	CALL wrk_dealloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[4161]	477
[825]	478	END SUBROUTINE lim_istate
	479
	480	SUBROUTINE lim_istate_init
	481	!!-------------------------------------------------------------------
	482	!! * ROUTINE lim_istate_init *
	483	!!
	484	!! ** Purpose : Definition of initial state of the ice
	485	!!
[4161]	486	!! ** Method : Read the namiceini namelist and check the parameter
	487	!! values called at the first timestep (nit000)
[825]	488	!!
[4161]	489	!! ** input :
	490	!! Namelist namiceini
	491	!!
	492	!! history :
	493	!! 8.5 ! 03-08 (C. Ethe) original code
	494	!! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization
[825]	495	!!-----------------------------------------------------------------------------
[6772]	496	!
	497	INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read
	498	INTEGER :: ji,jj
	499	INTEGER :: ifpr, ierror
	500	!
	501	CHARACTER(len=100) :: cn_dir ! Root directory for location of ice files
	502	TYPE(FLD_N) :: sn_hicif, sn_hsnif, sn_frld, sn_sist
	503	TYPE(FLD_N) :: sn_tbif1, sn_tbif2, sn_tbif3
	504	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	505	!
	506	NAMELIST/namiceini/ ln_iceini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, &
	507	& rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
	508	& rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, &
	509	& sn_hicif, sn_hsnif, sn_frld, sn_sist, &
	510	& sn_tbif1, sn_tbif2, sn_tbif3, cn_dir
[825]	511	!!-----------------------------------------------------------------------------
[2528]	512	!
[4147]	513	REWIND( numnam_ice_ref ) ! Namelist namiceini in reference namelist : Ice initial state
	514	READ ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
	515	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
	516
	517	REWIND( numnam_ice_cfg ) ! Namelist namiceini in configuration namelist : Ice initial state
	518	READ ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
	519	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
[4624]	520	IF(lwm) WRITE ( numoni, namiceini )
[4161]	521
[6772]	522	slf_i(jp_hicif) = sn_hicif ; slf_i(jp_hsnif) = sn_hsnif
	523	slf_i(jp_frld) = sn_frld ; slf_i(jp_sist) = sn_sist
	524	slf_i(jp_tbif1) = sn_tbif1 ; slf_i(jp_tbif2) = sn_tbif2 ; slf_i(jp_tbif3) = sn_tbif3
	525
[4161]	526	! Define the initial parameters
	527	! -------------------------
	528
	529	IF(lwp) THEN
[825]	530	WRITE(numout,*)
	531	WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
	532	WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5602]	533	WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini
[6772]	534	WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini_file = ', ln_limini_file
[5602]	535	WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst
	536	WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n
	537	WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s
	538	WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n
	539	WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s
	540	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n
	541	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s
	542	WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n
	543	WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s
	544	WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n
	545	WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s
[825]	546	ENDIF
[4161]	547
[6772]	548	IF( ln_limini_file ) THEN ! Ice initialization using input file
	549	!
	550	ierr = alloc_lim_istate_init()
	551	!
	552	! CALL iom_open( 'Ice_initialization.nc', inum_ice )
	553	! !
	554	! IF( inum_ice > 0 ) THEN
	555	! IF(lwp) WRITE(numout,*)
	556	! IF(lwp) WRITE(numout,*) ' ice state initialization with : Ice_initialization.nc'
	557	!
	558	! CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif_ini )
	559	! CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif_ini )
	560	! CALL iom_get( inum_ice, jpdom_data, 'frld' , frld_ini )
	561	! CALL iom_get( inum_ice, jpdom_data, 'ts' , sist_ini )
	562	! CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif_ini(1:nlci,1:nlcj,:), &
	563	! & kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,3 /) )
	564	! ! put some values in the extra-halo...
	565
	566	! set si structure
	567	ALLOCATE( si(jpfldi), STAT=ierror )
	568	IF( ierror > 0 ) THEN
	569	CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' ) ; RETURN
	570	ENDIF
	571
	572	DO ifpr= 1, jpfldi
	573	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
	574	ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
	575	END DO
	576
	577	! fill si with slf_i and control print
	578	CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
	579
	580	CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step
	581
	582	hicif_ini(:,:) = si(jp_hicif)%fnow(:,:,1)
	583	hsnif_ini(:,:) = si(jp_hsnif)%fnow(:,:,1)
	584	frld_ini(:,:) = si(jp_frld)%fnow(:,:,1)
	585	sist_ini(:,:) = si(jp_sist)%fnow(:,:,1)
	586	tbif_ini(:,:,1) = si(jp_tbif1)%fnow(:,:,1)
	587	tbif_ini(:,:,2) = si(jp_tbif2)%fnow(:,:,1)
	588	tbif_ini(:,:,3) = si(jp_tbif3)%fnow(:,:,1)
	589
	590	DO jj = nlcj+1, jpj ; tbif_ini(1:nlci,jj,:) = tbif_ini(1:nlci,nlej,:) ; END DO
	591	DO ji = nlci+1, jpi ; tbif_ini(ji ,: ,:) = tbif_ini(nlei ,: ,:) ; END DO
	592
	593	! CALL iom_close( inum_ice)
	594	! !
	595	! ENDIF
	596	ENDIF
	597
[825]	598	END SUBROUTINE lim_istate_init
	599
[6772]	600	SUBROUTINE limini_file
	601	!!-----------------------------------------------------------------------------
	602	!!
	603	!!
	604	!!
	605	!!
	606	!!-----------------------------------------------------------------------------
	607	INTEGER :: jl,ji,jj,jk
	608	INTEGER :: jl0
	609	INTEGER :: i_fill,jit,jjt
	610	REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv,zH
	611	REAL(wp) :: eps=1.e-6
	612	REAL(wp) :: zmin,zmax
	613	!rbb REAL(wp) :: epsi20,ztmelts,zdh
	614	REAL(wp) ::ztmelts,zdh
	615
	616	REAL(wp), POINTER, DIMENSION(:,:) :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini
	617	REAL(wp), POINTER, DIMENSION(:,:,:) :: zv_i_ini
	618	REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ini,za_i_ini
	619	REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator
	620	!-----------------------------------------------------------------------------
	621	IF(lwp)WRITE(numout,*)"limistate: read file : "
	622
	623	CALL wrk_alloc(jpl,jpi,jpj, zv_i_ini)
	624	CALL wrk_alloc( jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
	625	CALL wrk_alloc( jpl,jpi,jpj,zht_i_ini,za_i_ini)
	626	CALL wrk_alloc( jpi,jpj,zidto )
	627
	628	zhm_i_ini(:,:) = hicif_ini(:,:) ! ice thickness
	629	zat_i_ini(:,:) = 1._wp - frld_ini(:,:) ! ice concentration
	630	zvt_i_ini(:,:) = zhm_i_ini(:,:) * zat_i_ini(:,:) ! ice volume
	631	zhm_s_ini(:,:) = hsnif_ini(:,:) ! snow depth
	632
	633	zht_i_ini(:,:,:) = 0._wp
	634	za_i_ini(:,:,:) = 0._wp
	635	zv_i_ini(:,:,:) = 0._wp
	636
	637	zat_i_ini(:,:) = MIN( zat_i_ini(:,:) , 1.0_wp )
	638
	639
	640	DO ji = 1, jpi
	641	DO jj = 1, jpj
	642
	643	IF( zat_i_ini(ji,jj) .GT. 0._wp .AND. zhm_i_ini(ji,jj) .GT. 0._wp )THEN
	644
	645
	646	IF( gphit(ji,jj) .GE. 0._wp )THEN ; zsm_i_ini(ji,jj) = rn_smi_ini_n
	647	ELSE ; zsm_i_ini(ji,jj) = rn_smi_ini_s
	648	ENDIF
	649
	650	jl0 = 1
	651	DO jl = 2, jpl
	652	IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. &
	653	( zhm_i_ini(ji,jj) .LE. hi_max(jl) ) ) THEN
	654	jl0 = jl
	655	ENDIF
	656	END DO
	657
	658	IF( jl0==1 )THEN
	659
	660	zht_i_ini(1,ji,jj) = zhm_i_ini(ji,jj)
	661	za_i_ini(1,ji,jj) = zat_i_ini(ji,jj)
	662	zht_i_ini(2:jpl,ji,jj) = 0._wp
	663	za_i_ini(2:jpl,ji,jj) = 0._wp
	664
	665	ELSE ! jl0 ne 1
	666	ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
	667
	668	DO i_fill = jpl, 1, -1
	669	IF( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
	670
	671	!----------------------------
	672	! fill the i_fill categories
	673	!----------------------------
	674	! *** 1 category to fill
	675	IF( i_fill .EQ. 1 ) THEN
	676	zht_i_ini(1,ji,jj) = zhm_i_ini(ji,jj)
	677	za_i_ini(1,ji,jj) = zat_i_ini(ji,jj)
	678	zht_i_ini(2:jpl,ji,jj) = 0._wp
	679	za_i_ini(2:jpl,ji,jj) = 0._wp
	680	ELSE
	681
	682	! *** >1 categores to fill
	683	!--- Ice thicknesses in the i_fill - 1 first categories
	684	DO jl = 1, i_fill - 1
	685	zht_i_ini(jl,ji,jj) = 0.5 * ( hi_max(jl) + hi_max(jl-1) )
	686	END DO
	687
	688	!--- jl0: most likely index where cc will be maximum
	689	DO jl = 1, jpl
	690	IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. &
	691	( zhm_i_ini(ji,jj) .LE. hi_max(jl) ) ) THEN
	692	jl0 = jl
	693	ENDIF
	694	END DO
	695	jl0 = MIN(jl0, i_fill)
	696
	697	!--- Concentrations
	698	za_i_ini(jl0,ji,jj) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
	699	DO jl = 1, i_fill - 1
	700	IF ( jl .NE. jl0 ) THEN
	701	zsigma = 0.5 * zhm_i_ini(ji,jj)
	702	zarg = ( zht_i_ini(jl,ji,jj) - zhm_i_ini(ji,jj) ) / zsigma
	703	za_i_ini(jl,ji,jj) = za_i_ini(jl0,ji,jj) * EXP(-zarg**2)
	704	ENDIF
	705	END DO
	706
	707	zA = 0. ! sum of the areas in the jpl categories
	708	DO jl = 1, i_fill - 1
	709	zA = zA + za_i_ini(jl,ji,jj)
	710	END DO
	711	za_i_ini(i_fill,ji,jj) = zat_i_ini(ji,jj) - zA ! ice conc in the last category
	712	IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
	713
	714	!--- Ice thickness in the last category
	715	zV = 0. ! sum of the volumes of the N-1 categories
	716	DO jl = 1, i_fill - 1
	717	zV = zV + za_i_ini(jl,ji,jj)*zht_i_ini(jl,ji,jj)
	718	END DO
	719	zht_i_ini(i_fill,ji,jj) = ( zvt_i_ini(ji,jj) - zV ) /za_i_ini(i_fill,ji,jj)
	720	IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
	721
	722	!--- volumes
	723	zv_i_ini(:,ji,jj) = za_i_ini(:,ji,jj) * zht_i_ini(:,ji,jj)
	724	IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
	725
	726	ENDIF ! i_fill
	727
	728	!---------------------
	729	! Compatibility tests
	730	!---------------------
	731	! Test 1: area conservation
	732	zA_cons = SUM(za_i_ini(:,ji,jj)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons )
	733	IF ( zconv .LT. 1.0e-6 ) THEN
	734	ztest_1 = 1
	735	ELSE
	736	! this write is useful
	737	!WRITE(numout,) ' TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj)
	738	!WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
	739	!WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	740	!WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
	741	!WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
	742	!WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
	743	!WRITE(numout,*) ' hi_max ',hi_max
	744	!WRITE(numout,*) ' jl0 = ',jl0
	745	!WRITE(numout,*) ' vol = ',zvt_i_ini(ji,jj),SUM(zv_i_ini(:,ji,jj))
	746	ztest_1 = 0
	747	ENDIF
	748
	749	! Test 2: volume conservation
	750	zV_cons = SUM(zv_i_ini(:,ji,jj))
	751	zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
	752
	753	IF ( zconv .LT. 1.0e-6 ) THEN
	754	ztest_2 = 1
	755	ELSE
	756	! this write is useful
	757	!WRITE(numout,) ' TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
	758	! ' zvt_i_ini = ', zvt_i_ini(ji,jj)
	759	!WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
	760	!WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	761	!WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
	762	!WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
	763	!WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
	764	!WRITE(numout,*) ' hi_max ',hi_max
	765	!WRITE(numout,*) ' jl0 = ',jl0
	766	ztest_2 = 0
	767	ENDIF
	768
	769	! Test 3: thickness of the last category is in-bounds ?
	770	IF ( zht_i_ini(i_fill, ji,jj) .GT. hi_max(i_fill-1) ) THEN
	771	ztest_3 = 1
	772	ELSE
	773	! this write is useful
	774	!WRITE(numout,) ' TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,ji,jj) = ', &
	775	!zht_i_ini(i_fill,ji,jj), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
	776	!WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
	777	!WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	778	!WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
	779	!WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
	780	!WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
	781	!WRITE(numout,*) ' hi_max ',hi_max
	782	!WRITE(numout,*) ' jl0 = ',jl0
	783	ztest_3 = 0
	784	ENDIF
	785
	786	! Test 4: positivity of ice concentrations
	787	ztest_4 = 1
	788	DO jl = 1, jpl
	789	IF ( za_i_ini(jl,ji,jj) .LT. 0._wp ) THEN
	790	! this write is useful
	791	!WRITE(numout,) ' TEST 4 POSITIVITY NOT OK FOR CAT ', jl, 'WITH A = ', za_i_ini(jl,ji,jj)
	792	!WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
	793	!WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	794	!WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
	795	!WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
	796	!WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
	797	!WRITE(numout,*) ' hi_max ',hi_max
	798	!WRITE(numout,*) ' jl0 = ',jl0
	799	!WRITE(numout,*)
	800	ztest_4 = 0
	801	ENDIF
	802	END DO
	803
	804	ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
	805
	806	ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
	807
	808	END DO ! i_fill
	809
	810	!WRITE(numout,*) ' ztests : ', ztests
	811	!IF ( ztests .NE. 4 ) THEN
	812	!WRITE(numout,*)
	813	!WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
	814	!WRITE(numout,*) ' !!!! RED ALERT !!! '
	815	!WRITE(numout,*) ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!'
	816	!WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
	817	!WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
	818	!WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
	819	!WRITE(numout,) ' ** ztests is not equal to 4 '
	820	!WRITE(numout,) ' ** ztest_i (i=1,4) = ', ztest_1, ztest_2,ztest_3,ztest_4
	821	!WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	822	!WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
	823	!WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
	824	!WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
	825	!WRITE(numout,*) ' hi_max ',hi_max
	826	!ENDIF ! ztests .NE. 4
	827
	828	ENDIF ! jl0 ne 1
	829
	830	ENDIF ! zat_i_ini ne 0
	831	END DO ! jj
	832	END DO ! ji
	833
	834
	835	!---------------------------------------------------------------------
	836	! 3.3) Space-dependent arrays for ice state variables
	837	!---------------------------------------------------------------------
	838
	839	! Ice concentration, thickness and volume, ice salinity, ice age, surface
	840	! temperature
	841	DO jl = 1, jpl ! loop over categories
	842	DO jj = 1, jpj
	843	DO ji = 1, jpi
	844	a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,ji,jj) ! concentration
	845	ht_i(ji,jj,jl) = zswitch(ji,jj) * zht_i_ini(jl,ji,jj) !ice thickness
	846
	847	IF( zhm_i_ini( ji,jj ) .GT. 0_wp )THEN ; ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zhm_s_ini( ji,jj ) / zhm_i_ini( ji,jj ) )
	848	ELSE ; ht_s(ji,jj,jl) = 0._wp
	849	ENDIF
	850	sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ (1._wp - zswitch(ji,jj) ) * rn_simin ! salinity
	851	o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp + ( 1._wp -zswitch(ji,jj) ) ! age
	852	t_su(ji,jj,jl) = sist_ini(ji,jj)
	853
	854	! This case below should not be used if (ht_s/ht_i) is ok in
	855	! namelist
	856	! In case snow load is in excess that would lead to
	857	! transformation from snow to ice
	858	! Then, transfer the snow excess into the ice (different from
	859	! limthd_dh)
	860	zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) ht_i(ji,jj,jl) ) r1_rau0 )
	861	! recompute ht_i, ht_s avoiding out of bounds values
	862	ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
	863	ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic *r1_rhosn )
	864
	865	! ice volume, salt content, age content
	866	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) !ice volume
	867	v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) !snow volume
	868	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) *v_i(ji,jj,jl) ! salt content
	869	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) !age content
	870	END DO ! ji
	871	END DO ! jj
	872	END DO ! jl
	873
	874	!cbr
	875	DO jk = 1, nlay_s
	876	DO jl = 1, jpl ! loop over categories
	877	!rbb t_s(:,:,1,jl) = tbif_ini(:,:,1)
	878	t_s(:,:,1,jl) = tbif_ini(:,:,1)zswitch(:,:)+ ( 1._wp - zswitch(:,:) ) rt0
	879	END DO ! jl
	880	END DO ! jk
	881
	882	! Snow temperature and heat content
	883	DO jk = 1, nlay_s
	884	DO jl = 1, jpl ! loop over categories
	885	DO jj = 1, jpj
	886	DO ji = 1, jpi
	887	!cbr??? t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	888	! Snow energy of melting
	889	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
	890
	891	! Mutliply by volume, and divide by number of layers to get
	892	! heat content in J/m2
	893	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) *r1_nlay_s
	894	END DO ! ji
	895	END DO ! jj
	896	END DO ! jl
	897	END DO ! jk
	898
	899	! Ice salinity, temperature and heat content
	900	DO jk = 1, nlay_i
	901	DO jl = 1, jpl ! loop over categories
	902	DO jj = 1, jpj
	903	DO ji = 1, jpi
	904	!cbr??? t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	905	t_i(ji,jj,jk,jl) = tbif_ini(ji,jj,2)zswitch(ji,jj)+ ( 1._wp - zswitch(ji,jj) ) rt0
	906	s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
	907	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
	908
	909	! heat content per unit volume
	910	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
	911	+ lfus * ( 1._wp - (ztmelts-rt0) /MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
	912	- rcp * ( ztmelts - rt0 ) )
	913
	914	! Mutliply by ice volume, and divide by number of layers to
	915	! get heat content in J/m2
	916	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
	917	END DO ! ji
	918	END DO ! jj
	919	END DO ! jl
	920	END DO ! jk
	921
	922	!cbr tmp CALL wrk_dealloc(jpl,jpi,jpj, zht_i_ini, za_i_ini, zv_i_ini)
	923	CALL wrk_dealloc(jpl,jpi,jpj, zv_i_ini)
	924	CALL wrk_dealloc( jpl,jpi,jpj,zht_i_ini,za_i_ini)
	925	CALL wrk_dealloc( jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini,zsm_i_ini )
	926	CALL wrk_dealloc( jpi,jpj,zidto )
	927
	928	END SUBROUTINE limini_file
	929
	930
	931	INTEGER FUNCTION alloc_lim_istate_init()
	932	!!-----------------------------------------------------------------------------
	933	!!
	934	!!
	935	!!
	936	!!
	937	!!-----------------------------------------------------------------------------
	938	INTEGER :: ierr(1)
	939	!!-----------------------------------------------------------------------------
	940	ALLOCATE( hicif_ini(jpi,jpj) , hsnif_ini(jpi,jpj) , frld_ini(jpi,jpj) , sist_ini(jpi,jpj) , zswitch(jpi,jpj) , tbif_ini(jpi,jpj,3) , Stat=ierr(1) )
	941	alloc_lim_istate_init = MAXVAL(ierr)
	942	IF( alloc_lim_istate_init /= 0 ) CALL ctl_warn( 'lim_istate_init: failed to allocate arrays')
	943
	944	END FUNCTION alloc_lim_istate_init
[825]	945	#else
	946	!!----------------------------------------------------------------------
	947	!! Default option : Empty module NO LIM sea-ice model
	948	!!----------------------------------------------------------------------
	949	CONTAINS
	950	SUBROUTINE lim_istate ! Empty routine
	951	END SUBROUTINE lim_istate
	952	#endif
	953
	954	!!======================================================================
	955	END MODULE limistate

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

Download in other formats: