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.

limvar.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/limvar.F90 @ 8413

Last change on this file since 8413 was 8413, checked in by clem, 7 years ago
continue changing names (again)
Property svn:keywords set to `Id`
File size: 34.4 KB

Rev	Line
[825]	1	MODULE limvar
	2	!!======================================================================
	3	!! * MODULE limvar *
	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
[7646]	29	!! - tm_i (jpi,jpj) !mean ice temperature
[825]	30	!!======================================================================
[2715]	31	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
[5123]	32	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
[2715]	33	!!----------------------------------------------------------------------
[888]	34	#if defined key_lim3
[825]	35	!!----------------------------------------------------------------------
[2715]	36	!! 'key_lim3' LIM3 sea-ice model
	37	!!----------------------------------------------------------------------
[3625]	38	USE par_oce ! ocean parameters
	39	USE phycst ! physical constants (ocean directory)
[8316]	40	USE sbc_oce , ONLY : sss_m
[3625]	41	USE ice ! ice variables
	42	USE thd_ice ! ice variables (thermodynamics)
[8378]	43	!
[3625]	44	USE in_out_manager ! I/O manager
	45	USE lib_mpp ! MPP library
	46	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[921]	47
[825]	48	IMPLICIT NONE
	49	PRIVATE
	50
[5123]	51	PUBLIC lim_var_agg
	52	PUBLIC lim_var_glo2eqv
	53	PUBLIC lim_var_eqv2glo
	54	PUBLIC lim_var_salprof
	55	PUBLIC lim_var_bv
	56	PUBLIC lim_var_salprof1d
	57	PUBLIC lim_var_zapsmall
	58	PUBLIC lim_var_itd
[825]	59
	60	!!----------------------------------------------------------------------
[5123]	61	!! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
[1156]	62	!! $Id$
[2715]	63	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	64	!!----------------------------------------------------------------------
	65	CONTAINS
	66
[2715]	67	SUBROUTINE lim_var_agg( kn )
[921]	68	!!------------------------------------------------------------------
	69	!! * ROUTINE lim_var_agg *
[2715]	70	!!
	71	!! ** Purpose : aggregates ice-thickness-category variables to all-ice variables
	72	!! i.e. it turns VGLO into VAGG
[921]	73	!! ** Method :
	74	!!
	75	!! ** Arguments : n = 1, at_i vt_i only
	76	!! n = 2 everything
	77	!!
	78	!! note : you could add an argument when you need only at_i, vt_i
	79	!! and when you need everything
	80	!!------------------------------------------------------------------
[2715]	81	INTEGER, INTENT( in ) :: kn ! =1 at_i & vt only ; = what is needed
	82	!
	83	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	84	!!------------------------------------------------------------------
[825]	85
[7646]	86	! integrated values
[7753]	87	vt_i (:,:) = SUM( v_i, dim=3 )
	88	vt_s (:,:) = SUM( v_s, dim=3 )
	89	at_i (:,:) = SUM( a_i, dim=3 )
	90	et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
	91	et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
[825]	92
[8233]	93	! MV MP 2016
	94	IF ( ln_pnd ) THEN
	95	at_ip(:,:) = SUM( a_ip, dim=3 )
	96	vt_ip(:,:) = SUM( v_ip, dim=3 )
	97	ENDIF
	98	! END MP 2016
	99
[7646]	100	! open water fraction
[921]	101	DO jj = 1, jpj
	102	DO ji = 1, jpi
[7646]	103	ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )
[921]	104	END DO
	105	END DO
[825]	106
[2715]	107	IF( kn > 1 ) THEN
[7646]	108
	109	! mean ice/snow thickness
	110	DO jj = 1, jpj
	111	DO ji = 1, jpi
	112	rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
	113	htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
	114	htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
	115	ENDDO
	116	ENDDO
	117
	118	! mean temperature (K), salinity and age
[7753]	119	smt_i(:,:) = 0._wp
	120	tm_i(:,:) = 0._wp
	121	tm_su(:,:) = 0._wp
[8239]	122	tm_si(:,:) = 0._wp
[7753]	123	om_i (:,:) = 0._wp
[921]	124	DO jl = 1, jpl
[7646]	125
[921]	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
[7646]	128	rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
	129	tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
[8239]	130	tm_si(ji,jj) = tm_si(ji,jj) + rswitch * ( t_si(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
	131	om_i (ji,jj) = om_i (ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )
[921]	132	END DO
	133	END DO
[7646]	134
[921]	135	DO jk = 1, nlay_i
[7646]	136	DO jj = 1, jpj
	137	DO ji = 1, jpi
	138	rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
	139	tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) &
	140	& / MAX( vt_i(ji,jj) , epsi10 )
	141	smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &
	142	& / MAX( vt_i(ji,jj) , epsi10 )
	143	END DO
	144	END DO
[921]	145	END DO
	146	END DO
[7646]	147	tm_i = tm_i + rt0
	148	tm_su = tm_su + rt0
[8239]	149	tm_si = tm_si + rt0
[2715]	150	!
	151	ENDIF
	152	!
[921]	153	END SUBROUTINE lim_var_agg
[825]	154
	155
[921]	156	SUBROUTINE lim_var_glo2eqv
	157	!!------------------------------------------------------------------
[2715]	158	!! * ROUTINE lim_var_glo2eqv *
[921]	159	!!
[2715]	160	!! ** Purpose : computes equivalent variables as function of global variables
	161	!! i.e. it turns VGLO into VEQV
[921]	162	!!------------------------------------------------------------------
[2715]	163	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[8325]	164	REAL(wp) :: ze_i, zaaa, zbbb, zccc, zdiscrim ! local scalars
	165	REAL(wp) :: ztmelts, ze_s, zfac1, zfac2 ! - -
[2715]	166	!!------------------------------------------------------------------
[825]	167
	168	!-------------------------------------------------------
	169	! Ice thickness, snow thickness, ice salinity, ice age
	170	!-------------------------------------------------------
	171	DO jl = 1, jpl
	172	DO jj = 1, jpj
	173	DO ji = 1, jpi
[5167]	174	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
	175	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
[6416]	176	END DO
	177	END DO
	178	END DO
	179	! Force the upper limit of ht_i to always be < hi_max (99 m).
	180	DO jj = 1, jpj
	181	DO ji = 1, jpi
	182	rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
	183	ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
	184	a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
	185	END DO
	186	END DO
	187
	188	DO jl = 1, jpl
	189	DO jj = 1, jpj
	190	DO ji = 1, jpi
	191	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
[5167]	192	ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
	193	o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
[825]	194	END DO
	195	END DO
	196	END DO
[6416]	197
[5123]	198	IF( nn_icesal == 2 )THEN
[921]	199	DO jl = 1, jpl
	200	DO jj = 1, jpj
	201	DO ji = 1, jpi
[5167]	202	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
	203	sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch
[5202]	204	! ! bounding salinity
	205	sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin )
[921]	206	END DO
[825]	207	END DO
	208	END DO
	209	ENDIF
	210
[2715]	211	CALL lim_var_salprof ! salinity profile
[825]	212
	213	!-------------------
	214	! Ice temperatures
	215	!-------------------
	216	DO jl = 1, jpl
[921]	217	DO jk = 1, nlay_i
	218	DO jj = 1, jpj
	219	DO ji = 1, jpi
[8325]	220	! ! Energy of melting e(S,T) [J.m-3]
[5134]	221	rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes
[8325]	222	ze_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)
[5123]	223	ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature
[2715]	224	!
	225	zaaa = cpic ! Conversion q(S,T) -> T (second order equation)
[8325]	226	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + ze_i * r1_rhoic - lfus
[5123]	227	zccc = lfus * (ztmelts-rt0)
[2715]	228	zdiscrim = SQRT( MAX(zbbbzbbb - 4._wpzaaa*zccc , 0._wp) )
[5123]	229	t_i(ji,jj,jk,jl) = rt0 + rswitch ( - zbbb - zdiscrim ) / ( 2.0 zaaa )
[5202]	230	t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) ) ! -100 < t_i < ztmelts
[921]	231	END DO
[825]	232	END DO
[921]	233	END DO
[825]	234	END DO
	235
	236	!--------------------
	237	! Snow temperatures
	238	!--------------------
[2715]	239	zfac1 = 1._wp / ( rhosn * cpic )
[825]	240	zfac2 = lfus / cpic
	241	DO jl = 1, jpl
[921]	242	DO jk = 1, nlay_s
	243	DO jj = 1, jpj
	244	DO ji = 1, jpi
	245	!Energy of melting q(S,T) [J.m-3]
[5134]	246	rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes
[8325]	247	ze_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
[2715]	248	!
[8325]	249	t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * ze_s + zfac2 )
[5202]	250	t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) ) ! -100 < t_s < rt0
[921]	251	END DO
[825]	252	END DO
[921]	253	END DO
[825]	254	END DO
	255
[8233]	256	! integrated values
[7753]	257	vt_i (:,:) = SUM( v_i, dim=3 )
	258	vt_s (:,:) = SUM( v_s, dim=3 )
	259	at_i (:,:) = SUM( a_i, dim=3 )
	260
[8233]	261	! MV MP 2016
	262	! probably should resum for melt ponds ???
	263
[2715]	264	!
[825]	265	END SUBROUTINE lim_var_glo2eqv
	266
	267
	268	SUBROUTINE lim_var_eqv2glo
[921]	269	!!------------------------------------------------------------------
[2715]	270	!! * ROUTINE lim_var_eqv2glo *
	271	!!
	272	!! ** Purpose : computes global variables as function of equivalent variables
	273	!! i.e. it turns VEQV into VGLO
[921]	274	!! ** Method :
	275	!!
[2715]	276	!! ** History : (01-2006) Martin Vancoppenolle, UCL-ASTR
[921]	277	!!------------------------------------------------------------------
[2715]	278	!
[921]	279	v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
	280	v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
	281	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
[2715]	282	!
[921]	283	END SUBROUTINE lim_var_eqv2glo
[825]	284
	285
[921]	286	SUBROUTINE lim_var_salprof
	287	!!------------------------------------------------------------------
[2715]	288	!! * ROUTINE lim_var_salprof *
[921]	289	!!
[2715]	290	!! ** Purpose : computes salinity profile in function of bulk salinity
	291	!!
[5123]	292	!! ** Method : If bulk salinity greater than zsi1,
[921]	293	!! the profile is assumed to be constant (S_inf)
[5123]	294	!! If bulk salinity lower than zsi0,
[921]	295	!! the profile is linear with 0 at the surface (S_zero)
[5123]	296	!! If it is between zsi0 and zsi1, it is a
[921]	297	!! alpha-weighted linear combination of s_inf and s_zero
	298	!!
[5123]	299	!! ** References : Vancoppenolle et al., 2007
[921]	300	!!------------------------------------------------------------------
[2715]	301	INTEGER :: ji, jj, jk, jl ! dummy loop index
[5123]	302	REAL(wp) :: zfac0, zfac1, zsal
	303	REAL(wp) :: zswi0, zswi01, zargtemp , zs_zero
[8373]	304	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_slope_s, zalpha
[5123]	305	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	306	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[2715]	307	!!------------------------------------------------------------------
[825]	308
	309	!---------------------------------------
	310	! Vertically constant, constant in time
	311	!---------------------------------------
[6470]	312	IF( nn_icesal == 1 ) THEN
[7753]	313	s_i (:,:,:,:) = rn_icesal
	314	sm_i(:,:,:) = rn_icesal
[6470]	315	ENDIF
[825]	316
	317	!-----------------------------------
	318	! Salinity profile, varying in time
	319	!-----------------------------------
[5123]	320	IF( nn_icesal == 2 ) THEN
[2715]	321	!
[7753]	322	DO jk = 1, nlay_i
	323	s_i(:,:,jk,:) = sm_i(:,:,:)
[2715]	324	END DO
	325	!
	326	DO jl = 1, jpl ! Slope of the linear profile
[825]	327	DO jj = 1, jpj
	328	DO ji = 1, jpi
[5167]	329	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
	330	z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
[2715]	331	END DO
	332	END DO
	333	END DO
	334	!
[5123]	335	zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf
	336	zfac1 = zsi1 / ( zsi1 - zsi0 )
[3625]	337	!
[7753]	338	zalpha(:,:,:) = 0._wp
[825]	339	DO jl = 1, jpl
	340	DO jj = 1, jpj
	341	DO ji = 1, jpi
[5123]	342	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
	343	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) )
	344	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
	345	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) )
	346	! If 2.sm_i GE sss_m then rswitch = 1
[4333]	347	! this is to force a constant salinity profile in the Baltic Sea
[5123]	348	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
	349	zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
	350	zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
[825]	351	END DO
	352	END DO
	353	END DO
[4161]	354
[5123]	355	! Computation of the profile
[825]	356	DO jl = 1, jpl
	357	DO jk = 1, nlay_i
	358	DO jj = 1, jpj
	359	DO ji = 1, jpi
[2715]	360	! ! linear profile with 0 at the surface
[5123]	361	zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
[2715]	362	! ! weighting the profile
	363	s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
[5202]	364	! ! bounding salinity
	365	s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
[5134]	366	END DO
	367	END DO
	368	END DO
	369	END DO
[3625]	370	!
[5123]	371	ENDIF ! nn_icesal
[825]	372
	373	!-------------------------------------------------------
	374	! Vertically varying salinity profile, constant in time
	375	!-------------------------------------------------------
[921]	376
[5123]	377	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
[2715]	378	!
[7753]	379	sm_i(:,:,:) = 2.30_wp
[2715]	380	!
[825]	381	DO jl = 1, jpl
	382	DO jk = 1, nlay_i
[5123]	383	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
[2715]	384	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
[7753]	385	s_i(:,:,jk,jl) = zsal
[2715]	386	END DO
	387	END DO
[3625]	388	!
[5123]	389	ENDIF ! nn_icesal
[2715]	390	!
[825]	391	END SUBROUTINE lim_var_salprof
	392
	393
	394	SUBROUTINE lim_var_bv
[921]	395	!!------------------------------------------------------------------
[2715]	396	!! * ROUTINE lim_var_bv *
[921]	397	!!
[2715]	398	!! ** Purpose : computes mean brine volume (%) in sea ice
	399	!!
[921]	400	!! ** Method : e = - 0.054 * S (ppt) / T (C)
	401	!!
[2715]	402	!! References : Vancoppenolle et al., JGR, 2007
[921]	403	!!------------------------------------------------------------------
[2715]	404	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	405	!!------------------------------------------------------------------
	406	!
[7753]	407	bvm_i(:,:) = 0._wp
	408	bv_i (:,:,:) = 0._wp
[5167]	409	DO jl = 1, jpl
[921]	410	DO jk = 1, nlay_i
	411	DO jj = 1, jpj
	412	DO ji = 1, jpi
[7646]	413	rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) )
	414	bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i &
	415	& / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
[921]	416	END DO
	417	END DO
	418	END DO
[7646]	419
	420	DO jj = 1, jpj
	421	DO ji = 1, jpi
	422	rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
	423	bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
	424	END DO
	425	END DO
[921]	426	END DO
[2715]	427	!
[921]	428	END SUBROUTINE lim_var_bv
[825]	429
	430
[8342]	431	SUBROUTINE lim_var_salprof1d
[825]	432	!!-------------------------------------------------------------------
	433	!! * ROUTINE lim_thd_salprof1d *
	434	!!
	435	!! ** Purpose : 1d computation of the sea ice salinity profile
[2715]	436	!! Works with 1d vectors and is used by thermodynamic modules
[825]	437	!!-------------------------------------------------------------------
[2715]	438	INTEGER :: ji, jk ! dummy loop indices
[5123]	439	REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars
	440	REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - -
[2715]	441	!
[8373]	442	REAL(wp), DIMENSION(jpij) :: z_slope_s
[5123]	443	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	444	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[2715]	445	!!---------------------------------------------------------------------
[825]	446
	447	!---------------------------------------
	448	! Vertically constant, constant in time
	449	!---------------------------------------
[5123]	450	IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal
[825]	451
	452	!------------------------------------------------------
	453	! Vertically varying salinity profile, varying in time
	454	!------------------------------------------------------
	455
[5123]	456	IF( nn_icesal == 2 ) THEN
[2715]	457	!
[8342]	458	DO ji = 1, nidx ! Slope of the linear profile zs_zero
[5167]	459	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
	460	z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
[2715]	461	END DO
[825]	462
	463	! Weighting factor between zs_zero and zs_inf
	464	!---------------------------------------------
[5123]	465	zfac0 = 1._wp / ( zsi0 - zsi1 )
	466	zfac1 = zsi1 / ( zsi1 - zsi0 )
[825]	467	DO jk = 1, nlay_i
[8342]	468	DO ji = 1, nidx
[5123]	469	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
	470	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) )
	471	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
	472	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) )
	473	! if 2.sm_i GE sss_m then rswitch = 1
[4333]	474	! this is to force a constant salinity profile in the Baltic Sea
[8360]	475	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_1d(ji) ) )
[2715]	476	!
[5123]	477	zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch )
[2715]	478	!
[5123]	479	zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
[825]	480	! weighting the profile
[4872]	481	s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
[5202]	482	! bounding salinity
	483	s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
[5123]	484	END DO
	485	END DO
[825]	486
[5123]	487	ENDIF
[825]	488
	489	!-------------------------------------------------------
	490	! Vertically varying salinity profile, constant in time
	491	!-------------------------------------------------------
	492
[5123]	493	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
[2715]	494	!
[4872]	495	sm_i_1d(:) = 2.30_wp
[2715]	496	!
	497	DO jk = 1, nlay_i
[5123]	498	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	499	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
[8342]	500	DO ji = 1, nidx
[4872]	501	s_i_1d(ji,jk) = zsal
[2715]	502	END DO
	503	END DO
	504	!
	505	ENDIF
	506	!
[825]	507	END SUBROUTINE lim_var_salprof1d
	508
[5123]	509	SUBROUTINE lim_var_zapsmall
	510	!!-------------------------------------------------------------------
	511	!! * ROUTINE lim_var_zapsmall *
	512	!!
	513	!! ** Purpose : Remove too small sea ice areas and correct fluxes
	514	!!
	515	!! history : LIM3.5 - 01-2014 (C. Rousset) original code
	516	!!-------------------------------------------------------------------
	517	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[8233]	518	REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp
[5123]	519	!!-------------------------------------------------------------------
	520	DO jl = 1, jpl
	521
	522	!-----------------------------------------------------------------
	523	! Zap ice energy and use ocean heat to melt ice
	524	!-----------------------------------------------------------------
	525	DO jk = 1, nlay_i
	526	DO jj = 1 , jpj
	527	DO ji = 1 , jpi
	528	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
[5202]	529	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
	530	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
	531	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
[5123]	532	zei = e_i(ji,jj,jk,jl)
	533	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
	534	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
	535	! update exchanges with ocean
	536	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
	537	END DO
	538	END DO
	539	END DO
	540
	541	DO jj = 1 , jpj
	542	DO ji = 1 , jpi
	543	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
[5202]	544	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
	545	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
	546	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
[5123]	547	zsal = smv_i(ji,jj, jl)
	548	zvi = v_i (ji,jj, jl)
	549	zvs = v_s (ji,jj, jl)
	550	zes = e_s (ji,jj,1,jl)
[8233]	551	IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl)
[5123]	552	!-----------------------------------------------------------------
	553	! Zap snow energy
	554	!-----------------------------------------------------------------
	555	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
	556	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
	557
	558	!-----------------------------------------------------------------
	559	! zap ice and snow volume, add water and salt to ocean
	560	!-----------------------------------------------------------------
	561	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
	562	a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch
	563	v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch
	564	v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch
	565	t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
	566	oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
	567	smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
	568
[8360]	569	ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
	570	ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
	571
[8233]	572	! MV MP 2016
	573	IF ( ln_pnd ) THEN
	574	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
	575	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
	576	IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice
	577	ENDIF
	578	! END MV MP 2016
	579
[5123]	580	! update exchanges with ocean
	581	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
	582	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
[8341]	583	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
[5123]	584	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
	585	END DO
	586	END DO
	587	END DO
	588
	589	! to be sure that at_i is the sum of a_i(jl)
[7753]	590	at_i (:,:) = 0._wp
[5123]	591	DO jl = 1, jpl
[7753]	592	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[5123]	593	END DO
	594
	595	! open water = 1 if at_i=0
	596	DO jj = 1, jpj
	597	DO ji = 1, jpi
	598	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
	599	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
	600	END DO
	601	END DO
	602
	603	!
	604	END SUBROUTINE lim_var_zapsmall
	605
	606	SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
	607	!!------------------------------------------------------------------
	608	!! * ROUTINE lim_var_itd *
	609	!!
	610	!! ** Purpose : converting 1-cat ice to multiple ice categories
	611	!!
	612	!! ice thickness distribution follows a gaussian law
	613	!! around the concentration of the most likely ice thickness
[8413]	614	!! (similar as iceistate.F90)
[5123]	615	!!
	616	!! ** Method: Iterative procedure
	617	!!
	618	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	619	!!
	620	!! 2) Check whether the distribution conserves area and volume, positivity and
	621	!! category boundaries
	622	!!
	623	!! 3) If not (input ice is too thin), the last category is empty and
	624	!! the number of categories is reduced (jpl-1)
	625	!!
	626	!! 4) Iterate until ok (SUM(itest(:) = 4)
	627	!!
	628	!! ** Arguments : zhti: 1-cat ice thickness
	629	!! zhts: 1-cat snow depth
	630	!! zai : 1-cat ice concentration
	631	!!
	632	!! ** Output : jpl-cat
	633	!!
	634	!! (Example of application: BDY forcings when input are cell averaged)
	635	!!
	636	!!-------------------------------------------------------------------
	637	!! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code
	638	!! 2014 (C. Rousset) Rewriting
	639	!!-------------------------------------------------------------------
	640	!! Local variables
	641	INTEGER :: ji, jk, jl ! dummy loop indices
	642	INTEGER :: ijpij, i_fill, jl0
[7646]	643	REAL(wp) :: zarg, zV, zconv, zdh, zdv
[5123]	644	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables
	645	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables
[8373]	646	INTEGER , DIMENSION(4) :: itest
[5123]	647
	648	!--------------------------------------------------------------------
	649	! initialisation of variables
	650	!--------------------------------------------------------------------
	651	ijpij = SIZE(zhti,1)
	652	zht_i(1:ijpij,1:jpl) = 0._wp
	653	zht_s(1:ijpij,1:jpl) = 0._wp
	654	za_i (1:ijpij,1:jpl) = 0._wp
	655
	656	! ----------------------------------------
	657	! distribution over the jpl ice categories
	658	! ----------------------------------------
	659	DO ji = 1, ijpij
	660
	661	IF( zhti(ji) > 0._wp ) THEN
	662
[7646]	663	! find which category (jl0) the input ice thickness falls into
	664	jl0 = jpl
	665	DO jl = 1, jpl
	666	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
	667	jl0 = jl
	668	CYCLE
	669	ENDIF
	670	END DO
	671
	672	! initialisation of tests
	673	itest(:) = 0
[5123]	674
[7646]	675	i_fill = jpl + 1 !====================================
	676	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	677	! iteration !====================================
	678	i_fill = i_fill - 1
	679
	680	! initialisation of ice variables for each try
	681	zht_i(ji,1:jpl) = 0._wp
	682	za_i (ji,1:jpl) = 0._wp
[7813]	683	itest(:) = 0
[7646]	684
	685	! *** case very thin ice: fill only category 1
	686	IF ( i_fill == 1 ) THEN
	687	zht_i(ji,1) = zhti(ji)
	688	za_i (ji,1) = zai (ji)
	689
	690	! *** case ice is thicker: fill categories >1
	691	ELSE
[5123]	692
[7646]	693	! Fill ice thicknesses in the (i_fill-1) cat by hmean
	694	DO jl = 1, i_fill - 1
	695	zht_i(ji,jl) = hi_mean(jl)
	696	END DO
	697
	698	! Concentrations in the (i_fill-1) categories
	699	za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
	700	DO jl = 1, i_fill - 1
	701	IF ( jl /= jl0 ) THEN
	702	zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
	703	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
	704	ENDIF
	705	END DO
	706
	707	! Concentration in the last (i_fill) category
	708	za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
	709
	710	! Ice thickness in the last (i_fill) category
	711	zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
	712	zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
	713
	714	! clem: correction if concentration of upper cat is greater than lower cat
	715	! (it should be a gaussian around jl0 but sometimes it is not)
	716	IF ( jl0 /= jpl ) THEN
	717	DO jl = jpl, jl0+1, -1
	718	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
	719	zdv = zht_i(ji,jl) * za_i(ji,jl)
	720	zht_i(ji,jl ) = 0._wp
	721	za_i (ji,jl ) = 0._wp
	722	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
	723	END IF
	724	ENDDO
[5123]	725	ENDIF
	726
[7646]	727	ENDIF ! case ice is thick or thin
[7813]	728
[7646]	729	!---------------------
	730	! Compatibility tests
	731	!---------------------
	732	! Test 1: area conservation
	733	zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
	734	IF ( zconv < epsi06 ) itest(1) = 1
	735
	736	! Test 2: volume conservation
	737	zconv = ABS( zhti(ji)zai(ji) - SUM( za_i(ji,1:jpl)zht_i(ji,1:jpl) ) )
	738	IF ( zconv < epsi06 ) itest(2) = 1
[5123]	739
[7646]	740	! Test 3: thickness of the last category is in-bounds ?
	741	IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
[5123]	742
[7646]	743	! Test 4: positivity of ice concentrations
	744	itest(4) = 1
	745	DO jl = 1, i_fill
	746	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
	747	END DO
	748	! !============================
	749	END DO ! end iteration on categories
	750	! !============================
[5123]	751	ENDIF ! if zhti > 0
	752	END DO ! i loop
[7813]	753
[5123]	754	! ------------------------------------------------
	755	! Adding Snow in each category where za_i is not 0
	756	! ------------------------------------------------
	757	DO jl = 1, jpl
	758	DO ji = 1, ijpij
	759	IF( za_i(ji,jl) > 0._wp ) THEN
	760	zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
	761	! In case snow load is in excess that would lead to transformation from snow to ice
	762	! Then, transfer the snow excess into the ice (different from limthd_dh)
	763	zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )
	764	! recompute ht_i, ht_s avoiding out of bounds values
	765	zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
	766	zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
	767	ENDIF
	768	ENDDO
	769	ENDDO
	770	!
	771	END SUBROUTINE lim_var_itd
	772
	773
[825]	774	#else
[2715]	775	!!----------------------------------------------------------------------
	776	!! Default option Dummy module NO LIM3 sea-ice model
	777	!!----------------------------------------------------------------------
[825]	778	CONTAINS
	779	SUBROUTINE lim_var_agg ! Empty routines
	780	END SUBROUTINE lim_var_agg
	781	SUBROUTINE lim_var_glo2eqv ! Empty routines
	782	END SUBROUTINE lim_var_glo2eqv
	783	SUBROUTINE lim_var_eqv2glo ! Empty routines
	784	END SUBROUTINE lim_var_eqv2glo
	785	SUBROUTINE lim_var_salprof ! Empty routines
	786	END SUBROUTINE lim_var_salprof
	787	SUBROUTINE lim_var_bv ! Emtpy routines
[921]	788	END SUBROUTINE lim_var_bv
[825]	789	SUBROUTINE lim_var_salprof1d ! Emtpy routines
	790	END SUBROUTINE lim_var_salprof1d
[5123]	791	SUBROUTINE lim_var_zapsmall
	792	END SUBROUTINE lim_var_zapsmall
	793	SUBROUTINE lim_var_itd
	794	END SUBROUTINE lim_var_itd
[2715]	795	#endif
[825]	796
[2715]	797	!!======================================================================
[834]	798	END MODULE limvar

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

Download in other formats: