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.

limthd_2.F90 in branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90 @ 4859

Last change on this file since 4859 was 4859, checked in by smasson, 9 years ago
dev_4728_CNRS04_coupled_interface: cleaning of the coupling interface for OASIS3-MCT. 2
Property svn:keywords set to `Id`
File size: 34.7 KB

Line
1	MODULE limthd_2
2	!!======================================================================
3	!! * MODULE limthd_2 *
4	!! LIM thermo ice model : ice thermodynamic
5	!!======================================================================
6	!! History : 1.0 ! 2000-01 (LIM)
7	!! 2.0 ! 2002-07 (C. Ethe, G. Madec) F90
8	!! 2.0 ! 2003-08 (C. Ethe) add lim_thd_init
9	!! - ! 2008-2008 (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
10	!!---------------------------------------------------------------------
11	#if defined key_lim2
12	!!----------------------------------------------------------------------
13	!! 'key_lim2' : LIM 2.0 sea-ice model
14	!!----------------------------------------------------------------------
15	!! lim_thd_2 : thermodynamic of sea ice
16	!! lim_thd_init_2 : initialisation of sea-ice thermodynamic
17	!!----------------------------------------------------------------------
18	USE phycst ! physical constants
19	USE dom_oce ! ocean space and time domain variables
20	USE domvvl
21	USE lbclnk
22	USE in_out_manager ! I/O manager
23	USE lib_mpp
24	USE wrk_nemo ! work arrays
25	USE iom ! IOM library
26	USE ice_2 ! LIM sea-ice variables
27	USE sbc_oce !
28	USE sbc_ice !
29	USE thd_ice_2 ! LIM thermodynamic sea-ice variables
30	USE dom_ice_2 ! LIM sea-ice domain
31	USE limthd_zdf_2
32	USE limthd_lac_2
33	USE limtab_2
34	USE prtctl ! Print control
35	USE diaar5 , ONLY : lk_diaar5
36	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC lim_thd_2 ! called by lim_step
42
43	REAL(wp) :: epsi20 = 1.e-20 ! constant values
44	REAL(wp) :: epsi16 = 1.e-16 !
45	REAL(wp) :: epsi04 = 1.e-04 !
46	REAL(wp) :: rzero = 0.e0 !
47	REAL(wp) :: rone = 1.e0 !
48
49	!! * Substitutions
50	# include "domzgr_substitute.h90"
51	# include "vectopt_loop_substitute.h90"
52	!!-------- -------------------------------------------------------------
53	!! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
54	!! $Id$
55	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
56	!!----------------------------------------------------------------------
57	CONTAINS
58
59	SUBROUTINE lim_thd_2( kt )
60	!!-------------------------------------------------------------------
61	!! * ROUTINE lim_thd_2 *
62	!!
63	!! ** Purpose : This routine manages the ice thermodynamic.
64	!!
65	!! ** Action : - Initialisation of some variables
66	!! - Some preliminary computation (oceanic heat flux
67	!! at the ice base, snow acc.,heat budget of the leads)
68	!! - selection of the icy points and put them in an array
69	!! - call lim_vert_ther for vert ice thermodynamic
70	!! - back to the geographic grid
71	!! - selection of points for lateral accretion
72	!! - call lim_lat_acc for the ice accretion
73	!! - back to the geographic grid
74	!!
75	!! References : Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
76	!!---------------------------------------------------------------------
77	INTEGER, INTENT(in) :: kt ! number of iteration
78	!!
79	INTEGER :: ji, jj ! dummy loop indices
80	INTEGER :: nbpb ! nb of icy pts for thermo. cal.
81	INTEGER :: nbpac ! nb of pts for lateral accretion
82	CHARACTER (len=22) :: charout
83	REAL(wp) :: zfric_umin = 5e-03 ! lower bound for the friction velocity
84	REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity
85	REAL(wp) :: zinda ! switch for test. the val. of concen.
86	REAL(wp) :: zindb, zindg ! switches for test. the val of arg
87	REAL(wp) :: zfricp ! temporary scalar
88	REAL(wp) :: za , zh, zthsnice !
89	REAL(wp) :: zfric_u ! friction velocity
90	REAL(wp) :: zfntlat, zpareff ! test. the val. of lead heat budget
91
92	REAL(wp) :: zuice_m, zvice_m ! Sea-ice velocities at U & V-points
93	REAL(wp) :: zhice_u, zhice_v ! Sea-ice volume at U & V-points
94	REAL(wp) :: ztr_fram ! Sea-ice transport through Fram strait
95	REAL(wp) :: zrhoij, zrhoijm1 ! temporary scalars
96	REAL(wp) :: zztmp ! temporary scalars within a loop
97	REAL(wp), POINTER, DIMENSION(:,:) :: ztmp ! 2D workspace
98	REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif
99	REAL(wp), POINTER, DIMENSION(:,:) :: zlicegr ! link with lateral ice growth
100	!!$ REAL(wp), DIMENSION(:,:) :: firic ! IR flux over the ice (outputs only)
101	!!$ REAL(wp), DIMENSION(:,:) :: fcsic ! Sensible heat flux over the ice (outputs only)
102	!!$ REAL(wp), DIMENSION(:,:) :: fleic ! Latent heat flux over the ice (outputs only)
103	!!$ REAL(wp), DIMENSION(:,:) :: qlatic ! latent flux (outputs only)
104	REAL(wp), POINTER, DIMENSION(:,:) :: zdvosif ! Variation of volume at surface (outputs only)
105	REAL(wp), POINTER, DIMENSION(:,:) :: zdvobif ! Variation of ice volume at the bottom ice (outputs only)
106	REAL(wp), POINTER, DIMENSION(:,:) :: zdvolif ! Total variation of ice volume (outputs only)
107	REAL(wp), POINTER, DIMENSION(:,:) :: zdvonif ! Surface accretion Snow to Ice transformation (outputs only)
108	REAL(wp), POINTER, DIMENSION(:,:) :: zdvomif ! Bottom variation of ice volume due to melting (outputs only)
109	REAL(wp), POINTER, DIMENSION(:,:) :: zu_imasstr ! Sea-ice transport along i-axis at U-point (outputs only)
110	REAL(wp), POINTER, DIMENSION(:,:) :: zv_imasstr ! Sea-ice transport along j-axis at V-point (outputs only)
111	REAL(wp), POINTER, DIMENSION(:,:,:) :: zmsk ! 3D workspace
112	!!-------------------------------------------------------------------
113
114	CALL wrk_alloc( jpi, jpj, ztmp, zqlbsbq, zlicegr, zdvosif, zdvobif, zdvolif, zdvonif, zdvomif, zu_imasstr, zv_imasstr )
115	CALL wrk_alloc( jpi, jpj, jpk, zmsk )
116
117	IF( kt == nit000 ) CALL lim_thd_init_2 ! Initialization (first time-step only)
118
119	!-------------------------------------------!
120	! Initilization of diagnostic variables !
121	!-------------------------------------------!
122
123	!!gm needed? yes at least for some of these arrays
124	zdvosif(:,:) = 0.e0 ! variation of ice volume at surface
125	zdvobif(:,:) = 0.e0 ! variation of ice volume at bottom
126	zdvolif(:,:) = 0.e0 ! total variation of ice volume
127	zdvonif(:,:) = 0.e0 ! transformation of snow to sea-ice volume
128	zlicegr(:,:) = 0.e0 ! lateral variation of ice volume
129	zdvomif(:,:) = 0.e0 ! variation of ice volume at bottom due to melting only
130	ztr_fram = 0.e0 ! sea-ice transport through Fram strait
131	fstric (:,:) = 0.e0 ! part of solar radiation absorbing inside the ice
132	fscmbq (:,:) = 0.e0 ! linked with fstric
133	ffltbif(:,:) = 0.e0 ! linked with fstric
134	qfvbq (:,:) = 0.e0 ! linked with fstric
135	rdm_snw(:,:) = 0.e0 ! variation of snow mass over 1 time step
136	rdq_snw(:,:) = 0.e0 ! heat content associated with rdm_snw
137	rdm_ice(:,:) = 0.e0 ! variation of ice mass over 1 time step
138	rdq_ice(:,:) = 0.e0 ! heat content associated with rdm_ice
139	zmsk (:,:,:) = 0.e0
140	IF( ltrcdm2dc_ice ) fstric_daymean (:,:) = 0.e0 ! part of solar radiation absorbing inside the ice
141
142	! set to zero snow thickness smaller than epsi04
143	DO jj = 1, jpj
144	DO ji = 1, jpi
145	hsnif(ji,jj) = hsnif(ji,jj) * MAX( rzero, SIGN( rone , hsnif(ji,jj) - epsi04 ) )
146	END DO
147	END DO
148	!!gm better coded (do not use SIGN...)
149	! WHERE( hsnif(:,:) < epsi04 ) hsnif(:,:) = 0.e0
150	!!gm
151
152	IF(ln_ctl) CALL prt_ctl( tab2d_1=hsnif, clinfo1=' lim_thd: hsnif : ' )
153
154	!-----------------------------------!
155	! Treatment of particular cases !
156	!-----------------------------------!
157
158	DO jj = 1, jpj
159	DO ji = 1, jpi
160	! snow is transformed into ice if the original ice cover disappears.
161	zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -hicif(ji,jj) ) )
162	hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0
163	hsnif(ji,jj) = ( rone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
164	dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass
165
166	! the lead fraction, frld, must be little than or equal to amax (ice ridging).
167	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
168	zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) )
169	za = zindb * MIN( rone, ( 1.0 - frld(ji,jj) ) * uscomi )
170	hsnif (ji,jj) = hsnif(ji,jj) * za
171	hicif (ji,jj) = hicif(ji,jj) * za
172	qstoif(ji,jj) = qstoif(ji,jj) * za
173	frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za, epsi20 )
174
175	! the in situ ice thickness, hicif, must be equal to or greater than hiclim.
176	zh = MAX( rone , zindb * hiclim / MAX( hicif(ji,jj), epsi20 ) )
177	hsnif (ji,jj) = hsnif(ji,jj) * zh
178	hicif (ji,jj) = hicif(ji,jj) * zh
179	qstoif(ji,jj) = qstoif(ji,jj) * zh
180	frld (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh
181	END DO
182	END DO
183
184	IF(ln_ctl) THEN
185	CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ' )
186	CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ' )
187	CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ' )
188	CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ' )
189	CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ' )
190	ENDIF
191
192
193	!-------------------------------!
194	! Thermodynamics of sea ice !
195	!-------------------------------!
196
197	! Partial computation of forcing for the thermodynamic sea ice model.
198	!--------------------------------------------------------------------------
199
200	!CDIR NOVERRCHK
201	DO jj = 1, jpj
202	!CDIR NOVERRCHK
203	DO ji = 1, jpi
204	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
205	zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) )
206	pfrld(ji,jj) = frld(ji,jj)
207	zfricp = 1.0 - frld(ji,jj)
208	zinda = 1.0 - MAX( rzero , SIGN( rone , - zfricp ) )
209
210	! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
211	thcm(ji,jj) = 0.e0
212
213	! net downward heat flux from the ice to the ocean, expressed as a function of ocean
214	! temperature and turbulent mixing (McPhee, 1992)
215	zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity
216	fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) )
217	qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
218
219	! partial computation of the lead energy budget (qldif)
220	#if defined key_coupled
221	qldif(ji,jj) = tms(ji,jj) * rdt_ice &
222	& * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) ) &
223	& + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp ) &
224	& + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) ) )
225	#else
226	qldif(ji,jj) = tms(ji,jj) * rdt_ice * frld(ji,jj) &
227	& * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) &
228	& + qns(ji,jj) + fdtcn(ji,jj) &
229	& + ( 1.0 - zindb ) * fsbbq(ji,jj) )
230	#endif
231	! parlat : percentage of energy used for lateral ablation (0.0)
232	zfntlat = 1.0 - MAX( rzero , SIGN( rone , - qldif(ji,jj) ) )
233	zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat
234	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 )
235	qldif (ji,jj) = zpareff * qldif(ji,jj)
236	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
237
238	! energy needed to bring ocean surface layer until its freezing
239	qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
240
241	! calculate oceanic heat flux.
242	fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) )
243
244	! computation of the thermodynamic ice production (only needed for output)
245	hicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) )
246	END DO
247	END DO
248
249	! Select icy points and fulfill arrays for the vectorial grid.
250	!----------------------------------------------------------------------
251	nbpb = 0
252	DO jj = 1, jpj
253	DO ji = 1, jpi
254	IF ( frld(ji,jj) < 1.0 ) THEN
255	nbpb = nbpb + 1
256	npb(nbpb) = (jj - 1) * jpi + ji
257	ENDIF
258	END DO
259	END DO
260
261	IF(ln_ctl) THEN
262	CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld : ', tab2d_2=thcm , clinfo2=' thcm : ')
263	CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ')
264	CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif : ', tab2d_2=zqlbsbq, clinfo2=' zqlbsbq : ')
265	CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif : ', tab2d_2=fbif , clinfo2=' fbif : ')
266	zmsk(:,:,1) = tms(:,:)
267	CALL prt_ctl(tab2d_1=qcmif , clinfo1=' lim_thd: qcmif : ', mask1=zmsk)
268	CALL prt_ctl(tab2d_1=hicifp, clinfo1=' lim_thd: hicifp : ')
269	WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb
270	CALL prt_ctl_info(charout)
271	ENDIF
272
273
274	! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation
275	!------------------------------------------------------------------------------------
276
277	IF( nbpb > 0 ) THEN
278	!
279	! put the variable in a 1-D array for thermodynamics process
280	CALL tab_2d_1d_2( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) )
281	CALL tab_2d_1d_2( nbpb, h_ice_1d (1:nbpb) , hicif , jpi, jpj, npb(1:nbpb) )
282	CALL tab_2d_1d_2( nbpb, h_snow_1d (1:nbpb) , hsnif , jpi, jpj, npb(1:nbpb) )
283	CALL tab_2d_1d_2( nbpb, sist_1d (1:nbpb) , sist , jpi, jpj, npb(1:nbpb) )
284	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 1 ), tbif(:,:,1) , jpi, jpj, npb(1:nbpb) )
285	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 2 ), tbif(:,:,2) , jpi, jpj, npb(1:nbpb) )
286	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 3 ), tbif(:,:,3) , jpi, jpj, npb(1:nbpb) )
287	CALL tab_2d_1d_2( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice(:,:,1) , jpi, jpj, npb(1:nbpb) )
288	IF( ltrcdm2dc_ice ) &
289	& CALL tab_2d_1d_2( nbpb, qsr_ice_mean_1d (1:nbpb), qsr_ice_mean(:,:,1), jpi, jpj, npb(1:nbpb) )
290	CALL tab_2d_1d_2( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) )
291	CALL tab_2d_1d_2( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) )
292	CALL tab_2d_1d_2( nbpb, qns_ice_1d(1:nbpb) , qns_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
293	CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
294	IF( .NOT. lk_cpl ) THEN
295	CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
296	CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
297	ENDIF
298	CALL tab_2d_1d_2( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) )
299	CALL tab_2d_1d_2( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) )
300	CALL tab_2d_1d_2( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) )
301	CALL tab_2d_1d_2( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) )
302	CALL tab_2d_1d_2( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) )
303	CALL tab_2d_1d_2( nbpb, qstbif_1d (1:nbpb) , qstoif , jpi, jpj, npb(1:nbpb) )
304	CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb) , rdm_ice , jpi, jpj, npb(1:nbpb) )
305	CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb) , rdq_ice , jpi, jpj, npb(1:nbpb) )
306	CALL tab_2d_1d_2( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) )
307	CALL tab_2d_1d_2( nbpb, rdm_snw_1d (1:nbpb) , rdm_snw , jpi, jpj, npb(1:nbpb) )
308	CALL tab_2d_1d_2( nbpb, rdq_snw_1d (1:nbpb) , rdq_snw , jpi, jpj, npb(1:nbpb) )
309	CALL tab_2d_1d_2( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) )
310	!
311	CALL lim_thd_zdf_2( 1, nbpb ) ! compute ice growth
312	!
313	! back to the geographic grid.
314	CALL tab_1d_2d_2( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj )
315	CALL tab_1d_2d_2( nbpb, hicif , npb, h_ice_1d (1:nbpb) , jpi, jpj )
316	CALL tab_1d_2d_2( nbpb, hsnif , npb, h_snow_1d (1:nbpb) , jpi, jpj )
317	CALL tab_1d_2d_2( nbpb, sist , npb, sist_1d (1:nbpb) , jpi, jpj )
318	CALL tab_1d_2d_2( nbpb, tbif(:,:,1), npb, tbif_1d (1:nbpb , 1 ), jpi, jpj )
319	CALL tab_1d_2d_2( nbpb, tbif(:,:,2), npb, tbif_1d (1:nbpb , 2 ), jpi, jpj )
320	CALL tab_1d_2d_2( nbpb, tbif(:,:,3), npb, tbif_1d (1:nbpb , 3 ), jpi, jpj )
321	CALL tab_1d_2d_2( nbpb, fscmbq , npb, fscbq_1d (1:nbpb) , jpi, jpj )
322	CALL tab_1d_2d_2( nbpb, ffltbif , npb, fltbif_1d (1:nbpb) , jpi, jpj )
323	CALL tab_1d_2d_2( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
324	CALL tab_1d_2d_2( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
325	CALL tab_1d_2d_2( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
326	CALL tab_1d_2d_2( nbpb, qstoif , npb, qstbif_1d (1:nbpb) , jpi, jpj )
327	CALL tab_1d_2d_2( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj )
328	CALL tab_1d_2d_2( nbpb, rdq_ice , npb, rdq_ice_1d(1:nbpb) , jpi, jpj )
329	CALL tab_1d_2d_2( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
330	CALL tab_1d_2d_2( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj )
331	CALL tab_1d_2d_2( nbpb, rdq_snw , npb, rdq_snw_1d(1:nbpb) , jpi, jpj )
332	CALL tab_1d_2d_2( nbpb, zdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
333	CALL tab_1d_2d_2( nbpb, zdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
334	CALL tab_1d_2d_2( nbpb, zdvomif , npb, rdvomif_1d(1:nbpb) , jpi, jpj )
335	CALL tab_1d_2d_2( nbpb, zdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
336	CALL tab_1d_2d_2( nbpb, zdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
337	CALL tab_1d_2d_2( nbpb, qsr_ice(:,:,1), npb, qsr_ice_1d(1:nbpb) , jpi, jpj )
338	CALL tab_1d_2d_2( nbpb, qns_ice(:,:,1), npb, qns_ice_1d(1:nbpb) , jpi, jpj )
339	IF( ltrcdm2dc_ice )THEN
340	CALL tab_1d_2d_2( nbpb, fstric_daymean , npb, fstbif_daymean_1d (1:nbpb) , jpi, jpj )
341	CALL tab_1d_2d_2( nbpb, qsr_ice_mean(:,:,1), npb, qsr_ice_mean_1d(1:nbpb) , jpi, jpj )
342	ENDIF
343	IF( .NOT. lk_cpl ) CALL tab_1d_2d_2( nbpb, qla_ice(:,:,1), npb, qla_ice_1d(1:nbpb) , jpi, jpj )
344	!
345	ENDIF
346
347	! Up-date sea ice thickness
348	!--------------------------
349	DO jj = 1, jpj
350	DO ji = 1, jpi
351	phicif(ji,jj) = hicif(ji,jj)
352	hicif(ji,jj) = hicif(ji,jj) * ( rone - MAX( rzero, SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) )
353	END DO
354	END DO
355
356
357	! Tricky trick : add 2 to frld in the Southern Hemisphere
358	!--------------------------------------------------------
359	IF( fcor(1,1) < 0.e0 ) THEN
360	DO jj = 1, njeqm1
361	DO ji = 1, jpi
362	frld(ji,jj) = frld(ji,jj) + 2.0
363	END DO
364	END DO
365	ENDIF
366
367	CALL lbc_lnk( frld , 'T', 1. )
368
369	! Select points for lateral accretion (this occurs when heat exchange
370	! between ice and ocean is negative; ocean losing heat)
371	!-----------------------------------------------------------------
372	nbpac = 0
373	DO jj = 1, jpj
374	DO ji = 1, jpi
375	!i yes! IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
376	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
377	nbpac = nbpac + 1
378	npac( nbpac ) = (jj - 1) * jpi + ji
379	ENDIF
380	END DO
381	END DO
382
383	IF(ln_ctl) THEN
384	CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=hicif, clinfo2=' hicif : ')
385	WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac
386	CALL prt_ctl_info(charout)
387	ENDIF
388
389
390	! If ocean gains heat do nothing ; otherwise, one performs lateral accretion
391	!--------------------------------------------------------------------------------
392	IF( nbpac > 0 ) THEN
393	!
394	zlicegr(:,:) = rdm_ice(:,:) ! to output the lateral sea-ice growth
395	!...Put the variable in a 1-D array for lateral accretion
396	CALL tab_2d_1d_2( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) )
397	CALL tab_2d_1d_2( nbpac, h_snow_1d (1:nbpac) , hsnif , jpi, jpj, npac(1:nbpac) )
398	CALL tab_2d_1d_2( nbpac, h_ice_1d (1:nbpac) , hicif , jpi, jpj, npac(1:nbpac) )
399	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) )
400	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) )
401	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) )
402	CALL tab_2d_1d_2( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
403	CALL tab_2d_1d_2( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
404	CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac) , qstoif , jpi, jpj, npac(1:nbpac) )
405	CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice , jpi, jpj, npac(1:nbpac) )
406	CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac) , rdq_ice , jpi, jpj, npac(1:nbpac) )
407	CALL tab_2d_1d_2( nbpac, dvlbq_1d (1:nbpac) , zdvolif , jpi, jpj, npac(1:nbpac) )
408	CALL tab_2d_1d_2( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) )
409	!
410	CALL lim_thd_lac_2( 1 , nbpac ) ! lateral accretion routine.
411	!
412	! back to the geographic grid
413	CALL tab_1d_2d_2( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj )
414	CALL tab_1d_2d_2( nbpac, hsnif , npac(1:nbpac), h_snow_1d (1:nbpac) , jpi, jpj )
415	CALL tab_1d_2d_2( nbpac, hicif , npac(1:nbpac), h_ice_1d (1:nbpac) , jpi, jpj )
416	CALL tab_1d_2d_2( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d (1:nbpac , 1 ), jpi, jpj )
417	CALL tab_1d_2d_2( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d (1:nbpac , 2 ), jpi, jpj )
418	CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d (1:nbpac , 3 ), jpi, jpj )
419	CALL tab_1d_2d_2( nbpac, qstoif , npac(1:nbpac), qstbif_1d (1:nbpac) , jpi, jpj )
420	CALL tab_1d_2d_2( nbpac, rdm_ice , npac(1:nbpac), rdm_ice_1d(1:nbpac) , jpi, jpj )
421	CALL tab_1d_2d_2( nbpac, rdq_ice , npac(1:nbpac), rdq_ice_1d(1:nbpac) , jpi, jpj )
422	CALL tab_1d_2d_2( nbpac, zdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj )
423	!
424	ENDIF
425
426
427	! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
428	! Update daily thermodynamic ice production.
429	!------------------------------------------------------------------------------
430	DO jj = 1, jpj
431	DO ji = 1, jpi
432	frld (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) )
433	fr_i (ji,jj) = 1.0 - frld(ji,jj)
434	hicifp(ji,jj) = hicif(ji,jj) * fr_i(ji,jj) - hicifp(ji,jj)
435	END DO
436	END DO
437
438	! Outputs
439	!--------------------------------------------------------------------------------
440	ztmp(:,:) = 1. - pfrld(:,:) ! fraction of ice after the dynamic, before the thermodynamic
441	CALL iom_put( 'ist_cea', (sist(:,:) - rt0) * ztmp(:,:) ) ! Ice surface temperature [Celius]
442	CALL iom_put( 'qsr_ai_cea', qsr_ice(:,:,1) * ztmp(:,:) ) ! Solar flux over the ice [W/m2]
443	CALL iom_put( 'qns_ai_cea', qns_ice(:,:,1) * ztmp(:,:) ) ! Non-solar flux over the ice [W/m2]
444	IF( .NOT. lk_cpl ) CALL iom_put( 'qla_ai_cea', qla_ice(:,:,1) * ztmp(:,:) ) ! Latent flux over the ice [W/m2]
445	!
446	CALL iom_put( 'snowthic_cea', hsnif (:,:) * fr_i(:,:) ) ! Snow thickness [m]
447	CALL iom_put( 'icethic_cea' , hicif (:,:) * fr_i(:,:) ) ! Ice thickness [m]
448	zztmp = 1.0 / rdt_ice
449	CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp ) ! Ice produced [m/s]
450	IF( lk_diaar5 ) THEN
451	CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp ) ! Snow melt [kg/m2/s]
452	zztmp = rhoic / rdt_ice
453	CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp ) ! Snow to Ice transformation [kg/m2/s]
454	CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp ) ! Melt at Sea Ice top [kg/m2/s]
455	CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp ) ! Melt at Sea Ice bottom [kg/m2/s]
456	zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
457	CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp ) ! Lateral sea ice growth [kg/m2/s]
458	ENDIF
459	!
460	! Compute the Eastward & Northward sea-ice transport
461	zztmp = 0.25 * rhoic
462	DO jj = 1, jpjm1
463	DO ji = 1, jpim1 ! NO vector opt.
464	! Ice velocities, volume & transport at U & V-points
465	zuice_m = u_ice(ji+1,jj+1) + u_ice(ji+1,jj )
466	zvice_m = v_ice(ji+1,jj+1) + v_ice(ji ,jj+1)
467	zhice_u = hicif(ji,jj)e2t(ji,jj)fr_i(ji,jj) + hicif(ji+1,jj )e2t(ji+1,jj )fr_i(ji+1,jj )
468	zhice_v = hicif(ji,jj)e1t(ji,jj)fr_i(ji,jj) + hicif(ji ,jj+1)e1t(ji ,jj+1)fr_i(ji ,jj+1)
469	zu_imasstr(ji,jj) = zztmp * zhice_u * zuice_m
470	zv_imasstr(ji,jj) = zztmp * zhice_v * zvice_m
471	END DO
472	END DO
473	CALL lbc_lnk( zu_imasstr, 'U', -1. ) ; CALL lbc_lnk( zv_imasstr, 'V', -1. )
474	CALL iom_put( 'u_imasstr', zu_imasstr(:,:) ) ! Ice transport along i-axis at U-point [kg/s]
475	CALL iom_put( 'v_imasstr', zv_imasstr(:,:) ) ! Ice transport along j-axis at V-point [kg/s]
476
477	!! Fram Strait sea-ice transport (sea-ice + snow) (in ORCA2 = 5 points)
478	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
479	DO jj = mj0(137), mj1(137) ! B grid
480	IF( mj0(jj-1) >= nldj ) THEN
481	DO ji = MAX(mi0(134),nldi), MIN(mi1(138),nlei)
482	zrhoij = e1t(ji,jj ) * fr_i(ji,jj ) * ( rhoichicif(ji,jj ) + rhosnhsnif(ji,jj ) )
483	zrhoijm1 = e1t(ji,jj-1) * fr_i(ji,jj-1) * ( rhoichicif(ji,jj-1) + rhosnhsnif(ji,jj-1) )
484	ztr_fram = ztr_fram - 0.25 * ( v_ice(ji,jj)+ v_ice(ji+1,jj) ) * ( zrhoij + zrhoijm1 )
485	END DO
486	ENDIF
487	END DO
488	IF( lk_mpp ) CALL mpp_sum( ztr_fram )
489	CALL iom_put( 'fram_trans', ztr_fram ) ! Ice transport through Fram strait [kg/s]
490	ENDIF
491
492	!! ce ztmp(:,:) = 1. - AINT( frld(:,:), wp ) ! return 1 as soon as there is ice
493	!! ce A big warning because the model crashes on IDRIS/IBM SP6 with xlf 13.1.0.3, see ticket #761
494	!! ce We Unroll the loop and everything works fine
495	DO jj = 1, jpj
496	DO ji = 1, jpi
497	ztmp(ji,jj) = 1. - AINT( frld(ji,jj), wp ) ! return 1 as soon as there is ice
498	END DO
499	END DO
500	!
501	CALL iom_put( 'ice_pres' , ztmp ) ! Ice presence [-]
502	CALL iom_put( 'ist_ipa' , ( sist(:,:) - rt0 ) * ztmp(:,:) ) ! Ice surface temperature [Celius]
503	CALL iom_put( 'uice_ipa' , u_ice(:,:) * ztmp(:,:) ) ! Ice velocity along i-axis at I-point [m/s]
504	CALL iom_put( 'vice_ipa' , v_ice(:,:) * ztmp(:,:) ) ! Ice velocity along j-axis at I-point [m/s]
505
506	IF(ln_ctl) THEN
507	CALL prt_ctl_info(' lim_thd end ')
508	CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ' )
509	CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ' )
510	CALL prt_ctl( tab2d_1=phicif , clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ' )
511	CALL prt_ctl( tab2d_1=sist , clinfo1=' lim_thd: sist : ' )
512	CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ' )
513	CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ' )
514	CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ' )
515	CALL prt_ctl( tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ' )
516	CALL prt_ctl( tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ' )
517	ENDIF
518	!
519	CALL wrk_dealloc( jpi, jpj, ztmp, zqlbsbq, zlicegr, zdvosif, zdvobif, zdvolif, zdvonif, zdvomif, zu_imasstr, zv_imasstr )
520	CALL wrk_dealloc( jpi, jpj, jpk, zmsk )
521	!
522	END SUBROUTINE lim_thd_2
523
524
525	SUBROUTINE lim_thd_init_2
526	!!-------------------------------------------------------------------
527	!! * ROUTINE lim_thd_init_2 *
528	!!
529	!! ** Purpose : Physical constants and parameters linked to the ice
530	!! thermodynamics
531	!!
532	!! ** Method : Read the namicethd namelist and check the ice-thermo
533	!! parameter values called at the first timestep (nit000)
534	!!
535	!! ** input : Namelist namicether
536	!!-------------------------------------------------------------------
537	INTEGER :: ios ! Local integer output status for namelist read
538	NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , &
539	& swiqst, sbeta , parlat, hakspl, hibspl, exld, &
540	& hakdif, hnzst , thth , parsub, alphs
541	!!-------------------------------------------------------------------
542
543	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
544	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
545	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
546
547	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
548	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
549	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
550	IF(lwm) WRITE ( numoni, namicethd )
551
552	IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
553	!
554	IF(lwp) THEN ! control print
555	WRITE(numout,*)
556	WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation '
557	WRITE(numout,*)'~~~~~~~~~~~~~~'
558	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
559	WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit
560	WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin
561	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
562	WRITE(numout,*)' maximum lead fraction amax = ', amax
563	WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst
564	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
565	WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
566	WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat
567	WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl
568	WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl
569	WRITE(numout,*)' exponent for leads-closure rate exld = ', exld
570	WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif
571	WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth
572	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
573	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
574	WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs
575	ENDIF
576	!
577	uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction
578	rcdsn = hakdif * rcdsn
579	rcdic = hakdif * rcdic
580	!
581	IF( hsndif > 100.e0 .OR. hicdif > 100.e0 ) THEN
582	cnscg = 0.e0
583	ELSE
584	cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic
585	ENDIF
586	!
587	END SUBROUTINE lim_thd_init_2
588
589	#else
590	!!----------------------------------------------------------------------
591	!! Default option Dummy module NO LIM 2.0 sea-ice model
592	!!----------------------------------------------------------------------
593	CONTAINS
594	SUBROUTINE lim_thd_2 ! Dummy routine
595	END SUBROUTINE lim_thd_2
596	#endif
597
598	!!======================================================================
599	END MODULE limthd_2

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

Download in other formats: