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

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 8537

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

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

Download in other formats: