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.

limthd_dh.F90 in branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5449

Last change on this file since 5449 was 4924, checked in by mathiot, 10 years ago
UKM02_ice_shelves merged and SETTE tested with revision 4879 of trunk
Property svn:keywords set to `Id`
File size: 36.5 KB

Rev	Line
[825]	1	MODULE limthd_dh
[1572]	2	!!======================================================================
	3	!! * MODULE limthd_dh *
	4	!! LIM-3 : thermodynamic growth and decay of the ice
	5	!!======================================================================
	6	!! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D
	7	!! ! 2005-06 (M. Vancoppenolle) 3D version
[4924]	8	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
[3625]	9	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
	10	!! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes
[1572]	11	!!----------------------------------------------------------------------
[825]	12	#if defined key_lim3
[834]	13	!!----------------------------------------------------------------------
	14	!! 'key_lim3' LIM3 sea-ice model
	15	!!----------------------------------------------------------------------
[3625]	16	!! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice
[825]	17	!!----------------------------------------------------------------------
[3625]	18	USE par_oce ! ocean parameters
	19	USE phycst ! physical constants (OCE directory)
	20	USE sbc_oce ! Surface boundary condition: ocean fields
	21	USE ice ! LIM variables
	22	USE par_ice ! LIM parameters
	23	USE thd_ice ! LIM thermodynamics
	24	USE in_out_manager ! I/O manager
	25	USE lib_mpp ! MPP library
	26	USE wrk_nemo ! work arrays
	27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4924]	28	USE cpl_oasis3, ONLY : lk_cpl
	29
[825]	30	IMPLICIT NONE
	31	PRIVATE
	32
[1572]	33	PUBLIC lim_thd_dh ! called by lim_thd
[825]	34
[4333]	35	REAL(wp) :: epsi20 = 1.e-20 ! constant values
	36	REAL(wp) :: epsi10 = 1.e-10 !
[825]	37
	38	!!----------------------------------------------------------------------
[4161]	39	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]	40	!! $Id$
[2715]	41	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	42	!!----------------------------------------------------------------------
	43	CONTAINS
	44
[4924]	45	SUBROUTINE lim_thd_dh( kideb, kiut )
[921]	46	!!------------------------------------------------------------------
	47	!! * ROUTINE lim_thd_dh *
	48	!!
[1572]	49	!! ** Purpose : determines variations of ice and snow thicknesses.
[921]	50	!!
[1572]	51	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
	52	!! Bottom accretion/ablation arises from flux budget
	53	!! Snow thickness can increase by precipitation and decrease by sublimation
	54	!! If snow load excesses Archmiede limit, snow-ice is formed by
	55	!! the flooding of sea-water in the snow
[921]	56	!!
[1572]	57	!! 1) Compute available flux of heat for surface ablation
	58	!! 2) Compute snow and sea ice enthalpies
	59	!! 3) Surface ablation and sublimation
	60	!! 4) Bottom accretion/ablation
	61	!! 5) Case of Total ablation
	62	!! 6) Snow ice formation
[921]	63	!!
[1572]	64	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
	65	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
	66	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
	67	!! Vancoppenolle et al.,2009, Ocean Modelling
[921]	68	!!------------------------------------------------------------------
[1572]	69	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
	70	!!
	71	INTEGER :: ji , jk ! dummy loop indices
[4161]	72	INTEGER :: ii, ij ! 2D corresponding indices to ji
[1572]	73	INTEGER :: iter
[825]	74
[4924]	75	REAL(wp) :: ztmelts ! local scalar
	76	REAL(wp) :: zdh, zfdum !
[1572]	77	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
	78	REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads
[4924]	79	REAL(wp) :: zs_snic ! snow-ice salinity
[1572]	80	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
	81	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
	82	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
	83	REAL(wp) :: zgrr ! bottom growth rate
[4924]	84	REAL(wp) :: zt_i_new ! bottom formation temperature
	85
	86	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
	87	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
	88	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
	89	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
	90	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
	91	REAL(wp) :: zsstK ! SST in Kelvin
	92
[3294]	93	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness
[4924]	94	REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3)
	95	REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2)
	96	REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2)
	97	REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)
	98	REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
	99	REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
	100	INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting
[3294]	101
[3625]	102	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt
	103	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation
	104	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation
[3294]	105
	106	REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah
[4924]	107	REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness
[3294]	108
[4924]	109	REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2)
	110	REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2)
	111	REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3)
[3294]	112
[4161]	113	! mass and salt flux (clem)
[4924]	114	REAL(wp) :: zdvres, zswitch_sal, zswitch
[4161]	115
[3294]	116	! Heat conservation
[4924]	117	INTEGER :: num_iter_max
	118	REAL(wp) :: zinda, zindq, zindh
	119	REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug
	120
[1572]	121	!!------------------------------------------------------------------
[825]	122
[4924]	123	! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
	124	SELECT CASE( num_sal ) ! varying salinity or not
	125	CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile
	126	CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile
	127	END SELECT
[825]	128
[4924]	129	CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
	130	CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
	131	CALL wrk_alloc( jpij, zintermelt )
	132	CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i )
	133	CALL wrk_alloc( jpij, icount )
[4161]	134
[4924]	135	dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp
	136	dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp
	137
	138	zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp
	139	zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp
[2715]	140
[4924]	141	zh_s (:) = 0._wp
	142	zdh_s_pre(:) = 0._wp
	143	zdh_s_mel(:) = 0._wp
	144	zdh_s_sub(:) = 0._wp
	145	zqh_s (:) = 0._wp
	146	zqh_i (:) = 0._wp
[4161]	147
[4924]	148	zh_i (:,:) = 0._wp
	149	zdeltah (:,:) = 0._wp
	150	zintermelt(:) = 0._wp
	151	icount (:) = 0
	152
	153	! initialize layer thicknesses and enthalpies
	154	h_i_old (:,0:nlay_i+1) = 0._wp
	155	qh_i_old(:,0:nlay_i+1) = 0._wp
	156	DO jk = 1, nlay_i
	157	DO ji = kideb, kiut
	158	h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
	159	qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
	160	ENDDO
	161	ENDDO
[921]	162	!
	163	!------------------------------------------------------------------------------!
[4924]	164	! 1) Calculate available heat for surface and bottom ablation !
[921]	165	!------------------------------------------------------------------------------!
	166	!
[2715]	167	DO ji = kideb, kiut
[4924]	168	zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )
	169	ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt
[825]	170
[4924]	171	zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
	172	zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)
	173
	174	zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) )
	175	zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
	176	END DO
	177
[921]	178	!
	179	!------------------------------------------------------------------------------!
[4924]	180	! If snow temperature is above freezing point, then snow melts
	181	! (should not happen but sometimes it does)
[921]	182	!------------------------------------------------------------------------------!
[4924]	183	DO ji = kideb, kiut
	184	IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting
	185	! Contribution to heat flux to the ocean [W.m-2], < 0
	186	hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
	187	! Contribution to mass flux
	188	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
	189	! updates
	190	ht_s_1d(ji) = 0._wp
	191	q_s_1d (ji,1) = 0._wp
	192	t_s_1d (ji,1) = rtt
	193	END IF
	194	END DO
	195
	196	!------------------------------------------------------------!
	197	! 2) Computing layer thicknesses and enthalpies. !
	198	!------------------------------------------------------------!
[921]	199	!
[4924]	200	DO ji = kideb, kiut
	201	zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
[825]	202	END DO
[2715]	203	!
[825]	204	DO jk = 1, nlay_s
[2715]	205	DO ji = kideb, kiut
[4924]	206	zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)
[825]	207	END DO
	208	END DO
[2715]	209	!
[825]	210	DO jk = 1, nlay_i
[2715]	211	DO ji = kideb, kiut
[4924]	212	zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
	213	zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
[825]	214	END DO
	215	END DO
[921]	216	!
	217	!------------------------------------------------------------------------------\|
	218	! 3) Surface ablation and sublimation \|
	219	!------------------------------------------------------------------------------\|
	220	!
[834]	221	!-------------------------
	222	! 3.1 Snow precips / melt
	223	!-------------------------
[825]	224	! Snow accumulation in one thermodynamic time step
	225	! snowfall is partitionned between leads and ice
	226	! if snow fall was uniform, a fraction (1-at_i) would fall into leads
	227	! but because of the winds, more snow falls on leads than on sea ice
	228	! and a greater fraction (1-at_i)^beta of the total mass of snow
[834]	229	! (beta < 1) falls in leads.
[825]	230	! In reality, beta depends on wind speed,
	231	! and should decrease with increasing wind speed but here, it is
[834]	232	! considered as a constant. an average value is 0.66
[825]	233	! Martin Vancoppenolle, December 2006
	234
	235	DO ji = kideb, kiut
[4924]	236	!-----------
	237	! Snow fall
	238	!-----------
	239	! thickness change
	240	zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji)
[825]	241	zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
[4924]	242	! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
	243	zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )
	244	IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
	245	! heat flux from snow precip (>0, W.m-2)
	246	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
	247	! mass flux, <0
	248	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
	249	! update thickness
	250	ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
[825]	251
[4924]	252	!---------------------
	253	! Melt of falling snow
	254	!---------------------
	255	! thickness change
	256	IF( zdh_s_pre(ji) > 0._wp ) THEN
	257	zindq = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
	258	zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
	259	zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting
	260	! heat used to melt snow (W.m-2, >0)
	261	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
	262	! snow melting only = water into the ocean (then without snow precip), >0
	263	wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice
	264
	265	! updates available heat + thickness
	266	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )
	267	ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) )
	268	zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s )
	269
	270	ENDIF
[825]	271	END DO
	272
[4924]	273	! If heat still available, then melt more snow
	274	zdeltah(:,:) = 0._wp ! important
[825]	275	DO jk = 1, nlay_s
	276	DO ji = kideb, kiut
[4924]	277	! thickness change
	278	zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )
	279	zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) )
	280	zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
	281	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
	282	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
	283	! heat used to melt snow(W.m-2, >0)
	284	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice
	285	! snow melting only = water into the ocean (then without snow precip)
	286	wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[825]	287
[4924]	288	! updates available heat + thickness
	289	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
	290	ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
[825]	291
[1572]	292	END DO
	293	END DO
[825]	294
[4924]	295	!----------------------
	296	! 3.2 Snow sublimation
	297	!----------------------
	298	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
	299	! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)
	300	! clem comment: ice should also sublimate
	301	IF( lk_cpl ) THEN
	302	! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
	303	zdh_s_sub(:) = 0._wp
	304	ELSE
	305	! forced mode: snow thickness change due to sublimation
[921]	306	DO ji = kideb, kiut
[4924]	307	zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
	308	! Heat flux by sublimation [W.m-2], < 0
	309	! sublimate first snow that had fallen, then pre-existing snow
	310	zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + &
	311	& ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) &
	312	& * a_i_1d(ji) * r1_rdtice
	313	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff
	314	! Mass flux by sublimation
	315	wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
	316	! new snow thickness
	317	ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
[1572]	318	END DO
	319	ENDIF
[825]	320
[4924]	321	! --- Update snow diags --- !
[825]	322	DO ji = kideb, kiut
[4924]	323	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
	324	zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
	325	END DO ! ji
[825]	326
[4924]	327	!-------------------------------------------
	328	! 3.3 Update temperature, energy
	329	!-------------------------------------------
	330	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
	331	zq_s(:) = 0._wp
[825]	332	DO jk = 1, nlay_s
	333	DO ji = kideb,kiut
[4924]	334	zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) )
	335	q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) * &
	336	& ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + &
	337	& ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) )
	338	zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk)
[825]	339	END DO
	340	END DO
	341
[4924]	342	!--------------------------
	343	! 3.4 Surface ice ablation
	344	!--------------------------
	345	zdeltah(:,:) = 0._wp ! important
	346	DO jk = 1, nlay_i
	347	DO ji = kideb, kiut
	348	zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]
	349
	350	ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K]
	351
	352	zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0]
	353
	354	zdE = zEi - zEw ! Specific enthalpy difference < 0
	355
	356	zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
	357
	358	zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0]
	359
	360	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]
	361
	362	zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
	363
	364	dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt
	365
	366	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
	367
	368	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
	369
	370	! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
	371	sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
	372
	373	! Contribution to heat flux [W.m-2], < 0
	374	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
	375
	376	! Total heat flux used in this process [W.m-2], > 0
	377	hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
	378
	379	! Contribution to mass flux
	380	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
	381
	382	! record which layers have disappeared (for bottom melting)
	383	! => icount=0 : no layer has vanished
	384	! => icount=5 : 5 layers have vanished
	385	zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) )
	386	icount(ji) = icount(ji) + zindh
	387	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
	388
	389	! update heat content (J.m-2) and layer thickness
	390	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
	391	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]	392	END DO
[921]	393	END DO
[4924]	394	! update ice thickness
	395	DO ji = kideb, kiut
	396	ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
	397	END DO
[825]	398
[921]	399	!
	400	!------------------------------------------------------------------------------!
	401	! 4) Basal growth / melt !
	402	!------------------------------------------------------------------------------!
	403	!
[4924]	404	!------------------
	405	! 4.1 Basal growth
	406	!------------------
	407	! Basal growth is driven by heat imbalance at the ice-ocean interface,
	408	! between the inner conductive flux (fc_bo_i), from the open water heat flux
	409	! (fhld) and the turbulent ocean flux (fhtur).
	410	! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
[825]	411
[4924]	412	! If salinity varies in time, an iterative procedure is required, because
	413	! the involved quantities are inter-dependent.
	414	! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
	415	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
	416	! -> need for an iterative procedure, which converges quickly
	417
	418	IF ( num_sal == 2 ) THEN
	419	num_iter_max = 5
	420	ELSE
	421	num_iter_max = 1
[1572]	422	ENDIF
[825]	423
[4924]	424	!clem debug. Just to be sure that enthalpy at nlay_i+1 is null
	425	DO ji = kideb, kiut
	426	q_i_1d(ji,nlay_i+1) = 0._wp
	427	END DO
[825]	428
[4924]	429	! Iterative procedure
	430	DO iter = 1, num_iter_max
	431	DO ji = kideb, kiut
	432	IF( zf_tt(ji) < 0._wp ) THEN
[825]	433
[4924]	434	! New bottom ice salinity (Cox & Weeks, JGR88 )
	435	!--- zswi1 if dh/dt < 2.0e-8
	436	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
	437	!--- zswi2 if dh/dt > 3.6e-7
	438	zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
	439	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
	440	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
	441	zswi1 = 1. - zswi2 * zswi12
	442	zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
	443	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
[825]	444
[4924]	445	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
[825]	446
[4924]	447	s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity
	448	+ ( 1. - zswitch_sal ) * sm_i_1d(ji)
	449	! New ice growth
	450	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
[825]	451
[4924]	452	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
	453
	454	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
	455	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
	456	& + rcp * ( ztmelts-rtt )
	457
	458	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
	459
	460	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
	461
	462	dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
	463
	464	q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)
	465
	466	ENDIF ! fc_bo_i
	467	END DO ! ji
	468	END DO ! iter
	469
	470	! Contribution to Energy and Salt Fluxes
[825]	471	DO ji = kideb, kiut
[4924]	472	IF( zf_tt(ji) < 0._wp ) THEN
	473	! New ice growth
	474
	475	zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0)
	476
	477	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
	478
	479	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
	480
	481	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
	482	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
	483	& + rcp * ( ztmelts-rtt )
	484
	485	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
	486
	487	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
	488
	489	! Contribution to heat flux to the ocean [W.m-2], >0
	490	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
	491
	492	! Total heat flux used in this process [W.m-2], <0
	493	hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
	494
	495	! Contribution to salt flux, <0
	496	sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice
	497
	498	! Contribution to mass flux, <0
	499	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
	500
	501	! update heat content (J.m-2) and layer thickness
	502	qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1)
	503	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
[825]	504	ENDIF
	505	END DO
	506
[4924]	507	!----------------
	508	! 4.2 Basal melt
	509	!----------------
	510	zdeltah(:,:) = 0._wp ! important
[825]	511	DO jk = nlay_i, 1, -1
	512	DO ji = kideb, kiut
[4924]	513	IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting
[825]	514
[4924]	515	ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K)
[825]	516
[4924]	517	IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
	518	zintermelt(ji) = 1._wp
[825]	519
[4924]	520	zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
[825]	521
[4924]	522	!!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)
[4161]	523
[4924]	524	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
	525	! set up at 0 since no energy is needed to melt water...(it is already melted)
[4161]	526
[4924]	527	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	528	! this should normally not happen, but sometimes, heat diffusion leads to this
[825]	529
[4924]	530	dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk)
[825]	531
[4924]	532	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
[825]	533
[4924]	534	! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
	535	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
[825]	536
[4924]	537	! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
	538	sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
	539
	540	! Contribution to mass flux
	541	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
[825]	542
[4924]	543	! update heat content (J.m-2) and layer thickness
	544	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
	545	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
[825]	546
[4924]	547	ELSE !!! Basal melting
[825]	548
[4924]	549	zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
[825]	550
[4924]	551	zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
[825]	552
[4924]	553	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
	554
	555	zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)
	556
	557	zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change
	558
	559	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
	560
	561	zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
	562
	563	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt
	564
	565	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
	566
	567	zQm = zfmdt * zEw ! Heat exchanged with ocean
	568
	569	! Contribution to heat flux to the ocean [W.m-2], <0
	570	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
	571
	572	! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
	573	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
	574
	575	! Total heat flux used in this process [W.m-2], >0
	576	hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
	577
	578	! Contribution to mass flux
	579	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
	580
	581	! update heat content (J.m-2) and layer thickness
	582	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
	583	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	584	ENDIF
	585
	586	ENDIF
	587	END DO ! ji
	588	END DO ! jk
	589
	590	!------------------------------------------------------------------------------!
	591	! Excessive ablation in a 1-category model
	592	! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
	593	!------------------------------------------------------------------------------!
	594	! ??? keep ???
	595	! clem bug: I think this should be included above, so we would not have to
	596	! track heat/salt/mass fluxes backwards
	597	! IF( jpl == 1 ) THEN
	598	! DO ji = kideb, kiut
	599	! IF( zf_tt(ji) >= 0._wp ) THEN
	600	! zdh = MAX( hmelt , dh_i_bott(ji) )
	601	! zdvres = zdh - dh_i_bott(ji) ! >=0
	602	! dh_i_bott(ji) = zdh
	603	!
	604	! ! excessive energy is sent to lateral ablation
	605	! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) )
	606	! zq_1cat(ji) = zinda * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0
	607	!
	608	! ! correct salt and mass fluxes
	609	! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
	610	! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice
	611	! ENDIF
	612	! END DO
	613	! ENDIF
	614
[825]	615	!-------------------------------------------
[4924]	616	! Update temperature, energy
[825]	617	!-------------------------------------------
	618	DO ji = kideb, kiut
[4924]	619	ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
	620	END DO
[825]	621
[4924]	622	!-------------------------------------------
	623	! 5. What to do with remaining energy
	624	!-------------------------------------------
	625	! If heat still available for melting and snow remains, then melt more snow
	626	!-------------------------------------------
	627	zdeltah(:,:) = 0._wp ! important
	628	DO ji = kideb, kiut
	629	zq_rema(ji) = zq_su(ji) + zq_bo(ji)
	630	! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow
	631	! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
	632	! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
	633	! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
	634	! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)
	635	! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1)
	636	! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1)
	637	!
	638	! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2)
	639	! ! heat used to melt snow
	640	! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
	641	! ! Contribution to mass flux
	642	! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
	643	!
	644	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	645	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
	646	hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
[825]	647
[4924]	648	IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
	649	END DO
	650
[921]	651	!
	652	!------------------------------------------------------------------------------\|
	653	! 6) Snow-Ice formation \|
	654	!------------------------------------------------------------------------------\|
[1572]	655	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
	656	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
[825]	657	DO ji = kideb, kiut
[1572]	658	!
[4924]	659	dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) )
[825]	660
[4924]	661	ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji)
	662	ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji)
[825]	663
[4924]	664	! Salinity of snow ice
	665	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	666	zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
[825]	667
	668	! entrapment during snow ice formation
[4924]	669	! new salinity difference stored (to be used in limthd_ent.F90)
[4161]	670	IF ( num_sal == 2 ) THEN
[4924]	671	zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) )
[4161]	672	! salinity dif due to snow-ice formation
[4924]	673	dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch
[4161]	674	! salinity dif due to bottom growth
[4924]	675	IF ( zf_tt(ji) < 0._wp ) THEN
	676	dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch
[4161]	677	ENDIF
	678	ENDIF
[825]	679
[4924]	680	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
	681	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
	682	zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0
	683	zsstK = sst_m(ii,ij) + rt0
	684	zEw = rcp * ( zsstK - rt0 )
	685	zQm = zfmdt * zEw
	686
	687	! Contribution to heat flux
	688	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
[825]	689
[4924]	690	! Contribution to salt flux
	691	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice
	692
	693	! Contribution to mass flux
	694	! All snow is thrown in the ocean, and seawater is taken to replace the volume
	695	wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
	696	wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
[825]	697
[4924]	698	! update heat content (J.m-2) and layer thickness
	699	qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
	700	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
	701
	702	! Total ablation (to debug)
	703	IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp
[4161]	704
[825]	705	END DO !ji
[4924]	706
[2715]	707	!
[4924]	708	!-------------------------------------------
	709	! Update temperature, energy
	710	!-------------------------------------------
	711	!clem bug: we should take snow into account here
	712	DO ji = kideb, kiut
	713	zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )
	714	t_su_1d(ji) = zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt
	715	END DO ! ji
	716
	717	DO jk = 1, nlay_s
	718	DO ji = kideb,kiut
	719	! mask enthalpy
	720	zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) )
	721	q_s_1d(ji,jk) = ( 1.0 - zinda ) * q_s_1d(ji,jk)
	722	! recalculate t_s_1d from q_s_1d
	723	t_s_1d(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
	724	END DO
	725	END DO
	726
	727	CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
	728	CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
	729	CALL wrk_dealloc( jpij, zintermelt )
	730	CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i )
	731	CALL wrk_dealloc( jpij, icount )
[2715]	732	!
[4161]	733	!
[921]	734	END SUBROUTINE lim_thd_dh
[1572]	735
[825]	736	#else
[1572]	737	!!----------------------------------------------------------------------
	738	!! Default option NO LIM3 sea-ice model
	739	!!----------------------------------------------------------------------
[825]	740	CONTAINS
	741	SUBROUTINE lim_thd_dh ! Empty routine
	742	END SUBROUTINE lim_thd_dh
	743	#endif
[1572]	744
	745	!!======================================================================
[921]	746	END MODULE limthd_dh

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

Download in other formats: