source: trunk/SOURCES/deformation_mod_2lois.f90 @ 170

Last change on this file since 170 was 93, checked in by dumas, 8 years ago

First version with Schoof flux parameterisation at the grounding line. | New module furst_schoof_mod.f90 | New flag Schoof in grdline namelist (see in SOURCES/Fichiers-parametres/A-LBq15_param_list_Schoof.dat)
File size: 17.0 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)
113
114write(num_rep_42,*)'!___________________________________________________________' 
115write(num_rep_42,*)'! loi de deformation 1          module deformation_mod_2lois'
116write(num_rep_42,*)
117write(num_rep_42,*)
118write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)'
119write(num_rep_42,*)'! enhancement factor (sf)'
120write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :'
121write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)'
122write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :'
123write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)'
124write(num_rep_42,*)'!___________________________________________________________' 
125
126glen(1)   = exposant_1 
127ttrans(1) = temp_trans_1
128sf(1)     = enhanc_fact_1 
129Bat1(1)   = coef_cold_1
130Q1(1)     = Q_cold_1 
131Bat2(1)   = coef_warm_1
132Q2(1)     = Q_warm_1 
133
134
135! loi 2
136!--------------------------------------------------------------------------------
137 
138rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
139read(num_param,loidef_2)
140
141write(num_rep_42,*)'!___________________________________________________________' 
142write(num_rep_42,*)'! loi de deformation 2          module deformation_mod_2lois'
143write(num_rep_42,*)
144write(num_rep_42,*)
145write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)'
146write(num_rep_42,*)'! enhancement factor (sf)'
147write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :'
148write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)'
149write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :'
150write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)'
151write(num_rep_42,*)'!___________________________________________________________' 
152
153glen(2)   = exposant_2 
154ttrans(2) = temp_trans_2
155sf(2)     = enhanc_fact_2 
156Bat1(2)   = coef_cold_2
157Q1(2)     = Q_cold_2 
158Bat2(2)   = coef_warm_2
159Q2(2)     = Q_warm_2 
160
161! autre parametres ne changeant pas d'un run a l'autre
162RGAS=8.314
163
164! application des sf
165
166do iglen= n1poly,n2poly
167   bat1(iglen)=bat1(iglen)*sf(iglen)
168   bat2(iglen)=bat2(iglen)*sf(iglen)
169end do
170
171ddx(:,:,:)=0.
172ddy(:,:,:)=0.
173
174
175end subroutine init_deformation
176
177
178!--------------------------------------------------------------------
179subroutine flow_general
180
181  !$ USE OMP_LIB
182  implicit none
183
184  !$OMP PARALLEL
185  !$OMP WORKSHARE
186  teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:)
187  !$OMP END WORKSHARE
188
189  !$OMP DO COLLAPSE(2)
190  do k=nz-1,1,-1
191     do i=2,nx-1
192        do j=2,ny-1
193           ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* &
194                (273.15+(t(i,j,k+1)+t(i,j,k))/2.)
195        end do
196     end do
197  end do
198  !$OMP END DO
199  !$OMP END PARALLEL
200
201end subroutine flow_general
202!---------------------------------------------------------------------------------------
203subroutine flowlaw (iiglen)
204
205  !$ USE OMP_LIB
206
207  implicit none
208  integer ::  iiglen   !< compteur pour la boucle flowlaw
209  real :: a5 !< exposant de la loi de def +2
210  real :: a4 !< exposant de la loi de def +1
211  real :: ti1 !< utilise pour le calcul de btt a k=1
212  real :: bat !< prend la valeur de bat1 ou bat2 selon les cas
213  real :: q !< prend la valeur de q1 ou q2 selon les cas
214  real :: aat !< variable de travail =q*tetar(i,j,k)
215  real :: ssec !< variable de travail
216  real :: zetat !< variable de travail : profondeur ou t = trans(iiglen)
217  real,dimension(nx,ny,nz) :: si1 !< tableau de calcul
218  real,dimension(nx,ny,nz) :: si2 !< tableau de calcul
219  real,dimension(nx,ny) :: tab_mx
220  real,dimension(nx,ny) :: tab_my
221  real,dimension(nx,ny) :: tab2d
222
223!  real,dimension(nz) :: e          ! vertical coordinate in ice, scaled to h zeta
224  real :: de= 1./(nz-1)
225! variables openmp : 
226  !$  integer :: rang
227  !$ integer, dimension(:), allocatable :: tab_sync
228  !$ integer :: nb_taches
229
230
231  e(1)=0.
232  e(nz)=1.
233 
234  !$OMP PARALLEL PRIVATE(bat,q,aat,ssec,zetat)
235  !$ nb_taches = OMP_GET_NUM_THREADS()
236  !$OMP DO
237  do k=1,nz
238     if ((k.ne.1).and.(k.ne.nz)) e(k)=(k-1.)/(nz-1.)
239  end do
240  !$OMP END DO NOWAIT
241
242  ! exposant de la loi de deformation utilisee : glen(iiglen)
243  ! l'exposant correspondant a iiglen est defini dans deformation_mod
244  a5 = glen(iiglen) + 2 
245  a4 = glen(iiglen) + 1
246
247  !    boucle pour btt a k=1
248  ti1=273.15*273.15
249  !$OMP DO
250  do j=2,ny-1
251     do i=2,nx-1
252        if (t(i,j,1).le.ttrans(iiglen)) then
253           btt(i,j,1,iiglen)=bat1(iiglen)*exp(q1(iiglen)/rgas/ti1*t(i,j,1))
254        else
255           btt(i,j,1,iiglen)=bat2(iiglen)*exp(q2(iiglen)/rgas/ti1*t(i,j,1))
256        endif
257     end do
258  end do
259  !$OMP END DO
260
261  ! boucle pour tous les autres btt
262  !$OMP DO COLLAPSE(2)
263  do k=nz-1,1,-1
264        do j=2,ny-1 
265        do i=2,nx-1 
266           if (h(i,j).gt.1.) then
267              if ((teta(i,j,k).le.ttrans(iiglen)).and. &
268                   (teta(i,j,k+1).le.ttrans(iiglen))) then
269                 bat=bat1(iiglen)
270                 q=q1(iiglen)
271                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
272                 aat=max(-1.8,aat)
273                 aat=min(80.,aat)
274                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
275                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
276                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
277                      -(e(k+1)**(aat+a5))/(aat+a5) &
278                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
279                      * ssec/(aat+a4) 
280                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
281              else if ((teta(i,j,k).ge.ttrans(iiglen)).and. &
282                   (teta(i,j,k+1).ge.ttrans(iiglen))) then
283                 bat=bat2(iiglen)
284                 q=q2(iiglen)
285                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
286                 aat=max(-1.8,aat)
287                 aat=min(80.,aat)
288                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
289                 si1(i,j,k)=ssec/(aat+a4)*(e(k)**(aat+a4)-e(k+1)**(aat+a4))
290                 si2(i,j,k)=((e(k)**(aat+a5))/(aat+a5) &
291                      -(e(k+1)**(aat+a5))/(aat+a5) &
292                      -(e(k+1)**(aat+a4))*e(k)+e(k+1)**(aat+a5)) &
293                      * ssec/(aat+a4) 
294                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
295
296                 ! cas ou la temperature de la maille est en partie au dessus et au dessous
297                 ! de ttrans(iiglen)
298              else if ((teta(i,j,k).lt.ttrans(iiglen)).and. &
299                   (teta(i,j,k+1).gt.ttrans(iiglen))) then
300                 !         calcul de la profondeur de transition
301                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
302                    zetat=(e(k)+e(k+1))/2.
303                 else
304                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
305                         (teta(i,j,k)-teta(i,j,k+1))*de
306                 endif
307
308                 !         integration entre zeta2 et zetat, loi bat2
309                 bat=bat2(iiglen)
310                 q=q2(iiglen)
311                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
312                 aat=max(-1.8,aat)
313                 aat=min(80.,aat)
314                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
315                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
316                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
317                      -(e(k+1)**(aat+a5))/(aat+a5) &
318                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
319                      * ssec/(aat+a4)
320                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
321
322                 !         integration entre zetat et zeta1, loi bat1
323                 bat=bat1(iiglen)
324                 q=q1(iiglen)
325                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
326                 aat=max(-1.8,aat)
327                 aat=min(80.,aat)
328                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
329                 si1(i,j,k)=si1(i,j,k)+ &
330                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
331                 si2(i,j,k)=si2(i,j,k) &
332                      +((e(k)**(aat+a5))/(aat+a5) &
333                      -(zetat**(aat+a5))/(aat+a5) &
334                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
335                      * ssec/(aat+a4)
336
337
338                 ! deuxieme cas ou la temperature de la maille est en partie au dessus et
339                 ! au dessous de ttrans(iiglen)
340              else if ((teta(i,j,k).gt.ttrans(iiglen)).and. &
341                   (teta(i,j,k+1).lt.ttrans(iiglen))) then
342                 !         calcul de la profondeur de transition
343                 if (abs(teta(i,j,k)-teta(i,j,k+1)).lt.1.e-5) then
344                    zetat=(e(k)+e(k+1))/2.
345                 else
346                    zetat=e(k+1)-(ttrans(iiglen)-teta(i,j,k+1))/ &
347                         (teta(i,j,k)-teta(i,j,k+1))*de
348                 endif
349
350                 !         integration entre zeta2 et zetat, loi bat1
351                 bat=bat1(iiglen)
352                 q=q1(iiglen)
353                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
354                 aat=max(-1.8,aat)
355                 aat=min(80.,aat)
356                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
357                 si1(i,j,k)=ssec/(aat+a4)*(zetat**(aat+a4)-e(k+1)**(aat+a4))
358                 si2(i,j,k)=((zetat**(aat+a5))/(aat+a5) &
359                      -(e(k+1)**(aat+a5))/(aat+a5) &
360                      -(e(k+1)**(aat+a4))*zetat+e(k+1)**(aat+a5)) &
361                      * ssec/(aat+a4)
362                 btt(i,j,k+1,iiglen)=bat*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
363
364                 !         integration entre zetat et zeta1, loi bat2
365                 bat=bat2(iiglen)
366                 q=q2(iiglen)
367                 aat=q/rgas/ti2(i,j,k)*(teta(i,j,k+1)-teta(i,j,k))/de*e(k+1)
368                 aat=max(-1.8,aat)
369                 aat=min(80.,aat)
370                 ssec=bat/(e(k+1)**aat)*exp((q*teta(i,j,k+1))/(rgas*ti2(i,j,k)))
371                 si1(i,j,k)=si1(i,j,k)+ &
372                      ssec/(aat+a4)*(e(k)**(aat+a4)-zetat**(aat+a4))
373                 si2(i,j,k)=si2(i,j,k) &
374                      +((e(k)**(aat+a5))/(aat+a5) &
375                      -(zetat**(aat+a5))/(aat+a5) &
376                      -(zetat**(aat+a4))*e(k)+zetat**(aat+a5)) &
377                      * ssec/(aat+a4)
378              endif
379           endif
380        end do
381     end do
382  end do
383  !$OMP END DO NOWAIT
384 
385  !     integration de sa et s2a
386  !$OMP DO
387  do j=1,ny
388     do i=1,nx 
389        sa(i,j,nz,iiglen)=0.
390        s2a(i,j,nz,iiglen)=0.
391        if (h(i,j).le.1.) btt(i,j,nz,iiglen)=bat2(iiglen)
392     end do
393  end do
394  !$OMP END DO
395  !$OMP END PARALLEL
396 
397
398  ! Allocation et initialisation du tableau tab_sync
399  ! qui gere la synchronisation entre les differents threads
400  !$ allocate(tab_sync(0:nb_taches-1))
401  !$ tab_sync(0:nb_taches-1) = 1
402  !$OMP PARALLEL private(i,j,k,rang) shared(tab_sync)
403  !$ rang = OMP_GET_THREAD_NUM() 
404  do k=nz-1,1,-1
405       ! Synchronisation des threads
406     !$     if (rang /= 0) then
407     !$        do
408     !$OMP FLUSH(tab_sync)
409     !$           if (tab_sync(rang-1)>=tab_sync(rang)+1) exit
410     !$        enddo
411     !$OMP FLUSH(sa)
412     !$OMP FLUSH(s2a)
413     !$     endif
414     !$OMP DO SCHEDULE(STATIC)
415     do j=1,ny
416        do i=1,nx
417           if (h(i,j).gt.1.) then
418              sa(i,j,k,iiglen)=sa(i,j,k+1,iiglen)-si1(i,j,k)
419              s2a(i,j,k,iiglen)=s2a(i,j,k+1,iiglen)+si2(i,j,k)+ &
420                   sa(i,j,k+1,iiglen)*de
421           else
422              sa(i,j,k,iiglen)=bat2(iiglen)/a4*(1-e(k)**a4)
423              s2a(i,j,k,iiglen)=bat2(iiglen)/a4*(a4/a5-e(k)+e(k)**a5/a5)
424              btt(i,j,k,iiglen)=bat2(iiglen)
425           endif
426        end do
427     end do
428     !$OMP END DO NOWAIT
429!     ! Mis a jour du tableau tab_sync
430     !$     tab_sync(rang) = tab_sync(rang) + 1
431     !$OMP FLUSH(tab_sync,sa,s2a)
432  end do
433
434  !$OMP WORKSHARE
435  !     cas particulier des bords
436  sa(:,1,:,iiglen)=sa(:,2,:,iiglen)
437  s2a(:,1,:,iiglen)=s2a(:,2,:,iiglen)
438  btt(:,1,:,iiglen)=btt(:,2,:,iiglen)
439  sa(:,ny,:,iiglen)=sa(:,ny-1,:,iiglen)
440  s2a(:,ny,:,iiglen)=s2a(:,ny-1,:,iiglen)
441  btt(:,ny,:,iiglen)=btt(:,ny-1,:,iiglen)
442
443  sa(1,:,:,iiglen)=sa(2,:,:,iiglen)
444  s2a(1,:,:,iiglen)=s2a(2,:,:,iiglen)
445  btt(1,:,:,iiglen)=btt(2,:,:,iiglen)
446  sa(nx,:,:,iiglen)=sa(nx-1,:,:,iiglen)
447  s2a(nx,:,:,iiglen)=s2a(nx-1,:,:,iiglen)
448  btt(nx,:,:,iiglen)=btt(nx-1,:,:,iiglen)
449  !$OMP END WORKSHARE
450  !$OMP END PARALLEL
451 
452  ! calcul des moyennes sur les mailles staggered
453  do k=1,nz
454     tab2d=sa(:,:,k,iiglen)
455
456     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
457     sa_mx(:,:,k,iglen)=tab_mx
458     sa_my(:,:,k,iglen)=tab_my
459
460     tab2d=s2a(:,:,k,iiglen)
461     call moy_mxmy(nx,ny,tab2d,tab_mx,tab_my)
462     s2a_mx(:,:,k,iglen)=tab_mx
463     s2a_my(:,:,k,iglen)=tab_my
464  end do
465 
466! attribution de debug_3d pour les sorties s2a
467! loi 1 -> 190 dans description variable -> 90 dans debug_3d
468!debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen)
469
470end subroutine flowlaw
471
472
473end module deformation_mod_2lois
474
Note: See TracBrowser for help on using the repository browser.