Context Navigation

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

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

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

Download in other formats: