Changeset 74


Ignore:
Timestamp:
06/24/16 09:19:59 (8 years ago)
Author:
dumas
Message:

OpenMP parallelization in remplimat-shelves-tabTu, litho and taubed

Location:
trunk/SOURCES
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90

    r65 r74  
    4141!----------------------------------------------------------------------- 
    4242 
     43!$ USE OMP_LIB 
    4344use module3d_phy 
    4445use remplimat_declar        ! les tableaux tuij, .... 
    … …  
    104105 
    105106!----------------------------- 
     107!$OMP PARALLEL 
     108!$OMP WORKSHARE 
    106109Tu(:,:,:,:) = 0. ; Tv(:,:,:,:) = 0. ; Su(:,:,:,:) = 0. ; Sv(:,:,:,:) = 0. 
    107110opposx(:,:) = 0. ; opposy(:,:) = 0. 
    … …  
    117120Hmx_oppos(:,:) = Hmx(:,:) 
    118121Hmy_oppos(:,:) = Hmy(:,:) 
    119  
     122!$OMP END WORKSHARE 
     123!$OMP END PARALLEL 
    120124 
    121125! calcul de Hmx_oppos et Hmy_oppos dans le cas de fronts 
    … …  
    164168 call limit_frotm 
    165169else 
     170   !$OMP PARALLEL 
     171   !$OMP WORKSHARE 
    166172   frotmx(:,:)= betamx(:,:) 
    167173   frotmy(:,:)= betamy(:,:) 
     174   !$OMP END WORKSHARE 
     175   !$OMP END PARALLEL 
    168176end if 
    169177 
    … …  
    272280! pour avoir la diagonale=1 a faire imperativement APRES okmat 
    273281!-------------------------------------------------------------- 
    274  
     282!$OMP PARALLEL 
     283!$OMP DO PRIVATE(scal) 
    275284do j=ny1,ny2 
    276285   do i=nx1,nx2 
    … …  
    319328   end do 
    320329end do 
    321  
     330!$OMP END DO 
     331!$OMP END PARALLEL 
    322332 
    323333! mise a identite des noeuds fantomes 
    … …  
    329339 
    330340 
    331  
    332 debug_3d(:,:,35) = opposx(:,:) 
    333 debug_3d(:,:,36) = opposy(:,:) 
     341!debug_3d(:,:,35) = opposx(:,:) 
     342!debug_3d(:,:,36) = opposy(:,:) 
    334343 
    335344 
    … …  
    349358 
    350359! remise a 0 des noeuds fantomes 
    351  
     360!$OMP PARALLEL 
     361!$OMP WORKSHARE 
    352362where ((Hmx(nx1:nx2,ny1:ny2).le.1.).and.(flgzmx(nx1:nx2,ny1:ny2))) 
    353363   uxnew(nx1:nx2,ny1:ny2)=0. 
    … …  
    357367   uynew(nx1:nx2,ny1:ny2)=0. 
    358368end where 
    359  
     369!$OMP END WORKSHARE 
     370!$OMP END PARALLEL 
    360371 
    361372! call graphique_L2(kl,ku,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))    
    … …  
    404415 
    405416subroutine limit_frotm 
     417 
     418!$USE OMP_LIB 
    406419!------------------------- 
    407420! limite la valeur de frotmx a betamax 
    … …  
    411424if (itracebug.eq.1) write(Num_tracebug,*)'betamax = ',betamax 
    412425 
     426!$OMP PARALLEL 
     427!$OMP WORKSHARE 
    413428where (flgzmx(:,:)) 
    414429   frotmx(:,:)=min(abs(betamx(:,:)),betamax_2d(:,:)) 
    … …  
    422437   frotmy(:,:)=0 
    423438end where 
     439!$OMP END WORKSHARE 
     440!$OMP END PARALLEL 
    424441 
    425442end subroutine limit_frotm 
    … …  
    433450 
    434451subroutine rempli_Tuij                                    ! appelle tous les cas 
     452 
     453!$ USE OMP_LIB 
    435454 
    436455! imx(i,j)=0 -> ne pas traiter cette vitesse 
    … …  
    464483count_line=1                                 ! pour numeroter les lignes 
    465484!------------------------------------------------------------------ 
    466  
     485!$OMP PARALLEL 
     486!$OMP DO 
    467487lignes_UV: do j=ny1,ny2  
    468488   do i=nx1,nx2 
    … …  
    470490      if (i.gt.nx1) then                      ! Vitesses U 
    471491         ligu_L2(i,j)=count_line 
     492         !$OMP ATOMIC 
    472493         count_line=count_line+1 
    473494           
    … …  
    480501 
    481502         case(1)                              ! vitesse imposee 
    482             call vel_U_presc      
     503            call vel_U_presc(uxprec(i,j),Tu(i,j,0,0),opposx(i,j))     
    483504 
    484505         case(2)                              ! cas general 
    485             call vel_U_general 
     506            call vel_U_general(frotmx(i,j),dx2,pvi(i,j),pvi(i-1,j),pvm(i,j),pvm(i,j+1),sdx(i,j),rog,hmx_oppos(i,j), & 
     507            Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
    486508 
    487509         case(-1)                             ! bord sud 
    488             call vel_U_South           
     510            call vel_U_South(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))           
    489511 
    490512         case(-2)                             ! bord Est 
    491             call vel_U_East 
     513            call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
    492514        
    493515         case(-3)                             ! bord nord 
    494             call vel_U_North    
     516            call vel_U_North(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))    
    495517 
    496518         case(-4)                             ! bord West 
    497             call vel_U_West 
     519            call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j)) 
    498520 
    499521         case(-12)                            ! coin SE 
    500             call vel_U_SE 
     522            call vel_U_SE(Tu(i,j,:,:),opposx(i,j)) 
    501523 
    502524         case(-23)                            ! coin NE 
    503             call vel_U_NE 
     525            call vel_U_NE(Tu(i,j,:,:),opposx(i,j)) 
    504526 
    505527         case(-34)                            ! coin NW 
    506             call vel_U_NW 
     528            call vel_U_NW(Tu(i,j,:,:),opposx(i,j)) 
    507529 
    508530         case(-41)                            ! coin SW 
    509             call vel_U_SW 
     531            call vel_U_SW(Tu(i,j,:,:),opposx(i,j)) 
    510532 
    511533         end select case_imx 
    … …  
    515537 
    516538         ligv_L2(i,j)=count_line 
     539         !$OMP ATOMIC 
    517540         count_line=count_line+1 
    518541 
    … …  
    525548 
    526549         case(1)                              ! vitesse imposee 
    527             call vel_V_presc 
     550            call vel_V_presc(uyprec(i,j),Sv(i,j,0,0),opposy(i,j)) 
    528551             
    529552         case(2)                              ! cas general 
    530             call vel_V_general 
     553            call vel_V_general(frotmy(i,j),dx2,pvi(i,j),pvi(i,j-1),pvm(i,j),pvm(i+1,j),sdy(i,j),rog,hmy_oppos(i,j), & 
     554            Su(i,j,:,:),Sv(i,j,:,:),opposy(i,j)) 
    531555             
    532556         case(-1)                             ! bord sud 
    533             call vel_V_South              
     557            call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    534558             
    535559         case(-2)                             ! bord Est 
    536             call vel_V_East 
     560            call vel_V_East(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    537561             
    538562         case(-3)                             ! bord nord 
    539             call vel_V_North           
     563            call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    540564             
    541565         case(-4)                             ! bord West 
    542             call vel_V_West 
     566            call vel_V_West(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j)) 
    543567             
    544568         case(-12)                            ! coin SE 
    545             call vel_V_SE 
     569            call vel_V_SE(Sv(i,j,:,:),opposy(i,j))  
    546570             
    547571         case(-23)                            ! coin NE 
    548             call vel_V_NE 
     572            call vel_V_NE(Sv(i,j,:,:),opposy(i,j)) 
    549573             
    550574         case(-34)                            ! coin NW 
    551             call vel_V_NW 
     575            call vel_V_NW(Sv(i,j,:,:),opposy(i,j))  
    552576             
    553577         case(-41)                            ! coin SW 
    554578 
    555             call vel_V_SW 
     579            call vel_V_SW(Sv(i,j,:,:),opposy(i,j)) 
    556580             
    557581         end select case_imy 
    … …  
    559583   end do 
    560584end do lignes_UV 
     585!$OMP END DO 
     586!$OMP END PARALLEL 
    561587 
    562588end subroutine rempli_Tuij 
    … …  
    570596!------------------------------------------------------------------ 
    571597 
    572 subroutine vel_U_presc     !                     ! vitesse imposee 
    573  
    574   Tu(i,j,0,0)=1. 
    575   opposx(i,j)=uxprec(i,j) 
     598subroutine vel_U_presc(uxprec,Tu,opposx)     !                     ! vitesse imposee 
     599 
     600implicit none 
     601 
     602real,intent(in)  :: uxprec 
     603real,intent(out) :: Tu 
     604real,intent(out) :: opposx 
     605 
     606  Tu=1. 
     607  opposx=uxprec 
    576608 
    577609end subroutine vel_U_presc 
    578610!------------------------------------------------------------------ 
    579611                    
    580 subroutine vel_U_general                   ! cas general 
    581  
    582  
    583    beta=frotmx(i,j)*dx2                       !  Terme en u(i,j) 
    584  
    585  
    586  
    587    Tu(i,j,0,0)  = -4.*pvi(i,j)-4.*pvi(i-1,j)- pvm(i,j+1)-pvm(i,j)-beta 
    588  
    589    Tu(i,j,-1,0) =  4.*pvi(i-1,j)              !  Terme en u(i-1,j) 
    590  
    591    Tu(i,j,1,0)  =  4.*pvi(i,j)                !  Terme en u(i+1,j) 
    592  
    593    Tu(i,j,0,-1) =  pvm(i,j)                   !  Terme en u(i,j-1) 
    594  
    595    Tu(i,j,0,1)  =  pvm(i,j+1)                 !  Terme en u(i,j+1) 
    596  
    597    Tv(i,j,0,0)  = -2.*pvi(i,j)-pvm(i,j)       !  Terme en v(i,j) 
    598  
    599    Tv(i,j,-1,0) =  2.*pvi(i-1,j)+pvm(i,j)     !  Terme en v(i-1,j) 
    600  
    601    Tv(i,j,-1,1) = -2.*pvi(i-1,j)-pvm(i,j+1)   !  Terme en v(i-1,j+1) 
    602  
    603    Tv(i,j,0,1)  =  2.*pvi(i,j)+pvm(i,j+1)     !  Terme en v(i,j+1) 
    604  
    605    opposx(i,j)  =  rog*hmx_oppos(i,j)*sdx(i,j)*dx2  !  Terme en opposx(i,j) 
     612subroutine vel_U_general(frotmx,dx2,pvi,pvi_imoinsun,pvm,pvm_jplusun,sdx,rog,hmx_oppos,Tu,Tv,opposx)                   ! cas general 
     613 
     614implicit none 
     615 
     616real,intent(in) :: frotmx 
     617real,intent(in) :: dx2 
     618real,intent(in) :: pvi 
     619real,intent(in) :: pvi_imoinsun 
     620real,intent(in) :: pvm 
     621real,intent(in) :: pvm_jplusun 
     622real,intent(in) :: sdx 
     623real,intent(in) :: rog 
     624real,intent(in) :: hmx_oppos 
     625real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     626real,dimension(-1:1,-1:1),intent(inout) :: Tv 
     627real,intent(out) :: opposx 
     628 
     629   beta=frotmx*dx2                       !  Terme en u(i,j) 
     630 
     631   Tu(0,0)  = -4.*pvi - 4.*pvi_imoinsun - pvm_jplusun - pvm - beta 
     632 
     633   Tu(-1,0) =  4.*pvi_imoinsun              !  Terme en u(i-1,j) 
     634 
     635   Tu(1,0)  =  4.*pvi                !  Terme en u(i+1,j) 
     636 
     637   Tu(0,-1) =  pvm                   !  Terme en u(i,j-1) 
     638 
     639   Tu(0,1)  =  pvm_jplusun                 !  Terme en u(i,j+1) 
     640 
     641   Tv(0,0)  = -2.*pvi - pvm       !  Terme en v(i,j) 
     642 
     643   Tv(-1,0) =  2.*pvi_imoinsun + pvm     !  Terme en v(i-1,j) 
     644 
     645   Tv(-1,1) = -2.*pvi_imoinsun - pvm_jplusun   !  Terme en v(i-1,j+1) 
     646 
     647   Tv(0,1)  =  2.*pvi + pvm_jplusun     !  Terme en v(i,j+1) 
     648 
     649   opposx   =  rog*hmx_oppos*sdx*dx2  !  Terme en opposx(i,j) 
    606650 
    607651   return 
    … …  
    611655! voir explications dans vel_U_West 
    612656 
    613 subroutine vel_U_South                       ! bord sud non cisaillement            
    614  
    615   Tu(i,j,0,0)  =  1. 
    616   Tu(i,j,0,1)  = -1. 
    617   Tv(i,j,0,1)  = -1. 
    618   Tv(i,j,-1,1) =  1. 
    619   opposx(i,j)  =  0. 
     657subroutine vel_U_South(Tu,Tv,opposx)       ! bord sud non cisaillement            
     658 
     659implicit none                        
     660real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     661real,dimension(-1:1,-1:1),intent(inout) :: Tv 
     662real, intent(out) :: opposx 
     663 
     664  Tu(0,0)  =  1. 
     665  Tu(0,1)  = -1. 
     666  Tv(0,1)  = -1. 
     667  Tv(-1,1) =  1. 
     668  opposx   =  0. 
    620669 
    621670   return 
    … …  
    623672!------------------------------------------------------------------ 
    624673 
    625 subroutine vel_U_North                       ! bord nord   non cisaillement        
    626  
    627   Tu(i,j,0,0)  =  1. 
    628   Tu(i,j,0,-1) = -1. 
    629   Tv(i,j,0,0)  =  1. 
    630   Tv(i,j,-1,0) = -1. 
    631   opposx(i,j)  =  0. 
     674subroutine vel_U_North(Tu,Tv,opposx)     ! bord nord   non cisaillement        
     675 
     676implicit none                        
     677real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     678real,dimension(-1:1,-1:1),intent(inout) :: Tv 
     679real, intent(out) :: opposx 
     680 
     681  Tu(0,0)  =  1. 
     682  Tu(0,-1) = -1. 
     683  Tv(0,0)  =  1. 
     684  Tv(-1,0) = -1. 
     685  opposx   =  0. 
    632686 
    633687   return 
    … …  
    636690! voir explications dans vel_U_West 
    637691 
    638 subroutine vel_U_East                       ! bord Est   pression 
     692subroutine vel_U_East(pvi,pvi_imoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Tu,Tv,opposx) ! bord Est   pression 
     693 
     694implicit none 
     695real,intent(in) :: pvi 
     696real,intent(in) :: pvi_imoinsun 
     697real,intent(in) :: rog 
     698real,intent(in) :: H_imoinsun 
     699real,intent(in) :: rowg 
     700real,intent(in) :: slv_imoinsun 
     701real,intent(in) :: B_imoinsun 
     702real,intent(in) :: dx                        
     703real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     704real,dimension(-1:1,-1:1),intent(inout) :: Tv 
     705real, intent(out) :: opposx 
    639706        
    640   Tu(i,j,0,0)   =  4.*pvi(i-1,j) 
    641   Tu(i,j,-1,0)  = -4.*pvi(i-1,j) 
    642   Tv(i,j,0,1)   =  2.*pvi(i,j) 
    643   Tv(i,j,0,0)   = -2.*pvi(i,j) 
    644   opposx(i,j)   =  0.5*(rog*H(i-1,j)**2-rowg*(max(slv(i-1,j)-B(i-1,j),0.))**2)*dx 
     707  Tu(0,0)   =  4.*pvi_imoinsun 
     708  Tu(-1,0)  = -4.*pvi_imoinsun 
     709  Tv(0,1)   =  2.*pvi 
     710  Tv(0,0)   = -2.*pvi 
     711  opposx    =  0.5*(rog*H_imoinsun**2-rowg*(max(slv_imoinsun-B_imoinsun,0.))**2)*dx 
    645712 
    646713   return 
    … …  
    648715!------------------------------------------------------------------ 
    649716 
    650 subroutine vel_U_West                        ! bord West pression 
     717subroutine vel_U_West(pvi,pvi_imoinsun,rog,H,rowg,slv,B,dx,Tu,Tv,opposx) ! bord West pression 
    651718                                             ! tous les coef * -1 
     719implicit none 
     720real,intent(in) :: pvi 
     721real,intent(in) :: pvi_imoinsun 
     722real,intent(in) :: rog 
     723real,intent(in) :: H 
     724real,intent(in) :: rowg 
     725real,intent(in) :: slv 
     726real,intent(in) :: B 
     727real,intent(in) :: dx                        
     728real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     729real,dimension(-1:1,-1:1),intent(inout) :: Tv 
     730real, intent(out) :: opposx 
     731 
    652732 
    653733!  Tu(i,j,0,0)    =  4.*pvi(i-1,j) 
    654734!  Tu(i,j,1,0)    = -4.*pvi(i-1,j) 
    655735 
    656   Tu(i,j,0,0)    =  4.*pvi(i,j)              ! quand le bord est epais, il faut prendre  
    657   Tu(i,j,1,0)    = -4.*pvi(i,j)              ! le noeud du shelf pas celui en H=1 
     736  Tu(0,0)    =  4.*pvi              ! quand le bord est epais, il faut prendre  
     737  Tu(1,0)    = -4.*pvi              ! le noeud du shelf pas celui en H=1 
    658738                                             ! s'il est fin, pas de difference 
    659739 
    … …  
    661741!!$  Tv(i-1,j,-1,1) = -2.*pvi(i-1,j) 
    662742 
    663   Tv(i,j,-1,0) =  2.*pvi(i-1,j) 
    664   Tv(i,j,-1,1) = -2.*pvi(i-1,j) 
    665  opposx(i,j)    = -0.5*(rog*H(i,j)**2-rowg*(max(slv(i,j)-B(i,j),0.))**2)*dx 
     743  Tv(-1,0) =  2.*pvi_imoinsun 
     744  Tv(-1,1) = -2.*pvi_imoinsun 
     745 opposx    = -0.5*(rog*H**2-rowg*(max(slv-B,0.))**2)*dx 
    666746 
    667747  
    … …  
    674754!------------------------------------------------------------------ 
    675755 
    676 subroutine vel_V_presc                     ! vitesse imposee 
    677  
    678   Sv(i,j,0,0)=1. 
    679   opposy(i,j)=uyprec(i,j) 
     756subroutine vel_V_presc(uyprec,Sv,opposy)          ! vitesse imposee 
     757 
     758implicit none 
     759real,intent(in)  :: uyprec 
     760real,intent(out) :: Sv 
     761real,intent(out) :: opposy 
     762 
     763  Sv=1. 
     764  opposy=uyprec 
    680765 
    681766   return 
    … …  
    683768!------------------------------------------------------------------ 
    684769                     
    685 subroutine vel_V_general                      ! cas general 
    686  
    687  
    688    beta          =  frotmy(i,j)*dx2           !  Terme en v(i,j) 
    689    Sv(i,j,0,0)   = -4.*pvi(i,j)-4.*pvi(i,j-1)- pvm(i+1,j)-pvm(i,j)-beta 
    690  
    691    Sv(i,j,0,-1)  =  4.*pvi(i,j-1)             !  Terme en v(i,j-1) 
    692  
    693    Sv(i,j,0,1)   =  4.*pvi(i,j)               !  Terme en v(i,j+1) 
    694  
    695    Sv(i,j,-1,0)  =  pvm(i,j)                  !  Terme en v(i-1,j) 
    696  
    697    Sv(i,j,1,0)   =  pvm(i+1,j)                !  Terme en v(i+1,j) 
    698  
    699    Su(i,j,0,0)   = -2.*pvi(i,j)-pvm(i,j)      !  Terme en u(i,j) 
    700  
    701    Su(i,j,0,-1)  =  2.*pvi(i,j-1)+pvm(i,j)    !  Terme en u(i,j-1) 
    702  
    703    Su(i,j,1,-1)  = -2.*pvi(i,j-1)-pvm(i+1,j)  !  Terme en u(i+1,j-1) 
    704  
    705    Su(i,j,1,0)   =  2.*pvi(i,j)+pvm(i+1,j)    !  Terme en u(i+1,j) 
    706  
    707    opposy(i,j)   =  rog*hmy_oppos(i,j)*sdy(i,j)*dx2 !  Terme en opposy(i,j) 
     770subroutine vel_V_general(frotmy,dx2,pvi,pvi_jmoinsun,pvm,pvm_iplusun,sdy,rog,hmy_oppos,Su,Sv,opposy)! cas general 
     771 
     772implicit none 
     773real,intent(in) :: frotmy 
     774real,intent(in) :: dx2 
     775real,intent(in) :: pvi 
     776real,intent(in) :: pvi_jmoinsun 
     777real,intent(in) :: pvm 
     778real,intent(in) :: pvm_iplusun 
     779real,intent(in) :: sdy 
     780real,intent(in) :: rog 
     781real,intent(in) :: hmy_oppos 
     782real,dimension(-1:1,-1:1),intent(inout) :: Su 
     783real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     784real,intent(out) :: opposy 
     785 
     786real :: beta 
     787 
     788 
     789   beta          =  frotmy*dx2           !  Terme en v(i,j) 
     790   Sv(0,0)   = -4.*pvi - 4.*pvi_jmoinsun - pvm_iplusun - pvm-beta 
     791 
     792   Sv(0,-1)  =  4.*pvi_jmoinsun             !  Terme en v(i,j-1) 
     793 
     794   Sv(0,1)   =  4.*pvi               !  Terme en v(i,j+1) 
     795 
     796   Sv(-1,0)  =  pvm                  !  Terme en v(i-1,j) 
     797 
     798   Sv(1,0)   =  pvm_iplusun                !  Terme en v(i+1,j) 
     799 
     800   Su(0,0)   = -2.*pvi - pvm      !  Terme en u(i,j) 
     801 
     802   Su(0,-1)  =  2.*pvi_jmoinsun + pvm    !  Terme en u(i,j-1) 
     803 
     804   Su(1,-1)  = -2.*pvi_jmoinsun - pvm_iplusun  !  Terme en u(i+1,j-1) 
     805 
     806   Su(1,0)   =  2.*pvi + pvm_iplusun    !  Terme en u(i+1,j) 
     807 
     808   opposy    =  rog*hmy_oppos*sdy*dx2 !  Terme en opposy(i,j) 
    708809 
    709810   return 
    … …  
    712813!------------------------------------------------------------------ 
    713814 
    714 subroutine vel_V_West                       ! bord West non cisaillement            
    715  
    716   Sv(i,j,0,0)  =  1. 
    717   Sv(i,j,1,0)  = -1. 
    718   Su(i,j,1,0)  = -1. 
    719   Su(i,j,1,-1) =  1. 
    720   opposy(i,j)  =  0. 
     815subroutine vel_V_West(Sv,Su,opposy)     ! bord West non cisaillement            
     816 
     817implicit none 
     818real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     819real,dimension(-1:1,-1:1),intent(inout) :: Su 
     820real,intent(out) :: opposy 
     821 
     822  Sv(0,0)  =  1. 
     823  Sv(1,0)  = -1. 
     824  Su(1,0)  = -1. 
     825  Su(1,-1) =  1. 
     826  opposy   =  0. 
    721827 
    722828   return 
    … …  
    725831 
    726832 
    727 subroutine vel_V_East                       ! bord Est   non cisaillement        
    728  
    729   Sv(i,j,0,0)  =  1. 
    730   Sv(i,j,-1,0) = -1. 
    731   Su(i,j,0,0)  =  1. 
    732   Su(i,j,0,-1) = -1. 
    733   opposy(i,j)  =  0. 
     833subroutine vel_V_East(Sv,Su,opposy)         ! bord Est   non cisaillement        
     834 
     835implicit none 
     836real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     837real,dimension(-1:1,-1:1),intent(inout) :: Su 
     838real,intent(out) :: opposy 
     839 
     840  Sv(0,0)  =  1. 
     841  Sv(-1,0) = -1. 
     842  Su(0,0)  =  1. 
     843  Su(0,-1) = -1. 
     844  opposy   =  0. 
    734845 
    735846   return 
    … …  
    739850! voir explications dans vel_U_West 
    740851 
    741 subroutine vel_V_North                       ! bord Nord   pression 
     852subroutine vel_V_North(pvi,pvi_jmoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Sv,Su,opposy) ! bord Nord   pression 
     853 
     854implicit none 
     855real,intent(in) :: pvi 
     856real,intent(in) :: pvi_jmoinsun 
     857real,intent(in) :: rog 
     858real,intent(in) :: H_imoinsun 
     859real,intent(in) :: rowg 
     860real,intent(in) :: slv_imoinsun 
     861real,intent(in) :: B_imoinsun 
     862real,intent(in) :: dx 
     863real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     864real,dimension(-1:1,-1:1),intent(inout) :: Su 
     865real,intent(out) :: opposy 
    742866        
    743   Sv(i,j,0,0)  =  4.*pvi(i,j-1) 
    744   Sv(i,j,0,-1) = -4.*pvi(i,j-1) 
    745   Su(i,j,1,0)  =  2.*pvi(i,j) 
    746   Su(i,j,0,0)  = -2.*pvi(i,j) 
    747   opposy(i,j)  = 0.5*(rog*H(i-1,j)**2-rowg*(max(slv(i-1,j)-B(i-1,j),0.))**2)*dx 
     867  Sv(0,0)  =  4.*pvi_jmoinsun 
     868  Sv(0,-1) = -4.*pvi_jmoinsun 
     869  Su(1,0)  =  2.*pvi 
     870  Su(0,0)  = -2.*pvi 
     871  opposy   = 0.5*(rog*H_imoinsun**2-rowg*(max(slv_imoinsun-B_imoinsun,0.))**2)*dx 
    748872 
    749873   return 
    … …  
    752876! voir explications dans vel_U_West 
    753877 
    754 subroutine vel_V_South                        ! bord sud  pression 
     878subroutine vel_V_South(pvi,pvi_jmoinsun,rog,H,rowg,slv,B,dx,Sv,Su,opposy) ! bord sud  pression 
    755879                                              ! tous les coef * -1 
    756  
    757   Sv(i,j,0,0)  = 4.*pvi(i,j) 
    758   Sv(i,j,0,1)  =-4.*pvi(i,j) 
    759   Su(i,j,0,-1) = 2.*pvi(i,j-1) 
    760   Su(i,j,1,-1) =-2.*pvi(i,j-1) 
    761   opposy(i,j)  = -0.5*(rog*H(i,j)**2-rowg*(max(slv(i,j)-B(i,j),0.))**2)*dx 
     880implicit none 
     881real,intent(in) :: pvi 
     882real,intent(in) :: pvi_jmoinsun 
     883real,intent(in) :: rog 
     884real,intent(in) :: H 
     885real,intent(in) :: rowg 
     886real,intent(in) :: slv 
     887real,intent(in) :: B 
     888real,intent(in) :: dx 
     889real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     890real,dimension(-1:1,-1:1),intent(inout) :: Su 
     891real,intent(out) :: opposy 
     892 
     893  Sv(0,0)  = 4.*pvi 
     894  Sv(0,1)  =-4.*pvi 
     895  Su(0,-1) = 2.*pvi_jmoinsun 
     896  Su(1,-1) =-2.*pvi_jmoinsun 
     897  opposy   = -0.5*(rog*H**2-rowg*(max(slv-B,0.))**2)*dx 
    762898 
    763899   return 
    … …  
    780916!Coin SE 
    781917!---------- 
    782 subroutine vel_U_SE                         
    783  
    784   Tu(i,j,0,0) =  1.    !*pvi(i,j) 
    785   Tu(i,j,0,1) = -1.    !*pvi(i,j) 
    786   opposx(i,j) =  0. 
     918subroutine vel_U_SE(Tu,opposx)                         
     919 
     920implicit none 
     921real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     922real,intent(out) :: opposx 
     923 
     924  Tu(0,0) =  1.    !*pvi(i,j) 
     925  Tu(0,1) = -1.    !*pvi(i,j) 
     926  opposx  =  0. 
    787927 
    788928   return 
    … …  
    790930!------------------------------------------------------------------ 
    791931 
    792 subroutine vel_V_SE                         
    793  
    794   Sv(i,j,0,0) =  1.   !*pvi(i,j) 
    795   Sv(i,j,-1,0)= -1.   !*pvi(i,j) 
    796   opposy(i,j) =  0. 
     932subroutine vel_V_SE(Sv,opposy)                         
     933 
     934implicit none 
     935real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     936real,intent(out) :: opposy 
     937 
     938  Sv(0,0) =  1.   !*pvi(i,j) 
     939  Sv(-1,0)= -1.   !*pvi(i,j) 
     940  opposy  =  0. 
    797941 
    798942   return 
    … …  
    803947!Coin SW 
    804948!---------- 
    805 subroutine vel_U_SW                         
    806  
    807   Tu(i,j,0,0) =  1.   ! *pvi(i,j) 
    808   Tu(i,j,0,1) = -1.   ! *pvi(i,j) 
    809   opposx(i,j) =  0. 
     949subroutine vel_U_SW(Tu,opposx)                         
     950 
     951implicit none 
     952real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     953real,intent(out) :: opposx 
     954 
     955  Tu(0,0) =  1.   ! *pvi(i,j) 
     956  Tu(0,1) = -1.   ! *pvi(i,j) 
     957  opposx  =  0. 
    810958 
    811959   return 
    … …  
    813961!------------------------------------------------------------------ 
    814962 
    815 subroutine vel_V_SW                         
    816  
    817   Sv(i,j,0,0) =  1.   !*pvi(i,j) 
    818   Sv(i,j,1,0) = -1.   !*pvi(i,j) 
    819   opposy(i,j) =  0. 
     963subroutine vel_V_SW(Sv,opposy) 
     964 
     965implicit none 
     966real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     967real,intent(out) :: opposy                     
     968 
     969  Sv(0,0) =  1.   !*pvi(i,j) 
     970  Sv(1,0) = -1.   !*pvi(i,j) 
     971  opposy  =  0. 
    820972 
    821973   return 
    … …  
    825977!Coin NE 
    826978!---------- 
    827 subroutine vel_U_NE                         
    828  
    829   Tu(i,j,0,0)  =  1.  ! *pvi(i,j) 
    830   Tu(i,j,0,-1) = -1.  ! *pvi(i,j) 
    831   opposx(i,j)  =  0. 
     979subroutine vel_U_NE(Tu,opposx) 
     980 
     981implicit none 
     982real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     983real,intent(out) :: opposx                      
     984 
     985  Tu(0,0)  =  1.  ! *pvi(i,j) 
     986  Tu(0,-1) = -1.  ! *pvi(i,j) 
     987  opposx   =  0. 
    832988 
    833989   return 
    … …  
    835991!------------------------------------------------------------------ 
    836992 
    837 subroutine vel_V_NE                         
    838  
    839   Sv(i,j,0,0)  =  1.  ! *pvi(i,j) 
    840   Sv(i,j,-1,0) = -1.  ! *pvi(i,j) 
    841   opposy(i,j)  =  0. 
     993subroutine vel_V_NE(Sv,opposy)                         
     994 
     995implicit none 
     996real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     997real,intent(out) :: opposy 
     998 
     999  Sv(0,0)  =  1.  ! *pvi(i,j) 
     1000  Sv(-1,0) = -1.  ! *pvi(i,j) 
     1001  opposy   =  0. 
    8421002 
    8431003   return 
    … …  
    8471007!Coin NW 
    8481008!---------- 
    849 subroutine vel_U_NW                         
    850  
    851   Tu(i,j,0,0)  =  1.  ! *pvi(i,j) 
    852   Tu(i,j,0,-1) = -1.  ! *pvi(i,j) 
    853   opposx(i,j)  =  0. 
     1009subroutine vel_U_NW(Tu,opposx) 
     1010 
     1011implicit none                         
     1012real,dimension(-1:1,-1:1),intent(inout) :: Tu 
     1013real,intent(out) :: opposx  
     1014 
     1015  Tu(0,0)  =  1.  ! *pvi(i,j) 
     1016  Tu(0,-1) = -1.  ! *pvi(i,j) 
     1017  opposx   =  0. 
    8541018 
    8551019   return 
    … …  
    8571021!------------------------------------------------------------------ 
    8581022 
    859 subroutine vel_V_NW                         
    860  
    861   Sv(i,j,0,0)  =  1.  !*pvi(i,j) 
    862   Sv(i,j,1,0)  = -1.  !*pvi(i,j) 
    863   opposy(i,j)  =  0. 
     1023subroutine vel_V_NW(Sv,opposy)                         
     1024 
     1025implicit none                         
     1026real,dimension(-1:1,-1:1),intent(inout) :: Sv 
     1027real,intent(out) :: opposy  
     1028 
     1029  Sv(0,0)  =  1.  !*pvi(i,j) 
     1030  Sv(1,0)  = -1.  !*pvi(i,j) 
     1031  opposy   =  0. 
    8641032 
    8651033   return 

  • trunk/SOURCES/litho-0.4.f90

    r65 r74  
    3737 
    3838  
    39       subroutine litho 
     39subroutine litho 
     40  !$ USE OMP_LIB 
     41  USE module3D_phy 
     42  USE param_phy_mod 
     43  USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie 
    4044 
    41       USE module3D_phy 
    42       USE param_phy_mod 
    43       USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie 
     45  implicit none 
    4446 
    45       implicit none 
    46  
    47       INTEGER :: II,SOM1,SOM2 
    48       REAL, dimension(:,:), allocatable :: WLOC 
     47  INTEGER :: II,SOM1,SOM2 
     48  REAL, dimension(:,:), allocatable :: WLOC 
    4949 
    5050 
    51 !----- allocation de WLOC  et de croix ----------- 
     51  !----- allocation de WLOC  et de croix ----------- 
    5252 
    53       if (.not.allocated(WLOC)) THEN 
    54          allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err) 
    55          if (err/=0) then 
    56             print *,"Erreur a l'allocation du tableau WLOC",err 
    57             stop 4 
    58          end if 
    59       end if 
     53  if (.not.allocated(WLOC)) THEN 
     54     allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err) 
     55     if (err/=0) then 
     56        print *,"Erreur a l'allocation du tableau WLOC",err 
     57        stop 4 
     58     end if 
     59  end if 
    6060 
    61 !----- fin de l'allocation -------------- 
     61  !----- fin de l'allocation -------------- 
    6262 
    63 !     calcul de la deflexion par sommation des voisins appartenant 
    64 !     au bloc de taille 2*LBLOC 
    65        som1=0. 
    66        som2=0. 
     63  !     calcul de la deflexion par sommation des voisins appartenant 
     64  !     au bloc de taille 2*LBLOC 
     65  som1=0. 
     66  som2=0. 
    6767 
    68 !      On somme aussi les contributions des points exterieurs au domaine 
    69 !      lorsque la charge est due a l'ocean. On suppose alors  que 
    70 !      ces points ont la meme charge que les limites 
    71  
    72 boucleij : do J=1,NY 
    73            do I=1,NX 
     68  !      On somme aussi les contributions des points exterieurs au domaine 
     69  !      lorsque la charge est due a l'ocean. On suppose alors  que 
     70  !      ces points ont la meme charge que les limites 
     71  !$OMP PARALLEL 
     72  !$OMP DO PRIVATE(Wloc) 
     73  boucleij : do J=1,NY 
     74     do I=1,NX 
    7475 
    7576 
    76               W1(I,J)=0.  
    77               ii=0 
     77        W1(I,J)=0.  
     78 !   ii=0 
    7879 
    79 ! Wloc : enfoncement centre sur ij du a chaque charges autour 
    80       Wloc(:,:) = We(:,:) * charge(i-lbloc:i+lbloc,j-lbloc:j+lbloc) 
     80 ! Wloc : enfoncement centre sur ij du a chaque charges autour 
     81        Wloc(:,:) = We(:,:) * charge(i-lbloc:i+lbloc,j-lbloc:j+lbloc) 
    8182 
    8283 
    83       W1(i,j) = sum (Wloc(:,:))   ! effet en ij de toutes les charges autour (distance < lbloc) 
     84        W1(i,j) = sum (Wloc(:,:))   ! effet en ij de toutes les charges autour (distance < lbloc) 
    8485 
    85 !    ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC 
    86               som1=som1+W1(I,J) 
    87               som2=som2-charge(I,J)/ROMG 
     86 !    ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC 
     87 !   som1=som1+W1(I,J) 
     88 !   som2=som2-charge(I,J)/ROMG 
    8889 
    89       end do 
    90       end do  boucleij 
    91         
    92 !     write(6,*)'enfoncement total avec rigidite de la lithosphere:', 
    93 !    &  som1 
    94 !     write(6,*)'enfoncement total avec isostasie locale :', 
    95 !    &  som2 
     90     end do 
     91  end do  boucleij 
     92  !$OMP END DO 
     93  !$OMP END PARALLEL 
     94 
     95  !     write(6,*)'enfoncement total avec rigidite de la lithosphere:', 
     96  !    &  som1 
     97  !     write(6,*)'enfoncement total avec isostasie locale :', 
     98  !    &  som2 
    9699 
    97100 
    98     end subroutine litho 
    99        
     101end subroutine litho 
     102 

  • trunk/SOURCES/taubed-0.3.f90

    r4 r74  
    4343subroutine taubed() 
    4444 
     45  !$USE OMP_LIB  
    4546  USE module3D_phy 
    4647  USE param_phy_mod 
    … …  
    5253  ! NLITH est defini dans isostasie et permet le choix du modele d'isostasie 
    5354 
     55   
    5456  if (NLITH.eq.1) then 
    5557     !       avec rigidite de la lithosphere 
     58     !$OMP PARALLEL 
     59     !$OMP DO 
    5660     do J=1,NY  
    5761        do I=1,NX 
    … …  
    6569        end do 
    6670     end do 
    67  
     71     !$OMP END DO 
    6872 
    6973     ! il faut remplir CHARGE dans les parties a l'exterieur de la grille : 
    7074     ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord 
    71  
     75     !$OMP DO 
    7276     do J=1,NY 
    7377        CHARGE(1-LBLOC:0,J)=CHARGE(1,J)      ! bord W 
    7478        CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E 
    7579     end do 
     80     !$OMP END DO 
     81     !$OMP DO 
    7682     do I=1,NX 
    7783        CHARGE(I,1-LBLOC:0)=CHARGE(I,1)      ! bord S 
    7884        CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N 
    7985     end do 
     86     !$OMP END DO 
    8087 
    8188     ! valeur dans les quatres coins exterieurs au domaine        
     89          !$OMP WORKSHARE 
    8290     CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1)           ! coin SW 
    8391     CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY)      ! coin NW 
    8492     CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1)      ! coin SE 
    8593     CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE 
    86  
     94     !$OMP END WORKSHARE 
     95     !$OMP END PARALLEL  
    8796     call litho 
    8897 
    8998  else 
    9099     !     enfoncement local 
     100     !$OMP PARALLEL 
     101     !$OMP DO 
    91102     do J=1,NY 
    92103        do I=1,NX 
    … …  
    100111        end do 
    101112     end do 
     113     !$OMP END DO 
     114     !$OMP END PARALLEL 
    102115  endif 
    103116 
    104117  !     decroissance exponentielle de l'enfoncement 
     118  !$OMP PARALLEL 
     119  !$OMP DO 
    105120  do J=1,NY 
    106121     do I=1,NX 
    … …  
    109124     end do 
    110125  end do 
     126  !$OMP END DO 
     127  !$OMP END PARALLEL 
    111128 
    112129end subroutine taubed 
Note: See TracChangeset for help on using the changeset viewer.