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