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

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

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

Download in other formats: