Context Navigation

source: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_mod.f90 @ 159

Last change on this file since 159 was 70, checked in by dumas, 8 years ago
Bug fix in Qprod_icetemp (instability in ice temperature) : slowssa points heat = deformation SIA + sliding. Use fleuvemx and fleuvemy to identifie streams SSA and slow SSA points
File size: 18.5 KB

Line
1	!> \file dragging_prescr_beta_mod.f90
2	!! Module pour calculer le beta a partir d'une carte sur les noeuds centres.
3	!<
4
5	!> \namespace dragging_prescr_beta
6	!! Calcule le beta a partir de vitesses de bilan
7	!! @note Il faut partir d'un cptr
8	!! \author ...
9	!! \date ...
10	!! @note Used module
11	!! @note - use module3D_phy
12	!<
13
14	module dragging_prescr_beta
15
16
17	use module3d_phy
18	use interface_input
19	implicit none
20	real, dimension(nx,ny) :: Vcol_x !< uniquement pour compatibilite
21	real, dimension(nx,ny) :: Vcol_y !<
22	real, dimension(nx,ny) :: Vsl_x !<
23	real, dimension(nx,ny) :: Vsl_y !<
24
25
26	real, dimension(nx,ny) :: drag_centre !< beta on major node (average)
27	real, dimension(nx,ny) :: betamoy_x !< pour faire les moyennes glissantes
28	real, dimension(nx,ny) :: betamoy_y !<
29	integer, dimension(nx,ny) :: nbvoisx !<
30	integer, dimension(nx,ny) :: nbvoisy !<
31
32
33	real :: beta_limgz !< when beta gt beta_limgz -> not gzmx
34	real :: beta_min !< minimum value of beta
35	real :: beta_mult !< coefficient of beta field
36	integer :: ill,jll,nmoy
37
38	real :: maxi !< calcul correction dS a la grounding line
39	real :: mini
40
41	logical :: corr_def = .False. !< for deformation correction (compatibility)
42
43	contains
44
45	!-----------------------------------------------------------------------------------------
46	!> SUBROUTINE: init_dragging
47	!! Cette routine fait l'initialisation du dragging.
48	!<
49	subroutine init_dragging ! Cette routine fait l'initialisation du dragging.
50	! nouvelle version : lit les fichiers nc
51	implicit none
52	character(len=100) :: beta_c_file ! beta on centered grid
53	character(len=100) :: betamx_file ! beta on staggered grid mx
54	character(len=100) :: betamy_file ! beta on staggered grid mx
55
56	namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
57
58	if (itracebug.eq.1) call tracebug(' Prescr_beta subroutine init_dragging')
59
60	iter_beta=0
61
62	! lecture des parametres du run
63	!--------------------------------------------------------------------
64
65	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
66	428 format(A)
67	read(num_param,beta_prescr)
68
69	write(num_rep_42,428) '!___________________________________________________________'
70	write(num_rep_42,428) '! read beta on centered grid'
71	write(num_rep_42,beta_prescr)
72	write(num_rep_42,428) '! beta_file : nom des fichier qui contiennent les beta_c'
73	write(num_rep_42,428) '! above beta_limgz, gzmx is false'
74	write(num_rep_42,428) '! if grounded, beta_min minimum value '
75	write(num_rep_42,428) '! beta_mult : coefficient just after reading'
76	write(num_rep_42,*)
77	write(num_rep_42,428) '!___________________________________________________________'
78
79
80	beta_c_file = trim(dirnameinp)//trim(beta_c_file)
81	betamx_file = trim(dirnameinp)//trim(betamx_file)
82	betamy_file = trim(dirnameinp)//trim(betamy_file)
83	betamax = beta_limgz
84	betamax_2d(:,:) = betamax
85
86	! read the beta value on centered and staggered grids
87	! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
88	! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
89	! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
90
91
92	call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
93
94	! multiplication par le coefficient beta_mult
95
96	drag_centre(:,:) = drag_centre(:,:)*beta_mult
97
98	where (mk_init(:,:).eq.-2) ! iles a probleme
99	drag_centre(:,:) = 2.* abs(beta_limgz)
100	end where
101
102
103	where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes
104	drag_centre(:,:) = 1000. ! valeur relativement faible qui evite
105	end where ! les freinages
106
107
108
109	! call lect_input(1,'betamx',1,drag_mx,betamx_file,'')
110	! call lect_input(1,'betamy',1,drag_my,betamy_file,'')
111
112
113	call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
114
115	if (beta_limgz.gt.0.) then
116	beta_centre(:,:) = drag_centre(:,:)
117	betamx(:,:) = drag_mx(:,:)
118	betamy(:,:) = drag_my(:,:)
119
120	else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic
121
122	beta_centre(:,:) = - beta_limgz
123	betamx(:,:) = - beta_limgz
124	betamy(:,:) = - beta_limgz
125	drag_centre(:,:) = - beta_limgz
126	drag_mx(:,:) = - beta_limgz
127	drag_my(:,:) = - beta_limgz
128	beta_limgz = 1.
129	end if
130
131
132	! divers traitements qui peuvent etre appliques a beta : mis sous forme de routine
133
134	! call beta_plateau_slope ! give value beta_limgz to places where slope is >1.e-2
135	! call beta_stag2c ! average the staggered values on central ones
136	! call beta_plateau_S ! give value beta_limgz to places where S > 2500 m
137
138	call beta_min_value(beta_min) ! valeur mini que peut avoir beta (en Pa an m-1)
139	! on peut envisager une valeur ~ 10
140	! rappel : beta doit être positif.
141
142	! call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
143
144
145	! mstream attribution
146	!________________________
147	! mstream -> regions where streaming is allowed
148	! drag_mx -> value of beta if an ice stream is allowed
149	! betamx -> value of dragging (may depend on other conditions)
150
151	! to ensure having some "Dirichlet" points, there must be at least
152	! a mstream=0 in every island
153
154
155	! un tour de lissage
156	! call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre) ! average staggered
157	! ! on central
158
159	! call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! centered distributed
160	! ! on staggered
161	!
162
163
164	beta_centre(:,:) = drag_centre(:,:) ! surtout pour sorties
165
166
167
168	where (drag_mx(:,:).ge.beta_limgz)
169	mstream_mx(:,:) = 0
170	betamx(:,:) = beta_limgz
171	drag_mx(:,:) = beta_limgz
172	elsewhere
173	mstream_mx(:,:) = 1
174	betamx(:,:) = drag_mx(:,:)
175	endwhere
176
177	where (drag_my(:,:).ge.beta_limgz)
178	mstream_my(:,:) = 0
179	betamy(:,:) = beta_limgz
180	drag_my(:,:) = beta_limgz
181	elsewhere
182	mstream_my(:,:) = 1
183	betamy(:,:) = drag_my(:,:)
184	endwhere
185
186	if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture')
187	mstream_mx(:,:) = 0
188	mstream_my(:,:) = 0
189
190
191	call gzm_beta_prescr
192	! call write_datfile2(nx,ny,betamx,betamy,'beta-LBq15-08.dat')
193
194	! call make_gaussienne(n_g,sigma_beta,M_gauss)
195
196	return
197
198	end subroutine init_dragging
199	!________________________________________________________________________________
200
201	!> SUBROUTINE: dragging
202	!! Defini la localisation des streams et le frottement basal
203	!<
204	!-------------------------------------------------------------------------
205	subroutine dragging ! defini la localisation des streams et le frottement basal
206
207	implicit none
208
209
210
211	if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging')
212
213
214	! tentative d'iteration
215
216	! call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre) ! average staggered
217	! on central
218	! call iter_gauss(n_g,nx,ny,M_gauss,-1.,drag_centre,corrstream_centre,hdot) ! correct central
219	! call beta_c2stag(nx,ny,drag_mx,drag_my,corrstream_centre) ! centered distributed
220	! on staggered
221	! drag_centre(:,:) = corrstream_centre(:,:)
222	! beta_centre(:,:) = drag_centre(:,:) ! surtout pour sorties
223
224	! pour l'instant pas d'autre condition que mstream
225	where (mstream_mx(:,:).eq.1)
226	betamx(:,:) = drag_mx(:,:)
227	elsewhere
228	betamx(:,:)= beta_limgz
229	end where
230
231	where (mstream_my(:,:).eq.1)
232	betamy(:,:) = drag_my(:,:)
233	elsewhere
234	betamy(:,:)= beta_limgz
235	end where
236
237	betamx(:,:)=drag_mx(:,:)
238	betamy(:,:)=drag_my(:,:)
239
240	call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants
241
242
243	return
244	end subroutine dragging
245
246	!----------------------------------------------------------------------------------
247	!>SUBROUTINE: gzm_beta_prescr
248	!! Calcul de gzmx
249	!<
250
251	subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my)
252
253	logical :: test_Tx ! test if there is a basal melting point in the surrounding
254	logical :: test_Ty ! test if there is a basal melting point in the surrounding
255	logical :: test_beta_x ! test if there is a low drag point in the surrounding
256	logical :: test_beta_y ! test if there is a low drag point in the surrounding
257
258	real :: beta_forc_fleuve ! below this value -> ice stream
259
260	! beta_forc_fleuve = 5.e3
261	beta_forc_fleuve = beta_limgz
262
263	! calcul de gzmx
264
265
266	gzmx(:,:)=.true.
267	gzmy(:,:)=.true.
268	flgzmx(:,:)=.false.
269	flgzmy(:,:)=.false.
270
271
272	! points flottants : flgzmx mais pas gzmx
273	do j=2,ny
274	do i=2,nx
275	if (flot(i,j).and.(flot(i-1,j))) then
276	flotmx(i,j)=.true.
277	flgzmx(i,j)=.true.
278	gzmx(i,j)=.false.
279	betamx(i,j)=0.
280	end if
281	if (flot(i,j).and.(flot(i,j-1))) then
282	flotmy(i,j)=.true.
283	flgzmy(i,j)=.true.
284	gzmy(i,j)=.false.
285	betamy(i,j)=0.
286	end if
287	end do
288	end do
289
290	if (itracebug.eq.1) call tracebug(' apres criteres flot')
291
292	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
293
294	!--------- autres criteres
295
296	! points poses
297	gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points
298	gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta
299
300
301
302	! critere (lache) sur la temperature
303	do j = 2, ny-1
304	do i =2, nx-1
305
306	! test s'il y a un point tempere dans les environs
307	test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1)
308	test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1)
309
310	! test s'il y a un point low drag dans les environs
311	test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve )
312	test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve )
313
314	! critere pour gzmx
315
316	! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
317
318	! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x)
319	! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y)
320
321	end do
322	end do
323
324
325	if (itracebug.eq.1) call tracebug(' apres autres criteres ')
326	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
327
328	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
329	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
330
331	where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
332	betamx(:,:) = beta_limgz
333	end where
334
335	where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
336	betamy(:,:) = beta_limgz
337	end where
338
339	fleuvemx(:,:)=gzmx(:,:)
340	fleuvemy(:,:)=gzmy(:,:)
341
342	end subroutine gzm_beta_prescr
343
344
345	!----------------------------------------------------------------------------------
346	! Routines qui permettent des manip simples sur beta.
347
348
349	!______________________________________________________________________________________
350	!>SUBROUTINE: beta_plateau_slope
351	!! give value beta_limgz to places where slope is >1.e-2
352	!<
353
354	subroutine beta_plateau_slope ! give value beta_limgz to places where slope is >1.e-2
355
356	call slope_surf
357
358	where (abs(sdx(:,:)).gt.1.e-2)
359	drag_mx(:,:) = beta_limgz
360	end where
361
362	where (abs(sdy(:,:)).gt.1.e-2)
363	drag_my(:,:) = beta_limgz
364	end where
365	end subroutine beta_plateau_slope
366	!______________________________________________________________________________________
367	!> SUBROUTINE beta_plateau_S
368	!! give value beta_limgz to places where S > 2500 m
369	!<
370	subroutine beta_plateau_S ! give value beta_limgz to places where S > 2500 m
371
372	where (S(:,:).gt.2500.) ! moyen simple de nettoyer l'interieur. WAIS pas touche
373	drag_centre(:,:)=beta_limgz
374	drag_mx(:,:)=beta_limgz
375	drag_my(:,:)=beta_limgz
376	end where
377	end subroutine beta_plateau_S
378	!______________________________________________________________________________________
379	!>SUBROUTINE: beta_min_value
380	!! valeur mini que peut avoir beta (en Pa an m-1)
381	!<
382
383	subroutine beta_min_value(val_betamin)
384
385	real :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
386
387	drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
388	drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
389	drag_my(:,:) = max(drag_my(:,:),val_betamin)
390	where(flot(:,:))
391	drag_centre(:,:) = 0.
392	end where
393
394	end subroutine beta_min_value
395	!______________________________________________________________________________________
396	!>SUBROUTINE: beta_stag2c
397	!! average the staggered values on central ones
398	!<
399
400	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
401
402	implicit none
403	integer :: nxx,nyy ! dimension des tableaux
404	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
405	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
406	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
407
408	R_centre(:,:)=0.
409	do j=2,ny-1
410	do i=2,nx-1
411	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
412	+ (R_my(i,j)+R_my(i,j+1)))*0.25
413	end do
414	end do
415	end subroutine beta_stag2c
416	!______________________________________________________________________________________
417	!>SUBROUTINE: beta_c2stag_min
418	!! redistribute central values on staggered grid
419	!<
420
421	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
422
423	implicit none
424	integer :: nxx,nyy ! dimension des tableaux
425	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
426	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
427	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
428
429	do j=2,nyy
430	do i=2,nxx
431	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
432	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
433	end do
434	end do
435
436	end subroutine beta_c2stag
437	!______________________________________________________________________________________
438	!>SUBROUTINE: beta_c2stag_min
439	!! redistribute central values on staggered grid
440	!<
441
442	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
443
444	implicit none
445	integer :: nxx,nyy ! dimension des tableaux
446	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
447	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
448	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
449
450	do j=2,nyy
451	do i=2,nxx
452	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
453	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
454	end do
455	end do
456
457	end subroutine beta_c2stag_min
458	!______________________________________________________________________________________
459	!>SUBROUTINE: make_gaussienne
460	!! cree une matrice gaussienne normalisee a 1
461	!<
462
463	subroutine make_gaussienne(n,sigma,M)
464
465	implicit none
466	integer,intent(in) :: n ! demi taille du tableau gaussienne
467	real, intent(in) :: sigma ! sigma gaussienne exprime en dx
468	real,dimension(-n:n,-n:n),intent(inout) :: M ! tableau gaussienne
469
470
471	! variables locales
472	real :: dist2 ! distance au centre (au carre)
473	real :: som ! pour calibration
474	real :: sigma2 ! pour gaussienne
475
476	som = 0.
477	sigma2 = sigma * sigma * 2
478	do j = -n, n
479	do i = -n, n
480	dist2 = (ii+jj)
481	M(i,j) = exp(-dist2/sigma2)
482	som=som+M(i,j)
483	end do
484	end do
485
486	M(:,:) = M(:,:) / som
487
488	! do j = -n, n
489	! write(6,'(5(f0.4,1x))') (M(i,j),i=-n,n)
490	! end do
491
492	end subroutine make_gaussienne
493
494	!______________________________________________________________________________________
495	!>SUBROUTINE: iter_gauss
496	!! utilise hdot pour faire les iterations sur beta_centre
497	!<
498	subroutine iter_gauss(n,nxx,nyy,M,coef,beta_in,beta_out,Ecart)
499
500	integer, intent(in) :: n ! taille matrice gaussienne
501	integer, intent(in) :: nxx,nyy ! taille du domaine
502	real,dimension(-n:n,-n:n),intent(in) :: M ! matrice gaussienne
503	real, intent(in) :: coef ! coefficient multiplicateur
504	real,dimension(nx,ny), intent(in) :: beta_in ! la valeur initiale
505	real,dimension(nx,ny), intent(out) :: beta_out ! la valeur apres iteration
506	real,dimension(nx,ny), intent(in) :: Ecart ! l'ecart a l'obs
507
508	! variables locales
509	integer :: ii,jj
510
511	do j= 1 + n, nyy - n
512	do i = 1 + n, nxx - n
513	beta_out(i,j) = beta_in(i,j)
514
515	do jj = -n, n
516	do ii = -n, n
517	beta_out(i,j) = beta_out(i,j) + coef * M(ii,jj) * ecart(i+ii,j+jj)
518	end do
519	end do
520
521	end do
522	end do
523
524	beta_out(:,:) = max (beta_out(:,:),0.)
525	end subroutine iter_gauss
526
527	end module dragging_prescr_beta
528
529	!!$! add the nodes of topographical instability (on major nodes)
530	! a garder eventuellement pour runs instabilite
531	!!$beta_file = trim(dirnameinp)//'sensib_H_Gael-15km.dat'
532	!!$call lect_datfile(nx,ny,drag_instab_topo,1,beta_file)
533	!!$
534	!!$where (drag_instab_topo(:,:).gt.0.)
535	!!$ drag_centre(:,:)=min(drag_centre(:,:),400.)
536	!!$endwhere

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

Download in other formats: