Context Navigation

source: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_buoyency_mod.f90 @ 70

Last change on this file since 70 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: 17.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	!< Cette version est aussi fonction de la pression effective
4	!< calculee d'apres la hauteur au dessus de la flottaison
5	!<
6
7	!> \namespace dragging_prescr_beta
8	!! Calcule le beta a partir de vitesses de bilan
9	!! @note Il faut partir d'un cptr
10	!! \author ...
11	!! \date ...
12	!! @note Used module
13	!! @note - use module3D_phy
14	!<
15
16	module dragging_beta_buoy
17
18
19	use module3d_phy
20	use interface_input
21	implicit none
22	real, dimension(nx,ny) :: Vcol_x !< uniquement pour compatibilite
23	real, dimension(nx,ny) :: Vcol_y !<
24	real, dimension(nx,ny) :: Vsl_x !<
25	real, dimension(nx,ny) :: Vsl_y !<
26
27
28	real, dimension(nx,ny) :: drag_centre !< beta on major node (average)
29	!< beta_centre idem est la valeur courante declare dans 3D
30
31
32	real, dimension(nx,ny) :: H_buoy_init !< hauteur au dessus de la flottaison a initialisation
33	real, dimension(nx,ny) :: H_buoy_courant !< hauteur courante au dessus de la flottaison
34
35
36	real, dimension(nx,ny) :: betamoy_x !< pour faire les moyennes glissantes
37	real, dimension(nx,ny) :: betamoy_y !<
38	integer, dimension(nx,ny) :: nbvoisx !<
39	integer, dimension(nx,ny) :: nbvoisy !<
40
41
42	real :: beta_limgz !< when beta gt beta_limgz -> not gzmx
43	real :: beta_min !< minimum value of beta
44	real :: beta_mult !< coefficient of beta field
45	integer :: ill,jll,nmoy
46
47	real :: maxi !< calcul correction dS a la grounding line
48	real :: mini
49
50	logical :: corr_def = .False. !< for deformation correction (compatibility)
51
52	! declares dans le 3D (variables globales)
53	!-------------------------------------------
54	! beta_centre : valeur centree courante
55	! beta_mx et beta_my sont les valeurs courantes
56	! drag_mx, drag_my servent dans calc_beta (remplimat-shelves-tabTu.f90) avec un autre sens
57	! et de variable de calcul ici.
58
59
60
61	contains
62
63	!-----------------------------------------------------------------------------------------
64	!> SUBROUTINE: init_dragging
65	!! Cette routine fait l'initialisation du dragging.
66	!<
67	subroutine init_dragging ! Cette routine fait l'initialisation du dragging.
68	! nouvelle version : lit les fichiers nc
69	implicit none
70	character(len=100) :: beta_c_file ! beta on centered grid
71	character(len=100) :: betamx_file ! beta on staggered grid mx
72	character(len=100) :: betamy_file ! beta on staggered grid mx
73
74	namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
75
76	if (itracebug.eq.1) call tracebug(' Prescr_beta subroutine init_dragging')
77
78	iter_beta=0
79
80	! lecture des parametres du run
81	!--------------------------------------------------------------------
82
83	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
84	428 format(A)
85	read(num_param,beta_prescr)
86
87	write(num_rep_42,428) '!___________________________________________________________'
88	write(num_rep_42,428) '! read beta on centered grid'
89	write(num_rep_42,beta_prescr)
90	write(num_rep_42,428) '! beta_file : nom des fichier qui contiennent les beta_c'
91	write(num_rep_42,428) '! above beta_limgz, gzmx is false'
92	write(num_rep_42,428) '! if grounded, beta_min minimum value '
93	write(num_rep_42,428) '! beta_mult : coefficient just after reading'
94	write(num_rep_42,*)
95	write(num_rep_42,428) '!___________________________________________________________'
96
97
98	beta_c_file = trim(dirnameinp)//trim(beta_c_file)
99	betamx_file = trim(dirnameinp)//trim(betamx_file)
100	betamy_file = trim(dirnameinp)//trim(betamy_file)
101	betamax = beta_limgz
102	betamax_2d(:,:) = betamax
103
104	! read the beta value on centered and staggered grids
105	! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
106	! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
107	! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
108
109
110	call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
111
112	! calcule la valeur de reference de Hbuoy la hauteur au dessus de la flottaison
113	call calc_buoyency(nx,ny,H,B,H_buoy_init)
114
115
116	! multiplication du frottement par le coefficient beta_mult
117
118	drag_centre(:,:) = drag_centre(:,:)*beta_mult
119
120	where (mk_init(:,:).eq.-2) ! iles a probleme
121	drag_centre(:,:) = 2.* abs(beta_limgz)
122	end where
123
124
125	where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes
126	drag_centre(:,:) = 1000. ! valeur relativement faible qui evite
127	end where ! les freinages
128
129
130	call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
131
132	if (beta_limgz.gt.0.) then
133	beta_centre(:,:) = drag_centre(:,:)
134	betamx(:,:) = drag_mx(:,:)
135	betamy(:,:) = drag_my(:,:)
136
137	else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic
138
139	beta_centre(:,:) = - beta_limgz
140	betamx(:,:) = - beta_limgz
141	betamy(:,:) = - beta_limgz
142	drag_centre(:,:) = - beta_limgz
143	drag_mx(:,:) = - beta_limgz
144	drag_my(:,:) = - beta_limgz
145	beta_limgz = 1.
146	end if
147
148
149
150	! valeur mini que peut avoir beta (en Pa an m-1)
151	call beta_min_value(nx,ny,drag_mx,drag_my,drag_centre,beta_min)
152
153	! on peut envisager une valeur ~ 10
154	! rappel : beta doit être positif.
155
156
157
158	beta_centre(:,:) = drag_centre(:,:) ! utilise aussi dans schoof
159
160
161
162	where (drag_mx(:,:).ge.beta_limgz)
163	mstream_mx(:,:) = 0
164	betamx(:,:) = beta_limgz
165	drag_mx(:,:) = beta_limgz
166	elsewhere
167	mstream_mx(:,:) = 1
168	betamx(:,:) = drag_mx(:,:)
169	endwhere
170
171	where (drag_my(:,:).ge.beta_limgz)
172	mstream_my(:,:) = 0
173	betamy(:,:) = beta_limgz
174	drag_my(:,:) = beta_limgz
175	elsewhere
176	mstream_my(:,:) = 1
177	betamy(:,:) = drag_my(:,:)
178	endwhere
179
180	if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture')
181	mstream_mx(:,:) = 0
182	mstream_my(:,:) = 0
183	!
184
185	call gzm_beta_prescr
186
187	return
188
189	end subroutine init_dragging
190
191
192	!________________________________________________________________________________
193
194	!> SUBROUTINE: dragging
195	!! Defini la localisation des streams et le frottement basal
196	!<
197	!-------------------------------------------------------------------------
198	subroutine dragging ! defini la localisation des streams et le frottement basal
199
200	implicit none
201
202
203
204	if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging')
205
206	where (mstream_mx(:,:).eq.1)
207	betamx(:,:) = drag_mx(:,:)
208	elsewhere
209	betamx(:,:)= beta_limgz
210	end where
211
212	where (mstream_my(:,:).eq.1)
213	betamy(:,:) = drag_my(:,:)
214	elsewhere
215	betamy(:,:)= beta_limgz
216	end where
217
218	betamx(:,:)=drag_mx(:,:)
219	betamy(:,:)=drag_my(:,:)
220
221
222	! update la valeur de drag_mx et drag_my en fonction de la pression effective donnee par la hauteur au dessu de la flottaison
223	call update_buoy(nx,ny,drag_centre,H_buoy_init,betamx,betamy,beta_centre,beta_min)
224
225	debug_3D(:,:,63) = H_buoy_courant(:,:)
226	debug_3D(:,:,101) = drag_centre(:,:)
227	debug_3D(:,:,102) = H_buoy_init(:,:)
228
229
230
231	call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants
232
233
234	return
235	end subroutine dragging
236
237	!----------------------------------------------------------------------------------
238	!>SUBROUTINE: gzm_beta_prescr
239	!! Calcul de gzmx
240	!<
241
242	subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my)
243
244	logical :: test_Tx ! test if there is a basal melting point in the surrounding
245	logical :: test_Ty ! test if there is a basal melting point in the surrounding
246	logical :: test_beta_x ! test if there is a low drag point in the surrounding
247	logical :: test_beta_y ! test if there is a low drag point in the surrounding
248
249	real :: beta_forc_fleuve ! below this value -> ice stream
250
251	! beta_forc_fleuve = 5.e3
252	beta_forc_fleuve = beta_limgz
253
254	! calcul de gzmx
255
256
257	gzmx(:,:)=.true.
258	gzmy(:,:)=.true.
259	flgzmx(:,:)=.false.
260	flgzmy(:,:)=.false.
261
262
263	! points flottants : flgzmx mais pas gzmx
264	do j=2,ny
265	do i=2,nx
266	if (flot(i,j).and.(flot(i-1,j))) then
267	flotmx(i,j)=.true.
268	flgzmx(i,j)=.true.
269	gzmx(i,j)=.false.
270	betamx(i,j)=0.
271	end if
272	if (flot(i,j).and.(flot(i,j-1))) then
273	flotmy(i,j)=.true.
274	flgzmy(i,j)=.true.
275	gzmy(i,j)=.false.
276	betamy(i,j)=0.
277	end if
278	end do
279	end do
280
281	if (itracebug.eq.1) call tracebug(' apres criteres flot')
282
283	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
284
285	!--------- autres criteres
286
287	! points poses
288	gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points
289	gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta
290
291
292
293	! critere (lache) sur la temperature
294	do j = 2, ny-1
295	do i =2, nx-1
296
297	! test s'il y a un point tempere dans les environs
298	test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1)
299	test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1)
300
301	! test s'il y a un point low drag dans les environs
302	test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve )
303	test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve )
304
305	! critere pour gzmx
306
307	! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
308
309	! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x)
310	! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y)
311
312	end do
313	end do
314
315
316	if (itracebug.eq.1) call tracebug(' apres autres criteres ')
317	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
318
319	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
320	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
321
322	where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
323	betamx(:,:) = beta_limgz
324	end where
325
326	where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
327	betamy(:,:) = beta_limgz
328	end where
329
330	fleuvemx(:,:)=gzmx(:,:)
331	fleuvemy(:,:)=gzmy(:,:)
332
333
334	end subroutine gzm_beta_prescr
335
336
337	!----------------------------------------------------------------------------------
338	! Routines qui permettent des manip simples sur beta.
339
340
341
342	!______________________________________________________________________________________
343	!>SUBROUTINE: beta_min_value
344	!! valeur mini que peut avoir beta (en Pa an m-1)
345	!<
346
347	subroutine beta_min_value(nxx,nyy,R_mx,R_my,R_centre,val_betamin)
348	implicit none
349	integer :: nxx,nyy ! dimension des tableaux
350	real, dimension(nxx,nyy), intent(inout) :: R_mx ! valeur du beta sur maille mx
351	real, dimension(nxx,nyy), intent(inout) :: R_my ! valeur du beta sur maille my
352	real, dimension(nxx,nyy), intent(inout) :: R_centre ! valeur du beta sur maille centree
353	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
354
355	R_centre(:,:) = max(R_centre(:,:),val_betamin)
356	R_mx(:,:) = max(R_mx(:,:),val_betamin)
357	R_my(:,:) = max(R_my(:,:),val_betamin)
358
359	where(flot(:,:))
360	drag_centre(:,:) = 0.
361	end where
362
363	end subroutine beta_min_value
364	!______________________________________________________________________________________
365	!>SUBROUTINE: beta_stag2c
366	!! average the staggered values on central ones
367	!<
368
369	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
370
371	implicit none
372	integer :: nxx,nyy ! dimension des tableaux
373	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
374	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
375	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
376
377	R_centre(:,:)=0.
378	do j=2,ny-1
379	do i=2,nx-1
380	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
381	+ (R_my(i,j)+R_my(i,j+1)))*0.25
382	end do
383	end do
384	end subroutine beta_stag2c
385	!______________________________________________________________________________________
386	!>SUBROUTINE: beta_c2stag_min
387	!! redistribute central values on staggered grid
388	!<
389
390	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
391
392	implicit none
393	integer :: nxx,nyy ! dimension des tableaux
394	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
395	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
396	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
397
398	do j=2,nyy
399	do i=2,nxx
400	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
401	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
402	end do
403	end do
404
405	end subroutine beta_c2stag
406	!______________________________________________________________________________________
407	!>SUBROUTINE: beta_c2stag_min
408	!! redistribute central values on staggered grid
409	!<
410
411	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
412
413	implicit none
414	integer :: nxx,nyy ! dimension des tableaux
415	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
416	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
417	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
418
419	do j=2,nyy
420	do i=2,nxx
421	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
422	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
423	end do
424	end do
425
426	end subroutine beta_c2stag_min
427
428	!>SUBROUTINE: update_buoy
429	!! fait changer la valeur centree enfonction de la hauteur au dessus de la flottaison.
430	!< on pressuppose beta = beta0 * Hbuoy
431
432	subroutine update_buoy(nxx,nyy,Ref_centre,Ref_Hbuoy,R_mx,R_my,R_centre,val_betamin)
433	implicit none
434	integer :: nxx,nyy ! dimension des tableaux
435	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
436	real, dimension(nxx,nyy), intent(in) :: Ref_Hbuoy ! valeur de reference hauteur au dessus de la flottaison
437	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
438	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
439	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
440	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
441
442	call calc_buoyency(nxx,nyy,H,B,H_buoy_courant)
443
444
445
446	where (Ref_Hbuoy(:,:).gt.0.)
447	! R_centre(:,:) = Ref_centre(:,:)*max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)
448
449
450	! avec une dépendance non linéaire
451	R_centre(:,:) = Ref_centre(:,:)(max((max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)),0.01))*(1./3.)
452
453
454
455
456	elsewhere (Ref_Hbuoy(:,:).le.0.) ! flottant (0.) ou eventuellement pas de glace
457	R_centre(:,:) = Ref_centre(:,:)
458
459	end where
460
461	! valeur mini
462
463	R_centre(:,:) = max(R_centre(:,:),val_betamin)
464
465	! si flottaison drag = 0.
466
467	where(flot(:,:))
468	drag_centre(:,:) = 0.
469	end where
470
471	! remet sur les noeuds staggered
472
473	do j=2,nyy
474	do i=2,nxx
475	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
476	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
477	end do
478	end do
479
480	! valeur mini
481	R_mx(:,:) = max(R_mx(:,:),val_betamin)
482	R_my(:,:) = max(R_my(:,:),val_betamin)
483
484
485	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
486
487	end subroutine update_buoy
488
489
490	!>SUBROUTINE: calc_buoyency
491	!< calcule la hauteur au dessus de la flottaison.
492
493	subroutine calc_buoyency(nxx,nyy,H1,B1,H_buoy)
494
495	implicit none
496	integer :: nxx,nyy ! dimension des tableaux
497	real, dimension(nxx,nyy), intent(in) :: H1 ! epaisseur
498	real, dimension(nxx,nyy), intent(in) :: B1 ! base de la glace (attention : pas Bsoc)
499	real, dimension(nxx,nyy), intent(out) :: H_buoy
500
501
502	where ( sealevel-B1(:,:).le.0.) ! socle au dessus du niveau des mers
503	H_buoy(:,:) = H1(:,:)
504
505	elsewhere
506	H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
507
508	end where
509
510	end subroutine calc_buoyency
511
512
513	end module dragging_beta_buoy
514

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

Download in other formats: