Changeset 427 for branches/GRISLIv3


Ignore:
Timestamp:
04/25/23 15:00:21 (15 months ago)
Author:
dumas
Message:

Use only in reso_adv_diff_2D_vect relaxation_waterdif_mod and dat2netcdf

Location:
branches/GRISLIv3/SOURCES
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/readinput.f90

    r6 r427  
    1010 
    1111  use netcdf 
    12   use io_netcdf_grisli 
     12  use io_netcdf_grisli, only:write_ncdf_var 
    1313  implicit none 
    1414  integer,intent(in)  :: nxx                    !< dimension along x 
    … …  
    1616  integer,intent(in)  :: numcol                 !< column number  
    1717  character(len=*), intent(in) :: base_name     !< base name of the variable to write 
     18  character(len=*),intent(in) :: filename       !< name of the file 
     19  character(len=*),intent(in) :: fil_sortie     !< output filename 
     20 
     21  ! local variables 
    1822  character(len=20) :: base_name1               !< base name of the variable to write 
    19   character(len=*),intent(in) :: filename       !< name of the file 
    20   character(len=*)::fil_sortie 
    2123  real*8,dimension(:,:),pointer :: Tab1,Tab2,Tab3    !< the array that is read and returned 
    2224  integer  :: i,j                !< working integers 
    … …  
    109111subroutine init_ncdf(nxx,nyy,fil_sortie,dimnames2d) 
    110112  use netcdf 
    111   use io_netcdf_grisli 
     113  use io_netcdf_grisli, only: write_ncdf_dim 
    112114  implicit none  
    113115  integer,intent(in)  :: nxx                    !< dimension along x 
    114116  integer,intent(in)  :: nyy                    !< dimension along y 
    115   character(len=*)::fil_sortie 
    116   character(len=20),dimension(2) :: dimnames2d  !< dimensions pour netcdf pour tableau 2d 
     117  character(len=*),intent(in)::fil_sortie 
     118  character(len=20),dimension(2),intent(out) :: dimnames2d  !< dimensions pour netcdf pour tableau 2d 
    117119  ! initialisation 
    118120  call write_ncdf_dim('x',trim(fil_sortie),nxx)               ! dimensions des tableaux 
    … …  
    132134subroutine lect_ncfile(varname,Tab,filename) 
    133135  use netcdf 
    134   use io_netcdf_grisli 
     136  use io_netcdf_grisli, only: read_ncdf_var 
    135137  implicit none  
    136138 
    … …  
    140142  real*8, dimension(:,:), pointer :: tabvar 
    141143 
    142  
    143144  if(.not. associated(tabvar)) then 
    144145     allocate(tabvar(size(Tab,1),size(Tab,2)))  
    145146  else  ! to nullify the pointer if undefined  
    146     if (any(shape(tabvar).ne.(/size(Tab,1),size(Tab,2)/)) ) then  
     147     if (any(shape(tabvar).ne.(/size(Tab,1),size(Tab,2)/)) ) then  
    147148        tabvar => null() 
    148     end if 
     149     end if 
    149150  end if 
    150  
    151151 
    152152  call read_ncdf_var(varname,filename,tabvar) 
    153153  Tab(:,:)= tabvar(:,:) 
    154 !write(6,*) filename,varname,Tab(154,123) 
     154  !write(6,*) filename,varname,Tab(154,123) 
    155155end subroutine lect_ncfile 
    156156 
    … …  
    168168 
    169169subroutine lect_input (numcol,basename,col,tabvar,filename,ncfileout)  
    170   !use netcdf 
    171   !use io_netcdf_grisli 
     170 
    172171  implicit none 
    173172  integer,intent(in)  :: numcol                 !< column number 
    174173  character(len=*), intent(in) :: basename      !< base name of the variable to read 
    175174  integer,intent(in)  :: col                    !< column   1->average, 2->minval 3->maxval   
    176                                                 !          ??? double emploi avec numcol 
     175  !          ??? double emploi avec numcol 
    177176  real,dimension(:,:),intent(inout) :: tabvar   !< array to read in the variable 
    178177  character(len=*),intent(in) :: filename       !< name of the file to read 
    … …  
    197196     call lect_ncfile(trim(basename)//'_'//xcol,tabvar,ncfileout)                
    198197  else 
    199   
     198 
    200199     if(index(filename,'.grd') .ne.0) then   ! file in is a .grd 
    201   
     200 
    202201        call lect_ncfile('z',tabvar,filename) 
    203202     else 
    204                             ! file in is a .nc 
     203        ! file in is a .nc 
    205204        if(index(filename,'.nc') .ne.0) then 
    206205           call lect_ncfile(basename,tabvar,filename) 

  • branches/GRISLIv3/SOURCES/relaxation_water_diffusion.f90

    r196 r427  
    1717 
    1818 
    19 subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,HWATER) 
    20  
    21   !$ USE OMP_LIB 
    22  
    23   implicit none 
     19  subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,HWATER) 
     20 
     21    !$ USE OMP_LIB 
     22 
     23    implicit none 
    2424 
    2525 
    … …  
    5656    REAL,dimension(NXX,NYY) :: KMX, KMY 
    5757    !    REAL  :: RHO,RHOW,RHOG  !,SECYEAR!! ice density, water density, density*acceleration 
    58 !    write(6,*) 'entree relaxation' 
    59  
    60 ! write(166,*)' entree  relaxation waterdif' 
    61 !$OMP PARALLEL 
    62 !$OMP WORKSHARE 
     58    !    write(6,*) 'entree relaxation' 
     59 
     60    ! write(166,*)' entree  relaxation waterdif' 
     61    !$OMP PARALLEL 
     62    !$OMP WORKSHARE 
    6363    HWATER(:,:)= vieuxHWATER(:,:) 
    64 !$OMP END WORKSHARE 
     64    !$OMP END WORKSHARE 
    6565    !   calcul de kmx et kmx a partir de KOND 
    6666    !   conductivite hyrdraulique sur les noeuds mineurs 
    6767    !   moyenne harmonique 
    6868    !   ---------------------------------------- 
    69 !$OMP DO 
     69    !$OMP DO 
    7070    do j=2,nyy 
    7171       do i=2,nxx 
    … …  
    7979       end do 
    8080    end do 
    81 !$OMP END DO 
    82  
    83 !$OMP DO 
     81    !$OMP END DO 
     82 
     83    !$OMP DO 
    8484    do j=2,nyy 
    8585       do i=2,nxx 
    … …  
    9292       enddo 
    9393    enddo 
    94 !$OMP END DO 
     94    !$OMP END DO 
    9595 
    9696    !   attribution des coefficients  arelax .... 
    … …  
    106106    dtwdx2=dt/dx/dx 
    107107 
    108 !$OMP WORKSHARE 
    109   deltah(:,:)=0. 
    110   arelax(:,:)=0. 
    111   brelax(:,:)=0. 
    112   crelax(:,:)=1. 
    113   drelax(:,:)=0. 
    114   erelax(:,:)=0. 
    115   frelax(:,:)=limit_hw(:,:) 
    116 !$OMP END WORKSHARE 
    117   reste      =0. 
    118  
    119 !$OMP DO 
     108    !$OMP WORKSHARE 
     109    deltah(:,:)=0. 
     110    arelax(:,:)=0. 
     111    brelax(:,:)=0. 
     112    crelax(:,:)=1. 
     113    drelax(:,:)=0. 
     114    erelax(:,:)=0. 
     115    frelax(:,:)=limit_hw(:,:) 
     116    !$OMP END WORKSHARE 
     117    reste      =0. 
     118 
     119    !$OMP DO 
    120120    do J=2,NYY-1 
    121121       do I=2,NXX-1 
    122122 
    123123          if (klimit(i,j).eq.0) then 
    124 ! calcul du vecteur 
    125           FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT 
    126           frelax(i,j)=frelax(i,j)+(kmx(i,j)*pgmx(i,j)-kmx(i+1,j)*pgmx(i+1,j))*dtsrgdx 
    127           frelax(i,j)=frelax(i,j)+(kmy(i,j)*pgmy(i,j)-kmy(i,j+1)*pgmy(i,j+1))*dtsrgdx 
    128 ! calcul des diagonales       
    129           arelax(i,j)=-kmx(i,j)*dtwdx2        ! arelax : diagonale i-1,j  
    130  
    131           brelax(i,j)=-kmx(i+1,j)*dtwdx2      ! brelax : diagonale i+1,j  
    132  
    133           drelax(i,j)=-kmy(i,j)*dtwdx2        ! drelax : diagonale i,j-1  
    134            
    135           erelax(i,j)=-kmy(i,j+1)*dtwdx2      ! drelax : diagonale i,j+1  
    136  
    137           crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2  
    138                                               !crelax : diagonale i,j  
     124             ! calcul du vecteur 
     125             FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT 
     126             frelax(i,j)=frelax(i,j)+(kmx(i,j)*pgmx(i,j)-kmx(i+1,j)*pgmx(i+1,j))*dtsrgdx 
     127             frelax(i,j)=frelax(i,j)+(kmy(i,j)*pgmy(i,j)-kmy(i,j+1)*pgmy(i,j+1))*dtsrgdx 
     128             ! calcul des diagonales       
     129             arelax(i,j)=-kmx(i,j)*dtwdx2        ! arelax : diagonale i-1,j  
     130 
     131             brelax(i,j)=-kmx(i+1,j)*dtwdx2      ! brelax : diagonale i+1,j  
     132 
     133             drelax(i,j)=-kmy(i,j)*dtwdx2        ! drelax : diagonale i,j-1  
     134 
     135             erelax(i,j)=-kmy(i,j+1)*dtwdx2      ! drelax : diagonale i,j+1  
     136 
     137             crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2  
     138             !crelax : diagonale i,j  
    139139          else if (klimit(i,j).eq.1) then 
    140140             hwater(i,j)=limit_hw(i,j) 
    141 !             write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j) 
     141             !             write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j) 
    142142          endif 
    143143       end do 
    144144    end do 
    145 !$OMP END DO 
    146 !$OMP END PARALLEL 
     145    !$OMP END DO 
     146    !$OMP END PARALLEL 
    147147 
    148148    ! Boucle de relaxation : 
    … …  
    154154    Do  WHILE(.NOT.STOPP) 
    155155       ntour=ntour+1 
    156 !       write(6,*) 'boucle de relaxation numero',ntour 
     156       !       write(6,*) 'boucle de relaxation numero',ntour 
    157157       !$OMP PARALLEL 
    158158       !$OMP DO PRIVATE(reste) 
    … …  
    186186       Delh=0 
    187187       Vh=0 
    188         
     188 
    189189       !$OMP DO REDUCTION(+:Delh) 
    190190       DO j=2,NYY-1 
    191191          DO i=2,NXX-1 
    192               
     192 
    193193             !   write(166,*) I,J,delh,deltah(i,j) 
    194194             Delh=Delh+deltah(i,j)**2 
    … …  
    198198       !$OMP END DO 
    199199       !$OMP END PARALLEL 
    200 !       write(6,*) delh,maxval(deltah),minval(deltah) 
     200       !       write(6,*) delh,maxval(deltah),minval(deltah) 
    201201       !      testh=SQRT(Delh/Vh) 
    202202       if (delh.gt.0.) then 
    … …  
    206206       endif 
    207207       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000) 
    208 !       write(6,*) ntour,testh 
    209  
    210  
    211 362   continue   
    212 end do 
     208       !       write(6,*) ntour,testh 
     209 
     210 
     211362    continue   
     212    end do 
    213213  end subroutine relaxation_waterdif 
    214214end module relaxation_waterdif_mod 

  • branches/GRISLIv3/SOURCES/resol_adv_diff_2D-sept2009.f90

    r285 r427  
    1313module reso_adv_diff_2D_vect 
    1414 
    15 use module3D_phy 
    16  
    17  
    18  
    19 implicit none 
    20 double precision :: omega   !< parametre schema temporel de la resolution partie diffusion 
    21 double precision :: mu_adv  !< parametre schema temporel de la resolution advection 
    22 double precision :: upwind  !< schema spatial pour l'advection 
    23  
    24 ! tableaux de travail. resolution M H = Frelax 
    25 double precision,dimension(nx,ny) :: crelax      !< diagnonale de M 
    26 double precision,dimension(nx,ny) :: arelax      !< sous diagonale selon x 
    27 double precision,dimension(nx,ny) :: brelax      !< sur diagonale selon x 
    28 double precision,dimension(nx,ny) :: drelax      !< sous diagonale selon y 
    29 double precision,dimension(nx,ny) :: erelax      !< sur diagonale selon y 
    30 double precision,dimension(nx,ny) :: frelax      !< vecteur 
    31 double precision,dimension(nx,ny) :: c_west      !< sur demi mailles Ux 
    32 double precision,dimension(nx,ny) :: c_east      !< sur demi mailles Ux 
    33 double precision,dimension(nx,ny) :: c_north     !< sur demi mailles Uy 
    34 double precision,dimension(nx,ny) :: c_south     !< sur demi mailles Uy 
    35  
    36 double precision,dimension(nx,ny) :: bdx         !< pente socle 
    37 double precision,dimension(nx,ny) :: bdy         !< pente socle 
    38  
    39 double precision,dimension(nx,ny) :: hdx         !< pente epaisseur 
    40 double precision,dimension(nx,ny) :: hdy         !< pente epaisseur 
    41  
    42  
    43  
    44  
    45 integer, dimension(2) :: ijmax    ! position de maxval 
    46 integer :: iFAIL 
     15  use geography, only:nx,ny 
     16   
     17  implicit none 
     18 
     19  double precision :: omega   !< parametre schema temporel de la resolution partie diffusion 
     20  double precision :: mu_adv  !< parametre schema temporel de la resolution advection 
     21  double precision :: upwind  !< schema spatial pour l'advection 
     22 
     23  ! tableaux de travail. resolution M H = Frelax 
     24  double precision,dimension(nx,ny) :: crelax      !< diagnonale de M 
     25  double precision,dimension(nx,ny) :: arelax      !< sous diagonale selon x 
     26  double precision,dimension(nx,ny) :: brelax      !< sur diagonale selon x 
     27  double precision,dimension(nx,ny) :: drelax      !< sous diagonale selon y 
     28  double precision,dimension(nx,ny) :: erelax      !< sur diagonale selon y 
     29  double precision,dimension(nx,ny) :: frelax      !< vecteur 
     30  double precision,dimension(nx,ny) :: c_west      !< sur demi mailles Ux 
     31  double precision,dimension(nx,ny) :: c_east      !< sur demi mailles Ux 
     32  double precision,dimension(nx,ny) :: c_north     !< sur demi mailles Uy 
     33  double precision,dimension(nx,ny) :: c_south     !< sur demi mailles Uy 
     34  double precision,dimension(nx,ny) :: bdx         !< pente socle 
     35  double precision,dimension(nx,ny) :: bdy         !< pente socle 
     36  double precision,dimension(nx,ny) :: hdx         !< pente epaisseur 
     37  double precision,dimension(nx,ny) :: hdy         !< pente epaisseur 
     38  integer, dimension(2) :: ijmax    ! position de maxval 
     39  integer :: iFAIL 
    4740 
    4841contains  
    4942 
    50 !> initialise  init_reso_adv_diff_2D 
    51 !! definition des parametres qui gerent la resolution 
    52 !! @return omega, mu_adv, upwind 
    53  
    54 subroutine init_reso_adv_diff_2D 
    55  
    56  
    57 !write(num_rep_42,*)'partie diffusive' 
    58 !write(num_rep_42,*)'le type de schema temporel diffusif depend de omega' 
    59 !write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' 
    60 !write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite'  
    61  
    62 ! over implicite propose par Hindmasrsh 
    63 omega = 2.5 
    64 !write(num_rep_42,*)'omega = ',omega 
    65 !write(num_rep_42,*) 
    66 !write(num_rep_42,*)'partie advective' 
    67 !write(num_rep_42,*)' le schéma temporel advectif dépend de mu' 
    68 !write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' 
    69 !write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite' 
    70 mu_adv=1.  ! l'over implicite ne s'applique pas a l'equation advective 
    71 !write(num_rep_42,*)'mu_adv = ',mu_adv 
    72 !write(num_rep_42,*)' le schéma spatial dépend de upwind' 
    73 !write(num_rep_42,*)'1 -> schema upwind,  0.5 -> schema centre' 
    74 upwind=1 
    75 !write(num_rep_42,*)'upwind = ',upwind 
    76 !write(num_rep_42,*)'------------------------------------------------------' 
    77  
    78 return 
    79 end subroutine init_reso_adv_diff_2D 
    80 !------------------------------------------------------------------ 
    81  
    82 !> subroutine resol_adv_diff_2D_Hpresc 
    83 !! definition des parametres qui gerent la resolution 
    84 !! @param Dfx,Dfy             termes diffusifs 
    85 !! @param advx,advy           termes advectifs 
    86 !! @param vieuxH              ancienne valeur de H 
    87 !! @param H_presc             la valeur H prescrit pour certains noeuds 
    88 !! @param i_Hpresc            le masque ou la valeur de H est prescrite 
    89 !! @param bil                 le bilan pour la colonne  
    90 !! @return newH               la valeur de H apres le pas de temps 
    91  
    92 subroutine resol_adv_diff_2D_vect(Dfx,Dfy,advx,advy,H_presc,i_Hpresc,bil,vieuxH,newH) 
    93  
    94 !$ USE OMP_LIB 
    95 implicit none 
    96  
    97 double precision,dimension(nx,ny), intent(in) :: Dfx          !< terme diffusif selon x 
    98 double precision,dimension(nx,ny), intent(in) :: Dfy          !< terme diffusif selon y 
    99 double precision,dimension(nx,ny), intent(in) :: Advx         !< terme advectif selon x 
    100 double precision,dimension(nx,ny), intent(in) :: Advy         !< terme advectif selon y 
    101 double precision,dimension(nx,ny), intent(in) :: vieuxH       !< ancienne valeur de H 
    102 double precision,dimension(nx,ny), intent(out):: newH         !< nouvelle valeur de H 
    103  
    104 double precision,dimension(nx,ny), intent(in) :: bil          !< bilan de masse pour la colonne 
    105  
    106 double precision,dimension(nx,ny), intent(in) :: H_presc      !< H value if prescribed 
    107 integer,dimension(nx,ny), intent(in) :: i_Hpresc  !< 1 if H is prescribed on this node, else 0 
    108  
    109  
    110 ! tableaux de travail. resolution M H = Frelax 
    111 !!$real,dimension(nx,ny) :: crelax      !< diagnonale de M 
    112 !!$real,dimension(nx,ny) :: arelax      !< sous diagonale selon x 
    113 !!$real,dimension(nx,ny) :: brelax      !< sur diagonale selon x 
    114 !!$real,dimension(nx,ny) :: drelax      !< sous diagonale selon y 
    115 !!$real,dimension(nx,ny) :: erelax      !< sur diagonale selon y 
    116 !!$real,dimension(nx,ny) :: frelax      !< vecteur 
    117 !!$real,dimension(nx,ny) :: c_west      !< sur demi mailles Ux 
    118 !!$real,dimension(nx,ny) :: c_east      !< sur demi mailles Ux 
    119 !!$real,dimension(nx,ny) :: c_north     !< sur demi mailles Uy 
    120 !!$real,dimension(nx,ny) :: c_south     !<sur demi mailles Ux 
    121  
    122 !!$real,dimension(nx,ny) :: bdx         !< pente socle 
    123 !!$real,dimension(nx,ny) :: bdy         !< pente socle 
    124 !!$ 
    125 !!$real,dimension(nx,ny) :: hdx         !< pente epaisseur 
    126 !!$real,dimension(nx,ny) :: hdy         !< pente epaisseur 
    127  
    128 double precision :: frdx,frdy                    !< pour calcul frelax : termes diffusion 
    129 double precision :: fraxw,fraxe,frays,frayn      !< termes advection 
    130  
    131 double precision,dimension(nx,ny) :: deltah      ! dans calcul relax 
    132 double precision :: delh                         ! dans calcul relax 
    133 double precision :: testh                        ! dans calcul relax 
    134  
    135 logical :: stopp 
    136 integer :: ntour 
    137 double precision :: reste 
    138  
    139  
    140 if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion') 
    141  
    142 !$OMP PARALLEL 
    143 !$OMP WORKSHARE 
    144 ! calcul de bdx et hdx 
    145 hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0d0,dim=1)) 
    146 hdy(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0d0,dim=2)) 
    147 bdx(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0d0,dim=1)) 
    148 bdy(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0d0,dim=2)) 
    149  
    150  
    151 ! initialisations (qui feront aussi les conditions aux limites) 
    152 Arelax(:,:)=0. 
    153 Brelax(:,:)=0. 
    154 Drelax(:,:)=0. 
    155 Erelax(:,:)=0. 
    156 Crelax(:,:)=1. 
    157 Frelax(:,:)=0. 
    158 DeltaH(:,:)=0. 
    159 !$OMP END WORKSHARE 
    160  
    161 ! schema spatial 
    162  
    163 if (upwind.eq.1.) then                 !schema amont 
    164 !$OMP WORKSHARE 
    165    where (Advx(:,:).ge.0.) 
    166       c_west(:,:)=1. 
    167       c_east(:,:)=0. 
    168    elsewhere 
    169       c_west(:,:)=0. 
    170       c_east(:,:)=1. 
    171    end where 
    172  
    173    where (Advy(:,:).ge.0.) 
    174       c_south(:,:)=1. 
    175       c_north(:,:)=0. 
    176    elsewhere 
    177       c_south(:,:)=0. 
    178       c_north(:,:)=1. 
    179    end where 
    180 !$OMP END WORKSHARE 
    181 else if (upwind.lt.1.) then             ! schema centre 
    182 !$OMP WORKSHARE 
    183       c_west(:,:)=0.5 
    184       c_east(:,:)=0.5 
    185       c_south(:,:)=0.5 
    186       c_north(:,:)=0.5 
    187 !$OMP END WORKSHARE       
    188 end if 
    189  
    190  
    191 ! attribution des elements des diagonales 
    192 !$OMP DO PRIVATE(frdx,frdy,fraxw,fraxe,frays,frayn) 
    193 do j=2,ny-1 
    194   do i=2,nx-1 
    195  
    196 !  sous diagonale en x 
    197      arelax(i,j)=-omega*Dtdx2*Dfx(i,j)   &               ! partie diffusive en x 
    198           -mu_adv*dtdx*c_west(i,j)*Advx(i,j)             ! partie advective en x 
    199  
    200 !  sur diagonale en x 
    201      brelax(i,j)=-omega*Dtdx2*Dfx(i+1,j)  &              ! partie diffusive 
    202           +mu_adv*dtdx*c_east(i+1,j)*Advx(i+1,j)         ! partie advective 
    203  
    204 !  sous diagonale en y 
    205      drelax(i,j)=-omega*Dtdx2*Dfy(i,j)   &               ! partie diffusive en y 
    206           -mu_adv*dtdx*c_south(i,j)*Advy(i,j)            ! partie advective en y 
    207  
    208 !  sur diagonale en y 
    209      erelax(i,j)=-omega*Dtdx2*Dfy(i,j+1)  &              ! partie diffusive 
    210           +mu_adv*dtdx*c_north(i,j+1)*Advy(i,j+1)        ! partie advective 
    211  
    212  
    213  
    214 ! diagonale 
    215      crelax(i,j)=omega*Dtdx2*((Dfx(i+1,j)+Dfx(i,j))  &                         ! partie diffusive en x 
    216                              +(Dfy(i,j+1)+Dfy(i,j)))                           ! partie diffusive en y 
    217      crelax(i,j)=crelax(i,j)+mu_adv*dtdx* & 
    218                   ( (c_west(i+1,j)*Advx(i+1,j)-c_east(i,j)*Advx(i,j)) &        !partie advective en x 
    219                   +(c_south(i,j+1)*Advy(i,j+1)-c_north(i,j)*Advy(i,j)))        !partie advective en y 
    220      crelax(i,j)=1.+crelax(i,j)                                                ! partie temporelle 
    221  
    222  
    223 ! terme du vecteur 
    224  
    225      frdx= -Dfx(i,j) * (Bdx(i,j)  +(1.-omega)*Hdx(i,j))         &  ! partie diffusive en x 
    226            +Dfx(i+1,j)*(Bdx(i+1,j)+(1.-omega)*Hdx(i+1,j)) 
    227  
    228      frdy= -Dfy(i,j) * (Bdy(i,j)  +(1.-omega)*Hdy(i,j))         &  ! partie diffusive en y 
    229            +Dfy(i,j+1)*(Bdy(i,j+1)+(1.-omega)*Hdy(i,j+1)) 
    230  
    231      fraxw= -c_west(i,j)*  Advx(i,j) * vieuxH(i-1,j)           &  ! partie advective en x  
    232             +c_west(i+1,j)*Advx(i+1,j)*vieuxH(i,j)                ! venant de l'west 
    233  
    234      fraxe= -c_east(i,j) * Advx(i,j) * vieuxH(i,j)             &  ! partie advective en x  
    235             +c_east(i+1,j)*Advx(i+1,j)*vieuxH(i+1,j)              ! venant de l'est 
    236  
    237      frays= -c_south(i,j) * Advy(i,j) * vieuxH(i,j-1)          &  ! partie advective en y 
    238             +c_south(i,j+1)*Advy(i,j+1)*vieuxH(i,j)               ! venant du sud 
    239  
    240      frayn= -c_north(i,j) * Advy(i,j) * vieuxH(i,j)            &  ! partie advective en y 
    241             +c_north(i,j+1)*Advy(i,j+1)*vieuxH(i,j+1)             ! venant du nord 
    242  
    243  
    244  
    245      frelax(i,j) = vieuxH(i,j) + Dt*bil(i,j)    & 
    246                  + dtdx*(frdx+frdy)        &  ! termes lies a la diffusion (socle, omega) 
    247                  + (1.-mu_adv)*dtdx*((fraxw+fraxe)+(frays+frayn))  
    248                                               ! le dernier terme serait pour  
    249                                               ! de l'over implicite en diffusion 
    250                                               ! (en pratique pas utilise) 
    251      end do 
    252   end do 
    253 !$OMP END DO 
    254  
    255 !$OMP WORKSHARE 
    256 where (i_hpresc(:,:) .eq.1)          ! thickness prescribed 
    257    frelax(:,:) = H_presc(:,:) 
    258    arelax(:,:) = 0. 
    259    brelax(:,:) = 0. 
    260    crelax(:,:) = 1. 
    261    drelax(:,:) = 0. 
    262    erelax(:,:) = 0. 
    263 end where 
    264 !$OMP END WORKSHARE 
    265 !$OMP END PARALLEL 
    266  
    267 stopp = .false. 
    268 ntour=0 
    269  
    270 relax_loop: do while(.not.stopp) 
    271    ntour=ntour+1 
    272    !$OMP PARALLEL 
    273    !$OMP DO PRIVATE(reste) 
    274    do j=2,ny-1 
    275       do i=2,nx-1 
    276          reste = (((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) & 
    277                + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) & 
    278                + crelax(i,j)*newH(i,j))- frelax(i,j) 
    279  
    280          deltaH(i,j) = reste/crelax(i,j) 
    281       end do 
    282    end do 
    283    !$OMP END DO 
    284  
    285    !$OMP WORKSHARE 
    286    newH(:,:)=newH(:,:)-deltaH(:,:) 
    287    !$OMP END WORKSHARE 
    288    !$OMP END PARALLEL 
    289  
    290  
    291 ! critere d'arret: 
    292 ! ----------------          
    293    delh=0 
    294  
    295    !$OMP PARALLEL 
    296    !$OMP DO REDUCTION(+:delh) 
    297    do j=2,ny-1 
    298       do i=2,nx-1 
    299          delh=delh+deltaH(i,j)**2 
    300       end do 
    301    end do 
    302    !$OMP END DO 
    303    !$OMP END PARALLEL 
    304  
    305    if (delh.gt.0.) then 
    306       testh=sqrt(delh)/((nx-2)*(ny-2)) 
    307    else 
    308       testh=0. 
    309    endif 
    310    stopp = (testh.lt.1.e-4).or.(ntour.gt.100) 
    311  
    312 end do relax_loop 
    313  
    314 if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion') 
    315   
    316 return 
    317  
    318 end subroutine resol_adv_diff_2D_vect 
     43  !> initialise  init_reso_adv_diff_2D 
     44  !! definition des parametres qui gerent la resolution 
     45  !! @return omega, mu_adv, upwind 
     46 
     47  subroutine init_reso_adv_diff_2D 
     48 
     49    !write(num_rep_42,*)'partie diffusive' 
     50    !write(num_rep_42,*)'le type de schema temporel diffusif depend de omega' 
     51    !write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' 
     52    !write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite'  
     53 
     54    ! over implicite propose par Hindmasrsh 
     55    omega = 2.5 
     56    !write(num_rep_42,*)'omega = ',omega 
     57    !write(num_rep_42,*) 
     58    !write(num_rep_42,*)'partie advective' 
     59    !write(num_rep_42,*)' le schéma temporel advectif dépend de mu' 
     60    !write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' 
     61    !write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite' 
     62    mu_adv=1.  ! l'over implicite ne s'applique pas a l'equation advective 
     63    !write(num_rep_42,*)'mu_adv = ',mu_adv 
     64    !write(num_rep_42,*)' le schéma spatial dépend de upwind' 
     65    !write(num_rep_42,*)'1 -> schema upwind,  0.5 -> schema centre' 
     66    upwind=1 
     67    !write(num_rep_42,*)'upwind = ',upwind 
     68    !write(num_rep_42,*)'------------------------------------------------------' 
     69 
     70    return 
     71  end subroutine init_reso_adv_diff_2D 
     72  !------------------------------------------------------------------ 
     73 
     74  !> subroutine resol_adv_diff_2D_Hpresc 
     75  !! definition des parametres qui gerent la resolution 
     76  !! @param Dfx,Dfy             termes diffusifs 
     77  !! @param advx,advy           termes advectifs 
     78  !! @param vieuxH              ancienne valeur de H 
     79  !! @param H_presc             la valeur H prescrit pour certains noeuds 
     80  !! @param i_Hpresc            le masque ou la valeur de H est prescrite 
     81  !! @param bil                 le bilan pour la colonne  
     82  !! @return newH               la valeur de H apres le pas de temps 
     83 
     84  subroutine resol_adv_diff_2D_vect(Dfx,Dfy,advx,advy,H_presc,i_Hpresc,bil,vieuxH,newH) 
     85 
     86    use module3D_phy, only: dx1,B,dtdx2,dtdx,dt 
     87    use runparam, only: itracebug 
     88    !$ USE OMP_LIB 
     89    implicit none 
     90 
     91    double precision,dimension(nx,ny), intent(in) :: Dfx          !< terme diffusif selon x 
     92    double precision,dimension(nx,ny), intent(in) :: Dfy          !< terme diffusif selon y 
     93    double precision,dimension(nx,ny), intent(in) :: Advx         !< terme advectif selon x 
     94    double precision,dimension(nx,ny), intent(in) :: Advy         !< terme advectif selon y 
     95    double precision,dimension(nx,ny), intent(in) :: vieuxH       !< ancienne valeur de H 
     96    double precision,dimension(nx,ny), intent(out):: newH         !< nouvelle valeur de H 
     97    double precision,dimension(nx,ny), intent(in) :: bil          !< bilan de masse pour la colonne 
     98    double precision,dimension(nx,ny), intent(in) :: H_presc      !< H value if prescribed 
     99    integer,dimension(nx,ny), intent(in) :: i_Hpresc  !< 1 if H is prescribed on this node, else 0 
     100 
     101    ! local variables 
     102    double precision :: frdx,frdy                    !< pour calcul frelax : termes diffusion 
     103    double precision :: fraxw,fraxe,frays,frayn      !< termes advection 
     104    double precision,dimension(nx,ny) :: deltah      ! dans calcul relax 
     105    double precision :: delh                         ! dans calcul relax 
     106    double precision :: testh                        ! dans calcul relax 
     107    logical :: stopp 
     108    integer :: ntour 
     109    double precision :: reste 
     110    integer :: i,j 
     111 
     112 
     113    if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion') 
     114 
     115    !$OMP PARALLEL 
     116    !$OMP WORKSHARE 
     117    ! calcul de bdx et hdx 
     118    hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0d0,dim=1)) 
     119    hdy(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0d0,dim=2)) 
     120    bdx(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0d0,dim=1)) 
     121    bdy(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0d0,dim=2)) 
     122 
     123 
     124    ! initialisations (qui feront aussi les conditions aux limites) 
     125    Arelax(:,:)=0. 
     126    Brelax(:,:)=0. 
     127    Drelax(:,:)=0. 
     128    Erelax(:,:)=0. 
     129    Crelax(:,:)=1. 
     130    Frelax(:,:)=0. 
     131    DeltaH(:,:)=0. 
     132    !$OMP END WORKSHARE 
     133 
     134    ! schema spatial 
     135 
     136    if (upwind.eq.1.) then                 !schema amont 
     137       !$OMP WORKSHARE 
     138       where (Advx(:,:).ge.0.) 
     139          c_west(:,:)=1. 
     140          c_east(:,:)=0. 
     141       elsewhere 
     142          c_west(:,:)=0. 
     143          c_east(:,:)=1. 
     144       end where 
     145 
     146       where (Advy(:,:).ge.0.) 
     147          c_south(:,:)=1. 
     148          c_north(:,:)=0. 
     149       elsewhere 
     150          c_south(:,:)=0. 
     151          c_north(:,:)=1. 
     152       end where 
     153       !$OMP END WORKSHARE 
     154    else if (upwind.lt.1.) then             ! schema centre 
     155       !$OMP WORKSHARE 
     156       c_west(:,:)=0.5 
     157       c_east(:,:)=0.5 
     158       c_south(:,:)=0.5 
     159       c_north(:,:)=0.5 
     160       !$OMP END WORKSHARE       
     161    end if 
     162 
     163 
     164    ! attribution des elements des diagonales 
     165    !$OMP DO PRIVATE(frdx,frdy,fraxw,fraxe,frays,frayn) 
     166    do j=2,ny-1 
     167       do i=2,nx-1 
     168 
     169          !  sous diagonale en x 
     170          arelax(i,j)=-omega*Dtdx2*Dfx(i,j)   &               ! partie diffusive en x 
     171               -mu_adv*dtdx*c_west(i,j)*Advx(i,j)             ! partie advective en x 
     172 
     173          !  sur diagonale en x 
     174          brelax(i,j)=-omega*Dtdx2*Dfx(i+1,j)  &              ! partie diffusive 
     175               +mu_adv*dtdx*c_east(i+1,j)*Advx(i+1,j)         ! partie advective 
     176 
     177          !  sous diagonale en y 
     178          drelax(i,j)=-omega*Dtdx2*Dfy(i,j)   &               ! partie diffusive en y 
     179               -mu_adv*dtdx*c_south(i,j)*Advy(i,j)            ! partie advective en y 
     180 
     181          !  sur diagonale en y 
     182          erelax(i,j)=-omega*Dtdx2*Dfy(i,j+1)  &              ! partie diffusive 
     183               +mu_adv*dtdx*c_north(i,j+1)*Advy(i,j+1)        ! partie advective 
     184 
     185 
     186 
     187          ! diagonale 
     188          crelax(i,j)=omega*Dtdx2*((Dfx(i+1,j)+Dfx(i,j))  &                         ! partie diffusive en x 
     189               +(Dfy(i,j+1)+Dfy(i,j)))                           ! partie diffusive en y 
     190          crelax(i,j)=crelax(i,j)+mu_adv*dtdx* & 
     191               ( (c_west(i+1,j)*Advx(i+1,j)-c_east(i,j)*Advx(i,j)) &        !partie advective en x 
     192               +(c_south(i,j+1)*Advy(i,j+1)-c_north(i,j)*Advy(i,j)))        !partie advective en y 
     193          crelax(i,j)=1.+crelax(i,j)                                                ! partie temporelle 
     194 
     195 
     196          ! terme du vecteur 
     197 
     198          frdx= -Dfx(i,j) * (Bdx(i,j)  +(1.-omega)*Hdx(i,j))         &  ! partie diffusive en x 
     199               +Dfx(i+1,j)*(Bdx(i+1,j)+(1.-omega)*Hdx(i+1,j)) 
     200 
     201          frdy= -Dfy(i,j) * (Bdy(i,j)  +(1.-omega)*Hdy(i,j))         &  ! partie diffusive en y 
     202               +Dfy(i,j+1)*(Bdy(i,j+1)+(1.-omega)*Hdy(i,j+1)) 
     203 
     204          fraxw= -c_west(i,j)*  Advx(i,j) * vieuxH(i-1,j)           &  ! partie advective en x  
     205               +c_west(i+1,j)*Advx(i+1,j)*vieuxH(i,j)                ! venant de l'west 
     206 
     207          fraxe= -c_east(i,j) * Advx(i,j) * vieuxH(i,j)             &  ! partie advective en x  
     208               +c_east(i+1,j)*Advx(i+1,j)*vieuxH(i+1,j)              ! venant de l'est 
     209 
     210          frays= -c_south(i,j) * Advy(i,j) * vieuxH(i,j-1)          &  ! partie advective en y 
     211               +c_south(i,j+1)*Advy(i,j+1)*vieuxH(i,j)               ! venant du sud 
     212 
     213          frayn= -c_north(i,j) * Advy(i,j) * vieuxH(i,j)            &  ! partie advective en y 
     214               +c_north(i,j+1)*Advy(i,j+1)*vieuxH(i,j+1)             ! venant du nord 
     215 
     216 
     217 
     218          frelax(i,j) = vieuxH(i,j) + Dt*bil(i,j)    & 
     219               + dtdx*(frdx+frdy)        &  ! termes lies a la diffusion (socle, omega) 
     220               + (1.-mu_adv)*dtdx*((fraxw+fraxe)+(frays+frayn))  
     221          ! le dernier terme serait pour  
     222          ! de l'over implicite en diffusion 
     223          ! (en pratique pas utilise) 
     224       end do 
     225    end do 
     226    !$OMP END DO 
     227 
     228    !$OMP WORKSHARE 
     229    where (i_hpresc(:,:) .eq.1)          ! thickness prescribed 
     230       frelax(:,:) = H_presc(:,:) 
     231       arelax(:,:) = 0. 
     232       brelax(:,:) = 0. 
     233       crelax(:,:) = 1. 
     234       drelax(:,:) = 0. 
     235       erelax(:,:) = 0. 
     236    end where 
     237    !$OMP END WORKSHARE 
     238    !$OMP END PARALLEL 
     239 
     240    stopp = .false. 
     241    ntour=0 
     242 
     243    relax_loop: do while(.not.stopp) 
     244       ntour=ntour+1 
     245       !$OMP PARALLEL 
     246       !$OMP DO PRIVATE(reste) 
     247       do j=2,ny-1 
     248          do i=2,nx-1 
     249             reste = (((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) & 
     250                  + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) & 
     251                  + crelax(i,j)*newH(i,j))- frelax(i,j) 
     252 
     253             deltaH(i,j) = reste/crelax(i,j) 
     254          end do 
     255       end do 
     256       !$OMP END DO 
     257 
     258       !$OMP WORKSHARE 
     259       newH(:,:)=newH(:,:)-deltaH(:,:) 
     260       !$OMP END WORKSHARE 
     261       !$OMP END PARALLEL 
     262 
     263 
     264       ! critere d'arret: 
     265       ! ----------------          
     266       delh=0 
     267 
     268       !$OMP PARALLEL 
     269       !$OMP DO REDUCTION(+:delh) 
     270       do j=2,ny-1 
     271          do i=2,nx-1 
     272             delh=delh+deltaH(i,j)**2 
     273          end do 
     274       end do 
     275       !$OMP END DO 
     276       !$OMP END PARALLEL 
     277 
     278       if (delh.gt.0.) then 
     279          testh=sqrt(delh)/((nx-2)*(ny-2)) 
     280       else 
     281          testh=0. 
     282       endif 
     283       stopp = (testh.lt.1.e-4).or.(ntour.gt.100) 
     284 
     285    end do relax_loop 
     286 
     287    if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion') 
     288 
     289    return 
     290 
     291  end subroutine resol_adv_diff_2D_vect 
    319292 
    320293end module reso_adv_diff_2D_vect 
Note: See TracChangeset for help on using the changeset viewer.