Context Navigation

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

Last change on this file since 4 was 4, checked in by dumas, 10 years ago
initial import GRISLI trunk
File size: 17.3 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
103	! read the beta value on centered and staggered grids
104	! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
105	! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
106	! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
107
108
109	call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
110
111	! calcule la valeur de reference de Hbuoy la hauteur au dessus de la flottaison
112	call calc_buoyency(nx,ny,H,B,H_buoy_init)
113
114
115	! multiplication du frottement par le coefficient beta_mult
116
117	drag_centre(:,:) = drag_centre(:,:)*beta_mult
118
119	where (mk_init(:,:).eq.-2) ! iles a probleme
120	drag_centre(:,:) = 2.* abs(beta_limgz)
121	end where
122
123
124	where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes
125	drag_centre(:,:) = 1000. ! valeur relativement faible qui evite
126	end where ! les freinages
127
128
129	call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
130
131	if (beta_limgz.gt.0.) then
132	beta_centre(:,:) = drag_centre(:,:)
133	betamx(:,:) = drag_mx(:,:)
134	betamy(:,:) = drag_my(:,:)
135
136	else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic
137
138	beta_centre(:,:) = - beta_limgz
139	betamx(:,:) = - beta_limgz
140	betamy(:,:) = - beta_limgz
141	drag_centre(:,:) = - beta_limgz
142	drag_mx(:,:) = - beta_limgz
143	drag_my(:,:) = - beta_limgz
144	beta_limgz = 1.
145	end if
146
147
148
149	! valeur mini que peut avoir beta (en Pa an m-1)
150	call beta_min_value(nx,ny,drag_mx,drag_my,drag_centre,beta_min)
151
152	! on peut envisager une valeur ~ 10
153	! rappel : beta doit être positif.
154
155
156
157	beta_centre(:,:) = drag_centre(:,:) ! utilise aussi dans schoof
158
159
160
161	where (drag_mx(:,:).ge.beta_limgz)
162	mstream_mx(:,:) = 0
163	betamx(:,:) = beta_limgz
164	drag_mx(:,:) = beta_limgz
165	elsewhere
166	mstream_mx(:,:) = 1
167	betamx(:,:) = drag_mx(:,:)
168	endwhere
169
170	where (drag_my(:,:).ge.beta_limgz)
171	mstream_my(:,:) = 0
172	betamy(:,:) = beta_limgz
173	drag_my(:,:) = beta_limgz
174	elsewhere
175	mstream_my(:,:) = 1
176	betamy(:,:) = drag_my(:,:)
177	endwhere
178
179	if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture')
180	mstream_mx(:,:) = 0
181	mstream_my(:,:) = 0
182	!
183
184	call gzm_beta_prescr
185
186	return
187
188	end subroutine init_dragging
189
190
191	!________________________________________________________________________________
192
193	!> SUBROUTINE: dragging
194	!! Defini la localisation des streams et le frottement basal
195	!<
196	!-------------------------------------------------------------------------
197	subroutine dragging ! defini la localisation des streams et le frottement basal
198
199	implicit none
200
201
202
203	if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging')
204
205	where (mstream_mx(:,:).eq.1)
206	betamx(:,:) = drag_mx(:,:)
207	elsewhere
208	betamx(:,:)= beta_limgz
209	end where
210
211	where (mstream_my(:,:).eq.1)
212	betamy(:,:) = drag_my(:,:)
213	elsewhere
214	betamy(:,:)= beta_limgz
215	end where
216
217	betamx(:,:)=drag_mx(:,:)
218	betamy(:,:)=drag_my(:,:)
219
220
221	! update la valeur de drag_mx et drag_my en fonction de la pression effective donnee par la hauteur au dessu de la flottaison
222	call update_buoy(nx,ny,drag_centre,H_buoy_init,betamx,betamy,beta_centre,beta_min)
223
224	debug_3D(:,:,63) = H_buoy_courant(:,:)
225	debug_3D(:,:,101) = drag_centre(:,:)
226	debug_3D(:,:,102) = H_buoy_init(:,:)
227
228
229
230	call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants
231
232
233	return
234	end subroutine dragging
235
236	!----------------------------------------------------------------------------------
237	!>SUBROUTINE: gzm_beta_prescr
238	!! Calcul de gzmx
239	!<
240
241	subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my)
242
243	logical :: test_Tx ! test if there is a basal melting point in the surrounding
244	logical :: test_Ty ! test if there is a basal melting point in the surrounding
245	logical :: test_beta_x ! test if there is a low drag point in the surrounding
246	logical :: test_beta_y ! test if there is a low drag point in the surrounding
247
248	real :: beta_forc_fleuve ! below this value -> ice stream
249
250	! beta_forc_fleuve = 5.e3
251	beta_forc_fleuve = beta_limgz
252
253	! calcul de gzmx
254
255
256	gzmx(:,:)=.true.
257	gzmy(:,:)=.true.
258	flgzmx(:,:)=.false.
259	flgzmy(:,:)=.false.
260
261
262	! points flottants : flgzmx mais pas gzmx
263	do j=2,ny
264	do i=2,nx
265	if (flot(i,j).and.(flot(i-1,j))) then
266	flotmx(i,j)=.true.
267	flgzmx(i,j)=.true.
268	gzmx(i,j)=.false.
269	betamx(i,j)=0.
270	end if
271	if (flot(i,j).and.(flot(i,j-1))) then
272	flotmy(i,j)=.true.
273	flgzmy(i,j)=.true.
274	gzmy(i,j)=.false.
275	betamy(i,j)=0.
276	end if
277	end do
278	end do
279
280	if (itracebug.eq.1) call tracebug(' apres criteres flot')
281
282	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
283
284	!--------- autres criteres
285
286	! points poses
287	gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points
288	gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta
289
290
291
292	! critere (lache) sur la temperature
293	do j = 2, ny-1
294	do i =2, nx-1
295
296	! test s'il y a un point tempere dans les environs
297	test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1)
298	test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1)
299
300	! test s'il y a un point low drag dans les environs
301	test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve )
302	test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve )
303
304	! critere pour gzmx
305
306	! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
307
308	! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x)
309	! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y)
310
311	end do
312	end do
313
314
315	if (itracebug.eq.1) call tracebug(' apres autres criteres ')
316	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
317
318	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
319	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
320
321	where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
322	betamx(:,:) = beta_limgz
323	end where
324
325	where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
326	betamy(:,:) = beta_limgz
327	end where
328
329
330
331	end subroutine gzm_beta_prescr
332
333
334	!----------------------------------------------------------------------------------
335	! Routines qui permettent des manip simples sur beta.
336
337
338
339	!______________________________________________________________________________________
340	!>SUBROUTINE: beta_min_value
341	!! valeur mini que peut avoir beta (en Pa an m-1)
342	!<
343
344	subroutine beta_min_value(nxx,nyy,R_mx,R_my,R_centre,val_betamin)
345	implicit none
346	integer :: nxx,nyy ! dimension des tableaux
347	real, dimension(nxx,nyy), intent(inout) :: R_mx ! valeur du beta sur maille mx
348	real, dimension(nxx,nyy), intent(inout) :: R_my ! valeur du beta sur maille my
349	real, dimension(nxx,nyy), intent(inout) :: R_centre ! valeur du beta sur maille centree
350	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
351
352	R_centre(:,:) = max(R_centre(:,:),val_betamin)
353	R_mx(:,:) = max(R_mx(:,:),val_betamin)
354	R_my(:,:) = max(R_my(:,:),val_betamin)
355
356	where(flot(:,:))
357	drag_centre(:,:) = 0.
358	end where
359
360	end subroutine beta_min_value
361	!______________________________________________________________________________________
362	!>SUBROUTINE: beta_stag2c
363	!! average the staggered values on central ones
364	!<
365
366	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
367
368	implicit none
369	integer :: nxx,nyy ! dimension des tableaux
370	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
371	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
372	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
373
374	R_centre(:,:)=0.
375	do j=2,ny-1
376	do i=2,nx-1
377	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
378	+ (R_my(i,j)+R_my(i,j+1)))*0.25
379	end do
380	end do
381	end subroutine beta_stag2c
382	!______________________________________________________________________________________
383	!>SUBROUTINE: beta_c2stag_min
384	!! redistribute central values on staggered grid
385	!<
386
387	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
388
389	implicit none
390	integer :: nxx,nyy ! dimension des tableaux
391	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
392	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
393	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
394
395	do j=2,nyy
396	do i=2,nxx
397	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
398	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
399	end do
400	end do
401
402	end subroutine beta_c2stag
403	!______________________________________________________________________________________
404	!>SUBROUTINE: beta_c2stag_min
405	!! redistribute central values on staggered grid
406	!<
407
408	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
409
410	implicit none
411	integer :: nxx,nyy ! dimension des tableaux
412	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
413	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
414	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
415
416	do j=2,nyy
417	do i=2,nxx
418	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
419	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
420	end do
421	end do
422
423	end subroutine beta_c2stag_min
424
425	!>SUBROUTINE: update_buoy
426	!! fait changer la valeur centree enfonction de la hauteur au dessus de la flottaison.
427	!< on pressuppose beta = beta0 * Hbuoy
428
429	subroutine update_buoy(nxx,nyy,Ref_centre,Ref_Hbuoy,R_mx,R_my,R_centre,val_betamin)
430	implicit none
431	integer :: nxx,nyy ! dimension des tableaux
432	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
433	real, dimension(nxx,nyy), intent(in) :: Ref_Hbuoy ! valeur de reference hauteur au dessus de la flottaison
434	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
435	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
436	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
437	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
438
439	call calc_buoyency(nxx,nyy,H,B,H_buoy_courant)
440
441
442
443	where (Ref_Hbuoy(:,:).gt.0.)
444	! R_centre(:,:) = Ref_centre(:,:)*max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)
445
446
447	! avec une dépendance non linéaire
448	R_centre(:,:) = Ref_centre(:,:)(max((max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)),0.01))*(1./3.)
449
450
451
452
453	elsewhere (Ref_Hbuoy(:,:).le.0.) ! flottant (0.) ou eventuellement pas de glace
454	R_centre(:,:) = Ref_centre(:,:)
455
456	end where
457
458	! valeur mini
459
460	R_centre(:,:) = max(R_centre(:,:),val_betamin)
461
462	! si flottaison drag = 0.
463
464	where(flot(:,:))
465	drag_centre(:,:) = 0.
466	end where
467
468	! remet sur les noeuds staggered
469
470	do j=2,nyy
471	do i=2,nxx
472	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
473	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
474	end do
475	end do
476
477	! valeur mini
478	R_mx(:,:) = max(R_mx(:,:),val_betamin)
479	R_my(:,:) = max(R_my(:,:),val_betamin)
480
481
482	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
483
484	end subroutine update_buoy
485
486
487	!>SUBROUTINE: calc_buoyency
488	!< calcule la hauteur au dessus de la flottaison.
489
490	subroutine calc_buoyency(nxx,nyy,H1,B1,H_buoy)
491
492	implicit none
493	integer :: nxx,nyy ! dimension des tableaux
494	real, dimension(nxx,nyy), intent(in) :: H1 ! epaisseur
495	real, dimension(nxx,nyy), intent(in) :: B1 ! base de la glace (attention : pas Bsoc)
496	real, dimension(nxx,nyy), intent(out) :: H_buoy
497
498
499	where ( sealevel-B1(:,:).le.0.) ! socle au dessus du niveau des mers
500	H_buoy(:,:) = H1(:,:)
501
502	elsewhere
503	H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
504
505	end where
506
507	end subroutine calc_buoyency
508
509
510	end module dragging_beta_buoy
511

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

Download in other formats: