Changeset 71 for trunk/SOURCES/flottab2-0.7.f90

Timestamp:

06/15/16 17:13:33 (8 years ago)

Author:

dumas

Message:

First parallelization instructions with OpenMP tested in debug : exactly the same results | isostasy called only every 50 years

File:

• trunk/SOURCES/flottab2-0.7.f90 (modified) (51 diffs)

trunk/SOURCES/flottab2-0.7.f90

@@ -16 +16 @@
 module flottab_mod
   
+  !$ USE OMP_LIB
 USE module3D_phy
 use module_choix
@@ -104 +105 @@
 
 
+!~   integer :: t1,t2,ir
+!~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme  
+  
+!~    ! Temps CPU de calcul initial.
+!~    call cpu_time(t_cpu_0)
+!~    ! Temps elapsed de reference.
+!~    call system_clock(count=t1, count_rate=ir)
+
 
 if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab')
@@ -112 +121 @@
 ! cas particulier des runs paleo ou on impose un masque grounded
 
+!$OMP PARALLEL PRIVATE(archim,surnet)
  if (igrdline.eq.2) then
+    !$OMP WORKSHARE
     where ( mk_init(:,:).eq.1)                 ! pose
        flot(:,:) = .False.
@@ -120 +131 @@
        flot(:,:) = .True.
     end where
+    !$OMP END WORKSHARE
  end if
 
@@ -129 +141 @@
 
        appel_new_flot=.false.
+       !$OMP DO
        do j=1,ny
           do i=1,nx
@@ -136 +149 @@
           enddo
        enddo
-
+       !$OMP END DO
+       
 ! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m
 
+      !$OMP WORKSHARE
       ICE(:,:)=0
       front(:,:)=0
@@ -148 +163 @@
       cotemy(:,:)=.false.
       boost=.false.
+      !$OMP END WORKSHARE
 
 ! fin de l'initialisation
@@ -155 +171 @@
 ! -------------------------------------
 
-
+ !$OMP DO
  do j=1,ny
    do i=1,nx
@@ -223 +239 @@
     end do
  end do
-
+ !$OMP END DO
 
 !!$ do i=1,nx
@@ -240 +256 @@
 
 !-----------------------------------------------------------------------
-
+ !$OMP DO
 domain_x: do j=1,ny       
          do i=2,nx
@@ -280 +296 @@
           end do
        end do domain_x
-
-if (itracebug.eq.1)  call tracebug('  routine flottab apres domain_x')
+ !$OMP END DO
+!if (itracebug.eq.1)  call tracebug('  routine flottab apres domain_x')
 
 !     3_y    B- NOUVELLE DEFINITION DE FLOTMY
 !         --------------------------------
+ !$OMP DO
 domain_y: do j=2,ny       
          do i=1,nx
@@ -318 +335 @@
           end do
        end do domain_y
-
+ !$OMP END DO
 
 
@@ -351 +368 @@
 !      4- determination des iles
 !      -------------------------
- 
+          !$OMP WORKSHARE
           ilemx(:,:)=.false.
           ilemy(:,:)=.false.
-
-!       selon x 
+          !$OMP END WORKSHARE
+
+!       selon x
+  !$OMP DO
 ilesx:  do j=2,ny-1
            do i=3,nx-2
@@ -404 +423 @@
          end do
       end do ilesx
-
-!       selon y 
+   !$OMP END DO
+
+!       selon y
+   !$OMP DO
 ilesy: do j=3,ny-2
            do i=2,nx-1
@@ -452 +473 @@
          end do
       end do ilesy
+    !$OMP END DO
+    !$OMP END PARALLEL
 ! fin des iles
 
@@ -479 +502 @@
 
 !   6- calcule les vitesses des points qui sont devenus gzm
-         
+!$OMP PARALLEL
+!$OMP DO         
 do j=1,ny
    do i=2,nx-1
@@ -492 +516 @@
    end do
 end do
-
+!$OMP END DO
+
+!$OMP DO 
 do j=2,ny-1
    do i=1,nx
@@ -504 +530 @@
     end do
 end do
-
+!$OMP END DO
 
 
@@ -511 +537 @@
 
 if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage
+   !$OMP WORKSHARE
    uxbar(:,:)=uxs1(:,:)
-   uybar(:,:)=uys1(:,:) 
+   uybar(:,:)=uys1(:,:)
+   !$OMP END WORKSHARE
 endif
 
+!$OMP WORKSHARE
 flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:)))   &
      .or.(.not.marine.and.flotmx(:,:))
 flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:)))   &
      .or.(.not.marine.and.flotmy(:,:))
-
+!$OMP END WORKSHARE
 
 
@@ -526 +555 @@
 !       fbm est vrai si le point est flottant mais un des voisins est pose 
 !_________________________________________________________________________
+!$OMP DO
 do j=2,ny-1
   do i=2,nx-1
@@ -536 +566 @@
   end do
 end do
-
+!$OMP END DO
 
 
@@ -552 +582 @@
 !!$end do
 
+!$OMP WORKSHARE
 where (flot(:,:))
    where (H(:,:).gt.(1.1)) 
@@ -565 +596 @@
    end where
 end where
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 call DETERMIN_TACHE 
@@ -581 +614 @@
 !----------------------------------------------!
 !On determine les differents ice strean/shelf  !
-      call DETERMIN_TACHE                      !
+!      call DETERMIN_TACHE                      !
 !----------------------------------------------!
 
@@ -597 +630 @@
 
 !On compte comme englacé uniquement les calottes dont une partie est posée      
-
+!$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j)
+!$OMP DO
 do i=3,nx-2
    do j=3,ny-2
@@ -666 +700 @@
    end do
 end do
-
+!$OMP END DO
+!$OMP END PARALLEL
 
 !----------------------------------------------
@@ -691 +726 @@
 !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j)
 
-
+!~   ! Temps elapsed final
+!~   call system_clock(count=t2, count_rate=ir)
+!~   temps=real(t2 - t1,kind=4)/real(ir,kind=4)
+!~   ! Temps CPU de calcul final
+!~   call cpu_time(t_cpu_1)
+!~   t_cpu = t_cpu_1 - t_cpu_0
+ 
+!~   ! Impression du resultat.
+!~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              &
+!~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, &
+!~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, &
+!~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', &
+!~            nx,ny,temps,t_cpu,norme
 
 end subroutine flottab
@@ -700 +747 @@
 !>
 subroutine determin_tache
+
+!!$ USE OMP_LIB
 
 implicit none
@@ -720 +769 @@
 enddo
 !      table_in  = .false.
-
+!!$OMP PARALLEL
+!!$OMP WORKSHARE
 table_out(:,:) = 0
 iceberg(:)  = .true.
 icetrim (:) = .true.
 nb_pts_tache(:) = 0
-
+!!$OMP END WORKSHARE
+!!$OMP END PARALLEL
 !    open(unit=100,file="tache.data",status='replace')
 
 ! 2-reperage des taches
 !----------------------
+!!$OMP PARALLEL PRIVATE(mask,label,indice)
+!!$OMP DO
 do i=2,nx-1
   do j=2,ny-1
@@ -736 +789 @@
 
      IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
-
+                        
         if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then   !masque de 2 cases adjacentes
            !      un des voisins est deja en glace
@@ -747 +800 @@
               if (mask(indice).gt.0) label=min(label,mask(indice))
            enddo
+!cdc       label=min(label,minval(mask(:), mask=mask > 0))
 
            !on fixe la valeur de la tache voisine minimun au point etudie (via label)
@@ -798 +852 @@
   enddo
 enddo
-
+!!$OMP END DO
+!!$OMP END PARALLEL
 
 
 ! On reorganise compt en ecrivant le numero de la tache fondamentale
 ! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
-! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)        
+! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)
+
 do indice=1,label_max
    vartemp = compt(indice)
@@ -813 +869 @@
 enddo
 
+!!$OMP PARALLEL
+!!$OMP DO REDUCTION(+:nb_pts_tache)
 do i=1,nx
    do j=1,ny
@@ -821 +879 @@
    enddo
 enddo
-
-
+!!$OMP END DO
+!!$OMP END PARALLEL
 
 
@@ -853 +911 @@
 !>
 subroutine determin_front
-
+!!$ USE OMP_LIB
 integer :: i_moins1,i_plus1,i_plus2
 integer :: j_moins1,j_plus1,j_plus2
       
-      do i=3,nx-2
-       do j=3,ny-2
+      !!$OMP PARALLEL
+      !!$OMP DO
+      do j=3,ny-2
+        do i=3,nx-2
 
  surice:if  (ice(i,j).eq.0) then
@@ -910 +970 @@
        end do
       end do
-
+      !!$OMP END DO
 
 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
@@ -923 +983 @@
 
 !     print*,'dans remplissage baies',time
+       
 baies: do k=1,2
+!!$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2)
 do j=1,ny
    do i=1,nx
@@ -958 +1020 @@
    end do
 end do
+!!$OMP END DO
 end do baies
 
@@ -969 +1032 @@
 !!$end if
 
-
-do i=2,nx-1
-   do j=2,ny-1
+!!$OMP DO
+do j=2,ny-1
+   do i=2,nx-1
 
       if (ice(i,j).eq.1) then     !         test si ice=1
@@ -983 +1046 @@
    end do
 end do
+!!$OMP END DO
 
 ! traitement des bords. On considere que l'exterieur n'a pas de glace
 ! attention ce n'est vrai que sur la grande grille
 
-
+!!$OMP DO PRIVATE(i)
 do j=2,ny-1
    i=1
@@ -994 +1058 @@
    front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
 end do
-
+!!$OMP END DO 
+
+!!$OMP DO PRIVATE(j)
 do i=2,nx-1
    j=1 
@@ -1001 +1067 @@
    front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
 end do
-
+!!$OMP END DO 
 ! traitement des coins
 
@@ -1023 +1089 @@
 !   si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites
 
-
+!!$OMP DO 
 do j=1,ny
    do i=1,nx-1
@@ -1036 +1102 @@
    end do
 end do
-
+!!$OMP END DO
+
+!!$OMP DO 
 do j=1,ny-1
    do i=1,nx
@@ -1049 +1117 @@
    end do
 end do
+!!$OMP END DO
 
 !isolx signifie pas de voisins en x
@@ -1056 +1125 @@
 
 ! calcul de frontfacex et isolx
+!!$OMP DO
 do j=1,ny
    do i=2,nx-1
@@ -1075 +1145 @@
    end do
 end do
+!!$OMP END DO
 
 ! calcul de frontfacey et isoly
+!!$OMP DO
 do j=2,ny-1
    do i=1,nx
@@ -1096 +1168 @@
    end do
 end do
-
+!!$OMP END DO
 
 
@@ -1102 +1174 @@
 ! attention ce n'est vrai que sur la grande grille
 
-
+!!$OMP DO PRIVATE(i)
 do j=2,ny-1
    i=1
@@ -1121 +1193 @@
    end if
 end do
-
+!!$OMP END DO
+
+!!$OMP DO PRIVATE(j)
 do i=2,nx-1
    j=1 
@@ -1140 +1214 @@
    end if
 end do
-
-
+!!$OMP END DO
+!!OMP END PARALLEL
 
 return