1 | MODULE limthd |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE limthd *** |
---|
4 | !! LIM thermo ice model : ice thermodynamic |
---|
5 | !!====================================================================== |
---|
6 | #if defined key_ice_lim |
---|
7 | !!---------------------------------------------------------------------- |
---|
8 | !! 'key_ice_lim' : LIM sea-ice model |
---|
9 | !!---------------------------------------------------------------------- |
---|
10 | !! lim_thd : thermodynamic of sea ice |
---|
11 | !! lim_thd_init : initialisation of sea-ice thermodynamic |
---|
12 | !!---------------------------------------------------------------------- |
---|
13 | !! * Modules used |
---|
14 | USE phycst ! physical constants |
---|
15 | USE dom_oce ! ocean space and time domain variables |
---|
16 | USE lbclnk |
---|
17 | USE in_out_manager ! I/O manager |
---|
18 | USE ice ! LIM sea-ice variables |
---|
19 | USE ice_oce ! sea-ice/ocean variables |
---|
20 | USE flx_oce ! sea-ice/ocean forcings variables |
---|
21 | USE thd_ice ! LIM thermodynamic sea-ice variables |
---|
22 | USE dom_ice ! LIM sea-ice domain |
---|
23 | USE iceini |
---|
24 | USE limthd_zdf |
---|
25 | USE limthd_lac |
---|
26 | USE limtab |
---|
27 | USE prtctl ! Print control |
---|
28 | |
---|
29 | IMPLICIT NONE |
---|
30 | PRIVATE |
---|
31 | |
---|
32 | !! * Routine accessibility |
---|
33 | PUBLIC lim_thd ! called by lim_step |
---|
34 | |
---|
35 | !! * Module variables |
---|
36 | REAL(wp) :: & ! constant values |
---|
37 | epsi20 = 1.e-20 , & |
---|
38 | epsi16 = 1.e-16 , & |
---|
39 | epsi04 = 1.e-04 , & |
---|
40 | zzero = 0.e0 , & |
---|
41 | zone = 1.e0 |
---|
42 | |
---|
43 | !! * Substitutions |
---|
44 | # include "domzgr_substitute.h90" |
---|
45 | # include "vectopt_loop_substitute.h90" |
---|
46 | !!-------- ------------------------------------------------------------- |
---|
47 | !! LIM 2.0, UCL-LOCEAN-IPSL (2005) |
---|
48 | !! $Header$ |
---|
49 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
---|
50 | !!---------------------------------------------------------------------- |
---|
51 | |
---|
52 | CONTAINS |
---|
53 | |
---|
54 | SUBROUTINE lim_thd( kt ) |
---|
55 | !!------------------------------------------------------------------- |
---|
56 | !! *** ROUTINE lim_thd *** |
---|
57 | !! |
---|
58 | !! ** Purpose : This routine manages the ice thermodynamic. |
---|
59 | !! |
---|
60 | !! ** Action : - Initialisation of some variables |
---|
61 | !! - Some preliminary computation (oceanic heat flux |
---|
62 | !! at the ice base, snow acc.,heat budget of the leads) |
---|
63 | !! - selection of the icy points and put them in an array |
---|
64 | !! - call lim_vert_ther for vert ice thermodynamic |
---|
65 | !! - back to the geographic grid |
---|
66 | !! - selection of points for lateral accretion |
---|
67 | !! - call lim_lat_acc for the ice accretion |
---|
68 | !! - back to the geographic grid |
---|
69 | !! |
---|
70 | !! ** References : |
---|
71 | !! H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 |
---|
72 | !! |
---|
73 | !! History : |
---|
74 | !! 1.0 ! 00-01 (LIM) |
---|
75 | !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 |
---|
76 | !!--------------------------------------------------------------------- |
---|
77 | INTEGER, INTENT(in) :: kt ! number of iteration |
---|
78 | |
---|
79 | INTEGER :: ji, jj, & ! dummy loop indices |
---|
80 | nbpb , & ! nb of icy pts for thermo. cal. |
---|
81 | nbpac ! nb of pts for lateral accretion |
---|
82 | CHARACTER (len=22) :: charout |
---|
83 | REAL(wp) :: & |
---|
84 | zfric_umin = 5e-03 , & ! lower bound for the friction velocity |
---|
85 | zfric_umax = 2e-02 ! upper bound for the friction velocity |
---|
86 | REAL(wp) :: & |
---|
87 | zinda , & ! switch for test. the val. of concen. |
---|
88 | zindb, zindg , & ! switches for test. the val of arg |
---|
89 | za , zh, zthsnice , & |
---|
90 | zfric_u , & ! friction velocity |
---|
91 | zfnsol , & ! total non solar heat |
---|
92 | zfontn , & ! heat flux from snow thickness |
---|
93 | zfntlat, zpareff ! test. the val. of lead heat budget |
---|
94 | REAL(wp), DIMENSION(jpi,jpj) :: & |
---|
95 | zhicifp , & ! ice thickness for outputs |
---|
96 | zqlbsbq ! link with lead energy budget qldif |
---|
97 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: & |
---|
98 | zmsk ! working array |
---|
99 | !!------------------------------------------------------------------- |
---|
100 | |
---|
101 | IF( kt == nit000 ) CALL lim_thd_init ! Initialization (first time-step only) |
---|
102 | |
---|
103 | !-------------------------------------------! |
---|
104 | ! Initilization of diagnostic variables ! |
---|
105 | !-------------------------------------------! |
---|
106 | |
---|
107 | !i est-ce utile? oui au moins en partie |
---|
108 | rdvosif(:,:) = 0.e0 ! variation of ice volume at surface |
---|
109 | rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom |
---|
110 | fdvolif(:,:) = 0.e0 ! total variation of ice volume |
---|
111 | rdvonif(:,:) = 0.e0 ! lateral variation of ice volume |
---|
112 | fstric (:,:) = 0.e0 ! part of solar radiation absorbing inside the ice |
---|
113 | fscmbq (:,:) = 0.e0 ! linked with fstric |
---|
114 | ffltbif(:,:) = 0.e0 ! linked with fstric |
---|
115 | qfvbq (:,:) = 0.e0 ! linked with fstric |
---|
116 | rdmsnif(:,:) = 0.e0 ! variation of snow mass per unit area |
---|
117 | rdmicif(:,:) = 0.e0 ! variation of ice mass per unit area |
---|
118 | hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. |
---|
119 | zmsk (:,:,:) = 0.e0 |
---|
120 | |
---|
121 | DO jj = 1, jpj |
---|
122 | DO ji = 1, jpi |
---|
123 | hsnif(ji,jj) = hsnif(ji,jj) * MAX( zzero, SIGN( zone , hsnif(ji,jj) - epsi04 ) ) |
---|
124 | END DO |
---|
125 | END DO |
---|
126 | |
---|
127 | IF(ln_ctl) CALL prt_ctl(tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ') |
---|
128 | |
---|
129 | !-----------------------------------! |
---|
130 | ! Treatment of particular cases ! |
---|
131 | !-----------------------------------! |
---|
132 | |
---|
133 | DO jj = 1, jpj |
---|
134 | DO ji = 1, jpi |
---|
135 | ! snow is transformed into ice if the original ice cover disappears. |
---|
136 | zindg = tms(ji,jj) * MAX( zzero , SIGN( zone , -hicif(ji,jj) ) ) |
---|
137 | hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0 |
---|
138 | hsnif(ji,jj) = ( zone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn |
---|
139 | dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass |
---|
140 | |
---|
141 | ! the lead fraction, frld, must be little than or equal to amax (ice ridging). |
---|
142 | zthsnice = hsnif(ji,jj) + hicif(ji,jj) |
---|
143 | zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) |
---|
144 | za = zindb * MIN( zone, ( 1.0 - frld(ji,jj) ) * uscomi ) |
---|
145 | hsnif (ji,jj) = hsnif(ji,jj) * za |
---|
146 | hicif (ji,jj) = hicif(ji,jj) * za |
---|
147 | qstoif(ji,jj) = qstoif(ji,jj) * za |
---|
148 | frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za , epsi20 ) |
---|
149 | |
---|
150 | ! the in situ ice thickness, hicif, must be equal to or greater than hiclim. |
---|
151 | zh = MAX( zone , zindb * hiclim / MAX( hicif(ji,jj) , epsi20 ) ) |
---|
152 | hsnif (ji,jj) = hsnif(ji,jj) * zh |
---|
153 | hicif (ji,jj) = hicif(ji,jj) * zh |
---|
154 | qstoif(ji,jj) = qstoif(ji,jj) * zh |
---|
155 | frld (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh |
---|
156 | END DO |
---|
157 | END DO |
---|
158 | |
---|
159 | IF(ln_ctl) THEN |
---|
160 | CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif : ') |
---|
161 | CALL prt_ctl(tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ') |
---|
162 | CALL prt_ctl(tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ') |
---|
163 | CALL prt_ctl(tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ') |
---|
164 | CALL prt_ctl(tab2d_1=frld , clinfo1=' lim_thd: frld : ') |
---|
165 | ENDIF |
---|
166 | |
---|
167 | |
---|
168 | !-------------------------------! |
---|
169 | ! Thermodynamics of sea ice ! |
---|
170 | !-------------------------------! |
---|
171 | |
---|
172 | ! Partial computation of forcing for the thermodynamic sea ice model. |
---|
173 | !-------------------------------------------------------------------------- |
---|
174 | |
---|
175 | !CDIR NOVERRCHK |
---|
176 | DO jj = 1, jpj |
---|
177 | !CDIR NOVERRCHK |
---|
178 | DO ji = 1, jpi |
---|
179 | zthsnice = hsnif(ji,jj) + hicif(ji,jj) |
---|
180 | zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) |
---|
181 | pfrld(ji,jj) = frld(ji,jj) |
---|
182 | zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) |
---|
183 | |
---|
184 | ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget |
---|
185 | thcm(ji,jj) = 0.e0 |
---|
186 | |
---|
187 | ! net downward heat flux from the ice to the ocean, expressed as a function of ocean |
---|
188 | ! temperature and turbulent mixing (McPhee, 1992) |
---|
189 | zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity |
---|
190 | fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) ) |
---|
191 | qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice |
---|
192 | |
---|
193 | ! partial computation of the lead energy budget (qldif) |
---|
194 | zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting |
---|
195 | zfnsol = qnsr_oce(ji,jj) ! total non solar flux |
---|
196 | qldif(ji,jj) = tms(ji,jj) * ( qsr_oce(ji,jj) * ( 1.0 - thcm(ji,jj) ) & |
---|
197 | & + zfnsol + fdtcn(ji,jj) - zfontn & |
---|
198 | & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & |
---|
199 | & * frld(ji,jj) * rdt_ice |
---|
200 | ! parlat : percentage of energy used for lateral ablation (0.0) |
---|
201 | zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) ) |
---|
202 | zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat |
---|
203 | zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 ) |
---|
204 | qldif (ji,jj) = zpareff * qldif(ji,jj) |
---|
205 | qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) |
---|
206 | |
---|
207 | ! energy needed to bring ocean surface layer until its freezing |
---|
208 | qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda ) |
---|
209 | |
---|
210 | ! calculate oceanic heat flux. |
---|
211 | fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) ) |
---|
212 | |
---|
213 | ! computation of the daily thermodynamic ice production (only needed for output) |
---|
214 | zhicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) |
---|
215 | END DO |
---|
216 | END DO |
---|
217 | |
---|
218 | |
---|
219 | ! Select icy points and fulfill arrays for the vectorial grid. |
---|
220 | !---------------------------------------------------------------------- |
---|
221 | nbpb = 0 |
---|
222 | DO jj = 1, jpj |
---|
223 | DO ji = 1, jpi |
---|
224 | IF ( frld(ji,jj) < 1.0 ) THEN |
---|
225 | nbpb = nbpb + 1 |
---|
226 | npb(nbpb) = (jj - 1) * jpi + ji |
---|
227 | ENDIF |
---|
228 | END DO |
---|
229 | END DO |
---|
230 | |
---|
231 | IF(ln_ctl) THEN |
---|
232 | CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld : ', tab2d_2=thcm , clinfo2=' thcm : ') |
---|
233 | CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ') |
---|
234 | CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif : ', tab2d_2=zqlbsbq, clinfo2=' zqlbsbq : ') |
---|
235 | CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif : ', tab2d_2=fbif , clinfo2=' fbif : ') |
---|
236 | zmsk(:,:,1) = tms(:,:) |
---|
237 | CALL prt_ctl(tab2d_1=qcmif , clinfo1=' lim_thd: qcmif : ', mask1=zmsk) |
---|
238 | CALL prt_ctl(tab2d_1=zhicifp, clinfo1=' lim_thd: zhicifp : ') |
---|
239 | WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb |
---|
240 | CALL prt_ctl_info(charout) |
---|
241 | ENDIF |
---|
242 | |
---|
243 | |
---|
244 | ! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation |
---|
245 | !------------------------------------------------------------------------------------ |
---|
246 | |
---|
247 | IF ( nbpb > 0) THEN |
---|
248 | |
---|
249 | ! put the variable in a 1-D array for thermodynamics process |
---|
250 | CALL tab_2d_1d( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) ) |
---|
251 | CALL tab_2d_1d( nbpb, h_ice_1d (1:nbpb) , hicif , jpi, jpj, npb(1:nbpb) ) |
---|
252 | CALL tab_2d_1d( nbpb, h_snow_1d (1:nbpb) , hsnif , jpi, jpj, npb(1:nbpb) ) |
---|
253 | CALL tab_2d_1d( nbpb, sist_1d (1:nbpb) , sist , jpi, jpj, npb(1:nbpb) ) |
---|
254 | CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 1 ), tbif(:,:,1), jpi, jpj, npb(1:nbpb) ) |
---|
255 | CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 2 ), tbif(:,:,2), jpi, jpj, npb(1:nbpb) ) |
---|
256 | CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 3 ), tbif(:,:,3), jpi, jpj, npb(1:nbpb) ) |
---|
257 | CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice , jpi, jpj, npb(1:nbpb) ) |
---|
258 | CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) ) |
---|
259 | CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) ) |
---|
260 | CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb) , qnsr_ice , jpi, jpj, npb(1:nbpb) ) |
---|
261 | #if ! defined key_coupled |
---|
262 | CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb) , qla_ice , jpi, jpj, npb(1:nbpb) ) |
---|
263 | CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice , jpi, jpj, npb(1:nbpb) ) |
---|
264 | #endif |
---|
265 | CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice , jpi, jpj, npb(1:nbpb) ) |
---|
266 | CALL tab_2d_1d( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) ) |
---|
267 | CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) ) |
---|
268 | CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) ) |
---|
269 | CALL tab_2d_1d( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) ) |
---|
270 | CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) ) |
---|
271 | CALL tab_2d_1d( nbpb, qstbif_1d (1:nbpb) , qstoif , jpi, jpj, npb(1:nbpb) ) |
---|
272 | CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb) , rdmicif , jpi, jpj, npb(1:nbpb) ) |
---|
273 | CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) ) |
---|
274 | CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) ) |
---|
275 | |
---|
276 | CALL lim_thd_zdf( 1, nbpb ) ! compute ice growth |
---|
277 | |
---|
278 | ! back to the geographic grid. |
---|
279 | CALL tab_1d_2d( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj ) |
---|
280 | CALL tab_1d_2d( nbpb, hicif , npb, h_ice_1d (1:nbpb) , jpi, jpj ) |
---|
281 | CALL tab_1d_2d( nbpb, hsnif , npb, h_snow_1d (1:nbpb) , jpi, jpj ) |
---|
282 | CALL tab_1d_2d( nbpb, sist , npb, sist_1d (1:nbpb) , jpi, jpj ) |
---|
283 | CALL tab_1d_2d( nbpb, tbif(:,:,1), npb, tbif_1d (1:nbpb , 1 ), jpi, jpj ) |
---|
284 | CALL tab_1d_2d( nbpb, tbif(:,:,2), npb, tbif_1d (1:nbpb , 2 ), jpi, jpj ) |
---|
285 | CALL tab_1d_2d( nbpb, tbif(:,:,3), npb, tbif_1d (1:nbpb , 3 ), jpi, jpj ) |
---|
286 | CALL tab_1d_2d( nbpb, fscmbq , npb, fscbq_1d (1:nbpb) , jpi, jpj ) |
---|
287 | CALL tab_1d_2d( nbpb, ffltbif , npb, fltbif_1d (1:nbpb) , jpi, jpj ) |
---|
288 | CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj ) |
---|
289 | CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj ) |
---|
290 | CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj ) |
---|
291 | CALL tab_1d_2d( nbpb, qstoif , npb, qstbif_1d (1:nbpb) , jpi, jpj ) |
---|
292 | CALL tab_1d_2d( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb) , jpi, jpj ) |
---|
293 | CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj ) |
---|
294 | CALL tab_1d_2d( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb) , jpi, jpj ) |
---|
295 | CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj ) |
---|
296 | CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj ) |
---|
297 | CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj ) |
---|
298 | CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj ) |
---|
299 | |
---|
300 | |
---|
301 | ENDIF |
---|
302 | |
---|
303 | |
---|
304 | ! Up-date sea ice thickness. |
---|
305 | !--------------------------------- |
---|
306 | DO jj = 1, jpj |
---|
307 | DO ji = 1, jpi |
---|
308 | phicif(ji,jj) = hicif(ji,jj) |
---|
309 | hicif(ji,jj) = hicif(ji,jj) * ( zone - MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) ) |
---|
310 | END DO |
---|
311 | END DO |
---|
312 | |
---|
313 | |
---|
314 | ! Tricky trick : add 2 to frld in the Southern Hemisphere. |
---|
315 | !---------------------------------------------------------- |
---|
316 | IF( fcor(1,1) < 0.e0 ) THEN |
---|
317 | DO jj = 1, njeqm1 |
---|
318 | DO ji = 1, jpi |
---|
319 | frld(ji,jj) = frld(ji,jj) + 2.0 |
---|
320 | END DO |
---|
321 | END DO |
---|
322 | ENDIF |
---|
323 | |
---|
324 | |
---|
325 | ! Select points for lateral accretion (this occurs when heat exchange |
---|
326 | ! between ice and ocean is negative; ocean losing heat) |
---|
327 | !----------------------------------------------------------------- |
---|
328 | nbpac = 0 |
---|
329 | DO jj = 1, jpj |
---|
330 | DO ji = 1, jpi |
---|
331 | !i yes! IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN |
---|
332 | IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN |
---|
333 | nbpac = nbpac + 1 |
---|
334 | npac( nbpac ) = (jj - 1) * jpi + ji |
---|
335 | ENDIF |
---|
336 | END DO |
---|
337 | END DO |
---|
338 | |
---|
339 | IF(ln_ctl) THEN |
---|
340 | CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=hicif, clinfo2=' hicif : ') |
---|
341 | WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac |
---|
342 | CALL prt_ctl_info(charout) |
---|
343 | ENDIF |
---|
344 | |
---|
345 | |
---|
346 | ! |
---|
347 | ! If ocean gains heat do nothing ; otherwise, one performs lateral accretion |
---|
348 | !-------------------------------------------------------------------------------- |
---|
349 | |
---|
350 | IF( nbpac > 0 ) THEN |
---|
351 | |
---|
352 | !...Put the variable in a 1-D array for lateral accretion |
---|
353 | CALL tab_2d_1d( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) ) |
---|
354 | CALL tab_2d_1d( nbpac, h_snow_1d (1:nbpac) , hsnif , jpi, jpj, npac(1:nbpac) ) |
---|
355 | CALL tab_2d_1d( nbpac, h_ice_1d (1:nbpac) , hicif , jpi, jpj, npac(1:nbpac) ) |
---|
356 | CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) ) |
---|
357 | CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) ) |
---|
358 | CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) ) |
---|
359 | CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) |
---|
360 | CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) |
---|
361 | CALL tab_2d_1d( nbpac, qstbif_1d (1:nbpac) , qstoif , jpi, jpj, npac(1:nbpac) ) |
---|
362 | CALL tab_2d_1d( nbpac, rdmicif_1d(1:nbpac) , rdmicif , jpi, jpj, npac(1:nbpac) ) |
---|
363 | CALL tab_2d_1d( nbpac, dvlbq_1d (1:nbpac) , fdvolif , jpi, jpj, npac(1:nbpac) ) |
---|
364 | CALL tab_2d_1d( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) ) |
---|
365 | |
---|
366 | ! call lateral accretion routine. |
---|
367 | CALL lim_thd_lac( 1 , nbpac ) |
---|
368 | |
---|
369 | ! back to the geographic grid |
---|
370 | CALL tab_1d_2d( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj ) |
---|
371 | CALL tab_1d_2d( nbpac, hsnif , npac(1:nbpac), h_snow_1d (1:nbpac) , jpi, jpj ) |
---|
372 | CALL tab_1d_2d( nbpac, hicif , npac(1:nbpac), h_ice_1d (1:nbpac) , jpi, jpj ) |
---|
373 | CALL tab_1d_2d( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d (1:nbpac , 1 ), jpi, jpj ) |
---|
374 | CALL tab_1d_2d( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d (1:nbpac , 2 ), jpi, jpj ) |
---|
375 | CALL tab_1d_2d( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d (1:nbpac , 3 ), jpi, jpj ) |
---|
376 | CALL tab_1d_2d( nbpac, qstoif , npac(1:nbpac), qstbif_1d (1:nbpac) , jpi, jpj ) |
---|
377 | CALL tab_1d_2d( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d(1:nbpac) , jpi, jpj ) |
---|
378 | CALL tab_1d_2d( nbpac, fdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj ) |
---|
379 | |
---|
380 | ENDIF |
---|
381 | |
---|
382 | |
---|
383 | ! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick) |
---|
384 | ! Update daily thermodynamic ice production. |
---|
385 | !------------------------------------------------------------------------------ |
---|
386 | |
---|
387 | DO jj = 1, jpj |
---|
388 | DO ji = 1, jpi |
---|
389 | frld (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) ) |
---|
390 | hicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) - zhicifp(ji,jj) + hicifp(ji,jj) |
---|
391 | END DO |
---|
392 | END DO |
---|
393 | |
---|
394 | IF(ln_ctl) THEN |
---|
395 | CALL prt_ctl_info(' lim_thd end ') |
---|
396 | CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ') |
---|
397 | CALL prt_ctl(tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ') |
---|
398 | CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ') |
---|
399 | CALL prt_ctl(tab2d_1=sist , clinfo1=' lim_thd: sist : ') |
---|
400 | CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ') |
---|
401 | CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ') |
---|
402 | CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ') |
---|
403 | CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ') |
---|
404 | CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ') |
---|
405 | ENDIF |
---|
406 | |
---|
407 | END SUBROUTINE lim_thd |
---|
408 | |
---|
409 | |
---|
410 | SUBROUTINE lim_thd_init |
---|
411 | !!------------------------------------------------------------------- |
---|
412 | !! *** ROUTINE lim_thd_init *** |
---|
413 | !! |
---|
414 | !! ** Purpose : Physical constants and parameters linked to the ice |
---|
415 | !! thermodynamics |
---|
416 | !! |
---|
417 | !! ** Method : Read the namicethd namelist and check the ice-thermo |
---|
418 | !! parameter values called at the first timestep (nit000) |
---|
419 | !! |
---|
420 | !! ** input : Namelist namicether |
---|
421 | !! |
---|
422 | !! history : |
---|
423 | !! 8.5 ! 03-08 (C. Ethe) original code |
---|
424 | !!------------------------------------------------------------------- |
---|
425 | NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , & |
---|
426 | & swiqst, sbeta , parlat, hakspl, hibspl, exld, & |
---|
427 | & hakdif, hnzst , thth , parsub, alphs |
---|
428 | !!------------------------------------------------------------------- |
---|
429 | |
---|
430 | |
---|
431 | ! Define the initial parameters |
---|
432 | ! ------------------------- |
---|
433 | REWIND( numnam_ice ) |
---|
434 | READ ( numnam_ice , namicethd ) |
---|
435 | IF(lwp) THEN |
---|
436 | WRITE(numout,*) |
---|
437 | WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation ' |
---|
438 | WRITE(numout,*)'~~~~~~~~~~~~' |
---|
439 | WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt |
---|
440 | WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit |
---|
441 | WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin |
---|
442 | WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim |
---|
443 | WRITE(numout,*)' maximum lead fraction amax = ', amax |
---|
444 | WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst |
---|
445 | WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' |
---|
446 | WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta |
---|
447 | WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat |
---|
448 | WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl |
---|
449 | WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl |
---|
450 | WRITE(numout,*)' exponent for leads-closure rate exld = ', exld |
---|
451 | WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif |
---|
452 | WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth |
---|
453 | WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst |
---|
454 | WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub |
---|
455 | WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs |
---|
456 | ENDIF |
---|
457 | |
---|
458 | uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction |
---|
459 | rcdsn = hakdif * rcdsn |
---|
460 | rcdic = hakdif * rcdic |
---|
461 | |
---|
462 | IF ( ( hsndif > 100.e0 ) .OR. ( hicdif > 100.e0 ) ) THEN |
---|
463 | cnscg = 0.e0 |
---|
464 | ELSE |
---|
465 | cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic |
---|
466 | ENDIF |
---|
467 | |
---|
468 | END SUBROUTINE lim_thd_init |
---|
469 | |
---|
470 | #else |
---|
471 | !!---------------------------------------------------------------------- |
---|
472 | !! Default option Dummy module NO LIM sea-ice model |
---|
473 | !!---------------------------------------------------------------------- |
---|
474 | CONTAINS |
---|
475 | SUBROUTINE lim_thd ! Dummy routine |
---|
476 | END SUBROUTINE lim_thd |
---|
477 | #endif |
---|
478 | |
---|
479 | !!====================================================================== |
---|
480 | END MODULE limthd |
---|