Context Navigation

source: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_mod.f90 @ 150

Last change on this file since 150 was 93, checked in by dumas, 8 years ago
First version with Schoof flux parameterisation at the grounding line. \| New module furst_schoof_mod.f90 \| New flag Schoof in grdline namelist (see in SOURCES/Fichiers-parametres/A-LBq15_param_list_Schoof.dat)
File size: 24.2 KB

Line
1	!> \file dragging_prescr_beta_nolin_mod.f90
2	!< Module pour calculer le beta a partir d'une carte sur les noeuds centres.
3	!< Cette version utilise le beta inverse lineaire comme entree
4	!< mais ensuite suppose une loi non lineaire
5	!< copie sur la formulation basee sur la hauteur au dessus de la flottaison
6	!< dont on garde les variables, ....
7	!<
8
9	!> \namespace dragging_prescr_beta
10	!! Calcule le beta a partir de vitesses de bilan
11	!! @note Il faut partir d'un cptr
12	!! \author ...
13	!! \date ...
14	!! @note Used module
15	!! @note - use module3D_phy
16	!<
17
18	module dragging_beta_nolin
19
20
21	use module3d_phy
22	use interface_input
23	implicit none
24	real, dimension(nx,ny) :: Vcol_x !< uniquement pour compatibilite
25	real, dimension(nx,ny) :: Vcol_y !<
26	real, dimension(nx,ny) :: Vsl_x !<
27	real, dimension(nx,ny) :: Vsl_y !<
28
29
30	real, dimension(nx,ny) :: drag_centre !< beta on major node (average)
31	!< beta_centre idem est la valeur courante declare dans 3D
32
33	! real, dimension(nx,ny) :: coef_drag !< coefficient de la loi de friction non lineaire : depend de la valeur de alpha_drag
34	! passe dans 3D !< si alpha_drag = 1, coef_drag = drag_centre
35
36
37	real, dimension(nx,ny) :: H_buoy_init !< hauteur au dessus de la flottaison a initialisation
38	real, dimension(nx,ny) :: H_buoy_courant !< hauteur courante au dessus de la flottaison
39	real, dimension(nx,ny) :: Uslid_init !< vitesse de glissement sur noeud centre a initialisation
40	real, dimension(nx,ny) :: Uslid_centre !< vitesse de glissement sur noeud centre courant
41
42	real, dimension(nx,ny) :: betamoy_x !< pour faire les moyennes glissantes
43	real, dimension(nx,ny) :: betamoy_y !<
44	integer, dimension(nx,ny) :: nbvoisx !<
45	integer, dimension(nx,ny) :: nbvoisy !<
46
47
48	real :: beta_limgz !< when beta gt beta_limgz -> not gzmx
49	real :: beta_min !< minimum value of beta
50	real :: beta_mult !< coefficient of beta field
51	! real :: alpha_drag != 3 !< exposant non linearite, dans la namelist
52	! real :: exp_alpha !< coefficient dans la loi
53
54	integer :: ill,jll,nmoy
55
56	real :: maxi !< calcul correction dS a la grounding line
57	real :: mini
58
59	logical :: corr_def = .False. !< for deformation correction (compatibility)
60
61	! declares dans le 3D (variables globales)
62	!-------------------------------------------
63	! beta_centre : valeur centree courante
64	! beta_mx et beta_my sont les valeurs courantes
65	! drag_mx, drag_my servent dans calc_beta (remplimat-shelves-tabTu.f90) avec un autre sens
66	! et de variable de calcul ici.
67
68
69
70	contains
71
72	!-----------------------------------------------------------------------------------------
73	!> SUBROUTINE: init_dragging
74	!! Cette routine fait l'initialisation du dragging.
75	!<
76	subroutine init_dragging ! Cette routine fait l'initialisation du dragging.
77	! nouvelle version : lit les fichiers nc
78	implicit none
79	character(len=100) :: beta_c_file ! beta on centered grid
80	character(len=100) :: betamx_file ! beta on staggered grid mx
81	character(len=100) :: betamy_file ! beta on staggered grid mx
82
83	! pour la lecture dans le fichier parametre
84	character(len=4) :: char_num ! pour la lecture
85	integer :: lenrun ! longueur du run name
86	integer :: num_job ! numero du job
87	integer :: num ! pour la lecture
88	integer :: nlenght
89	character(len=200) :: file_exp ! nom du fichier experience
90	character(len=300) :: string_lect ! pour lire toute la ligne
91
92
93	namelist/beta_prescr_nolin/beta_c_file,beta_limgz,beta_min,beta_mult,alpha_drag
94
95	if (itracebug.eq.1) call tracebug(' Prescr_beta subroutine init_dragging')
96
97
98
99	! lecture des parametres du run
100	!--------------------------------------------------------------------
101
102	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
103	428 format(A)
104	read(num_param,beta_prescr_nolin)
105
106	write(num_rep_42,428) '!___________________________________________________________'
107	write(num_rep_42,428) '! read beta on centered grid'
108	write(num_rep_42,beta_prescr_nolin)
109	write(num_rep_42,428) '! beta_file : nom des fichier qui contiennent les beta_c'
110	write(num_rep_42,428) '! above beta_limgz, gzmx is false'
111	write(num_rep_42,428) '! if grounded, beta_min minimum value '
112	write(num_rep_42,428) '! beta_mult : coefficient just after reading'
113	write(num_rep_42,*)
114	write(num_rep_42,428) '!___________________________________________________________'
115
116
117	if (alpha_drag.lt.1.e-5) then ! mettre 0 dans le fichier parametre pour lire l'exposant ailleurs
118
119
120	! bloc specifique pour les experiences nolin de recul de ligne d'echouage janvier 2013------------
121
122	! recupere le numero du job dans le runname
123
124	lenrun = len(trim(runname))
125	char_num = runname(lenrun-3:lenrun)
126
127	read(char_num,*) num_job
128
129	! fichier d'experience : attention ! le fichier d'experience n'est lu qu'apres, ici c est juste la premiere et la derniere colonne
130	! file_exp = 'antarctic_stability_1000_samples-nolin.dat'
131	file_exp = 'antarctic_stab_1000_samples-inv_nolin.dat'
132	file_exp = trim(dirnameinp)//trim(file_exp)
133	write(6,*) file_exp
134
135	! Read the experiment list and keep the line corresponding to the job number
136
137	open(98,file=file_exp)
138	! read(98,*) ! lit la ligne de titre
139
140	do i=1,2000 ! each line is a parameter set
141	read(98,201,end=200) string_lect
142	read(string_lect,*) num
143	! write(6,*) num
144	! write(6,*) string_lect
145	201 format(A300)
146
147	if (num.eq.num_job) then ! recherche de la derniere colonne
148	! write(6,*) string_lect
149	! nlenght=len_trim(adjustr(string_lect))
150	nlenght=len_trim(string_lect)
151	string_lect = string_lect(nlenght-4:nlenght)
152	! write(6,*) nlenght, string_lect
153	read(string_lect,*) alpha_drag
154	write(6,*) 'run',num,'exposant du frottement', alpha_drag
155	write(42,*) '! exposant du frottement =', alpha_drag
156	exit ! stop at the job number
157	end if
158
159	end do
160	200 continue
161	close(98)
162	end if
163	!--------------------------------------- fin du bloc ---------------------------------------------------------------------
164	iter_beta=0
165	exp_alpha = 1/alpha_drag - 1.
166
167	beta_c_file = trim(dirnameinp)//trim(beta_c_file)
168	betamx_file = trim(dirnameinp)//trim(betamx_file)
169	betamy_file = trim(dirnameinp)//trim(betamy_file)
170	betamax = beta_limgz
171	betamax_2d(:,:)=betamax
172
173	! read the beta value on centered and staggered grids
174	! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
175	! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
176	! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
177
178
179	call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
180
181	! calcule la valeur de reference de Hbuoy la hauteur au dessus de la flottaison
182	call calc_buoyency(nx,ny,H,B,H_buoy_init)
183
184
185	! calcule la valeur de reference de Uslid
186	call calc_Uslid_maj(nx,ny,ux(:,:,nz),uy(:,:,nz),Uslid_init)
187
188	! multiplication du frottement par le coefficient beta_mult
189
190	drag_centre(:,:) = drag_centre(:,:)*beta_mult
191
192	where (mk_init(:,:).eq.-2) ! iles a probleme
193	drag_centre(:,:) = 2.* abs(beta_limgz)
194	end where
195
196
197	where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes
198	drag_centre(:,:) = 1000. ! valeur relativement faible qui evite
199	end where ! les freinages
200
201
202	call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
203
204	if (beta_limgz.gt.0.) then
205	beta_centre(:,:) = drag_centre(:,:)
206	betamx(:,:) = drag_mx(:,:)
207	betamy(:,:) = drag_my(:,:)
208
209	else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic
210
211	beta_centre(:,:) = - beta_limgz
212	betamx(:,:) = - beta_limgz
213	betamy(:,:) = - beta_limgz
214	drag_centre(:,:) = - beta_limgz
215	drag_mx(:,:) = - beta_limgz
216	drag_my(:,:) = - beta_limgz
217	beta_limgz = 1.
218	end if
219
220
221
222	! valeur mini que peut avoir beta (en Pa an m-1)
223	call beta_min_value(nx,ny,drag_mx,drag_my,drag_centre,beta_min)
224
225	! on peut envisager une valeur ~ 10
226	! rappel : beta doit être positif.
227
228
229
230	beta_centre(:,:) = drag_centre(:,:) ! utilise aussi dans schoof
231
232	! initialisation de coef_drag. depend de l'exposant de la loi de dragging
233
234	where (abs(Uslid_init(:,:)).gt.1.)
235	coef_drag(:,:)= drag_centre(:,:)(1./abs(Uslid_init(:,:)))*exp_alpha
236	elsewhere
237	coef_drag(:,:)= drag_centre(:,:)
238	end where
239	debug_3D(:,:,105) = coef_drag(:,:)
240
241
242	where (drag_mx(:,:).ge.beta_limgz)
243	mstream_mx(:,:) = 0
244	betamx(:,:) = beta_limgz
245	drag_mx(:,:) = beta_limgz
246	elsewhere
247	mstream_mx(:,:) = 1
248	betamx(:,:) = drag_mx(:,:)
249	endwhere
250
251	where (drag_my(:,:).ge.beta_limgz)
252	mstream_my(:,:) = 0
253	betamy(:,:) = beta_limgz
254	drag_my(:,:) = beta_limgz
255	elsewhere
256	mstream_my(:,:) = 1
257	betamy(:,:) = drag_my(:,:)
258	endwhere
259
260	if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture')
261	mstream_mx(:,:) = 0
262	mstream_my(:,:) = 0
263	!
264
265	call gzm_beta_prescr
266
267	return
268
269	end subroutine init_dragging
270
271
272	!________________________________________________________________________________
273
274	!> SUBROUTINE: dragging
275	!! Defini la localisation des streams et le frottement basal
276	!<
277	!-------------------------------------------------------------------------
278	subroutine dragging ! defini la localisation des streams et le frottement basal
279
280	implicit none
281
282	real :: som1 ! pour calcul test_iter_diag
283	real :: som2
284
285
286
287	if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging')
288
289	where (mstream_mx(:,:).eq.1)
290	betamx(:,:) = drag_mx(:,:)
291	elsewhere
292	betamx(:,:)= beta_limgz
293	end where
294
295	where (mstream_my(:,:).eq.1)
296	betamy(:,:) = drag_my(:,:)
297	elsewhere
298	betamy(:,:)= beta_limgz
299	end where
300
301	betamx(:,:)=drag_mx(:,:)
302	betamy(:,:)=drag_my(:,:)
303
304
305	! update la valeur de drag_mx et drag_my en fonction de la pression effective donnee par la hauteur au dessu de la flottaison
306	! call update_buoy(nx,ny,drag_centre,H_buoy_init,betamx,betamy,beta_centre,beta_min)
307
308	! update la valeur de beta_mx, beta_my, ... en fonction de la vitesse basale
309
310	call update_nolin(nx,ny,drag_centre,Uslid_init,betamx,betamy,beta_centre,beta_min)
311
312	debug_3D(:,:,63) = H_buoy_courant(:,:)
313	debug_3D(:,:,101) = drag_centre(:,:)
314	debug_3D(:,:,102) = H_buoy_init(:,:)
315	debug_3D(:,:,103) = Uslid_init(:,:)
316	debug_3D(:,:,104) = Uslid_centre(:,:)
317
318	! calcul du critere de convergence
319	som1 = 0.
320	som2 = 0.
321
322	do j=2,ny
323	do i = 2,nx
324	if (.not.flot(i,j)) then
325	som1 = (Uslid_centre(i,j)-Uiter_centre(i,j))**2+som1
326	som2 = Uslid_centre(i,j)**2+som2
327	end if
328	end do
329	end do
330	Uiter_centre(:,:) = Uslid_centre(:,:)
331	test_iter_diag = som1/som2
332
333
334
335
336
337	call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants
338
339
340	return
341	end subroutine dragging
342
343	!----------------------------------------------------------------------------------
344	!>SUBROUTINE: gzm_beta_prescr
345	!! Calcul de gzmx
346	!<
347
348	subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my)
349
350	logical :: test_Tx ! test if there is a basal melting point in the surrounding
351	logical :: test_Ty ! test if there is a basal melting point in the surrounding
352	logical :: test_beta_x ! test if there is a low drag point in the surrounding
353	logical :: test_beta_y ! test if there is a low drag point in the surrounding
354
355	real :: beta_forc_fleuve ! below this value -> ice stream
356
357	! beta_forc_fleuve = 5.e3
358	beta_forc_fleuve = beta_limgz
359
360	! calcul de gzmx
361
362
363	gzmx(:,:)=.true.
364	gzmy(:,:)=.true.
365	flgzmx(:,:)=.false.
366	flgzmy(:,:)=.false.
367
368
369	! points flottants : flgzmx mais pas gzmx
370	do j=2,ny
371	do i=2,nx
372	if (flot(i,j).and.(flot(i-1,j))) then
373	flotmx(i,j)=.true.
374	flgzmx(i,j)=.true.
375	gzmx(i,j)=.false.
376	betamx(i,j)=0.
377	end if
378	if (flot(i,j).and.(flot(i,j-1))) then
379	flotmy(i,j)=.true.
380	flgzmy(i,j)=.true.
381	gzmy(i,j)=.false.
382	betamy(i,j)=0.
383	end if
384	end do
385	end do
386
387	if (itracebug.eq.1) call tracebug(' apres criteres flot')
388
389	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
390
391	!--------- autres criteres
392
393	! points poses
394	gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points
395	gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta
396
397
398
399	! critere (lache) sur la temperature
400	do j = 2, ny-1
401	do i =2, nx-1
402
403	! test s'il y a un point tempere dans les environs
404	test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1)
405	test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1)
406
407	! test s'il y a un point low drag dans les environs
408	test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve )
409	test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve )
410
411	! critere pour gzmx
412
413	! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
414
415	! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x)
416	! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y)
417
418	end do
419	end do
420
421
422	if (itracebug.eq.1) call tracebug(' apres autres criteres ')
423	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
424
425	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
426	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
427
428	where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
429	betamx(:,:) = beta_limgz
430	end where
431
432	where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
433	betamy(:,:) = beta_limgz
434	end where
435
436	fleuvemx(:,:)=gzmx(:,:)
437	fleuvemy(:,:)=gzmy(:,:)
438
439	end subroutine gzm_beta_prescr
440
441
442	!----------------------------------------------------------------------------------
443	! Routines qui permettent des manip simples sur beta.
444
445
446
447	!______________________________________________________________________________________
448	!>SUBROUTINE: beta_min_value
449	!! valeur mini que peut avoir beta (en Pa an m-1)
450	!<
451
452	subroutine beta_min_value(nxx,nyy,R_mx,R_my,R_centre,val_betamin)
453	implicit none
454	integer :: nxx,nyy ! dimension des tableaux
455	real, dimension(nxx,nyy), intent(inout) :: R_mx ! valeur du beta sur maille mx
456	real, dimension(nxx,nyy), intent(inout) :: R_my ! valeur du beta sur maille my
457	real, dimension(nxx,nyy), intent(inout) :: R_centre ! valeur du beta sur maille centree
458	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
459
460	R_centre(:,:) = max(R_centre(:,:),val_betamin)
461	R_mx(:,:) = max(R_mx(:,:),val_betamin)
462	R_my(:,:) = max(R_my(:,:),val_betamin)
463
464	where(flot(:,:))
465	drag_centre(:,:) = 0.
466	end where
467
468	end subroutine beta_min_value
469	!______________________________________________________________________________________
470	!>SUBROUTINE: beta_stag2c
471	!! average the staggered values on central ones
472	!<
473
474	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
475
476	implicit none
477	integer :: nxx,nyy ! dimension des tableaux
478	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
479	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
480	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
481
482	R_centre(:,:)=0.
483	do j=2,ny-1
484	do i=2,nx-1
485	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
486	+ (R_my(i,j)+R_my(i,j+1)))*0.25
487	end do
488	end do
489	end subroutine beta_stag2c
490	!______________________________________________________________________________________
491	!>SUBROUTINE: beta_c2stag_min
492	!! redistribute central values on staggered grid
493	!<
494
495	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
496
497	implicit none
498	integer :: nxx,nyy ! dimension des tableaux
499	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
500	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
501	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
502
503	do j=2,nyy
504	do i=2,nxx
505	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
506	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
507	end do
508	end do
509
510	end subroutine beta_c2stag
511	!______________________________________________________________________________________
512	!>SUBROUTINE: beta_c2stag_min
513	!! redistribute central values on staggered grid
514	!<
515
516	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
517
518	implicit none
519	integer :: nxx,nyy ! dimension des tableaux
520	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
521	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
522	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
523
524	do j=2,nyy
525	do i=2,nxx
526	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
527	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
528	end do
529	end do
530
531	end subroutine beta_c2stag_min
532
533	!>SUBROUTINE: update_buoy
534	!! fait changer la valeur centree enfonction de la hauteur au dessus de la flottaison.
535	!< on pressuppose beta = beta0 * Hbuoy
536
537	subroutine update_buoy(nxx,nyy,Ref_centre,Ref_Hbuoy,R_mx,R_my,R_centre,val_betamin)
538	implicit none
539	integer :: nxx,nyy ! dimension des tableaux
540	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
541	real, dimension(nxx,nyy), intent(in) :: Ref_Hbuoy ! valeur de reference hauteur au dessus de la flottaison
542	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
543	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
544	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
545	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
546
547	call calc_buoyency(nxx,nyy,H,B,H_buoy_courant)
548
549
550
551	where (Ref_Hbuoy(:,:).gt.0.)
552	! R_centre(:,:) = Ref_centre(:,:)*max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)
553
554
555	! avec une dépendance non linéaire
556	R_centre(:,:) = Ref_centre(:,:)(max((max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)),0.01))*(1./3.)
557
558
559
560
561	elsewhere (Ref_Hbuoy(:,:).le.0.) ! flottant (0.) ou eventuellement pas de glace
562	R_centre(:,:) = Ref_centre(:,:)
563
564	end where
565
566	! valeur mini
567
568	R_centre(:,:) = max(R_centre(:,:),val_betamin)
569
570	! si flottaison drag = 0.
571
572	where(flot(:,:))
573	drag_centre(:,:) = 0.
574	end where
575
576	! remet sur les noeuds staggered
577
578	do j=2,nyy
579	do i=2,nxx
580	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
581	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
582	end do
583	end do
584
585	! valeur mini
586	R_mx(:,:) = max(R_mx(:,:),val_betamin)
587	R_my(:,:) = max(R_my(:,:),val_betamin)
588
589
590	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
591
592	end subroutine update_buoy
593
594
595	!>SUBROUTINE: calc_buoyency
596	!< calcule la hauteur au dessus de la flottaison.
597
598	subroutine calc_buoyency(nxx,nyy,H1,B1,H_buoy)
599
600	implicit none
601	integer :: nxx,nyy ! dimension des tableaux
602	real, dimension(nxx,nyy), intent(in) :: H1 ! epaisseur
603	real, dimension(nxx,nyy), intent(in) :: B1 ! base de la glace (attention : pas Bsoc)
604	real, dimension(nxx,nyy), intent(out) :: H_buoy
605
606
607	where ( sealevel-B1(:,:).le.0.) ! socle au dessus du niveau des mers
608	H_buoy(:,:) = H1(:,:)
609
610	elsewhere
611	H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
612
613	end where
614
615	end subroutine calc_buoyency
616
617	!>SUBROUTINE: calc_buoyency
618	!< calcule la hauteur au dessus de la flottaison.
619
620	subroutine calc_Uslid_maj(nxx,nyy,ux_stag,uy_stag,u_centre)
621
622	implicit none
623	integer :: nxx,nyy ! dimension des tableaux
624	real, dimension(nxx,nyy), intent(in) :: ux_stag ! vitesse selon ux
625	real, dimension(nxx,nyy), intent(in) :: uy_stag ! vitesse selon uy
626	real, dimension(nxx,nyy), intent(out) :: u_centre ! vitesse sur noeud majeur
627
628	do j=2,nyy-1
629	do i=2,nxx-1
630	u_centre(i,j) = 0.5* ((ux_stag(i,j)+ux_stag(i+1,j))2 +(uy_stag(i,j)+uy_stag(i,j+1))2)**0.5
631	end do
632	end do
633
634	end subroutine calc_Uslid_maj
635
636	subroutine update_nolin(nxx,nyy,Ref_centre,Ref_Uslid,R_mx,R_my,R_centre,val_betamin)
637
638	implicit none
639	integer :: nxx,nyy ! dimension des tableaux
640	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
641	real, dimension(nxx,nyy), intent(in) :: Ref_Uslid ! valeur de reference vitesse de glissement
642	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
643	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
644	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
645	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
646
647	call calc_Uslid_maj(nxx,nyy,ux(:,:,nz),uy(:,:,nz),Uslid_centre)
648
649	where ((Ref_Uslid(:,:).gt.1.).and.(Uslid_centre(:,:).gt.1.))
650	R_centre(:,:) = Ref_centre(:,:)(Uslid_centre(:,:)/Ref_Uslid(:,:))*exp_alpha
651	! coef_drag(:,:)= Ref_centre(:,:)(1./Ref_Uslid(:,:))*exp_alpha
652	elsewhere
653	R_centre(:,:) = Ref_centre(:,:)
654	! coef_drag(:,:)= Ref_centre(:,:)
655	end where
656
657	! valeur mini
658	R_centre(:,:) = max(R_centre(:,:),val_betamin)
659
660	! si flottaison drag = 0.
661
662	where(flot(:,:))
663	R_centre(:,:) = 0.
664	end where
665
666	! remet sur les noeuds staggered
667
668	do j=2,nyy
669	do i=2,nxx
670	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
671	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
672	end do
673	end do
674
675	! valeur mini
676	R_mx(:,:) = max(R_mx(:,:),val_betamin)
677	R_my(:,:) = max(R_my(:,:),val_betamin)
678
679
680	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
681
682	end subroutine update_nolin
683
684
685	end module dragging_beta_nolin
686

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

Download in other formats: