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/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/ICE/icethd_dh.F90 @ 14933

Last change on this file since 14933 was 14495, checked in by mcastril, 3 years ago
Further changes to make the mixed precision branch compliant with cpp -traditional and GNU compiler
Property svn:keywords set to `Id`
File size: 35.9 KB

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

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

Download in other formats: