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.

iceistate.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90 @ 8518

Last change on this file since 8518 was 8518, checked in by clem, 7 years ago
changes in style - part6 - commits of the day
File size: 29.6 KB

Rev	Line
[8514]	1	MODULE iceistate
	2	!!======================================================================
	3	!! * MODULE iceistate *
	4	!! Initialisation of diagnostics ice variables
	5	!!======================================================================
	6	!! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code
	7	!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
	8	!! 3.0 ! 2009-11 (M. Vancoppenolle) Enhanced version for ice cats
	9	!! 3.0 ! 2011-02 (G. Madec) dynamical allocation
	10	!! - ! 2014 (C. Rousset) add N/S initializations
	11	!!----------------------------------------------------------------------
	12	#if defined key_lim3
	13	!!----------------------------------------------------------------------
	14	!! 'key_lim3' LIM3 sea-ice model
	15	!!----------------------------------------------------------------------
	16	!! ice_istate : initialization of diagnostics ice variables
	17	!! ice_istate_init : initialization of ice state and namelist read
	18	!!----------------------------------------------------------------------
[8518]	19	USE par_oce ! ocean parameters
[8514]	20	USE phycst ! physical constant
	21	USE oce ! dynamics and tracers variables
	22	USE dom_oce ! ocean domain
[8518]	23	USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
	24	USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
[8514]	25	USE eosbn2 ! equation of state
[8518]	26	USE domvvl ! Variable volume
[8514]	27	USE ice ! sea-ice variables
	28	USE icevar ! ice_var_salprof
	29	!
	30	USE in_out_manager ! I/O manager
	31	USE lib_mpp ! MPP library
	32	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
	33	USE fldread ! read input fields
	34	USE iom
	35
	36	IMPLICIT NONE
	37	PRIVATE
	38
[8515]	39	PUBLIC ice_istate ! called by icestp.F90
	40	PUBLIC ice_istate_init ! called by icestp.F90
[8514]	41
	42	INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read
	43	INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) at T-point
	44	INTEGER , PARAMETER :: jp_hts = 2 ! index of snow thicknes (m) at T-point
	45	INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (%) at T-point
	46	INTEGER , PARAMETER :: jp_tsu = 4 ! index of ice surface temp (K) at T-point
	47	INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temp at T-point
	48	INTEGER , PARAMETER :: jp_smi = 6 ! index of ice sali at T-point
	49	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	50	!
	51	! namelist (namice_ini)
	52	LOGICAL :: ln_iceini ! initialization or not
	53	LOGICAL :: ln_iceini_file ! Ice initialization state from 2D netcdf file
	54	REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice
	55	REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north
	56	REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south
	57	REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north
	58	REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south
	59	REAL(wp) :: rn_ati_ini_n ! initial leads area in the north
	60	REAL(wp) :: rn_ati_ini_s ! initial leads area in the south
	61	REAL(wp) :: rn_smi_ini_n ! initial salinity
	62	REAL(wp) :: rn_smi_ini_s ! initial salinity
	63	REAL(wp) :: rn_tmi_ini_n ! initial temperature
	64	REAL(wp) :: rn_tmi_ini_s ! initial temperature
	65
	66	!!----------------------------------------------------------------------
	67	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
	68	!! $Id: iceistate.F90 8378 2017-07-26 13:55:59Z clem $
	69	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	70	!!----------------------------------------------------------------------
	71	CONTAINS
	72
	73	SUBROUTINE ice_istate
	74	!!-------------------------------------------------------------------
	75	!! * ROUTINE ice_istate *
	76	!!
	77	!! ** Purpose : defined the sea-ice initial state
	78	!!
	79	!! ** Method : This routine will put some ice where ocean
	80	!! is at the freezing point, then fill in ice
	81	!! state variables using prescribed initial
	82	!! values in the namelist
	83	!!
	84	!! ** Steps : 1) Read namelist
	85	!! 2) Basal temperature; ice and hemisphere masks
	86	!! 3) Fill in the ice thickness distribution using gaussian
	87	!! 4) Fill in space-dependent arrays for state variables
	88	!! 5) Diagnostic arrays
	89	!! 6) Lateral boundary conditions
	90	!!
	91	!! ** Notes : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
	92	!! where there is no ice (clem: I do not know why, is it mandatory?)
	93	!!--------------------------------------------------------------------
[8518]	94	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	95	REAL(wp) :: ztmelts, zdh
	96	INTEGER :: i_hemis, i_fill, jl0
	97	REAL(wp) :: zarg, zV, zconv, zdv
	98	INTEGER , DIMENSION(4) :: itest
	99	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	100	REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
	101	REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file
	102	REAL(wp), DIMENSION(jpi,jpj) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
	103	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini !data by cattegories to fill
[8514]	104	!--------------------------------------------------------------------
	105
	106	IF(lwp) WRITE(numout,*)
[8515]	107	IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
	108	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[8514]	109
	110	!--------------------------------------------------------------------
	111	! 1) Read namelist
	112	!--------------------------------------------------------------------
	113	!
	114	! init surface temperature
	115	DO jl = 1, jpl
	116	t_su (:,:,jl) = rt0 * tmask(:,:,1)
	117	tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
	118	END DO
	119
	120	! init basal temperature (considered at freezing point)
	121	CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
	122	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
	123
	124	!--------------------------------------------------------------------
	125	! 2) Initialization of sea ice state variables
	126	!--------------------------------------------------------------------
	127	IF( ln_iceini ) THEN
	128	!
	129	IF( ln_iceini_file )THEN
	130	!
	131	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
	132	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
	133	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
	134	zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
	135	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
	136	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
	137	!
	138	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	139	ELSEWHERE ; zswitch(:,:) = 0._wp
	140	END WHERE
	141	!
	142	ELSE ! ln_iceini_file = F
	143
	144	!--------------------------------------------------------------------
	145	! 3) Basal temperature, ice mask
	146	!--------------------------------------------------------------------
	147	! no ice if sst <= t-freez + ttest
	148	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	149	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
	150	END WHERE
	151
	152	!-----------------------------
	153	! 3.1) Hemisphere-dependent arrays
	154	!-----------------------------
	155	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
[8515]	156	WHERE( ff_t(:,:) >= 0._wp )
	157	zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
	158	zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
	159	zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
	160	zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
	161	zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
	162	ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
	163	ELSEWHERE
	164	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
	165	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
	166	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
	167	zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
	168	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
	169	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
	170	END WHERE
[8514]	171	!
	172	ENDIF ! ln_iceini_file
	173
	174	zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume
	175	!---------------------------------------------------------------------
	176	! 3.2) Distribute ice concentration and thickness into the categories
	177	!---------------------------------------------------------------------
	178	! a gaussian distribution for ice concentration is used
	179	! then we check whether the distribution fullfills
	180	! volume and area conservation, positivity and ice categories bounds
	181	zh_i_ini(:,:,:) = 0._wp
	182	za_i_ini(:,:,:) = 0._wp
	183	!
	184	DO jj = 1, jpj
	185	DO ji = 1, jpi
	186	!
	187	IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
	188
	189	!--- jl0: most likely index where cc will be maximum
	190	jl0 = jpl
	191	DO jl = 1, jpl
	192	IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
	193	jl0 = jl
	194	CYCLE
	195	ENDIF
	196	END DO
	197	!
	198	! initialisation of tests
	199	itest(:) = 0
	200
	201	i_fill = jpl + 1 !====================================
	202	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	203	! iteration !====================================
	204	i_fill = i_fill - 1
	205
	206	! initialisation of ice variables for each try
	207	zh_i_ini(ji,jj,:) = 0._wp
	208	za_i_ini(ji,jj,:) = 0._wp
	209	itest(:) = 0
	210	!
	211	! *** case very thin ice: fill only category 1
	212	IF ( i_fill == 1 ) THEN
	213	zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
	214	za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
	215
	216	! *** case ice is thicker: fill categories >1
	217	ELSE
	218
	219	! Fill ice thicknesses in the (i_fill-1) cat by hmean
	220	DO jl = 1, i_fill-1
	221	zh_i_ini(ji,jj,jl) = hi_mean(jl)
	222	END DO
	223	!
	224	!--- Concentrations
	225	za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
	226	DO jl = 1, i_fill - 1
	227	IF( jl /= jl0 )THEN
	228	zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
	229	za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
	230	ENDIF
	231	END DO
	232	!
	233	! Concentration in the last (i_fill) category
	234	za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
	235
	236	! Ice thickness in the last (i_fill) category
	237	zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
	238	zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )
	239
	240	! clem: correction if concentration of upper cat is greater than lower cat
	241	! (it should be a gaussian around jl0 but sometimes it is not)
	242	IF ( jl0 /= jpl ) THEN
	243	DO jl = jpl, jl0+1, -1
	244	IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
	245	zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
	246	zh_i_ini(ji,jj,jl ) = 0._wp
	247	za_i_ini(ji,jj,jl ) = 0._wp
	248	za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) &
	249	& + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
	250	END IF
	251	ENDDO
	252	ENDIF
	253	!
	254	ENDIF ! case ice is thick or thin
	255
	256	!---------------------
	257	! Compatibility tests
	258	!---------------------
	259	! Test 1: area conservation
	260	zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
	261	IF ( zconv < epsi06 ) itest(1) = 1
	262
	263	! Test 2: volume conservation
	264	zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) &
	265	& - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
	266	IF ( zconv < epsi06 ) itest(2) = 1
	267
	268	! Test 3: thickness of the last category is in-bounds ?
	269	IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
	270
	271	! Test 4: positivity of ice concentrations
	272	itest(4) = 1
	273	DO jl = 1, i_fill
	274	IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
	275	END DO
	276	! !============================
	277	END DO ! end iteration on categories
	278	! !============================
	279	!
	280	IF( lwp .AND. SUM(itest) /= 4 ) THEN
	281	WRITE(numout,*)
	282	WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! '
	283	WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
	284	WRITE(numout,*)
	285	WRITE(numout,) ' ** itest_i (i=1,4) = ', itest(:)
	286	WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	287	WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
	288	ENDIF
	289
	290	ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
	291	!
	292	END DO
	293	END DO
	294
	295	!---------------------------------------------------------------------
	296	! 3.3) Space-dependent arrays for ice state variables
	297	!---------------------------------------------------------------------
	298
	299	! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
	300	DO jl = 1, jpl ! loop over categories
	301	DO jj = 1, jpj
	302	DO ji = 1, jpi
	303	a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration
	304	ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness
	305	sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity
	306	o_i(ji,jj,jl) = 0._wp ! age (0 day)
	307	t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
	308
	309	IF( zht_i_ini(ji,jj) > 0._wp )THEN
	310	ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth
	311	ELSE
	312	ht_s(ji,jj,jl)= 0._wp
	313	ENDIF
	314
	315	! This case below should not be used if (ht_s/ht_i) is ok in namelist
	316	! In case snow load is in excess that would lead to transformation from snow to ice
	317	! Then, transfer the snow excess into the ice (different from icethd_dh)
	318	zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
	319	! recompute ht_i, ht_s avoiding out of bounds values
	320	ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
	321	ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
	322
	323	! ice volume, salt content, age content
	324	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume
	325	v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume
	326	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
	327	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content
	328	END DO
	329	END DO
	330	END DO
	331
	332	! for constant salinity in time
	333	IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
	334	CALL ice_var_salprof
	335	smv_i = sm_i * v_i
	336	ENDIF
	337
	338	! Snow temperature and heat content
	339	DO jk = 1, nlay_s
	340	DO jl = 1, jpl ! loop over categories
	341	DO jj = 1, jpj
	342	DO ji = 1, jpi
	343	t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	344	! Snow energy of melting
	345	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
	346
	347	! Mutliply by volume, and divide by number of layers to get heat content in J/m2
	348	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
	349	END DO
	350	END DO
	351	END DO
	352	END DO
	353
	354	! Ice salinity, temperature and heat content
	355	DO jk = 1, nlay_i
	356	DO jl = 1, jpl ! loop over categories
	357	DO jj = 1, jpj
	358	DO ji = 1, jpi
	359	t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	360	s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
	361	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
	362
	363	! heat content per unit volume
	364	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
	365	+ lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
	366	- rcp * ( ztmelts - rt0 ) )
	367
	368	! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
	369	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
	370	END DO
	371	END DO
	372	END DO
	373	END DO
	374
	375	tn_ice (:,:,:) = t_su (:,:,:)
	376
	377	! MV MP 2016
	378	! Melt pond volume and fraction
	379	IF ( ln_pnd ) THEN
	380	DO jl = 1, jpl
	381	a_ip_frac(:,:,jl) = 0.2 * zswitch(:,:)
	382	h_ip (:,:,jl) = 0.05 * zswitch(:,:)
	383	a_ip(:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)
	384	v_ip(:,:,jl) = h_ip (:,:,jl) * a_ip(:,:,jl)
	385	END DO
	386	ELSE
	387	a_ip(:,:,:) = 0._wp
	388	v_ip(:,:,:) = 0._wp
	389	a_ip_frac(:,:,:) = 0._wp
	390	h_ip (:,:,:) = 0._wp
	391	ENDIF
	392	! END MV MP 2016
	393
	394	ELSE ! if ln_iceini=false
	395	a_i (:,:,:) = 0._wp
	396	v_i (:,:,:) = 0._wp
	397	v_s (:,:,:) = 0._wp
	398	smv_i(:,:,:) = 0._wp
	399	oa_i (:,:,:) = 0._wp
	400	ht_i (:,:,:) = 0._wp
	401	ht_s (:,:,:) = 0._wp
	402	sm_i (:,:,:) = 0._wp
	403	o_i (:,:,:) = 0._wp
	404
	405	e_i(:,:,:,:) = 0._wp
	406	e_s(:,:,:,:) = 0._wp
	407
	408	DO jl = 1, jpl
	409	DO jk = 1, nlay_i
	410	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
	411	END DO
	412	DO jk = 1, nlay_s
	413	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
	414	END DO
	415	END DO
	416
	417	a_ip(:,:,:) = 0._wp
	418	v_ip(:,:,:) = 0._wp
	419	a_ip_frac(:,:,:) = 0._wp
	420	h_ip (:,:,:) = 0._wp
	421
	422	ENDIF ! ln_iceini
	423
	424	at_i (:,:) = 0.0_wp
	425	DO jl = 1, jpl
	426	at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
	427	END DO
	428
	429	!--------------------------------------------------------------------
	430	! 4) Global ice variables for output diagnostics \|
	431	!--------------------------------------------------------------------
	432	u_ice (:,:) = 0._wp
	433	v_ice (:,:) = 0._wp
[8518]	434	!
[8514]	435	!--------------------------------------------------------------------
[8518]	436	! Snow-ice mass (case ice is fully embedded) \|
[8514]	437	!--------------------------------------------------------------------
[8518]	438	snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) ! snow+ice mass
	439	snwice_mass_b(:,:) = snwice_mass(:,:)
	440	!
	441	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[8514]	442
[8518]	443	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	444	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[8514]	445
[8518]	446	IF( .NOT.ln_linssh ) THEN
	447
	448	WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
	449	ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
	450
	451	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	452	e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
	453	e3t_b(:,:,jk) = e3t_n(:,:,jk)
	454	e3t_a(:,:,jk) = e3t_n(:,:,jk)
	455	END DO
	456
	457	! Reconstruction of all vertical scale factors at now and before time-steps
	458	! =========================================================================
	459	! Horizontal scale factor interpolations
	460	! --------------------------------------
	461	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	462	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	463	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	464	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	465	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	466	! Vertical scale factor interpolations
	467	! ------------------------------------
	468	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	469	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	470	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	471	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	472	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	473	! t- and w- points depth
	474	! ----------------------
	475	!!gm not sure of that....
	476	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	477	gdepw_n(:,:,1) = 0.0_wp
	478	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	479	DO jk = 2, jpk
	480	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk )
	481	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	482	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	483	END DO
	484	ENDIF
	485	ENDIF
	486
[8514]	487	!------------------------------------
	488	! 6) store fields at before time-step
	489	!------------------------------------
	490	! it is only necessary for the 1st interpolation by Agrif
	491	a_i_b (:,:,:) = a_i (:,:,:)
	492	e_i_b (:,:,:,:) = e_i (:,:,:,:)
	493	v_i_b (:,:,:) = v_i (:,:,:)
	494	v_s_b (:,:,:) = v_s (:,:,:)
	495	e_s_b (:,:,:,:) = e_s (:,:,:,:)
	496	smv_i_b(:,:,:) = smv_i(:,:,:)
	497	oa_i_b (:,:,:) = oa_i (:,:,:)
	498	u_ice_b(:,:) = u_ice(:,:)
	499	v_ice_b(:,:) = v_ice(:,:)
	500
	501	!!!clem
	502	!! ! Output the initial state and forcings
	503	!! CALL dia_wri_state( 'output.init', nit000 )
	504	!!!
	505
	506	END SUBROUTINE ice_istate
	507
	508	SUBROUTINE ice_istate_init
	509	!!-------------------------------------------------------------------
	510	!! * ROUTINE ice_istate_init *
	511	!!
	512	!! ** Purpose : Definition of initial state of the ice
	513	!!
	514	!! ** Method : Read the namice_ini namelist and check the parameter
	515	!! values called at the first timestep (nit000)
	516	!!
	517	!! ** input :
	518	!! Namelist namice_ini
	519	!!
	520	!! history :
	521	!! 8.5 ! 03-08 (C. Ethe) original code
	522	!! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization
	523	!!-----------------------------------------------------------------------------
	524	!
	525	INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read
	526	INTEGER :: ji,jj
	527	INTEGER :: ifpr, ierror
	528	!
	529	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
	530	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
	531	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	532	!
	533	NAMELIST/namice_ini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, &
	534	& rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
	535	& rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, &
	536	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
	537	!!-----------------------------------------------------------------------------
	538	!
	539	REWIND( numnam_ice_ref ) ! Namelist namice_ini in reference namelist : Ice initial state
	540	READ ( numnam_ice_ref, namice_ini, IOSTAT = ios, ERR = 901)
	541	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in reference namelist', lwp )
	542
	543	REWIND( numnam_ice_cfg ) ! Namelist namice_ini in configuration namelist : Ice initial state
	544	READ ( numnam_ice_cfg, namice_ini, IOSTAT = ios, ERR = 902 )
	545	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in configuration namelist', lwp )
	546	IF(lwm) WRITE ( numoni, namice_ini )
	547
	548	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
	549	slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu
	550	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi
[8515]	551	!
	552	!
	553	IF(lwp) THEN ! control print
[8514]	554	WRITE(numout,*)
	555	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
	556	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	557	WRITE(numout,*) ' Namelist namice_ini'
[8515]	558	WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini
	559	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
	560	WRITE(numout,*) ' max delta ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
	561	WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n
	562	WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s
	563	WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n
	564	WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s
	565	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n
	566	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s
	567	WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n
	568	WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s
	569	WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n
	570	WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s
[8514]	571	ENDIF
	572
	573	IF( ln_iceini_file ) THEN ! Ice initialization using input file
	574	!
	575	! set si structure
	576	ALLOCATE( si(jpfldi), STAT=ierror )
	577	IF( ierror > 0 ) THEN
	578	CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' ) ; RETURN
	579	ENDIF
	580
	581	DO ifpr = 1, jpfldi
	582	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
	583	ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
	584	END DO
	585
	586	! fill si with slf_i and control print
	587	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' )
	588
	589	CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step
	590
	591	ENDIF
	592
	593	END SUBROUTINE ice_istate_init
	594
	595	#else
	596	!!----------------------------------------------------------------------
	597	!! Default option : Empty module NO LIM sea-ice model
	598	!!----------------------------------------------------------------------
	599	#endif
	600
	601	!!======================================================================
	602	END MODULE iceistate

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

Download in other formats: