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

Last change on this file since 465 was 465, checked in by aquiquet, 4 months ago

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