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

Last change on this file since 8316 was 8316, checked in by clem, 7 years ago
STEP2 (3): remove obsolete features
Property svn:keywords set to `Id`
File size: 35.1 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)
	43	USE in_out_manager ! I/O manager
	44	USE lib_mpp ! MPP library
	45	USE wrk_nemo ! work arrays
	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
	164	REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars
[4990]	165	REAL(wp) :: ztmelts, zq_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
[2715]	220	! ! Energy of melting q(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
[5123]	222	zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)
	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)
[5123]	226	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus
	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
[5123]	247	zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
[2715]	248	!
[5123]	249	t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_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
	304	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha
	305	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	306	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[2715]	307	!!------------------------------------------------------------------
[825]	308
[3294]	309	CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
[825]	310
	311	!---------------------------------------
	312	! Vertically constant, constant in time
	313	!---------------------------------------
[6470]	314	IF( nn_icesal == 1 ) THEN
[7753]	315	s_i (:,:,:,:) = rn_icesal
	316	sm_i(:,:,:) = rn_icesal
[6470]	317	ENDIF
[825]	318
	319	!-----------------------------------
	320	! Salinity profile, varying in time
	321	!-----------------------------------
[5123]	322	IF( nn_icesal == 2 ) THEN
[2715]	323	!
[7753]	324	DO jk = 1, nlay_i
	325	s_i(:,:,jk,:) = sm_i(:,:,:)
[2715]	326	END DO
	327	!
	328	DO jl = 1, jpl ! Slope of the linear profile
[825]	329	DO jj = 1, jpj
	330	DO ji = 1, jpi
[5167]	331	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
	332	z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
[2715]	333	END DO
	334	END DO
	335	END DO
	336	!
[5123]	337	zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf
	338	zfac1 = zsi1 / ( zsi1 - zsi0 )
[3625]	339	!
[7753]	340	zalpha(:,:,:) = 0._wp
[825]	341	DO jl = 1, jpl
	342	DO jj = 1, jpj
	343	DO ji = 1, jpi
[5123]	344	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
	345	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) )
	346	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
	347	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) )
	348	! If 2.sm_i GE sss_m then rswitch = 1
[4333]	349	! this is to force a constant salinity profile in the Baltic Sea
[5123]	350	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
	351	zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
	352	zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
[825]	353	END DO
	354	END DO
	355	END DO
[4161]	356
[5123]	357	! Computation of the profile
[825]	358	DO jl = 1, jpl
	359	DO jk = 1, nlay_i
	360	DO jj = 1, jpj
	361	DO ji = 1, jpi
[2715]	362	! ! linear profile with 0 at the surface
[5123]	363	zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
[2715]	364	! ! weighting the profile
	365	s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
[5202]	366	! ! bounding salinity
	367	s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
[5134]	368	END DO
	369	END DO
	370	END DO
	371	END DO
[3625]	372	!
[5123]	373	ENDIF ! nn_icesal
[825]	374
	375	!-------------------------------------------------------
	376	! Vertically varying salinity profile, constant in time
	377	!-------------------------------------------------------
[921]	378
[5123]	379	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
[2715]	380	!
[7753]	381	sm_i(:,:,:) = 2.30_wp
[2715]	382	!
[825]	383	DO jl = 1, jpl
	384	DO jk = 1, nlay_i
[5123]	385	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
[2715]	386	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
[7753]	387	s_i(:,:,jk,jl) = zsal
[2715]	388	END DO
	389	END DO
[3625]	390	!
[5123]	391	ENDIF ! nn_icesal
[2715]	392	!
[3294]	393	CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
[2715]	394	!
[825]	395	END SUBROUTINE lim_var_salprof
	396
	397
	398	SUBROUTINE lim_var_bv
[921]	399	!!------------------------------------------------------------------
[2715]	400	!! * ROUTINE lim_var_bv *
[921]	401	!!
[2715]	402	!! ** Purpose : computes mean brine volume (%) in sea ice
	403	!!
[921]	404	!! ** Method : e = - 0.054 * S (ppt) / T (C)
	405	!!
[2715]	406	!! References : Vancoppenolle et al., JGR, 2007
[921]	407	!!------------------------------------------------------------------
[2715]	408	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	409	!!------------------------------------------------------------------
	410	!
[7753]	411	bvm_i(:,:) = 0._wp
	412	bv_i (:,:,:) = 0._wp
[5167]	413	DO jl = 1, jpl
[921]	414	DO jk = 1, nlay_i
	415	DO jj = 1, jpj
	416	DO ji = 1, jpi
[7646]	417	rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) )
	418	bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i &
	419	& / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
[921]	420	END DO
	421	END DO
	422	END DO
[7646]	423
	424	DO jj = 1, jpj
	425	DO ji = 1, jpi
	426	rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
	427	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 )
	428	END DO
	429	END DO
[921]	430	END DO
[2715]	431	!
[921]	432	END SUBROUTINE lim_var_bv
[825]	433
	434
[2715]	435	SUBROUTINE lim_var_salprof1d( kideb, kiut )
[825]	436	!!-------------------------------------------------------------------
	437	!! * ROUTINE lim_thd_salprof1d *
	438	!!
	439	!! ** Purpose : 1d computation of the sea ice salinity profile
[2715]	440	!! Works with 1d vectors and is used by thermodynamic modules
[825]	441	!!-------------------------------------------------------------------
[2715]	442	INTEGER, INTENT(in) :: kideb, kiut ! thickness category index
	443	!
	444	INTEGER :: ji, jk ! dummy loop indices
[5123]	445	INTEGER :: ii, ij ! local integers
	446	REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars
	447	REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - -
[2715]	448	!
	449	REAL(wp), POINTER, DIMENSION(:) :: z_slope_s
[5123]	450	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	451	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
[2715]	452	!!---------------------------------------------------------------------
[825]	453
[3294]	454	CALL wrk_alloc( jpij, z_slope_s )
[825]	455
	456	!---------------------------------------
	457	! Vertically constant, constant in time
	458	!---------------------------------------
[5123]	459	IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal
[825]	460
	461	!------------------------------------------------------
	462	! Vertically varying salinity profile, varying in time
	463	!------------------------------------------------------
	464
[5123]	465	IF( nn_icesal == 2 ) THEN
[2715]	466	!
	467	DO ji = kideb, kiut ! Slope of the linear profile zs_zero
[5167]	468	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
	469	z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
[2715]	470	END DO
[825]	471
	472	! Weighting factor between zs_zero and zs_inf
	473	!---------------------------------------------
[5123]	474	zfac0 = 1._wp / ( zsi0 - zsi1 )
	475	zfac1 = zsi1 / ( zsi1 - zsi0 )
[825]	476	DO jk = 1, nlay_i
	477	DO ji = kideb, kiut
[4161]	478	ii = MOD( npb(ji) - 1 , jpi ) + 1
	479	ij = ( npb(ji) - 1 ) / jpi + 1
[5123]	480	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
	481	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) )
	482	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
	483	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) )
	484	! if 2.sm_i GE sss_m then rswitch = 1
[4333]	485	! this is to force a constant salinity profile in the Baltic Sea
[5123]	486	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
[2715]	487	!
[5123]	488	zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch )
[2715]	489	!
[5123]	490	zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
[825]	491	! weighting the profile
[4872]	492	s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
[5202]	493	! bounding salinity
	494	s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
[5123]	495	END DO
	496	END DO
[825]	497
[5123]	498	ENDIF
[825]	499
	500	!-------------------------------------------------------
	501	! Vertically varying salinity profile, constant in time
	502	!-------------------------------------------------------
	503
[5123]	504	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
[2715]	505	!
[4872]	506	sm_i_1d(:) = 2.30_wp
[2715]	507	!
	508	DO jk = 1, nlay_i
[5123]	509	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	510	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
[2715]	511	DO ji = kideb, kiut
[4872]	512	s_i_1d(ji,jk) = zsal
[2715]	513	END DO
	514	END DO
	515	!
	516	ENDIF
	517	!
[3294]	518	CALL wrk_dealloc( jpij, z_slope_s )
[2715]	519	!
[825]	520	END SUBROUTINE lim_var_salprof1d
	521
[5123]	522	SUBROUTINE lim_var_zapsmall
	523	!!-------------------------------------------------------------------
	524	!! * ROUTINE lim_var_zapsmall *
	525	!!
	526	!! ** Purpose : Remove too small sea ice areas and correct fluxes
	527	!!
	528	!! history : LIM3.5 - 01-2014 (C. Rousset) original code
	529	!!-------------------------------------------------------------------
	530	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[8233]	531	REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp
[5123]	532	!!-------------------------------------------------------------------
[7753]	533	at_i (:,:) = 0._wp
[5123]	534	DO jl = 1, jpl
[7753]	535	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[5123]	536	END DO
	537
	538	DO jl = 1, jpl
	539
	540	!-----------------------------------------------------------------
	541	! Zap ice energy and use ocean heat to melt ice
	542	!-----------------------------------------------------------------
	543	DO jk = 1, nlay_i
	544	DO jj = 1 , jpj
	545	DO ji = 1 , jpi
	546	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
	547	rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch
[5202]	548	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
	549	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
	550	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
[5123]	551	zei = e_i(ji,jj,jk,jl)
	552	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
	553	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
	554	! update exchanges with ocean
	555	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
	556	END DO
	557	END DO
	558	END DO
	559
	560	DO jj = 1 , jpj
	561	DO ji = 1 , jpi
	562	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
	563	rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch
[5202]	564	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
	565	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
	566	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
[5123]	567	zsal = smv_i(ji,jj, jl)
	568	zvi = v_i (ji,jj, jl)
	569	zvs = v_s (ji,jj, jl)
	570	zes = e_s (ji,jj,1,jl)
[8233]	571	IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl)
[5123]	572	!-----------------------------------------------------------------
	573	! Zap snow energy
	574	!-----------------------------------------------------------------
	575	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
	576	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
	577
	578	!-----------------------------------------------------------------
	579	! zap ice and snow volume, add water and salt to ocean
	580	!-----------------------------------------------------------------
	581	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
	582	a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch
	583	v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch
	584	v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch
	585	t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
	586	oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
	587	smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
	588
[8233]	589	! MV MP 2016
	590	IF ( ln_pnd ) THEN
	591	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
	592	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
	593	IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice
	594	ENDIF
	595	! END MV MP 2016
	596
[5123]	597	! update exchanges with ocean
	598	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
	599	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	600	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	601	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
	602	END DO
	603	END DO
	604	END DO
	605
	606	! to be sure that at_i is the sum of a_i(jl)
[7753]	607	at_i (:,:) = 0._wp
[5123]	608	DO jl = 1, jpl
[7753]	609	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[5123]	610	END DO
	611
	612	! open water = 1 if at_i=0
	613	DO jj = 1, jpj
	614	DO ji = 1, jpi
	615	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
	616	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
	617	END DO
	618	END DO
	619
	620	!
	621	END SUBROUTINE lim_var_zapsmall
	622
	623	SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
	624	!!------------------------------------------------------------------
	625	!! * ROUTINE lim_var_itd *
	626	!!
	627	!! ** Purpose : converting 1-cat ice to multiple ice categories
	628	!!
	629	!! ice thickness distribution follows a gaussian law
	630	!! around the concentration of the most likely ice thickness
	631	!! (similar as limistate.F90)
	632	!!
	633	!! ** Method: Iterative procedure
	634	!!
	635	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	636	!!
	637	!! 2) Check whether the distribution conserves area and volume, positivity and
	638	!! category boundaries
	639	!!
	640	!! 3) If not (input ice is too thin), the last category is empty and
	641	!! the number of categories is reduced (jpl-1)
	642	!!
	643	!! 4) Iterate until ok (SUM(itest(:) = 4)
	644	!!
	645	!! ** Arguments : zhti: 1-cat ice thickness
	646	!! zhts: 1-cat snow depth
	647	!! zai : 1-cat ice concentration
	648	!!
	649	!! ** Output : jpl-cat
	650	!!
	651	!! (Example of application: BDY forcings when input are cell averaged)
	652	!!
	653	!!-------------------------------------------------------------------
	654	!! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code
	655	!! 2014 (C. Rousset) Rewriting
	656	!!-------------------------------------------------------------------
	657	!! Local variables
	658	INTEGER :: ji, jk, jl ! dummy loop indices
	659	INTEGER :: ijpij, i_fill, jl0
[7646]	660	REAL(wp) :: zarg, zV, zconv, zdh, zdv
[5123]	661	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables
	662	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables
	663	INTEGER , POINTER, DIMENSION(:) :: itest
	664
	665	CALL wrk_alloc( 4, itest )
	666	!--------------------------------------------------------------------
	667	! initialisation of variables
	668	!--------------------------------------------------------------------
	669	ijpij = SIZE(zhti,1)
	670	zht_i(1:ijpij,1:jpl) = 0._wp
	671	zht_s(1:ijpij,1:jpl) = 0._wp
	672	za_i (1:ijpij,1:jpl) = 0._wp
	673
	674	! ----------------------------------------
	675	! distribution over the jpl ice categories
	676	! ----------------------------------------
	677	DO ji = 1, ijpij
	678
	679	IF( zhti(ji) > 0._wp ) THEN
	680
[7646]	681	! find which category (jl0) the input ice thickness falls into
	682	jl0 = jpl
	683	DO jl = 1, jpl
	684	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
	685	jl0 = jl
	686	CYCLE
	687	ENDIF
	688	END DO
	689
	690	! initialisation of tests
	691	itest(:) = 0
[5123]	692
[7646]	693	i_fill = jpl + 1 !====================================
	694	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	695	! iteration !====================================
	696	i_fill = i_fill - 1
	697
	698	! initialisation of ice variables for each try
	699	zht_i(ji,1:jpl) = 0._wp
	700	za_i (ji,1:jpl) = 0._wp
[7813]	701	itest(:) = 0
[7646]	702
	703	! *** case very thin ice: fill only category 1
	704	IF ( i_fill == 1 ) THEN
	705	zht_i(ji,1) = zhti(ji)
	706	za_i (ji,1) = zai (ji)
	707
	708	! *** case ice is thicker: fill categories >1
	709	ELSE
[5123]	710
[7646]	711	! Fill ice thicknesses in the (i_fill-1) cat by hmean
	712	DO jl = 1, i_fill - 1
	713	zht_i(ji,jl) = hi_mean(jl)
	714	END DO
	715
	716	! Concentrations in the (i_fill-1) categories
	717	za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
	718	DO jl = 1, i_fill - 1
	719	IF ( jl /= jl0 ) THEN
	720	zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
	721	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
	722	ENDIF
	723	END DO
	724
	725	! Concentration in the last (i_fill) category
	726	za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
	727
	728	! Ice thickness in the last (i_fill) category
	729	zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
	730	zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
	731
	732	! clem: correction if concentration of upper cat is greater than lower cat
	733	! (it should be a gaussian around jl0 but sometimes it is not)
	734	IF ( jl0 /= jpl ) THEN
	735	DO jl = jpl, jl0+1, -1
	736	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
	737	zdv = zht_i(ji,jl) * za_i(ji,jl)
	738	zht_i(ji,jl ) = 0._wp
	739	za_i (ji,jl ) = 0._wp
	740	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
	741	END IF
	742	ENDDO
[5123]	743	ENDIF
	744
[7646]	745	ENDIF ! case ice is thick or thin
[7813]	746
[7646]	747	!---------------------
	748	! Compatibility tests
	749	!---------------------
	750	! Test 1: area conservation
	751	zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
	752	IF ( zconv < epsi06 ) itest(1) = 1
	753
	754	! Test 2: volume conservation
	755	zconv = ABS( zhti(ji)zai(ji) - SUM( za_i(ji,1:jpl)zht_i(ji,1:jpl) ) )
	756	IF ( zconv < epsi06 ) itest(2) = 1
[5123]	757
[7646]	758	! Test 3: thickness of the last category is in-bounds ?
	759	IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
[5123]	760
[7646]	761	! Test 4: positivity of ice concentrations
	762	itest(4) = 1
	763	DO jl = 1, i_fill
	764	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
	765	END DO
	766	! !============================
	767	END DO ! end iteration on categories
	768	! !============================
[5123]	769	ENDIF ! if zhti > 0
	770	END DO ! i loop
[7813]	771
[5123]	772	! ------------------------------------------------
	773	! Adding Snow in each category where za_i is not 0
	774	! ------------------------------------------------
	775	DO jl = 1, jpl
	776	DO ji = 1, ijpij
	777	IF( za_i(ji,jl) > 0._wp ) THEN
	778	zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
	779	! In case snow load is in excess that would lead to transformation from snow to ice
	780	! Then, transfer the snow excess into the ice (different from limthd_dh)
	781	zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )
	782	! recompute ht_i, ht_s avoiding out of bounds values
	783	zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
	784	zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
	785	ENDIF
	786	ENDDO
	787	ENDDO
	788
	789	CALL wrk_dealloc( 4, itest )
	790	!
	791	END SUBROUTINE lim_var_itd
	792
	793
[825]	794	#else
[2715]	795	!!----------------------------------------------------------------------
	796	!! Default option Dummy module NO LIM3 sea-ice model
	797	!!----------------------------------------------------------------------
[825]	798	CONTAINS
	799	SUBROUTINE lim_var_agg ! Empty routines
	800	END SUBROUTINE lim_var_agg
	801	SUBROUTINE lim_var_glo2eqv ! Empty routines
	802	END SUBROUTINE lim_var_glo2eqv
	803	SUBROUTINE lim_var_eqv2glo ! Empty routines
	804	END SUBROUTINE lim_var_eqv2glo
	805	SUBROUTINE lim_var_salprof ! Empty routines
	806	END SUBROUTINE lim_var_salprof
	807	SUBROUTINE lim_var_bv ! Emtpy routines
[921]	808	END SUBROUTINE lim_var_bv
[825]	809	SUBROUTINE lim_var_salprof1d ! Emtpy routines
	810	END SUBROUTINE lim_var_salprof1d
[5123]	811	SUBROUTINE lim_var_zapsmall
	812	END SUBROUTINE lim_var_zapsmall
	813	SUBROUTINE lim_var_itd
	814	END SUBROUTINE lim_var_itd
[2715]	815	#endif
[825]	816
[2715]	817	!!======================================================================
[834]	818	END MODULE limvar

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

Download in other formats: