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 NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago
Add neighborhood collectives calls in the NEMO src - ticket #2496
Property svn:keywords set to `Id`
File size: 34.8 KB

Line
1	MODULE icethd
2	!!======================================================================
3	!! * MODULE icethd *
4	!! sea-ice : master routine for thermodynamics
5	!!======================================================================
6	!! History : 1.0 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) original code 1D
7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
8	!!----------------------------------------------------------------------
9	#if defined key_si3
10	!!----------------------------------------------------------------------
11	!! 'key_si3' SI3 sea-ice model
12	!!----------------------------------------------------------------------
13	!! ice_thd : thermodynamics of sea ice
14	!! ice_thd_init : initialisation of sea-ice thermodynamics
15	!!----------------------------------------------------------------------
16	USE phycst ! physical constants
17	USE dom_oce ! ocean space and time domain variables
18	USE ice ! sea-ice: variables
19	!!gm list trop longue ==>>> why not passage en argument d'appel ?
20	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
21	USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
22	& qml_ice, qcn_ice, qtr_ice_top
23	USE ice1D ! sea-ice: thermodynamics variables
24	USE icethd_zdf ! sea-ice: vertical heat diffusion
25	USE icethd_dh ! sea-ice: ice-snow growth and melt
26	USE icethd_da ! sea-ice: lateral melting
27	USE icethd_sal ! sea-ice: salinity
28	USE icethd_ent ! sea-ice: enthalpy redistribution
29	USE icethd_do ! sea-ice: growth in open water
30	USE icethd_pnd ! sea-ice: melt ponds
31	USE iceitd ! sea-ice: remapping thickness distribution
32	USE icetab ! sea-ice: 1D <==> 2D transformation
33	USE icevar ! sea-ice: operations
34	USE icectl ! sea-ice: control print
35	!
36	USE in_out_manager ! I/O manager
37	USE iom ! I/O manager library
38	USE lib_mpp ! MPP library
39	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
40	USE lbclnk ! lateral boundary conditions (or mpp links)
41	USE timing ! Timing
42
43	IMPLICIT NONE
44	PRIVATE
45
46	PUBLIC ice_thd ! called by limstp module
47	PUBLIC ice_thd_init ! called by ice_init
48
49	!! namelist (namthd)
50	LOGICAL :: ln_icedH ! activate ice thickness change from growing/melting (T) or not (F)
51	LOGICAL :: ln_icedA ! activate lateral melting param. (T) or not (F)
52	LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F)
53	LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F)
54	LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean
55
56	!! for convergence tests
57	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztice_cvgerr, ztice_cvgstp
58
59	!! * Substitutions
60	# include "do_loop_substitute.h90"
61	!!----------------------------------------------------------------------
62	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
63	!! $Id$
64	!! Software governed by the CeCILL license (see ./LICENSE)
65	!!----------------------------------------------------------------------
66	CONTAINS
67
68	SUBROUTINE ice_thd( kt )
69	!!-------------------------------------------------------------------
70	!! * ROUTINE ice_thd *
71	!!
72	!! ** Purpose : This routine manages ice thermodynamics
73	!!
74	!! ** Action : - computation of oceanic sensible heat flux at the ice base
75	!! energy budget in the leads
76	!! net fluxes on top of ice and of ocean
77	!! - selection of grid cells with ice
78	!! - call ice_thd_zdf for vertical heat diffusion
79	!! - call ice_thd_dh for vertical ice growth and melt
80	!! - call ice_thd_pnd for melt ponds
81	!! - call ice_thd_ent for enthalpy remapping
82	!! - call ice_thd_sal for ice desalination
83	!! - call ice_thd_temp to retrieve temperature from ice enthalpy
84	!! - call ice_thd_mono for extra lateral ice melt if active virtual thickness distribution
85	!! - call ice_thd_da for lateral ice melt
86	!! - back to the geographic grid
87	!! - call ice_thd_rem for remapping thickness distribution
88	!! - call ice_thd_do for ice growth in leads
89	!!-------------------------------------------------------------------
90	INTEGER, INTENT(in) :: kt ! number of iteration
91	!
92	INTEGER :: ji, jj, jk, jl ! dummy loop indices
93	REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
94	REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
95	REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
96	REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
97	!
98	!!-------------------------------------------------------------------
99	! controls
100	IF( ln_timing ) CALL timing_start('icethd') ! timing
101	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
102	IF( ln_icediachk ) CALL ice_cons2D (0, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
103
104	IF( kt == nit000 .AND. lwp ) THEN
105	WRITE(numout,*)
106	WRITE(numout,*) 'ice_thd: sea-ice thermodynamics'
107	WRITE(numout,*) '~~~~~~~'
108	ENDIF
109
110	! convergence tests
111	IF( ln_zdf_chkcvg ) THEN
112	ALLOCATE( ztice_cvgerr(jpi,jpj,jpl) , ztice_cvgstp(jpi,jpj,jpl) )
113	ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp
114	ENDIF
115
116	!---------------------------------------------!
117	! computation of friction velocity at T points
118	!---------------------------------------------!
119	IF( ln_icedyn ) THEN
120	zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
121	zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
122	DO_2D( 0, 0, 0, 0 )
123	zfric(ji,jj) = rn_cio * ( 0.5_wp * &
124	& ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
125	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
126	END_2D
127	ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean
128	DO_2D( 0, 0, 0, 0 )
129	zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * &
130	& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
131	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
132	END_2D
133	ENDIF
134	#if defined key_mpi3
135	CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T', 1.0_wp )
136	#else
137	CALL lbc_lnk( 'icethd', zfric, 'T', 1.0_wp )
138	#endif
139	!
140	!--------------------------------------------------------------------!
141	! Partial computation of forcing for the thermodynamic sea ice model
142	!--------------------------------------------------------------------!
143	DO_2D( 1, 1, 1, 1 )
144	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
145	!
146	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
147	! ! practically no "direct lateral ablation"
148	!
149	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
150	! ! temperature and turbulent mixing (McPhee, 1992)
151	!
152	! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
153	zqld = tmask(ji,jj,1) * rDt_ice * &
154	& ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &
155	& ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
156
157	! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- !
158	zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
159	zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0
160
161	! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2)
162	zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
163	qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
164
165	qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )
166	! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach
167	! the freezing point, so that we do not have SST < T_freeze
168	! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
169
170	!-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice
171	qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr )
172
173	! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting
174	! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
175	IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
176	IF( ln_leadhfx ) THEN ; fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
177	ELSE ; fhld(ji,jj) = 0._wp
178	ENDIF
179	qlead(ji,jj) = 0._wp
180	ELSE
181	fhld (ji,jj) = 0._wp
182	ENDIF
183	!
184	! Net heat flux on top of the ice-ocean [W.m-2]
185	! ---------------------------------------------
186	qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
187	END_2D
188
189	! In case we bypass open-water ice formation
190	IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp
191	! In case we bypass growing/melting from top and bottom
192	IF( .NOT. ln_icedH ) THEN
193	qsb_ice_bot(:,:) = 0._wp
194	fhld (:,:) = 0._wp
195	ENDIF
196
197	! ---------------------------------------------------------------------
198	! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
199	! ---------------------------------------------------------------------
200	! First step here : non solar + precip - qlead - qsensible
201	! Second step in icethd_dh : heat remaining if total melt (zq_rema)
202	! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar
203	qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean
204	& - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation
205	& - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux
206	& - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt
207	! ! (fhld should be 0 while bott growth)
208	!-------------------------------------------------------------------------------------------!
209	! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
210	!-------------------------------------------------------------------------------------------!
211	DO jl = 1, jpl
212
213	! select ice covered grid points
214	npti = 0 ; nptidx(:) = 0
215	DO_2D( 1, 1, 1, 1 )
216	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
217	npti = npti + 1
218	nptidx(npti) = (jj - 1) * jpi + ji
219	ENDIF
220	END_2D
221
222	IF( npti > 0 ) THEN ! If there is no ice, do nothing.
223	!
224	CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- !
225	! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
226	!
227	s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here)
228	dh_i_sum (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm (1:npti) = 0._wp
229	dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
230	dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
231	!
232	CALL ice_thd_zdf ! --- Ice-Snow temperature --- !
233	!
234	IF( ln_icedH ) THEN ! --- Growing/Melting --- !
235	CALL ice_thd_dh ! Ice-Snow thickness
236	CALL ice_thd_pnd ! Melt ponds formation
237	CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping
238	ENDIF
239	CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- !
240	!
241	CALL ice_thd_temp ! --- Temperature update --- !
242	!
243	IF( ln_icedH .AND. ln_virtual_itd ) &
244	& CALL ice_thd_mono ! --- Extra lateral melting if virtual_itd --- !
245	!
246	IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- !
247	!
248	CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
249	! ! --- & Move to 2D arrays --- !
250	ENDIF
251	!
252	END DO
253	!
254	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
255	IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
256	!
257	IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
258	!
259	IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- !
260	!
261	! convergence tests
262	IF( ln_zdf_chkcvg ) THEN
263	CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr )
264	CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp )
265	ENDIF
266	!
267	! controls
268	IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
269	IF( sn_cfctl%l_prtctl ) &
270	& CALL ice_prt3D ('icethd') ! prints
271	IF( ln_timing ) CALL timing_stop('icethd') ! timing
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, zbbb, zccc ! local scalar
286	!!-------------------------------------------------------------------
287	! Recover ice temperature
288	DO jk = 1, nlay_i
289	DO ji = 1, npti
290	ztmelts = -rTmlt * sz_i_1d(ji,jk)
291	! Conversion q(S,T) -> T (second order equation)
292	zbbb = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus
293	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) )
294	t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi
295
296	! mask temperature
297	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
298	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
299	END DO
300	END DO
301	!
302	END SUBROUTINE ice_thd_temp
303
304
305	SUBROUTINE ice_thd_mono
306	!!-----------------------------------------------------------------------
307	!! * ROUTINE ice_thd_mono *
308	!!
309	!! ** Purpose : Lateral melting in case virtual_itd
310	!! ( dA = A/2h dh )
311	!!-----------------------------------------------------------------------
312	INTEGER :: ji ! dummy loop indices
313	REAL(wp) :: zhi_bef ! ice thickness before thermo
314	REAL(wp) :: zdh_mel, zda_mel ! net melting
315	REAL(wp) :: zvi, zvs ! ice/snow volumes
316	!!-----------------------------------------------------------------------
317	!
318	DO ji = 1, npti
319	zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) )
320	IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN
321	zvi = a_i_1d(ji) * h_i_1d(ji)
322	zvs = a_i_1d(ji) * h_s_1d(ji)
323	! lateral melting = concentration change
324	zhi_bef = h_i_1d(ji) - zdh_mel
325	rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
326	zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
327	a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel )
328	! adjust thickness
329	h_i_1d(ji) = zvi / a_i_1d(ji)
330	h_s_1d(ji) = zvs / a_i_1d(ji)
331	! retrieve total concentration
332	at_i_1d(ji) = a_i_1d(ji)
333	END IF
334	END DO
335	!
336	END SUBROUTINE ice_thd_mono
337
338
339	SUBROUTINE ice_thd_1d2d( kl, kn )
340	!!-----------------------------------------------------------------------
341	!! * ROUTINE ice_thd_1d2d *
342	!!
343	!! ** Purpose : move arrays from 1d to 2d and the reverse
344	!!-----------------------------------------------------------------------
345	INTEGER, INTENT(in) :: kl ! index of the ice category
346	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
347	!
348	INTEGER :: jk ! dummy loop indices
349	!!-----------------------------------------------------------------------
350	!
351	SELECT CASE( kn )
352	! !---------------------!
353	CASE( 1 ) !== from 2D to 1D ==!
354	! !---------------------!
355	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i )
356	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) )
357	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) )
358	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) )
359	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) )
360	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) )
361	DO jk = 1, nlay_s
362	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
363	CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
364	END DO
365	DO jk = 1, nlay_i
366	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) )
367	CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) )
368	CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
369	END DO
370	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )
371	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
372	CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )
373	!
374	CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice )
375	CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d (1:npti), qsr_ice (:,:,kl) )
376	CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) )
377	CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d (1:npti), evap_ice(:,:,kl) )
378	CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d (1:npti), dqns_ice(:,:,kl) )
379	CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti), t_bo )
380	CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d (1:npti), sprecip )
381	CALL tab_2d_1d( npti, nptidx(1:npti), qsb_ice_bot_1d(1:npti), qsb_ice_bot )
382	CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d (1:npti), fhld )
383
384	CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d (1:npti), qml_ice (:,:,kl) )
385	CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d (1:npti), qcn_ice (:,:,kl) )
386	CALL tab_2d_1d( npti, nptidx(1:npti), qtr_ice_top_1d(1:npti), qtr_ice_top(:,:,kl) )
387	!
388	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
389	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
390	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub )
391	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
392	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
393	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
394	!
395	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog )
396	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom )
397	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
398	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni )
399	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res )
400	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr )
401	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam )
402	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
403	!
404	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog )
405	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom )
406	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum )
407	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni )
408	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri )
409	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res )
410	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub )
411	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam )
412	!
413	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd )
414	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr )
415	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum )
416	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom )
417	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog )
418	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif )
419	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw )
420	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw )
421	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub )
422	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res )
423	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
424	CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )
425	!
426	! ocean surface fields
427	CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
428	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
429	!
430	! to update ice age
431	CALL tab_2d_1d( npti, nptidx(1:npti), o_i_1d (1:npti), o_i (:,:,kl) )
432	CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
433	!
434	! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
435	DO jk = 1, nlay_i
436	WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
437	END DO
438	DO jk = 1, nlay_s
439	WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
440	END DO
441	!
442	! !---------------------!
443	CASE( 2 ) !== from 1D to 2D ==!
444	! !---------------------!
445	! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
446	DO jk = 1, nlay_i
447	e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * h_i_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_i
448	END DO
449	DO jk = 1, nlay_s
450	e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) * h_s_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_s
451	END DO
452	!
453	! Change thickness to volume (replaces routine ice_var_eqv2glo)
454	v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
455	v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
456	sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
457	v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
458	v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
459	oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
460
461	CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i )
462	CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) )
463	CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) )
464	CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) )
465	CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) )
466	CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) )
467	DO jk = 1, nlay_s
468	CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
469	CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
470	END DO
471	DO jk = 1, nlay_i
472	CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) )
473	CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) )
474	CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
475	END DO
476	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )
477	CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
478	CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )
479	!
480	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
481	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
482	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub )
483	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
484	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
485	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
486	!
487	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog )
488	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom )
489	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
490	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni )
491	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res )
492	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr )
493	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam )
494	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
495	!
496	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog )
497	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom )
498	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum )
499	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni )
500	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri )
501	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res )
502	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub )
503	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam )
504	!
505	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd )
506	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr )
507	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum )
508	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom )
509	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog )
510	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif )
511	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw )
512	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw )
513	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub )
514	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res )
515	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
516	CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )
517	!
518	CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) )
519	CALL tab_1d_2d( npti, nptidx(1:npti), qtr_ice_bot_1d(1:npti), qtr_ice_bot(:,:,kl) )
520	! effective conductivity and 1st layer temperature (ln_cndflx=T)
521	CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
522	CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
523	! SIMIP diagnostics
524	CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d (1:npti), t_si (:,:,kl) )
525	CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) )
526	CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) )
527	! extensive variables
528	CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
529	CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
530	CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
531	CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
532	CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )
533	CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
534	! check convergence of heat diffusion scheme
535	IF( ln_zdf_chkcvg ) THEN
536	CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgerr_1d(1:npti), ztice_cvgerr(:,:,kl) )
537	CALL tab_1d_2d( npti, nptidx(1:npti), tice_cvgstp_1d(1:npti), ztice_cvgstp(:,:,kl) )
538	ENDIF
539	!
540	END SELECT
541	!
542	END SUBROUTINE ice_thd_1d2d
543
544
545	SUBROUTINE ice_thd_init
546	!!-------------------------------------------------------------------
547	!! * ROUTINE ice_thd_init *
548	!!
549	!! ** Purpose : Physical constants and parameters associated with
550	!! ice thermodynamics
551	!!
552	!! ** Method : Read the namthd namelist and check the parameters
553	!! called at the first timestep (nit000)
554	!!
555	!! ** input : Namelist namthd
556	!!-------------------------------------------------------------------
557	INTEGER :: ios ! Local integer output status for namelist read
558	!!
559	NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS, ln_leadhfx
560	!!-------------------------------------------------------------------
561	!
562	READ ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
563	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in reference namelist' )
564	READ ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
565	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd in configuration namelist' )
566	IF(lwm) WRITE( numoni, namthd )
567	!
568	IF(lwp) THEN ! control print
569	WRITE(numout,*)
570	WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
571	WRITE(numout,*) '~~~~~~~~~~~~'
572	WRITE(numout,*) ' Namelist namthd:'
573	WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH
574	WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA
575	WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO
576	WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS
577	WRITE(numout,*) ' heat in the leads is used to melt sea-ice before warming the ocean ln_leadhfx = ', ln_leadhfx
578	ENDIF
579	!
580	CALL ice_thd_zdf_init ! set ice heat diffusion parameters
581	IF( ln_icedA ) CALL ice_thd_da_init ! set ice lateral melting parameters
582	IF( ln_icedO ) CALL ice_thd_do_init ! set ice growth in open water parameters
583	CALL ice_thd_sal_init ! set ice salinity parameters
584	CALL ice_thd_pnd_init ! set melt ponds parameters
585	!
586	END SUBROUTINE ice_thd_init
587
588	#else
589	!!----------------------------------------------------------------------
590	!! Default option Dummy module NO SI3 sea-ice model
591	!!----------------------------------------------------------------------
592	#endif
593
594	!!======================================================================
595	END MODULE icethd

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

Download in other formats: