Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

limvar.F90 @ 7910

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

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7910

Download in other formats: