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_lac.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 26.0 KB

Line
1	MODULE limthd_lac
2	!!======================================================================
3	!! * MODULE limthd_lac *
4	!! lateral thermodynamic growth of the ice
5	!!======================================================================
6	!! History : LIM ! 2005-12 (M. Vancoppenolle) Original code
7	!! - ! 2006-01 (M. Vancoppenolle) add ITD
8	!! 3.0 ! 2007-07 (M. Vancoppenolle) Mass and energy conservation tested
9	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
10	!!----------------------------------------------------------------------
11	#if defined key_lim3
12	!!----------------------------------------------------------------------
13	!! 'key_lim3' LIM3 sea-ice model
14	!!----------------------------------------------------------------------
15	!! lim_lat_acr : lateral accretion of ice
16	!!----------------------------------------------------------------------
17	USE par_oce ! ocean parameters
18	USE dom_oce ! domain variables
19	USE phycst ! physical constants
20	USE sbc_oce ! Surface boundary condition: ocean fields
21	USE sbc_ice ! Surface boundary condition: ice fields
22	USE thd_ice ! LIM thermodynamics
23	USE ice ! LIM variables
24	USE limtab ! LIM 2D <==> 1D
25	USE limcons ! LIM conservation
26	USE in_out_manager ! I/O manager
27	USE lib_mpp ! MPP library
28	USE wrk_nemo ! work arrays
29	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
30	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
31	USE limthd_ent
32	USE limvar
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC lim_thd_lac ! called by lim_thd
38
39	!!----------------------------------------------------------------------
40	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
41	!! $Id$
42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
43	!!----------------------------------------------------------------------
44	CONTAINS
45
46	SUBROUTINE lim_thd_lac
47	!!-------------------------------------------------------------------
48	!! * ROUTINE lim_thd_lac *
49	!!
50	!! ** Purpose : Computation of the evolution of the ice thickness and
51	!! concentration as a function of the heat balance in the leads.
52	!! It is only used for lateral accretion
53	!!
54	!! ** Method : Ice is formed in the open water when ocean lose heat
55	!! (heat budget of open water Bl is negative) .
56	!! Computation of the increase of 1-A (ice concentration) fol-
57	!! lowing the law :
58	!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
59	!! where - h0 is the thickness of ice created in the lead
60	!! - a is a minimum fraction for leads
61	!! - F is a monotonic non-increasing function defined as:
62	!! F(X)=( 1 - Xexld )(1.0/exld)
63	!! - exld is the exponent closure rate (=2 default val.)
64	!!
65	!! ** Action : - Adjustment of snow and ice thicknesses and heat
66	!! content in brine pockets
67	!! - Updating ice internal temperature
68	!! - Computation of variation of ice volume and mass
69	!! - Computation of frldb after lateral accretion and
70	!! update ht_s_1d, ht_i_1d and tbif_1d(:,:)
71	!!------------------------------------------------------------------------
72	INTEGER :: ji,jj,jk,jl ! dummy loop indices
73	INTEGER :: nbpac ! local integers
74	INTEGER :: ii, ij, iter ! - -
75	REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars
76	REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - -
77	REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - -
78	CHARACTER (len = 15) :: fieldid
79
80	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
81	REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg)
82	REAL(wp) :: zEw ! seawater specific enthalpy (J/kg)
83	REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean)
84
85	REAL(wp) :: zv_newfra
86
87	INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows
88	REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not
89
90	REAL(wp), POINTER, DIMENSION(:) :: zv_newice ! volume of accreted ice
91	REAL(wp), POINTER, DIMENSION(:) :: za_newice ! fractional area of accreted ice
92	REAL(wp), POINTER, DIMENSION(:) :: zh_newice ! thickness of accreted ice
93	REAL(wp), POINTER, DIMENSION(:) :: ze_newice ! heat content of accreted ice
94	REAL(wp), POINTER, DIMENSION(:) :: zs_newice ! salinity of accreted ice
95	REAL(wp), POINTER, DIMENSION(:) :: zo_newice ! age of accreted ice
96	REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget
97	REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget
98	REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction
99	REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom
100	REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector)
101
102	REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl
103	REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl
104	REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i
105	REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i
106	REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i
107
108	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i
109
110	REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity
111
112	REAL(wp) :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used)
113	!!-----------------------------------------------------------------------!
114
115	CALL wrk_alloc( jpij, jcat ) ! integer
116	CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
117	CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
118	CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
119	CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d )
120	CALL wrk_alloc( jpi,jpj, zvrel )
121
122	CALL lim_var_agg(1)
123	CALL lim_var_glo2eqv
124	!------------------------------------------------------------------------------\|
125	! 2) Convert units for ice internal energy
126	!------------------------------------------------------------------------------\|
127	!$OMP PARALLEL
128	DO jl = 1, jpl
129	DO jk = 1, nlay_i
130	!$OMP DO schedule(static) private(jj,ji,rswitch)
131	DO jj = 1, jpj
132	DO ji = 1, jpi
133	!Energy of melting q(S,T) [J.m-3]
134	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice
135	e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )
136	END DO
137	END DO
138	END DO
139	END DO
140
141	!------------------------------------------------------------------------------!
142	! 3) Collection thickness of ice formed in leads and polynyas
143	!------------------------------------------------------------------------------!
144	! hicol is the thickness of new ice formed in open water
145	! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
146	! Frazil ice forms in open water, is transported by wind
147	! accumulates at the edge of the consolidated ice edge
148	! where it forms aggregates of a specific thickness called
149	! collection thickness.
150
151	! Note : the following algorithm currently breaks vectorization
152	!
153
154	!$OMP DO schedule(static) private(jj,ji)
155	DO jj = 1, jpj
156	DO ji = 1, jpi
157	zvrel(ji,jj) = 0._wp
158	END DO
159	END DO
160
161	!$OMP DO schedule(static) private(jj,ji)
162	DO jj = 1, jpj
163	DO ji = 1, jpi
164	! Default new ice thickness
165	IF( qlead(ji,jj) < 0._wp ) THEN ; hicol(ji,jj) = rn_hnewice
166	ELSE ; hicol(ji,jj) = 0._wp
167	END IF
168	END DO
169	END DO
170	!$OMP END PARALLEL
171
172	IF( ln_frazil ) THEN
173
174	!--------------------
175	! Physical constants
176	!--------------------
177
178	zhicrit = 0.04 ! frazil ice thickness
179	ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
180	zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
181	zgamafr = 0.03
182
183	!$OMP PARALLEL
184	!$OMP DO schedule(static) private(jj,ji)
185	DO jj = 1, jpj
186	DO ji = 1, jpi
187	hicol(ji,jj) = 0._wp
188	END DO
189	END DO
190
191	!$OMP DO schedule(static) private(jj,ji,ztaux,ztauy,ztenagm,rswitch,zvfrx,zvfry,zvgx,zvgy,zvrel2,iter,zf,zfp)
192	DO jj = 2, jpjm1
193	DO ji = 2, jpim1
194	IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
195	!-------------
196	! Wind stress
197	!-------------
198	! C-grid wind stress components
199	ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) &
200	& + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp
201	ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) &
202	& + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp
203	! Square root of wind stress
204	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
205
206	!---------------------
207	! Frazil ice velocity
208	!---------------------
209	rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
210	zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
211	zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
212
213	!-------------------
214	! Pack ice velocity
215	!-------------------
216	! C-grid ice velocity
217	zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
218	zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
219
220	!-----------------------------------
221	! Relative frazil/pack ice velocity
222	!-----------------------------------
223	! absolute relative velocity
224	rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
225	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) &
226	& + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
227	zvrel(ji,jj) = SQRT( zvrel2 )
228
229	!---------------------
230	! Iterative procedure
231	!---------------------
232	hicol(ji,jj) = zhicrit + ( zhicrit + 0.1 ) &
233	& / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2
234
235	iter = 1
236	DO WHILE ( iter < 20 )
237	zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) - &
238	& hicol(ji,jj) * zhicrit * ztwogp * zvrel2
239	zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
240
241	hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 )
242	iter = iter + 1
243	END DO
244
245	ENDIF ! end of selection of pixels where ice forms
246
247	END DO
248	END DO
249	!$OMP END PARALLEL
250	!
251	CALL lbc_lnk( zvrel, 'T', 1. )
252	CALL lbc_lnk( hicol, 'T', 1. )
253
254	ENDIF ! End of computation of frazil ice collection thickness
255
256	!------------------------------------------------------------------------------!
257	! 4) Identify grid points where new ice forms
258	!------------------------------------------------------------------------------!
259
260	!-------------------------------------
261	! Select points for new ice formation
262	!-------------------------------------
263	! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
264	nbpac = 0
265	npac(:) = 0
266	!
267	DO jj = 1, jpj
268	DO ji = 1, jpi
269	IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
270	nbpac = nbpac + 1
271	npac( nbpac ) = (jj - 1) * jpi + ji
272	ENDIF
273	END DO
274	END DO
275
276	! debug point to follow
277	jiindex_1d = 0
278	IF( ln_limctl ) THEN
279	DO ji = mi0(iiceprt), mi1(iiceprt)
280	DO jj = mj0(jiceprt), mj1(jiceprt)
281	IF ( qlead(ji,jj) < 0._wp ) THEN
282	jiindex_1d = (jj - 1) * jpi + ji
283	ENDIF
284	END DO
285	END DO
286	ENDIF
287
288	IF( ln_limctl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
289
290	!------------------------------
291	! Move from 2-D to 1-D vectors
292	!------------------------------
293	! If ocean gains heat do nothing. Otherwise compute new ice formation
294
295	IF ( nbpac > 0 ) THEN
296
297	CALL tab_2d_1d( nbpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )
298	DO jl = 1, jpl
299	CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
300	CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
301	CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
302	DO jk = 1, nlay_i
303	CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
304	END DO
305	END DO
306
307	CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) )
308	CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )
309	CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw , jpi, jpj, npac(1:nbpac) )
310	CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw , jpi, jpj, npac(1:nbpac) )
311	CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) )
312	CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )
313
314	CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd , jpi, jpj, npac(1:nbpac) )
315	CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw , jpi, jpj, npac(1:nbpac) )
316	CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac) , rn_amax_2d, jpi, jpj, npac(1:nbpac) )
317
318	!------------------------------------------------------------------------------!
319	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
320	!------------------------------------------------------------------------------!
321
322	!-----------------------------------------
323	! Keep old ice areas and volume in memory
324	!-----------------------------------------
325	zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:)
326	za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
327
328	!----------------------
329	! Thickness of new ice
330	!----------------------
331	zh_newice(1:nbpac) = hicol_1d(1:nbpac)
332
333	!----------------------
334	! Salinity of new ice
335	!----------------------
336	SELECT CASE ( nn_icesal )
337	CASE ( 1 ) ! Sice = constant
338	zs_newice(1:nbpac) = rn_icesal
339	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
340	DO ji = 1, nbpac
341	ii = MOD( npac(ji) - 1 , jpi ) + 1
342	ij = ( npac(ji) - 1 ) / jpi + 1
343	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) )
344	END DO
345	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
346	zs_newice(1:nbpac) = 2.3
347	END SELECT
348
349	!-------------------------
350	! Heat content of new ice
351	!-------------------------
352	! We assume that new ice is formed at the seawater freezing point
353	DO ji = 1, nbpac
354	ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K)
355	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) &
356	& + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) &
357	& - rcp * ( ztmelts - rt0 ) )
358	END DO
359
360	!----------------
361	! Age of new ice
362	!----------------
363	DO ji = 1, nbpac
364	zo_newice(ji) = 0._wp
365	END DO
366
367	!-------------------
368	! Volume of new ice
369	!-------------------
370	DO ji = 1, nbpac
371
372	zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg]
373
374	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg]
375	! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
376
377	zdE = zEi - zEw ! specific enthalpy difference [J/kg]
378
379	zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0)
380	! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point
381	zv_newice(ji) = - zfmdt * r1_rhoic
382
383	zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux
384
385	! Contribution to heat flux to the ocean [W.m-2], >0
386	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
387	! Total heat flux used in this process [W.m-2]
388	hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
389	! mass flux
390	wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
391	! salt flux
392	sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
393	END DO
394
395	zv_frazb(:) = 0._wp
396	IF( ln_frazil ) THEN
397	! A fraction zfrazb of frazil ice is accreted at the ice bottom
398	DO ji = 1, nbpac
399	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
400	zfrazb = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
401	zv_frazb(ji) = zfrazb * zv_newice(ji)
402	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
403	END DO
404	END IF
405
406	!-----------------
407	! Area of new ice
408	!-----------------
409	DO ji = 1, nbpac
410	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
411	END DO
412
413	!------------------------------------------------------------------------------!
414	! 6) Redistribute new ice area and volume into ice categories !
415	!------------------------------------------------------------------------------!
416
417	!------------------------
418	! 6.1) lateral ice growth
419	!------------------------
420	! If lateral ice growth gives an ice concentration gt 1, then
421	! we keep the excessive volume in memory and attribute it later to bottom accretion
422	DO ji = 1, nbpac
423	IF ( za_newice(ji) > ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
424	zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
425	zdv_res(ji) = zda_res (ji) * zh_newice(ji)
426	za_newice(ji) = za_newice(ji) - zda_res (ji)
427	zv_newice(ji) = zv_newice(ji) - zdv_res (ji)
428	ELSE
429	zda_res(ji) = 0._wp
430	zdv_res(ji) = 0._wp
431	ENDIF
432	END DO
433
434	! find which category to fill
435	zat_i_1d(:) = 0._wp
436	DO jl = 1, jpl
437	DO ji = 1, nbpac
438	IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
439	za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
440	zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
441	jcat (ji) = jl
442	ENDIF
443	zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl)
444	END DO
445	END DO
446
447	! Heat content
448	DO ji = 1, nbpac
449	jl = jcat(ji) ! categroy in which new ice is put
450	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice
451	END DO
452
453	DO jk = 1, nlay_i
454	!$OMP PARALLEL DO schedule(static) private(ji,jl,rswitch)
455	DO ji = 1, nbpac
456	jl = jcat(ji)
457	rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
458	ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &
459	& ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) &
460	& * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
461	END DO
462	END DO
463
464	!------------------------------------------------
465	! 6.2) bottom ice growth + ice enthalpy remapping
466	!------------------------------------------------
467	DO jl = 1, jpl
468
469	! for remapping
470	h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
471	qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
472	DO jk = 1, nlay_i
473	!$OMP PARALLEL DO schedule(static) private(ji)
474	DO ji = 1, nbpac
475	h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
476	qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
477	END DO
478	END DO
479
480	! new volumes including lateral/bottom accretion + residual
481	!$OMP PARALLEL DO schedule(static) private(ji,rswitch,zv_newfra)
482	DO ji = 1, nbpac
483	rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
484	zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
485	za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)
486	zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
487	! for remapping
488	h_i_old (ji,nlay_i+1) = zv_newfra
489	qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
490	ENDDO
491	! --- Ice enthalpy remapping --- !
492	CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) )
493	ENDDO
494
495	!-----------------
496	! Update salinity
497	!-----------------
498	DO jl = 1, jpl
499	!$OMP PARALLEL DO schedule(static) private(ji,zdv)
500	DO ji = 1, nbpac
501	zdv = zv_i_1d(ji,jl) - zv_b(ji,jl)
502	zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
503	END DO
504	END DO
505
506	!------------------------------------------------------------------------------!
507	! 7) Change 2D vectors to 1D vectors
508	!------------------------------------------------------------------------------!
509	DO jl = 1, jpl
510	CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
511	CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
512	CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
513	DO jk = 1, nlay_i
514	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj )
515	END DO
516	END DO
517	CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
518	CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
519
520	CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
521	CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj )
522	!
523	ENDIF ! nbpac > 0
524
525	!------------------------------------------------------------------------------!
526	! 8) Change units for e_i
527	!------------------------------------------------------------------------------!
528	DO jl = 1, jpl
529	DO jk = 1, nlay_i
530	!$OMP PARALLEL DO schedule(static) private(jj,ji)
531	DO jj = 1, jpj
532	DO ji = 1, jpi
533	! heat content in J/m2
534	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
535	END DO
536	END DO
537	END DO
538	END DO
539
540	!
541	CALL wrk_dealloc( jpij, jcat ) ! integer
542	CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
543	CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
544	CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
545	CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d )
546	CALL wrk_dealloc( jpi,jpj, zvrel )
547	!
548	END SUBROUTINE lim_thd_lac
549
550	#else
551	!!----------------------------------------------------------------------
552	!! Default option NO LIM3 sea-ice model
553	!!----------------------------------------------------------------------
554	CONTAINS
555	SUBROUTINE lim_thd_lac ! Empty routine
556	END SUBROUTINE lim_thd_lac
557	#endif
558
559	!!======================================================================
560	END MODULE limthd_lac

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

Download in other formats: