Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

limthd.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 34.5 KB

Line
1	MODULE limthd
2	!!======================================================================
3	!! * MODULE limthd *
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 lim_thd_glohec, lim_thd_con_dh and lim_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	!! lim_thd : thermodynamic of sea ice
20	!! lim_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	USE sbc_oce ! Surface boundary condition: ocean fields
26	USE sbc_ice ! Surface boundary condition: ice fields
27	USE thd_ice ! thermodynamic sea-ice variables
28	USE limthd_dif ! vertical diffusion
29	USE limthd_dh ! ice-snow growth and melt
30	USE limthd_da ! lateral melting
31	USE limthd_sal ! ice salinity
32	USE limthd_ent ! ice enthalpy redistribution
33	USE limthd_lac ! lateral accretion
34	USE limitd_th ! remapping thickness distribution
35	USE limtab ! 1D <==> 2D transformation
36	USE limvar !
37	USE limcons ! conservation tests
38	USE limctl ! control print
39	!
40	USE in_out_manager ! I/O manager
41	USE lbclnk ! lateral boundary condition - MPP links
42	USE lib_mpp ! MPP library
43	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
44	USE timing ! Timing
45
46	IMPLICIT NONE
47	PRIVATE
48
49	PUBLIC lim_thd ! called by limstp module
50	PUBLIC lim_thd_init ! called by sbc_lim_init
51
52	!! * Substitutions
53	# include "vectopt_loop_substitute.h90"
54	!!----------------------------------------------------------------------
55	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
56	!! $Id$
57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
58	!!----------------------------------------------------------------------
59	CONTAINS
60
61	SUBROUTINE lim_thd( kt )
62	!!-------------------------------------------------------------------
63	!! * ROUTINE lim_thd *
64	!!
65	!! ** Purpose : This routine manages ice thermodynamics
66	!!
67	!! ** Action : - Initialisation of some variables
68	!! - Some preliminary computation (oceanic heat flux
69	!! at the ice base, snow acc.,heat budget of the leads)
70	!! - selection of the icy points and put them in an array
71	!! - call lim_thd_dif for vertical heat diffusion
72	!! - call lim_thd_dh for vertical ice growth and melt
73	!! - call lim_thd_ent for enthalpy remapping
74	!! - call lim_thd_sal for ice desalination
75	!! - call lim_thd_temp to retrieve temperature from ice enthalpy
76	!! - back to the geographic grid
77	!!
78	!! ** References :
79	!!---------------------------------------------------------------------
80	INTEGER, INTENT(in) :: kt ! number of iteration
81	!
82	INTEGER :: ji, jj, jk, jl ! dummy loop indices
83	INTEGER :: nbpb ! nb of icy pts for vertical thermo calculations
84	REAL(wp) :: zfric_u, zqld, zqfr
85	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
86	REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
87	REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
88	REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
89	!
90	!!-------------------------------------------------------------------
91
92	IF( nn_timing == 1 ) CALL timing_start('limthd')
93
94
95	IF( kt == nit000 .AND. lwp ) THEN
96	WRITE(numout,*)''
97	WRITE(numout,*)' lim_thd '
98	WRITE(numout,*)' ~~~~~~~~'
99	ENDIF
100
101	! conservation test
102	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
103
104	CALL lim_var_glo2eqv
105
106	!---------------------------------------------!
107	! computation of friction velocity at T points
108	!---------------------------------------------!
109	IF( ln_limdyn ) THEN
110	zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
111	zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
112	DO jj = 2, jpjm1
113	DO ji = fs_2, fs_jpim1
114	zfric(ji,jj) = rn_cio * ( 0.5_wp * &
115	& ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
116	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
117	END DO
118	END DO
119	ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean
120	DO jj = 2, jpjm1
121	DO ji = fs_2, fs_jpim1
122	zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp * &
123	& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
124	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
125	END DO
126	END DO
127	ENDIF
128	CALL lbc_lnk( zfric, 'T', 1. )
129	!
130	!----------------------------------!
131	! Initialization and units change
132	!----------------------------------!
133	ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice
134
135	! Change the units of heat content; from J/m2 to J/m3
136	DO jl = 1, jpl
137	DO jk = 1, nlay_i
138	DO jj = 1, jpj
139	DO ji = 1, jpi
140	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) )
141	!Energy of melting q(S,T) [J.m-3]
142	e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
143	END DO
144	END DO
145	END DO
146	DO jk = 1, nlay_s
147	DO jj = 1, jpj
148	DO ji = 1, jpi
149	rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )
150	!Energy of melting q(S,T) [J.m-3]
151	e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
152	END DO
153	END DO
154	END DO
155	END DO
156
157	!--------------------------------------------------------------------!
158	! Partial computation of forcing for the thermodynamic sea ice model
159	!--------------------------------------------------------------------!
160	DO jj = 1, jpj
161	DO ji = 1, jpi
162	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
163	!
164	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
165	! ! practically no "direct lateral ablation"
166	!
167	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
168	! ! temperature and turbulent mixing (McPhee, 1992)
169	!
170	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
171	zqld = tmask(ji,jj,1) * rdt_ice * &
172	& ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
173
174	! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
175	zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
176
177	! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
178	zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
179	fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
180	fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
181	! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
182	! the freezing point, so that we do not have SST < T_freeze
183	! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
184
185	!-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
186	qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
187
188	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
189	IF( zqld > 0._wp ) THEN
190	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 limthd_dh.F90
191	qlead(ji,jj) = 0._wp
192	ELSE
193	fhld (ji,jj) = 0._wp
194	ENDIF
195	!
196	! Net heat flux on top of the ice-ocean [W.m-2]
197	! ---------------------------------------------
198	hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
199	END DO
200	END DO
201
202	! In case we bypass open-water ice formation
203	IF( .NOT. ln_limdO ) qlead(:,:) = 0._wp
204	! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
205	IF( .NOT. ln_limdH ) hfx_in(:,:) = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
206	IF( .NOT. ln_limdH ) fhtur (:,:) = 0._wp ; fhld (:,:) = 0._wp
207
208	! ---------------------------------------------------------------------
209	! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
210	! ---------------------------------------------------------------------
211	! First step here : non solar + precip - qlead - qturb
212	! Second step in limthd_dh : heat remaining if total melt (zq_rema)
213	! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar
214	DO jj = 1, jpj
215	DO ji = 1, jpi
216	hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean
217	& - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation
218	& - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence
219	& - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt
220	! (fhld should be 0 while bott growth)
221	END DO
222	END DO
223
224	!------------------------------------------------------------------------------!
225	! Thermodynamic computation (only on grid points covered by ice)
226	!------------------------------------------------------------------------------!
227
228	DO jl = 1, jpl !loop over ice categories
229
230	! select ice covered grid points
231	nbpb = 0
232	DO jj = 1, jpj
233	DO ji = 1, jpi
234	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
235	nbpb = nbpb + 1
236	npb(nbpb) = (jj - 1) * jpi + ji
237	ENDIF
238	END DO
239	END DO
240
241	! debug point to follow
242	jiindex_1d = 0
243	IF( ln_limctl ) THEN
244	DO ji = mi0(iiceprt), mi1(iiceprt)
245	DO jj = mj0(jiceprt), mj1(jiceprt)
246	jiindex_1d = (jj - 1) * jpi + ji
247	WRITE(numout,*) ' lim_thd : Category no : ', jl
248	END DO
249	END DO
250	ENDIF
251
252	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
253
254	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
255	!
256	s_i_new (:) = 0._wp ; dh_s_tot (:) = 0._wp ! --- some init --- !
257	dh_i_surf (:) = 0._wp ; dh_i_bott(:) = 0._wp
258	dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp
259
260	CALL lim_thd_1d2d( nbpb, jl, 1 ) ! --- Move to 1D arrays --- !
261	!
262	IF( ln_limdH ) CALL lim_thd_dif( 1, nbpb ) ! --- Ice/Snow Temperature profile --- !
263	!
264	IF( ln_limdH ) CALL lim_thd_dh( 1, nbpb ) ! --- Ice/Snow thickness --- !
265	!
266	IF( ln_limdH ) CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) ! --- Ice enthalpy remapping --- !
267	!
268	CALL lim_thd_sal( 1, nbpb ) ! --- Ice salinity --- !
269	!
270	CALL lim_thd_temp( 1, nbpb ) ! --- temperature update --- !
271	!
272	IF( ln_limdH ) THEN
273	IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
274	CALL lim_thd_lam( 1, nbpb ) ! --- extra lateral melting if monocat --- !
275	END IF
276	END IF
277	!
278	CALL lim_thd_1d2d( nbpb, jl, 2 ) ! --- Move to 2D arrays --- !
279	!
280	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
281	ENDIF
282	!
283	END DO !jl
284
285	IF( ln_limdA) CALL lim_thd_da ! --- lateral melting --- !
286
287	! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
288	DO jl = 1, jpl
289	DO jk = 1, nlay_i
290	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
291	END DO
292	DO jk = 1, nlay_s
293	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
294	END DO
295	END DO
296
297	! Change thickness to volume
298	v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
299	v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
300	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
301
302	! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat)
303	DO jl = 1, jpl
304	DO jj = 1, jpj
305	DO ji = 1, jpi
306	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) )
307	oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 )
308	END DO
309	END DO
310	END DO
311
312	CALL lim_var_zapsmall
313
314	! control checks
315	IF( ln_limctl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' ) ! control print
316	!
317	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
318
319	!------------------------------------------------!
320	! Transport ice between thickness categories
321	!------------------------------------------------!
322	! Given thermodynamic growth rates, transport ice between thickness categories.
323	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
324
325	IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )
326
327	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
328
329	!------------------------------------------------!
330	! Add frazil ice growing in leads
331	!------------------------------------------------!
332	IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
333
334	IF( ln_limdO ) CALL lim_thd_lac
335
336	! conservation test
337	IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
338
339	! Control print
340	IF( ln_ctl ) CALL lim_prt3D( 'limthd' )
341	!
342	!
343	IF( nn_timing == 1 ) CALL timing_stop('limthd')
344
345	END SUBROUTINE lim_thd
346
347
348	SUBROUTINE lim_thd_temp( kideb, kiut )
349	!!-----------------------------------------------------------------------
350	!! * ROUTINE lim_thd_temp *
351	!!
352	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
353	!!
354	!! ** Method : Formula (Bitz and Lipscomb, 1999)
355	!!-------------------------------------------------------------------
356	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
357	!
358	INTEGER :: ji, jk ! dummy loop indices
359	REAL(wp) :: ztmelts, zaaa, zbbb, zccc, zdiscrim ! local scalar
360	!!-------------------------------------------------------------------
361	! Recover ice temperature
362	DO jk = 1, nlay_i
363	DO ji = kideb, kiut
364	ztmelts = -tmut * s_i_1d(ji,jk) + rt0
365	! Conversion q(S,T) -> T (second order equation)
366	zaaa = cpic
367	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
368	zccc = lfus * ( ztmelts - rt0 )
369	zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
370	t_i_1d(ji,jk) = rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
371
372	! mask temperature
373	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )
374	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
375	END DO
376	END DO
377	!
378	END SUBROUTINE lim_thd_temp
379
380
381	SUBROUTINE lim_thd_lam( kideb, kiut )
382	!!-----------------------------------------------------------------------
383	!! * ROUTINE lim_thd_lam *
384	!!
385	!! ** Purpose : Lateral melting in case monocategory
386	!! ( dA = A/2h dh )
387	!!-----------------------------------------------------------------------
388	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
389	!
390	INTEGER :: ji ! dummy loop indices
391	REAL(wp) :: zhi_bef ! ice thickness before thermo
392	REAL(wp) :: zdh_mel, zda_mel ! net melting
393	REAL(wp) :: zvi, zvs ! ice/snow volumes
394	!!-----------------------------------------------------------------------
395	!
396	DO ji = kideb, kiut
397	zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
398	IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN
399	zvi = a_i_1d(ji) * ht_i_1d(ji)
400	zvs = a_i_1d(ji) * ht_s_1d(ji)
401	! lateral melting = concentration change
402	zhi_bef = ht_i_1d(ji) - zdh_mel
403	rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
404	zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
405	a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel )
406	! adjust thickness
407	ht_i_1d(ji) = zvi / a_i_1d(ji)
408	ht_s_1d(ji) = zvs / a_i_1d(ji)
409	! retrieve total concentration
410	at_i_1d(ji) = a_i_1d(ji)
411	END IF
412	END DO
413	!
414	END SUBROUTINE lim_thd_lam
415
416
417	SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
418	!!-----------------------------------------------------------------------
419	!! * ROUTINE lim_thd_1d2d *
420	!!
421	!! ** Purpose : move arrays from 1d to 2d and the reverse
422	!!-----------------------------------------------------------------------
423	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
424	INTEGER, INTENT(in) :: nbpb ! size of 1D arrays
425	INTEGER, INTENT(in) :: jl ! ice cat
426	!
427	INTEGER :: jk ! dummy loop indices
428	!!-----------------------------------------------------------------------
429	!
430	SELECT CASE( kn )
431	!
432	CASE( 1 ) ! from 2D to 1D
433	!
434	CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
435	CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
436	CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
437	CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
438	!
439	CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
440	CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
441	DO jk = 1, nlay_s
442	CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
443	CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
444	END DO
445	DO jk = 1, nlay_i
446	CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
447	CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
448	CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
449	END DO
450	!
451	CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
452	CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
453	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
454	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
455	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
456	CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
457	CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
458	CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
459	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
460	CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
461	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
462	CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) )
463	CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) )
464	CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) )
465	!
466	CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) )
467	CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) )
468	!
469	CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) )
470	CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) )
471	CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) )
472	CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) )
473	CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) )
474	CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) )
475	!
476	CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) )
477	CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) )
478	CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) )
479	CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) )
480	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
481	CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) )
482	CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub , jpi, jpj,npb(1:nbpb) )
483	!
484	CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) )
485	CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) )
486	CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) )
487	CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) )
488	CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) )
489	CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) )
490	CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) )
491	CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) )
492	CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) )
493	CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) )
494	CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) )
495	CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
496	CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
497	!
498	CASE( 2 ) ! from 1D to 2D
499	!
500	CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj )
501	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj )
502	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj )
503	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj )
504	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj )
505	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj )
506	DO jk = 1, nlay_s
507	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj)
508	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj)
509	END DO
510	DO jk = 1, nlay_i
511	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj)
512	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj)
513	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj)
514	END DO
515	CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj )
516	!
517	CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj )
518	CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj )
519	!
520	CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj )
521	CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj )
522	CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj )
523	CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj )
524	CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj )
525	CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj )
526	!
527	CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj )
528	CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj )
529	CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj )
530	CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj )
531	CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj )
532	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
533	CALL tab_1d_2d( nbpb, sfx_sub , npb, sfx_sub_1d(1:nbpb) , jpi, jpj )
534	!
535	CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj )
536	CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj )
537	CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj )
538	CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj )
539	CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj )
540	CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj )
541	CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj )
542	CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj )
543	CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj )
544	CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj )
545	CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj )
546	CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
547	CALL tab_1d_2d( nbpb, hfx_err_dif , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
548	!
549	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
550	CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
551	!
552	END SELECT
553	!
554	END SUBROUTINE lim_thd_1d2d
555
556
557	SUBROUTINE lim_thd_init
558	!!-----------------------------------------------------------------------
559	!! * ROUTINE lim_thd_init *
560	!!
561	!! ** Purpose : Physical constants and parameters linked to the ice
562	!! thermodynamics
563	!!
564	!! ** Method : Read the namicethd namelist and check the ice-thermo
565	!! parameter values called at the first timestep (nit000)
566	!!
567	!! ** input : Namelist namicether
568	!!-------------------------------------------------------------------
569	INTEGER :: ios ! Local integer output status for namelist read
570	NAMELIST/namicethd/ rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, ln_it_qnsice, nn_monocat, rn_cdsn, &
571	& ln_limdH, rn_betas, &
572	& ln_limdA, rn_beta, rn_dmin, &
573	& ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin
574	!!-------------------------------------------------------------------
575	!
576	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
577	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
578	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
579
580	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
581	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
582	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
583	IF(lwm) WRITE ( numoni, namicethd )
584	!
585	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
586	nn_monocat = 0
587	IF(lwp) WRITE(numout,*)
588	IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
589	ENDIF
590	!
591	IF(lwp) THEN ! control print
592	WRITE(numout,*) 'lim_thd_init : Ice Thermodynamics'
593	WRITE(numout,*) '~~~~~~~~~~~~~'
594	WRITE(numout,*)' -- limthd_dif --'
595	WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i
596	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif
597	WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif
598	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon
599	WRITE(numout,*)' iterate the surface non-solar flux (T) or not (F) ln_it_qnsice = ', ln_it_qnsice
600	WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat
601	WRITE(numout,*)' thermal conductivity of the snow rn_cdsn = ', rn_cdsn
602	WRITE(numout,*)' -- limthd_dh --'
603	WRITE(numout,*)' activate ice thick change from top/bot (T) or not (F) ln_limdH = ', ln_limdH
604	WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas
605	WRITE(numout,*)' -- limthd_da --'
606	WRITE(numout,*)' activate lateral melting (T) or not (F) ln_limdA = ', ln_limdA
607	WRITE(numout,*)' Coef. beta for lateral melting param. rn_beta = ', rn_beta
608	WRITE(numout,*)' Minimum floe diameter for lateral melting param. rn_dmin = ', rn_dmin
609	WRITE(numout,*)' -- limthd_lac --'
610	WRITE(numout,*)' activate ice growth in open-water (T) or not (F) ln_limdO = ', ln_limdO
611	WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice
612	WRITE(numout,*)' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
613	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb
614	WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb
615	WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb
616	WRITE(numout,*)' -- limitd_th --'
617	WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin
618	WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i
619	ENDIF
620	!
621	END SUBROUTINE lim_thd_init
622
623	#else
624	!!----------------------------------------------------------------------
625	!! Default option Dummy module NO LIM3 sea-ice model
626	!!----------------------------------------------------------------------
627	#endif
628
629	!!======================================================================
630	END MODULE limthd

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 7910

Download in other formats: