Context Navigation

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

Last change on this file since 59 was 59, checked in by aquiquet, 8 years ago
Introduction of Betamax_2d for all geoms / call to sliding removed / Drastic change in dragging module based on effective pressure (no curvature but slow_ssa zones)
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
340
341	end subroutine gzm_beta_prescr
342
343
344	!----------------------------------------------------------------------------------
345	! Routines qui permettent des manip simples sur beta.
346
347
348	!______________________________________________________________________________________
349	!>SUBROUTINE: beta_plateau_slope
350	!! give value beta_limgz to places where slope is >1.e-2
351	!<
352
353	subroutine beta_plateau_slope ! give value beta_limgz to places where slope is >1.e-2
354
355	call slope_surf
356
357	where (abs(sdx(:,:)).gt.1.e-2)
358	drag_mx(:,:) = beta_limgz
359	end where
360
361	where (abs(sdy(:,:)).gt.1.e-2)
362	drag_my(:,:) = beta_limgz
363	end where
364	end subroutine beta_plateau_slope
365	!______________________________________________________________________________________
366	!> SUBROUTINE beta_plateau_S
367	!! give value beta_limgz to places where S > 2500 m
368	!<
369	subroutine beta_plateau_S ! give value beta_limgz to places where S > 2500 m
370
371	where (S(:,:).gt.2500.) ! moyen simple de nettoyer l'interieur. WAIS pas touche
372	drag_centre(:,:)=beta_limgz
373	drag_mx(:,:)=beta_limgz
374	drag_my(:,:)=beta_limgz
375	end where
376	end subroutine beta_plateau_S
377	!______________________________________________________________________________________
378	!>SUBROUTINE: beta_min_value
379	!! valeur mini que peut avoir beta (en Pa an m-1)
380	!<
381
382	subroutine beta_min_value(val_betamin)
383
384	real :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
385
386	drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
387	drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
388	drag_my(:,:) = max(drag_my(:,:),val_betamin)
389	where(flot(:,:))
390	drag_centre(:,:) = 0.
391	end where
392
393	end subroutine beta_min_value
394	!______________________________________________________________________________________
395	!>SUBROUTINE: beta_stag2c
396	!! average the staggered values on central ones
397	!<
398
399	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
400
401	implicit none
402	integer :: nxx,nyy ! dimension des tableaux
403	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
404	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
405	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
406
407	R_centre(:,:)=0.
408	do j=2,ny-1
409	do i=2,nx-1
410	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
411	+ (R_my(i,j)+R_my(i,j+1)))*0.25
412	end do
413	end do
414	end subroutine beta_stag2c
415	!______________________________________________________________________________________
416	!>SUBROUTINE: beta_c2stag_min
417	!! redistribute central values on staggered grid
418	!<
419
420	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
421
422	implicit none
423	integer :: nxx,nyy ! dimension des tableaux
424	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
425	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
426	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
427
428	do j=2,nyy
429	do i=2,nxx
430	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
431	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
432	end do
433	end do
434
435	end subroutine beta_c2stag
436	!______________________________________________________________________________________
437	!>SUBROUTINE: beta_c2stag_min
438	!! redistribute central values on staggered grid
439	!<
440
441	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
442
443	implicit none
444	integer :: nxx,nyy ! dimension des tableaux
445	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
446	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
447	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
448
449	do j=2,nyy
450	do i=2,nxx
451	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
452	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
453	end do
454	end do
455
456	end subroutine beta_c2stag_min
457	!______________________________________________________________________________________
458	!>SUBROUTINE: make_gaussienne
459	!! cree une matrice gaussienne normalisee a 1
460	!<
461
462	subroutine make_gaussienne(n,sigma,M)
463
464	implicit none
465	integer,intent(in) :: n ! demi taille du tableau gaussienne
466	real, intent(in) :: sigma ! sigma gaussienne exprime en dx
467	real,dimension(-n:n,-n:n),intent(inout) :: M ! tableau gaussienne
468
469
470	! variables locales
471	real :: dist2 ! distance au centre (au carre)
472	real :: som ! pour calibration
473	real :: sigma2 ! pour gaussienne
474
475	som = 0.
476	sigma2 = sigma * sigma * 2
477	do j = -n, n
478	do i = -n, n
479	dist2 = (ii+jj)
480	M(i,j) = exp(-dist2/sigma2)
481	som=som+M(i,j)
482	end do
483	end do
484
485	M(:,:) = M(:,:) / som
486
487	! do j = -n, n
488	! write(6,'(5(f0.4,1x))') (M(i,j),i=-n,n)
489	! end do
490
491	end subroutine make_gaussienne
492
493	!______________________________________________________________________________________
494	!>SUBROUTINE: iter_gauss
495	!! utilise hdot pour faire les iterations sur beta_centre
496	!<
497	subroutine iter_gauss(n,nxx,nyy,M,coef,beta_in,beta_out,Ecart)
498
499	integer, intent(in) :: n ! taille matrice gaussienne
500	integer, intent(in) :: nxx,nyy ! taille du domaine
501	real,dimension(-n:n,-n:n),intent(in) :: M ! matrice gaussienne
502	real, intent(in) :: coef ! coefficient multiplicateur
503	real,dimension(nx,ny), intent(in) :: beta_in ! la valeur initiale
504	real,dimension(nx,ny), intent(out) :: beta_out ! la valeur apres iteration
505	real,dimension(nx,ny), intent(in) :: Ecart ! l'ecart a l'obs
506
507	! variables locales
508	integer :: ii,jj
509
510	do j= 1 + n, nyy - n
511	do i = 1 + n, nxx - n
512	beta_out(i,j) = beta_in(i,j)
513
514	do jj = -n, n
515	do ii = -n, n
516	beta_out(i,j) = beta_out(i,j) + coef * M(ii,jj) * ecart(i+ii,j+jj)
517	end do
518	end do
519
520	end do
521	end do
522
523	beta_out(:,:) = max (beta_out(:,:),0.)
524	end subroutine iter_gauss
525
526	end module dragging_prescr_beta
527
528	!!$! add the nodes of topographical instability (on major nodes)
529	! a garder eventuellement pour runs instabilite
530	!!$beta_file = trim(dirnameinp)//'sensib_H_Gael-15km.dat'
531	!!$call lect_datfile(nx,ny,drag_instab_topo,1,beta_file)
532	!!$
533	!!$where (drag_instab_topo(:,:).gt.0.)
534	!!$ drag_centre(:,:)=min(drag_centre(:,:),400.)
535	!!$endwhere

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

Download in other formats: