Changeset 76 for trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90

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:

• trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90 (modified) (17 diffs)

trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90

@@ -80 +80 @@
 subroutine icethick3
 
+!$ USE OMP_LIB
+
 implicit none
 real,dimension(nx,ny) :: Dminx,Dminy
@@ -92 +94 @@
 ! Copie de H sur vieuxH        
 ! --------------------
+!$OMP PARALLEL
+!$OMP WORKSHARE
 where (mk0(:,:).eq.0)
    vieuxH(:,:)=0.
@@ -99 +103 @@
    vieuxH(:,:)=H(:,:)
 end where
+!$OMP END WORKSHARE
 
 ! mk0 est le masque à ne jamais dépasser
@@ -110 +115 @@
 jmaxdia=0
 
-! attention limitation ! 
+! attention limitation !
+!$OMP WORKSHARE 
 uxbar(:,:)=max(uxbar(:,:),-V_limit)
 uxbar(:,:)=min(uxbar(:,:),V_limit)
 uybar(:,:)=max(uybar(:,:),-V_limit)
 uybar(:,:)=min(uybar(:,:),V_limit)
+!$OMP END WORKSHARE
 
 ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante
-do i=2,nx-1
-   do j=2,ny-1
+!$OMP END PARALLEL
+!$OMP PARALLEL
+!$OMP DO PRIVATE(aux)
+!do i=2,nx-1
+!   do j=2,ny-1
+do j=2,ny-1
+  do i=2,nx-1
       aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) &
            + (abs(uybar(i,j)) + abs(uybar(i,j+1)))
       aux=aux*dx1*0.5
+      !$OMP CRITICAL
       if (aux.gt.maxdia) then
          maxdia=aux
-         imaxdia=i
-         jmaxdia=j
+!         imaxdia=i
+!         jmaxdia=j
       endif
+      !$OMP END CRITICAL
    end do
 end do
-
+!$OMP END DO
+!$OMP END PARALLEL
 
 if (maxdia.ge.(testdiag/dtmax)) then
@@ -141 +156 @@
 ! avec le pas de temps long (dtt)
 
- call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
+call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
 
 
 ! calcul diffusivite - vitesse 
 !----------------------------------------------------------------- 
-
+!$OMP PARALLEL
+!$OMP WORKSHARE
 Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif
 Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar
@@ -172 +188 @@
   uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot
 end where
-
+!$OMP END WORKSHARE
 
 
@@ -179 +195 @@
 
 if (adv_frac.gt.1.) then ! advection seulement
+   !$OMP WORKSHARE
    advmx(:,:)=Uxbar(:,:)
    advmy(:,:)=Uybar(:,:)
    Dmx(:,:)=0.
    Dmy(:,:)=0.
-
+   !$OMP END WORKSHARE
 else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement
+   !$OMP WORKSHARE
    advmx(:,:)=0.
    advmy(:,:)=0.
-
+   !$OMP END WORKSHARE
    ! D reste la valeur au dessus)
 
@@ -197 +215 @@
 ! diffusion = ce qui peut être exprime en diffusion 
 !             + une partie (1-adv_frac) de l'advection
-
+   !$OMP WORKSHARE
    where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general
        advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion
@@ -236 +254 @@
           dmy(:,:)=dminy(:,:)*H(:,:)
        end where
-
-
-
    elsewhere   ! zones trop plates ou trop fines
       Dmy(:,:)=0.
       advmy(:,:)=Uybar(:,:)
    end where
-
+   !$OMP END WORKSHARE
 
 else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs
-
+   !$OMP WORKSHARE
    zonemx(:,:)=flgzmx(:,:)
    zonemy(:,:)=flgzmy(:,:)
@@ -265 +280 @@
     advmy(:,:)=0.
    end where
-
+   !$OMP END WORKSHARE
 end if
 
-
+!$OMP WORKSHARE
 where (flot(:,:) )
    tabtest(:,:)=1.
@@ -274 +289 @@
    tabtest(:,:)=0.
 end where
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 !------------------------------------------------------------------detect
 !!$ tabtest(:,:)=dmx(:,:)
@@ -375 +391 @@
 
 
-debug_3D(:,:,45)=dmx(:,:)
-debug_3D(:,:,46)=dmy(:,:)
-debug_3D(:,:,47)=advmx(:,:)
-debug_3D(:,:,48)=advmy(:,:)
-
-
+!debug_3D(:,:,45)=dmx(:,:)
+!debug_3D(:,:,46)=dmy(:,:)
+!debug_3D(:,:,47)=advmx(:,:)
+!debug_3D(:,:,48)=advmy(:,:)
+
+!$OMP PARALLEL
+!$OMP WORKSHARE
 bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 ! diverses precription de l'epaisseur en fonction de l'objectif
@@ -387 +406 @@
 
 
-  call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions)
-
-
-   
+call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions)
 
 
@@ -436 +452 @@
 !!$end if
 
-debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)
-
-debug_3D(:,:,86)=Hp(:,:)
-debug_3D(:,:,85)=i_Hp(:,:)
+!debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)
+
+!debug_3D(:,:,86)=Hp(:,:)
+!debug_3D(:,:,85)=i_Hp(:,:)
 
 !!$where (i_hp(:,:).eq.1)
@@ -463 +479 @@
 call resol_adv_diff_2D_vect (Dmx,Dmy,advmx,advmy,Hp,i_Hp,bilmass,vieuxH,H)
 
-! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord 
+! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord
+!$OMP PARALLEL
+!$OMP WORKSHARE 
 where (H(:,:).lt.0.)
    ablbord(:,:)=H(:,:)/dt
@@ -489 +507 @@
 endwhere
 
-
 ! eventuellement retirer apres spinup ou avoir un cas serac
 Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt
+!$OMP END WORKSHARE
 
 if (igrdline.ne.3) then
+   !$OMP WORKSHARE
    Hdot(:,:)=min(Hdot(:,:),10.)
    Hdot(:,:)=max(Hdot(:,:),-10.)
+   !$OMP END WORKSHARE
 endif
 
+!$OMP WORKSHARE
 where (i_hp(:,:).ne.1)
    H(:,:)=vieuxH(:,:)+Hdot(:,:)*dt
@@ -511 +532 @@
    Hdot(:,:)=0.
 end where
-
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 !calul de l'ablation sur les bords (pourrait n'être appelé que pour les sorties courtes)