Context Navigation

source: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_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: 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	betamax2d(:,:) = betamax
172
173
174	! read the beta value on centered and staggered grids
175	! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
176	! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
177	! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
178
179
180	call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
181
182	! calcule la valeur de reference de Hbuoy la hauteur au dessus de la flottaison
183	call calc_buoyency(nx,ny,H,B,H_buoy_init)
184
185
186	! calcule la valeur de reference de Uslid
187	call calc_Uslid_maj(nx,ny,ux(:,:,nz),uy(:,:,nz),Uslid_init)
188
189	! multiplication du frottement par le coefficient beta_mult
190
191	drag_centre(:,:) = drag_centre(:,:)*beta_mult
192
193	where (mk_init(:,:).eq.-2) ! iles a probleme
194	drag_centre(:,:) = 2.* abs(beta_limgz)
195	end where
196
197
198	where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes
199	drag_centre(:,:) = 1000. ! valeur relativement faible qui evite
200	end where ! les freinages
201
202
203	call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid
204
205	if (beta_limgz.gt.0.) then
206	beta_centre(:,:) = drag_centre(:,:)
207	betamx(:,:) = drag_mx(:,:)
208	betamy(:,:) = drag_my(:,:)
209
210	else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic
211
212	beta_centre(:,:) = - beta_limgz
213	betamx(:,:) = - beta_limgz
214	betamy(:,:) = - beta_limgz
215	drag_centre(:,:) = - beta_limgz
216	drag_mx(:,:) = - beta_limgz
217	drag_my(:,:) = - beta_limgz
218	beta_limgz = 1.
219	end if
220
221
222
223	! valeur mini que peut avoir beta (en Pa an m-1)
224	call beta_min_value(nx,ny,drag_mx,drag_my,drag_centre,beta_min)
225
226	! on peut envisager une valeur ~ 10
227	! rappel : beta doit être positif.
228
229
230
231	beta_centre(:,:) = drag_centre(:,:) ! utilise aussi dans schoof
232
233	! initialisation de coef_drag. depend de l'exposant de la loi de dragging
234
235	where (abs(Uslid_init(:,:)).gt.1.)
236	coef_drag(:,:)= drag_centre(:,:)(1./abs(Uslid_init(:,:)))*exp_alpha
237	elsewhere
238	coef_drag(:,:)= drag_centre(:,:)
239	end where
240	debug_3D(:,:,105) = coef_drag(:,:)
241
242
243	where (drag_mx(:,:).ge.beta_limgz)
244	mstream_mx(:,:) = 0
245	betamx(:,:) = beta_limgz
246	drag_mx(:,:) = beta_limgz
247	elsewhere
248	mstream_mx(:,:) = 1
249	betamx(:,:) = drag_mx(:,:)
250	endwhere
251
252	where (drag_my(:,:).ge.beta_limgz)
253	mstream_my(:,:) = 0
254	betamy(:,:) = beta_limgz
255	drag_my(:,:) = beta_limgz
256	elsewhere
257	mstream_my(:,:) = 1
258	betamy(:,:) = drag_my(:,:)
259	endwhere
260
261	if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture')
262	mstream_mx(:,:) = 0
263	mstream_my(:,:) = 0
264	!
265
266	call gzm_beta_prescr
267
268	return
269
270	end subroutine init_dragging
271
272
273	!________________________________________________________________________________
274
275	!> SUBROUTINE: dragging
276	!! Defini la localisation des streams et le frottement basal
277	!<
278	!-------------------------------------------------------------------------
279	subroutine dragging ! defini la localisation des streams et le frottement basal
280
281	implicit none
282
283	real :: som1 ! pour calcul test_iter_diag
284	real :: som2
285
286
287
288	if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging')
289
290	where (mstream_mx(:,:).eq.1)
291	betamx(:,:) = drag_mx(:,:)
292	elsewhere
293	betamx(:,:)= beta_limgz
294	end where
295
296	where (mstream_my(:,:).eq.1)
297	betamy(:,:) = drag_my(:,:)
298	elsewhere
299	betamy(:,:)= beta_limgz
300	end where
301
302	betamx(:,:)=drag_mx(:,:)
303	betamy(:,:)=drag_my(:,:)
304
305
306	! update la valeur de drag_mx et drag_my en fonction de la pression effective donnee par la hauteur au dessu de la flottaison
307	! call update_buoy(nx,ny,drag_centre,H_buoy_init,betamx,betamy,beta_centre,beta_min)
308
309	! update la valeur de beta_mx, beta_my, ... en fonction de la vitesse basale
310
311	call update_nolin(nx,ny,drag_centre,Uslid_init,betamx,betamy,beta_centre,beta_min)
312
313	debug_3D(:,:,63) = H_buoy_courant(:,:)
314	debug_3D(:,:,101) = drag_centre(:,:)
315	debug_3D(:,:,102) = H_buoy_init(:,:)
316	debug_3D(:,:,103) = Uslid_init(:,:)
317	debug_3D(:,:,104) = Uslid_centre(:,:)
318
319	! calcul du critere de convergence
320	som1 = 0.
321	som2 = 0.
322
323	do j=2,ny
324	do i = 2,nx
325	if (.not.flot(i,j)) then
326	som1 = (Uslid_centre(i,j)-Uiter_centre(i,j))**2+som1
327	som2 = Uslid_centre(i,j)**2+som2
328	end if
329	end do
330	end do
331	Uiter_centre(:,:) = Uslid_centre(:,:)
332	test_iter_diag = som1/som2
333
334
335
336
337
338	call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants
339
340
341	return
342	end subroutine dragging
343
344	!----------------------------------------------------------------------------------
345	!>SUBROUTINE: gzm_beta_prescr
346	!! Calcul de gzmx
347	!<
348
349	subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my)
350
351	logical :: test_Tx ! test if there is a basal melting point in the surrounding
352	logical :: test_Ty ! test if there is a basal melting point in the surrounding
353	logical :: test_beta_x ! test if there is a low drag point in the surrounding
354	logical :: test_beta_y ! test if there is a low drag point in the surrounding
355
356	real :: beta_forc_fleuve ! below this value -> ice stream
357
358	! beta_forc_fleuve = 5.e3
359	beta_forc_fleuve = beta_limgz
360
361	! calcul de gzmx
362
363
364	gzmx(:,:)=.true.
365	gzmy(:,:)=.true.
366	flgzmx(:,:)=.false.
367	flgzmy(:,:)=.false.
368
369
370	! points flottants : flgzmx mais pas gzmx
371	do j=2,ny
372	do i=2,nx
373	if (flot(i,j).and.(flot(i-1,j))) then
374	flotmx(i,j)=.true.
375	flgzmx(i,j)=.true.
376	gzmx(i,j)=.false.
377	betamx(i,j)=0.
378	end if
379	if (flot(i,j).and.(flot(i,j-1))) then
380	flotmy(i,j)=.true.
381	flgzmy(i,j)=.true.
382	gzmy(i,j)=.false.
383	betamy(i,j)=0.
384	end if
385	end do
386	end do
387
388	if (itracebug.eq.1) call tracebug(' apres criteres flot')
389
390	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
391
392	!--------- autres criteres
393
394	! points poses
395	gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points
396	gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta
397
398
399
400	! critere (lache) sur la temperature
401	do j = 2, ny-1
402	do i =2, nx-1
403
404	! test s'il y a un point tempere dans les environs
405	test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1)
406	test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1)
407
408	! test s'il y a un point low drag dans les environs
409	test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve )
410	test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve )
411
412	! critere pour gzmx
413
414	! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
415
416	! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x)
417	! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y)
418
419	end do
420	end do
421
422
423	if (itracebug.eq.1) call tracebug(' apres autres criteres ')
424	if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)
425
426	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
427	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
428
429	where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
430	betamx(:,:) = beta_limgz
431	end where
432
433	where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
434	betamy(:,:) = beta_limgz
435	end where
436
437	fleuvemx(:,:)=gzmx(:,:)
438	fleuvemy(:,:)=gzmy(:,:)
439
440	end subroutine gzm_beta_prescr
441
442
443	!----------------------------------------------------------------------------------
444	! Routines qui permettent des manip simples sur beta.
445
446
447
448	!______________________________________________________________________________________
449	!>SUBROUTINE: beta_min_value
450	!! valeur mini que peut avoir beta (en Pa an m-1)
451	!<
452
453	subroutine beta_min_value(nxx,nyy,R_mx,R_my,R_centre,val_betamin)
454	implicit none
455	integer :: nxx,nyy ! dimension des tableaux
456	real, dimension(nxx,nyy), intent(inout) :: R_mx ! valeur du beta sur maille mx
457	real, dimension(nxx,nyy), intent(inout) :: R_my ! valeur du beta sur maille my
458	real, dimension(nxx,nyy), intent(inout) :: R_centre ! valeur du beta sur maille centree
459	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
460
461	R_centre(:,:) = max(R_centre(:,:),val_betamin)
462	R_mx(:,:) = max(R_mx(:,:),val_betamin)
463	R_my(:,:) = max(R_my(:,:),val_betamin)
464
465	where(flot(:,:))
466	drag_centre(:,:) = 0.
467	end where
468
469	end subroutine beta_min_value
470	!______________________________________________________________________________________
471	!>SUBROUTINE: beta_stag2c
472	!! average the staggered values on central ones
473	!<
474
475	subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones
476
477	implicit none
478	integer :: nxx,nyy ! dimension des tableaux
479	real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx
480	real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my
481	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
482
483	R_centre(:,:)=0.
484	do j=2,ny-1
485	do i=2,nx-1
486	R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
487	+ (R_my(i,j)+R_my(i,j+1)))*0.25
488	end do
489	end do
490	end subroutine beta_stag2c
491	!______________________________________________________________________________________
492	!>SUBROUTINE: beta_c2stag_min
493	!! redistribute central values on staggered grid
494	!<
495
496	subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
497
498	implicit none
499	integer :: nxx,nyy ! dimension des tableaux
500	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
501	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
502	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
503
504	do j=2,nyy
505	do i=2,nxx
506	R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
507	R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
508	end do
509	end do
510
511	end subroutine beta_c2stag
512	!______________________________________________________________________________________
513	!>SUBROUTINE: beta_c2stag_min
514	!! redistribute central values on staggered grid
515	!<
516
517	subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid
518
519	implicit none
520	integer :: nxx,nyy ! dimension des tableaux
521	real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx
522	real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my
523	real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree
524
525	do j=2,nyy
526	do i=2,nxx
527	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
528	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
529	end do
530	end do
531
532	end subroutine beta_c2stag_min
533
534	!>SUBROUTINE: update_buoy
535	!! fait changer la valeur centree enfonction de la hauteur au dessus de la flottaison.
536	!< on pressuppose beta = beta0 * Hbuoy
537
538	subroutine update_buoy(nxx,nyy,Ref_centre,Ref_Hbuoy,R_mx,R_my,R_centre,val_betamin)
539	implicit none
540	integer :: nxx,nyy ! dimension des tableaux
541	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
542	real, dimension(nxx,nyy), intent(in) :: Ref_Hbuoy ! valeur de reference hauteur au dessus de la flottaison
543	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
544	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
545	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
546	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
547
548	call calc_buoyency(nxx,nyy,H,B,H_buoy_courant)
549
550
551
552	where (Ref_Hbuoy(:,:).gt.0.)
553	! R_centre(:,:) = Ref_centre(:,:)*max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)
554
555
556	! avec une dépendance non linéaire
557	R_centre(:,:) = Ref_centre(:,:)(max((max(H_buoy_courant(:,:),1.)/max(Ref_Hbuoy(:,:),1.)),0.01))*(1./3.)
558
559
560
561
562	elsewhere (Ref_Hbuoy(:,:).le.0.) ! flottant (0.) ou eventuellement pas de glace
563	R_centre(:,:) = Ref_centre(:,:)
564
565	end where
566
567	! valeur mini
568
569	R_centre(:,:) = max(R_centre(:,:),val_betamin)
570
571	! si flottaison drag = 0.
572
573	where(flot(:,:))
574	drag_centre(:,:) = 0.
575	end where
576
577	! remet sur les noeuds staggered
578
579	do j=2,nyy
580	do i=2,nxx
581	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
582	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
583	end do
584	end do
585
586	! valeur mini
587	R_mx(:,:) = max(R_mx(:,:),val_betamin)
588	R_my(:,:) = max(R_my(:,:),val_betamin)
589
590
591	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
592
593	end subroutine update_buoy
594
595
596	!>SUBROUTINE: calc_buoyency
597	!< calcule la hauteur au dessus de la flottaison.
598
599	subroutine calc_buoyency(nxx,nyy,H1,B1,H_buoy)
600
601	implicit none
602	integer :: nxx,nyy ! dimension des tableaux
603	real, dimension(nxx,nyy), intent(in) :: H1 ! epaisseur
604	real, dimension(nxx,nyy), intent(in) :: B1 ! base de la glace (attention : pas Bsoc)
605	real, dimension(nxx,nyy), intent(out) :: H_buoy
606
607
608	where ( sealevel-B1(:,:).le.0.) ! socle au dessus du niveau des mers
609	H_buoy(:,:) = H1(:,:)
610
611	elsewhere
612	H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
613
614	end where
615
616	end subroutine calc_buoyency
617
618	!>SUBROUTINE: calc_buoyency
619	!< calcule la hauteur au dessus de la flottaison.
620
621	subroutine calc_Uslid_maj(nxx,nyy,ux_stag,uy_stag,u_centre)
622
623	implicit none
624	integer :: nxx,nyy ! dimension des tableaux
625	real, dimension(nxx,nyy), intent(in) :: ux_stag ! vitesse selon ux
626	real, dimension(nxx,nyy), intent(in) :: uy_stag ! vitesse selon uy
627	real, dimension(nxx,nyy), intent(out) :: u_centre ! vitesse sur noeud majeur
628
629	do j=2,nyy-1
630	do i=2,nxx-1
631	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
632	end do
633	end do
634
635	end subroutine calc_Uslid_maj
636
637	subroutine update_nolin(nxx,nyy,Ref_centre,Ref_Uslid,R_mx,R_my,R_centre,val_betamin)
638
639	implicit none
640	integer :: nxx,nyy ! dimension des tableaux
641	real, dimension(nxx,nyy), intent(in) :: Ref_centre ! valeur de reference dragging centre
642	real, dimension(nxx,nyy), intent(in) :: Ref_Uslid ! valeur de reference vitesse de glissement
643	real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx
644	real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my
645	real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree
646	real , intent(in) :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1)
647
648	call calc_Uslid_maj(nxx,nyy,ux(:,:,nz),uy(:,:,nz),Uslid_centre)
649
650	where ((Ref_Uslid(:,:).gt.1.).and.(Uslid_centre(:,:).gt.1.))
651	R_centre(:,:) = Ref_centre(:,:)(Uslid_centre(:,:)/Ref_Uslid(:,:))*exp_alpha
652	! coef_drag(:,:)= Ref_centre(:,:)(1./Ref_Uslid(:,:))*exp_alpha
653	elsewhere
654	R_centre(:,:) = Ref_centre(:,:)
655	! coef_drag(:,:)= Ref_centre(:,:)
656	end where
657
658	! valeur mini
659	R_centre(:,:) = max(R_centre(:,:),val_betamin)
660
661	! si flottaison drag = 0.
662
663	where(flot(:,:))
664	R_centre(:,:) = 0.
665	end where
666
667	! remet sur les noeuds staggered
668
669	do j=2,nyy
670	do i=2,nxx
671	R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
672	R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
673	end do
674	end do
675
676	! valeur mini
677	R_mx(:,:) = max(R_mx(:,:),val_betamin)
678	R_my(:,:) = max(R_my(:,:),val_betamin)
679
680
681	! par contre, je ne reteste pas sur beta_limgz pour eviter d'avoir des points qui sautent
682
683	end subroutine update_nolin
684
685
686	end module dragging_beta_nolin
687

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

Download in other formats: