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/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 6584

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

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

Download in other formats: