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.

limvar.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7702

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

Line
1	MODULE limvar
2	!!======================================================================
3	!! * MODULE limvar *
4	!! Different sets of ice model variables
5	!! how to switch from one to another
6	!!
7	!! There are three sets of variables
8	!! VGLO : global variables of the model
9	!! - v_i (jpi,jpj,jpl)
10	!! - v_s (jpi,jpj,jpl)
11	!! - a_i (jpi,jpj,jpl)
12	!! - t_s (jpi,jpj,jpl)
13	!! - e_i (jpi,jpj,nlay_i,jpl)
14	!! - smv_i(jpi,jpj,jpl)
15	!! - oa_i (jpi,jpj,jpl)
16	!! VEQV : equivalent variables sometimes used in the model
17	!! - ht_i(jpi,jpj,jpl)
18	!! - ht_s(jpi,jpj,jpl)
19	!! - t_i (jpi,jpj,nlay_i,jpl)
20	!! ...
21	!! VAGG : aggregate variables, averaged/summed over all
22	!! thickness categories
23	!! - vt_i(jpi,jpj)
24	!! - vt_s(jpi,jpj)
25	!! - at_i(jpi,jpj)
26	!! - et_s(jpi,jpj) !total snow heat content
27	!! - et_i(jpi,jpj) !total ice thermal content
28	!! - smt_i(jpi,jpj) !mean ice salinity
29	!! - tm_i (jpi,jpj) !mean ice temperature
30	!!======================================================================
31	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
32	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
33	!!----------------------------------------------------------------------
34	#if defined key_lim3
35	!!----------------------------------------------------------------------
36	!! 'key_lim3' LIM3 sea-ice model
37	!!----------------------------------------------------------------------
38	USE par_oce ! ocean parameters
39	USE phycst ! physical constants (ocean directory)
40	USE sbc_oce ! Surface boundary condition: ocean fields
41	USE ice ! ice variables
42	USE thd_ice ! ice variables (thermodynamics)
43	USE in_out_manager ! I/O manager
44	USE lib_mpp ! MPP library
45	USE wrk_nemo ! work arrays
46	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
47
48	IMPLICIT NONE
49	PRIVATE
50
51	PUBLIC lim_var_agg
52	PUBLIC lim_var_glo2eqv
53	PUBLIC lim_var_eqv2glo
54	PUBLIC lim_var_salprof
55	PUBLIC lim_var_bv
56	PUBLIC lim_var_salprof1d
57	PUBLIC lim_var_zapsmall
58	PUBLIC lim_var_itd
59
60	!!----------------------------------------------------------------------
61	!! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
62	!! $Id$
63	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
64	!!----------------------------------------------------------------------
65	CONTAINS
66
67	SUBROUTINE lim_var_agg( kn )
68	!!------------------------------------------------------------------
69	!! * ROUTINE lim_var_agg *
70	!!
71	!! ** Purpose : aggregates ice-thickness-category variables to all-ice variables
72	!! i.e. it turns VGLO into VAGG
73	!! ** Method :
74	!!
75	!! ** Arguments : n = 1, at_i vt_i only
76	!! n = 2 everything
77	!!
78	!! note : you could add an argument when you need only at_i, vt_i
79	!! and when you need everything
80	!!------------------------------------------------------------------
81	INTEGER, INTENT( in ) :: kn ! =1 at_i & vt only ; = what is needed
82	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_s, ze_i
83	!
84	INTEGER :: ji, jj, jk, jl ! dummy loop indices
85	!!------------------------------------------------------------------
86
87	CALL wrk_alloc( jpi, jpj, nlay_s, ze_s )
88	CALL wrk_alloc( jpi, jpj, nlay_i, ze_i )
89	! integrated values
90	!$OMP PARALLEL
91	!$OMP DO schedule(static) private(jj, ji)
92	DO jj = 1, jpj
93	DO ji = 1, jpi
94	vt_i (ji,jj) = 0._wp
95	vt_s (ji,jj) = 0._wp
96	at_i (ji,jj) = 0._wp
97	et_s(ji,jj) = 0._wp
98	et_i(ji,jj) = 0._wp
99	END DO
100	END DO
101	DO jl = 1, jpl
102	!$OMP DO schedule(static) private(jj, ji)
103	DO jj = 1, jpj
104	DO ji = 1, jpi
105	vt_i (ji,jj) = vt_i (ji,jj) + v_i (ji,jj,jl)
106	vt_s (ji,jj) = vt_s (ji,jj) + v_s (ji,jj,jl)
107	at_i (ji,jj) = at_i (ji,jj) + a_i (ji,jj,jl)
108	END DO
109	END DO
110	END DO
111	DO jk = 1, nlay_s
112	!$OMP DO schedule(static) private(jj, ji)
113	DO jj = 1, jpj
114	DO ji = 1, jpi
115	ze_s(ji,jj,jk) = 0._wp
116	END DO
117	END DO
118	END DO
119	DO jk = 1, nlay_i
120	!$OMP DO schedule(static) private(jj, ji)
121	DO jj = 1, jpj
122	DO ji = 1, jpi
123	ze_i(ji,jj,jk) = 0._wp
124	END DO
125	END DO
126	END DO
127	DO jl = 1, jpl
128	DO jk = 1, nlay_s
129	!$OMP DO schedule(static) private(jj, ji)
130	DO jj = 1, jpj
131	DO ji = 1, jpi
132	ze_s(ji,jj,jk) = ze_s(ji,jj,jk) + e_s(ji,jj,jk,jl)
133	END DO
134	END DO
135	END DO
136	END DO
137	DO jl = 1, jpl
138	DO jk = 1, nlay_i
139	!$OMP DO schedule(static) private(jj, ji)
140	DO jj = 1, jpj
141	DO ji = 1, jpi
142	ze_i(ji,jj,jk) = ze_i(ji,jj,jk) + e_i(ji,jj,jk,jl)
143	END DO
144	END DO
145	END DO
146	END DO
147	DO jk = 1, nlay_s
148	!$OMP DO schedule(static) private(jj, ji)
149	DO jj = 1, jpj
150	DO ji = 1, jpi
151	et_s(ji,jj) = et_s(ji,jj) + ze_s(ji,jj,jk)
152	END DO
153	END DO
154	END DO
155	DO jk = 1, nlay_i
156	!$OMP DO schedule(static) private(jj, ji)
157	DO jj = 1, jpj
158	DO ji = 1, jpi
159	et_i(ji,jj) = et_i(ji,jj) + ze_i(ji,jj,jk)
160	END DO
161	END DO
162	END DO
163
164	! open water fraction
165	!$OMP DO schedule(static) private(jj, ji)
166	DO jj = 1, jpj
167	DO ji = 1, jpi
168	ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )
169	END DO
170	END DO
171	!$OMP END PARALLEL
172
173	IF( kn > 1 ) THEN
174
175	!$OMP PARALLEL
176	! mean ice/snow thickness
177	!$OMP DO schedule(static) private(jj,ji,rswitch)
178	DO jj = 1, jpj
179	DO ji = 1, jpi
180	rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
181	htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
182	htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
183	ENDDO
184	ENDDO
185
186	! mean temperature (K), salinity and age
187	!$OMP DO schedule(static) private(jj,ji)
188	DO jj = 1, jpj
189	DO ji = 1, jpi
190	smt_i(ji,jj) = 0._wp
191	tm_i(ji,jj) = 0._wp
192	tm_su(ji,jj) = 0._wp
193	om_i (ji,jj) = 0._wp
194	ENDDO
195	ENDDO
196	DO jl = 1, jpl
197
198	!$OMP DO schedule(static) private(jj,ji,rswitch)
199	DO jj = 1, jpj
200	DO ji = 1, jpi
201	rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
202	tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
203	om_i (ji,jj) = om_i (ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
204	END DO
205	END DO
206
207	DO jk = 1, nlay_i
208	!$OMP DO schedule(static) private(jj,ji,rswitch)
209	DO jj = 1, jpj
210	DO ji = 1, jpi
211	rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
212	tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) &
213	& / MAX( vt_i(ji,jj) , epsi10 )
214	smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &
215	& / MAX( vt_i(ji,jj) , epsi10 )
216	END DO
217	END DO
218	END DO
219	END DO
220	!$OMP END PARALLEL
221	tm_i = tm_i + rt0
222	tm_su = tm_su + rt0
223	!
224	ENDIF
225	CALL wrk_dealloc( jpi, jpj, nlay_s, ze_s )
226	CALL wrk_dealloc( jpi, jpj, nlay_i, ze_i )
227	!
228	END SUBROUTINE lim_var_agg
229
230
231	SUBROUTINE lim_var_glo2eqv
232	!!------------------------------------------------------------------
233	!! * ROUTINE lim_var_glo2eqv *
234	!!
235	!! ** Purpose : computes equivalent variables as function of global variables
236	!! i.e. it turns VGLO into VEQV
237	!!------------------------------------------------------------------
238	INTEGER :: ji, jj, jk, jl ! dummy loop indices
239	REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars
240	REAL(wp) :: ztmelts, zq_s, zfac1, zfac2 ! - -
241	!!------------------------------------------------------------------
242
243	!-------------------------------------------------------
244	! Ice thickness, snow thickness, ice salinity, ice age
245	!-------------------------------------------------------
246	!$OMP PARALLEL
247	DO jl = 1, jpl
248	!$OMP DO schedule(static) private(jj,ji,rswitch)
249	DO jj = 1, jpj
250	DO ji = 1, jpi
251	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
252	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
253	END DO
254	END DO
255	END DO
256	! Force the upper limit of ht_i to always be < hi_max (99 m).
257	!$OMP DO schedule(static) private(jj,ji,rswitch)
258	DO jj = 1, jpj
259	DO ji = 1, jpi
260	rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
261	ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
262	a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
263	END DO
264	END DO
265
266	DO jl = 1, jpl
267	!$OMP DO schedule(static) private(jj,ji,rswitch)
268	DO jj = 1, jpj
269	DO ji = 1, jpi
270	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
271	ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
272	o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
273	END DO
274	END DO
275	END DO
276
277	IF( nn_icesal == 2 )THEN
278	DO jl = 1, jpl
279	!$OMP DO schedule(static) private(jj,ji,rswitch)
280	DO jj = 1, jpj
281	DO ji = 1, jpi
282	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes
283	sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch
284	! ! bounding salinity
285	sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin )
286	END DO
287	END DO
288	END DO
289	ENDIF
290	!$OMP END PARALLEL
291
292	CALL lim_var_salprof ! salinity profile
293
294	!-------------------
295	! Ice temperatures
296	!-------------------
297	!$OMP PARALLEL
298	DO jl = 1, jpl
299	DO jk = 1, nlay_i
300	!$OMP DO schedule(static) private(jj,ji,rswitch,zq_i,ztmelts,zaaa,zbbb,zccc,zdiscrim)
301	DO jj = 1, jpj
302	DO ji = 1, jpi
303	! ! Energy of melting q(S,T) [J.m-3]
304	rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes
305	zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)
306	ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature
307	!
308	zaaa = cpic ! Conversion q(S,T) -> T (second order equation)
309	zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus
310	zccc = lfus * (ztmelts-rt0)
311	zdiscrim = SQRT( MAX(zbbbzbbb - 4._wpzaaa*zccc , 0._wp) )
312	t_i(ji,jj,jk,jl) = rt0 + rswitch ( - zbbb - zdiscrim ) / ( 2.0 zaaa )
313	t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) ) ! -100 < t_i < ztmelts
314	END DO
315	END DO
316	END DO
317	END DO
318
319	!--------------------
320	! Snow temperatures
321	!--------------------
322	zfac1 = 1._wp / ( rhosn * cpic )
323	zfac2 = lfus / cpic
324	DO jl = 1, jpl
325	DO jk = 1, nlay_s
326	!$OMP DO schedule(static) private(jj,ji,rswitch,zq_s)
327	DO jj = 1, jpj
328	DO ji = 1, jpi
329	!Energy of melting q(S,T) [J.m-3]
330	rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes
331	zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
332	!
333	t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 )
334	t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) ) ! -100 < t_s < rt0
335	END DO
336	END DO
337	END DO
338	END DO
339
340	! integrated values
341	!$OMP DO schedule(static) private(jj, ji)
342	DO jj = 1, jpj
343	DO ji = 1, jpi
344	vt_i (ji,jj) = 0._wp
345	vt_s (ji,jj) = 0._wp
346	at_i (ji,jj) = 0._wp
347	END DO
348	END DO
349	DO jl = 1, jpl
350	!$OMP DO schedule(static) private(jj, ji)
351	DO jj = 1, jpj
352	DO ji = 1, jpi
353	vt_i (ji,jj) = vt_i (ji,jj) + v_i (ji,jj,jl)
354	vt_s (ji,jj) = vt_s (ji,jj) + v_s (ji,jj,jl)
355	at_i (ji,jj) = at_i (ji,jj) + a_i (ji,jj,jl)
356	END DO
357	END DO
358	END DO
359	!$OMP END PARALLEL
360	!
361	END SUBROUTINE lim_var_glo2eqv
362
363
364	SUBROUTINE lim_var_eqv2glo
365	!!------------------------------------------------------------------
366	!! * ROUTINE lim_var_eqv2glo *
367	!!
368	!! ** Purpose : computes global variables as function of equivalent variables
369	!! i.e. it turns VEQV into VGLO
370	!! ** Method :
371	!!
372	!! ** History : (01-2006) Martin Vancoppenolle, UCL-ASTR
373	!!------------------------------------------------------------------
374	!
375	v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
376	v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
377	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
378	!
379	END SUBROUTINE lim_var_eqv2glo
380
381
382	SUBROUTINE lim_var_salprof
383	!!------------------------------------------------------------------
384	!! * ROUTINE lim_var_salprof *
385	!!
386	!! ** Purpose : computes salinity profile in function of bulk salinity
387	!!
388	!! ** Method : If bulk salinity greater than zsi1,
389	!! the profile is assumed to be constant (S_inf)
390	!! If bulk salinity lower than zsi0,
391	!! the profile is linear with 0 at the surface (S_zero)
392	!! If it is between zsi0 and zsi1, it is a
393	!! alpha-weighted linear combination of s_inf and s_zero
394	!!
395	!! ** References : Vancoppenolle et al., 2007
396	!!------------------------------------------------------------------
397	INTEGER :: ji, jj, jk, jl ! dummy loop index
398	REAL(wp) :: zfac0, zfac1, zsal
399	REAL(wp) :: zswi0, zswi01, zargtemp , zs_zero
400	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha
401	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
402	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
403	!!------------------------------------------------------------------
404
405	CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
406
407	!---------------------------------------
408	! Vertically constant, constant in time
409	!---------------------------------------
410	IF( nn_icesal == 1 ) THEN
411	!$OMP PARALLEL
412	DO jl = 1, jpl
413	DO jk = 1, nlay_i
414	!$OMP DO schedule(static) private(jj, ji)
415	DO jj = 1, jpj
416	DO ji = 1, jpi
417	s_i (ji,jj,jk,jl) = rn_icesal
418	END DO
419	END DO
420	END DO
421	END DO
422	DO jl = 1, jpl
423	!$OMP DO schedule(static) private(jj, ji)
424	DO jj = 1, jpj
425	DO ji = 1, jpi
426	sm_i(ji,jj,jl) = rn_icesal
427	END DO
428	END DO
429	END DO
430	!$OMP END PARALLEL
431	ENDIF
432
433	!-----------------------------------
434	! Salinity profile, varying in time
435	!-----------------------------------
436	IF( nn_icesal == 2 ) THEN
437	!
438	!$OMP PARALLEL
439	DO jl = 1, jpl
440	DO jk = 1, nlay_i
441	!$OMP DO schedule(static) private(jj, ji)
442	DO jj = 1, jpj
443	DO ji = 1, jpi
444	s_i(ji,jj,jk,jl) = sm_i(ji,jj,jl)
445	END DO
446	END DO
447	!$OMP END DO NOWAIT
448	END DO
449	END DO
450	!
451	DO jl = 1, jpl ! Slope of the linear profile
452	!$OMP DO schedule(static) private(jj,ji,rswitch)
453	DO jj = 1, jpj
454	DO ji = 1, jpi
455	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
456	z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
457	END DO
458	END DO
459	!$OMP END DO NOWAIT
460	END DO
461	!
462	zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf
463	zfac1 = zsi1 / ( zsi1 - zsi0 )
464	!
465	DO jl = 1, jpl
466	!$OMP DO schedule(static) private(jj, ji)
467	DO jj = 1, jpj
468	DO ji = 1, jpi
469	zalpha(ji,jj,jl) = 0._wp
470	END DO
471	END DO
472	END DO
473	DO jl = 1, jpl
474	!$OMP DO schedule(static) private(jj,ji,zswi0,zswi01,rswitch)
475	DO jj = 1, jpj
476	DO ji = 1, jpi
477	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
478	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) )
479	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
480	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) )
481	! If 2.sm_i GE sss_m then rswitch = 1
482	! this is to force a constant salinity profile in the Baltic Sea
483	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
484	zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
485	zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
486	END DO
487	END DO
488	END DO
489
490	! Computation of the profile
491	DO jl = 1, jpl
492	DO jk = 1, nlay_i
493	!$OMP DO schedule(static) private(jj,ji,zs_zero)
494	DO jj = 1, jpj
495	DO ji = 1, jpi
496	! ! linear profile with 0 at the surface
497	zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
498	! ! weighting the profile
499	s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
500	! ! bounding salinity
501	s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
502	END DO
503	END DO
504	END DO
505	END DO
506	!$OMP END PARALLEL
507	!
508	ENDIF ! nn_icesal
509
510	!-------------------------------------------------------
511	! Vertically varying salinity profile, constant in time
512	!-------------------------------------------------------
513
514	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
515	!
516	!$OMP PARALLEL
517	DO jl = 1, jpl
518	!$OMP DO schedule(static) private(jj,ji)
519	DO jj = 1, jpj
520	DO ji = 1, jpi
521	sm_i(ji,jj,jl) = 2.30_wp
522	END DO
523	END DO
524	!$OMP END DO NOWAIT
525	END DO
526	!
527	DO jl = 1, jpl
528	DO jk = 1, nlay_i
529	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
530	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
531	!$OMP DO schedule(static) private(jj,ji)
532	DO jj = 1, jpj
533	DO ji = 1, jpi
534	s_i(ji,jj,jk,jl) = zsal
535	END DO
536	END DO
537	END DO
538	END DO
539	!$OMP END PARALLEL
540	!
541	ENDIF ! nn_icesal
542	!
543	CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
544	!
545	END SUBROUTINE lim_var_salprof
546
547
548	SUBROUTINE lim_var_bv
549	!!------------------------------------------------------------------
550	!! * ROUTINE lim_var_bv *
551	!!
552	!! ** Purpose : computes mean brine volume (%) in sea ice
553	!!
554	!! ** Method : e = - 0.054 * S (ppt) / T (C)
555	!!
556	!! References : Vancoppenolle et al., JGR, 2007
557	!!------------------------------------------------------------------
558	INTEGER :: ji, jj, jk, jl ! dummy loop indices
559	!!------------------------------------------------------------------
560	!
561	!$OMP PARALLEL
562	!$OMP DO schedule(static) private(jj,ji)
563	DO jj = 1, jpj
564	DO ji = 1, jpi
565	bvm_i(ji,jj) = 0._wp
566	END DO
567	END DO
568	DO jl = 1, jpl
569	!$OMP DO schedule(static) private(jj,ji)
570	DO jj = 1, jpj
571	DO ji = 1, jpi
572	bv_i (ji,jj,jl) = 0._wp
573	END DO
574	END DO
575	END DO
576	DO jl = 1, jpl
577	DO jk = 1, nlay_i
578	!$OMP DO schedule(static) private(jj,ji,rswitch)
579	DO jj = 1, jpj
580	DO ji = 1, jpi
581	rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) )
582	bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i &
583	& / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
584	END DO
585	END DO
586	END DO
587
588	!$OMP DO schedule(static) private(jj,ji,rswitch)
589	DO jj = 1, jpj
590	DO ji = 1, jpi
591	rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
592	bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
593	END DO
594	END DO
595	END DO
596	!$OMP END PARALLEL
597	!
598	END SUBROUTINE lim_var_bv
599
600
601	SUBROUTINE lim_var_salprof1d( kideb, kiut )
602	!!-------------------------------------------------------------------
603	!! * ROUTINE lim_thd_salprof1d *
604	!!
605	!! ** Purpose : 1d computation of the sea ice salinity profile
606	!! Works with 1d vectors and is used by thermodynamic modules
607	!!-------------------------------------------------------------------
608	INTEGER, INTENT(in) :: kideb, kiut ! thickness category index
609	!
610	INTEGER :: ji, jk ! dummy loop indices
611	INTEGER :: ii, ij ! local integers
612	REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars
613	REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - -
614	!
615	REAL(wp), POINTER, DIMENSION(:) :: z_slope_s
616	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
617	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
618	!!---------------------------------------------------------------------
619
620	CALL wrk_alloc( jpij, z_slope_s )
621
622	!---------------------------------------
623	! Vertically constant, constant in time
624	!---------------------------------------
625	IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal
626
627	!------------------------------------------------------
628	! Vertically varying salinity profile, varying in time
629	!------------------------------------------------------
630
631	IF( nn_icesal == 2 ) THEN
632	!
633	DO ji = kideb, kiut ! Slope of the linear profile zs_zero
634	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
635	z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
636	END DO
637
638	! Weighting factor between zs_zero and zs_inf
639	!---------------------------------------------
640	zfac0 = 1._wp / ( zsi0 - zsi1 )
641	zfac1 = zsi1 / ( zsi1 - zsi0 )
642	DO jk = 1, nlay_i
643	DO ji = kideb, kiut
644	ii = MOD( npb(ji) - 1 , jpi ) + 1
645	ij = ( npb(ji) - 1 ) / jpi + 1
646	! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
647	zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) )
648	! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
649	zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) )
650	! if 2.sm_i GE sss_m then rswitch = 1
651	! this is to force a constant salinity profile in the Baltic Sea
652	rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
653	!
654	zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch )
655	!
656	zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
657	! weighting the profile
658	s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
659	! bounding salinity
660	s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
661	END DO
662	END DO
663
664	ENDIF
665
666	!-------------------------------------------------------
667	! Vertically varying salinity profile, constant in time
668	!-------------------------------------------------------
669
670	IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
671	!
672	sm_i_1d(:) = 2.30_wp
673	!
674	DO jk = 1, nlay_i
675	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
676	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
677	DO ji = kideb, kiut
678	s_i_1d(ji,jk) = zsal
679	END DO
680	END DO
681	!
682	ENDIF
683	!
684	CALL wrk_dealloc( jpij, z_slope_s )
685	!
686	END SUBROUTINE lim_var_salprof1d
687
688	SUBROUTINE lim_var_zapsmall
689	!!-------------------------------------------------------------------
690	!! * ROUTINE lim_var_zapsmall *
691	!!
692	!! ** Purpose : Remove too small sea ice areas and correct fluxes
693	!!
694	!! history : LIM3.5 - 01-2014 (C. Rousset) original code
695	!!-------------------------------------------------------------------
696	INTEGER :: ji, jj, jl, jk ! dummy loop indices
697	REAL(wp) :: zsal, zvi, zvs, zei, zes
698	!!-------------------------------------------------------------------
699	!$OMP PARALLEL
700	!$OMP DO schedule(static) private(jj,ji)
701	DO jj = 1, jpj
702	DO ji = 1, jpi
703	at_i (ji,jj) = 0._wp
704	END DO
705	END DO
706	DO jl = 1, jpl
707	!$OMP DO schedule(static) private(jj,ji)
708	DO jj = 1, jpj
709	DO ji = 1, jpi
710	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
711	END DO
712	END DO
713	END DO
714
715	DO jl = 1, jpl
716
717	!-----------------------------------------------------------------
718	! Zap ice energy and use ocean heat to melt ice
719	!-----------------------------------------------------------------
720	DO jk = 1, nlay_i
721	!$OMP DO schedule(static) private(jj,ji,rswitch,zei)
722	DO jj = 1 , jpj
723	DO ji = 1 , jpi
724	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
725	rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch
726	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
727	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
728	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
729	zei = e_i(ji,jj,jk,jl)
730	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
731	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
732	! update exchanges with ocean
733	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
734	END DO
735	END DO
736	END DO
737
738	!$OMP DO schedule(static) private(jj,ji,rswitch,zsal,zvi,zvs,zes)
739	DO jj = 1 , jpj
740	DO ji = 1 , jpi
741	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
742	rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch
743	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
744	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
745	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
746	zsal = smv_i(ji,jj, jl)
747	zvi = v_i (ji,jj, jl)
748	zvs = v_s (ji,jj, jl)
749	zes = e_s (ji,jj,1,jl)
750	!-----------------------------------------------------------------
751	! Zap snow energy
752	!-----------------------------------------------------------------
753	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
754	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
755
756	!-----------------------------------------------------------------
757	! zap ice and snow volume, add water and salt to ocean
758	!-----------------------------------------------------------------
759	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
760	a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch
761	v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch
762	v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch
763	t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
764	oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
765	smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
766
767	! update exchanges with ocean
768	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
769	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
770	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
771	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
772	END DO
773	END DO
774	END DO
775
776	! to be sure that at_i is the sum of a_i(jl)
777	!$OMP DO schedule(static) private(jj,ji)
778	DO jj = 1, jpj
779	DO ji = 1, jpi
780	at_i (ji,jj) = 0._wp
781	END DO
782	END DO
783	DO jl = 1, jpl
784	!$OMP DO schedule(static) private(jj,ji)
785	DO jj = 1, jpj
786	DO ji = 1, jpi
787	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
788	END DO
789	END DO
790	END DO
791
792	! open water = 1 if at_i=0
793	!$OMP DO schedule(static) private(jj,ji,rswitch)
794	DO jj = 1, jpj
795	DO ji = 1, jpi
796	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
797	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
798	END DO
799	END DO
800	!$OMP END PARALLEL
801
802	!
803	END SUBROUTINE lim_var_zapsmall
804
805	SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
806	!!------------------------------------------------------------------
807	!! * ROUTINE lim_var_itd *
808	!!
809	!! ** Purpose : converting 1-cat ice to multiple ice categories
810	!!
811	!! ice thickness distribution follows a gaussian law
812	!! around the concentration of the most likely ice thickness
813	!! (similar as limistate.F90)
814	!!
815	!! ** Method: Iterative procedure
816	!!
817	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
818	!!
819	!! 2) Check whether the distribution conserves area and volume, positivity and
820	!! category boundaries
821	!!
822	!! 3) If not (input ice is too thin), the last category is empty and
823	!! the number of categories is reduced (jpl-1)
824	!!
825	!! 4) Iterate until ok (SUM(itest(:) = 4)
826	!!
827	!! ** Arguments : zhti: 1-cat ice thickness
828	!! zhts: 1-cat snow depth
829	!! zai : 1-cat ice concentration
830	!!
831	!! ** Output : jpl-cat
832	!!
833	!! (Example of application: BDY forcings when input are cell averaged)
834	!!
835	!!-------------------------------------------------------------------
836	!! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code
837	!! 2014 (C. Rousset) Rewriting
838	!!-------------------------------------------------------------------
839	!! Local variables
840	INTEGER :: ji, jk, jl ! dummy loop indices
841	INTEGER :: ijpij, i_fill, jl0
842	REAL(wp) :: zarg, zV, zconv, zdh, zdv
843	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables
844	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables
845	INTEGER , POINTER, DIMENSION(:) :: itest
846
847	CALL wrk_alloc( 4, itest )
848	!--------------------------------------------------------------------
849	! initialisation of variables
850	!--------------------------------------------------------------------
851	ijpij = SIZE(zhti,1)
852	zht_i(1:ijpij,1:jpl) = 0._wp
853	zht_s(1:ijpij,1:jpl) = 0._wp
854	za_i (1:ijpij,1:jpl) = 0._wp
855
856	! ----------------------------------------
857	! distribution over the jpl ice categories
858	! ----------------------------------------
859	DO ji = 1, ijpij
860
861	IF( zhti(ji) > 0._wp ) THEN
862
863	! find which category (jl0) the input ice thickness falls into
864	jl0 = jpl
865	DO jl = 1, jpl
866	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
867	jl0 = jl
868	CYCLE
869	ENDIF
870	END DO
871
872	! initialisation of tests
873	itest(:) = 0
874
875	i_fill = jpl + 1 !====================================
876	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
877	! iteration !====================================
878	i_fill = i_fill - 1
879
880	! initialisation of ice variables for each try
881	zht_i(ji,1:jpl) = 0._wp
882	za_i (ji,1:jpl) = 0._wp
883	itest(:) = 0
884
885	! *** case very thin ice: fill only category 1
886	IF ( i_fill == 1 ) THEN
887	zht_i(ji,1) = zhti(ji)
888	za_i (ji,1) = zai (ji)
889
890	! *** case ice is thicker: fill categories >1
891	ELSE
892
893	! Fill ice thicknesses in the (i_fill-1) cat by hmean
894	DO jl = 1, i_fill - 1
895	zht_i(ji,jl) = hi_mean(jl)
896	END DO
897
898	! Concentrations in the (i_fill-1) categories
899	za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
900	DO jl = 1, i_fill - 1
901	IF ( jl /= jl0 ) THEN
902	zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
903	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
904	ENDIF
905	END DO
906
907	! Concentration in the last (i_fill) category
908	za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
909
910	! Ice thickness in the last (i_fill) category
911	zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
912	zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
913
914	! clem: correction if concentration of upper cat is greater than lower cat
915	! (it should be a gaussian around jl0 but sometimes it is not)
916	IF ( jl0 /= jpl ) THEN
917	DO jl = jpl, jl0+1, -1
918	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
919	zdv = zht_i(ji,jl) * za_i(ji,jl)
920	zht_i(ji,jl ) = 0._wp
921	za_i (ji,jl ) = 0._wp
922	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
923	END IF
924	ENDDO
925	ENDIF
926
927	ENDIF ! case ice is thick or thin
928
929	!---------------------
930	! Compatibility tests
931	!---------------------
932	! Test 1: area conservation
933	zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
934	IF ( zconv < epsi06 ) itest(1) = 1
935
936	! Test 2: volume conservation
937	zconv = ABS( zhti(ji)zai(ji) - SUM( za_i(ji,1:jpl)zht_i(ji,1:jpl) ) )
938	IF ( zconv < epsi06 ) itest(2) = 1
939
940	! Test 3: thickness of the last category is in-bounds ?
941	IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
942
943	! Test 4: positivity of ice concentrations
944	itest(4) = 1
945	DO jl = 1, i_fill
946	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
947	END DO
948	! !============================
949	END DO ! end iteration on categories
950	! !============================
951	ENDIF ! if zhti > 0
952	END DO ! i loop
953
954	! ------------------------------------------------
955	! Adding Snow in each category where za_i is not 0
956	! ------------------------------------------------
957	DO jl = 1, jpl
958	DO ji = 1, ijpij
959	IF( za_i(ji,jl) > 0._wp ) THEN
960	zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
961	! In case snow load is in excess that would lead to transformation from snow to ice
962	! Then, transfer the snow excess into the ice (different from limthd_dh)
963	zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )
964	! recompute ht_i, ht_s avoiding out of bounds values
965	zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
966	zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
967	ENDIF
968	ENDDO
969	ENDDO
970
971	CALL wrk_dealloc( 4, itest )
972	!
973	END SUBROUTINE lim_var_itd
974
975
976	#else
977	!!----------------------------------------------------------------------
978	!! Default option Dummy module NO LIM3 sea-ice model
979	!!----------------------------------------------------------------------
980	CONTAINS
981	SUBROUTINE lim_var_agg ! Empty routines
982	END SUBROUTINE lim_var_agg
983	SUBROUTINE lim_var_glo2eqv ! Empty routines
984	END SUBROUTINE lim_var_glo2eqv
985	SUBROUTINE lim_var_eqv2glo ! Empty routines
986	END SUBROUTINE lim_var_eqv2glo
987	SUBROUTINE lim_var_salprof ! Empty routines
988	END SUBROUTINE lim_var_salprof
989	SUBROUTINE lim_var_bv ! Emtpy routines
990	END SUBROUTINE lim_var_bv
991	SUBROUTINE lim_var_salprof1d ! Emtpy routines
992	END SUBROUTINE lim_var_salprof1d
993	SUBROUTINE lim_var_zapsmall
994	END SUBROUTINE lim_var_zapsmall
995	SUBROUTINE lim_var_itd
996	END SUBROUTINE lim_var_itd
997	#endif
998
999	!!======================================================================
1000	END MODULE limvar

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

Download in other formats: