1 | !> \file tab-litho-0.3.f90 |
---|
2 | !! Repartition des enfoncement en fonction de la rigidite de la lithosphere. |
---|
3 | !< |
---|
4 | |
---|
5 | !> SUBROUTINE: tab_litho |
---|
6 | !! \author ... |
---|
7 | !! \date ... |
---|
8 | !! @note Subroutine qui donne la repartition des enfoncements |
---|
9 | !! en fonction de la rigidite de la lithosphere. |
---|
10 | !! @note variables en entree |
---|
11 | !! - ROM masse volumique du manteau |
---|
12 | !! - RL rayon de rigidite relative |
---|
13 | !! - DL rigidite flexural |
---|
14 | !! |
---|
15 | !! @note variables en sortie |
---|
16 | !! - WE deflection due a une charge unitaire |
---|
17 | !! |
---|
18 | !! @note Used modules: |
---|
19 | !! @note - use module3D_phy |
---|
20 | !! @note - use param_phy_mod |
---|
21 | !! @note - use ISO_DECLAR |
---|
22 | !< |
---|
23 | |
---|
24 | subroutine tab_litho |
---|
25 | |
---|
26 | USE module3D_phy |
---|
27 | USE param_phy_mod |
---|
28 | USE ISO_DECLAR ! module de declaration pour l'isostasie |
---|
29 | |
---|
30 | |
---|
31 | implicit none |
---|
32 | |
---|
33 | |
---|
34 | INTEGER, PARAMETER :: nk=1000 |
---|
35 | REAL, dimension(0:nk) :: kei |
---|
36 | REAL :: stepk,AL,XL,DIST |
---|
37 | REAL :: som |
---|
38 | |
---|
39 | ! pour la lithosphere |
---|
40 | STEPK=100. |
---|
41 | AL=-RL*RL/(2.*PI*DL)*DX*DY |
---|
42 | |
---|
43 | ! fonction de kelvin |
---|
44 | ! lecture de la table kei qui est tous les 0.01 entre 0 et 10 |
---|
45 | ! STEPK=100=1/ecart |
---|
46 | |
---|
47 | open(num_kelvin,file=trim(dirnameinp)//'../kelvin.res') |
---|
48 | read(num_kelvin,*) |
---|
49 | do K=0,NK |
---|
50 | read(num_kelvin,*) XL,kei(K) |
---|
51 | end do |
---|
52 | close(num_kelvin) |
---|
53 | |
---|
54 | |
---|
55 | do I=-LBLOC,LBLOC |
---|
56 | do J=-LBLOC,LBLOC |
---|
57 | DIST=dx*sqrt(1.*(I*I+J*J)) |
---|
58 | XL=DIST/RL*STEPK |
---|
59 | K=int(XL) |
---|
60 | if ((K.gt.834).or.(DIST.gt.dx*LBLOC)) then |
---|
61 | We(I,J)=0. |
---|
62 | else |
---|
63 | We(I,J)=kei(k)+(kei(k+1)-kei(k))*(XL-K) |
---|
64 | if (K.eq.834) We(I,J)=max(We(I,J),0.) |
---|
65 | endif |
---|
66 | We(I,J)=We(I,J)*Al |
---|
67 | end do |
---|
68 | end do |
---|
69 | |
---|
70 | |
---|
71 | |
---|
72 | !open(55,file='iso-cat') |
---|
73 | !write(55,*) 'WE avec LBLOC sur 480 km' |
---|
74 | !do i=-lbloc,lbloc |
---|
75 | ! write(55,'25(i4,1x)') (int(WE(i,j)*1.e08),j=-lbloc,lbloc) |
---|
76 | !end do |
---|
77 | !close(55) |
---|
78 | |
---|
79 | ! sommation des elements de WE dans som pour normaliser |
---|
80 | ! la somme des enfoncements a 1 |
---|
81 | som = SUM(WE(:,:)) |
---|
82 | |
---|
83 | |
---|
84 | ! normalisation |
---|
85 | WE(:,:) = WE(:,:)/(som*ROM*g) |
---|
86 | |
---|
87 | |
---|
88 | |
---|
89 | ! ATTENTION DESSOUS VALABLE SEULEMENT EN 40 km |
---|
90 | |
---|
91 | ! la partie ci-dessous doit etre mise dans une routine specifique Philippe |
---|
92 | ! pour faire exactement la meme isostasie que Philippe |
---|
93 | !if ((icouple.eq.2).or.(icouple.eq.4)) then |
---|
94 | |
---|
95 | ! LBLOC=10 |
---|
96 | |
---|
97 | ! do I=-LBLOC,LBLOC |
---|
98 | ! do J=-LBLOC,LBLOC |
---|
99 | ! We(i,j)=0. ! mise a 0 partout |
---|
100 | ! end do |
---|
101 | ! end do |
---|
102 | ! open(num_phil,file='../INPUT/iso-philippe') |
---|
103 | |
---|
104 | ! read(num_phil,*) |
---|
105 | ! read(num_phil,*) |
---|
106 | ! read(num_phil,*) |
---|
107 | ! read(num_phil,*) |
---|
108 | ! do J=-10,10 |
---|
109 | ! read(num_phil,*)(We(i,j),i=-10,10) |
---|
110 | ! end do |
---|
111 | ! close(num_phil) |
---|
112 | ! do i=-LBLOC,LBLOC |
---|
113 | ! do j=-LBLOC,LBLOC |
---|
114 | ! we(i,j)=we(i,j)/11030/g/rom |
---|
115 | ! end do |
---|
116 | ! end do |
---|
117 | !endif |
---|
118 | |
---|
119 | return |
---|
120 | end subroutine tab_litho |
---|