Context Navigation

icevar.F90

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

Download in other formats: