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

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

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