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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 6416

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

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

Download in other formats: