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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 2633

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

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

Download in other formats: