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

Last change on this file since 8517 was 8517, checked in by clem, 7 years ago
changes in style - part6 - one more round
File size: 32.7 KB

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

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

Download in other formats: