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 @ 8498

Last change on this file since 8498 was 8498, checked in by clem, 7 years ago
changes in style - part2 -
File size: 35.9 KB

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

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

Download in other formats: