Context Navigation

ice_sal_adv.f @ 59

Last change on this file since 59 was 20, checked in by vancop, 8 years ago
Ludivine source files
File size: 13.2 KB

Rev	Line
[10]	1	SUBROUTINE ice_sal_adv(nlay_i,kideb,kiut)
	2
	3	!!------------------------------------------------------------------
	4	!! * ROUTINE ice_sal_adv *
	5	!!
	6	!! ** Purpose :
	7	!! This routine computes new salinities in the ice
	8	!!
	9	!! ** Method : Vertical salinity profile computation
	10	!! Resolves brine transport equation
	11	!!
	12	!! ** Steps
	13	!!
	14	!! ** Arguments
	15	!!
	16	!! ** Inputs / Outputs
	17	!!
	18	!! ** External
	19	!!
	20	!! ** References : Vancop. et al., 2008
	21	!!
	22	!! ** History :
	23	!! (06-2003) Martin Vancop. LIM1D
	24	!! (06-2008) Martin Vancop. BIO-LIM
	25	!! (09-2008) Martin Vancop. Explicit gravity drainage
	26	!!
	27	!!------------------------------------------------------------------
	28
	29	USE lib_fortran
	30
	31	INCLUDE 'type.com'
	32	INCLUDE 'para.com'
	33	INCLUDE 'const.com'
	34	INCLUDE 'ice.com'
	35	INCLUDE 'thermo.com'
	36
	37	REAL(8), DIMENSION(nlay_i) ::
	38	& z_ms_i , !: mass of salt times thickness
	39	& z_sbr_i !: brine salinity
	40
	41	REAL(8), DIMENSION(nlay_i) :: !: dummy factors for tracer equation
	42	& za , !: all
	43	& zb , !: gravity drainage
	44	& zc , !: upward advective flow
	45	& zRae , !: effective Ra
	46	& ze , !: downward advective flow
	47	& zind , !: independent term in the tridiag system
	48	& zindtbis , !:
	49	& zdiagbis !:
	50
[20]	51
[10]	52	REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation
	53	& ztrid !: tridiagonal matrix
	54
	55	REAL(8) ::
	56	& zdummy1 , !: dummy factors
	57	& zdummy2 , !:
	58	& zdummy3 , !:
	59	& zswitchs , !: switch for summer drainage
	60	& zeps = 1.0e-20 !: numerical limit
	61
	62	! Rayleigh number computation
	63	REAL(8) ::
	64	& ze_i_min , !: minimum brine volume
	65	& zcp , !: temporary scalar for sea ice specific heat
	66	& zk , !: temporary scalar for sea ice thermal conductivity
	67	& zalphara !: multiplicator for diffusivity
	68
	69	REAL(8), DIMENSION(nlay_i) ::
	70	& zsigma , !: brine salinity at layer interfaces
	71	& zperm , !: permeability
	72	& zpermin , !: minimum permeability
	73	& zperm_eff , !: minimum permeability
	74	& zrhodiff , !: density difference
	75	& zlevel , !: height of the water column
	76	& zthdiff !: thermal diffusivity
[20]	77
[10]	78	REAL(8), DIMENSION(nlay_i+1) ::
	79	& z_sbr_int !: brine salinity at layer interfaces
	80
	81	INTEGER ::
	82	& layer2 , !: layer loop index
	83	& indtr !: index of tridiagonal system
	84
	85	CHARACTER(len=4) ::
	86	& bc = 'conc' !: Boundary condition 'conc' or 'flux'
	87
	88	REAL(8) ::
	89	& z_ms_i_ini , !: initial mass of salt
	90	& z_ms_i_fin , !: final mass of salt
	91	& z_fs_b , !: basal flux of salt
	92	& z_fs_su , !: surface flux of salt
	93	& z_dms_i !: mass variation
	94
	95	LOGICAL ::
	96	& ln_write ,
	97	& ln_con ,
	98	& ln_sal
	99
	100	ln_write = .TRUE. ! write outputs
	101	ln_con = .TRUE. ! conservation check
	102	ln_sal = .TRUE. ! compute salinity variations or not
	103
	104	IF ( ln_write ) THEN
	105	WRITE(numout,*)
	106	WRITE(numout,) ' * ice_sal_adv : '
	107	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ '
	108	WRITE(numout,*) ' ln_sal = ', ln_sal
	109	WRITE(numout,*) ' ln_grd = ', ln_grd
	110	WRITE(numout,*) ' ln_flu = ', ln_flu
	111	WRITE(numout,*) ' ln_flo = ', ln_flo
	112	WRITE(numout,*) ' c_gravdr = ', c_gravdr
	113	WRITE(numout,*) ' c_sbr = ', c_sbr
	114	WRITE(numout,*) ' c_perm = ', c_perm
	115	WRITE(numout,*) ' c_permeff= ', c_permeff
	116	ENDIF
	117
	118	ji = 1
	119
	120	IF ( ln_sal ) THEN
	121	!
	122	!------------------------------------------------------------------------------\|
	123	! 1) Initialization
	124	!------------------------------------------------------------------------------\|
	125	!
	126	IF ( ln_write ) THEN
	127	WRITE(numout,*) ' - Initialization ... '
	128	ENDIF
	129
	130	! brine diffusivity
	131	diff_br(:) = 0.0
	132
	133	! Darcy velocity
	134	w_adv_br(:) = 0.0
	135
	136	!--------------------
	137	! Conservation check
	138	!--------------------
	139	IF ( ln_con ) THEN
	140	CALL ice_sal_column( kideb , kiut , z_ms_i_ini ,
	141	& s_i_b(1,1:nlay_i),
	142	& deltaz_i_phy, nlay_i, .FALSE. )
	143	ENDIF ! ln_con
	144
	145	IF ( ln_write ) THEN
	146	WRITE(numout,*) ' '
	147	WRITE(numout,*) ' nlay_i : ', nlay_i
	148	WRITE(numout,*) ' kideb : ', kideb
	149	WRITE(numout,*) ' kiut : ', kiut
	150	WRITE(numout,*) ' '
	151	WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer),
	152	& layer = 1, nlay_i )
	153	WRITE(numout,*) ' z_i_phy : ', ( z_i_phy(layer),
	154	& layer = 1, nlay_i )
	155	WRITE(numout,*) ' s_i_b : ', ( s_i_b (ji,layer),
	156	& layer = 1, nlay_i )
	157	WRITE(numout,*) ' t_i_b : ', ( t_i_b (ji,layer),
	158	& layer = 1, nlay_i )
	159	WRITE(numout,*) ' t_i_int : ', ( t_i_int (ji,layer),
	160	& layer = 1, nlay_i+1 )
	161	WRITE(numout,*)
	162	ENDIF ! ln_write
	163
	164	!
	165	!------------------------------------------------------------------------------\|
	166	! 2) Brine salinity at layer mid points and interfaces, brine fraction
	167	!------------------------------------------------------------------------------\|
	168	!
	169	DO layer = 1, nlay_i
	170	e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer)
	171	& - tpw )
	172	END DO
	173
	174	IF ( c_sbr .EQ. 'LIN' ) THEN ! Linear liquidus
	175
	176	DO layer = 1, nlay_i + 1
	177	zTc = t_i_int(ji,layer) - tpw
	178	z_sbr_int(layer) = - zTc / tmut !--- interfacial value
	179	END DO
	180
	181	DO layer = 1, nlay_i
	182	zTc = t_i_b(ji,layer) - tpw
	183	z_sbr_i(layer) = - zTc / tmut !--- mid-point value
	184	END DO
	185
	186	ENDIF
	187
	188	IF ( c_sbr .EQ. 'WEA' ) THEN ! Weast (1971) 3rd order liquidus
	189
	190	DO layer = 1, nlay_i + 1
	191	zTc = t_i_int(ji,layer) - tpw
	192	z_sbr_int(layer) = -17.6zTc -0.389zTc*2. -0.00362zTc**3.
	193	END DO
	194
	195	DO layer = 1, nlay_i
	196	zTc = t_i_b(ji,layer) - tpw
	197	z_sbr_i(layer) = -17.6zTc -0.389zTc*2. -0.00362zTc**3.
	198	END DO
	199
	200	ENDIF
	201
	202	IF ( c_sbr .EQ. 'NTZ' ) THEN ! Notz (2005) 3rd order liquidus
	203
	204	DO layer = 1, nlay_i + 1
	205	zTc = t_i_int(ji,layer) - tpw
	206	z_sbr_int(layer) = -21.4zTc - 0.886zTc*2. - 0.017zTc**3.
	207	END DO
	208
	209	DO layer = 1, nlay_i
	210	zTc = t_i_b(ji,layer) - tpw
	211	z_sbr_i(layer) = -21.4zTc - 0.886zTc*2. - 0.017zTc**3.
	212	END DO
	213
	214	ENDIF
	215
	216	IF ( ln_write ) THEN
	217	WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer),
	218	& layer = 1, nlay_i )
	219	WRITE(numout,*) ' z_sbr_int : ', ( z_sbr_int (layer),
	220	& layer = 1, nlay_i + 1 )
	221	WRITE(numout,*) ' e_i_b : ', ( e_i_b (layer),
	222	& layer = 1, nlay_i )
	223	ENDIF ! ln_write
	224	!
	225	!------------------------------------------------------------------------------\|
	226	! 3) Effective permeability at layer mid-points
	227	!------------------------------------------------------------------------------\|
	228	!
	229	IF ( c_permeff .EQ. 'HAR' ) THEN ! Harmonic mean
	230	WRITE(numout,*)
	231	ENDIF
	232
	233	IF ( c_permeff .EQ. 'MIN' ) THEN ! Minimum permeability
	234
	235	DO layer = 1, nlay_i
	236	ze_i_min = 99999.0
	237	DO layer2 = layer, nlay_i
	238	ze_i_min = MIN( ze_i_min , e_i_b(layer2) )
	239	IF ( c_perm .EQ. 'FRE' ) ! Freitag (1999)
	240	& zperm_eff(layer) = 1.995e-8 * ze_i_min**3.1
	241	IF ( c_perm .EQ. 'RJW' ) ! Rees-Jones and Worster (2014)
	242	& zperm_eff(layer) = 1.0e-8 * ze_i_min**3
	243	END DO
	244	END DO ! layer
	245
	246	ENDIF
	247
	248	!
	249	!------------------------------------------------------------------------------\|
	250	! 4) Rayleigh number at mid-points (see Vancoppenolle et al TCD2013)
	251	!------------------------------------------------------------------------------\|
	252	!
	253	z1 = gpes / ( thdiff_br * visc_br );
	254
	255	DO layer = 1, nlay_i
	256
	257	z2 = beta_ocs * ( z_sbr_i(layer) - oce_sal ) ! Delta rho
	258	z3 = ht_i_b(ji) - z_i_phy(layer)
	259	z4 = zperm_eff(layer);
	260
	261	rayleigh(layer) = MAX( z1 * z2 * z3 * z4, 0.0)
	262
	263	END DO
	264
	265	!
	266	!------------------------------------------------------------------------------\|
	267	! 5) Brine diffusivity
	268	!------------------------------------------------------------------------------\|
	269	!
	270
	271	!-------------------
	272	! Brine Diffusivity
	273	!-------------------
	274	DO layer = 1, nlay_i
	275	zalphara = ( TANH( ra_smooth * ( rayleigh(layer) - ra_c ) )
	276	& + 1 ) / 2.0
	277	diff_br(layer) = ( 1.0 - zalphara ) * d_br_mol +
	278	& zalphara * ( d_br_tur )
	279	IF ( .NOT. ln_grd ) diff_br(layer) = 0.
	280	END DO
	281
	282
	283	IF ( ln_write ) THEN
	284	WRITE(numout,*)
	285	WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer),
	286	& layer = 1, nlay_i )
	287	WRITE(numout,*) ' diff_br : ', ( diff_br(layer),
	288	& layer = 1, nlay_i )
	289	WRITE(numout,*)
	290	ENDIF
	291
	292	!
	293	!------------------------------------------------------------------------------\|
	294	! 6) Brine velocity
	295	!------------------------------------------------------------------------------\|
	296	!
	297	! c_wadv = 'GN' or 'RW'
	298
	299	alpha_GN = 1.56e-3
	300	rho_br_GN = 1020.
	301	Rc_GN = 1.01
[20]	302	!alpha_GN = 1.0e-3
	303	!Rc_GN = 1.0
[10]	304
	305	w_adv_br(:) = 0.0
	306	zRae(:) = 0.
	307	DO layer = 1, nlay_i
	308	zRae(layer) = MAX( rayleigh(layer) - Rc_GN, 0.) ! correction from the orig scheme
	309	END DO
	310
	311	DO layer = 1, nlay_i
	312	w_adv_br(layer) = - alpha_GN / rho_br_GN *
	313	& SUM ( zRae(1:layer)*deltaz_i_phy(1:layer) )
	314	END DO
	315
	316	IF ( ln_write ) THEN
	317	WRITE(numout,*)
	318	WRITE(numout,*) ' w_adv_br : ', ( w_adv_br(layer),
	319	& layer = 1, nlay_i )
	320	ENDIF
	321	!
	322	!------------------------------------------------------------------------------\|
	323	! 7) New salinities
	324	!------------------------------------------------------------------------------\|
	325	!
[20]	326
[10]	327	DO layer = 1, nlay_i
	328	za(layer) = w_adv_br(layer) * ddtb / deltaz_i_phy(layer)
	329	END DO
	330
	331	! first layer
	332	sn_i_b(1) = z_sbr_i(1) * ( e_i_b(1) + za(1) ) +
	333	& z_sbr_i(2) * ( - za(1) )
	334
[20]	335
[10]	336	! inner layers
	337	DO layer = 2, nlay_i - 1
	338	sn_i_b(layer) = z_sbr_i(layer-1) * ( za(layer)/2. ) +
	339	& z_sbr_i(layer) * e_i_b(layer) +
	340	& z_sbr_i(layer+1) * ( - za(layer)/2. )
	341	END DO
	342
	343	! lowermost layer
	344	sn_i_b(nlay_i) = z_sbr_i(nlay_i-1) * ( za(nlay_i)/2. ) +
	345	& z_sbr_i(nlay_i) * ( e_i_b(nlay_i) +
	346	& za(nlay_i)/2. ) - za(nlay_i) * oce_sal
[20]	347
[10]	348
	349	IF ( ln_write ) THEN
	350	WRITE(numout,*)
	351	WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) ,
	352	& layer = 1, nlay_i )
	353	ENDIF
	354
	355	!
	356	!-----------------------------------------------------------------------
	357	! 8) Mass of salt conserved ?
	358	!-----------------------------------------------------------------------
	359	!
	360	! Final mass of salt
	361	CALL ice_sal_column( kideb , kiut , z_ms_i_fin ,
	362	& sn_i_b(1:nlay_i),
	363	& deltaz_i_phy, nlay_i, .FALSE. )
	364
	365	! Bottom flux ( positive upwards for conservation routine )
	366	zfb = - e_i_b(nlay_i) *
	367	& ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) *
	368	& ( z_sbr_i(nlay_i) - oce_sal ) + w_flood * ( z_flood *
	369	& oce_sal + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) )
	370	& - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i)
	371	! & - qsummer * z_sbr_i(nlay_i) / ddtb
	372
	373	fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards
	374	IF ( ln_write ) THEN
	375	WRITE(numout,*) ' fsb : ', fsb
	376	WRITE(numout,*)
	377	ENDIF
	378
	379	! Surface flux of salt
	380	zfsu = 0.0
	381
	382	! conservation check
	383	zerror = 1.0e-15
	384	CALL ice_sal_conserv(kideb,kiut,'ice_sal_adv : ',zerror,
	385	& z_ms_i_ini,z_ms_i_fin,
	386	& zfb , zfsu , ddtb)
	387
	388	ENDIF ! ln_sal
[20]	389
	390
[10]	391	!
	392	!------------------------------------------------------------------------------\|
	393	! End of la sous-routine
	394	WRITE(numout,*)
	395	END SUBROUTINE

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

Context Navigation

source: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_sal_adv.f @ 59

Download in other formats: