Changeset 76


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

Location:
trunk/SOURCES
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc

    r75 r76  
    4040        IFORT= ifort 
    4141 
    42         ARITHMi = -O2 -fp-model precise 
     42        ARITHMi = -O2 -fp-model precise  #-openmp 
    4343                                                                      # (normalement reproductible) 
    4444    ifeq ($(debug), 1) 

  • trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90

    r24 r76  
    1515 
    1616Subroutine Thermal_prop_icetemp 
     17!$ USE OMP_LIB 
    1718  Use Icetemp_declar  
    1819 
    … …  
    2122  If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Thermal_prop_icetemp' 
    2223 
     24!$OMP PARALLEL 
     25!$OMP DO COLLAPSE(2) 
    2326  Do K=1,Nz   
    2427     Do J=1,Ny 
    … …  
    2932           !  Attention Pour La Conductivite C'Est La Formule De Yin-Chao Yen 
    3033 
    31            Ct(I,J,K)=6.62E7 
    32            Cp(I,J,K)=2009.  
     34!           Ct(I,J,K)=6.62E7 
     35!           Cp(I,J,K)=2009.  
    3336           Tpmp(I,J,K)=-0.00087*(K-1)*1./(Nz-1)*H(I,J)!De=1/(Nz_m-1) 
    3437           Cp(I,J,K)=(2115.3+7.79293*T(I,J,K))  
    … …  
    3841     End Do 
    3942  End Do 
    40  
     43!$OMP END DO 
     44!$OMP END PARALLEL 
    4145End Subroutine Thermal_prop_icetemp 
    4246 

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

    r65 r76  
    8080subroutine icethick3 
    8181 
     82!$ USE OMP_LIB 
     83 
    8284implicit none 
    8385real,dimension(nx,ny) :: Dminx,Dminy 
    … …  
    9294! Copie de H sur vieuxH         
    9395! -------------------- 
     96!$OMP PARALLEL 
     97!$OMP WORKSHARE 
    9498where (mk0(:,:).eq.0) 
    9599   vieuxH(:,:)=0. 
    … …  
    99103   vieuxH(:,:)=H(:,:) 
    100104end where 
     105!$OMP END WORKSHARE 
    101106 
    102107! mk0 est le masque à ne jamais dépasser 
    … …  
    110115jmaxdia=0 
    111116 
    112 ! attention limitation !  
     117! attention limitation ! 
     118!$OMP WORKSHARE  
    113119uxbar(:,:)=max(uxbar(:,:),-V_limit) 
    114120uxbar(:,:)=min(uxbar(:,:),V_limit) 
    115121uybar(:,:)=max(uybar(:,:),-V_limit) 
    116122uybar(:,:)=min(uybar(:,:),V_limit) 
     123!$OMP END WORKSHARE 
    117124 
    118125! le pas de temps est choisi pour rendre la matrice advective diagonale dominante 
    119 do i=2,nx-1 
    120    do j=2,ny-1 
     126!$OMP END PARALLEL 
     127!$OMP PARALLEL 
     128!$OMP DO PRIVATE(aux) 
     129!do i=2,nx-1 
     130!   do j=2,ny-1 
     131do j=2,ny-1 
     132  do i=2,nx-1 
    121133      aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) & 
    122134           + (abs(uybar(i,j)) + abs(uybar(i,j+1))) 
    123135      aux=aux*dx1*0.5 
     136      !$OMP CRITICAL 
    124137      if (aux.gt.maxdia) then 
    125138         maxdia=aux 
    126          imaxdia=i 
    127          jmaxdia=j 
     139!         imaxdia=i 
     140!         jmaxdia=j 
    128141      endif 
     142      !$OMP END CRITICAL 
    129143   end do 
    130144end do 
    131  
     145!$OMP END DO 
     146!$OMP END PARALLEL 
    132147 
    133148if (maxdia.ge.(testdiag/dtmax)) then 
    … …  
    141156! avec le pas de temps long (dtt) 
    142157 
    143  call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
     158call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
    144159 
    145160 
    146161! calcul diffusivite - vitesse  
    147162!-----------------------------------------------------------------  
    148  
     163!$OMP PARALLEL 
     164!$OMP WORKSHARE 
    149165Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif 
    150166Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar 
    … …  
    172188  uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot 
    173189end where 
    174  
     190!$OMP END WORKSHARE 
    175191 
    176192 
    … …  
    179195 
    180196if (adv_frac.gt.1.) then ! advection seulement 
     197   !$OMP WORKSHARE 
    181198   advmx(:,:)=Uxbar(:,:) 
    182199   advmy(:,:)=Uybar(:,:) 
    183200   Dmx(:,:)=0. 
    184201   Dmy(:,:)=0. 
    185  
     202   !$OMP END WORKSHARE 
    186203else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement 
     204   !$OMP WORKSHARE 
    187205   advmx(:,:)=0. 
    188206   advmy(:,:)=0. 
    189  
     207   !$OMP END WORKSHARE 
    190208   ! D reste la valeur au dessus) 
    191209 
    … …  
    197215! diffusion = ce qui peut être exprime en diffusion  
    198216!             + une partie (1-adv_frac) de l'advection 
    199  
     217   !$OMP WORKSHARE 
    200218   where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general 
    201219       advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion 
    … …  
    236254          dmy(:,:)=dminy(:,:)*H(:,:) 
    237255       end where 
    238  
    239  
    240  
    241256   elsewhere   ! zones trop plates ou trop fines 
    242257      Dmy(:,:)=0. 
    243258      advmy(:,:)=Uybar(:,:) 
    244259   end where 
    245  
     260   !$OMP END WORKSHARE 
    246261 
    247262else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs 
    248  
     263   !$OMP WORKSHARE 
    249264   zonemx(:,:)=flgzmx(:,:) 
    250265   zonemy(:,:)=flgzmy(:,:) 
    … …  
    265280    advmy(:,:)=0. 
    266281   end where 
    267  
     282   !$OMP END WORKSHARE 
    268283end if 
    269284 
    270  
     285!$OMP WORKSHARE 
    271286where (flot(:,:) ) 
    272287   tabtest(:,:)=1. 
    … …  
    274289   tabtest(:,:)=0. 
    275290end where 
    276  
     291!$OMP END WORKSHARE 
     292!$OMP END PARALLEL 
    277293!------------------------------------------------------------------detect 
    278294!!$ tabtest(:,:)=dmx(:,:) 
    … …  
    375391 
    376392 
    377 debug_3D(:,:,45)=dmx(:,:) 
    378 debug_3D(:,:,46)=dmy(:,:) 
    379 debug_3D(:,:,47)=advmx(:,:) 
    380 debug_3D(:,:,48)=advmy(:,:) 
    381  
    382  
     393!debug_3D(:,:,45)=dmx(:,:) 
     394!debug_3D(:,:,46)=dmy(:,:) 
     395!debug_3D(:,:,47)=advmx(:,:) 
     396!debug_3D(:,:,48)=advmy(:,:) 
     397 
     398!$OMP PARALLEL 
     399!$OMP WORKSHARE 
    383400bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance 
     401!$OMP END WORKSHARE 
     402!$OMP END PARALLEL 
    384403 
    385404! diverses precription de l'epaisseur en fonction de l'objectif 
    … …  
    387406 
    388407 
    389   call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions) 
    390  
    391  
    392     
     408call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions) 
    393409 
    394410 
    … …  
    436452!!$end if 
    437453 
    438 debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88) 
    439  
    440 debug_3D(:,:,86)=Hp(:,:) 
    441 debug_3D(:,:,85)=i_Hp(:,:) 
     454!debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88) 
     455 
     456!debug_3D(:,:,86)=Hp(:,:) 
     457!debug_3D(:,:,85)=i_Hp(:,:) 
    442458 
    443459!!$where (i_hp(:,:).eq.1) 
    … …  
    463479call resol_adv_diff_2D_vect (Dmx,Dmy,advmx,advmy,Hp,i_Hp,bilmass,vieuxH,H) 
    464480 
    465 ! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord  
     481! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord 
     482!$OMP PARALLEL 
     483!$OMP WORKSHARE  
    466484where (H(:,:).lt.0.) 
    467485   ablbord(:,:)=H(:,:)/dt 
    … …  
    489507endwhere 
    490508 
    491  
    492509! eventuellement retirer apres spinup ou avoir un cas serac 
    493510Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt 
     511!$OMP END WORKSHARE 
    494512 
    495513if (igrdline.ne.3) then 
     514   !$OMP WORKSHARE 
    496515   Hdot(:,:)=min(Hdot(:,:),10.) 
    497516   Hdot(:,:)=max(Hdot(:,:),-10.) 
     517   !$OMP END WORKSHARE 
    498518endif 
    499519 
     520!$OMP WORKSHARE 
    500521where (i_hp(:,:).ne.1) 
    501522   H(:,:)=vieuxH(:,:)+Hdot(:,:)*dt 
    … …  
    511532   Hdot(:,:)=0. 
    512533end where 
    513  
    514  
     534!$OMP END WORKSHARE 
     535!$OMP END PARALLEL 
    515536 
    516537!calul de l'ablation sur les bords (pourrait n'être appelé que pour les sorties courtes) 

  • trunk/SOURCES/diffusiv-polyn-0.6.f90

    r68 r76  
    3434 
    3535!     =============================================================== 
    36  
     36!$ USE OMP_LIB 
    3737USE module3D_phy 
    3838 
    … …  
    5858 
    5959! initialisation de la diffusion 
     60!$OMP PARALLEL 
     61!$OMP WORKSHARE 
    6062diffmx(:,:)=0. 
    6163diffmy(:,:)=0. 
    62  
     64!$OMP END WORKSHARE 
     65!$OMP END PARALLEL 
    6366 
    6467! calcul de SDX, SDY et de la pente au carre en mx et en my : 
    … …  
    7982! Calcul de la vitesse de glissement, qu'elle soit calculee par diagno ou par sliding 
    8083!------------------------------------------------------------------------------------ 
    81 !  
     84!$OMP PARALLEL 
     85!$OMP WORKSHARE 
    8286where (flgzmx(:,:)) 
    8387    ubx(:,:) = uxflgz(:,:) 
    … …  
    9195    uby(:,:) = ddby(:,:)* (-sdy(:,:)) 
    9296endwhere 
    93  
    94  
     97!$OMP END WORKSHARE 
     98!$OMP END PARALLEL 
    9599if (itracebug.eq.1) call tracebug(' diffusiv  apres calcul glissement') 
    96100 
    … …  
    99103! Deformation SIA : Calcul de ddy  et ddx pour chaque valeur de iglen 
    100104!------------------------------------------------------------------ 
    101  
    102105do iglen=n1poly,n2poly 
    103106   glenexp=max(0.,(glen(iglen)-1.)/2.) 
    104  
    105  
     107  !$OMP PARALLEL  
     108  !$OMP DO 
    106109   do j=1,ny 
    107110      do i=1,nx 
    108  
    109111         if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams 
    110112            ddy(i,j,iglen)=((slope2my(i,j)**  &  ! slope2my calcule en debut de routine 
    111113                 glenexp)*(rog)**glen(iglen)) & 
    112114                 *hmy(i,j)**(glen(iglen)+1) 
    113  
    114  
    115115         endif 
    116116      end do 
    117117   end do 
     118   !$OMP END DO 
     119   !$OMP END PARALLEL 
    118120end do 
     121 
     122 
    119123 
    120124do iglen=n1poly,n2poly 
    121125   glenexp=max(0.,(glen(iglen)-1.)/2.) 
    122  
    123  
     126   !$OMP PARALLEL 
     127   !$OMP DO 
    124128   do j=1,ny   
    125129      do i=1,nx 
    126  
    127130         if (.not.flotmx(i,j)) then      ! bug y->x corrige le 5 avril 2008 ( 
    128                                      
    129131            ddx(i,j,iglen)=((slope2mx(i,j)**  & ! slope2mx calcule en debut de routine 
    130132                 glenexp)*(rog)**glen(iglen)) & 
    … …  
    133135      end do 
    134136   end do 
     137   !$OMP END DO 
     138   !$OMP END PARALLEL 
    135139end do 
    136140 
     141 
    137142! vitesse due a la déformation----------------------------------------------- 
    138  
    139  
     143!$OMP PARALLEL 
     144!$OMP DO 
    140145ydef: do j=2,ny 
    141146       do i=1,nx  
    … …  
    164169    end do 
    165170   end do ydef 
    166  
     171!$OMP END DO 
     172 
     173!$OMP DO 
    167174xdef: do j=1,ny 
    168175       do i=2,nx  
    … …  
    192199    end do 
    193200   end do xdef 
     201!$OMP END DO 
    194202 
    195203if (itracebug.eq.1) call tracebug(' diffusiv  avant limit') 
    … …  
    199207 
    200208! la limitation selon x 
     209!$OMP WORKSHARE 
    201210where(.not.flgzmx(:,:)) 
    202211   uxbar(:,:)=min(uxbar(:,:),ultot) 
    … …  
    209218   uybar(:,:)=max(uybar(:,:),-ultot) 
    210219end where 
     220!$OMP END WORKSHARE 
    211221 
    212222!    cas particulier des sommets ou il ne reste plus de neige. 
     223!$OMP DO 
    213224do j=2,ny-1 
    214225   do i=2,nx-1 
    … …  
    221232   end do 
    222233end do 
    223  
     234!$OMP END DO 
     235!$OMP END PARALLEL 
    224236if (itracebug.eq.1)  call tracebug(' fin diffusiv ') 
    225237 

  • trunk/SOURCES/dragging_neff_slope_mod.f90

    r66 r76  
    168168!------------------------------------------------------------------------- 
    169169subroutine dragging   ! defini la localisation des streams et le frottement basal 
    170  
     170!$ USE OMP_LIB 
    171171 
    172172!         les iles n'ont pas de condition neff mais ont la condition toblim 
    173173!         (ce bloc etait avant dans flottab) 
    174174 
    175  
     175!$OMP PARALLEL 
     176!$OMP DO 
    176177do j=2,ny 
    177178   do i=2,nx 
    … …  
    180181   end do 
    181182end do 
     183!$OMP END DO 
    182184 
    183185! le calcul des fleuves se fait sur les mailles majeures en partant de la cote 
    … …  
    193195! Attention : ce systeme ne permet pas d'avoir des fleuves convergents 
    194196!             et les fleuves n'ont une épaisseur que de 1 noeud. 
    195  
     197!$OMP WORKSHARE 
    196198fleuvemx(:,:)=.false. 
    197199fleuvemy(:,:)=.false. 
    198200fleuve(:,:)=.false. 
    199201cote(:,:)=.false. 
     202!$OMP END WORKSHARE 
    200203 
    201204! detection des cotes 
     205!$OMP DO 
    202206do  j=2,ny-1   
    203207   do i=2,nx-1   
    … …  
    208212   end do 
    209213end do 
     214!$OMP END DO 
    210215 
    211216! calcul de neff (pression effective sur noeuds majeurs) 
     217!$OMP DO 
    212218do j=1,ny-1 
    213219   do i=1,nx-1 
    … …  
    215221   enddo 
    216222enddo 
     223!$OMP END DO 
    217224!aurel, for the borders: 
     225!$OMP WORKSHARE 
    218226neff(:,ny)=1.e5 
    219227neff(nx,:)=1.e5 
     228 
    220229 
    221230!!$ 
    … …  
    277286! aurel, we add the neff threshold: 
    278287where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true. 
    279  
     288!$OMP END WORKSHARE 
     289 
     290!$OMP DO 
    280291do j=1,ny-1 
    281292   do i=1,nx-1 
    … …  
    288299   end do 
    289300end do 
     301!$OMP END DO 
    290302 
    291303! we look for the warm based points that will not be treated as stream (ub from SSA): 
     304!$OMP WORKSHARE 
    292305slowssa(:,:)=.false. 
    293306slowssamx(:,:)=.false. 
    294307slowssamy(:,:)=.false. 
     308!$OMP END WORKSHARE 
     309!$OMP DO 
    295310do  j=1,ny 
    296311   do i=1,nx 
    … …  
    299314   end do 
    300315end do 
     316!$OMP END DO 
     317!$OMP DO 
    301318do j=1,ny-1 
    302319   do i=1,nx-1 
    … …  
    309326   end do 
    310327end do 
    311  
     328!$OMP END DO 
     329!$OMP DO 
    312330do j=1,ny 
    313331   do i=1,nx 
    … …  
    354372   end do 
    355373end do 
    356  
     374!$OMP END DO 
     375 
     376!$OMP DO 
    357377do j=1,ny 
    358378   do i=1,nx 
    … …  
    399419   end do 
    400420end do 
    401  
    402  
     421!$OMP END DO 
     422 
     423!$OMP DO 
    403424do j=2,ny-1 
    404425   do i=2,nx-1 
    … …  
    407428   end do 
    408429end do 
    409  
    410  
    411 where (fleuve(:,:)) 
    412    debug_3D(:,:,1)=1 
    413 elsewhere (flot(:,:)) 
    414    debug_3D(:,:,1)=2 
    415 elsewhere 
    416    debug_3D(:,:,1)=0 
    417 endwhere 
    418  
    419 where (cote(:,:)) 
    420    debug_3D(:,:,2)=1 
    421 elsewhere 
    422    debug_3D(:,:,2)=0 
    423 endwhere 
    424  
    425 where (fleuvemx(:,:)) 
    426    debug_3D(:,:,3)=1 
    427 elsewhere 
    428    debug_3D(:,:,3)=0 
    429 endwhere 
    430  
    431 where (flotmx(:,:)) 
    432    debug_3D(:,:,3)=-1 
    433 endwhere 
     430!$OMP END DO 
     431!$OMP END PARALLEL 
     432 
     433!~ where (fleuve(:,:)) 
     434!~    debug_3D(:,:,1)=1 
     435!~ elsewhere (flot(:,:)) 
     436!~    debug_3D(:,:,1)=2 
     437!~ elsewhere 
     438!~    debug_3D(:,:,1)=0 
     439!~ endwhere 
     440!~  
     441!~ where (cote(:,:)) 
     442!~    debug_3D(:,:,2)=1 
     443!~ elsewhere 
     444!~    debug_3D(:,:,2)=0 
     445!~ endwhere 
     446!~  
     447!~ where (fleuvemx(:,:)) 
     448!~    debug_3D(:,:,3)=1 
     449!~ elsewhere 
     450!~    debug_3D(:,:,3)=0 
     451!~ endwhere 
     452!~  
     453!~ where (flotmx(:,:)) 
     454!~    debug_3D(:,:,3)=-1 
     455!~ endwhere 
    434456 
    435457!_____________________ 
    436458 
    437459 
    438 where (fleuvemy(:,:)) 
    439    debug_3D(:,:,4)=1 
    440 elsewhere 
    441    debug_3D(:,:,4)=0 
    442 endwhere 
    443  
    444 where (flotmy(:,:)) 
    445    debug_3D(:,:,4)=-1 
    446 end where 
    447  
    448 debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5)  
    449 debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5)  
     460!~ where (fleuvemy(:,:)) 
     461!~    debug_3D(:,:,4)=1 
     462!~ elsewhere 
     463!~    debug_3D(:,:,4)=0 
     464!~ endwhere 
     465!~  
     466!~ where (flotmy(:,:)) 
     467!~    debug_3D(:,:,4)=-1 
     468!~ end where 
     469!~  
     470!~ debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5)  
     471!~ debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5)  
    450472 
    451473return 

  • 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 

  • trunk/SOURCES/relaxation_water_diffusion.f90

    r65 r76  
    1818 
    1919subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER) 
    20     implicit none 
     20 
     21  !$ USE OMP_LIB 
     22 
     23  implicit none 
    2124 
    2225 
    … …  
    5760 
    5861! write(166,*)' entree  relaxation waterdif' 
     62!$OMP PARALLEL 
     63!$OMP WORKSHARE 
    5964    HWATER(:,:)= vieuxHWATER(:,:) 
    60  
     65!$OMP END WORKSHARE 
    6166    !   calcul de kmx et kmx a partir de KOND 
    6267    !   conductivite hyrdraulique sur les noeuds mineurs 
    6368    !   moyenne harmonique 
    6469    !   ---------------------------------------- 
    65  
     70!$OMP DO 
    6671    do j=2,nyy 
    6772       do i=2,nxx 
    … …  
    7580       end do 
    7681    end do 
    77  
     82!$OMP END DO 
     83 
     84!$OMP DO 
    7885    do j=2,nyy 
    7986       do i=2,nxx 
    … …  
    8693       enddo 
    8794    enddo 
    88  
     95!$OMP END DO 
    8996 
    9097    !   attribution des coefficients  arelax .... 
    … …  
    100107    dtwdx2=dt/dx/dx 
    101108 
    102 arelax(:,:)=0. 
    103 brelax(:,:)=0. 
    104 crelax(:,:)=1. 
    105 drelax(:,:)=0. 
    106 erelax(:,:)=0. 
    107 frelax(:,:)=limit_hw(:,:) 
    108  
    109  
    110  
     109!$OMP WORKSHARE 
     110  arelax(:,:)=0. 
     111  brelax(:,:)=0. 
     112  crelax(:,:)=1. 
     113  drelax(:,:)=0. 
     114  erelax(:,:)=0. 
     115  frelax(:,:)=limit_hw(:,:) 
     116!$OMP END WORKSHARE 
     117 
     118!$OMP DO 
    111119    do J=2,NYY-1 
    112120       do I=2,NXX-1 
    113121 
    114122          if (klimit(i,j).eq.0) then 
    115  
    116123! calcul du vecteur 
    117  
    118124          FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT 
    119125          frelax(i,j)=frelax(i,j)+(kmx(i,j)*pgmx(i,j)-kmx(i+1,j)*pgmx(i+1,j))*dtsrgdx 
    120126          frelax(i,j)=frelax(i,j)+(kmy(i,j)*pgmy(i,j)-kmy(i,j+1)*pgmy(i,j+1))*dtsrgdx 
    121  
    122127! calcul des diagonales       
    123128          arelax(i,j)=-kmx(i,j)*dtwdx2        ! arelax : diagonale i-1,j  
    … …  
    131136          crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2  
    132137                                              !crelax : diagonale i,j  
    133  
    134138          else if (klimit(i,j).eq.1) then 
    135139             hwater(i,j)=limit_hw(i,j) 
    136140!             write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j) 
    137  
    138141          endif 
    139  
    140142       end do 
    141143    end do 
    142  
    143  
     144!$OMP END DO 
     145!$OMP END PARALLEL 
    144146 
    145147    ! Boucle de relaxation : 
    … …  
    152154       ntour=ntour+1 
    153155!       write(6,*) 'boucle de relaxation numero',ntour 
    154  
     156       !$OMP PARALLEL 
     157       !$OMP DO PRIVATE(reste) 
    155158       do j=2,NYY-1 
    156159          do i=2,NXX-1 
    … …  
    161164 
    162165             DELTAH(I,J) = RESTE/CRELAX(I,J)              
    163  
    164166          end do 
    165167       end do 
    166  
    167 deltah(:,:)=min(deltah(:,:),10.) 
    168 deltah(:,:)=max(deltah(:,:),-10.) 
    169       
    170  
    171  
     168       !$OMP END DO 
     169 
     170       !$OMP WORKSHARE 
     171       deltah(:,:)=min(deltah(:,:),10.) 
     172       deltah(:,:)=max(deltah(:,:),-10.) 
     173       !$OMP END WORKSHARE 
    172174       ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction 
    173175       ! de hwater sur les points voisins. 
     176       !$OMP DO 
    174177       do j=2,NYY-1 
    175178          do i=2,NXX-1 
    176179             HWATER(I,J) = HWATER(I,J) - DELTAH(I,J) 
    177  
    178180          end do 
    179181       end do 
    180  
     182       !$OMP END DO 
    181183 
    182184       ! critere d'arret: 
    183  
    184        
    185  
    186185       Delh=0 
    187186       Vh=0 
    188  
     187        
     188       !$OMP DO REDUCTION(+:Delh) 
    189189       DO j=2,NYY-1 
    190190          DO i=2,NXX-1 
    191191              
    192192             !   write(166,*) I,J,delh,deltah(i,j) 
    193                          Delh=Delh+deltah(i,j)**2 
     193             Delh=Delh+deltah(i,j)**2 
    194194             !            Vh=Vh+h(i,j)**2. 
    195195          END DO 
    196196       END DO 
    197  
     197       !$OMP END DO 
     198       !$OMP END PARALLEL 
    198199!       write(6,*) delh,maxval(deltah),minval(deltah) 
    199200       !      testh=SQRT(Delh/Vh) 
Note: See TracChangeset for help on using the changeset viewer.