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.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90 @ 8522

Last change on this file since 8522 was 8522, checked in by clem, 7 years ago
changes in style - part6 - ice diffusion (hope the scheme still converges)
File size: 36.0 KB

Rev	Line
[8422]	1	MODULE icethd
	2	!!======================================================================
	3	!! * MODULE icethd *
	4	!! LIM-3 : ice thermodynamic
	5	!!======================================================================
	6	!! History : LIM ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
	7	!! 2.0 ! 2002-07 (C. Ethe, G. Madec) LIM-2 (F90 rewriting)
	8	!! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations
	9	!! - ! 2007-04 (M. Vancoppenolle) add ice_thd_glohec, ice_thd_con_dh and ice_thd_con_dif
	10	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
	11	!! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
	12	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	13	!! - ! 2012-05 (C. Rousset) add penetration solar flux
	14	!!----------------------------------------------------------------------
	15	#if defined key_lim3
	16	!!----------------------------------------------------------------------
[8486]	17	!! 'key_lim3' LIM3 sea-ice model
[8422]	18	!!----------------------------------------------------------------------
	19	!! ice_thd : thermodynamic of sea ice
	20	!! ice_thd_init : initialisation of sea-ice thermodynamic
	21	!!----------------------------------------------------------------------
	22	USE phycst ! physical constants
	23	USE dom_oce ! ocean space and time domain variables
	24	USE ice ! sea-ice variables
[8486]	25	!!gm list trop longue ==>>> why not passage en argument d'appel ?
[8422]	26	USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
	27	USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
[8512]	28	& fr1_i0, fr2_i0
[8422]	29	USE ice1D ! thermodynamic sea-ice variables
	30	USE icethd_dif ! vertical diffusion
	31	USE icethd_dh ! ice-snow growth and melt
	32	USE icethd_da ! lateral melting
	33	USE icethd_sal ! ice salinity
	34	USE icethd_ent ! ice enthalpy redistribution
	35	USE icethd_lac ! lateral accretion
	36	USE iceitd ! remapping thickness distribution
	37	USE icetab ! 1D <==> 2D transformation
[8424]	38	USE icevar !
[8422]	39	USE icectl ! control print
	40	!
	41	USE in_out_manager ! I/O manager
	42	USE lbclnk ! lateral boundary condition - MPP links
	43	USE lib_mpp ! MPP library
	44	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
	45	USE timing ! Timing
	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
	50	PUBLIC ice_thd ! called by limstp module
	51	PUBLIC ice_thd_init ! called by ice_init
	52
	53	!! * Substitutions
	54	# include "vectopt_loop_substitute.h90"
	55	!!----------------------------------------------------------------------
[8486]	56	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8422]	57	!! $Id: icethd.F90 8420 2017-08-08 12:18:46Z clem $
	58	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	59	!!----------------------------------------------------------------------
	60	CONTAINS
	61
	62	SUBROUTINE ice_thd( kt )
	63	!!-------------------------------------------------------------------
	64	!! * ROUTINE ice_thd *
	65	!!
	66	!! ** Purpose : This routine manages ice thermodynamics
	67	!!
	68	!! ** Action : - Initialisation of some variables
	69	!! - Some preliminary computation (oceanic heat flux
	70	!! at the ice base, snow acc.,heat budget of the leads)
	71	!! - selection of the icy points and put them in an array
	72	!! - call ice_thd_dif for vertical heat diffusion
	73	!! - call ice_thd_dh for vertical ice growth and melt
	74	!! - call ice_thd_ent for enthalpy remapping
	75	!! - call ice_thd_sal for ice desalination
	76	!! - call ice_thd_temp to retrieve temperature from ice enthalpy
	77	!! - back to the geographic grid
	78	!!---------------------------------------------------------------------
	79	INTEGER, INTENT(in) :: kt ! number of iteration
	80	!
	81	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	82	REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
	83	REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
	84	REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
	85	REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
	86	!
	87	!!-------------------------------------------------------------------
[8517]	88	! controls
	89	IF( nn_timing == 1 ) CALL timing_start('icethd') ! timing
	90	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[8422]	91
	92	IF( kt == nit000 .AND. lwp ) THEN
[8426]	93	WRITE(numout,*)
[8512]	94	WRITE(numout,*) 'ice_thd: sea-ice thermodynamics'
	95	WRITE(numout,*) '~~~~~~~'
[8422]	96	ENDIF
	97
[8424]	98	CALL ice_var_glo2eqv
[8422]	99
	100	!---------------------------------------------!
	101	! computation of friction velocity at T points
	102	!---------------------------------------------!
[8512]	103	IF( ln_icedyn ) THEN
[8422]	104	zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
	105	zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
	106	DO jj = 2, jpjm1
	107	DO ji = fs_2, fs_jpim1
	108	zfric(ji,jj) = rn_cio * ( 0.5_wp * &
	109	& ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
	110	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
	111	END DO
	112	END DO
	113	ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean
	114	DO jj = 2, jpjm1
	115	DO ji = fs_2, fs_jpim1
	116	zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp * &
	117	& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
	118	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
	119	END DO
	120	END DO
	121	ENDIF
	122	CALL lbc_lnk( zfric, 'T', 1. )
	123	!
	124	ftr_ice(:,:,:) = 0._wp ! initialization (part of solar radiation transmitted through the ice)
	125
	126	!--------------------------------------------------------------------!
	127	! Partial computation of forcing for the thermodynamic sea ice model
	128	!--------------------------------------------------------------------!
	129	DO jj = 1, jpj
	130	DO ji = 1, jpi
	131	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
	132	!
	133	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	134	! ! practically no "direct lateral ablation"
	135	!
	136	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	137	! ! temperature and turbulent mixing (McPhee, 1992)
	138	!
	139	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
	140	zqld = tmask(ji,jj,1) * rdt_ice * &
	141	& ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &
	142	& ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
	143
	144	! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
	145	! includes supercooling potential energy (>0) or "above-freezing" energy (<0)
	146	zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
	147
	148	! --- Above-freezing sensible heat content (J/m2 grid)
	149	zqfr_neg = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp )
	150
	151	! --- Sensible ocean-to-ice heat flux (W/m2)
	152	zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
	153	fhtur(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
	154
	155	fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
	156	! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
	157	! the freezing point, so that we do not have SST < T_freeze
	158	! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
	159
	160	!-- Energy Budget of the leads (J.m-2), source of lateral accretion. Must be < 0 to form ice
	161	qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
	162
	163	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
	164	IF( zqld > 0._wp ) THEN
	165	fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
	166	qlead(ji,jj) = 0._wp
	167	ELSE
	168	fhld (ji,jj) = 0._wp
	169	ENDIF
	170	!
	171	! Net heat flux on top of the ice-ocean [W.m-2]
	172	! ---------------------------------------------
	173	hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
	174	END DO
	175	END DO
	176
	177	! In case we bypass open-water ice formation
[8514]	178	IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp
[8422]	179	! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
[8514]	180	IF( .NOT. ln_icedH ) THEN
[8483]	181	hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
	182	fhtur (:,:) = 0._wp
	183	fhld (:,:) = 0._wp
	184	ENDIF
[8422]	185
	186	! ---------------------------------------------------------------------
	187	! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
	188	! ---------------------------------------------------------------------
	189	! First step here : non solar + precip - qlead - qturb
	190	! Second step in icethd_dh : heat remaining if total melt (zq_rema)
	191	! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar
[8498]	192	hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean
	193	& - qlead(:,:) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation
	194	& - at_i (:,:) * fhtur(:,:) & ! heat flux taken by turbulence
	195	& - at_i (:,:) * fhld(:,:) ! heat flux taken during bottom growth/melt
	196	! (fhld should be 0 while bott growth)
[8422]	197	!-------------------------------------------------------------------------------------------!
	198	! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
	199	!-------------------------------------------------------------------------------------------!
	200	DO jl = 1, jpl
	201
	202	! select ice covered grid points
	203	nidx = 0 ; idxice(:) = 0
	204	DO jj = 1, jpj
	205	DO ji = 1, jpi
	206	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
	207	nidx = nidx + 1
	208	idxice(nidx) = (jj - 1) * jpi + ji
	209	ENDIF
	210	END DO
	211	END DO
	212
	213	IF( lk_mpp ) CALL mpp_ini_ice( nidx , numout )
	214
	215	IF( nidx > 0 ) THEN ! If there is no ice, do nothing.
	216	!
	217	CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- !
	218	! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
	219	!
	220	s_i_new (1:nidx) = 0._wp ; dh_s_tot (1:nidx) = 0._wp ! --- some init --- ! (important to have them here)
	221	dh_i_surf (1:nidx) = 0._wp ; dh_i_bott(1:nidx) = 0._wp
	222	dh_snowice(1:nidx) = 0._wp ; dh_i_sub (1:nidx) = 0._wp
	223	!
[8514]	224	IF( ln_icedH ) THEN ! --- growing/melting --- !
[8486]	225	CALL ice_thd_dif ! Ice/Snow Temperature profile
	226	CALL ice_thd_dh ! Ice/Snow thickness
	227	CALL ice_thd_ent( e_i_1d(1:nidx,:) ) ! Ice enthalpy remapping
	228	ENDIF
[8422]	229	!
	230	CALL ice_thd_sal ! --- Ice salinity --- !
	231	!
	232	CALL ice_thd_temp ! --- temperature update --- !
	233	!
[8486]	234	!!gm please create a new logical (l_thd_lam or a better explicit name) set one for all in icestp.F90 module
[8514]	235	!!gm l_thd_lam = ln_icedH .AND. ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 )
[8486]	236	!!gm by the way, the different options associated with nn_monocat =1 to 4 are quite impossible to identify
	237	!!gm more comment to add when ready the namelist, with an explicit print in the ocean.output
[8514]	238	IF( ln_icedH ) THEN
[8422]	239	IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
	240	CALL ice_thd_lam ! --- extra lateral melting if monocat --- !
	241	END IF
	242	END IF
	243	!
[8514]	244	IF( ln_icedA ) CALL ice_thd_da ! --- lateral melting --- !
[8422]	245	!
	246	CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
	247	! ! --- & Move to 2D arrays --- !
	248	!
	249	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	250	ENDIF
	251	!
	252	END DO
	253	! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting)
	254	oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
	255
[8517]	256	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[8422]	257	!
[8486]	258	CALL ice_var_zapsmall ! --- remove very small ice concentration (<1e-10) --- !
	259	! ! & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
[8422]	260	!
[8486]	261	IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
[8422]	262	!
[8514]	263	IF( ln_icedO ) CALL ice_thd_lac ! --- frazil ice growing in leads --- !
[8422]	264	!
[8517]	265	! controls
	266	IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
	267	IF( ln_ctl ) CALL ice_prt3D ('icethd') ! prints
	268	IF( nn_timing == 1 ) CALL timing_stop('icethd') ! timing
[8422]	269	!
	270	END SUBROUTINE ice_thd
	271
	272
	273	SUBROUTINE ice_thd_temp
	274	!!-----------------------------------------------------------------------
	275	!! * ROUTINE ice_thd_temp *
	276	!!
	277	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
	278	!!
	279	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	280	!!-------------------------------------------------------------------
	281	INTEGER :: ji, jk ! dummy loop indices
[8522]	282	REAL(wp) :: ztmelts, zbbb, zccc ! local scalar
[8422]	283	!!-------------------------------------------------------------------
	284	! Recover ice temperature
	285	DO jk = 1, nlay_i
	286	DO ji = 1, nidx
[8498]	287	ztmelts = -tmut * s_i_1d(ji,jk)
[8422]	288	! Conversion q(S,T) -> T (second order equation)
[8498]	289	zbbb = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus
	290	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) )
[8522]	291	t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic
[8422]	292
	293	! mask temperature
[8498]	294	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )
	295	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
[8422]	296	END DO
	297	END DO
	298	!
	299	END SUBROUTINE ice_thd_temp
	300
	301
	302	SUBROUTINE ice_thd_lam
	303	!!-----------------------------------------------------------------------
	304	!! * ROUTINE ice_thd_lam *
	305	!!
	306	!! ** Purpose : Lateral melting in case monocategory
	307	!! ( dA = A/2h dh )
	308	!!-----------------------------------------------------------------------
	309	INTEGER :: ji ! dummy loop indices
	310	REAL(wp) :: zhi_bef ! ice thickness before thermo
	311	REAL(wp) :: zdh_mel, zda_mel ! net melting
	312	REAL(wp) :: zvi, zvs ! ice/snow volumes
	313	!!-----------------------------------------------------------------------
	314	!
	315	DO ji = 1, nidx
	316	zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
	317	IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN
	318	zvi = a_i_1d(ji) * ht_i_1d(ji)
	319	zvs = a_i_1d(ji) * ht_s_1d(ji)
	320	! lateral melting = concentration change
	321	zhi_bef = ht_i_1d(ji) - zdh_mel
	322	rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
	323	zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
	324	a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel )
	325	! adjust thickness
	326	ht_i_1d(ji) = zvi / a_i_1d(ji)
	327	ht_s_1d(ji) = zvs / a_i_1d(ji)
	328	! retrieve total concentration
	329	at_i_1d(ji) = a_i_1d(ji)
	330	END IF
	331	END DO
	332	!
	333	END SUBROUTINE ice_thd_lam
	334
	335
[8486]	336	SUBROUTINE ice_thd_1d2d( kl, kn )
[8422]	337	!!-----------------------------------------------------------------------
	338	!! * ROUTINE ice_thd_1d2d *
	339	!!
	340	!! ** Purpose : move arrays from 1d to 2d and the reverse
	341	!!-----------------------------------------------------------------------
[8486]	342	INTEGER, INTENT(in) :: kl ! index of the ice category
	343	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
[8422]	344	!
[8486]	345	INTEGER :: jk ! dummy loop indices
[8422]	346	!!-----------------------------------------------------------------------
	347	!
	348	SELECT CASE( kn )
[8486]	349	! !---------------------!
	350	CASE( 1 ) !== from 2D to 1D ==!
	351	! !---------------------!
[8422]	352	CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i )
[8486]	353	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl) )
	354	CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d(1:nidx), ht_i(:,:,kl) )
	355	CALL tab_2d_1d( nidx, idxice(1:nidx), ht_s_1d(1:nidx), ht_s(:,:,kl) )
	356	CALL tab_2d_1d( nidx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl) )
	357	CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,kl) )
[8422]	358	DO jk = 1, nlay_s
[8486]	359	CALL tab_2d_1d( nidx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl) )
	360	CALL tab_2d_1d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl) )
[8422]	361	END DO
	362	DO jk = 1, nlay_i
[8486]	363	CALL tab_2d_1d( nidx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl) )
	364	CALL tab_2d_1d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl) )
	365	CALL tab_2d_1d( nidx, idxice(1:nidx), s_i_1d(1:nidx,jk), s_i(:,:,jk,kl) )
[8422]	366	END DO
	367	!
	368	CALL tab_2d_1d( nidx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice )
[8486]	369	CALL tab_2d_1d( nidx, idxice(1:nidx), qsr_ice_1d (1:nidx), qsr_ice (:,:,kl) )
[8422]	370	CALL tab_2d_1d( nidx, idxice(1:nidx), fr1_i0_1d (1:nidx), fr1_i0 )
	371	CALL tab_2d_1d( nidx, idxice(1:nidx), fr2_i0_1d (1:nidx), fr2_i0 )
[8486]	372	CALL tab_2d_1d( nidx, idxice(1:nidx), qns_ice_1d (1:nidx), qns_ice (:,:,kl) )
	373	CALL tab_2d_1d( nidx, idxice(1:nidx), ftr_ice_1d (1:nidx), ftr_ice (:,:,kl) )
	374	CALL tab_2d_1d( nidx, idxice(1:nidx), evap_ice_1d (1:nidx), evap_ice(:,:,kl) )
	375	CALL tab_2d_1d( nidx, idxice(1:nidx), dqns_ice_1d (1:nidx), dqns_ice(:,:,kl) )
[8422]	376	CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d (1:nidx), t_bo )
	377	CALL tab_2d_1d( nidx, idxice(1:nidx), sprecip_1d (1:nidx), sprecip )
	378	CALL tab_2d_1d( nidx, idxice(1:nidx), fhtur_1d (1:nidx), fhtur )
	379	CALL tab_2d_1d( nidx, idxice(1:nidx), fhld_1d (1:nidx), fhld )
	380	!
	381	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni )
	382	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum )
	383	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sub_1d (1:nidx), wfx_sub )
	384	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub )
	385	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub )
	386	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub )
	387	!
	388	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog )
	389	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom )
	390	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum )
	391	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni )
	392	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res )
	393	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr )
	394	CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam )
	395	!
	396	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog )
	397	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom )
	398	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum )
	399	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni )
	400	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri )
	401	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res )
	402	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub )
	403	CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam )
	404	!
	405	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd )
	406	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr )
	407	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum )
	408	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom )
	409	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog )
	410	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif )
	411	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw )
	412	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw )
	413	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub )
	414	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res )
	415	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )
	416	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem )
	417	CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out )
	418	!
	419	! SIMIP diagnostics
	420	CALL tab_2d_1d( nidx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo )
	421	CALL tab_2d_1d( nidx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su )
	422	! ocean surface fields
	423	CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m )
	424	CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m )
	425
	426	! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
	427	DO jk = 1, nlay_i
	428	WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i
[8486]	429	END DO
[8422]	430	DO jk = 1, nlay_s
	431	WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s
[8486]	432	END DO
[8422]	433	!
[8486]	434	! !---------------------!
	435	CASE( 2 ) !== from 1D to 2D ==!
	436	! !---------------------!
[8422]	437	! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
	438	DO jk = 1, nlay_i
	439	e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i
[8486]	440	END DO
[8422]	441	DO jk = 1, nlay_s
	442	e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s
[8486]	443	END DO
[8422]	444	!
[8506]	445	! Change thickness to volume (replaces routine ice_var_eqv2glo)
[8422]	446	v_i_1d(1:nidx) = ht_i_1d(1:nidx) * a_i_1d(1:nidx)
	447	v_s_1d(1:nidx) = ht_s_1d(1:nidx) * a_i_1d(1:nidx)
	448	smv_i_1d(1:nidx) = sm_i_1d(1:nidx) * v_i_1d(1:nidx)
	449
	450	CALL tab_1d_2d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i )
[8486]	451	CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl) )
	452	CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d(1:nidx), ht_i(:,:,kl) )
	453	CALL tab_1d_2d( nidx, idxice(1:nidx), ht_s_1d(1:nidx), ht_s(:,:,kl) )
	454	CALL tab_1d_2d( nidx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl) )
	455	CALL tab_1d_2d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,kl) )
[8422]	456	DO jk = 1, nlay_s
[8486]	457	CALL tab_1d_2d( nidx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl) )
	458	CALL tab_1d_2d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl) )
[8422]	459	END DO
	460	DO jk = 1, nlay_i
[8486]	461	CALL tab_1d_2d( nidx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl) )
	462	CALL tab_1d_2d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl) )
	463	CALL tab_1d_2d( nidx, idxice(1:nidx), s_i_1d(1:nidx,jk), s_i(:,:,jk,kl) )
[8422]	464	END DO
	465	!
[8486]	466	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni )
	467	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum )
	468	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sub_1d (1:nidx), wfx_sub )
	469	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub )
	470	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub )
	471	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub )
[8422]	472	!
[8486]	473	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog )
	474	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom )
	475	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum )
	476	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni )
	477	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res )
	478	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr )
	479	CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam )
[8422]	480	!
[8486]	481	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog )
	482	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom )
	483	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum )
	484	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni )
	485	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri )
	486	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res )
	487	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub )
	488	CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam )
[8422]	489	!
[8486]	490	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd )
	491	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr )
	492	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum )
	493	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom )
	494	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog )
	495	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif )
	496	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw )
	497	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw )
	498	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub )
	499	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res )
	500	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )
	501	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem )
	502	CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out )
[8422]	503	!
[8486]	504	CALL tab_1d_2d( nidx, idxice(1:nidx), qns_ice_1d(1:nidx), qns_ice(:,:,kl) )
	505	CALL tab_1d_2d( nidx, idxice(1:nidx), ftr_ice_1d(1:nidx), ftr_ice(:,:,kl) )
[8422]	506	!
	507	! SIMIP diagnostics
[8486]	508	CALL tab_1d_2d( nidx, idxice(1:nidx), t_si_1d (1:nidx), t_si(:,:,kl) )
	509	CALL tab_1d_2d( nidx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo )
	510	CALL tab_1d_2d( nidx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su )
[8422]	511	! extensive variables
[8486]	512	CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i (:,:,kl) )
	513	CALL tab_1d_2d( nidx, idxice(1:nidx), v_s_1d (1:nidx), v_s (:,:,kl) )
	514	CALL tab_1d_2d( nidx, idxice(1:nidx), smv_i_1d(1:nidx), smv_i(:,:,kl) )
	515	!
[8422]	516	END SELECT
	517	!
	518	END SUBROUTINE ice_thd_1d2d
	519
	520
	521	SUBROUTINE ice_thd_init
	522	!!-----------------------------------------------------------------------
	523	!! * ROUTINE ice_thd_init *
	524	!!
	525	!! ** Purpose : Physical constants and parameters linked to the ice
	526	!! thermodynamics
	527	!!
[8512]	528	!! ** Method : Read the namice_thd namelist and check the ice-thermo
[8422]	529	!! parameter values called at the first timestep (nit000)
	530	!!
	531	!! ** input : Namelist namicether
	532	!!-------------------------------------------------------------------
[8486]	533	INTEGER :: ios ! Local integer output status for namelist read
	534	!!
[8517]	535	NAMELIST/namice_thd/ rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s, &
[8515]	536	& ln_icedH, rn_blow_s, &
	537	& ln_icedA, rn_beta, rn_dmin, &
	538	& ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, &
	539	& nn_iceflx
[8422]	540	!!-------------------------------------------------------------------
	541	!
[8512]	542	REWIND( numnam_ice_ref ) ! Namelist namice_thd in reference namelist : Ice thermodynamics
	543	READ ( numnam_ice_ref, namice_thd, IOSTAT = ios, ERR = 901)
	544	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in reference namelist', lwp )
[8422]	545
[8512]	546	REWIND( numnam_ice_cfg ) ! Namelist namice_thd in configuration namelist : Ice thermodynamics
	547	READ ( numnam_ice_cfg, namice_thd, IOSTAT = ios, ERR = 902 )
	548	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in configuration namelist', lwp )
	549	IF(lwm) WRITE ( numoni, namice_thd )
[8422]	550	!
	551	!
	552	IF(lwp) THEN ! control print
	553	WRITE(numout,*) 'ice_thd_init : Ice Thermodynamics'
	554	WRITE(numout,*) '~~~~~~~~~~~~~'
[8512]	555	WRITE(numout,*) ' Namelist namice_thd'
	556	WRITE(numout,*) ' -- icethd_dif --'
	557	WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i
[8514]	558	WRITE(numout,*) ' thermal conductivity in the ice (Untersteiner 1964) ln_cndi_U64 = ', ln_cndi_U64
	559	WRITE(numout,*) ' thermal conductivity in the ice (Pringle et al 2007) ln_cndi_P07 = ', ln_cndi_P07
	560	WRITE(numout,*) ' change the surface non-solar flux with Tsu or not ln_dqns_i = ', ln_dqns_i
	561	WRITE(numout,*) ' thermal conductivity of the snow rn_cnd_s = ', rn_cnd_s
[8512]	562	WRITE(numout,*) ' -- icethd_dh --'
[8514]	563	WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH
	564	WRITE(numout,*) ' coefficient for ice-lead partition of snowfall rn_blow_s = ', rn_blow_s
[8512]	565	WRITE(numout,*) ' -- icethd_da --'
[8514]	566	WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA
[8512]	567	WRITE(numout,*) ' Coef. beta for lateral melting param. rn_beta = ', rn_beta
	568	WRITE(numout,*) ' Minimum floe diameter for lateral melting param. rn_dmin = ', rn_dmin
	569	WRITE(numout,*) ' -- icethd_lac --'
[8514]	570	WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO
	571	WRITE(numout,*) ' ice thickness for lateral accretion rn_hinew = ', rn_hinew
[8512]	572	WRITE(numout,*) ' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
[8514]	573	WRITE(numout,*) ' Maximum proportion of frazil ice collecting at bottom rn_maxfraz = ', rn_maxfraz
	574	WRITE(numout,*) ' Threshold relative drift speed for collection of frazil rn_vfraz = ', rn_vfraz
	575	WRITE(numout,*) ' Squeezing coefficient for collection of frazil rn_Cfraz = ', rn_Cfraz
[8512]	576	WRITE(numout,*) ' -- icestp --'
[8514]	577	WRITE(numout,*) ' Multicategory heat flux formulation nn_iceflx = ', nn_iceflx
[8422]	578	ENDIF
	579	!
[8514]	580	IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
	581	CALL ctl_stop( 'ice_thd_init: choose one and only one formulation for thermal conductivity (ln_cndi_U64 or ln_cndi_P07)' )
	582	ENDIF
[8422]	583	!
[8518]	584	IF ( rn_hinew < rn_himin ) CALL ctl_stop( 'ice_thd_init : rn_hinew should be >= rn_himin' )
[8514]	585	!
[8422]	586	IF(lwp) WRITE(numout,*)
[8514]	587	SELECT CASE( nn_iceflx ) ! LIM3 Multi-category heat flux formulation
[8486]	588	CASE( -1 )
[8514]	589	IF(lwp) WRITE(numout,*) ' LIM3: use per-category fluxes (nn_iceflx = -1) '
	590	IF( ln_cpl ) CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
[8486]	591	CASE( 0 )
[8514]	592	IF(lwp) WRITE(numout,*) ' LIM3: use average per-category fluxes (nn_iceflx = 0) '
[8486]	593	CASE( 1 )
[8514]	594	IF(lwp) WRITE(numout,*) ' LIM3: use average then redistribute per-category fluxes (nn_iceflx = 1) '
	595	IF( ln_cpl ) CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
[8486]	596	CASE( 2 )
[8514]	597	IF(lwp) WRITE(numout,*) ' LIM3: Redistribute a single flux over categories (nn_iceflx = 2) '
	598	IF( .NOT. ln_cpl ) CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in forced mode cannot be 2' )
[8422]	599	CASE DEFAULT
[8514]	600	CALL ctl_stop( 'ice_thd_init: LIM3 option, nn_iceflx, should be between -1 and 2' )
[8422]	601	END SELECT
	602	!
	603	END SUBROUTINE ice_thd_init
	604
	605	#else
	606	!!----------------------------------------------------------------------
	607	!! Default option Dummy module NO LIM3 sea-ice model
	608	!!----------------------------------------------------------------------
	609	#endif
	610
	611	!!======================================================================
	612	END MODULE icethd

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

Download in other formats: