Changeset 395


Ignore:
Timestamp:
03/24/23 11:55:20 (13 months ago)
Author:
dumas
Message:

use only in module eau_basale

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/eaubasale-0.5_mod.f90

    r237 r395  
    1515module eau_basale 
    1616 
    17 use module3d_phy 
    18 use param_phy_mod 
    19 use relaxation_waterdif_mod ! module qui contient la routine relaxation  
    20                            ! pour l'eau avec interface explicite 
    21  
    22  
    23 implicit none 
    24  
    25 ! dimensionnements 
    26 !------------------------------------------------------------------------  
    27     LOGICAL :: ecoulement_eau 
    28  
    29     REAL :: KONDMAX 
    30     real :: kond0 
    31     REAL :: INFILTR 
    32     REAL :: hmax_till    !< épaisseur de la couche de till 
    33     REAL :: hmax_wat     !< épaisseur de la couche d'eau dans le till 
    34     REAL :: poro_till    !< porosité du till 
    35     REAL :: Zflot        !< hauteur d'eau donnant la flottaison 
    36     real :: keffmax      !< kondmax*hmax_wat 
    37     real :: nefflocal 
    38  
    39     REAL,dimension(NX,NY) :: limit_hw    !< conditions aux limites 
    40     integer,dimension(NX,NY) :: klimit    !< ou appliquer les conditions  
    41     REAL,dimension(NX,NY) :: pot_w             !< pour calculer le gradient de pression  
    42     REAL,dimension(NX,NY) :: pot_f             !< pour les points flottants 
    43     REAL,dimension(NX,NY) :: hw                !< hauteur d'eau dans le till, limite a hmax_wat 
    44                                                  !< c'est la hauteur sur laquelle peut se faire  
    45                                                  !< l'ecoulement de l'eau 
    46     REAL,dimension(nx,ny) :: keff              !< conductivite effective : Kond*hw 
    47  
    48  
    49     REAL,dimension(NX,NY) :: bmelt_w          !< fusion (terme source) exprimé en m d'eau 
    50     REAL,dimension(NX,NY) :: vieuxhwater      !< valeur de hwater au debut de l'appel 
    51  
    52  
    53     INTEGER :: NXlocal,NYlocal !< pour passer NX et NY en parametre a la routine relaxation 
     17  use module3d_phy, only:nx,ny 
     18  use param_phy_mod, only: 
     19 
     20  implicit none 
     21 
     22  ! dimensionnements 
     23  !------------------------------------------------------------------------  
     24  LOGICAL :: ecoulement_eau 
     25 
     26  REAL :: KONDMAX 
     27  real :: kond0 
     28  REAL :: INFILTR 
     29  REAL :: hmax_till    !< épaisseur de la couche de till 
     30  REAL :: hmax_wat     !< épaisseur de la couche d'eau dans le till 
     31  REAL :: poro_till    !< porosité du till 
     32  REAL :: Zflot        !< hauteur d'eau donnant la flottaison 
     33  real :: keffmax      !< kondmax*hmax_wat 
     34  real :: nefflocal 
     35 
     36  REAL,dimension(NX,NY) :: limit_hw    !< conditions aux limites 
     37  integer,dimension(NX,NY) :: klimit    !< ou appliquer les conditions  
     38  REAL,dimension(NX,NY) :: pot_w             !< pour calculer le gradient de pression  
     39  REAL,dimension(NX,NY) :: pot_f             !< pour les points flottants 
     40  REAL,dimension(NX,NY) :: hw                !< hauteur d'eau dans le till, limite a hmax_wat 
     41  !< c'est la hauteur sur laquelle peut se faire  
     42  !< l'ecoulement de l'eau 
     43  REAL,dimension(nx,ny) :: keff              !< conductivite effective : Kond*hw 
     44 
     45 
     46  REAL,dimension(NX,NY) :: bmelt_w          !< fusion (terme source) exprimé en m d'eau 
     47  REAL,dimension(NX,NY) :: vieuxhwater      !< valeur de hwater au debut de l'appel 
     48 
     49 
     50  INTEGER :: NXlocal,NYlocal !< pour passer NX et NY en parametre a la routine relaxation 
    5451 
    5552 
    5653contains 
    57 !> SUBROUTINE: init_eaubasale 
    58 !!Initialisation du block eaubasale 
    59 !> 
    60 subroutine init_eaubasale 
    61  
    62 namelist/eaubasale1/ecoulement_eau,hwatermax,infiltr  
    63 namelist/param_hydr/hmax_till,poro_till,kond0 
    64  
    65 ! formats pour les ecritures dans 42 
     54  !> SUBROUTINE: init_eaubasale 
     55  !!Initialisation du block eaubasale 
     56  !> 
     57  subroutine init_eaubasale 
     58 
     59    use module3d_phy, only:hwatermax,num_param,num_rep_42,kond,secyear,hdotwater,pgx,pgy 
     60     
     61    namelist/eaubasale1/ecoulement_eau,hwatermax,infiltr  
     62    namelist/param_hydr/hmax_till,poro_till,kond0 
     63 
     64    ! formats pour les ecritures dans 42 
    6665428 format(A) 
    6766 
    68 ! lecture des parametres du run                      block eaubasale1 
    69 !-------------------------------------------------------------------- 
    70  
    71  
    72 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    73 read(num_param,eaubasale1) 
    74 write(num_rep_42,428)'!___________________________________________________________'  
    75 write(num_rep_42,428) '&eaubasale1              ! nom du premier bloc eau basale ' 
    76 write(num_rep_42,*) 
    77 write(num_rep_42,*) 'ecoulement_eau = ',ecoulement_eau 
    78 write(num_rep_42,*) 'hwatermax      = ',hwatermax 
    79 write(num_rep_42,*) 'infiltr        = ', infiltr 
    80 write(num_rep_42,*)'/'                             
    81 write(num_rep_42,428) '! ecoulement eau : .false. -> modele bucket, sinon equ. diffusion' 
    82 write(num_rep_42,428) '! hwatermax :  hauteur d eau basale maximum dans le sediment (m)' 
    83 write(num_rep_42,428) '! infiltr est la quantite d eau qui peut s infiltrer dans le sol (m/an)' 
    84 write(num_rep_42,*) 
    85  
    86  
    87 !valeurs numeriques des parametres hydrauliques 
    88  
    89 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    90 read(num_param,param_hydr) 
    91  
    92 write(num_rep_42,428)'!___________________________________________________________'  
    93 write(num_rep_42,428) '&param_hydr             ! nom du  bloc parametres hydrauliques ' 
    94 write(num_rep_42,*) 
    95 write(num_rep_42,*) 'hmax_till      = ',hmax_till 
    96 write(num_rep_42,*) 'poro_till      = ',poro_till 
    97 write(num_rep_42,*) 'kond0          = ',kond0  
    98 write(num_rep_42,*)'/'                             
    99 write(num_rep_42,428) '! hmax_till (m) : epaisseur max du sediment ' 
    100 write(num_rep_42,428) '! poro_till : porosite du sediment ' 
    101 write(num_rep_42,428) '! conductivite du sediment :  kond0 (m/s)' 
    102 write(num_rep_42,*) 
    103  
    104  
    105  
    106  
    107 hmax_wat=hmax_till*poro_till ! hauteur maxi que peut atteindre l'eau dans la couche de till 
    108  
    109  
    110 ! Conductivite hydraulique : cond passée en m/an ( car le dt est en années) 
    111 kond(:,:)=kond0 
    112 kond(:,:) = kond(:,:)*SECYEAR  
    113 kondmax   = 1.*SECYEAR  
    114 hdotwater(:,:)=0. 
    115 keff(:,:) = kond(:,:) 
    116 pgx(:,:) = 0. 
    117 pgy(:,:) = 0. 
    118  
    119 NXlocal=NX 
    120 NYlocal=NY 
    121  
    122 end subroutine init_eaubasale 
    123  
    124 !> SUBROUTINE: eaubasale 
    125 !! to do 
    126 !> 
    127 subroutine eaubasale !(pwater)   version correspondant à la thèse de Vincent 
    128   ! A terme, il faudrait en faire un module  
    129 !$ USE OMP_LIB 
    130 !~   integer :: t1,t2,ir 
    131 !~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme   
    132 !~    
    133 !~    ! Temps CPU de calcul initial. 
    134 !~    call cpu_time(t_cpu_0) 
    135 !~    ! Temps elapsed de reference. 
    136 !~    call system_clock(count=t1, count_rate=ir) 
    137  
    138  
    139  
    140  
    141   if (itracebug.eq.1)  call tracebug(' Entree dans routine eau basale') 
    142  
    143 !$OMP PARALLEL 
    144 !$OMP WORKSHARE 
    145   vieuxhwater(:,:) = hwater(:,:) 
    146   kond(:,:) = kond0*SECYEAR  
    147  
    148  
    149   ! conditions aux limites 
    150   klimit(:,:)=0 
    151   limit_hw(:,:)=-9999. 
    152 !$OMP END WORKSHARE 
    153 !$OMP DO   
    154   do j=1,ny 
    155      do i=1,nx 
    156  
    157  
    158         if (flot(i,j))then  ! points flottants 
    159            klimit(i,j)=1 
    160            limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg 
    161  
    162         else if (IBASE(I,J).eq.1) then ! base froide 
    163            klimit(i,j)=1 
    164            limit_hw(i,j)=0. 
    165  
    166         else if ((.not.flot(i,j)).and.(H(i,j).lt.1.)) then  ! bord de la calotte 
    167            klimit(i,j)=1 
    168            limit_hw(i,j)=10.    ! riviere de 10 m de profondeur 
    169         endif 
    170  
    171      end do 
    172   end do 
    173 !$OMP END DO 
    174  
    175 !$OMP DO 
    176 !  do j=2,ny-1 
    177   do j=1,ny 
    178      do i=1,nx 
    179         bmelt_w(i,j)=bmelt(i,j)*rofresh/ro     
    180         hwater(i,j)=max(hwater(i,j),0.) 
    181         hw(i,j)=min(hwater(i,j),hmax_wat)  
    182      enddo 
    183   enddo 
    184 !$OMP END DO 
    185 !$OMP END PARALLEL 
    186  
    187   ecoul:  if (ecoulement_eau) then 
    188      !  print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5' 
    189      !   write(6,*) 'entree ecoulement_eau' 
    190      !$OMP PARALLEL 
    191      !$OMP DO 
    192      do j=1,ny 
    193         do i=1,nx               
    194  
    195            !   calcul des potentiels 
    196            pot_w(I,J)=rofreshg*B(I,J)  ! potentiel de gravite 
    197  
    198            ! on ajoute pression glace mais pas la pression d'eau qui est traitée comme diffusion 
    199  
    200            pot_w(I,J)=pot_w(I,J)+rog*H(I,J) 
    201  
    202            pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
    203         enddo 
    204      enddo 
    205 !$OMP END DO 
    206  
    207      ! sorties debug 17 juillet 2007 
    208      debug_3D(:,:,5)=pot_w(:,:) 
    209      debug_3D(:,:,6)=pot_f(:,:) 
    210      debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:) 
    211      debug_3D(:,:,8)=hwater(:,:) 
    212 !$OMP DO 
    213      do j=2,ny 
    214         do i=2,nx 
    215            if (H(i,j).gt.25.) then !afq: pourquoi GT 25 ??? a tester avec des valeurs plus petites 
    216               !           calcul du gradient de pression 
    217               if (flotmx(i,j)) then  
    218                  pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx 
    219               else 
    220                  pgx(i,j)=(pot_w(i-1,j)-pot_w(i,j))/dx 
    221               endif 
    222  
    223               if (flotmy(i,j)) then  
    224                  pgy(i,j)=(pot_f(i,j-1)-pot_f(i,j))/dy 
    225               else 
    226                  pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy 
    227               endif 
    228            endif 
    229            pgx(i,j)=pgx(i,j)/rofreshg   ! pour passer pgx sans unité 
    230            pgy(i,j)=pgy(i,j)/rofreshg       
    231         end do 
    232      end do 
    233 !$OMP END DO 
    234 !$OMP END PARALLEL 
    235  
    236      if (nt.gt.0) then 
    237         if (dt.gt.0.)  then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0 
    238            !------------------------------------------------------------------------- 
    239            ! les points de la grounding line ont une conductivité hydraulique élevée 
    240            ! même si la base est froide.  modif catritz 18 janvier 2005 
    241            !open(166,file='erreur_eau') 
    242 !$OMP PARALLEL 
    243 !$OMP DO 
    244            do j=2,Ny-1 
    245               do i=2,Nx-1 
    246                  ! cond=0 si glace froide et pas sur la grounding line 
    247                  if ((IBASE(i,j).eq.1).and.   & 
    248                       (.not.flot(i,j-1)).and.(.not.flot(i,j+1)).and.  & 
    249                       (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(i,j)=0.! 1.0e-5 
    250  
    251                  ! cond infinie quand epaisseur faible et glace flottante 
    252                  if (flot(i,j).or.H(i,j).le.1.5)  kond(i,j)= kondmax 
    253  
    254                  ! conductivite forte lorsque N est faible (croit à partir de 100 bars)  
    255                  nefflocal=0.91*H(i,j)-hwater(i,j) 
    256                  nefflocal=max(100.,nefflocal) 
    257                  if (nefflocal.le.1000.) then 
    258                     kond(i,j)=kond(i,j)*1000./nefflocal 
    259                  endif 
    260                  kond(i,j)=min(kondmax,kond(i,j))   
    261  
    262                  ! conductivite effective (conductivité * taille du tuyau en m2/an)  
    263                  keff(i,j)=kond(i,j)*hw(i,j) 
    264               end do 
    265            end do 
    266 !$OMP END DO 
    267            ! condition sur les bords de la grille 
    268 !$OMP WORKSHARE 
    269            kond(1,:)=kondmax 
    270            kond(nx,:)=kondmax 
    271            kond(:,1)=kondmax 
    272            kond(:,ny)=kondmax 
    273            vieuxhwater(:,:) = hwater(:,:) 
    274 !$OMP END WORKSHARE 
    275 !$OMP END PARALLEL 
     67    ! lecture des parametres du run                      block eaubasale1 
     68    !-------------------------------------------------------------------- 
     69 
     70 
     71    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     72    read(num_param,eaubasale1) 
     73    write(num_rep_42,428)'!___________________________________________________________'  
     74    write(num_rep_42,428) '&eaubasale1              ! nom du premier bloc eau basale ' 
     75    write(num_rep_42,*) 
     76    write(num_rep_42,*) 'ecoulement_eau = ',ecoulement_eau 
     77    write(num_rep_42,*) 'hwatermax      = ',hwatermax 
     78    write(num_rep_42,*) 'infiltr        = ', infiltr 
     79    write(num_rep_42,*)'/'                             
     80    write(num_rep_42,428) '! ecoulement eau : .false. -> modele bucket, sinon equ. diffusion' 
     81    write(num_rep_42,428) '! hwatermax :  hauteur d eau basale maximum dans le sediment (m)' 
     82    write(num_rep_42,428) '! infiltr est la quantite d eau qui peut s infiltrer dans le sol (m/an)' 
     83    write(num_rep_42,*) 
     84 
     85 
     86    !valeurs numeriques des parametres hydrauliques 
     87 
     88    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     89    read(num_param,param_hydr) 
     90 
     91    write(num_rep_42,428)'!___________________________________________________________'  
     92    write(num_rep_42,428) '&param_hydr             ! nom du  bloc parametres hydrauliques ' 
     93    write(num_rep_42,*) 
     94    write(num_rep_42,*) 'hmax_till      = ',hmax_till 
     95    write(num_rep_42,*) 'poro_till      = ',poro_till 
     96    write(num_rep_42,*) 'kond0          = ',kond0  
     97    write(num_rep_42,*)'/'                             
     98    write(num_rep_42,428) '! hmax_till (m) : epaisseur max du sediment ' 
     99    write(num_rep_42,428) '! poro_till : porosite du sediment ' 
     100    write(num_rep_42,428) '! conductivite du sediment :  kond0 (m/s)' 
     101    write(num_rep_42,*) 
     102 
     103 
     104 
     105 
     106    hmax_wat=hmax_till*poro_till ! hauteur maxi que peut atteindre l'eau dans la couche de till 
     107 
     108 
     109    ! Conductivite hydraulique : cond passée en m/an ( car le dt est en années) 
     110    kond(:,:)=kond0 
     111    kond(:,:) = kond(:,:)*SECYEAR  
     112    kondmax   = 1.*SECYEAR  
     113    hdotwater(:,:)=0. 
     114    keff(:,:) = kond(:,:) 
     115    pgx(:,:) = 0. 
     116    pgy(:,:) = 0. 
     117 
     118    NXlocal=NX 
     119    NYlocal=NY 
     120 
     121  end subroutine init_eaubasale 
     122 
     123  !> SUBROUTINE: eaubasale 
     124  !! to do 
     125  !> 
     126  subroutine eaubasale !(pwater)   version correspondant à la thèse de Vincent 
     127 
     128    use module3d_phy, only:hwater,kond,secyear,flot,sealevel_2D,Bsoc,rofreshg,ibase,S,H,B,bmelt,rofresh,& 
     129         debug_3d,flotmx,flotmy,pgx,pgy,phiwx,phiwy,isynchro,dtt,dx,dy,dt,hdotwater,pwater,hwatermax 
     130    use param_phy_mod, only:rowg,ro,rog 
     131    use runparam, only:itracebug,nt 
     132    use relaxation_waterdif_mod, only:relaxation_waterdif 
     133     
     134    !$ USE OMP_LIB 
     135    !~   integer :: t1,t2,ir 
     136    !~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme   
     137    !~    
     138    !~    ! Temps CPU de calcul initial. 
     139    !~    call cpu_time(t_cpu_0) 
     140    !~    ! Temps elapsed de reference. 
     141    !~    call system_clock(count=t1, count_rate=ir) 
     142 
     143    implicit none 
     144 
     145    integer :: i,j 
     146 
     147    if (itracebug.eq.1)  call tracebug(' Entree dans routine eau basale') 
     148 
     149    !$OMP PARALLEL 
     150    !$OMP WORKSHARE 
     151    vieuxhwater(:,:) = hwater(:,:) 
     152    kond(:,:) = kond0*SECYEAR  
     153 
     154 
     155    ! conditions aux limites 
     156    klimit(:,:)=0 
     157    limit_hw(:,:)=-9999. 
     158    !$OMP END WORKSHARE 
     159    !$OMP DO   
     160    do j=1,ny 
     161       do i=1,nx 
     162 
     163 
     164          if (flot(i,j))then  ! points flottants 
     165             klimit(i,j)=1 
     166             limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg 
     167 
     168          else if (IBASE(I,J).eq.1) then ! base froide 
     169             klimit(i,j)=1 
     170             limit_hw(i,j)=0. 
     171 
     172          else if ((.not.flot(i,j)).and.(H(i,j).lt.1.)) then  ! bord de la calotte 
     173             klimit(i,j)=1 
     174             limit_hw(i,j)=10.    ! riviere de 10 m de profondeur 
     175          endif 
     176 
     177       end do 
     178    end do 
     179    !$OMP END DO 
     180 
     181    !$OMP DO 
     182    !  do j=2,ny-1 
     183    do j=1,ny 
     184       do i=1,nx 
     185          bmelt_w(i,j)=bmelt(i,j)*rofresh/ro     
     186          hwater(i,j)=max(hwater(i,j),0.) 
     187          hw(i,j)=min(hwater(i,j),hmax_wat)  
     188       enddo 
     189    enddo 
     190    !$OMP END DO 
     191    !$OMP END PARALLEL 
     192 
     193    ecoul:  if (ecoulement_eau) then 
     194       !  print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5' 
     195       !   write(6,*) 'entree ecoulement_eau' 
     196       !$OMP PARALLEL 
     197       !$OMP DO 
     198       do j=1,ny 
     199          do i=1,nx               
     200 
     201             !   calcul des potentiels 
     202             pot_w(I,J)=rofreshg*B(I,J)  ! potentiel de gravite 
     203 
     204             ! on ajoute pression glace mais pas la pression d'eau qui est traitée comme diffusion 
     205 
     206             pot_w(I,J)=pot_w(I,J)+rog*H(I,J) 
     207 
     208             pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf 
     209          enddo 
     210       enddo 
     211       !$OMP END DO 
     212 
     213       ! sorties debug 17 juillet 2007 
     214       debug_3D(:,:,5)=pot_w(:,:) 
     215       debug_3D(:,:,6)=pot_f(:,:) 
     216       debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:) 
     217       debug_3D(:,:,8)=hwater(:,:) 
     218       !$OMP DO 
     219       do j=2,ny 
     220          do i=2,nx 
     221             if (H(i,j).gt.25.) then !afq: pourquoi GT 25 ??? a tester avec des valeurs plus petites 
     222                !           calcul du gradient de pression 
     223                if (flotmx(i,j)) then  
     224                   pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx 
     225                else 
     226                   pgx(i,j)=(pot_w(i-1,j)-pot_w(i,j))/dx 
     227                endif 
     228 
     229                if (flotmy(i,j)) then  
     230                   pgy(i,j)=(pot_f(i,j-1)-pot_f(i,j))/dy 
     231                else 
     232                   pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy 
     233                endif 
     234             endif 
     235             pgx(i,j)=pgx(i,j)/rofreshg   ! pour passer pgx sans unité 
     236             pgy(i,j)=pgy(i,j)/rofreshg       
     237          end do 
     238       end do 
     239       !$OMP END DO 
     240       !$OMP END PARALLEL 
     241 
     242       if (nt.gt.0) then 
     243          if (dt.gt.0.)  then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0 
     244             !------------------------------------------------------------------------- 
     245             ! les points de la grounding line ont une conductivité hydraulique élevée 
     246             ! même si la base est froide.  modif catritz 18 janvier 2005 
     247             !open(166,file='erreur_eau') 
     248             !$OMP PARALLEL 
     249             !$OMP DO 
     250             do j=2,Ny-1 
     251                do i=2,Nx-1 
     252                   ! cond=0 si glace froide et pas sur la grounding line 
     253                   if ((IBASE(i,j).eq.1).and.   & 
     254                        (.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 
     256 
     257                   ! cond infinie quand epaisseur faible et glace flottante 
     258                   if (flot(i,j).or.H(i,j).le.1.5)  kond(i,j)= kondmax 
     259 
     260                   ! conductivite forte lorsque N est faible (croit à partir de 100 bars)  
     261                   nefflocal=0.91*H(i,j)-hwater(i,j) 
     262                   nefflocal=max(100.,nefflocal) 
     263                   if (nefflocal.le.1000.) then 
     264                      kond(i,j)=kond(i,j)*1000./nefflocal 
     265                   endif 
     266                   kond(i,j)=min(kondmax,kond(i,j))   
     267 
     268                   ! conductivite effective (conductivité * taille du tuyau en m2/an)  
     269                   keff(i,j)=kond(i,j)*hw(i,j) 
     270                end do 
     271             end do 
     272             !$OMP END DO 
     273             ! condition sur les bords de la grille 
     274             !$OMP WORKSHARE 
     275             kond(1,:)=kondmax 
     276             kond(nx,:)=kondmax 
     277             kond(:,1)=kondmax 
     278             kond(:,ny)=kondmax 
     279             vieuxhwater(:,:) = hwater(:,:) 
     280             !$OMP END WORKSHARE 
     281             !$OMP END PARALLEL 
    276282!!$OMP SINGLE 
    277            call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,hwater) 
     283             call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,hwater) 
    278284!!$OMP END SINGLE 
    279         else 
    280            !         print*,'dt=',dt,'pas de relax_water' 
    281         endif!!!!!!!!!!!!!!!!!! fin relax_water si dt>0 
    282  
    283      endif 
    284  
    285  
    286  
    287  
    288      ! Apres relaxation, boundary conditions, extremum values 
    289      ! ================, ===================, =============== 
    290  
    291      !   ------------variation d'epaisseur entre 2 pas de temps ------------ 
    292  
    293      ! on le fait avant les butoirs pour justement pouvoir les estimer 
    294 !$OMP PARALLEL 
    295      if (dt.gt.0.) then 
    296         !         print*,'dt=',dt,'pas de relax_water' 
    297         !$OMP DO 
    298         do j=1,ny 
    299            do i=1,nx 
    300               hdotwater(i,j)=(hwater(i,j)-vieuxhwater(i,j))/dt  
    301            end do 
    302         end do 
    303         !$OMP END DO 
    304      endif 
    305  
    306 !$OMP DO PRIVATE(Zflot) 
    307      do i=1,nx 
    308         do j=1,ny 
    309  
    310  
    311            !  ______________ valeurs de hwater et pwater _____________________________ 
    312            !  
    313            if (flot(I,J).or.H(I,J).le.1.5) then ! we are outside of the ice sheet 
    314  
    315               if (flot(i,j)) then  ! if flot > hwater=hwater in ocean  
    316                  hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j) 
    317                  !             hwater(i,j)= max(0.,hwater(i,j)) 
    318                  if (hwater(i,j).lt.0.) hwater(i,j)=0. 
    319                  pwater(i,j)= hwater(i,j)*rowg 
    320               else 
    321                  hwater(i,j) = 0. ! bare grounded land > no hwater 
    322                  pwater(i,j)=0. 
    323               endif 
    324  
    325            else          !            sous la  calotte ---------------------- 
    326  
    327               ! Attention le bloc suivant est pour le run 20  
    328               !           Zflot=row/ro*(sealevel-Bsoc(i,j))-10. 
    329               Zflot=H(i,j)*rog/rofreshg 
    330  
    331               if (hwater(i,j).le.0.) then  
    332                  hwater(i,j)=0. 
    333  
    334               else if (hwater(i,j).gt.zflot) then 
    335                  hwater(i,j)=zflot 
    336                  hw(i,j)=min(hwater(i,j),hmax_wat) 
    337                  pwater(i,j)=rofreshg*Hwater(i,j) 
    338               endif 
    339               !            if ((i.eq.60).and.(j.eq.60)) write(6,*) hw(i,j),hwater(i,j) 
    340  
    341               ! sinon 
    342  
    343               !          hw(i,j)=min(hw(i,j),hmax_wat)   
    344               !          Hwater(i,j)=hw(i,j) 
    345               !          pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5) 
    346            endif 
    347  
    348            ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression 
    349            ! ----------------------------------------------------------------------------- 
    350            !            if (HWATER(i,j).gt.poro_till*hmax_till) then 
    351            !             pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till) 
    352            !            endif 
    353         end do 
    354      end do 
    355 !$OMP END DO 
    356  
    357      !   ************ valeurs de HWATER pour les coins ******** 
    358  
    359      hwater(1,1)=(hwater(1,2)+hwater(2,1))/2. 
    360      hwater(1,ny)=(hwater(1,ny-1)+hwater(2,ny))/2. 
    361      hwater(nx,1)=(hwater(nx,2)+hwater(nx-1,1))/2. 
    362      hwater(nx,ny)=(hwater(nx,ny-1)+hwater(nx-1,ny))/2. 
    363  
    364      ! pour les sorties de flux d'eau 
    365      !$OMP DO 
    366      do j=2,ny 
    367         do i=2,nx 
    368            if  (keff(i,j)==0..or.keff(i-1,j)==0.) then 
    369               phiwx(i,j)=0. ! to avoid division by o              
    370            else 
    371               phiwx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) 
    372  
    373               phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j)) 
    374            endif 
    375            ! ligne pour sortir les pgx 
    376            pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) 
    377         end do 
    378      end do 
    379      !$OMP END DO 
    380      !$OMP DO 
    381      do j=2,ny 
    382         do i=2,nx 
    383            if  (keff(i,j)==0..or.keff(i,j-1)==0.) then 
    384               phiWy(i,j)=0. ! to avoid division by o 
    385            else 
    386               phiWy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) 
    387               phiWy(i,j)=phiWy(i,j)*2*(keff(i,j)*keff(i,j-1))/(keff(i,j)+keff(i,j-1)) 
    388            endif 
    389            pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) 
    390         enddo 
    391      enddo 
    392      !$OMP END DO 
    393      !$OMP END PARALLEL 
    394  
    395   else  ! partie originellement dans icetemp à mettre dans un autre module 
    396      ! (système module choix)     hauteur d'eau cumulee 
    397  
    398      if (ISYNCHRO.eq.1) then 
    399         !$OMP PARALLEL 
    400         !$OMP DO 
    401         do j=1,ny 
    402            do i=1,nx 
    403               if (.not.flot(i,j)) then 
    404                  hwater(i,j)=hwater(i,j)+(bmelt(i,j)*dtt) 
    405                  hwater(i,j)=hwater(i,j)-(infiltr*dtt) 
    406                  hwater(i,j)=max(hwater(i,j),0.) 
    407                  hwater(i,j)=min(hwater(i,j),hwatermax) 
    408               else 
    409                  hwater(i,j)=hwatermax 
    410               endif 
    411               !        if (IBASE(I,J).eq.1) HWATER(I,J)=0. 
    412            end do 
    413         end do 
    414         !$OMP END DO 
    415         !$OMP END PARALLEL 
    416      end if 
    417       
    418   endif ecoul 
    419  
    420 !~   ! Temps elapsed final 
    421 !~   call system_clock(count=t2, count_rate=ir) 
    422 !~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
    423 !~   ! Temps CPU de calcul final 
    424 !~   call cpu_time(t_cpu_1) 
    425 !~   t_cpu = t_cpu_1 - t_cpu_0 
    426  
    427 !~   ! Impression du resultat. 
    428 !~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
    429 !~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
    430 !~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
    431 !~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
    432 !~            nx,ny,temps,t_cpu,norme 
    433  
    434   return 
    435 end subroutine eaubasale 
    436  end module eau_basale 
     285          else 
     286             !         print*,'dt=',dt,'pas de relax_water' 
     287          endif!!!!!!!!!!!!!!!!!! fin relax_water si dt>0 
     288 
     289       endif 
     290 
     291 
     292 
     293 
     294       ! Apres relaxation, boundary conditions, extremum values 
     295       ! ================, ===================, =============== 
     296 
     297       !   ------------variation d'epaisseur entre 2 pas de temps ------------ 
     298 
     299       ! on le fait avant les butoirs pour justement pouvoir les estimer 
     300       !$OMP PARALLEL 
     301       if (dt.gt.0.) then 
     302          !         print*,'dt=',dt,'pas de relax_water' 
     303          !$OMP DO 
     304          do j=1,ny 
     305             do i=1,nx 
     306                hdotwater(i,j)=(hwater(i,j)-vieuxhwater(i,j))/dt  
     307             end do 
     308          end do 
     309          !$OMP END DO 
     310       endif 
     311 
     312       !$OMP DO PRIVATE(Zflot) 
     313       do i=1,nx 
     314          do j=1,ny 
     315 
     316 
     317             !  ______________ valeurs de hwater et pwater _____________________________ 
     318             !  
     319             if (flot(I,J).or.H(I,J).le.1.5) then ! we are outside of the ice sheet 
     320 
     321                if (flot(i,j)) then  ! if flot > hwater=hwater in ocean  
     322                   hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j) 
     323                   !             hwater(i,j)= max(0.,hwater(i,j)) 
     324                   if (hwater(i,j).lt.0.) hwater(i,j)=0. 
     325                   pwater(i,j)= hwater(i,j)*rowg 
     326                else 
     327                   hwater(i,j) = 0. ! bare grounded land > no hwater 
     328                   pwater(i,j)=0. 
     329                endif 
     330 
     331             else          !            sous la  calotte ---------------------- 
     332 
     333                ! Attention le bloc suivant est pour le run 20  
     334                !           Zflot=row/ro*(sealevel-Bsoc(i,j))-10. 
     335                Zflot=H(i,j)*rog/rofreshg 
     336 
     337                if (hwater(i,j).le.0.) then  
     338                   hwater(i,j)=0. 
     339 
     340                else if (hwater(i,j).gt.zflot) then 
     341                   hwater(i,j)=zflot 
     342                   hw(i,j)=min(hwater(i,j),hmax_wat) 
     343                   pwater(i,j)=rofreshg*Hwater(i,j) 
     344                endif 
     345                !            if ((i.eq.60).and.(j.eq.60)) write(6,*) hw(i,j),hwater(i,j) 
     346 
     347                ! sinon 
     348 
     349                !          hw(i,j)=min(hw(i,j),hmax_wat)   
     350                !          Hwater(i,j)=hw(i,j) 
     351                !          pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5) 
     352             endif 
     353 
     354             ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression 
     355             ! ----------------------------------------------------------------------------- 
     356             !            if (HWATER(i,j).gt.poro_till*hmax_till) then 
     357             !             pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till) 
     358             !            endif 
     359          end do 
     360       end do 
     361       !$OMP END DO 
     362 
     363       !   ************ valeurs de HWATER pour les coins ******** 
     364 
     365       hwater(1,1)=(hwater(1,2)+hwater(2,1))/2. 
     366       hwater(1,ny)=(hwater(1,ny-1)+hwater(2,ny))/2. 
     367       hwater(nx,1)=(hwater(nx,2)+hwater(nx-1,1))/2. 
     368       hwater(nx,ny)=(hwater(nx,ny-1)+hwater(nx-1,ny))/2. 
     369 
     370       ! pour les sorties de flux d'eau 
     371       !$OMP DO 
     372       do j=2,ny 
     373          do i=2,nx 
     374             if  (keff(i,j)==0..or.keff(i-1,j)==0.) then 
     375                phiwx(i,j)=0. ! to avoid division by o              
     376             else 
     377                phiwx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) 
     378 
     379                phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j)) 
     380             endif 
     381             ! ligne pour sortir les pgx 
     382             pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) 
     383          end do 
     384       end do 
     385       !$OMP END DO 
     386       !$OMP DO 
     387       do j=2,ny 
     388          do i=2,nx 
     389             if  (keff(i,j)==0..or.keff(i,j-1)==0.) then 
     390                phiWy(i,j)=0. ! to avoid division by o 
     391             else 
     392                phiWy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) 
     393                phiWy(i,j)=phiWy(i,j)*2*(keff(i,j)*keff(i,j-1))/(keff(i,j)+keff(i,j-1)) 
     394             endif 
     395             pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) 
     396          enddo 
     397       enddo 
     398       !$OMP END DO 
     399       !$OMP END PARALLEL 
     400 
     401    else  ! partie originellement dans icetemp à mettre dans un autre module 
     402       ! (système module choix)     hauteur d'eau cumulee 
     403 
     404       if (ISYNCHRO.eq.1) then 
     405          !$OMP PARALLEL 
     406          !$OMP DO 
     407          do j=1,ny 
     408             do i=1,nx 
     409                if (.not.flot(i,j)) then 
     410                   hwater(i,j)=hwater(i,j)+(bmelt(i,j)*dtt) 
     411                   hwater(i,j)=hwater(i,j)-(infiltr*dtt) 
     412                   hwater(i,j)=max(hwater(i,j),0.) 
     413                   hwater(i,j)=min(hwater(i,j),hwatermax) 
     414                else 
     415                   hwater(i,j)=hwatermax 
     416                endif 
     417                !        if (IBASE(I,J).eq.1) HWATER(I,J)=0. 
     418             end do 
     419          end do 
     420          !$OMP END DO 
     421          !$OMP END PARALLEL 
     422       end if 
     423 
     424    endif ecoul 
     425 
     426    !~   ! Temps elapsed final 
     427    !~   call system_clock(count=t2, count_rate=ir) 
     428    !~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
     429    !~   ! Temps CPU de calcul final 
     430    !~   call cpu_time(t_cpu_1) 
     431    !~   t_cpu = t_cpu_1 - t_cpu_0 
     432 
     433    !~   ! Impression du resultat. 
     434    !~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
     435    !~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
     436    !~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
     437    !~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
     438    !~            nx,ny,temps,t_cpu,norme 
     439 
     440    return 
     441  end subroutine eaubasale 
     442end module eau_basale 
Note: See TracChangeset for help on using the changeset viewer.