Ignore:
Timestamp:
06/29/16 16:21:13 (8 years ago)
Author:
dumas
Message:

OpenMP parallelization in conserv-mass-adv-diff_sept2009_mod.f90, diffusiv-polyn-0.6.f90, dragging_neff_slope_mod.f90, eaubasale-0.5_mod.f90 and relaxation_water_diffusion.f90

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/eaubasale-0.5_mod.f90

    r4 r76  
    133133subroutine eaubasale !(pwater)   version correspondant à la thèse de Vincent 
    134134  ! A terme, il faudrait en faire un module  
    135  
    136  
     135!$ USE OMP_LIB 
     136!~   integer :: t1,t2,ir 
     137!~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme   
     138!~    
     139!~    ! Temps CPU de calcul initial. 
     140!~    call cpu_time(t_cpu_0) 
     141!~    ! Temps elapsed de reference. 
     142!~    call system_clock(count=t1, count_rate=ir) 
    137143 
    138144 
    … …  
    141147  if (itracebug.eq.1)  call tracebug(' Entree dans routine eau basale') 
    142148 
     149!$OMP PARALLEL 
     150!$OMP WORKSHARE 
    143151  vieuxhwater(:,:) = hwater(:,:) 
    144152  kond(:,:) = kond0*SECYEAR  
    … …  
    148156  klimit(:,:)=0 
    149157  limit_hw(:,:)=-9999. 
     158!$OMP END WORKSHARE 
     159!$OMP DO   
    150160  do j=1,ny 
    151161     do i=1,nx 
    … …  
    167177     end do 
    168178  end do 
    169  
    170  
    171  
    172   do i=1,nx 
    173      do j=2,ny-1 
     179!$OMP END DO 
     180 
     181!$OMP DO 
     182!  do j=2,ny-1 
     183  do j=1,ny 
     184     do i=1,nx 
    174185        bmelt_w(i,j)=bmelt(i,j)*rofresh/ro     
    175186        hwater(i,j)=max(hwater(i,j),0.) 
    176187        hw(i,j)=min(hwater(i,j),hmax_wat)  
    177  
    178188     enddo 
    179189  enddo 
    180  
     190!$OMP END DO 
     191!$OMP END PARALLEL 
    181192 
    182193  ecoul:  if (ecoulement_eau) then 
    183194     !  print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5' 
    184195     !   write(6,*) 'entree ecoulement_eau' 
    185  
    186      do J=1,NY 
    187         do I=1,NX 
     196     !$OMP PARALLEL 
     197     !$OMP DO 
     198     do j=1,ny 
     199        do i=1,nx 
    188200           tetar(i,j)=(xlong(i,j))*PI/180. ! pourrait etre fait une fois pour toute 
    189201           PGX(I,J)=101*sin(tetar(i,j))*1.e-2                  
    … …  
    200212        enddo 
    201213     enddo 
     214!$OMP END DO 
    202215 
    203216     ! sorties debug 17 juillet 2007 
    204      debug_3D(:,:,5)=pot_w(:,:) 
    205      debug_3D(:,:,6)=pot_f(:,:) 
    206      debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:) 
    207      debug_3D(:,:,8)=hwater(:,:) 
    208  
     217!     debug_3D(:,:,5)=pot_w(:,:) 
     218!     debug_3D(:,:,6)=pot_f(:,:) 
     219!     debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:) 
     220!     debug_3D(:,:,8)=hwater(:,:) 
     221!$OMP DO 
    209222     do j=2,ny 
    210223        do i=2,nx 
    211224           if (H(i,j).gt.25.) then 
    212  
    213225              !           calcul du gradient de pression 
    214               ! _______________________________________ 
    215  
    216226              if (flotmx(i,j)) then  
    217227                 pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx 
    … …  
    225235                 pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy 
    226236              endif 
    227  
    228  
    229237           endif 
    230  
    231  
    232238           pgx(i,j)=pgx(i,j)/rofreshg   ! pour passer pgx sans unité 
    233239           pgy(i,j)=pgy(i,j)/rofreshg       
    234240        end do 
    235241     end do 
    236  
    237  
     242!$OMP END DO 
     243!$OMP END PARALLEL 
    238244 
    239245     if (nt.gt.0) then 
    240246        if (dt.gt.0.)  then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0 
    241  
    242  
    243247           !------------------------------------------------------------------------- 
    244248           ! les points de la grounding line ont une conductivité hydraulique élevée 
    245249           ! même si la base est froide.  modif catritz 18 janvier 2005 
    246250           !open(166,file='erreur_eau') 
    247  
    248            do J=2,NY-1 
    249               do I=2,NX-1 
    250  
     251!$OMP PARALLEL 
     252!$OMP DO 
     253           do j=2,Ny-1 
     254              do i=2,Nx-1 
    251255                 ! cond=0 si glace froide et pas sur la grounding line 
    252  
    253                  if ((IBASE(I,J).eq.1).and.   & 
     256                 if ((IBASE(i,j).eq.1).and.   & 
    254257                      (.not.flot(i,j-1)).and.(.not.flot(i,j+1)).and.  & 
    255                       (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(I,J)=0.! 1.0e-5 
     258                      (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(i,j)=0.! 1.0e-5 
    256259 
    257260                 ! cond infinie quand epaisseur faible et glace flottante 
    … …  
    267270 
    268271                 ! conductivite effective (conductivité * taille du tuyau en m2/an)  
    269  
    270272                 keff(i,j)=kond(i,j)*hw(i,j) 
    271273              end do 
    272274           end do 
    273  
     275!$OMP END DO 
    274276           ! condition sur les bords de la grille 
    275  
     277!$OMP WORKSHARE 
    276278           kond(1,:)=kondmax 
    277279           kond(nx,:)=kondmax 
    278280           kond(:,1)=kondmax 
    279281           kond(:,ny)=kondmax 
    280  
    281            ! fin de la modif 18 janvier 2005--------------------------------------- 
    282  
    283282           vieuxhwater(:,:) = hwater(:,:) 
    284  
    285  
     283!$OMP END WORKSHARE 
     284!$OMP END PARALLEL 
     285!!$OMP SINGLE 
    286286           call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,keffmax,hwater) 
    287  
    288  
    289  
     287!!$OMP END SINGLE 
    290288        else 
    291289           !         print*,'dt=',dt,'pas de relax_water' 
    … …  
    303301 
    304302     ! on le fait avant les butoirs pour justement pouvoir les estimer 
    305  
     303!$OMP PARALLEL 
    306304     if (dt.gt.0.) then 
    307305        !         print*,'dt=',dt,'pas de relax_water' 
     306        !$OMP DO 
    308307        do j=1,ny 
    309308           do i=1,nx 
    … …  
    311310           end do 
    312311        end do 
     312        !$OMP END DO 
    313313     endif 
    314314 
    315  
     315!$OMP DO PRIVATE(Zflot) 
    316316     do i=1,nx 
    317317        do j=1,ny 
    … …  
    353353              !          Hwater(i,j)=hw(i,j) 
    354354              !          pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5) 
    355  
    356  
    357  
    358  
    359355           endif 
    360  
    361356 
    362357           ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression 
    … …  
    365360           !             pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till) 
    366361           !            endif 
    367  
    368  
    369  
    370  
    371362        end do 
    372363     end do 
    373  
     364!$OMP END DO 
    374365 
    375366     !   ************ valeurs de HWATER pour les coins ******** 
    … …  
    381372 
    382373     ! pour les sorties de flux d'eau 
     374     !$OMP DO 
    383375     do j=2,ny 
    384376        do i=2,nx 
    … …  
    390382              phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j)) 
    391383           endif 
    392  
    393384           ! ligne pour sortir les pgx 
    394  
    395385           pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) 
    396386        end do 
    397387     end do 
    398  
     388     !$OMP END DO 
     389     !$OMP DO 
    399390     do j=2,ny 
    400391        do i=2,nx 
    … …  
    406397           endif 
    407398           pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) 
    408  
    409399        enddo 
    410400     enddo 
    411  
    412  
     401     !$OMP END DO 
     402     !$OMP END PARALLEL 
    413403 
    414404  else  ! partie originellement dans icetemp à mettre dans un autre module 
    … …  
    416406 
    417407     if (ISYNCHRO.eq.1) then 
     408        !$OMP PARALLEL 
     409        !$OMP DO 
    418410        do j=1,ny 
    419411           do i=1,nx 
    … …  
    429421           end do 
    430422        end do 
     423        !$OMP END DO 
     424        !$OMP END PARALLEL 
    431425     end if 
     426      
    432427  endif ecoul 
     428 
     429!~   ! Temps elapsed final 
     430!~   call system_clock(count=t2, count_rate=ir) 
     431!~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
     432!~   ! Temps CPU de calcul final 
     433!~   call cpu_time(t_cpu_1) 
     434!~   t_cpu = t_cpu_1 - t_cpu_0 
     435 
     436!~   ! Impression du resultat. 
     437!~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
     438!~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
     439!~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
     440!~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
     441!~            nx,ny,temps,t_cpu,norme 
    433442 
    434443  return 
Note: See TracChangeset for help on using the changeset viewer.