source: branches/GRISLIv3/SOURCES/deformation_mod_2lois.f90 @ 397

Last change on this file since 397 was 397, checked in by dumas, 13 months ago

use only in module furst_schoof_mod
File size: 17.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
29  use module3d_phy, only:nx,ny,nz!,e,num_param,num_rep_42,rgas,T,iglen,tpmp,H
30
31  implicit none
32
33
34  ! declarations ne dependant pas du nombre de lois
35
36  real, dimension(nx,ny,nz)   :: teta           !< teta = t - tpmp
37  real, dimension(nx,ny,nz-1) :: ti2            !< calcule dans flow_general
38  real, dimension(nx,ny,nz)   :: visc           !< viscosite de la glace (pour les shelves)
39  real, dimension(nz)         :: edecal         !< travail (decalage de E de 1 indice)
40
41  ! declaration des lois
42  ! la valeur de n1poly et de n2poly determine le nombre de lois additionnees
43
44  integer, parameter :: n1poly=1     !< 2 lois numerotees de 1 a 2
45  integer, parameter :: n2poly=2
46
47
48  ! Tous les parametres de la loi sous forme de tableau (n1poly:n2poly)
49
50  real, dimension(n1poly:n2poly) :: glen             !< l'exposant
51  real, dimension(n1poly:n2poly) :: ttrans           !< la temperature de transition
52  real, dimension(n1poly:n2poly) :: sf               !< softening factor for flow law
53
54  ! Pour chaque loi (meme exposant), il y a deux domaines avec une temperature de transition
55
56  ! 1- temperature inferieure a Ttrans
57  real, dimension(n1poly:n2poly) :: Q1               !< energies d'activation (temp < ttrans) 
58  real, dimension(n1poly:n2poly) :: Bat1             !< coefficient           (temp < ttrans)
59
60  ! 2- temperature superieure a Ttrans
61  real, dimension(n1poly:n2poly) :: Q2               !< energies d'activation (temp > ttrans)
62  real, dimension(n1poly:n2poly) :: Bat2             !< coefficient           (temp > ttrans)
63
64
65  ! declaration des tableaux utilises dans le calcul des vitesses, viscosites, ...
66
67  real, dimension(nx,ny,nz,n1poly:n2poly) :: Btt     !< coefficient au point considere
68  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa      !< sert dans l'integration des vitesses
69  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_mx   !< Sa sur demi mailles >
70  real, dimension(nx,ny,nz,n1poly:n2poly) :: Sa_my   !< Sa sur demi mailles ^
71  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a     !< sert dans l'integration des vitesses
72  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_mx  !< S2a sur demi mailles >
73  real, dimension(nx,ny,nz,n1poly:n2poly) :: S2a_my  !< S2a sur demi mailles ^
74  real, dimension(nx,ny,n1poly:n2poly) :: ddx     !< sert dans le calcul de conserv. masse
75  real, dimension(nx,ny,n1poly:n2poly) :: ddy     !< sert dans le calcul de conserv. masse
76
77contains
78
79  !-------------------------------------------------------------------------------------------
80
81  !> SUBROUTINE: init_deformation
82  !! Routine qui lit les variables de deformation
83  !>
84  subroutine init_deformation
85   
86    use module3d_phy, only:num_param,num_rep_42,rgas,iglen
87
88    implicit none
89   
90    real :: exposant_1
91    real :: temp_trans_1
92    real :: enhanc_fact_1
93    real :: coef_cold_1
94    real :: Q_cold_1
95    real :: coef_warm_1
96    real :: Q_warm_1
97
98    real :: exposant_2
99    real :: temp_trans_2
100    real :: enhanc_fact_2
101    real :: coef_cold_2
102    real :: Q_cold_2
103    real :: coef_warm_2
104    real :: Q_warm_2
105
106    namelist/loidef_1/exposant_1,temp_trans_1,enhanc_fact_1,   & 
107         coef_cold_1,Q_cold_1,coef_warm_1,Q_warm_1
108
109    namelist/loidef_2/exposant_2,temp_trans_2,enhanc_fact_2,   & 
110         coef_cold_2,Q_cold_2,coef_warm_2,Q_warm_2
111
112
113    ! loi 1
114    !--------------------------------------------------------------------------------
115
116    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
117    read(num_param,loidef_1)
118    write(num_rep_42,loidef_1)
119
120    glen(1)   = exposant_1 
121    ttrans(1) = temp_trans_1
122    sf(1)     = enhanc_fact_1 
123    Bat1(1)   = coef_cold_1
124    Q1(1)     = Q_cold_1 
125    Bat2(1)   = coef_warm_1
126    Q2(1)     = Q_warm_1 
127
128
129    ! loi 2
130    !--------------------------------------------------------------------------------
131
132    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
133    read(num_param,loidef_2)
134    write(num_rep_42,loidef_2)
135
136    write(num_rep_42,*)'!___________________________________________________________' 
137    write(num_rep_42,*)'! loi de deformation            module deformation_mod_2lois'
138    write(num_rep_42,*)'! exposant (glen), temperature de transition (ttrans)'
139    write(num_rep_42,*)'! enhancement factor (sf)'
140    write(num_rep_42,*)'! pour les temperatures inf. a Temp_trans :'
141    write(num_rep_42,*)'!                            coef_cold (Bat1) et Q_cold (Q1)'
142    write(num_rep_42,*)'! pour les temperatures sup. a Temp_trans :'
143    write(num_rep_42,*)'!                            coef_warm (Bat2) et Q_warm (Q2)'
144    write(num_rep_42,*)'!___________________________________________________________' 
145
146    glen(2)   = exposant_2 
147    ttrans(2) = temp_trans_2
148    sf(2)     = enhanc_fact_2 
149    Bat1(2)   = coef_cold_2
150    Q1(2)     = Q_cold_2 
151    Bat2(2)   = coef_warm_2
152    Q2(2)     = Q_warm_2 
153
154    ! autre parametres ne changeant pas d'un run a l'autre
155    RGAS=8.314
156
157    ! application des sf
158
159    do iglen= n1poly,n2poly
160       bat1(iglen)=bat1(iglen)*sf(iglen)
161       bat2(iglen)=bat2(iglen)*sf(iglen)
162    end do
163
164    ddx(:,:,:)=0.
165    ddy(:,:,:)=0.
166
167
168  end subroutine init_deformation
169
170
171  !--------------------------------------------------------------------
172  subroutine flow_general
173   
174    use module3d_phy, only:T,tpmp
175    !$ USE OMP_LIB
176   
177    implicit none
178
179    integer :: i,j,k
180
181    !$OMP PARALLEL
182    !$OMP WORKSHARE
183    teta(:,:,:) = t(:,:,1:nz)-tpmp(:,:,:)
184    !$OMP END WORKSHARE
185
186    !$OMP DO COLLAPSE(2)
187    do k=nz-1,1,-1
188       do i=2,nx-1
189          do j=2,ny-1
190             ti2(i,j,k) = (273.15+(tpmp(i,j,k+1)+tpmp(i,j,k))/2.)* &
191                  (273.15+(t(i,j,k+1)+t(i,j,k))/2.)
192          end do
193       end do
194    end do
195    !$OMP END DO
196    !$OMP END PARALLEL
197
198  end subroutine flow_general
199  !---------------------------------------------------------------------------------------
200  subroutine flowlaw (iiglen)
201   
202    use module3d_phy, only:e,T,rgas,H,iglen
203    !$ USE OMP_LIB
204
205    implicit none
206
207    integer :: i,j,k
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
470  end subroutine flowlaw
471
472
473end module deformation_mod_2lois
474
Note: See TracBrowser for help on using the repository browser.