Context Navigation

source: branches/GRISLIv3/SOURCES/spinup_mod.f90 @ 474

Last change on this file since 474 was 467, checked in by aquiquet, 6 months ago
Cleaning branch: continuing module3D cleaning
File size: 25.3 KB

Line
1	!> \file spinup_mod.f90
2	!! This module allows various spinup tasks mostly based on balance velocities.
3	!<
4
5	!> \namespace spinup_vitbil
6	!! This module allows various spinup tasks mostly based on balance velocities.
7	!! @note Could be transformed to use observed surface velocities
8	!! @note must be used as a dragging module
9	!! @note in module "choix" other dragging module must be switched off
10	!! \author ...
11	!! \date ...
12	!! @note Used modules:
13	!! @note - use module3D_phy
14	!! @note - use param_phy_mod
15	!! @note - use deformation_mod_2lois
16	!<
17
18
19	module spinup_vitbil
20
21	use geography, only: nx,ny,nz
22
23	implicit none
24
25	real, dimension(nx,ny) :: Vcol_x !< vertically averaged velocity x direction (balance)
26	real, dimension(nx,ny) :: Vcol_y !< vertically averaged velocity y direction (balance)
27	real, dimension(nx,ny) :: Ux_deformation !< vertically averaged deformation velocity x direction
28	real, dimension(nx,ny) :: Uy_deformation !< vertically averaged deformation velocity y direction
29	real, dimension(nx,ny) :: coef_defmx !< rescaling coefficient of Sa_mx and s2a_mx
30	real, dimension(nx,ny) :: coef_defmy !< rescaling coefficient of Sa_my and s2a_my
31	real, dimension(nx,ny) :: init_coefmx !< rescalling coefficient before limitation
32	real, dimension(nx,ny) :: init_coefmy !< rescalling coefficient before limitation
33
34	! compatibilites avec remplimat
35	logical :: corr_def = .False. !< for deformation correction (compatibility)
36
37	real, dimension(nx,ny) :: Vsl_x !<
38	real, dimension(nx,ny) :: Vsl_y !<
39
40	integer :: type_vitbil ! defini le type de fichier a lire : lect_vitbil ou lect_vitbil_Lebrocq
41
42	contains
43	!--------------------------------------------------------------------------------
44
45	!> SUBROUTINE: init_spinup
46	!! Initialize spinup
47	!<
48	subroutine init_spinup
49
50	use module3D_phy, only: ispinup,num_param,num_rep_42
51	use runparam, only : itracebug
52
53	namelist/spinup/ispinup,type_vitbil
54
55	! put here which type of spinup
56	! ispinup = 0 ! run standard
57	! ispinup = 1 ! ispinup = 1 -> saute icethick3, diagnoshelf,
58	! diffusiv, isostasie pour equilibre temperature
59	! en prenant les vitesses calculees aux premiers
60	! pas de temps
61	! ispinup = 2 ! ispinup = 2 -> fait la conservation
62	! ! de la masse avec les vitesses de bilan
63	! ispinup = 3 ! fait le calcul des temperatures
64	! ! avec les vitesses de bilan
65
66
67
68
69	! lecture des parametres
70
71	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
72	read(num_param,spinup)
73
74	write(num_rep_42,*)'!___________________________________________________________'
75	write(num_rep_42,*)'! spinup module spinup_vitbil '
76	write(num_rep_42,spinup) ! pour ecrire les parametres lus
77	write(num_rep_42,*)
78	write(num_rep_42,*)
79	write(num_rep_42,*)'! ispinup = 0 run standard voir no_spinup'
80	write(num_rep_42,*)'! ispinup = 1 equilibre temperature avec vitesses grisli voir no_spinup'
81	write(num_rep_42,*)
82	write(num_rep_42,*)'! ispinup >1 use spinup_vitbil'
83	write(num_rep_42,*)'! ispinup = 2 conservation de la masse avec vitesses bilan '
84	write(num_rep_42,*)'! ispinup = 3 equilibre temperature avec vitesses bilan'
85	write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
86	write(num_rep_42,*)
87
88
89	if (ispinup.eq.0) then
90	write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
91	write(6,*) 'et il faut rajouter un dragging'
92	stop
93	endif
94
95	if (ispinup.eq.3) then
96	select case (type_vitbil)
97	case (1)
98	!cdc Methode Cat lecture fichier selon x et y sur grille stag
99	call lect_vitbil
100	case(2)
101	!cdc Methode Tof lecture fichier Lebrocq vitesse et direction
102	call lect_vitbil_Lebrocq
103	case default
104	print*,'type_vitbil valeur invalide dans spinup_vitbil'
105	stop
106	end select
107	endif
108
109	if (itracebug.eq.1) call tracebug(' fin routine init_spinup de spinup_vitbil')
110	end subroutine init_spinup
111
112	subroutine lect_vitbil
113
114	use module3D_phy, only: debug_3D,num_rep_42,num_param
115	use io_netcdf_grisli, only: Read_Ncdf_var
116	use geography, only: dirnameinp
117	use runparam, only : itracebug
118
119	character(len=100) :: balance_Ux_file ! balance velocity file Ux
120	character(len=100) :: balance_Uy_file ! balance velocity file Uy
121	real*8, dimension(:,:), pointer :: tab !< tableau 2d real pointer
122
123	namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file
124
125	if (itracebug.eq.1) call tracebug(' Subroutine lect_vitbil')
126
127	! lecture des parametres du run
128	!--------------------------------------------------------------------
129
130	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
131	428 format(A)
132	read(num_param,vitbil_upwind)
133
134	write(num_rep_42,428) '!___________________________________________________________'
135	write(num_rep_42,428) '!read balance velocities on staggered grid'
136	write(num_rep_42,vitbil_upwind)
137	write(num_rep_42,428) '! balance_Ux_file : nom du fichier qui contient les vit. bilan Ux'
138	write(num_rep_42,428) '! balance_Uy_file : nom du fichier qui contient les vit. bilan Uy'
139	write(num_rep_42,*)
140
141	! read balance velocities on staggered nodes
142
143	balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file)
144	balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file)
145
146	call Read_Ncdf_var('Vcol_x',balance_Ux_file,tab)
147	Vcol_x(:,:)=tab(:,:)
148	call Read_Ncdf_var('Vcol_y',balance_Uy_file,tab)
149	Vcol_y(:,:)=tab(:,:)
150
151	! call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc') !read Vcol_x
152	! call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc') !read Vcol_y
153	!call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file) ! read Vcol_x
154	!call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file) ! read Vcol_y
155
156	debug_3D(:,:,59)=Vcol_x(:,:)
157	debug_3D(:,:,60)=Vcol_y(:,:)
158
159	end subroutine lect_vitbil
160
161	subroutine lect_vitbil_Lebrocq
162
163	use module3D_phy, only: num_rep_42,mstream_mx,num_param,pi
164	use io_netcdf_grisli, only: Read_Ncdf_var
165	use geography, only: dirnameinp
166	use runparam, only : itracebug
167
168	character(len=100) :: vit_balance_file ! balance velocity file
169	character(len=100) :: flowdir_balance_file ! balance flow direction file
170	real*8, dimension(:,:), pointer :: tab !< tableau 2d real pointer
171	real, dimension(nx,ny) :: V_Lebrocq !< vertically averaged velocity (balance)
172	real, dimension(nx,ny) :: flowdir_Lebrocq !< flow direction (degree)
173	real, dimension(nx,ny) :: Vx_Lebrocq, Vy_Lebrocq ! vitesses selon x et y
174	integer, dimension(nx,ny) :: v_good !points avec vitesse Lebrocq OK
175	integer :: i,j
176
177	namelist/vitbil_Lebrocq/vit_balance_file, flowdir_balance_file
178
179	if (itracebug.eq.1) call tracebug(' Subroutine lect_vitbil_Lebrocq')
180
181	! lecture des parametres du run
182	!--------------------------------------------------------------------
183
184	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
185	428 format(A)
186	read(num_param,vitbil_Lebrocq)
187
188	write(num_rep_42,428) '!___________________________________________________________'
189	write(num_rep_42,428) '!read balance velocities from Lebrocq code'
190	write(num_rep_42,vitbil_Lebrocq)
191	write(num_rep_42,428) '! vit_balance_file : balance velocity file'
192	write(num_rep_42,428) '! flowdir_balance_file : balance flow direction file'
193	write(num_rep_42,*)
194
195	! read balance velocities on staggered nodes
196
197	vit_balance_file = trim(dirnameinp)//trim(vit_balance_file)
198	flowdir_balance_file = trim(dirnameinp)//trim(flowdir_balance_file)
199
200	call Read_Ncdf_var('z',vit_balance_file,tab)
201	V_Lebrocq(:,:)=tab(:,:)
202	call Read_Ncdf_var('z',flowdir_balance_file,tab)
203	flowdir_Lebrocq(:,:)=tab(:,:)*PI/180.
204
205	where (V_Lebrocq(:,:).ge.0.)
206	v_good(:,:)=1
207	elsewhere
208	v_good(:,:)=0
209	endwhere
210
211	do j=2,ny-1
212	do i=2,nx-1
213	if (V_Lebrocq(i,j).lt.0.) then
214
215	V_Lebrocq(i,j)= (V_Lebrocq(i-1,j)v_good(i-1,j) + V_Lebrocq(i+1,j)v_good(i+1,j) + &
216	V_Lebrocq(i,j-1)v_good(i,j-1) + V_Lebrocq(i,j+1)v_good(i,j+1)) / &
217	(v_good(i-1,j)+v_good(i+1,j)+v_good(i,j-1)+v_good(i,j-1))
218	endif
219	enddo
220	enddo
221
222
223	where (V_Lebrocq(:,:).ge.100.)
224	V_Lebrocq(:,:)=100.
225	endwhere
226
227	! calcul des vitesses selon x et selon y :
228	where (V_Lebrocq(:,:).ge.0.)
229	Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
230	Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
231	elsewhere
232	Vx_Lebrocq(:,:) = 0.
233	Vy_Lebrocq(:,:) = 0.
234	endwhere
235
236
237	! calcul des vitesses sur les points staggered
238
239	do j=2,ny
240	do i=2,nx
241	Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
242	Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
243	enddo
244	enddo
245
246	end subroutine lect_vitbil_Lebrocq
247	!____________________________________________________________________________________
248
249	!> SUBROUTINE: init_dragging
250	!! @ Cette routine fait l'initialisation du dragging.
251	!<
252	subroutine init_dragging ! Cette routine fait l'initialisation du dragging.
253
254	use module3D_phy, only: mstream_mx,mstream_my,drag_mx,drag_my,gzmx,gzmy,flgzmx,flgzmy,flotmx,flotmy
255
256	implicit none
257
258	mstream_mx(:,:)=0 ! pas de dragging a l'initialisation
259	mstream_my(:,:)=0
260
261	! coefficient permettant de modifier le basal drag.
262	drag_mx(:,:)=1.
263	drag_my(:,:)=1.
264	gzmx(:,:)=.false.
265	gzmy(:,:)=.false.
266	flgzmx(:,:)=flotmx(:,:)
267	flgzmy(:,:)=flotmy(:,:)
268
269	end subroutine init_dragging
270
271
272	subroutine dragging
273
274	end subroutine dragging
275	!----------------------------------------------------------------------
276
277
278	subroutine force_balance_vel
279
280	use module3D_phy, only: uxbar, uybar,debug_3D
281	use runparam, only : itracebug
282	if (itracebug.eq.1) call tracebug(' Subroutine force_balance_vel')
283
284	Uxbar(:,:)=Vcol_x(:,:)
285	Uybar(:,:)=Vcol_y(:,:)
286	debug_3D(:,:,59)=Uxbar(:,:)
287	debug_3D(:,:,60)=Uybar(:,:)
288	end subroutine force_balance_vel
289
290
291
292	!> SUBROUTINE: calc_coef_vitbil
293	!! @ Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
294	!! @note Must be called after flowlaw and before velocities
295	!<
296	subroutine calc_coef_vitbil
297
298	use module3D_phy, only: ddbx,ddby,gzmx,gzmy,flgzmx,flgzmy,fleuvemx,fleuvemy,flot,flotmx,flotmy,&
299	coefmxbmelt,coefmybmelt,sdx,sdy,uxbar,uybar,uxdef,uydef,ubx,uby,debug_3D,iglen
300	use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
301	use runparam, only : itracebug
302
303	integer :: i,j,k
304
305	ddbx(:,:)=0.
306	ddby(:,:)=0.
307	gzmx(:,:)=.true.
308	gzmy(:,:)=.true.
309	flgzmx(:,:)=.false.
310	flgzmy(:,:)=.false.
311	fleuvemx(:,:)=.false.
312	fleuvemy(:,:)=.false.
313
314	do j=2,ny
315	do i=2,nx
316	if (flot(i,j).and.(flot(i-1,j))) then
317	flotmx(i,j)=.true.
318	flgzmx(i,j)=.true.
319	gzmx(i,j)=.false.
320	end if
321	if (flot(i,j).and.(flot(i,j-1))) then
322	flotmy(i,j)=.true.
323	flgzmy(i,j)=.true.
324	gzmy(i,j)=.false.
325	end if
326	end do
327	end do
328
329	where (coefmxbmelt(:,:).le.0.)
330	gzmx(:,:)=.false.
331	endwhere
332	where (coefmybmelt(:,:).le.0.)
333	gzmy(:,:)=.false.
334	endwhere
335	flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
336	flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
337	fleuvemx(:,:)=gzmx(:,:)
338	fleuvemy(:,:)=gzmy(:,:)
339
340	if (itracebug.eq.1) call tracebug(' Subroutine calc_coef_vitbil')
341	call slope_surf
342
343	where (abs(sdx(:,:)).lt.1.e-6)
344	sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
345	end where
346	where (abs(sdy(:,:)).lt.1.e-6)
347	sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
348	end where
349
350
351	call calc_ubar_def
352
353	if (itracebug.eq.1) call tracebug(' apres calc_ubar_def')
354
355	!---------------compute coef_defmx -----------------------------------
356
357	do j=1,ny
358	do i=1,nx
359
360	flotx: if (flotmx(i,j)) then
361
362	coef_defmx(i,j) = 1.
363	Uxbar(i,j) = Vcol_x(i,j)
364	! dmr below is a cast from a real4 to logical4
365	! dmr cannot be implicit in gfortran
366	! dmr flgzmx(i,j) = Vcol_x(i,j)
367	flgzmx(i,j) = (nint(Vcol_x(i,j)).ne.0)
368	uxdef(i,j) = 0.
369	Ubx(i,j) = Vcol_x(i,j)
370
371	else
372	coldx: if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then !base froide
373	!meme test que dans sliding Bindshadler
374	! sans doute trop fort
375	ddbx(i,j) = 0.
376	Ubx(i,j) = 0.
377
378	if (abs(Ux_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def
379	coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
380
381	else
382	coef_defmx(i,j)=1.
383	end if
384	else ! base temperee
385	if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
386	coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
387	ddbx(i,j) = 0
388	Ubx(i,j) = 0.
389	else
390	coef_defmx(i,j)=1. ! no rescaling of deformation
391	ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j) ! pb si directions opposees
392	Ubx(i,j) = Vcol_x(i,j)-Ux_deformation(i,j)
393	end if
394	end if coldx
395
396	end if flotx
397
398	end do
399	end do
400
401
402	do j=1,ny
403	do i=1,nx
404
405	floty: if (flotmy(i,j)) then
406	coef_defmy(i,j) = 1.
407	Uybar(i,j) = Vcol_y(i,j)
408	! dmr below is a cast from a real4 to logical4
409	! dmr cannot be implicit in gfortran
410	! dmr flgzmy(i,j) = Vcol_y(i,j)
411	flgzmy(i,j) = (nint(Vcol_y(i,j)).ne.0)
412
413	uydef(i,j) = 0.
414	Uby(i,j) = Vcol_y(i,j)
415
416	if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test0')
417	else
418	coldy: if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then !base froide
419	!meme test que dans sliding Bindshadler
420	ddby(i,j) = 0.
421	Uby(i,j) = 0.
422
423	if (abs(Uy_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def
424	coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
425	if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test1')
426	else
427	if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test2')
428	coef_defmy(i,j)=1.
429	end if
430	else ! base temperee
431	if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! que deformation
432	coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
433	ddby(i,j) = 0.
434	Uby(i,j) = 0.
435	if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test3')
436	else
437	coef_defmy(i,j)=1. ! pas de rescaling de la deformation
438	ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j) ! pb si directions opposees
439	Uby(i,j) = Vcol_y(i,j)-Uy_deformation(i,j)
440	if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144)) call tracebug(' test4')
441
442	end if
443	end if coldy
444
445	end if floty
446
447	end do
448	end do
449	if (itracebug.eq.1) call tracebug(' apres test floty')
450
451	! remarque coef peut être <1 a cause de la pente locale, mais il joue toujours son role
452	debug_3D(:,:,61)=coef_defmx(:,:)
453	debug_3D(:,:,62)=coef_defmy(:,:)
454
455
456	calib: do iglen=n1poly,n2poly
457	do k=1,nz
458	SA_mx(:,:,k,iglen) = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
459	S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
460	SA_my(:,:,k,iglen) = SA_my(:,:,k,iglen)*coef_defmy(:,:)
461	S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
462	end do
463	end do calib
464
465	end subroutine calc_coef_vitbil
466
467	!> SUBROUTINE: limit_coef_vitbil
468	!! Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
469	!! @note The difference with calc_coef_vitbil is that the coefficient is limited to the range 0.5-2
470	!! @note - if > 2, sliding is assumed.
471	!! @note - if < 0.5, it means that ice is actually colder than computed (to be implemented)
472	!! @note Must be called after flowlaw and before velocities
473	!<
474
475	subroutine limit_coef_vitbil
476
477	use module3D_phy, only: sdx,sdy,flotmx,flotmy,coefmxbmelt,coefmybmelt,gzmx,gzmy,ddbx,ddby,debug_3D,&
478	iglen
479	use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
480	use runparam, only : itracebug
481
482	integer :: i,j,k
483
484	call slope_surf
485
486	where (abs(sdx(:,:)).lt.1.e-6)
487	sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
488	end where
489	where (abs(sdy(:,:)).lt.1.e-6)
490	sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
491	end where
492
493
494	call calc_ubar_def
495
496	if (itracebug.eq.1) call tracebug(' Entree dans routine calc_coef_vitbil')
497	!---------------compute coef_defmx -----------------------------------
498
499	do j=1,ny
500	do i=1,nx
501
502	flotx: if (flotmx(i,j)) then
503	coef_defmx(i,j)=1.
504	init_coefmx(i,j)= coef_defmx(i,j)
505	else
506	coldx: if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then !base froide
507	!meme test que dans sliding Bindshadler
508	! sans doute trop fort
509	ddbx(i,j)=0.
510	if (abs(Ux_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def
511	coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
512	init_coefmx(i,j)= coef_defmx(i,j)
513
514	if (coef_defmx(i,j).gt.1.) then
515	coef_defmx(i,j)=1.
516	ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)
517	end if
518	else
519	coef_defmx(i,j)=1.
520	init_coefmx(i,j)= coef_defmx(i,j)
521	end if
522	else ! base temperee
523	if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
524	coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
525	init_coefmx(i,j)= coef_defmx(i,j)
526	ddbx(i,j) = 0
527	else
528	coef_defmx(i,j)=1. ! no rescaling of deformation
529	init_coefmx(i,j)= coef_defmx(i,j)
530	ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j) ! pb si directions opposees
531	end if
532	end if coldx
533
534	end if flotx
535
536	end do
537	end do
538
539
540	do j=1,ny
541	do i=1,nx
542
543	floty: if (flotmy(i,j)) then
544	coef_defmy(i,j)=1.
545	init_coefmy(i,j)= coef_defmy(i,j)
546	else
547	coldy: if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then !base froide
548	!meme test que dans sliding Bindshadler
549	! sans doute trop fort
550	ddby(i,j)=0.
551	if (abs(Uy_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def
552	coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
553	init_coefmy(i,j)= coef_defmy(i,j)
554
555	if (coef_defmy(i,j).gt.1.) then
556	coef_defmy(i,j)=1.
557	ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)
558	end if
559	else
560	coef_defmy(i,j)=1.
561	init_coefmy(i,j)= coef_defmy(i,j)
562	end if
563	else ! base temperee
564	if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! only deformation
565	coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
566	init_coefmy(i,j)= coef_defmy(i,j)
567	ddby(i,j) = 0
568	else
569	coef_defmy(i,j)=1. ! no rescaling of deformation
570	init_coefmy(i,j)= coef_defmy(i,j)
571	ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j) ! pb si directions opposees
572	end if
573	end if coldy
574
575	end if floty
576
577	end do
578	end do
579
580
581	! remarque coef peut être <1 a cause de la pente locale, mais il joue toujours son role
582	debug_3D(:,:,73)=init_coefmx(:,:)
583	debug_3D(:,:,74)=init_coefmy(:,:)
584
585	debug_3D(:,:,61)=coef_defmx(:,:)
586	debug_3D(:,:,62)=coef_defmy(:,:)
587
588
589	calib: do iglen=n1poly,n2poly
590	do k=1,nz
591	SA_mx(:,:,k,iglen) = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
592	S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
593	SA_my(:,:,k,iglen) = SA_my(:,:,k,iglen)*coef_defmy(:,:)
594	S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
595	end do
596	end do calib
597
598	! call ajust_ghf ATTENTION NE PAS ACTIVER SAUF TESTS
599
600	! debug_3D(:,:,73)= init_coefmx(:,:)
601	! debug_3D(:,:,74)= init_coefmy(:,:)
602
603
604	end subroutine limit_coef_vitbil
605
606
607	!> SUBROUTINE: calc_ubar_def
608	!! Calculate velocity due to deformation before calibration
609	!<
610
611	subroutine calc_ubar_def ! calculate velocity due to deformation before calibration
612
613	use module3D_phy, only: iglen,flotmx,flotmy,hmx,hmy,sdx,sdy,slope2mx,slope2my
614	use deformation_mod_2lois, only: n1poly,n2poly,glen,s2a_mx,s2a_my,ddx,ddy
615	use geography, only: dx,dy
616	use param_phy_mod, only: rog
617	use runparam, only: itracebug
618
619	implicit none ! extrait de diffusiv pour calculer la partie deformation
620
621	real :: glenexp
622	real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
623	real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
624	integer :: i,j
625
626	inv_4dx=1./(4*dx)
627	inv_4dy=1./(4*dy)
628
629
630	do iglen=n1poly,n2poly
631	glenexp=max(0.,(glen(iglen)-1.)/2.)
632
633	do j=1,ny
634	do i=1,nx
635
636	if (.not.flotmy(i,j)) then
637	ddy(i,j,iglen)=((slope2my(i,j)** & ! SLOPE2my calcule dans slope_surf
638	glenexp)(rog)*glen(iglen)) &
639	Hmy(i,j)*(glen(iglen)+1)
640
641	endif
642	if (.not.flotmx(i,j)) then
643	ddx(i,j,iglen)=((slope2mx(i,j)** & ! slope2mx calcule en debut de routine
644	glenexp)(rog)*glen(iglen)) &
645	Hmx(i,j)*(glen(iglen)+1)
646	endif
647
648	end do
649	end do
650	end do
651
652	do j=2,ny-1
653	do i=2,nx-1
654	ux_deformation(i,j)=0.
655	Uy_deformation(i,j)=0.
656	do iglen=n1poly,n2poly
657	ux_deformation(i,j)=ux_deformation(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen)
658	Uy_deformation(i,j)=Uy_deformation(i,j)+ddy(i,j,iglen)*s2a_my(i,j,1,iglen)
659	end do
660	end do
661	end do
662
663	do j=2,ny-1
664	do i=2,nx-1
665	ux_deformation(i,j)=ux_deformation(i,j)*(-sdx(i,j))
666	Uy_deformation(i,j)=Uy_deformation(i,j)*(-sdy(i,j))
667	end do
668	end do
669	!debug_3D(:,:,63)=ux_deformation(:,:)
670	!debug_3D(:,:,64)=uy_deformation(:,:)
671
672	if (itracebug.eq.1) call tracebug(' fin de calc_ubar_def ')
673
674
675	end subroutine calc_ubar_def
676
677	!> SUBROUTINE: ajust_ghf
678	!! Ajuste le flux geothermique pour avoir une temperature coherente
679	!! avec les vitesses de bilan
680	!<
681
682	subroutine ajust_ghf
683
684	use module3D_phy, only: debug_3D,flot,ibase,ghf,secyear
685
686	implicit none
687
688	real,dimension(nx,ny) :: ghf0 !< geothermal heat flux J/m2/a 'o'
689	real,dimension(nx,ny) :: coefdef_maj !< coefficient deformation
690	real :: increment_ghf
691	real :: ghf_min
692	integer :: i,j
693
694	increment_ghf=0.00001 !< exprime en W/m2
695	ghf_min=0.025
696
697	do j=2,ny-1
698	do i=2,nx-1
699	if (.not. flot(i,j)) then
700	coefdef_maj(i,j)= (init_coefmx(i,j)+init_coefmx(i+1,j))+ &
701	(init_coefmy(i,j)+init_coefmy(i,j+1))
702	coefdef_maj(i,j)=0.25*coefdef_maj(i,j)
703
704	! un coef trop grand indique eventuellement une base trop froide
705	if ((coefdef_maj(i,j).gt.1.5).and.(ibase(i,j).eq.1)) then ! en base froide
706	ghf(i,j)=ghf(i,j)-SECYEAR*increment_ghf ! attention ghf est negatif
707
708	! un coef trop petit indique une base trop chaude
709	else if ((coefdef_maj(i,j).lt.0.7).and.(ghf(i,j).lt.-secyear*ghf_min)) then
710	ghf(i,j)=ghf(i,j)+SECYEAR*increment_ghf ! attention ghf est negatif
711	end if
712	end if
713	end do
714	end do
715
716	debug_3D(:,:,75)=(ghf0(:,:)-ghf(:,:))*1000./secyear
717
718	end subroutine ajust_ghf
719
720
721
722	end module spinup_vitbil
723
724
725
726

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

Download in other formats: