Context Navigation

icethd.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

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

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

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icethd.F90 @ 12340

Download in other formats: