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.

icethd_dh.F90 in NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/ICE/icethd_dh.F90 @ 9900

Last change on this file since 9900 was 9900, checked in by davestorkey, 6 years ago
UKMO dev_r9888_GO6_mixing branch: clear SVN keywords.
File size: 35.5 KB

Rev	Line
[8586]	1	MODULE icethd_dh
	2	!!======================================================================
	3	!! * MODULE icethd_dh *
	4	!! seaice : thermodynamic growth and melt
	5	!!======================================================================
[9604]	6	!! History : ! 2003-05 (M. Vancoppenolle) Original code in 1D
	7	!! ! 2005-06 (M. Vancoppenolle) 3D version
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8586]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8586]	13	!!----------------------------------------------------------------------
	14	!! ice_thd_dh : vertical sea-ice growth and melt
	15	!! ice_thd_snwblow : distribute snow fall between ice and ocean
	16	!!----------------------------------------------------------------------
	17	USE dom_oce ! ocean space and time domain
	18	USE phycst ! physical constants
	19	USE ice ! sea-ice: variables
	20	USE ice1D ! sea-ice: thermodynamics variables
	21	USE icethd_sal ! sea-ice: salinity profiles
	22	!
	23	USE in_out_manager ! I/O manager
	24	USE lib_mpp ! MPP library
	25	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
	30	PUBLIC ice_thd_dh ! called by ice_thd
[9274]	31	PUBLIC ice_thd_snwblow ! called in sbcblk/sbccpl and here
[8586]	32
	33	INTERFACE ice_thd_snwblow
	34	MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
	35	END INTERFACE
	36
	37	!!----------------------------------------------------------------------
[9598]	38	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[9900]	39	!! $Id$
[9598]	40	!! Software governed by the CeCILL licence (./LICENSE)
[8586]	41	!!----------------------------------------------------------------------
	42	CONTAINS
	43
	44	SUBROUTINE ice_thd_dh
	45	!!------------------------------------------------------------------
	46	!! * ROUTINE ice_thd_dh *
	47	!!
[9274]	48	!! ** Purpose : compute ice and snow thickness changes due to growth/melting
[8586]	49	!!
	50	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
[9274]	51	!! Bottom accretion/ablation arises from flux budget
	52	!! Snow thickness can increase by precipitation and decrease by sublimation
	53	!! If snow load excesses Archmiede limit, snow-ice is formed by
	54	!! the flooding of sea-water in the snow
[8586]	55	!!
[9274]	56	!! - Compute available flux of heat for surface ablation
	57	!! - Compute snow and sea ice enthalpies
	58	!! - Surface ablation and sublimation
	59	!! - Bottom accretion/ablation
	60	!! - Snow ice formation
[8586]	61	!!
	62	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
	63	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
	64	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
	65	!! Vancoppenolle et al.,2009, Ocean Modelling
	66	!!------------------------------------------------------------------
	67	INTEGER :: ji, jk ! dummy loop indices
	68	INTEGER :: iter ! local integer
	69
	70	REAL(wp) :: ztmelts ! local scalar
	71	REAL(wp) :: zdum
	72	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
	73	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
	74	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
	75	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
	76	REAL(wp) :: zgrr ! bottom growth rate
	77	REAL(wp) :: zt_i_new ! bottom formation temperature
	78	REAL(wp) :: z1_rho ! 1/(rhosn+rau0-rhoic)
	79
	80	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
	81	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
	82	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
	83	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
	84	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
	85
	86	REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)
	87	REAL(wp), DIMENSION(jpij) :: zq_su ! heat for surface ablation (J.m-2)
	88	REAL(wp), DIMENSION(jpij) :: zq_bo ! heat for bottom ablation (J.m-2)
	89	REAL(wp), DIMENSION(jpij) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
	90	REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
	91	REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2)
	92
	93	REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt
	94	REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation
	95	REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation
	96
[9271]	97	REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness
	98	REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness
[8586]	99	REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah
	100	INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting
	101
	102	REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing
	103
	104	REAL(wp) :: zswitch_sal
	105
	106	INTEGER :: num_iter_max ! Heat conservation
	107	!!------------------------------------------------------------------
	108
	109	! Discriminate between time varying salinity and constant
	110	SELECT CASE( nn_icesal ) ! varying salinity or not
	111	CASE( 1 , 3 ) ; zswitch_sal = 0._wp ! prescribed salinity profile
	112	CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
	113	END SELECT
	114
	115	! initialize layer thicknesses and enthalpies
	116	h_i_old (1:npti,0:nlay_i+1) = 0._wp
	117	eh_i_old(1:npti,0:nlay_i+1) = 0._wp
	118	DO jk = 1, nlay_i
	119	DO ji = 1, npti
	120	h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
	121	eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
	122	END DO
	123	END DO
	124	!
[9604]	125	! ! ============================================== !
	126	! ! Available heat for surface and bottom ablation !
	127	! ! ============================================== !
[9274]	128	SELECT CASE( nice_jules )
[8813]	129	!
[9274]	130	CASE( np_jules_ACTIVE )
[8813]	131	!
	132	DO ji = 1, npti
[9274]	133	zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
[8813]	134	END DO
	135	!
[9274]	136	CASE( np_jules_EMULE )
[8813]	137	!
	138	DO ji = 1, npti
[9274]	139	zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji)
	140	qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
	141	zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
[8813]	142	END DO
	143	!
[9274]	144	CASE( np_jules_OFF )
[8813]	145	!
	146	DO ji = 1, npti
[9274]	147	zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji)
[9767]	148	qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
	149	zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
[8813]	150	END DO
	151	!
	152	END SELECT
	153	!
[8586]	154	DO ji = 1, npti
[9274]	155	zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)
	156	zq_bo(ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
[8586]	157	END DO
	158
[9274]	159	! Ice and snow layer thicknesses
	160	!-------------------------------
[9271]	161	DO jk = 1, nlay_i
	162	DO ji = 1, npti
	163	zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
	164	END DO
	165	END DO
	166
	167	DO jk = 1, nlay_s
	168	DO ji = 1, npti
	169	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
	170	END DO
	171	END DO
	172
[9274]	173	! ! ============ !
	174	! ! Snow !
	175	! ! ============ !
	176	!
	177	! Internal melting
	178	! ----------------
	179	! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
[9271]	180	DO jk = 1, nlay_s
[8586]	181	DO ji = 1, npti
[9274]	182	IF( t_s_1d(ji,jk) > rt0 ) THEN
	183	hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0
	184	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux
[9271]	185	! updates
[9274]	186	dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)
	187	h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)
	188	zh_s (ji,jk) = 0._wp
	189	e_s_1d (ji,jk) = 0._wp
	190	t_s_1d (ji,jk) = rt0
[9271]	191	END IF
[8586]	192	END DO
[9274]	193	END DO
[8586]	194
[9274]	195	! Snow precipitation
	196	!-------------------
	197	CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
[8586]	198
	199	zdeltah(1:npti,:) = 0._wp
	200	DO ji = 1, npti
[9274]	201	IF( sprecip_1d(ji) > 0._wp ) THEN
	202	!
	203	! --- precipitation ---
	204	zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) ! thickness change
	205	zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3)
	206	!
	207	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat flux from snow precip (>0, W.m-2)
	208	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! mass flux, <0
	209
	210	! --- melt of falling snow ---
	211	rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
	212	zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change
	213	zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
	214	hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat used to melt snow (W.m-2, >0)
	215	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0
	216
	217	! updates available heat + precipitations after melting
	218	dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
	219	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )
	220	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
	221
	222	! update thickness
	223	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
	224	!
	225	ELSE
	226	!
	227	zdh_s_pre(ji) = 0._wp
	228	zqprec (ji) = 0._wp
	229	!
	230	ENDIF
[8586]	231	END DO
	232
[9271]	233	! recalculate snow layers
	234	DO jk = 1, nlay_s
	235	DO ji = 1, npti
	236	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
	237	END DO
	238	END DO
	239
[9274]	240	! Snow melting
	241	! ------------
[8586]	242	! If heat still available (zq_su > 0), then melt more snow
	243	zdeltah(1:npti,:) = 0._wp
	244	zdh_s_mel(1:npti) = 0._wp
	245	DO jk = 1, nlay_s
	246	DO ji = 1, npti
[9274]	247	IF( zh_s(ji,jk) > 0._wp .AND. zq_su(ji) > 0._wp ) THEN
	248	!
	249	rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
	250	zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change
	251	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting
	252	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
	253
	254	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice ! heat used to melt snow(W.m-2, >0)
	255	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip)
	256
	257	! updates available heat + thickness
	258	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk)
	259	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
	260	h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
	261	zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
	262	!
	263	ENDIF
[8586]	264	END DO
	265	END DO
	266
[9274]	267	! Snow sublimation
	268	!-----------------
[8586]	269	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
[8885]	270	! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
[8586]	271	zdeltah(1:npti,:) = 0._wp
	272	DO ji = 1, npti
[9274]	273	IF( evap_ice_1d(ji) > 0._wp ) THEN
	274	!
	275	zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
	276	zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn ! remaining evap in kg.m-2 (used for ice melting later on)
	277	zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
	278
	279	hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
	280	& ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) &
	281	& * a_i_1d(ji) * r1_rdtice
	282	wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! Mass flux by sublimation
	283
	284	! new snow thickness
	285	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
	286	! update precipitations after sublimation and correct sublimation
	287	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
	288	zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
	289	!
	290	ELSE
	291	!
	292	zdh_s_sub (ji) = 0._wp
	293	zevap_rema(ji) = 0._wp
	294	!
	295	ENDIF
[8586]	296	END DO
	297
	298	! --- Update snow diags --- !
	299	DO ji = 1, npti
	300	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
	301	END DO
	302
[9274]	303	! Update temperature, energy
	304	!---------------------------
[8586]	305	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
	306	DO jk = 1, nlay_s
	307	DO ji = 1,npti
	308	rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
[9274]	309	e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * &
	310	& ( ( zdh_s_pre(ji) ) * zqprec(ji) + &
	311	& ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
[8586]	312	END DO
	313	END DO
[9274]	314
	315	! ! ============ !
	316	! ! Ice !
	317	! ! ============ !
[8586]	318
[9274]	319	! Surface ice melting
	320	!--------------------
[8586]	321	zdeltah(1:npti,:) = 0._wp ! important
	322	DO jk = 1, nlay_i
	323	DO ji = 1, npti
[9274]	324	ztmelts = - tmut * sz_i_1d(ji,jk) ! Melting point of layer k [C]
[8586]	325
[9274]	326	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	327
	328	zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]
[9274]	329	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
[8586]	330	! set up at 0 since no energy is needed to melt water...(it is already melted)
	331	zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	332	! this should normally not happen, but sometimes, heat diffusion leads to this
	333	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
	334
[9750]	335	dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting
[8586]	336
	337	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
	338
[9274]	339	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	340	! ice enthalpy zEi is "sent" to the ocean
	341	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
	342	! using s_i_1d and not sz_i_1d(jk) is ok
	343	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	344
[9274]	345	ELSE !-- Surface melting
[8586]	346
	347	zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]
	348	zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0]
	349	zdE = zEi - zEw ! Specific enthalpy difference < 0
	350
	351	zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
	352
	353	zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0]
	354
	355	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]
	356
	357	zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
	358
[9750]	359	dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt
[8586]	360
	361	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
	362
	363	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
	364
[9274]	365	sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux >0
	366	! using s_i_1d and not sz_i_1d(jk) is ok)
	367	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux [W.m-2], < 0
	368	hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], > 0
	369	!
	370	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	371
	372	END IF
[9274]	373
	374	! Ice sublimation
	375	! ---------------
	376	zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
	377	zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
	378	dh_i_sub(ji) = dh_i_sub(ji) + zdum
	379
	380	sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0
	381	! clem: flux is sent to the ocean for simplicity
	382	! but salt should remain in the ice except
	383	! if all ice is melted. => must be corrected
	384	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0
[8586]	385
[9274]	386	wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice ! Mass flux > 0
	387
[8586]	388	! update remaining mass flux
	389	zevap_rema(ji) = zevap_rema(ji) + zdum * rhoic
	390
	391	! record which layers have disappeared (for bottom melting)
	392	! => icount=0 : no layer has vanished
	393	! => icount=5 : 5 layers have vanished
	394	rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )
	395	icount(ji,jk) = NINT( rswitch )
	396	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
	397
	398	! update heat content (J.m-2) and layer thickness
	399	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	400	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	401	END DO
	402	END DO
[9274]	403
[8586]	404	! update ice thickness
	405	DO ji = 1, npti
[9750]	406	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
[8586]	407	END DO
	408
	409	! remaining "potential" evap is sent to ocean
	410	DO ji = 1, npti
	411	wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice ! <=0 (net evap for the ocean in kg.m-2.s-1)
	412	END DO
	413
	414
[9274]	415	! Ice Basal growth
[8586]	416	!------------------
	417	! Basal growth is driven by heat imbalance at the ice-ocean interface,
	418	! between the inner conductive flux (fc_bo_i), from the open water heat flux
	419	! (fhld) and the turbulent ocean flux (fhtur).
	420	! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
	421
	422	! If salinity varies in time, an iterative procedure is required, because
	423	! the involved quantities are inter-dependent.
[9750]	424	! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
	425	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
[8586]	426	! -> need for an iterative procedure, which converges quickly
	427
	428	num_iter_max = 1
	429	IF( nn_icesal == 2 ) num_iter_max = 5 ! salinity varying in time
	430
	431	DO ji = 1, npti
	432	IF( zf_tt(ji) < 0._wp ) THEN
[9274]	433	DO iter = 1, num_iter_max ! iterations
[8586]	434
	435	! New bottom ice salinity (Cox & Weeks, JGR88 )
	436	!--- zswi1 if dh/dt < 2.0e-8
	437	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
	438	!--- zswi2 if dh/dt > 3.6e-7
[9750]	439	zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) )
[9274]	440	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
	441	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
	442	zswi1 = 1. - zswi2 * zswi12
	443	zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
	444	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
[8586]	445
[9274]	446	s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity
[8586]	447
[9274]	448	ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C)
	449
	450	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
[8586]	451
[9274]	452	zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)
	453	& - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts
[8586]	454
[9274]	455	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
[8586]	456
[9274]	457	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	458
[9750]	459	dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
[8586]	460
	461	END DO
	462	! Contribution to Energy and Salt Fluxes
[9750]	463	zfmdt = - rhoic * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)
[8586]	464
[9274]	465	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], >0
	466	hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], <0
[8586]	467
[9750]	468	sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice ! Salt flux, <0
[8586]	469
[9750]	470	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice ! Mass flux, <0
[8586]	471
	472	! update heat content (J.m-2) and layer thickness
[9750]	473	eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoic)
	474	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
[8586]	475
	476	ENDIF
	477
	478	END DO
	479
[9274]	480	! Ice Basal melt
	481	!---------------
[8586]	482	zdeltah(1:npti,:) = 0._wp ! important
	483	DO jk = nlay_i, 1, -1
	484	DO ji = 1, npti
	485	IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting
	486
	487	ztmelts = - tmut * sz_i_1d(ji,jk) ! Melting point of layer jk (C)
	488
[9274]	489	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
[8586]	490
	491	zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
	492	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
[9274]	493	! set up at 0 since no energy is needed to melt water...(it is already melted)
[8586]	494	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
	495	! this should normally not happen, but sometimes, heat diffusion leads to this
	496
[9750]	497	dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk)
[8586]	498
	499	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
	500
[9274]	501	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	502	! ice enthalpy zEi is "sent" to the ocean
	503	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
	504	! using s_i_1d and not sz_i_1d(jk) is ok
	505	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	506
	507	! update heat content (J.m-2) and layer thickness
	508	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	509	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	510
[9274]	511	ELSE !-- Basal melting
[8586]	512
[9274]	513	zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
	514	zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)
	515	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
[8586]	516
[9274]	517	zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)
[8586]	518
[9274]	519	zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change
[8586]	520
[9274]	521	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
[8586]	522
[9274]	523	zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
[8586]	524
[9750]	525	dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt
[8586]	526
[9274]	527	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
[8586]	528
[9274]	529	zQm = zfmdt * zEw ! Heat exchanged with ocean
[8586]	530
[9274]	531	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], <0
	532	hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat used in this process [W.m-2], >0
[8586]	533
[9274]	534	sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux
	535	! using s_i_1d and not sz_i_1d(jk) is ok
[8586]	536
[9274]	537	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux
[8586]	538
	539	! update heat content (J.m-2) and layer thickness
	540	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
	541	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
	542	ENDIF
	543
	544	ENDIF
	545	END DO
	546	END DO
	547
	548	! Update temperature, energy
[9274]	549	! --------------------------
[8586]	550	DO ji = 1, npti
[9750]	551	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
[8586]	552	END DO
	553
[9274]	554	! If heat still available then melt more snow
[8586]	555	!-------------------------------------------
	556	zdeltah(1:npti,:) = 0._wp ! important
	557	DO ji = 1, npti
[9274]	558	zq_rema (ji) = zq_su(ji) + zq_bo(ji)
	559	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow
	560	rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
	561	zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
	562	zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
	563	dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1)
	564	h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1)
[8586]	565
[9274]	566	zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2)
	567	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! Heat used to melt snow, W.m-2 (>0)
	568	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! Mass flux
[8637]	569	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1)
[8586]	570	!
	571	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
[9274]	572	hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
[8586]	573
	574	IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
	575	END DO
	576
	577	!
[9274]	578	! Snow-Ice formation
	579	! ------------------
[8586]	580	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
	581	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
	582	z1_rho = 1._wp / ( rhosn+rau0-rhoic )
	583	DO ji = 1, npti
	584	!
	585	dh_snowice(ji) = MAX( 0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho )
	586
	587	h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji)
	588	h_s_1d(ji) = h_s_1d(ji) - dh_snowice(ji)
	589
	590	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
	591	zfmdt = ( rhosn - rhoic ) * dh_snowice(ji) ! <0
	592	zEw = rcp * sst_1d(ji)
	593	zQm = zfmdt * zEw
	594
[9274]	595	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux
[8586]	596
[9274]	597	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux
[8586]	598
	599	! Case constant salinity in time: virtual salt flux to keep salinity constant
[8637]	600	IF( nn_icesal /= 2 ) THEN
[8586]	601	sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_rdtice & ! put back sss_m into the ocean
	602	& - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice ! and get rn_icesal from the ocean
	603	ENDIF
	604
[9274]	605	! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
[8586]	606	wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
	607	wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
	608
	609	! update heat content (J.m-2) and layer thickness
	610	eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
	611	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
	612
	613	END DO
	614
	615	!
	616	! Update temperature, energy
[9274]	617	! --------------------------
[8586]	618	DO ji = 1, npti
[9274]	619	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
	620	t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
[8586]	621	END DO
	622
	623	DO jk = 1, nlay_s
	624	DO ji = 1,npti
	625	! mask enthalpy
	626	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - h_s_1d(ji) ) )
	627	e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
	628	! recalculate t_s_1d from e_s_1d
	629	t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic )
	630	END DO
	631	END DO
	632
	633	! --- ensure that a_i = 0 where h_i = 0 ---
	634	WHERE( h_i_1d(1:npti) == 0._wp ) a_i_1d(1:npti) = 0._wp
	635	!
	636	END SUBROUTINE ice_thd_dh
	637
	638
	639	!!--------------------------------------------------------------------------
	640	!! INTERFACE ice_thd_snwblow
[9274]	641	!!
[8586]	642	!! ** Purpose : Compute distribution of precip over the ice
[9274]	643	!!
	644	!! Snow accumulation in one thermodynamic time step
	645	!! snowfall is partitionned between leads and ice.
	646	!! If snow fall was uniform, a fraction (1-at_i) would fall into leads
	647	!! but because of the winds, more snow falls on leads than on sea ice
	648	!! and a greater fraction (1-at_i)^beta of the total mass of snow
	649	!! (beta < 1) falls in leads.
	650	!! In reality, beta depends on wind speed,
	651	!! and should decrease with increasing wind speed but here, it is
	652	!! considered as a constant. an average value is 0.66
[8586]	653	!!--------------------------------------------------------------------------
	654	!!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
	655	SUBROUTINE ice_thd_snwblow_2d( pin, pout )
	656	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b )
	657	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
	658	pout = ( 1._wp - ( pin )**rn_blow_s )
	659	END SUBROUTINE ice_thd_snwblow_2d
	660
	661	SUBROUTINE ice_thd_snwblow_1d( pin, pout )
	662	REAL(wp), DIMENSION(:), INTENT(in ) :: pin
	663	REAL(wp), DIMENSION(:), INTENT(inout) :: pout
	664	pout = ( 1._wp - ( pin )**rn_blow_s )
	665	END SUBROUTINE ice_thd_snwblow_1d
	666
	667	#else
	668	!!----------------------------------------------------------------------
[9570]	669	!! Default option NO SI3 sea-ice model
[8586]	670	!!----------------------------------------------------------------------
	671	#endif
	672
	673	!!======================================================================
	674	END MODULE icethd_dh

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

Download in other formats: