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_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90 @ 4846

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

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

Download in other formats: