source: trunk/SOURCES/deformation_mod_2lois.f90 @ 334

Last change on this file since 334 was 224, checked in by dumas, 5 years ago

format of the parameter file created by the code, possibility to use it in reading
File size: 16.4 KB
Line 
1!> \file deformation_mod_2lois.f90
2!! C'est ce module qui doit etre selectionne pour faire le calcul de la deformation de la glace 
3!! en utilisant une loi de deformation polynomiale
4!! a deux composantes (habituellement n=3 + n=1)
5!!
6!<
7
8!> \namespace deformation_mod_2lois
9!! C'est ce module qui doit etre selectionne pour faire le calcul de la
10!! deformation de la glace en utilisant une loi de deformation polynomiale
11!! a deux composantes (habituellement n=3 + n=1)
12!! \author CatRitz
13!! \date decmebre 2008
14!! @note Used modules
15!! @note  - use module3D_phy
16!! @note  Contient :
17!! @note  - les declarations des tableaux (etait avant dans deform_declar)
18!! en fait la declaration est independante du nombre d'elements de la loi
19!! et peut etre simplement copie pour un autre nombre
20!! @note  - lecture des valeurs utilisees par namelist
21!! N'est valable que pour la loi a 2 composants. Pour un nombre different de composants,
22!! le recopier et modifier
23!! @note  -  les parameters n1poly et n2poly
24!! @note  -  init_deformation en ajoutant ou supprimant des blocs de namelist
25!!
26!<
27module deformation_mod_2lois
28
29use module3d_phy
30implicit none
31
32
33! declarations ne dependant pas du nombre de lois
34
35real, dimension(nx,ny,nz)   :: teta           !< teta = t - tpmp
36real, dimension(nx,ny,nz-1) :: ti2            !< calcule dans flow_general
37real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves)
38real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice)
39
40! declaration des lois
41! la valeur de n1poly et de n2poly determine le nombre de lois additionnees
42
43integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2
44integer, parameter :: n2poly=2
45
46
47! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly)
48
49real, dimension(n1poly:n2poly) :: glen             !< l'exposant
50real, dimension(n1poly:n2poly) :: ttrans           !< la temperature de transition
51real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law
52
53! Pour chaque loi (meme exposant), il y a deux domaines avec une temperature de transition
54
55! 1- temperature inferieure a Ttrans
56real, dimension(n1poly:n2poly) :: Q1               !< energies d'activation (temp < ttrans) 
57real, dimension(n1poly:n2poly) :: Bat1             !< coefficient           (temp < ttrans)
58         
59! 2- temperature superieure a Ttrans
60real, dimension(n1poly:n2poly) :: Q2               !< energies d'activation (temp > ttrans)
61real, dimension(n1poly:n2poly) :: Bat2             !< coefficient           (temp > ttrans)
62
63
64! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ...
65
66real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere
67real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses
68real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles >
69real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^
70real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses
71real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles >
72real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^
73real, dimension(nx,ny,n1poly:n2poly) :: ddx     !< sert dans le calcul de conserv. masse
74real, dimension(nx,ny,n1poly:n2poly) :: ddy     !< sert dans le calcul de conserv. masse
75
76contains
77
78!-------------------------------------------------------------------------------------------
79
80!> SUBROUTINE: init_deformation
81!! Routine qui lit les variables de deformation
82!>
83subroutine init_deformation 
84
85real :: exposant_1
86real :: temp_trans_1
87real :: enhanc_fact_1
88real :: coef_cold_1
89real :: Q_cold_1
90real :: coef_warm_1
91real :: Q_warm_1
92
93real :: exposant_2
94real :: temp_trans_2
95real :: enhanc_fact_2
96real :: coef_cold_2
97real :: Q_cold_2
98real :: coef_warm_2
99real :: Q_warm_2
100
101namelist/loidef_1/exposant_1,temp_trans_1,enhanc_fact_1,   & 
102                  coef_cold_1,Q_cold_1,coef_warm_1,Q_warm_1
103
104namelist/loidef_2/exposant_2,temp_trans_2,enhanc_fact_2,   & 
105                  coef_cold_2,Q_cold_2,coef_warm_2,Q_warm_2
106
107
108! loi 1
109!--------------------------------------------------------------------------------
110
111rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
112read(num_param,loidef_1)
113write(num_rep_42,loidef_1)
114
115glen(1)   = exposant_1 
116ttrans(1) = temp_trans_1
117sf(1)     = enhanc_fact_1 
118Bat1(1)   = coef_cold_1
119Q1(1)     = Q_cold_1 
120Bat2(1)   = coef_warm_1
121Q2(1)     = Q_warm_1 
122
123
124! loi 2
125!--------------------------------------------------------------------------------
126 
127rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
128read(num_param,loidef_2)
129write(num_rep_42,loidef_2)
130
131write(num_rep_42,*)'!___________________________________________________________' 
132write(num_rep_42,*)'! loi de deformation            module deformation_mod_2lois'
133write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)'
134write(num_rep_42,*)'! enhancement factor (sf)'
135write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :'
136write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)'
137write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :'
138write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)'
139write(num_rep_42,*)'!___________________________________________________________' 
140
141glen(2)   = exposant_2 
142ttrans(2) = temp_trans_2
143sf(2)     = enhanc_fact_2 
144Bat1(2)   = coef_cold_2
145Q1(2)     = Q_cold_2 
146Bat2(2)   = coef_warm_2
147Q2(2)     = Q_warm_2 
148
149! autre parametres ne changeant pas d'un run a l'autre
150RGAS=8.314
151
152! application des sf
153
154do iglen= n1poly,n2poly
155   bat1(iglen)=bat1(iglen)*sf(iglen)
156   bat2(iglen)=bat2(iglen)*sf(iglen)
157end do
158
159ddx(:,:,:)=0.
160ddy(:,:,:)=0.
161
162
163end subroutine init_deformation
164
165
166!--------------------------------------------------------------------
167subroutine flow_general
168
169  !$ USE OMP_LIB
170  implicit none
171
172  !$OMP PARALLEL
173  !$OMP WORKSHARE
174  teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:)
175  !$OMP END WORKSHARE
176
177  !$OMP DO COLLAPSE(2)
178  do k=nz-1,1,-1
179     do i=2,nx-1
180        do j=2,ny-1
181           ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* &
182                (273.15+(t(i,j,k+1)+t(i,j,k))/2.)
183        end do
184     end do
185  end do
186  !$OMP END DO
187  !$OMP END PARALLEL
188
189end subroutine flow_general
190!---------------------------------------------------------------------------------------
191subroutine flowlaw (iiglen)
192
193  !$ USE OMP_LIB
194
195  implicit none
196  integer ::  iiglen   !< compteur pour la boucle flowlaw
197  real :: a5 !< exposant de la loi de def +2
198  real :: a4 !< exposant de la loi de def +1
199  real :: ti1 !< utilise pour le calcul de btt a k=1
200  real :: bat !< prend la valeur de bat1 ou bat2 selon les cas
201  real :: q !< prend la valeur de q1 ou q2 selon les cas
202  real :: aat !< variable de travail =q*tetar(i,j,k)
203  real :: ssec !< variable de travail
204  real :: zetat !< variable de travail : profondeur ou t = trans(iiglen)
205  real,dimension(nx,ny,nz) :: si1 !< tableau de calcul
206  real,dimension(nx,ny,nz) :: si2 !< tableau de calcul
207  real,dimension(nx,ny) :: tab_mx
208  real,dimension(nx,ny) :: tab_my
209  real,dimension(nx,ny) :: tab2d
210
211!  real,dimension(nz) :: e          ! vertical coordinate in ice, scaled to h zeta
212  real :: de= 1./(nz-1)
213! variables openmp : 
214  !$  integer :: rang
215  !$ integer, dimension(:), allocatable :: tab_sync
216  !$ integer :: nb_taches
217
218
219  e(1)=0.
220  e(nz)=1.
221 
222  !$OMP PARALLEL PRIVATE(bat,q,aat,ssec,zetat)
223  !$ nb_taches = OMP_GET_NUM_THREADS()
224  !$OMP DO
225  do k=1,nz
226     if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.)
227  end do
228  !$OMP END DO NOWAIT
229
230  ! exposant de la loi de deformation utilisee : glen(iiglen)
231  ! l'exposant correspondant a iiglen est defini dans deformation_mod
232  a5 = glen(iiglen) + 2 
233  a4 = glen(iiglen) + 1
234
235  !    boucle pour btt a k=1
236  ti1=273.15*273.15
237  !$OMP DO
238  do j=2,ny-1
239     do i=2,nx-1
240        if (t(i,j,1).le.ttrans(iiglen)) then
241           btt(i,j,1,iiglen)=bat1(iiglen)*exp(q1(iiglen)/rgas/ti1*t(i,j,1))
242        else
243           btt(i,j,1,iiglen)=bat2(iiglen)*exp(q2(iiglen)/rgas/ti1*t(i,j,1))
244        endif
245     end do
246  end do
247  !$OMP END DO
248
249  ! boucle pour tous les autres btt
250  !$OMP DO COLLAPSE(2)
251  do k=nz-1,1,-1
252        do j=2,ny-1 
253        do i=2,nx-1 
254           if (h(i,j).gt.1.) then
255              if ((teta(i,j,k).le.ttrans(iiglen)).and. &
256                   (teta(i,j,k+1).le.ttrans(iiglen))) then
257                 bat=bat1(iiglen)
258                 q=q1(iiglen)
259                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
260                 aat=max(-1.8,aat)
261                 aat=min(80.,aat)
262                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
263                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
264                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
265                      -(e(k+1)**(aat+a5))/(aat+a5) &
266                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
267                      * ssec/(aat+a4) 
268                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
269              else if ((teta(i,j,k).ge.ttrans(iiglen)).and. &
270                   (teta(i,j,k+1).ge.ttrans(iiglen))) then
271                 bat=bat2(iiglen)
272                 q=q2(iiglen)
273                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
274                 aat=max(-1.8,aat)
275                 aat=min(80.,aat)
276                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
277                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
278                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
279                      -(e(k+1)**(aat+a5))/(aat+a5) &
280                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
281                      * ssec/(aat+a4) 
282                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
283
284                 ! cas ou la temperature de la maille est en partie au dessus et au dessous
285                 ! de ttrans(iiglen)
286              else if ((teta(i,j,k).lt.ttrans(iiglen)).and. &
287                   (teta(i,j,k+1).gt.ttrans(iiglen))) then
288                 !         calcul de la profondeur de transition
289                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
290                    zetat=(e(k)+e(k+1))/2.
291                 else
292                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
293                         (teta(i,j,k)-teta(i,j,k+1))*de
294                 endif
295
296                 !         integration entre zeta2 et zetat, loi bat2
297                 bat=bat2(iiglen)
298                 q=q2(iiglen)
299                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
300                 aat=max(-1.8,aat)
301                 aat=min(80.,aat)
302                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
303                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
304                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
305                      -(e(k+1)**(aat+a5))/(aat+a5) &
306                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
307                      * ssec/(aat+a4)
308                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
309
310                 !         integration entre zetat et zeta1, loi bat1
311                 bat=bat1(iiglen)
312                 q=q1(iiglen)
313                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
314                 aat=max(-1.8,aat)
315                 aat=min(80.,aat)
316                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
317                 si1(i,j,k)=si1(i,j,k)+ &
318                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
319                 si2(i,j,k)=si2(i,j,k) &
320                      +((e(k)**(aat+a5))/(aat+a5) &
321                      -(zetat**(aat+a5))/(aat+a5) &
322                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
323                      * ssec/(aat+a4)
324
325
326                 ! deuxieme cas ou la temperature de la maille est en partie au dessus et
327                 ! au dessous de ttrans(iiglen)
328              else if ((teta(i,j,k).gt.ttrans(iiglen)).and. &
329                   (teta(i,j,k+1).lt.ttrans(iiglen))) then
330                 !         calcul de la profondeur de transition
331                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
332                    zetat=(e(k)+e(k+1))/2.
333                 else
334                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
335                         (teta(i,j,k)-teta(i,j,k+1))*de
336                 endif
337
338                 !         integration entre zeta2 et zetat, loi bat1
339                 bat=bat1(iiglen)
340                 q=q1(iiglen)
341                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
342                 aat=max(-1.8,aat)
343                 aat=min(80.,aat)
344                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
345                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
346                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
347                      -(e(k+1)**(aat+a5))/(aat+a5) &
348                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
349                      * ssec/(aat+a4)
350                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
351
352                 !         integration entre zetat et zeta1, loi bat2
353                 bat=bat2(iiglen)
354                 q=q2(iiglen)
355                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
356                 aat=max(-1.8,aat)
357                 aat=min(80.,aat)
358                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
359                 si1(i,j,k)=si1(i,j,k)+ &
360                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
361                 si2(i,j,k)=si2(i,j,k) &
362                      +((e(k)**(aat+a5))/(aat+a5) &
363                      -(zetat**(aat+a5))/(aat+a5) &
364                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
365                      * ssec/(aat+a4)
366              endif
367           endif
368        end do
369     end do
370  end do
371  !$OMP END DO NOWAIT
372 
373  !     integration de sa et s2a
374  !$OMP DO
375  do j=1,ny
376     do i=1,nx 
377        sa(i,j,nz,iiglen)=0.
378        s2a(i,j,nz,iiglen)=0.
379        if (h(i,j).le.1.) btt(i,j,nz,iiglen)=bat2(iiglen)
380     end do
381  end do
382  !$OMP END DO
383  !$OMP END PARALLEL
384 
385
386  ! Allocation et initialisation du tableau tab_sync
387  ! qui gere la synchronisation entre les differents threads
388  !$ allocate(tab_sync(0:nb_taches-1))
389  !$ tab_sync(0:nb_taches-1) = 1
390  !$OMP PARALLEL private(i,j,k,rang) shared(tab_sync)
391  !$ rang = OMP_GET_THREAD_NUM() 
392  do k=nz-1,1,-1
393       ! Synchronisation des threads
394     !$     if (rang /= 0) then
395     !$        do
396     !$OMP FLUSH(tab_sync)
397     !$           if (tab_sync(rang-1)>=tab_sync(rang)+1) exit
398     !$        enddo
399     !$OMP FLUSH(sa)
400     !$OMP FLUSH(s2a)
401     !$     endif
402     !$OMP DO SCHEDULE(STATIC)
403     do j=1,ny
404        do i=1,nx
405           if (h(i,j).gt.1.) then
406              sa(i,j,k,iiglen)=sa(i,j,k+1,iiglen)-si1(i,j,k)
407              s2a(i,j,k,iiglen)=s2a(i,j,k+1,iiglen)+si2(i,j,k)+ &
408                   sa(i,j,k+1,iiglen)*de
409           else
410              sa(i,j,k,iiglen)=bat2(iiglen)/a4*(1-e(k)**a4)
411              s2a(i,j,k,iiglen)=bat2(iiglen)/a4*(a4/a5-e(k)+e(k)**a5/a5)
412              btt(i,j,k,iiglen)=bat2(iiglen)
413           endif
414        end do
415     end do
416     !$OMP END DO NOWAIT
417!     ! Mis a jour du tableau tab_sync
418     !$     tab_sync(rang) = tab_sync(rang) + 1
419     !$OMP FLUSH(tab_sync,sa,s2a)
420  end do
421
422  !$OMP WORKSHARE
423  !     cas particulier des bords
424  sa(:,1,:,iiglen)=sa(:,2,:,iiglen)
425  s2a(:,1,:,iiglen)=s2a(:,2,:,iiglen)
426  btt(:,1,:,iiglen)=btt(:,2,:,iiglen)
427  sa(:,ny,:,iiglen)=sa(:,ny-1,:,iiglen)
428  s2a(:,ny,:,iiglen)=s2a(:,ny-1,:,iiglen)
429  btt(:,ny,:,iiglen)=btt(:,ny-1,:,iiglen)
430
431  sa(1,:,:,iiglen)=sa(2,:,:,iiglen)
432  s2a(1,:,:,iiglen)=s2a(2,:,:,iiglen)
433  btt(1,:,:,iiglen)=btt(2,:,:,iiglen)
434  sa(nx,:,:,iiglen)=sa(nx-1,:,:,iiglen)
435  s2a(nx,:,:,iiglen)=s2a(nx-1,:,:,iiglen)
436  btt(nx,:,:,iiglen)=btt(nx-1,:,:,iiglen)
437  !$OMP END WORKSHARE
438  !$OMP END PARALLEL
439 
440  ! calcul des moyennes sur les mailles staggered
441  do k=1,nz
442     tab2d=sa(:,:,k,iiglen)
443
444     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
445     sa_mx(:,:,k,iglen)=tab_mx
446     sa_my(:,:,k,iglen)=tab_my
447
448     tab2d=s2a(:,:,k,iiglen)
449     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
450     s2a_mx(:,:,k,iglen)=tab_mx
451     s2a_my(:,:,k,iglen)=tab_my
452  end do
453 
454! attribution de debug_3d pour les sorties s2a
455! loi 1 -> 190 dans description variable -> 90 dans debug_3d
456!debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen)
457
458end subroutine flowlaw
459
460
461end module deformation_mod_2lois
462
Note: See TracBrowser for help on using the repository browser.