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.

icevar.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/icevar.F90 @ 8587

Last change on this file since 8587 was 8587, checked in by clem, 7 years ago
bug fixes
File size: 32.0 KB

Rev	Line
[8424]	1	MODULE icevar
	2	!!======================================================================
	3	!! * MODULE icevar *
[8534]	4	!! sea-ice: Different sets of ice model variables
[8424]	5	!! how to switch from one to another
	6	!!
	7	!! There are three sets of variables
	8	!! VGLO : global variables of the model
	9	!! - v_i (jpi,jpj,jpl)
	10	!! - v_s (jpi,jpj,jpl)
	11	!! - a_i (jpi,jpj,jpl)
	12	!! - t_s (jpi,jpj,jpl)
	13	!! - e_i (jpi,jpj,nlay_i,jpl)
[8564]	14	!! - sv_i(jpi,jpj,jpl)
	15	!! - oa_i(jpi,jpj,jpl)
[8424]	16	!! VEQV : equivalent variables sometimes used in the model
[8563]	17	!! - h_i(jpi,jpj,jpl)
	18	!! - h_s(jpi,jpj,jpl)
	19	!! - t_i(jpi,jpj,nlay_i,jpl)
[8424]	20	!! ...
	21	!! VAGG : aggregate variables, averaged/summed over all
	22	!! thickness categories
	23	!! - vt_i(jpi,jpj)
	24	!! - vt_s(jpi,jpj)
	25	!! - at_i(jpi,jpj)
	26	!! - et_s(jpi,jpj) !total snow heat content
	27	!! - et_i(jpi,jpj) !total ice thermal content
[8564]	28	!! - sm_i(jpi,jpj) !mean ice salinity
[8424]	29	!! - tm_i (jpi,jpj) !mean ice temperature
	30	!!======================================================================
	31	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
	32	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
[8486]	33	!! 3.5 ! 2012 (M. Vancoppenolle) add ice_var_itd
	34	!! 3.6 ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
[8424]	35	!!----------------------------------------------------------------------
	36	#if defined key_lim3
	37	!!----------------------------------------------------------------------
[8534]	38	!! 'key_lim3' ESIM sea-ice model
[8424]	39	!!----------------------------------------------------------------------
[8486]	40	!! ice_var_agg : integrate variables over layers and categories
	41	!! ice_var_glo2eqv : transform from VGLO to VEQV
	42	!! ice_var_eqv2glo : transform from VEQV to VGLO
	43	!! ice_var_salprof : salinity profile in the ice
	44	!! ice_var_salprof1d : salinity profile in the ice 1D
	45	!! ice_var_zapsmall : remove very small area and volume
	46	!! ice_var_itd : convert 1-cat to multiple cat
[8559]	47	!! ice_var_bv : brine volume
[8486]	48	!!----------------------------------------------------------------------
[8534]	49	USE dom_oce ! ocean space and time domain
[8424]	50	USE phycst ! physical constants (ocean directory)
	51	USE sbc_oce , ONLY : sss_m
[8534]	52	USE ice ! sea-ice: variables
	53	USE ice1D ! sea-ice: thermodynamics variables
[8424]	54	!
	55	USE in_out_manager ! I/O manager
	56	USE lib_mpp ! MPP library
[8534]	57	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
[8424]	58
	59	IMPLICIT NONE
	60	PRIVATE
	61
	62	PUBLIC ice_var_agg
	63	PUBLIC ice_var_glo2eqv
	64	PUBLIC ice_var_eqv2glo
	65	PUBLIC ice_var_salprof
	66	PUBLIC ice_var_salprof1d
	67	PUBLIC ice_var_zapsmall
	68	PUBLIC ice_var_itd
[8559]	69	PUBLIC ice_var_bv
[8424]	70
	71	!!----------------------------------------------------------------------
[8486]	72	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8424]	73	!! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
	74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	75	!!----------------------------------------------------------------------
	76	CONTAINS
	77
	78	SUBROUTINE ice_var_agg( kn )
[8534]	79	!!-------------------------------------------------------------------
[8424]	80	!! * ROUTINE ice_var_agg *
	81	!!
[8486]	82	!! ** Purpose : aggregates ice-thickness-category variables to
	83	!! all-ice variables, i.e. it turns VGLO into VAGG
[8534]	84	!!-------------------------------------------------------------------
[8486]	85	INTEGER, INTENT( in ) :: kn ! =1 state variables only
	86	! ! >1 state variables + others
[8424]	87	!
[8486]	88	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	89	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i
[8534]	90	!!-------------------------------------------------------------------
[8486]	91	!
	92	! ! integrated values
	93	vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 )
	94	vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 )
	95	at_i(:,:) = SUM( a_i(:,:,:) , dim=3 )
[8424]	96	et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
	97	et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
	98
	99	! MV MP 2016
[8486]	100	IF ( ln_pnd ) THEN ! Melt pond
	101	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
	102	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
[8424]	103	ENDIF
	104	! END MP 2016
	105
[8517]	106	ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction
[8424]	107
	108	IF( kn > 1 ) THEN
	109	!
[8486]	110	ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
[8500]	111	WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:)
[8486]	112	ELSEWHERE ; z1_at_i(:,:) = 0._wp
	113	END WHERE
[8500]	114	WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:)
[8486]	115	ELSEWHERE ; z1_vt_i(:,:) = 0._wp
	116	END WHERE
	117	!
	118	! ! mean ice/snow thickness
[8564]	119	hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
	120	hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
[8486]	121	!
	122	! ! mean temperature (K), salinity and age
[8496]	123	tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
	124	tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
	125	om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
[8486]	126	!
[8564]	127	tm_i(:,:) = 0._wp
	128	sm_i(:,:) = 0._wp
[8488]	129	DO jl = 1, jpl
	130	DO jk = 1, nlay_i
[8564]	131	tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
	132	sm_i(:,:) = sm_i(:,:) + r1_nlay_i * sz_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
[8488]	133	END DO
	134	END DO
	135	!
[8587]	136	! ! put rt0 where there is no ice
	137	WHERE( at_i(:,:)<=epsi20 )
	138	tm_su(:,:) = rt0
	139	tm_si(:,:) = rt0
	140	tm_i (:,:) = rt0
	141	END WHERE
	142
[8486]	143	DEALLOCATE( z1_at_i , z1_vt_i )
[8424]	144	ENDIF
	145	!
	146	END SUBROUTINE ice_var_agg
	147
	148
	149	SUBROUTINE ice_var_glo2eqv
[8534]	150	!!-------------------------------------------------------------------
[8424]	151	!! * ROUTINE ice_var_glo2eqv *
	152	!!
[8486]	153	!! ** Purpose : computes equivalent variables as function of
	154	!! global variables, i.e. it turns VGLO into VEQV
[8534]	155	!!-------------------------------------------------------------------
[8424]	156	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[8522]	157	REAL(wp) :: ze_i ! local scalars
[8498]	158	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
[8564]	159	REAL(wp) :: zhmax, z1_zhmax ! - -
[8498]	160	REAL(wp) :: zlay_i, zlay_s ! - -
[8486]	161	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
[8534]	162	!!-------------------------------------------------------------------
[8424]	163
[8486]	164	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
	165	!! ice area and ice volume ?
	166	!! the idea is to be able to define one for all at the begining of this routine
	167	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
	168
[8424]	169	!-------------------------------------------------------
	170	! Ice thickness, snow thickness, ice salinity, ice age
	171	!-------------------------------------------------------
[8486]	172	! !--- inverse of the ice area
	173	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
	174	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
	175	END WHERE
	176	!
[8559]	177	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
	178	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
	179	END WHERE
	180	!
[8563]	181	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness
[8424]	182
[8498]	183	zhmax = hi_max(jpl)
[8486]	184	z1_zhmax = 1._wp / hi_max(jpl)
[8563]	185	WHERE( h_i(:,:,jpl) > zhmax ) !--- bound h_i by hi_max (i.e. 99 m) with associated update of ice area
	186	h_i (:,:,jpl) = zhmax
[8486]	187	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
[8563]	188	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) ! NB: v_i always /=0 as h_i > hi_max
[8486]	189	END WHERE
	190
[8563]	191	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness
[8424]	192
[8486]	193	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age
	194
[8514]	195	IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere)
[8564]	196	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
	197	ELSEWHERE ; s_i(:,:,:) = rn_simin
[8486]	198	END WHERE
[8424]	199	ENDIF
	200
	201	CALL ice_var_salprof ! salinity profile
	202
	203	!-------------------
[8522]	204	! Ice temperature [K] (with a minimum value (rt0 - 100.))
[8424]	205	!-------------------
[8486]	206	zlay_i = REAL( nlay_i , wp ) ! number of layers
[8424]	207	DO jl = 1, jpl
	208	DO jk = 1, nlay_i
	209	DO jj = 1, jpj
	210	DO ji = 1, jpi
[8486]	211	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
	212	!
[8564]	213	ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
	214	ztmelts = - sz_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C]
[8498]	215	! Conversion q(S,T) -> T (second order equation)
	216	zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
	217	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
[8522]	218	t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
[8486]	219	!
	220	ELSE !--- no ice
[8522]	221	t_i(ji,jj,jk,jl) = rt0
[8486]	222	ENDIF
[8424]	223	END DO
	224	END DO
	225	END DO
	226	END DO
	227
	228	!--------------------
[8522]	229	! Snow temperature [K] (with a minimum value (rt0 - 100.))
[8424]	230	!--------------------
[8486]	231	zlay_s = REAL( nlay_s , wp )
	232	DO jk = 1, nlay_s
	233	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
[8522]	234	t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
[8486]	235	ELSEWHERE !--- no ice
[8522]	236	t_s(:,:,jk,:) = rt0
[8486]	237	END WHERE
[8424]	238	END DO
	239
	240	! integrated values
	241	vt_i (:,:) = SUM( v_i, dim=3 )
	242	vt_s (:,:) = SUM( v_s, dim=3 )
	243	at_i (:,:) = SUM( a_i, dim=3 )
	244
[8486]	245	! MV MP 2016
	246	! probably should resum for melt ponds ???
[8424]	247
	248	!
	249	END SUBROUTINE ice_var_glo2eqv
	250
	251
	252	SUBROUTINE ice_var_eqv2glo
[8534]	253	!!-------------------------------------------------------------------
[8424]	254	!! * ROUTINE ice_var_eqv2glo *
	255	!!
[8486]	256	!! ** Purpose : computes global variables as function of
	257	!! equivalent variables, i.e. it turns VEQV into VGLO
[8534]	258	!!-------------------------------------------------------------------
[8424]	259	!
[8564]	260	v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:)
	261	v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:)
	262	sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:)
[8424]	263	!
	264	END SUBROUTINE ice_var_eqv2glo
	265
	266
	267	SUBROUTINE ice_var_salprof
[8534]	268	!!-------------------------------------------------------------------
[8424]	269	!! * ROUTINE ice_var_salprof *
	270	!!
	271	!! ** Purpose : computes salinity profile in function of bulk salinity
	272	!!
	273	!! ** Method : If bulk salinity greater than zsi1,
	274	!! the profile is assumed to be constant (S_inf)
	275	!! If bulk salinity lower than zsi0,
	276	!! the profile is linear with 0 at the surface (S_zero)
	277	!! If it is between zsi0 and zsi1, it is a
	278	!! alpha-weighted linear combination of s_inf and s_zero
	279	!!
	280	!! ** References : Vancoppenolle et al., 2007
[8534]	281	!!-------------------------------------------------------------------
[8424]	282	INTEGER :: ji, jj, jk, jl ! dummy loop index
[8486]	283	REAL(wp) :: zsal, z1_dS
	284	REAL(wp) :: zargtemp , zs0, zs
	285	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
[8424]	286	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	287	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[8534]	288	!!-------------------------------------------------------------------
[8424]	289
[8486]	290	!!gm Question: Remove the option 3 ? How many years since it last use ?
	291
	292	SELECT CASE ( nn_icesal )
	293	!
[8514]	294	! !---------------------------------------!
	295	CASE( 1 ) ! constant salinity in time and space !
	296	! !---------------------------------------!
[8564]	297	sz_i(:,:,:,:) = rn_icesal
	298	s_i(:,:,:) = rn_icesal
[8424]	299	!
[8514]	300	! !---------------------------------------------!
	301	CASE( 2 ) ! time varying salinity with linear profile !
	302	! !---------------------------------------------!
[8486]	303	!
	304	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
	305	!
[8559]	306	DO jl = 1, jpl
	307	DO jk = 1, nlay_i
[8564]	308	sz_i(:,:,jk,jl) = s_i(:,:,jl)
[8559]	309	END DO
[8424]	310	END DO
[8486]	311	! ! Slope of the linear profile
[8564]	312	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
[8486]	313	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
	314	END WHERE
[8424]	315	!
[8486]	316	z1_dS = 1._wp / ( zsi1 - zsi0 )
[8424]	317	DO jl = 1, jpl
	318	DO jj = 1, jpj
	319	DO ji = 1, jpi
[8564]	320	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
[8486]	321	! ! force a constant profile when SSS too low (Baltic Sea)
[8564]	322	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
[8424]	323	END DO
	324	END DO
	325	END DO
	326
	327	! Computation of the profile
	328	DO jl = 1, jpl
	329	DO jk = 1, nlay_i
	330	DO jj = 1, jpj
	331	DO ji = 1, jpi
[8486]	332	! ! linear profile with 0 surface value
[8563]	333	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
[8564]	334	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
	335	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
[8424]	336	END DO
	337	END DO
	338	END DO
	339	END DO
	340	!
[8486]	341	DEALLOCATE( z_slope_s , zalpha )
[8424]	342	!
[8514]	343	! !-------------------------------------------!
	344	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	345	! !-------------------------------------------! (mean = 2.30)
[8486]	346	!
[8564]	347	s_i(:,:,:) = 2.30_wp
[8486]	348	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
	349	!! a array S_prof(1:nlay_i) containing the calculation and just do:
	350	! DO jk = 1, nlay_i
[8564]	351	! sz_i(:,:,jk,:) = S_prof(jk)
[8486]	352	! END DO
	353	!!gm end
[8424]	354	!
	355	DO jl = 1, jpl
	356	DO jk = 1, nlay_i
	357	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
[8564]	358	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
[8424]	359	END DO
	360	END DO
	361	!
[8486]	362	END SELECT
[8424]	363	!
	364	END SUBROUTINE ice_var_salprof
	365
	366	SUBROUTINE ice_var_salprof1d
	367	!!-------------------------------------------------------------------
	368	!! * ROUTINE ice_var_salprof1d *
	369	!!
	370	!! ** Purpose : 1d computation of the sea ice salinity profile
	371	!! Works with 1d vectors and is used by thermodynamic modules
	372	!!-------------------------------------------------------------------
	373	INTEGER :: ji, jk ! dummy loop indices
[8486]	374	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
[8559]	375	REAL(wp) :: zs, zs0 ! - -
[8424]	376	!
[8559]	377	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
[8424]	378	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	379	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[8534]	380	!!-------------------------------------------------------------------
[8486]	381	!
	382	SELECT CASE ( nn_icesal )
	383	!
[8514]	384	! !---------------------------------------!
	385	CASE( 1 ) ! constant salinity in time and space !
	386	! !---------------------------------------!
[8565]	387	sz_i_1d(1:npti,:) = rn_icesal
[8424]	388	!
[8514]	389	! !---------------------------------------------!
	390	CASE( 2 ) ! time varying salinity with linear profile !
	391	! !---------------------------------------------!
[8486]	392	!
[8559]	393	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
[8486]	394	!
[8559]	395	! ! Slope of the linear profile
[8565]	396	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
	397	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
[8559]	398	END WHERE
	399
	400	z1_dS = 1._wp / ( zsi1 - zsi0 )
[8565]	401	DO ji = 1, npti
[8564]	402	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
[8559]	403	! ! force a constant profile when SSS too low (Baltic Sea)
[8564]	404	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
[8424]	405	END DO
[8559]	406	!
	407	! Computation of the profile
[8424]	408	DO jk = 1, nlay_i
[8565]	409	DO ji = 1, npti
[8559]	410	! ! linear profile with 0 surface value
[8563]	411	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
[8564]	412	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
	413	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
[8486]	414	END DO
	415	END DO
	416	!
[8559]	417	DEALLOCATE( z_slope_s, zalpha )
[8424]	418
[8514]	419	! !-------------------------------------------!
	420	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	421	! !-------------------------------------------! (mean = 2.30)
[8424]	422	!
[8565]	423	s_i_1d(1:npti) = 2.30_wp
[8424]	424	!
[8486]	425	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
[8424]	426	DO jk = 1, nlay_i
	427	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	428	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
[8565]	429	DO ji = 1, npti
[8564]	430	sz_i_1d(ji,jk) = zsal
[8424]	431	END DO
	432	END DO
	433	!
[8486]	434	END SELECT
[8424]	435	!
	436	END SUBROUTINE ice_var_salprof1d
	437
[8486]	438
[8424]	439	SUBROUTINE ice_var_zapsmall
	440	!!-------------------------------------------------------------------
	441	!! * ROUTINE ice_var_zapsmall *
	442	!!
	443	!! ** Purpose : Remove too small sea ice areas and correct fluxes
	444	!!-------------------------------------------------------------------
	445	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[8559]	446	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
[8424]	447	!!-------------------------------------------------------------------
[8486]	448	!
	449	DO jl = 1, jpl !== loop over the categories ==!
	450	!
[8424]	451	!-----------------------------------------------------------------
	452	! Zap ice energy and use ocean heat to melt ice
	453	!-----------------------------------------------------------------
[8563]	454	WHERE( a_i(:,:,jl) > epsi20 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
	455	ELSEWHERE ; h_i(:,:,jl) = 0._wp
[8559]	456	END WHERE
	457	!
[8563]	458	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
[8564]	459	ELSEWHERE ; zswitch(:,:) = 1._wp
[8559]	460	END WHERE
	461
[8424]	462	DO jk = 1, nlay_i
	463	DO jj = 1 , jpj
	464	DO ji = 1 , jpi
	465	! update exchanges with ocean
[8559]	466	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
	467	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
	468	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
[8424]	469	END DO
	470	END DO
	471	END DO
	472
	473	DO jj = 1 , jpj
	474	DO ji = 1 , jpi
[8559]	475	! update exchanges with ocean
[8564]	476	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_rdtice
	477	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
	478	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
	479	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
[8424]	480	!-----------------------------------------------------------------
	481	! Zap snow energy
	482	!-----------------------------------------------------------------
[8559]	483	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
	484	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
[8424]	485
	486	!-----------------------------------------------------------------
	487	! zap ice and snow volume, add water and salt to ocean
	488	!-----------------------------------------------------------------
[8564]	489	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
	490	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
	491	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
	492	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
	493	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
	494	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
	495	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
[8424]	496
[8563]	497	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
	498	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
[8424]	499
[8559]	500	END DO
	501	END DO
[8424]	502
[8559]	503	IF( ln_pnd ) THEN
	504	DO jj = 1 , jpj
	505	DO ji = 1 , jpi
	506	IF( ln_pnd_fw ) &
	507	& wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
	508	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
	509	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
	510	END DO
[8424]	511	END DO
[8559]	512	ENDIF
	513
[8424]	514	END DO
	515
	516	! to be sure that at_i is the sum of a_i(jl)
[8517]	517	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
	518	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
[8424]	519
[8517]	520	! open water = 1 if at_i=0
	521	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
[8424]	522	!
	523	END SUBROUTINE ice_var_zapsmall
	524
[8486]	525
[8563]	526	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
[8534]	527	!!-------------------------------------------------------------------
[8424]	528	!! * ROUTINE ice_var_itd *
	529	!!
	530	!! ** Purpose : converting 1-cat ice to multiple ice categories
	531	!!
	532	!! ice thickness distribution follows a gaussian law
	533	!! around the concentration of the most likely ice thickness
[8514]	534	!! (similar as iceistate.F90)
[8424]	535	!!
	536	!! ** Method: Iterative procedure
	537	!!
	538	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	539	!!
	540	!! 2) Check whether the distribution conserves area and volume, positivity and
	541	!! category boundaries
	542	!!
	543	!! 3) If not (input ice is too thin), the last category is empty and
	544	!! the number of categories is reduced (jpl-1)
	545	!!
	546	!! 4) Iterate until ok (SUM(itest(:) = 4)
	547	!!
	548	!! ** Arguments : zhti: 1-cat ice thickness
	549	!! zhts: 1-cat snow depth
[8563]	550	!! zati: 1-cat ice concentration
[8424]	551	!!
	552	!! ** Output : jpl-cat
	553	!!
	554	!! (Example of application: BDY forcings when input are cell averaged)
	555	!!-------------------------------------------------------------------
	556	INTEGER :: ji, jk, jl ! dummy loop indices
	557	INTEGER :: ijpij, i_fill, jl0
	558	REAL(wp) :: zarg, zV, zconv, zdh, zdv
[8563]	559	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
	560	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
[8424]	561	INTEGER , DIMENSION(4) :: itest
[8486]	562	!!-------------------------------------------------------------------
[8534]	563	!
	564	! ----------------------------------------
	565	! distribution over the jpl ice categories
	566	! ----------------------------------------
	567	! a gaussian distribution for ice concentration is used
	568	! then we check whether the distribution fullfills
	569	! volume and area conservation, positivity and ice categories bounds
[8486]	570	ijpij = SIZE( zhti , 1 )
[8563]	571	zh_i(1:ijpij,1:jpl) = 0._wp
	572	zh_s(1:ijpij,1:jpl) = 0._wp
[8424]	573	za_i (1:ijpij,1:jpl) = 0._wp
	574
	575	DO ji = 1, ijpij
	576
	577	IF( zhti(ji) > 0._wp ) THEN
	578
	579	! find which category (jl0) the input ice thickness falls into
	580	jl0 = jpl
	581	DO jl = 1, jpl
	582	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
	583	jl0 = jl
	584	CYCLE
	585	ENDIF
	586	END DO
	587
[8534]	588	itest(:) = 0
	589	i_fill = jpl + 1 !------------------------------------
	590	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	591	! !------------------------------------
[8424]	592	i_fill = i_fill - 1
[8534]	593	!
[8563]	594	zh_i(ji,1:jpl) = 0._wp
[8424]	595	za_i (ji,1:jpl) = 0._wp
	596	itest(:) = 0
	597
[8534]	598	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
[8563]	599	zh_i(ji,1) = zhti(ji)
	600	za_i (ji,1) = zati (ji)
[8534]	601	ELSE !-- case ice is thicker: fill categories >1
	602	! thickness
[8424]	603	DO jl = 1, i_fill - 1
[8563]	604	zh_i(ji,jl) = hi_mean(jl)
[8424]	605	END DO
	606
[8534]	607	! concentration
[8563]	608	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
[8424]	609	DO jl = 1, i_fill - 1
	610	IF ( jl /= jl0 ) THEN
[8563]	611	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
[8424]	612	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
	613	ENDIF
	614	END DO
	615
[8534]	616	! last category
[8563]	617	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
	618	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
	619	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
[8424]	620
	621	! clem: correction if concentration of upper cat is greater than lower cat
	622	! (it should be a gaussian around jl0 but sometimes it is not)
	623	IF ( jl0 /= jpl ) THEN
	624	DO jl = jpl, jl0+1, -1
	625	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
[8563]	626	zdv = zh_i(ji,jl) * za_i(ji,jl)
	627	zh_i(ji,jl ) = 0._wp
[8424]	628	za_i (ji,jl ) = 0._wp
	629	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
	630	END IF
	631	ENDDO
	632	ENDIF
	633
[8534]	634	ENDIF
[8424]	635
	636	! Compatibility tests
[8563]	637	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
[8534]	638	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
[8424]	639
[8563]	640	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
[8534]	641	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
[8424]	642
[8563]	643	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
[8424]	644
	645	itest(4) = 1
	646	DO jl = 1, i_fill
[8534]	647	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
[8424]	648	END DO
[8534]	649	! !----------------------------
[8424]	650	END DO ! end iteration on categories
[8534]	651	! !----------------------------
	652	ENDIF
	653	END DO
[8424]	654
[8534]	655	! Add Snow in each category where za_i is not 0
[8424]	656	DO jl = 1, jpl
	657	DO ji = 1, ijpij
	658	IF( za_i(ji,jl) > 0._wp ) THEN
[8563]	659	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
[8424]	660	! In case snow load is in excess that would lead to transformation from snow to ice
	661	! Then, transfer the snow excess into the ice (different from icethd_dh)
[8563]	662	zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
	663	! recompute h_i, h_s avoiding out of bounds values
	664	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
	665	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
[8424]	666	ENDIF
[8486]	667	END DO
	668	END DO
[8424]	669	!
	670	END SUBROUTINE ice_var_itd
	671
[8559]	672
	673	SUBROUTINE ice_var_bv
	674	!!-------------------------------------------------------------------
	675	!! * ROUTINE ice_var_bv *
	676	!!
	677	!! ** Purpose : computes mean brine volume (%) in sea ice
	678	!!
	679	!! ** Method : e = - 0.054 * S (ppt) / T (C)
	680	!!
	681	!! References : Vancoppenolle et al., JGR, 2007
	682	!!-------------------------------------------------------------------
	683	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	684	!!-------------------------------------------------------------------
	685	!
	686	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
	687	!! instead of setting everything to zero as just below
	688	bv_i (:,:,:) = 0._wp
	689	DO jl = 1, jpl
	690	DO jk = 1, nlay_i
	691	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
[8564]	692	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
[8559]	693	END WHERE
	694	END DO
	695	END DO
	696	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
	697	ELSEWHERE ; bvm_i(:,:) = 0._wp
	698	END WHERE
	699	!
	700	END SUBROUTINE ice_var_bv
	701
	702
[8424]	703	#else
	704	!!----------------------------------------------------------------------
[8534]	705	!! Default option Dummy module NO ESIM sea-ice model
[8424]	706	!!----------------------------------------------------------------------
	707	#endif
	708
	709	!!======================================================================
	710	END MODULE icevar

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

Download in other formats: