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

Context Navigation

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 8179

Last change on this file since 8179 was 8106, checked in by vancop, 7 years ago
Follow up on melt ponds and albedo
Property svn:keywords set to `Id`
File size: 27.8 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
[6469]	25	USE limvar ! lim_var_salprof
[2528]	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
[6515]	30	USE fldread ! read input fields
	31	USE iom
[825]	32
[6853]	33	!!!clem
	34	!! USE diawri
	35	!!!
	36
[825]	37	IMPLICIT NONE
	38	PRIVATE
	39
[2528]	40	PUBLIC lim_istate ! routine called by lim_init.F90
[825]	41
[6515]	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)
[825]	50	!!----------------------------------------------------------------------
[4161]	51	!! LIM 3.0, UCL-LOCEAN-IPSL (2008)
[1156]	52	!! $Id$
[4161]	53	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]	54	!!----------------------------------------------------------------------
	55	CONTAINS
	56
	57	SUBROUTINE lim_istate
	58	!!-------------------------------------------------------------------
	59	!! * ROUTINE lim_istate *
	60	!!
	61	!! ** Purpose : defined the sea-ice initial state
	62	!!
[4161]	63	!! ** Method :
	64	!! This routine will put some ice where ocean
	65	!! is at the freezing point, then fill in ice
	66	!! state variables using prescribed initial
	67	!! values in the namelist
	68	!!
	69	!! ** Steps :
	70	!! 1) Read namelist
	71	!! 2) Basal temperature; ice and hemisphere masks
	72	!! 3) Fill in the ice thickness distribution using gaussian
	73	!! 4) Fill in space-dependent arrays for state variables
	74	!! 5) Diagnostic arrays
	75	!! 6) Lateral boundary conditions
	76	!!
[4333]	77	!! ** Notes : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
[4990]	78	!! where there is no ice (clem: I do not know why, is it mandatory?)
[4333]	79	!!
[4161]	80	!! History :
	81	!! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code
	82	!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
	83	!! 4.0 ! 09-11 (M. Vancoppenolle) Enhanced version for ice cats
	84	!!--------------------------------------------------------------------
	85
	86	!! * Local variables
	87	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[5123]	88	REAL(wp) :: ztmelts, zdh
[8085]	89	INTEGER :: i_hemis, i_fill, jl0
[8098]	90	REAL(wp) :: zarg, zV, zconv, zdv
[4688]	91	REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator
[6515]	92	REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file
	93	REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
[6853]	94	REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini !data by cattegories to fill
	95	INTEGER , POINTER, DIMENSION(:) :: itest
[825]	96	!--------------------------------------------------------------------
[921]	97
[6853]	98	CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini )
[6515]	99	CALL wrk_alloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
	100	CALL wrk_alloc( jpi, jpj, zswitch )
[6853]	101	Call wrk_alloc( 4, itest )
[2715]	102
[4161]	103	IF(lwp) WRITE(numout,*)
	104	IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
	105	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
	106
[825]	107	!--------------------------------------------------------------------
[4161]	108	! 1) Read namelist
[825]	109	!--------------------------------------------------------------------
[7060]	110	CALL lim_istate_init
[825]	111
[7060]	112	! init surface temperature
	113	DO jl = 1, jpl
[5123]	114	t_su (:,:,jl) = rt0 * tmask(:,:,1)
	115	tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
[4688]	116	END DO
	117
[7060]	118	! init basal temperature (considered at freezing point)
[5540]	119	CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
	120	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
[4990]	121
[4688]	122
[7060]	123	!--------------------------------------------------------------------
	124	! 2) Initialization of sea ice state variables
	125	!--------------------------------------------------------------------
[6515]	126	IF( ln_limini ) THEN
[4990]	127
[6515]	128	IF( ln_limini_file )THEN
[1037]	129
[6515]	130	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
	131	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
	132	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
	133	zts_u_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
	134	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
	135	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
[825]	136
[7060]	137	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	138	ELSEWHERE ; zswitch(:,:) = 0._wp
	139	END WHERE
	140
[6515]	141	ELSE ! ln_limini_file = F
[825]	142
[7060]	143	!--------------------------------------------------------------------
	144	! 3) Basal temperature, ice mask
	145	!--------------------------------------------------------------------
	146	! no ice if sst <= t-freez + ttest
	147	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	148	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
	149	END WHERE
	150
[6515]	151	!-----------------------------
	152	! 3.1) Hemisphere-dependent arrays
	153	!-----------------------------
	154	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
	155	DO jj = 1, jpj
	156	DO ji = 1, jpi
[6746]	157	IF( ff(ji,jj) >= 0._wp ) THEN
[6853]	158	zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
	159	zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
	160	zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
	161	zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
	162	zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
	163	ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
[6515]	164	ELSE
[6853]	165	zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
	166	zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
	167	zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
	168	zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
	169	zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
	170	ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
[6515]	171	ENDIF
	172	END DO
	173	END DO
[825]	174
[6515]	175	ENDIF ! ln_limini_file
	176
	177	zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume
	178	!---------------------------------------------------------------------
	179	! 3.2) Distribute ice concentration and thickness into the categories
	180	!---------------------------------------------------------------------
	181	! a gaussian distribution for ice concentration is used
	182	! then we check whether the distribution fullfills
	183	! volume and area conservation, positivity and ice categories bounds
	184	zh_i_ini(:,:,:) = 0._wp
	185	za_i_ini(:,:,:) = 0._wp
[4688]	186
[6515]	187	DO jj = 1, jpj
	188	DO ji = 1, jpi
[825]	189
[6515]	190	IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
[825]	191
[6853]	192	!--- jl0: most likely index where cc will be maximum
	193	jl0 = jpl
	194	DO jl = 1, jpl
	195	IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
	196	jl0 = jl
	197	CYCLE
	198	ENDIF
	199	END DO
[825]	200
[6853]	201	! initialisation of tests
	202	itest(:) = 0
	203
	204	i_fill = jpl + 1 !====================================
	205	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	206	! iteration !====================================
	207	i_fill = i_fill - 1
[825]	208
[6853]	209	! initialisation of ice variables for each try
	210	zh_i_ini(ji,jj,:) = 0._wp
	211	za_i_ini(ji,jj,:) = 0._wp
	212	itest(:) = 0
[6515]	213
[6853]	214	! *** case very thin ice: fill only category 1
	215	IF ( i_fill == 1 ) THEN
	216	zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
	217	za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
	218
	219	! *** case ice is thicker: fill categories >1
	220	ELSE
	221
	222	! Fill ice thicknesses in the (i_fill-1) cat by hmean
	223	DO jl = 1, i_fill-1
	224	zh_i_ini(ji,jj,jl) = hi_mean(jl)
	225	END DO
[4688]	226
[6853]	227	!--- Concentrations
	228	za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
	229	DO jl = 1, i_fill - 1
	230	IF( jl /= jl0 )THEN
	231	zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
	232	za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
	233	ENDIF
	234	END DO
[4688]	235
[6853]	236	! Concentration in the last (i_fill) category
	237	za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
[825]	238
[6853]	239	! Ice thickness in the last (i_fill) category
	240	zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
	241	zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )
[921]	242
[6853]	243	! clem: correction if concentration of upper cat is greater than lower cat
	244	! (it should be a gaussian around jl0 but sometimes it is not)
	245	IF ( jl0 /= jpl ) THEN
	246	DO jl = jpl, jl0+1, -1
	247	IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
	248	zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
	249	zh_i_ini(ji,jj,jl ) = 0._wp
	250	za_i_ini(ji,jj,jl ) = 0._wp
	251	za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) &
	252	& + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
	253	END IF
	254	ENDDO
[6515]	255	ENDIF
[825]	256
[6853]	257	ENDIF ! case ice is thick or thin
[825]	258
[6853]	259	!---------------------
	260	! Compatibility tests
	261	!---------------------
	262	! Test 1: area conservation
	263	zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
	264	IF ( zconv < epsi06 ) itest(1) = 1
	265
	266	! Test 2: volume conservation
	267	zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) &
	268	& - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
	269	IF ( zconv < epsi06 ) itest(2) = 1
	270
	271	! Test 3: thickness of the last category is in-bounds ?
	272	IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
	273
	274	! Test 4: positivity of ice concentrations
	275	itest(4) = 1
	276	DO jl = 1, i_fill
	277	IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
	278	END DO
	279	! !============================
	280	END DO ! end iteration on categories
	281	! !============================
[825]	282
[6853]	283	IF( lwp .AND. SUM(itest) /= 4 ) THEN
	284	WRITE(numout,*)
	285	WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! '
	286	WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
	287	WRITE(numout,*)
	288	WRITE(numout,) ' ** itest_i (i=1,4) = ', itest(:)
	289	WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
	290	WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
[6515]	291	ENDIF
[6853]	292
	293	ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
	294
	295	ENDDO
	296	ENDDO
[825]	297
[6515]	298	!---------------------------------------------------------------------
	299	! 3.3) Space-dependent arrays for ice state variables
	300	!---------------------------------------------------------------------
[4333]	301
[6515]	302	! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
	303	DO jl = 1, jpl ! loop over categories
	304	DO jj = 1, jpj
	305	DO ji = 1, jpi
	306	a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini(ji,jj,jl) ! concentration
	307	ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(ji,jj,jl) ! ice thickness
	308	sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) ! salinity
	309	o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day)
	310	t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
[4333]	311
[6515]	312	IF( zht_i_ini(ji,jj) > 0._wp )THEN
	313	ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) ) ! snow depth
	314	ELSE
	315	ht_s(ji,jj,jl)= 0._wp
	316	ENDIF
	317
	318	! This case below should not be used if (ht_s/ht_i) is ok in namelist
	319	! In case snow load is in excess that would lead to transformation from snow to ice
	320	! Then, transfer the snow excess into the ice (different from limthd_dh)
	321	zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
	322	! recompute ht_i, ht_s avoiding out of bounds values
	323	ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
	324	ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
	325
	326	! ice volume, salt content, age content
	327	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume
	328	v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume
	329	smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
	330	oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content
	331	END DO
[5202]	332	END DO
	333	END DO
[825]	334
[6515]	335	! for constant salinity in time
	336	IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
	337	CALL lim_var_salprof
	338	smv_i = sm_i * v_i
	339	ENDIF
	340
	341	! Snow temperature and heat content
	342	DO jk = 1, nlay_s
	343	DO jl = 1, jpl ! loop over categories
	344	DO jj = 1, jpj
	345	DO ji = 1, jpi
	346	t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	347	! Snow energy of melting
	348	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
[6469]	349
[6515]	350	! Mutliply by volume, and divide by number of layers to get heat content in J/m2
	351	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
	352	END DO
[5202]	353	END DO
	354	END DO
	355	END DO
[825]	356
[6515]	357	! Ice salinity, temperature and heat content
	358	DO jk = 1, nlay_i
	359	DO jl = 1, jpl ! loop over categories
	360	DO jj = 1, jpj
	361	DO ji = 1, jpi
	362	t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
	363	s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
	364	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
[825]	365
[6515]	366	! heat content per unit volume
	367	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &
	368	+ lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
	369	- rcp * ( ztmelts - rt0 ) )
[825]	370
[6515]	371	! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
	372	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
	373	END DO
[5202]	374	END DO
	375	END DO
	376	END DO
[825]	377
[6515]	378	tn_ice (:,:,:) = t_su (:,:,:)
[4688]	379
[8085]	380	! MV MP 2016
	381	! Melt pond volume and fraction
[8098]	382	IF ( ln_pnd ) THEN
	383
	384	DO jl = 1, jpl
[8106]	385
	386	a_ip_frac(:,:,jl) = 0.2 * zswitch(:,:)
	387	h_ip (:,:,jl) = 0.05 * zswitch(:,:)
	388	a_ip(:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)
	389	v_ip(:,:,jl) = h_ip (:,:,jl) * a_ip(:,:,jl)
	390
[8098]	391	END DO
	392
	393	ELSE
	394
[8106]	395	a_ip(:,:,:) = 0._wp
	396	v_ip(:,:,:) = 0._wp
	397	a_ip_frac(:,:,:) = 0._wp
	398	h_ip (:,:,:) = 0._wp
[8098]	399
	400	ENDIF
[8085]	401	! END MV MP 2016
	402
[6515]	403	ELSE ! if ln_limini=false
[4688]	404	a_i (:,:,:) = 0._wp
	405	v_i (:,:,:) = 0._wp
	406	v_s (:,:,:) = 0._wp
	407	smv_i(:,:,:) = 0._wp
	408	oa_i (:,:,:) = 0._wp
	409	ht_i (:,:,:) = 0._wp
	410	ht_s (:,:,:) = 0._wp
	411	sm_i (:,:,:) = 0._wp
	412	o_i (:,:,:) = 0._wp
	413
	414	e_i(:,:,:,:) = 0._wp
	415	e_s(:,:,:,:) = 0._wp
	416
	417	DO jl = 1, jpl
	418	DO jk = 1, nlay_i
[5123]	419	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]	420	END DO
	421	DO jk = 1, nlay_s
[5123]	422	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]	423	END DO
	424	END DO
[6515]	425
[8085]	426	a_ip(:,:,:) = 0._wp
	427	v_ip(:,:,:) = 0._wp
[8106]	428	a_ip_frac(:,:,:) = 0._wp
	429	h_ip (:,:,:) = 0._wp
[8085]	430
[6515]	431	ENDIF ! ln_limini
[4688]	432
	433	at_i (:,:) = 0.0_wp
	434	DO jl = 1, jpl
	435	at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
	436	END DO
[7293]	437
[825]	438	!--------------------------------------------------------------------
[4161]	439	! 4) Global ice variables for output diagnostics \|
[825]	440	!--------------------------------------------------------------------
[4161]	441	u_ice (:,:) = 0._wp
	442	v_ice (:,:) = 0._wp
	443	stress1_i(:,:) = 0._wp
	444	stress2_i(:,:) = 0._wp
	445	stress12_i(:,:) = 0._wp
[825]	446
	447	!--------------------------------------------------------------------
[4161]	448	! 5) Moments for advection
[825]	449	!--------------------------------------------------------------------
	450
[4161]	451	sxopw (:,:) = 0._wp
	452	syopw (:,:) = 0._wp
	453	sxxopw(:,:) = 0._wp
	454	syyopw(:,:) = 0._wp
	455	sxyopw(:,:) = 0._wp
[3610]	456
[4161]	457	sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp
	458	syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp
	459	sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp
	460	syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp
	461	sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp
[825]	462
[4161]	463	sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp
	464	syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp
	465	sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp
	466	syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp
	467	sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp
[825]	468
[4161]	469	sxsal (:,:,:) = 0._wp
	470	sysal (:,:,:) = 0._wp
	471	sxxsal (:,:,:) = 0._wp
	472	syysal (:,:,:) = 0._wp
	473	sxysal (:,:,:) = 0._wp
[825]	474
[4335]	475	sxage (:,:,:) = 0._wp
	476	syage (:,:,:) = 0._wp
	477	sxxage (:,:,:) = 0._wp
	478	syyage (:,:,:) = 0._wp
	479	sxyage (:,:,:) = 0._wp
	480
[7293]	481	! MV MP 2016
[8098]	482	IF ( nn_pnd_scheme > 0 ) THEN
[7306]	483	sxap (:,:,:) = 0._wp ; sxvp (:,:,:) = 0._wp
	484	syap (:,:,:) = 0._wp ; syvp (:,:,:) = 0._wp
	485	sxxap (:,:,:) = 0._wp ; sxxvp (:,:,:) = 0._wp
	486	syyap (:,:,:) = 0._wp ; syyvp (:,:,:) = 0._wp
	487	sxyap (:,:,:) = 0._wp ; sxyvp (:,:,:) = 0._wp
[7293]	488	ENDIF
	489	! END MV MP 2016
	490
[6853]	491	!!!clem
	492	!! ! Output the initial state and forcings
	493	!! CALL dia_wri_state( 'output.init', nit000 )
[7060]	494	!!!
	495
[6853]	496	CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini )
[6515]	497	CALL wrk_dealloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
	498	CALL wrk_dealloc( jpi, jpj, zswitch )
[6853]	499	Call wrk_dealloc( 4, itest )
[4161]	500
[825]	501	END SUBROUTINE lim_istate
	502
	503	SUBROUTINE lim_istate_init
	504	!!-------------------------------------------------------------------
	505	!! * ROUTINE lim_istate_init *
	506	!!
	507	!! ** Purpose : Definition of initial state of the ice
	508	!!
[4161]	509	!! ** Method : Read the namiceini namelist and check the parameter
	510	!! values called at the first timestep (nit000)
[825]	511	!!
[4161]	512	!! ** input :
	513	!! Namelist namiceini
	514	!!
	515	!! history :
	516	!! 8.5 ! 03-08 (C. Ethe) original code
	517	!! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization
[825]	518	!!-----------------------------------------------------------------------------
[6515]	519	!
	520	INTEGER :: ios,ierr,inum_ice ! Local integer output status for namelist read
	521	INTEGER :: ji,jj
	522	INTEGER :: ifpr, ierror
	523	!
	524	CHARACTER(len=100) :: cn_dir ! Root directory for location of ice files
	525	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
	526	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	527	!
	528	NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, &
	529	& rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
	530	& rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, &
	531	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
[825]	532	!!-----------------------------------------------------------------------------
[2528]	533	!
[4147]	534	REWIND( numnam_ice_ref ) ! Namelist namiceini in reference namelist : Ice initial state
	535	READ ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
	536	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
	537
	538	REWIND( numnam_ice_cfg ) ! Namelist namiceini in configuration namelist : Ice initial state
	539	READ ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
	540	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
[4624]	541	IF(lwm) WRITE ( numoni, namiceini )
[4161]	542
[6515]	543	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
	544	slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu
	545	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi
	546
[4161]	547	! Define the initial parameters
	548	! -------------------------
	549
	550	IF(lwp) THEN
[825]	551	WRITE(numout,*)
	552	WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
	553	WRITE(numout,*) '~~~~~~~~~~~~~~~'
[6515]	554	WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini = ', ln_limini
	555	WRITE(numout,*) ' ice initialization from a netcdf file ln_limini_file = ', ln_limini_file
[5123]	556	WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst
	557	WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n
	558	WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s
	559	WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n
	560	WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s
	561	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n
	562	WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s
	563	WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n
	564	WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s
	565	WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n
	566	WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s
[825]	567	ENDIF
[4161]	568
[6515]	569	IF( ln_limini_file ) THEN ! Ice initialization using input file
	570	!
	571	! set si structure
	572	ALLOCATE( si(jpfldi), STAT=ierror )
	573	IF( ierror > 0 ) THEN
	574	CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' ) ; RETURN
	575	ENDIF
	576
	577	DO ifpr = 1, jpfldi
	578	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
	579	ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
	580	END DO
	581
	582	! fill si with slf_i and control print
	583	CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
	584
	585	CALL fld_read( nit000, 1, si ) ! input fields provided at the current time-step
	586
	587	ENDIF
	588
[825]	589	END SUBROUTINE lim_istate_init
	590
	591	#else
	592	!!----------------------------------------------------------------------
	593	!! Default option : Empty module NO LIM sea-ice model
	594	!!----------------------------------------------------------------------
	595	CONTAINS
	596	SUBROUTINE lim_istate ! Empty routine
	597	END SUBROUTINE lim_istate
	598	#endif
	599
	600	!!======================================================================
	601	END MODULE limistate

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

Download in other formats: