source: trunk/SOURCES/flowlaw-0.3.f90 @ 25

Last change on this file since 25 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk
File size: 9.9 KB
Line 
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
31SUBROUTINE 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
259debug_3d(:,:,90+iiglen-1) = S2a(:,:,1,iiglen)
260
261END SUBROUTINE FLOWLAW
Note: See TracBrowser for help on using the repository browser.