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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 3544

Last change on this file since 3544 was 3349, checked in by rblod, 12 years ago
Fix unitialized variables in LIM3, see ticket #935
Property svn:keywords set to `Id`
File size: 25.4 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
[2715]	7	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[2528]	8	!!----------------------------------------------------------------------
[825]	9	#if defined key_lim3
	10	!!----------------------------------------------------------------------
[834]	11	!! 'key_lim3' : LIM3 sea-ice model
[825]	12	!!----------------------------------------------------------------------
	13	!! lim_istate : Initialisation of diagnostics ice variables
	14	!! lim_istate_init : initialization of ice state and namelist read
	15	!!----------------------------------------------------------------------
[2528]	16	USE phycst ! physical constant
	17	USE oce ! dynamics and tracers variables
	18	USE dom_oce ! ocean domain
	19	USE sbc_oce ! Surface boundary condition: ocean fields
	20	USE eosbn2 ! equation of state
	21	USE ice ! sea-ice variables
	22	USE par_ice ! ice parameters
	23	USE dom_ice ! sea-ice domain
	24	USE in_out_manager ! I/O manager
	25	USE lbclnk ! lateral boundary condition - MPP exchanges
[2715]	26	USE lib_mpp ! MPP library
[3294]	27	USE wrk_nemo ! work arrays
[825]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC lim_istate ! routine called by lim_init.F90
[825]	33
[2528]	34	! !! init namelist (namiceini)
	35	REAL(wp) :: ttest = 2.0_wp ! threshold water temperature for initial sea ice
	36	REAL(wp) :: hninn = 0.5_wp ! initial snow thickness in the north
	37	REAL(wp) :: hginn_u = 2.5_wp ! initial ice thickness in the north
	38	REAL(wp) :: aginn_u = 0.7_wp ! initial leads area in the north
	39	REAL(wp) :: hginn_d = 5.0_wp ! initial ice thickness in the north
	40	REAL(wp) :: aginn_d = 0.25_wp ! initial leads area in the north
	41	REAL(wp) :: hnins = 0.1_wp ! initial snow thickness in the south
	42	REAL(wp) :: hgins_u = 1.0_wp ! initial ice thickness in the south
	43	REAL(wp) :: agins_u = 0.7_wp ! initial leads area in the south
	44	REAL(wp) :: hgins_d = 2.0_wp ! initial ice thickness in the south
	45	REAL(wp) :: agins_d = 0.2_wp ! initial leads area in the south
	46	REAL(wp) :: sinn = 6.301_wp ! initial salinity
	47	REAL(wp) :: sins = 6.301_wp !
[825]	48
	49	!!----------------------------------------------------------------------
[2715]	50	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	51	!! $Id$
[2528]	52	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	53	!!----------------------------------------------------------------------
	54	CONTAINS
	55
	56	SUBROUTINE lim_istate
	57	!!-------------------------------------------------------------------
	58	!! * ROUTINE lim_istate *
	59	!!
	60	!! ** Purpose : defined the sea-ice initial state
	61	!!
	62	!! ** Method : restart from a state defined in a binary file
	63	!! or from arbitrary sea-ice conditions
[2528]	64	!!-------------------------------------------------------------------
[834]	65	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[2528]	66	REAL(wp) :: zeps6, zeps, ztmelts, epsi06 ! local scalars
[2715]	67	REAL(wp) :: zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs
[3294]	68	REAL(wp), POINTER, DIMENSION(:) :: zgfactorn, zhin
	69	REAL(wp), POINTER, DIMENSION(:) :: zgfactors, zhis
	70	REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator
[825]	71	!--------------------------------------------------------------------
[921]	72
[3294]	73	CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis )
	74	CALL wrk_alloc( jpi, jpj, zidto )
[2715]	75
[825]	76	!--------------------------------------------------------------------
	77	! 1) Preliminary things
	78	!--------------------------------------------------------------------
[2528]	79	epsi06 = 1.e-6_wp
[825]	80
	81	CALL lim_istate_init ! reading the initials parameters of the ice
	82
[1037]	83	!!gm in lim2 the initialisation if only done if required in the namelist :
	84	!!gm IF( .NOT. ln_limini ) THEN
	85	!!gm this should be added in lim3 namelist...
[825]	86
	87	!--------------------------------------------------------------------
	88	! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s) \|
	89	!--------------------------------------------------------------------
	90
[1037]	91	IF(lwp) WRITE(numout,*)
	92	IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
	93	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
[825]	94
[3294]	95	t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius]
[825]	96
[1037]	97	DO jj = 1, jpj ! ice if sst <= t-freez + ttest
	98	DO ji = 1, jpi
[3294]	99	IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0.e0 ! no ice
	100	ELSE ; zidto(ji,jj) = 1.e0 ! ice
[1037]	101	ENDIF
	102	END DO
	103	END DO
	104
	105	t_bo(:,:) = t_bo(:,:) + rt0 ! t_bo converted from Celsius to Kelvin (rt0 over land)
	106
[825]	107	! constants for heat contents
[2528]	108	zeps = 1.e-20_wp
	109	zeps6 = 1.e-06_wp
[825]	110
	111	! zgfactor for initial ice distribution
[2528]	112	zgfactorn(:) = 0._wp
	113	zgfactors(:) = 0._wp
[825]	114
	115	! first ice type
	116	DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
[2528]	117	zhin (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
	118	zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)(zhin(1)-hginn_u) 0.5_wp )
	119	zhis (1) = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
	120	zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)(zhis(1)-hgins_u) 0.5_wp )
[825]	121	END DO ! jl
	122	zgfactorn(1) = aginn_u / zgfactorn(1)
	123	zgfactors(1) = agins_u / zgfactors(1)
	124
	125	! -------------
	126	! new distribution, polynom of second order, conserving area and volume
[2528]	127	zh1 = 0._wp
	128	zh2 = 0._wp
	129	zh3 = 0._wp
[825]	130	DO jl = 1, jpl
[2528]	131	zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
[825]	132	zh1 = zh1 + zh
[2528]	133	zh2 = zh2 + zh * zh
	134	zh3 = zh3 + zh * zh * zh
[825]	135	END DO
[1037]	136	IF(lwp) WRITE(numout,*) ' zh1 : ', zh1
	137	IF(lwp) WRITE(numout,*) ' zh2 : ', zh2
	138	IF(lwp) WRITE(numout,*) ' zh3 : ', zh3
[825]	139
[2528]	140	zvol = aginn_u * hginn_u
[825]	141	zare = aginn_u
[2528]	142	IF( jpl >= 2 ) THEN
[825]	143	zbn = ( zvolzh2 - zarezh3 ) / ( zh2zh2 - zh1zh3)
	144	zan = ( zare - zbn*zh1 ) / zh2
	145	ENDIF
	146
[1037]	147	IF(lwp) WRITE(numout,*) ' zvol: ', zvol
	148	IF(lwp) WRITE(numout,*) ' zare: ', zare
	149	IF(lwp) WRITE(numout,*) ' zbn : ', zbn
	150	IF(lwp) WRITE(numout,*) ' zan : ', zan
[825]	151
[2528]	152	zvol = agins_u * hgins_u
[825]	153	zare = agins_u
[2528]	154	IF( jpl >= 2 ) THEN
[825]	155	zbs = ( zvolzh2 - zarezh3 ) / ( zh2zh2 - zh1zh3)
	156	zas = ( zare - zbs*zh1 ) / zh2
	157	ENDIF
	158
[1037]	159	IF(lwp) WRITE(numout,*) ' zvol: ', zvol
	160	IF(lwp) WRITE(numout,*) ' zare: ', zare
	161	IF(lwp) WRITE(numout,*) ' zbn : ', zbn
	162	IF(lwp) WRITE(numout,*) ' zan : ', zan
[825]	163
	164	!end of new lines
	165	! -------------
	166	!!!
[921]	167	! retour a LIMA_MEC
	168	! ! second ice type
	169	! zdummy = hi_max(ice_cat_bounds(2,1)-1)
	170	! hi_max(ice_cat_bounds(2,1)-1) = 0.0
[825]	171
[921]	172	! ! here to change !!!!
	173	! jm = 2
	174	! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	175	! zhin (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
	176	! zhin (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + &
	177	! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0
	178	! zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0)
	179	! zhis (2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
	180	! zhis (2) = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm ) + &
	181	! hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm) ) / 2.0
	182	! zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0)
	183	! END DO ! jl
	184	! zgfactorn(2) = aginn_d / zgfactorn(2)
	185	! zgfactors(2) = agins_d / zgfactors(2)
	186	! hi_max(ice_cat_bounds(2,1)-1) = zdummy
	187	! END retour a LIMA_MEC
[825]	188	!!!
[1037]	189
	190	!!gm optimisation : loop over the ice categories inside the ji, jj loop !!!
	191
[825]	192	DO jj = 1, jpj
	193	DO ji = 1, jpi
	194
	195	!--- Northern hemisphere
	196	!----------------------------------------------------------------
[2528]	197	IF( fcor(ji,jj) >= 0._wp ) THEN
[825]	198
	199	!-----------------------
	200	! Ice area / thickness
	201	!-----------------------
	202
	203	IF ( jpl .EQ. 1) THEN ! one category
	204
	205	DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
[1037]	206	a_i(ji,jj,jl) = zidto(ji,jj) * aginn_u
	207	ht_i(ji,jj,jl) = zidto(ji,jj) * hginn_u
[825]	208	v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
	209	END DO
	210
	211	ELSE ! several categories
	212
	213	DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
	214	zhin(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
[1037]	215	a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* &
[921]	216	(zhin(1)-hginn_u)/2.0) , epsi06)
[825]	217	! new line
[1037]	218	a_i(ji,jj,jl) = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) )
	219	ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(1)
[825]	220	v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
	221	END DO
	222
	223	ENDIF
	224
	225
	226	!!!
[921]	227	! retour a LIMA_MEC
	228	! !ridged ice
	229	! zdummy = hi_max(ice_cat_bounds(2,1)-1)
	230	! hi_max(ice_cat_bounds(2,1)-1) = 0.0
	231	! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories
	232	! zhin(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
[1037]	233	! a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* &
[921]	234	! (zhin(2)-hginn_d)/2.0) , epsi06)
[1037]	235	! ht_i(ji,jj,jl) = zidto(ji,jj) * zhin(2)
[921]	236	! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
	237	! END DO
	238	! hi_max(ice_cat_bounds(2,1)-1) = zdummy
[825]	239
[921]	240	! !rafted ice
	241	! jl = 6
	242	! a_i(ji,jj,jl) = 0.0
	243	! ht_i(ji,jj,jl) = 0.0
	244	! v_i(ji,jj,jl) = 0.0
	245	! END retour a LIMA_MEC
[825]	246	!!!
	247
	248	DO jl = 1, jpl
	249
	250	!-------------
	251	! Snow depth
	252	!-------------
[1037]	253	ht_s(ji,jj,jl) = zidto(ji,jj) * hninn
[825]	254	v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
	255
	256	!---------------
	257	! Ice salinity
	258	!---------------
[1037]	259	sm_i(ji,jj,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
[888]	260	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
[825]	261
	262	!----------
	263	! Ice age
	264	!----------
[1037]	265	o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) )
[825]	266	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)
[921]	267
[825]	268	!------------------------------
	269	! Sea ice surface temperature
	270	!------------------------------
	271
[1037]	272	t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
[825]	273
	274	!------------------------------------
	275	! Snow temperature and heat content
	276	!------------------------------------
	277
	278	DO jk = 1, nlay_s
[1037]	279	t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
[825]	280	! Snow energy of melting
[1103]	281	e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
[825]	282	! Change dimensions
	283	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
	284	! Multiply by volume, so that heat content in 10^9 Joules
	285	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
[921]	286	v_s(ji,jj,jl) / nlay_s
[825]	287	END DO !jk
	288
	289	!-----------------------------------------------
	290	! Ice salinities, temperature and heat content
	291	!-----------------------------------------------
	292
	293	DO jk = 1, nlay_i
[1037]	294	t_i(ji,jj,jk,jl) = zidto(ji,jj)270.00 + ( 1.0 - zidto(ji,jj) ) rtt
	295	s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
[825]	296	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
[921]	297
	298	! heat content per unit volume
[1037]	299	e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
[921]	300	( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
	301	+ lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
	302	- rcp * ( ztmelts - rtt ) &
	303	)
[825]	304
[921]	305	! Correct dimensions to avoid big values
[825]	306	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
	307
[921]	308	! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
[825]	309	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
[921]	310	area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
	311	nlay_i
[825]	312	END DO ! jk
	313
	314	END DO ! jl
	315
	316	ELSE ! on fcor
	317
[921]	318	!--- Southern hemisphere
	319	!----------------------------------------------------------------
[825]	320
	321	!-----------------------
	322	! Ice area / thickness
	323	!-----------------------
	324
	325	IF ( jpl .EQ. 1) THEN ! one category
	326
	327	DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
[1037]	328	a_i(ji,jj,jl) = zidto(ji,jj) * agins_u
	329	ht_i(ji,jj,jl) = zidto(ji,jj) * hgins_u
[825]	330	v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
	331	END DO
	332
	333	ELSE ! several categories
[921]	334
	335	!level ice
[825]	336	DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories
	337
	338	zhis(1) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
[1037]	339	a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * &
[921]	340	(zhis(1)-hgins_u)/2.0) , epsi06 )
[825]	341	! new line square distribution volume conserving
[1037]	342	a_i(ji,jj,jl) = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) )
	343	ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(1)
[825]	344	v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
[921]	345
[825]	346	END DO ! jl
	347
	348	ENDIF
	349
	350	!!!
[921]	351	! retour a LIMA_MEC
	352	! !ridged ice
	353	! zdummy = hi_max(ice_cat_bounds(2,1)-1)
	354	! hi_max(ice_cat_bounds(2,1)-1) = 0.0
	355	! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories
	356	! zhis(2) = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
[1037]	357	! a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2) &
	358	! & * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 )
	359	! ht_i(ji,jj,jl) = zidto(ji,jj) * zhis(2)
[921]	360	! v_i(ji,jj,jl) = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
	361	! END DO
	362	! hi_max(ice_cat_bounds(2,1)-1) = zdummy
[825]	363
[921]	364	! !rafted ice
	365	! jl = 6
	366	! a_i(ji,jj,jl) = 0.0
	367	! ht_i(ji,jj,jl) = 0.0
	368	! v_i(ji,jj,jl) = 0.0
	369	! END retour a LIMA_MEC
[825]	370	!!!
	371
	372	DO jl = 1, jpl !over thickness categories
	373
	374	!---------------
	375	! Snow depth
	376	!---------------
	377
[1037]	378	ht_s(ji,jj,jl) = zidto(ji,jj) * hnins
[825]	379	v_s(ji,jj,jl) = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
	380
	381	!---------------
	382	! Ice salinity
	383	!---------------
	384
[1037]	385	sm_i(ji,jj,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
[888]	386	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
[825]	387
	388	!----------
	389	! Ice age
	390	!----------
	391
[1037]	392	o_i(ji,jj,jl) = zidto(ji,jj) * 1.0 + ( 1.0 - zidto(ji,jj) )
[825]	393	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)
	394
	395	!------------------------------
	396	! Sea ice surface temperature
	397	!------------------------------
	398
[1037]	399	t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
[825]	400
	401	!----------------------------------
	402	! Snow temperature / heat content
	403	!----------------------------------
	404
	405	DO jk = 1, nlay_s
[1037]	406	t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
[825]	407	! Snow energy of melting
[1037]	408	e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
[825]	409	! Change dimensions
	410	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
	411	! Multiply by volume, so that heat content in 10^9 Joules
	412	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
[921]	413	v_s(ji,jj,jl) / nlay_s
[825]	414	END DO
	415
	416	!---------------------------------------------
	417	! Ice temperature, salinity and heat content
	418	!---------------------------------------------
	419
	420	DO jk = 1, nlay_i
[1037]	421	t_i(ji,jj,jk,jl) = zidto(ji,jj)270.00 + ( 1.0 - zidto(ji,jj) ) rtt
	422	s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
[825]	423	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
[921]	424
	425	! heat content per unit volume
[1037]	426	e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
[921]	427	( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
	428	+ lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
	429	- rcp * ( ztmelts - rtt ) &
	430	)
[825]	431
[921]	432	! Correct dimensions to avoid big values
[825]	433	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
	434
[921]	435	! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
[825]	436	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
[921]	437	area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
	438	nlay_i
[825]	439	END DO !jk
	440
	441	END DO ! jl
	442
	443	ENDIF ! on fcor
	444
[2528]	445	END DO
	446	END DO
[825]	447
	448	!--------------------------------------------------------------------
	449	! 3) Global ice variables for output diagnostics \|
	450	!--------------------------------------------------------------------
	451
	452	fsbbq (:,:) = 0.e0
	453	u_ice (:,:) = 0.e0
	454	v_ice (:,:) = 0.e0
	455	stress1_i(:,:) = 0.0
	456	stress2_i(:,:) = 0.0
	457	stress12_i(:,:) = 0.0
	458
	459	!--------------------------------------------------------------------
	460	! 4) Moments for advection
	461	!--------------------------------------------------------------------
	462
[3349]	463	sxopw (:,:) = 0.e0
	464	syopw (:,:) = 0.e0
	465	sxxopw(:,:) = 0.e0
	466	syyopw(:,:) = 0.e0
	467	sxyopw(:,:) = 0.e0
	468
[825]	469	sxice (:,:,:) = 0.e0 ; sxsn (:,:,:) = 0.e0 ; sxa (:,:,:) = 0.e0
	470	syice (:,:,:) = 0.e0 ; sysn (:,:,:) = 0.e0 ; sya (:,:,:) = 0.e0
	471	sxxice(:,:,:) = 0.e0 ; sxxsn(:,:,:) = 0.e0 ; sxxa (:,:,:) = 0.e0
	472	syyice(:,:,:) = 0.e0 ; syysn(:,:,:) = 0.e0 ; syya (:,:,:) = 0.e0
	473	sxyice(:,:,:) = 0.e0 ; sxysn(:,:,:) = 0.e0 ; sxya (:,:,:) = 0.e0
	474
	475	sxc0 (:,:,:) = 0.e0 ; sxe (:,:,:,:)= 0.e0
	476	syc0 (:,:,:) = 0.e0 ; sye (:,:,:,:)= 0.e0
	477	sxxc0 (:,:,:) = 0.e0 ; sxxe (:,:,:,:)= 0.e0
	478	syyc0 (:,:,:) = 0.e0 ; syye (:,:,:,:)= 0.e0
	479	sxyc0 (:,:,:) = 0.e0 ; sxye (:,:,:,:)= 0.e0
	480
	481	sxsal (:,:,:) = 0.e0
	482	sysal (:,:,:) = 0.e0
	483	sxxsal (:,:,:) = 0.e0
	484	syysal (:,:,:) = 0.e0
	485	sxysal (:,:,:) = 0.e0
	486
	487	!--------------------------------------------------------------------
	488	! 5) Lateral boundary conditions \|
	489	!--------------------------------------------------------------------
	490
	491	DO jl = 1, jpl
	492	CALL lbc_lnk( a_i(:,:,jl) , 'T', 1. )
	493	CALL lbc_lnk( v_i(:,:,jl) , 'T', 1. )
	494	CALL lbc_lnk( v_s(:,:,jl) , 'T', 1. )
	495	CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
	496	CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
[2528]	497	!
[825]	498	CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
	499	CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
	500	CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
	501	CALL lbc_lnk( o_i(:,:,jl) , 'T', 1. )
	502	CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
	503	DO jk = 1, nlay_s
	504	CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
[869]	505	CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
[825]	506	END DO
	507	DO jk = 1, nlay_i
	508	CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
	509	CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
	510	END DO
[2528]	511	!
	512	a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
[825]	513	END DO
	514
	515	CALL lbc_lnk( at_i , 'T', 1. )
	516	at_i(:,:) = tms(:,:) * at_i(:,:) ! put 0 over land
[2528]	517	!
[825]	518	CALL lbc_lnk( fsbbq , 'T', 1. )
[2528]	519	!
[3294]	520	CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis )
	521	CALL wrk_dealloc( jpi, jpj, zidto )
[2715]	522	!
[825]	523	END SUBROUTINE lim_istate
	524
[2528]	525
[825]	526	SUBROUTINE lim_istate_init
	527	!!-------------------------------------------------------------------
	528	!! * ROUTINE lim_istate_init *
	529	!!
	530	!! ** Purpose : Definition of initial state of the ice
	531	!!
[2528]	532	!! ** Method : Read the namiceini namelist and check the parameter
	533	!! values called at the first timestep (nit000)
[825]	534	!!
[2528]	535	!! ** input : namelist namiceini
[825]	536	!!-----------------------------------------------------------------------------
[2528]	537	NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins, &
	538	& hgins_u, agins_u, hgins_d, agins_d, sinn, sins
[825]	539	!!-----------------------------------------------------------------------------
[2528]	540	!
	541	REWIND ( numnam_ice ) ! Read Namelist namiceini
[825]	542	READ ( numnam_ice , namiceini )
[2528]	543	!
	544	IF(lwp) THEN ! control print
[825]	545	WRITE(numout,*)
	546	WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
	547	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	548	WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest
	549	WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn
	550	WRITE(numout,*) ' initial undef ice thickness in the north hginn_u = ', hginn_u
	551	WRITE(numout,*) ' initial undef ice concentr. in the north aginn_u = ', aginn_u
	552	WRITE(numout,*) ' initial def ice thickness in the north hginn_d = ', hginn_d
	553	WRITE(numout,*) ' initial def ice concentr. in the north aginn_d = ', aginn_d
	554	WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins
	555	WRITE(numout,*) ' initial undef ice thickness in the north hgins_u = ', hgins_u
	556	WRITE(numout,*) ' initial undef ice concentr. in the north agins_u = ', agins_u
	557	WRITE(numout,*) ' initial def ice thickness in the north hgins_d = ', hgins_d
	558	WRITE(numout,*) ' initial def ice concentr. in the north agins_d = ', agins_d
	559	WRITE(numout,*) ' initial ice salinity in the north sinn = ', sinn
	560	WRITE(numout,*) ' initial ice salinity in the south sins = ', sins
	561	ENDIF
[2528]	562	!
[825]	563	END SUBROUTINE lim_istate_init
	564
	565	#else
	566	!!----------------------------------------------------------------------
	567	!! Default option : Empty module NO LIM sea-ice model
	568	!!----------------------------------------------------------------------
	569	CONTAINS
	570	SUBROUTINE lim_istate ! Empty routine
	571	END SUBROUTINE lim_istate
	572	#endif
	573
	574	!!======================================================================
	575	END MODULE limistate

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

Download in other formats: