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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 5078

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

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

Download in other formats: