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/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5079

Last change on this file since 5079 was 5079, checked in by clem, 9 years ago
LIM3: solve ticket #1421
Property svn:keywords set to `Id`
File size: 32.4 KB

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

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

Download in other formats: