1 | !> \file flowlaw-0.3.f90 |
---|
2 | !!Subroutine qui calcul les parametres pour chaque loi de deformation |
---|
3 | !< |
---|
4 | |
---|
5 | !> SUBROUTINE: FLOWLAW |
---|
6 | !! \author ... |
---|
7 | !! \date 25 Novembre 1999 |
---|
8 | !! |
---|
9 | !!@note subroutine qui calcul SA et S2A et BTT pour chaque loi de deformation (exposant variable) |
---|
10 | !! |
---|
11 | !! variables en entree : - iglen : precise sur quel exposant de la loi de deformation |
---|
12 | !! flowlaw doit travailler. |
---|
13 | !! - glen : tableau contenant les exposants des lois |
---|
14 | !! - BAT1 et BAT2 |
---|
15 | !! - TTRANS |
---|
16 | !! - Q1 et Q2 |
---|
17 | !! - T |
---|
18 | !! - RGAS |
---|
19 | !! - TI2, TETAR et TETA : tableaux de travail calcule dans |
---|
20 | !! la routine flow_general |
---|
21 | !! - E |
---|
22 | !! |
---|
23 | !! variables en sortie : - SA et S2A |
---|
24 | !! - BTT |
---|
25 | !! |
---|
26 | !! @note Used modules: |
---|
27 | !! @note - use module3D_phy |
---|
28 | !! @note - use module_choix |
---|
29 | !> |
---|
30 | |
---|
31 | SUBROUTINE FLOWLAW (iiglen) |
---|
32 | |
---|
33 | USE module3D_phy |
---|
34 | ! USE DEFORM_DECLAR |
---|
35 | USE module_choix |
---|
36 | |
---|
37 | implicit none |
---|
38 | INTEGER :: iiglen !< compteur pour la boucle flowlaw |
---|
39 | REAL :: a5 !< exposant de la loi de def +2 |
---|
40 | REAL :: a4 !< exposant de la loi de def +1 |
---|
41 | REAL :: TI1 !< utilise pour le calcul de btt a k=1 |
---|
42 | REAL :: BAT !< prend la valeur de BAT1 ou BAT2 selon les cas |
---|
43 | REAL :: Q !< prend la valeur de Q1 ou Q2 selon les cas |
---|
44 | REAL :: AAT !< variable de travail =Q*TETAR(I,J,K) |
---|
45 | REAL :: SSEC !< variable de travail |
---|
46 | REAL :: ZETAT !< variable de travail : profondeur ou T = TRANS(iiglen) |
---|
47 | REAL,dimension(NX,NY,NZ) :: SI1 !< tableau de calcul |
---|
48 | REAL,dimension(NX,NY,NZ) :: SI2 !< tableau de calcul |
---|
49 | real,dimension(nx,ny) :: tab_mx |
---|
50 | real,dimension(nx,ny) :: tab_my |
---|
51 | real,dimension(nx,ny) :: tab2D |
---|
52 | |
---|
53 | REAL,dimension(NZ) :: E ! vertical coordinate in ice, scaled to H zeta |
---|
54 | REAL :: De= 1./(NZ-1) |
---|
55 | |
---|
56 | E(1)=0. |
---|
57 | E(NZ)=1. |
---|
58 | do K=1,NZ |
---|
59 | if ((K.ne.1).and.(K.ne.NZ)) E(K)=(K-1.)/(NZ-1.) |
---|
60 | end do |
---|
61 | |
---|
62 | ! exposant de la loi de deformation utilisee : glen(iiglen) |
---|
63 | ! l'exposant correspondant a iiglen est defini dans DEFORMATION_MOD |
---|
64 | a5 = glen(iiglen) + 2 |
---|
65 | a4 = glen(iiglen) + 1 |
---|
66 | |
---|
67 | ! boucle pour btt a k=1 |
---|
68 | TI1=273.15*273.15 |
---|
69 | do J=2,NY-1 |
---|
70 | do I=2,NX-1 |
---|
71 | if (T(i,j,1).le.TTRANS(iiglen)) then |
---|
72 | BTT(I,J,1,iiglen)=BAT1(iiglen)*exp(Q1(iiglen)/RGAS/TI1*T(I,J,1)) |
---|
73 | else |
---|
74 | BTT(I,J,1,iiglen)=BAT2(iiglen)*exp(Q2(iiglen)/RGAS/TI1*T(I,J,1)) |
---|
75 | endif |
---|
76 | end do |
---|
77 | end do |
---|
78 | |
---|
79 | ! boucle pour tous les autres BTT |
---|
80 | do K=NZ-1,1,-1 |
---|
81 | do J=2,NY-1 |
---|
82 | do I=2,NX-1 |
---|
83 | if (H(i,j).gt.1.) then |
---|
84 | if ((TETA(I,J,K).le.TTRANS(iiglen)).and. & |
---|
85 | (TETA(I,J,K+1).le.TTRANS(iiglen))) then |
---|
86 | BAT=BAT1(iiglen) |
---|
87 | Q=Q1(iiglen) |
---|
88 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
89 | AAT=max(-1.8,AAT) |
---|
90 | AAT=min(80.,AAT) |
---|
91 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
92 | SI1(I,J,K)=SSEC/(AAT+a4)*(E(K)**(AAT+a4)-E(K+1)**(AAT+a4)) |
---|
93 | SI2(I,J,K)=((E(K)**(AAT+a5))/(AAT+a5) & |
---|
94 | -(E(K+1)**(AAT+a5))/(AAT+a5) & |
---|
95 | -(E(K+1)**(AAT+a4))*E(K)+E(K+1)**(AAT+a5)) & |
---|
96 | * SSEC/(AAT+a4) |
---|
97 | BTT(I,J,K+1,iiglen)=BAT*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
98 | else if ((TETA(I,J,K).ge.TTRANS(iiglen)).and. & |
---|
99 | (TETA(I,J,K+1).ge.TTRANS(iiglen))) then |
---|
100 | BAT=BAT2(iiglen) |
---|
101 | Q=Q2(iiglen) |
---|
102 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
103 | AAT=max(-1.8,AAT) |
---|
104 | AAT=min(80.,AAT) |
---|
105 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
106 | SI1(I,J,K)=SSEC/(AAT+a4)*(E(K)**(AAT+a4)-E(K+1)**(AAT+a4)) |
---|
107 | SI2(I,J,K)=((E(K)**(AAT+a5))/(AAT+a5) & |
---|
108 | -(E(K+1)**(AAT+a5))/(AAT+a5) & |
---|
109 | -(E(K+1)**(AAT+a4))*E(K)+E(K+1)**(AAT+a5)) & |
---|
110 | * SSEC/(AAT+a4) |
---|
111 | BTT(I,J,K+1,iiglen)=BAT*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
112 | |
---|
113 | ! cas ou la temperature de la maille est en partie au dessus et au dessous |
---|
114 | ! de TTRANS(iiglen) |
---|
115 | else if ((TETA(I,J,K).LT.TTRANS(iiglen)).and. & |
---|
116 | (TETA(I,J,K+1).GT.TTRANS(iiglen))) then |
---|
117 | ! calcul de la profondeur de transition |
---|
118 | if (abs(TETA(I,J,K)-TETA(I,J,K+1)).lt.1.e-5) then |
---|
119 | ZETAT=(E(K)+E(K+1))/2. |
---|
120 | else |
---|
121 | ZETAT=E(K+1)-(TTRANS(iiglen)-TETA(I,J,K+1))/ & |
---|
122 | (TETA(I,J,K)-TETA(I,J,K+1))*DE |
---|
123 | endif |
---|
124 | |
---|
125 | ! integration entre ZETA2 et ZETAT, loi bat2 |
---|
126 | BAT=BAT2(iiglen) |
---|
127 | Q=Q2(iiglen) |
---|
128 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
129 | AAT=max(-1.8,AAT) |
---|
130 | AAT=min(80.,AAT) |
---|
131 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
132 | SI1(I,J,K)=SSEC/(AAT+a4)*(ZETAT**(AAT+a4)-E(K+1)**(AAT+a4)) |
---|
133 | SI2(I,J,K)=((ZETAT**(AAT+a5))/(AAT+a5) & |
---|
134 | -(E(K+1)**(AAT+a5))/(AAT+a5) & |
---|
135 | -(E(K+1)**(AAT+a4))*ZETAT+E(K+1)**(AAT+a5)) & |
---|
136 | * SSEC/(AAT+a4) |
---|
137 | BTT(I,J,K+1,iiglen)=BAT*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
138 | |
---|
139 | ! integration entre ZETAT et ZETA1, loi bat1 |
---|
140 | BAT=BAT1(iiglen) |
---|
141 | Q=Q1(iiglen) |
---|
142 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
143 | AAT=max(-1.8,AAT) |
---|
144 | AAT=min(80.,AAT) |
---|
145 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
146 | SI1(I,J,K)=SI1(I,J,K)+ & |
---|
147 | SSEC/(AAT+a4)*(E(K)**(AAT+a4)-ZETAT**(AAT+a4)) |
---|
148 | SI2(I,J,K)=SI2(I,J,K) & |
---|
149 | +((E(K)**(AAT+a5))/(AAT+a5) & |
---|
150 | -(ZETAT**(AAT+a5))/(AAT+a5) & |
---|
151 | -(ZETAT**(AAT+a4))*E(K)+ZETAT**(AAT+a5)) & |
---|
152 | * SSEC/(AAT+a4) |
---|
153 | |
---|
154 | |
---|
155 | ! deuxieme cas ou la temperature de la maille est en partie au dessus et |
---|
156 | ! au dessous de TTRANS(iiglen) |
---|
157 | else if ((TETA(I,J,K).GT.TTRANS(iiglen)).and. & |
---|
158 | (TETA(I,J,K+1).LT.TTRANS(iiglen))) then |
---|
159 | ! calcul de la profondeur de transition |
---|
160 | if (abs(TETA(I,J,K)-TETA(I,J,K+1)).lt.1.e-5) then |
---|
161 | ZETAT=(E(K)+E(K+1))/2. |
---|
162 | else |
---|
163 | ZETAT=E(K+1)-(TTRANS(iiglen)-TETA(I,J,K+1))/ & |
---|
164 | (TETA(I,J,K)-TETA(I,J,K+1))*DE |
---|
165 | endif |
---|
166 | |
---|
167 | ! integration entre ZETA2 et ZETAT, loi bat1 |
---|
168 | BAT=BAT1(iiglen) |
---|
169 | Q=Q1(iiglen) |
---|
170 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
171 | AAT=max(-1.8,AAT) |
---|
172 | AAT=min(80.,AAT) |
---|
173 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
174 | SI1(I,J,K)=SSEC/(AAT+a4)*(ZETAT**(AAT+a4)-E(K+1)**(AAT+a4)) |
---|
175 | SI2(I,J,K)=((ZETAT**(AAT+a5))/(AAT+a5) & |
---|
176 | -(E(K+1)**(AAT+a5))/(AAT+a5) & |
---|
177 | -(E(K+1)**(AAT+a4))*ZETAT+E(K+1)**(AAT+a5)) & |
---|
178 | * SSEC/(AAT+a4) |
---|
179 | BTT(I,J,K+1,iiglen)=BAT*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
180 | |
---|
181 | ! integration entre ZETAT et ZETA1, loi bat2 |
---|
182 | BAT=BAT2(iiglen) |
---|
183 | Q=Q2(iiglen) |
---|
184 | AAT=Q/RGAS/TI2(I,J,K)*(TETA(I,J,K+1)-TETA(I,J,K))/DE*E(K+1) |
---|
185 | AAT=max(-1.8,AAT) |
---|
186 | AAT=min(80.,AAT) |
---|
187 | SSEC=BAT/(E(K+1)**AAT)*exp((Q*TETA(I,J,K+1))/(RGAS*TI2(I,J,K))) |
---|
188 | SI1(I,J,K)=SI1(I,J,K)+ & |
---|
189 | SSEC/(AAT+a4)*(E(K)**(AAT+a4)-ZETAT**(AAT+a4)) |
---|
190 | SI2(I,J,K)=SI2(I,J,K) & |
---|
191 | +((E(K)**(AAT+a5))/(AAT+a5) & |
---|
192 | -(ZETAT**(AAT+a5))/(AAT+a5) & |
---|
193 | -(ZETAT**(AAT+a4))*E(K)+ZETAT**(AAT+a5)) & |
---|
194 | * SSEC/(AAT+a4) |
---|
195 | endif |
---|
196 | endif |
---|
197 | end do |
---|
198 | end do |
---|
199 | end do |
---|
200 | |
---|
201 | ! integration de SA et S2A |
---|
202 | |
---|
203 | do J=1,NY |
---|
204 | do I=1,NX |
---|
205 | SA(I,J,NZ,iiglen)=0. |
---|
206 | S2A(I,J,NZ,iiglen)=0. |
---|
207 | if (H(I,J).le.1.) BTT(I,J,NZ,iiglen)=BAT2(iiglen) |
---|
208 | end do |
---|
209 | end do |
---|
210 | |
---|
211 | do K=NZ-1,1,-1 |
---|
212 | do J=1,NY |
---|
213 | do I=1,NX |
---|
214 | if (H(i,j).gt.1.) then |
---|
215 | SA(I,J,K,iiglen)=SA(I,J,K+1,iiglen)-SI1(I,J,K) |
---|
216 | S2A(I,J,K,iiglen)=S2A(I,J,K+1,iiglen)+SI2(I,J,K)+ & |
---|
217 | SA(I,J,K+1,iiglen)*DE |
---|
218 | else |
---|
219 | SA(I,J,K,iiglen)=BAT2(iiglen)/a4*(1-E(K)**a4) |
---|
220 | S2A(I,J,K,iiglen)=BAT2(iiglen)/a4*(a4/a5-E(K)+E(K)**a5/a5) |
---|
221 | BTT(I,J,K,iiglen)=BAT2(iiglen) |
---|
222 | endif |
---|
223 | end do |
---|
224 | end do |
---|
225 | end do |
---|
226 | |
---|
227 | ! cas particulier des bords |
---|
228 | SA(:,1,:,iiglen)=SA(:,2,:,iiglen) |
---|
229 | S2A(:,1,:,iiglen)=S2A(:,2,:,iiglen) |
---|
230 | BTT(:,1,:,iiglen)=BTT(:,2,:,iiglen) |
---|
231 | SA(:,NY,:,iiglen)=SA(:,NY-1,:,iiglen) |
---|
232 | S2A(:,NY,:,iiglen)=S2A(:,NY-1,:,iiglen) |
---|
233 | BTT(:,NY,:,iiglen)=BTT(:,NY-1,:,iiglen) |
---|
234 | |
---|
235 | SA(1,:,:,iiglen)=SA(2,:,:,iiglen) |
---|
236 | S2A(1,:,:,iiglen)=S2A(2,:,:,iiglen) |
---|
237 | BTT(1,:,:,iiglen)=BTT(2,:,:,iiglen) |
---|
238 | SA(NX,:,:,iiglen)=SA(NX-1,:,:,iiglen) |
---|
239 | S2A(NX,:,:,iiglen)=S2A(NX-1,:,:,iiglen) |
---|
240 | BTT(NX,:,:,iiglen)=BTT(NX-1,:,:,iiglen) |
---|
241 | |
---|
242 | ! calcul des moyennes sur les mailles staggered |
---|
243 | |
---|
244 | do k=1,nz |
---|
245 | tab2D=SA(:,:,k,iiglen) |
---|
246 | |
---|
247 | call moy_mxmy(nx,ny,tab2D,tab_mx,tab_my) |
---|
248 | SA_mx(:,:,k,iglen)=tab_mx |
---|
249 | SA_my(:,:,k,iglen)=tab_my |
---|
250 | |
---|
251 | tab2D=S2A(:,:,k,iiglen) |
---|
252 | call moy_mxmy(nx,ny,tab2D,tab_mx,tab_my) |
---|
253 | S2A_mx(:,:,k,iglen)=tab_mx |
---|
254 | S2A_my(:,:,k,iglen)=tab_my |
---|
255 | end do |
---|
256 | |
---|
257 | ! attribution de debug_3D pour les sorties S2a |
---|
258 | ! loi 1 -> 190 dans description variable -> 90 dans debug_3D |
---|
259 | debug_3d(:,:,90+iiglen-1) = S2a(:,:,1,iiglen) |
---|
260 | |
---|
261 | END SUBROUTINE FLOWLAW |
---|