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 @ 8559

Last change on this file since 8559 was 8559, checked in by clem, 7 years ago
almost useless commits
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)
	14	!! - smv_i(jpi,jpj,jpl)
	15	!! - oa_i (jpi,jpj,jpl)
	16	!! VEQV : equivalent variables sometimes used in the model
	17	!! - ht_i(jpi,jpj,jpl)
	18	!! - ht_s(jpi,jpj,jpl)
	19	!! - t_i (jpi,jpj,nlay_i,jpl)
	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
	28	!! - smt_i(jpi,jpj) !mean ice salinity
	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
	119	htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
	120	htm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
	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	!
[8488]	127	tm_i (:,:) = 0._wp
	128	smt_i(:,:) = 0._wp
	129	DO jl = 1, jpl
	130	DO jk = 1, nlay_i
	131	tm_i (:,:) = tm_i (:,:) + r1_nlay_i * t_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
	132	smt_i(:,:) = smt_i(:,:) + r1_nlay_i * s_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
	133	END DO
	134	END DO
	135	!
[8486]	136	!!gm QUESTION 1 : why salinity is named smt_i and not just sm_i ? since the 4D field is named s_i. (NB for temp: tm_i, t_i)
	137	!
	138	DEALLOCATE( z1_at_i , z1_vt_i )
[8424]	139	ENDIF
	140	!
	141	END SUBROUTINE ice_var_agg
	142
	143
	144	SUBROUTINE ice_var_glo2eqv
[8534]	145	!!-------------------------------------------------------------------
[8424]	146	!! * ROUTINE ice_var_glo2eqv *
	147	!!
[8486]	148	!! ** Purpose : computes equivalent variables as function of
	149	!! global variables, i.e. it turns VGLO into VEQV
[8534]	150	!!-------------------------------------------------------------------
[8424]	151	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[8522]	152	REAL(wp) :: ze_i ! local scalars
[8498]	153	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
	154	REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - -
	155	REAL(wp) :: zlay_i, zlay_s ! - -
[8486]	156	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
[8534]	157	!!-------------------------------------------------------------------
[8424]	158
[8486]	159	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
	160	!! ice area and ice volume ?
	161	!! the idea is to be able to define one for all at the begining of this routine
	162	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
	163
[8424]	164	!-------------------------------------------------------
	165	! Ice thickness, snow thickness, ice salinity, ice age
	166	!-------------------------------------------------------
[8486]	167	! !--- inverse of the ice area
	168	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
	169	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
	170	END WHERE
	171	!
[8559]	172	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
	173	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
	174	END WHERE
	175	!
[8486]	176	ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness
[8424]	177
[8498]	178	zhmax = hi_max(jpl)
[8486]	179	z1_zhmax = 1._wp / hi_max(jpl)
	180	WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
	181	ht_i (:,:,jpl) = zhmax
	182	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
[8559]	183	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max
[8486]	184	END WHERE
	185
	186	ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness
[8424]	187
[8486]	188	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age
	189
[8514]	190	IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere)
[8559]	191	WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) )
[8486]	192	ELSEWHERE ; sm_i(:,:,:) = rn_simin
	193	END WHERE
[8424]	194	ENDIF
	195
	196	CALL ice_var_salprof ! salinity profile
	197
	198	!-------------------
[8522]	199	! Ice temperature [K] (with a minimum value (rt0 - 100.))
[8424]	200	!-------------------
[8486]	201	zlay_i = REAL( nlay_i , wp ) ! number of layers
[8424]	202	DO jl = 1, jpl
	203	DO jk = 1, nlay_i
	204	DO jj = 1, jpj
	205	DO ji = 1, jpi
[8486]	206	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
	207	!
[8559]	208	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]
[8498]	209	ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C]
	210	! Conversion q(S,T) -> T (second order equation)
	211	zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
	212	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
[8522]	213	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]	214	!
	215	ELSE !--- no ice
[8522]	216	t_i(ji,jj,jk,jl) = rt0
[8486]	217	ENDIF
[8424]	218	END DO
	219	END DO
	220	END DO
	221	END DO
	222
	223	!--------------------
[8522]	224	! Snow temperature [K] (with a minimum value (rt0 - 100.))
[8424]	225	!--------------------
[8486]	226	zlay_s = REAL( nlay_s , wp )
	227	DO jk = 1, nlay_s
	228	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
[8522]	229	t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
[8486]	230	ELSEWHERE !--- no ice
[8522]	231	t_s(:,:,jk,:) = rt0
[8486]	232	END WHERE
[8424]	233	END DO
	234
	235	! integrated values
	236	vt_i (:,:) = SUM( v_i, dim=3 )
	237	vt_s (:,:) = SUM( v_s, dim=3 )
	238	at_i (:,:) = SUM( a_i, dim=3 )
	239
[8486]	240	! MV MP 2016
	241	! probably should resum for melt ponds ???
[8424]	242
	243	!
	244	END SUBROUTINE ice_var_glo2eqv
	245
	246
	247	SUBROUTINE ice_var_eqv2glo
[8534]	248	!!-------------------------------------------------------------------
[8424]	249	!! * ROUTINE ice_var_eqv2glo *
	250	!!
[8486]	251	!! ** Purpose : computes global variables as function of
	252	!! equivalent variables, i.e. it turns VEQV into VGLO
[8534]	253	!!-------------------------------------------------------------------
[8424]	254	!
[8486]	255	v_i (:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
	256	v_s (:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
[8424]	257	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
	258	!
	259	END SUBROUTINE ice_var_eqv2glo
	260
	261
	262	SUBROUTINE ice_var_salprof
[8534]	263	!!-------------------------------------------------------------------
[8424]	264	!! * ROUTINE ice_var_salprof *
	265	!!
	266	!! ** Purpose : computes salinity profile in function of bulk salinity
	267	!!
	268	!! ** Method : If bulk salinity greater than zsi1,
	269	!! the profile is assumed to be constant (S_inf)
	270	!! If bulk salinity lower than zsi0,
	271	!! the profile is linear with 0 at the surface (S_zero)
	272	!! If it is between zsi0 and zsi1, it is a
	273	!! alpha-weighted linear combination of s_inf and s_zero
	274	!!
	275	!! ** References : Vancoppenolle et al., 2007
[8534]	276	!!-------------------------------------------------------------------
[8424]	277	INTEGER :: ji, jj, jk, jl ! dummy loop index
[8486]	278	REAL(wp) :: zsal, z1_dS
	279	REAL(wp) :: zargtemp , zs0, zs
	280	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
[8424]	281	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	282	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[8534]	283	!!-------------------------------------------------------------------
[8424]	284
[8486]	285	!!gm Question: Remove the option 3 ? How many years since it last use ?
	286
	287	SELECT CASE ( nn_icesal )
	288	!
[8514]	289	! !---------------------------------------!
	290	CASE( 1 ) ! constant salinity in time and space !
	291	! !---------------------------------------!
[8424]	292	s_i (:,:,:,:) = rn_icesal
	293	sm_i(:,:,:) = rn_icesal
	294	!
[8514]	295	! !---------------------------------------------!
	296	CASE( 2 ) ! time varying salinity with linear profile !
	297	! !---------------------------------------------!
[8486]	298	!
	299	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
	300	!
[8559]	301	DO jl = 1, jpl
	302	DO jk = 1, nlay_i
	303	s_i(:,:,jk,jl) = sm_i(:,:,jl)
	304	END DO
[8424]	305	END DO
[8486]	306	! ! Slope of the linear profile
	307	WHERE( ht_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:)
	308	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
	309	END WHERE
[8424]	310	!
[8486]	311	z1_dS = 1._wp / ( zsi1 - zsi0 )
[8424]	312	DO jl = 1, jpl
	313	DO jj = 1, jpj
	314	DO ji = 1, jpi
[8486]	315	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
	316	! ! force a constant profile when SSS too low (Baltic Sea)
	317	IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
[8424]	318	END DO
	319	END DO
	320	END DO
	321
	322	! Computation of the profile
	323	DO jl = 1, jpl
	324	DO jk = 1, nlay_i
	325	DO jj = 1, jpj
	326	DO ji = 1, jpi
[8486]	327	! ! linear profile with 0 surface value
	328	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
	329	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) ! weighting the profile
	330	s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
[8424]	331	END DO
	332	END DO
	333	END DO
	334	END DO
	335	!
[8486]	336	DEALLOCATE( z_slope_s , zalpha )
[8424]	337	!
[8514]	338	! !-------------------------------------------!
	339	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	340	! !-------------------------------------------! (mean = 2.30)
[8486]	341	!
[8424]	342	sm_i(:,:,:) = 2.30_wp
[8486]	343	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
	344	!! a array S_prof(1:nlay_i) containing the calculation and just do:
	345	! DO jk = 1, nlay_i
	346	! s_i(:,:,jk,:) = S_prof(jk)
	347	! END DO
	348	!!gm end
[8424]	349	!
	350	DO jl = 1, jpl
	351	DO jk = 1, nlay_i
	352	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
[8486]	353	s_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
[8424]	354	END DO
	355	END DO
	356	!
[8486]	357	END SELECT
[8424]	358	!
	359	END SUBROUTINE ice_var_salprof
	360
	361	SUBROUTINE ice_var_salprof1d
	362	!!-------------------------------------------------------------------
	363	!! * ROUTINE ice_var_salprof1d *
	364	!!
	365	!! ** Purpose : 1d computation of the sea ice salinity profile
	366	!! Works with 1d vectors and is used by thermodynamic modules
	367	!!-------------------------------------------------------------------
	368	INTEGER :: ji, jk ! dummy loop indices
[8486]	369	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
[8559]	370	REAL(wp) :: zs, zs0 ! - -
[8424]	371	!
[8559]	372	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
[8424]	373	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	374	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[8534]	375	!!-------------------------------------------------------------------
[8486]	376	!
	377	SELECT CASE ( nn_icesal )
	378	!
[8514]	379	! !---------------------------------------!
	380	CASE( 1 ) ! constant salinity in time and space !
	381	! !---------------------------------------!
[8559]	382	s_i_1d(1:nidx,:) = rn_icesal
[8424]	383	!
[8514]	384	! !---------------------------------------------!
	385	CASE( 2 ) ! time varying salinity with linear profile !
	386	! !---------------------------------------------!
[8486]	387	!
[8559]	388	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
[8486]	389	!
[8559]	390	! ! Slope of the linear profile
	391	WHERE( ht_i_1d(1:nidx) > epsi20 ) ; z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx)
	392	ELSEWHERE ; z_slope_s(1:nidx) = 0._wp
	393	END WHERE
	394
	395	z1_dS = 1._wp / ( zsi1 - zsi0 )
	396	DO ji = 1, nidx
	397	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) )
	398	! ! force a constant profile when SSS too low (Baltic Sea)
	399	IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
[8424]	400	END DO
[8559]	401	!
	402	! Computation of the profile
[8424]	403	DO jk = 1, nlay_i
	404	DO ji = 1, nidx
[8559]	405	! ! linear profile with 0 surface value
	406	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
	407	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji)
	408	s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
[8486]	409	END DO
	410	END DO
	411	!
[8559]	412	DEALLOCATE( z_slope_s, zalpha )
[8424]	413
[8514]	414	! !-------------------------------------------!
	415	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	416	! !-------------------------------------------! (mean = 2.30)
[8424]	417	!
[8559]	418	sm_i_1d(1:nidx) = 2.30_wp
[8424]	419	!
[8486]	420	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
[8424]	421	DO jk = 1, nlay_i
	422	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	423	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
	424	DO ji = 1, nidx
	425	s_i_1d(ji,jk) = zsal
	426	END DO
	427	END DO
	428	!
[8486]	429	END SELECT
[8424]	430	!
	431	END SUBROUTINE ice_var_salprof1d
	432
[8486]	433
[8424]	434	SUBROUTINE ice_var_zapsmall
	435	!!-------------------------------------------------------------------
	436	!! * ROUTINE ice_var_zapsmall *
	437	!!
	438	!! ** Purpose : Remove too small sea ice areas and correct fluxes
	439	!!-------------------------------------------------------------------
	440	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[8559]	441	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
[8424]	442	!!-------------------------------------------------------------------
[8486]	443	!
	444	DO jl = 1, jpl !== loop over the categories ==!
	445	!
[8424]	446	!-----------------------------------------------------------------
	447	! Zap ice energy and use ocean heat to melt ice
	448	!-----------------------------------------------------------------
[8559]	449	WHERE( a_i(:,:,jl) > epsi20 ) ; ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
	450	ELSEWHERE ; ht_i(:,:,jl) = 0._wp
	451	END WHERE
	452	!
	453	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
	454	ELSEWHERE ; zswitch(:,:) = 1._wp
	455	END WHERE
	456
[8424]	457	DO jk = 1, nlay_i
	458	DO jj = 1 , jpj
	459	DO ji = 1 , jpi
	460	! update exchanges with ocean
[8559]	461	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
	462	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
	463	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
[8424]	464	END DO
	465	END DO
	466	END DO
	467
	468	DO jj = 1 , jpj
	469	DO ji = 1 , jpi
[8559]	470	! update exchanges with ocean
	471	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl) * rhoic * r1_rdtice
	472	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
	473	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
	474	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]	475	!-----------------------------------------------------------------
	476	! Zap snow energy
	477	!-----------------------------------------------------------------
[8559]	478	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
	479	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
[8424]	480
	481	!-----------------------------------------------------------------
	482	! zap ice and snow volume, add water and salt to ocean
	483	!-----------------------------------------------------------------
[8559]	484	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
	485	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
	486	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
	487	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
	488	t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
	489	oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj)
	490	smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj)
[8424]	491
[8559]	492	ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj)
	493	ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj)
[8424]	494
[8559]	495	END DO
	496	END DO
[8424]	497
[8559]	498	IF( ln_pnd ) THEN
	499	DO jj = 1 , jpj
	500	DO ji = 1 , jpi
	501	IF( ln_pnd_fw ) &
	502	& wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
	503	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
	504	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
	505	END DO
[8424]	506	END DO
[8559]	507	ENDIF
	508
[8424]	509	END DO
	510
	511	! to be sure that at_i is the sum of a_i(jl)
[8517]	512	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
	513	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
[8424]	514
[8517]	515	! open water = 1 if at_i=0
	516	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
[8424]	517	!
	518	END SUBROUTINE ice_var_zapsmall
	519
[8486]	520
[8424]	521	SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
[8534]	522	!!-------------------------------------------------------------------
[8424]	523	!! * ROUTINE ice_var_itd *
	524	!!
	525	!! ** Purpose : converting 1-cat ice to multiple ice categories
	526	!!
	527	!! ice thickness distribution follows a gaussian law
	528	!! around the concentration of the most likely ice thickness
[8514]	529	!! (similar as iceistate.F90)
[8424]	530	!!
	531	!! ** Method: Iterative procedure
	532	!!
	533	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	534	!!
	535	!! 2) Check whether the distribution conserves area and volume, positivity and
	536	!! category boundaries
	537	!!
	538	!! 3) If not (input ice is too thin), the last category is empty and
	539	!! the number of categories is reduced (jpl-1)
	540	!!
	541	!! 4) Iterate until ok (SUM(itest(:) = 4)
	542	!!
	543	!! ** Arguments : zhti: 1-cat ice thickness
	544	!! zhts: 1-cat snow depth
	545	!! zai : 1-cat ice concentration
	546	!!
	547	!! ** Output : jpl-cat
	548	!!
	549	!! (Example of application: BDY forcings when input are cell averaged)
	550	!!-------------------------------------------------------------------
	551	INTEGER :: ji, jk, jl ! dummy loop indices
	552	INTEGER :: ijpij, i_fill, jl0
	553	REAL(wp) :: zarg, zV, zconv, zdh, zdv
	554	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables
	555	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables
	556	INTEGER , DIMENSION(4) :: itest
[8486]	557	!!-------------------------------------------------------------------
[8534]	558	!
	559	! ----------------------------------------
	560	! distribution over the jpl ice categories
	561	! ----------------------------------------
	562	! a gaussian distribution for ice concentration is used
	563	! then we check whether the distribution fullfills
	564	! volume and area conservation, positivity and ice categories bounds
[8486]	565	ijpij = SIZE( zhti , 1 )
[8424]	566	zht_i(1:ijpij,1:jpl) = 0._wp
	567	zht_s(1:ijpij,1:jpl) = 0._wp
	568	za_i (1:ijpij,1:jpl) = 0._wp
	569
	570	DO ji = 1, ijpij
	571
	572	IF( zhti(ji) > 0._wp ) THEN
	573
	574	! find which category (jl0) the input ice thickness falls into
	575	jl0 = jpl
	576	DO jl = 1, jpl
	577	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
	578	jl0 = jl
	579	CYCLE
	580	ENDIF
	581	END DO
	582
[8534]	583	itest(:) = 0
	584	i_fill = jpl + 1 !------------------------------------
	585	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	586	! !------------------------------------
[8424]	587	i_fill = i_fill - 1
[8534]	588	!
[8424]	589	zht_i(ji,1:jpl) = 0._wp
	590	za_i (ji,1:jpl) = 0._wp
	591	itest(:) = 0
	592
[8534]	593	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
[8424]	594	zht_i(ji,1) = zhti(ji)
	595	za_i (ji,1) = zai (ji)
[8534]	596	ELSE !-- case ice is thicker: fill categories >1
	597	! thickness
[8424]	598	DO jl = 1, i_fill - 1
	599	zht_i(ji,jl) = hi_mean(jl)
	600	END DO
	601
[8534]	602	! concentration
[8424]	603	za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
	604	DO jl = 1, i_fill - 1
	605	IF ( jl /= jl0 ) THEN
	606	zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
	607	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
	608	ENDIF
	609	END DO
	610
[8534]	611	! last category
[8424]	612	za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
	613	zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
	614	zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
	615
	616	! clem: correction if concentration of upper cat is greater than lower cat
	617	! (it should be a gaussian around jl0 but sometimes it is not)
	618	IF ( jl0 /= jpl ) THEN
	619	DO jl = jpl, jl0+1, -1
	620	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
	621	zdv = zht_i(ji,jl) * za_i(ji,jl)
	622	zht_i(ji,jl ) = 0._wp
	623	za_i (ji,jl ) = 0._wp
	624	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
	625	END IF
	626	ENDDO
	627	ENDIF
	628
[8534]	629	ENDIF
[8424]	630
	631	! Compatibility tests
[8534]	632	zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
	633	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
[8424]	634
	635	zconv = ABS( zhti(ji)zai(ji) - SUM( za_i(ji,1:jpl)zht_i(ji,1:jpl) ) )
[8534]	636	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
[8424]	637
[8534]	638	IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
[8424]	639
	640	itest(4) = 1
	641	DO jl = 1, i_fill
[8534]	642	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
[8424]	643	END DO
[8534]	644	! !----------------------------
[8424]	645	END DO ! end iteration on categories
[8534]	646	! !----------------------------
	647	ENDIF
	648	END DO
[8424]	649
[8534]	650	! Add Snow in each category where za_i is not 0
[8424]	651	DO jl = 1, jpl
	652	DO ji = 1, ijpij
	653	IF( za_i(ji,jl) > 0._wp ) THEN
	654	zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
	655	! In case snow load is in excess that would lead to transformation from snow to ice
	656	! Then, transfer the snow excess into the ice (different from icethd_dh)
	657	zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )
	658	! recompute ht_i, ht_s avoiding out of bounds values
	659	zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
	660	zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
	661	ENDIF
[8486]	662	END DO
	663	END DO
[8424]	664	!
	665	END SUBROUTINE ice_var_itd
	666
[8559]	667
	668	SUBROUTINE ice_var_bv
	669	!!-------------------------------------------------------------------
	670	!! * ROUTINE ice_var_bv *
	671	!!
	672	!! ** Purpose : computes mean brine volume (%) in sea ice
	673	!!
	674	!! ** Method : e = - 0.054 * S (ppt) / T (C)
	675	!!
	676	!! References : Vancoppenolle et al., JGR, 2007
	677	!!-------------------------------------------------------------------
	678	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	679	!!-------------------------------------------------------------------
	680	!
	681	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
	682	!! instead of setting everything to zero as just below
	683	bv_i (:,:,:) = 0._wp
	684	DO jl = 1, jpl
	685	DO jk = 1, nlay_i
	686	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
	687	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
	688	END WHERE
	689	END DO
	690	END DO
	691	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
	692	ELSEWHERE ; bvm_i(:,:) = 0._wp
	693	END WHERE
	694	!
	695	END SUBROUTINE ice_var_bv
	696
	697
[8424]	698	#else
	699	!!----------------------------------------------------------------------
[8534]	700	!! Default option Dummy module NO ESIM sea-ice model
[8424]	701	!!----------------------------------------------------------------------
	702	#endif
	703
	704	!!======================================================================
	705	END MODULE icevar

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

Download in other formats: