1 | !> \file litho-0.3.f90 |
---|
2 | !! Donne la repartition des enfoncements |
---|
3 | !! en fonction de la rigidite de la lithosphere. |
---|
4 | !< |
---|
5 | |
---|
6 | !> SUBROUTINE: litho |
---|
7 | !!Petit routine qui donne la repartition des enfoncements |
---|
8 | !!en fonction de la rigidite de la lithosphere. |
---|
9 | !! \author .... |
---|
10 | !! \date 10 Novembre 1999 |
---|
11 | !! @note - En entree |
---|
12 | !! @note WE(-LBLOC:LBLOC,-LBLOC:LBLOC) : deflection due a une charge unitaire |
---|
13 | !! defini dans tab-litho |
---|
14 | !! @note LBLOC : relie a la distance : distance en noeud autour de laquelle |
---|
15 | !! la flexure de la lithosphere est calculee |
---|
16 | !! @note CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) : poids par unite de surface |
---|
17 | !! (unite ?) au temps time, calcule avant 'appel a litho |
---|
18 | !! dans taubed ou initial2 |
---|
19 | !! @note - En sortie |
---|
20 | !! @note W1(nx,ny) : enfoncement courant |
---|
21 | !! @note - tableaux de travail. |
---|
22 | !! @note WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC) : enfoncement au point i,j du a la charge |
---|
23 | !! au point IP,JP |
---|
24 | !! @note LPX et LPY sont les coordonnees entre (-LBLOC:LBLOC,-LBLOC:LBLOC) autour |
---|
25 | !! du point I,J pour calculer tous les WLOC |
---|
26 | !! @note IP et JP sont les coordonnees de CHARGE dans un repere : |
---|
27 | !! (1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) |
---|
28 | !! @note Used modules: |
---|
29 | !! @note - use module3D_phy |
---|
30 | !! @note - use ISO_DECLAR |
---|
31 | !< |
---|
32 | subroutine litho |
---|
33 | |
---|
34 | USE module3D_phy |
---|
35 | USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie |
---|
36 | |
---|
37 | implicit none |
---|
38 | |
---|
39 | INTEGER :: IP,JP,LPX,LPY,II,SOM1,SOM2 |
---|
40 | REAL, dimension(:,:), allocatable :: WLOC |
---|
41 | REAL, dimension(:), allocatable :: croix |
---|
42 | |
---|
43 | !----- allocation de WLOC et de croix ----------- |
---|
44 | |
---|
45 | if (.not.allocated(WLOC)) THEN |
---|
46 | allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err) |
---|
47 | if (err/=0) then |
---|
48 | print *,"Erreur a l'allocation du tableau WLOC",err |
---|
49 | stop 4 |
---|
50 | end if |
---|
51 | end if |
---|
52 | |
---|
53 | if (.not.allocated(croix)) THEN |
---|
54 | allocate(croix(0:LBLOC),stat=err) |
---|
55 | if (err/=0) then |
---|
56 | print *,"Erreur a l'allocation du tableau croix",err |
---|
57 | stop 4 |
---|
58 | end if |
---|
59 | end if |
---|
60 | !----- fin de l'allocation -------------- |
---|
61 | |
---|
62 | |
---|
63 | ! calcul de la deflexion par sommation des voisins appartenant |
---|
64 | ! au bloc de taille 2*LBLOC |
---|
65 | som1=0. |
---|
66 | som2=0. |
---|
67 | |
---|
68 | ! On somme aussi les contributions des points exterieurs au domaine |
---|
69 | ! lorsque la charge est due a l'ocean. On suppose alors que |
---|
70 | ! ces points ont la meme charge que les limites |
---|
71 | |
---|
72 | boucleij : do J=1,NY |
---|
73 | do I=1,NX |
---|
74 | |
---|
75 | |
---|
76 | W1(I,J)=0. |
---|
77 | ii=0 |
---|
78 | |
---|
79 | ! IP,JP compteurs de la sommation dans le repere |
---|
80 | ! (1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) |
---|
81 | |
---|
82 | ! LPX,LPY compteurs dans le repere -LBLOC:LBLOC |
---|
83 | |
---|
84 | ! do IP=I-LBLOC,I+LBLOC |
---|
85 | ! do JP=J-LBLOC,J+LBLOC |
---|
86 | ! LPX = IP-I |
---|
87 | ! LPY = JP-J |
---|
88 | ! WLOC(LPX,LPY)=We(LPX,LPY)*CHARGE(IP,JP) |
---|
89 | ! ii=ii+1 |
---|
90 | ! end do |
---|
91 | ! end do |
---|
92 | ! ecriture tableau f90 de la boucle ci-dessus : |
---|
93 | |
---|
94 | WLOC(:,:) = WE(:,:) * CHARGE(I-LBLOC:I+LBLOC,J-LBLOC:J+LBLOC) |
---|
95 | |
---|
96 | |
---|
97 | |
---|
98 | ! fin de l'attribution des wloc |
---|
99 | |
---|
100 | |
---|
101 | ! sommation de tous les effets (tentative pour |
---|
102 | ! eviter les erreurs d'arrondi) |
---|
103 | W1(I,J)=WLOC(0,0) |
---|
104 | |
---|
105 | ! dans croix on calcul la somme des effets WLOC puis on met le resultat dans W1 |
---|
106 | do LPX=1,LBLOC |
---|
107 | LPY=0 |
---|
108 | ! print*,i,j,LPX,WLOC(LPX,LPY),WLOC(-LPX,LPY),WLOC(LPY,LPX),WLOC(LPY,-LPX) |
---|
109 | ! print*,(WLOC(LPX,LPY)+WLOC(-LPX,LPY)) |
---|
110 | ! print*,(WLOC(LPY,LPX)+WLOC(LPY,-LPX)) |
---|
111 | ! CROIX(LPY)=(WLOC(LPX,LPY)+WLOC(-LPX,LPY)) |
---|
112 | |
---|
113 | ! CROIX(LPY)=(WLOC(LPY,LPX)+WLOC(LPY,-LPX)) + CROIX(LPY) |
---|
114 | |
---|
115 | |
---|
116 | ! write(6,*) i,j,lpy,lpx,lbloc,WLOC(LPY,LPX) |
---|
117 | ! write(6,*) WLOC(LPY,-LPX),WLOC(LPX,LPY),WLOC(-LPX,LPY) |
---|
118 | ! write(6,*) (WLOC(LPY,LPX)+WLOC(LPY,-LPX))+(WLOC(LPX,LPY)+WLOC(-LPX,LPY)) |
---|
119 | CROIX(LPY)=(WLOC(LPY,LPX)+WLOC(LPY,-LPX))+(WLOC(LPX,LPY)+WLOC(-LPX,LPY)) |
---|
120 | |
---|
121 | |
---|
122 | ! print*,'croix ',croix(lpy) |
---|
123 | LPY=LPX |
---|
124 | CROIX(LPY)= ((WLOC(LPX,LPY)+WLOC(LPX,-LPY)) & |
---|
125 | +(WLOC(-LPX,LPY)+WLOC(-LPX,-LPY))) |
---|
126 | do LPY=1,LPX-1 |
---|
127 | CROIX(LPY)= (((WLOC(LPX,LPY)+WLOC(LPX,-LPY)) & |
---|
128 | +(WLOC(-LPX,LPY)+WLOC(-LPX,-LPY))) & |
---|
129 | +((WLOC(LPY,LPX)+WLOC(LPY,-LPX)) & |
---|
130 | +(WLOC(-LPY,LPX)+WLOC(-LPY,-LPX)))) |
---|
131 | end do |
---|
132 | do LPY=0,LPX |
---|
133 | W1(I,J)=W1(I,J)+CROIX(LPY) ! sommation dans W1 |
---|
134 | end do |
---|
135 | |
---|
136 | end do |
---|
137 | |
---|
138 | |
---|
139 | ! if (ii.ne.(2*LBLOC+1)*(2*LBLOC+1)) then |
---|
140 | ! write(6,*)'attention mauvais compteur pour i,j',I,J,II |
---|
141 | ! endif |
---|
142 | ! ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC |
---|
143 | som1=som1+W1(I,J) |
---|
144 | som2=som2-charge(I,J)/ROMG |
---|
145 | |
---|
146 | end do |
---|
147 | end do boucleij |
---|
148 | |
---|
149 | ! write(6,*)'enfoncement total avec rigidite de la lithosphere:', |
---|
150 | ! & som1 |
---|
151 | ! write(6,*)'enfoncement total avec isostasie locale :', |
---|
152 | ! & som2 |
---|
153 | |
---|
154 | deallocate(WLOC,stat=err) |
---|
155 | if (err/=0) then |
---|
156 | print *,"Erreur a l'allocation du tableau WLOC",err |
---|
157 | stop 4 |
---|
158 | end if |
---|
159 | |
---|
160 | deallocate(croix,stat=err) |
---|
161 | if (err/=0) then |
---|
162 | print *,"Erreur a l'allocation du tableau croix",err |
---|
163 | stop 4 |
---|
164 | end if |
---|
165 | return |
---|
166 | end subroutine litho |
---|
167 | |
---|