| 1 | !> \file conserv-mass-adv-diff_mod.F90 |
|---|
| 2 | !! Module qui calcule l evolution de l epaisseur en resolvant |
|---|
| 3 | !! une equation d'advection diffusion |
|---|
| 4 | !< |
|---|
| 5 | |
|---|
| 6 | !> \namespace equat_adv_diff_2D |
|---|
| 7 | !! Module qui calcule l evolution de l epaisseur en resolvant |
|---|
| 8 | !! une equation d'advection diffusion |
|---|
| 9 | !! \author Cat |
|---|
| 10 | !! \date Fevrier 2007 |
|---|
| 11 | !! @note Used modules |
|---|
| 12 | !! @note - use module3D_phy |
|---|
| 13 | !! @note - use reso_adv_diff_2D |
|---|
| 14 | !< |
|---|
| 15 | module equat_adv_diff_2D ! Cat Fevrier 2007 |
|---|
| 16 | |
|---|
| 17 | use module3D_phy |
|---|
| 18 | use reso_adv_diff_2D |
|---|
| 19 | |
|---|
| 20 | |
|---|
| 21 | implicit none |
|---|
| 22 | |
|---|
| 23 | real, dimension(nx,ny) :: advmx !< partie advective en x |
|---|
| 24 | real, dimension(nx,ny) :: advmy !< partie advective en y |
|---|
| 25 | real, dimension(nx,ny) :: dmx !< partie advective en x |
|---|
| 26 | real, dimension(nx,ny) :: dmy !< partie advective en y |
|---|
| 27 | real, dimension(nx,ny) :: VieuxH !< ancienne valeur de H |
|---|
| 28 | logical, dimension(nx,ny) :: zonemx !< pour separer advection-diffusion |
|---|
| 29 | logical, dimension(nx,ny) :: zonemy !< pour separer advection-diffusion |
|---|
| 30 | real :: adv_frac !< fraction du flux traitee en advection |
|---|
| 31 | integer :: itesti |
|---|
| 32 | integer :: itour |
|---|
| 33 | |
|---|
| 34 | contains |
|---|
| 35 | !> SUBROUTINE: init_icethick |
|---|
| 36 | !! Routine qui permet d'initialiser les parametres icethick |
|---|
| 37 | !! @note Definition du parametre qui gerent la repartition advection diffusion |
|---|
| 38 | !> |
|---|
| 39 | subroutine init_icethick |
|---|
| 40 | |
|---|
| 41 | |
|---|
| 42 | write(num_rep_42,*)'Conservation de la masse avec equation advection-diffusion ' |
|---|
| 43 | write(num_rep_42,*)'la repartition depend de adv_frac' |
|---|
| 44 | write(num_rep_42,*)'>1 -> advection seule' |
|---|
| 45 | write(num_rep_42,*)'0 -> diffusion seule' |
|---|
| 46 | write(num_rep_42,*)'0<*<1 -> fraction de l advection' |
|---|
| 47 | write(num_rep_42,*) '-1 -> zones diffusion + zones advecttion' |
|---|
| 48 | |
|---|
| 49 | !varname='adv_frac' |
|---|
| 50 | |
|---|
| 51 | adv_frac=2. ! 0.8 |
|---|
| 52 | |
|---|
| 53 | !adv_frac=0 |
|---|
| 54 | write(num_rep_42,*)'adv_frac=',adv_frac |
|---|
| 55 | |
|---|
| 56 | call init_reso_adv_diff_2D |
|---|
| 57 | |
|---|
| 58 | Dx1=1./Dx |
|---|
| 59 | |
|---|
| 60 | itour=0 |
|---|
| 61 | |
|---|
| 62 | end subroutine init_icethick |
|---|
| 63 | |
|---|
| 64 | !------------------------------------------------------------------ |
|---|
| 65 | !> SUBROUTINE: icethick3 |
|---|
| 66 | !! Routine de calcul de l'epaisseur de glace |
|---|
| 67 | !> |
|---|
| 68 | subroutine icethick3 |
|---|
| 69 | |
|---|
| 70 | implicit none |
|---|
| 71 | real,dimension(nx,ny) :: Dminx,Dminy |
|---|
| 72 | real,dimension(nx,ny) :: Uxdiff,Uydiff ! vitesse due a la diffusion |
|---|
| 73 | integer :: it1,it2,it3 |
|---|
| 74 | |
|---|
| 75 | real aux ! pour le calcul du critere |
|---|
| 76 | real maxdia ! sur le pas de temps |
|---|
| 77 | integer imaxdia,jmaxdia |
|---|
| 78 | |
|---|
| 79 | if (itracebug.eq.1) call tracebug(' Entree dans routine icethick') |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | ! Copie de H sur vieuxH |
|---|
| 83 | ! -------------------- |
|---|
| 84 | where (mk0(:,:).eq.0) |
|---|
| 85 | vieuxH(:,:)=0. |
|---|
| 86 | elsewhere |
|---|
| 87 | vieuxH(:,:)=H(:,:) |
|---|
| 88 | end where |
|---|
| 89 | |
|---|
| 90 | ! mk0 est le masque à ne jamais dépasser |
|---|
| 91 | ! mk0=0 -> pas de glace |
|---|
| 92 | ! mk0=-1 -> on garde la même epaisseur |
|---|
| 93 | |
|---|
| 94 | Maxdia = -1.0 |
|---|
| 95 | imaxdia=0 |
|---|
| 96 | jmaxdia=0 |
|---|
| 97 | |
|---|
| 98 | ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante |
|---|
| 99 | do i=2,nx-1 |
|---|
| 100 | do j=2,ny-1 |
|---|
| 101 | aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) & |
|---|
| 102 | + (abs(uybar(i,j)) + abs(uybar(i,j+1))) |
|---|
| 103 | aux=aux*dx1*0.5 |
|---|
| 104 | if (aux.gt.maxdia) then |
|---|
| 105 | maxdia=aux |
|---|
| 106 | imaxdia=i |
|---|
| 107 | jmaxdia=j |
|---|
| 108 | endif |
|---|
| 109 | end do |
|---|
| 110 | end do |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | if (maxdia.ge.(testdiag/dtmax)) then |
|---|
| 114 | dt = testdiag/Maxdia |
|---|
| 115 | dt=max(dt,dtmin) |
|---|
| 116 | else |
|---|
| 117 | dt = dtmax |
|---|
| 118 | end if |
|---|
| 119 | |
|---|
| 120 | ! next_time ajuste le pas de temps pour faciliter la synchronisation |
|---|
| 121 | ! avec le pas de temps long (dtt) |
|---|
| 122 | |
|---|
| 123 | call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) |
|---|
| 124 | |
|---|
| 125 | |
|---|
| 126 | !!$write(166,*)'--------------------------------------------------------' |
|---|
| 127 | !!$write(166,*)'apres old calcul de dt',time,dt |
|---|
| 128 | |
|---|
| 129 | ! calcul diffusivite - vitesse |
|---|
| 130 | !----------------------------------------------------------------- |
|---|
| 131 | |
|---|
| 132 | Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:)) ! vitesse qui peut s'exprimer en terme diffusif |
|---|
| 133 | Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:)) ! dans les where qui suivent elle est limitee a uxbar |
|---|
| 134 | ! pour eviter des problemes dans le cas 0< adv_frac <1 |
|---|
| 135 | dmx(:,:)=diffmx(:,:) |
|---|
| 136 | dmy(:,:)=diffmy(:,:) |
|---|
| 137 | |
|---|
| 138 | ! schema amont pour la diffusion en x |
|---|
| 139 | where (Uxbar(:,:).ge.0.) |
|---|
| 140 | dmx(:,:)=diffmx(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=1) |
|---|
| 141 | uxdiff(:,:)=min(uxdiff(:,:),uxbar(:,:)) ! pour tenir compte limitation utot |
|---|
| 142 | elsewhere |
|---|
| 143 | dmx(:,:)=diffmx(:,:)*H(:,:) |
|---|
| 144 | uxdiff(:,:)=max(uxdiff(:,:),uxbar(:,:)) ! pour tenir compte limitation utot |
|---|
| 145 | end where |
|---|
| 146 | |
|---|
| 147 | |
|---|
| 148 | |
|---|
| 149 | ! schema amont pour la diffusion en y |
|---|
| 150 | where (uybar(:,:).ge.0.) |
|---|
| 151 | dmy(:,:)=diffmy(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=2) |
|---|
| 152 | uydiff(:,:)=min(uydiff(:,:),uybar(:,:)) ! pour tenir compte limitation utot |
|---|
| 153 | elsewhere |
|---|
| 154 | dmy(:,:)=dmy(:,:)*H(:,:) |
|---|
| 155 | uydiff(:,:)=max(uydiff(:,:),uybar(:,:)) ! pour tenir compte limitation utot |
|---|
| 156 | end where |
|---|
| 157 | |
|---|
| 158 | |
|---|
| 159 | |
|---|
| 160 | ! tests pour la répartition advection - diffusion |
|---|
| 161 | !------------------------------------------------ |
|---|
| 162 | |
|---|
| 163 | if (adv_frac.gt.1.) then ! advection seulement |
|---|
| 164 | advmx(:,:)=Uxbar(:,:) |
|---|
| 165 | advmy(:,:)=Uybar(:,:) |
|---|
| 166 | Dmx(:,:)=0. |
|---|
| 167 | Dmy(:,:)=0. |
|---|
| 168 | |
|---|
| 169 | else if (abs(adv_frac).lt.1.e-8) then ! diffusion seulement |
|---|
| 170 | advmx(:,:)=0. |
|---|
| 171 | advmy(:,:)=0. |
|---|
| 172 | |
|---|
| 173 | ! D reste la valeur au dessus) |
|---|
| 174 | |
|---|
| 175 | else if ((adv_frac.ge.1.e-8).and.(adv_frac.le.1.)) then ! advection-diffusion) |
|---|
| 176 | |
|---|
| 177 | ! selon x -------------------------------- |
|---|
| 178 | |
|---|
| 179 | ! advection = adv_frac* tout ce qui n'est pas diffusion |
|---|
| 180 | ! diffusion = ce qui peut être exprime en diffusion |
|---|
| 181 | ! + une partie (1-adv_frac) de l'advection |
|---|
| 182 | |
|---|
| 183 | where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.)) ! cas general |
|---|
| 184 | advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:)) ! tout ce qui n'est pas diffusion |
|---|
| 185 | advmx(:,:)=advmx(:,:)*adv_frac ! la fraction adv_frac du precedent |
|---|
| 186 | Dminx(:,:)=-(Uxbar(:,:)-advmx(:,:))/sdx(:,:) ! ce qui reste exprime en diffusivite |
|---|
| 187 | ! a multiplier par H |
|---|
| 188 | |
|---|
| 189 | ! schema amont pour la diffusivite |
|---|
| 190 | where (uxbar(:,:).ge.0.) |
|---|
| 191 | dmx(:,:)=dminx(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=1) |
|---|
| 192 | elsewhere |
|---|
| 193 | dmx(:,:)=dminx(:,:)*H(:,:) |
|---|
| 194 | end where |
|---|
| 195 | |
|---|
| 196 | |
|---|
| 197 | elsewhere ! zones trop plates ou trop fines |
|---|
| 198 | Dmx(:,:)=0. |
|---|
| 199 | advmx(:,:)=Uxbar(:,:) |
|---|
| 200 | end where |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | ! selon y -------------------------------- |
|---|
| 204 | |
|---|
| 205 | ! advection = adv_frac* tout ce qui n'est pas diffusion |
|---|
| 206 | ! diffusion = ce qui peut être exprime en diffusion |
|---|
| 207 | ! + une partie (1-adv_frac) de l'advection |
|---|
| 208 | |
|---|
| 209 | where ((abs(sdy(:,:)).gt.1.e-8).and.(Hmy(:,:).gt.2.)) ! cas general |
|---|
| 210 | advmy(:,:)=(Uybar(:,:)-Uydiff(:,:)) ! tout ce qui n'est pas diffusion |
|---|
| 211 | advmy(:,:)=advmy(:,:)*adv_frac ! la fraction adv_frac du precedent |
|---|
| 212 | Dminy(:,:)=-(Uybar(:,:)-advmy(:,:))/sdy(:,:) ! ce qui reste exprime en diffusivite |
|---|
| 213 | ! a multiplier par H |
|---|
| 214 | |
|---|
| 215 | ! schema amont pour la diffusivite |
|---|
| 216 | where (uybar(:,:).ge.0.) |
|---|
| 217 | dmy(:,:)=dminy(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=2) |
|---|
| 218 | elsewhere |
|---|
| 219 | dmy(:,:)=dminy(:,:)*H(:,:) |
|---|
| 220 | end where |
|---|
| 221 | |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | elsewhere ! zones trop plates ou trop fines |
|---|
| 225 | Dmy(:,:)=0. |
|---|
| 226 | advmy(:,:)=Uybar(:,:) |
|---|
| 227 | end where |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | else if (adv_frac.lt.0) then ! diffusif dans zones SIA, advectif ailleurs |
|---|
| 231 | |
|---|
| 232 | zonemx(:,:)=flgzmx(:,:) |
|---|
| 233 | zonemy(:,:)=flgzmy(:,:) |
|---|
| 234 | |
|---|
| 235 | ! selon x -------------- |
|---|
| 236 | where (zonemx(:,:)) |
|---|
| 237 | advmx(:,:)=Uxbar(:,:) |
|---|
| 238 | Dmx(:,:)=0. |
|---|
| 239 | elsewhere |
|---|
| 240 | advmx(:,:)=0. |
|---|
| 241 | end where |
|---|
| 242 | |
|---|
| 243 | ! selon y -------------- |
|---|
| 244 | where (zonemy(:,:)) |
|---|
| 245 | advmy(:,:)=Uybar(:,:) |
|---|
| 246 | Dmy(:,:)=0. |
|---|
| 247 | elsewhere |
|---|
| 248 | advmy(:,:)=0. |
|---|
| 249 | end where |
|---|
| 250 | |
|---|
| 251 | end if |
|---|
| 252 | |
|---|
| 253 | |
|---|
| 254 | where (flot(:,:) ) |
|---|
| 255 | tabtest(:,:)=1. |
|---|
| 256 | elsewhere |
|---|
| 257 | tabtest(:,:)=0. |
|---|
| 258 | end where |
|---|
| 259 | |
|---|
| 260 | !------------------------------------------------------------------detect |
|---|
| 261 | !!$ tabtest(:,:)=dmx(:,:) |
|---|
| 262 | !!$ |
|---|
| 263 | !!$call detect_assym(nx,ny,0,41,1,0,1,0,tabtest,itesti) |
|---|
| 264 | !!$ |
|---|
| 265 | !!$if (itesti.gt.0) then |
|---|
| 266 | !!$ write(6,*) 'asymetrie sur dmx avant resol pour time=',time,itesti |
|---|
| 267 | !!$ stop |
|---|
| 268 | !!$else |
|---|
| 269 | !!$ write(6,*) ' pas d asymetrie sur dmx avant resol pour time=',time,itesti |
|---|
| 270 | !!$end if |
|---|
| 271 | !----------------------------------------------------------- fin detect |
|---|
| 272 | |
|---|
| 273 | |
|---|
| 274 | !------------------------------------------------------------------detect |
|---|
| 275 | !!$ tabtest(:,:)=dmy(:,:) |
|---|
| 276 | !!$ |
|---|
| 277 | !!$call detect_assym(nx,ny,0,41,1,0,1,1,tabtest,itesti) |
|---|
| 278 | !!$ |
|---|
| 279 | !!$if (itesti.gt.0) then |
|---|
| 280 | !!$ write(6,*) 'asymetrie sur dmy avant resol pour time=',time,itesti |
|---|
| 281 | !!$ stop |
|---|
| 282 | !!$else |
|---|
| 283 | !!$ write(6,*) ' pas d asymetrie sur dmy avant resol pour time=',time,itesti |
|---|
| 284 | !!$end if |
|---|
| 285 | !----------------------------------------------------------- fin detect |
|---|
| 286 | |
|---|
| 287 | !------------------------------------------------------------------detect |
|---|
| 288 | !!$ tabtest(:,:)=advmx(:,:) |
|---|
| 289 | !!$ |
|---|
| 290 | !!$call detect_assym(nx,ny,0,41,1,0,1,0,tabtest,itesti) |
|---|
| 291 | !!$ |
|---|
| 292 | !!$if (itesti.gt.0) then |
|---|
| 293 | !!$ write(6,*) 'asymetrie sur advmx avant resol pour time=',time,itesti |
|---|
| 294 | !!$ stop |
|---|
| 295 | !!$else |
|---|
| 296 | !!$ write(6,*) ' pas d asymetrie sur advdmx avant resol pour time=',time,itesti |
|---|
| 297 | !!$end if |
|---|
| 298 | !----------------------------------------------------------- fin detect |
|---|
| 299 | |
|---|
| 300 | !------------------------------------------------------------------detect |
|---|
| 301 | !!$ tabtest(:,:)=advmy(:,:) |
|---|
| 302 | !!$ |
|---|
| 303 | !!$call detect_assym(nx,ny,0,41,1,0,-1,1,tabtest,itesti) |
|---|
| 304 | !!$ |
|---|
| 305 | !!$if (itesti.gt.0) then |
|---|
| 306 | !!$ write(6,*) 'asymetrie sur advmy avant resol pour time=',time,itesti |
|---|
| 307 | !!$ stop |
|---|
| 308 | !!$else |
|---|
| 309 | !!$ write(6,*) ' pas d asymetrie sur advmy avant resol pour time=',time,itesti |
|---|
| 310 | !!$end if |
|---|
| 311 | !----------------------------------------------------------- fin detect |
|---|
| 312 | |
|---|
| 313 | ! nouveau calcul de dt |
|---|
| 314 | aux=maxval( (abs(advmx(:,:))+abs(advmy(:,:)))/dx+(abs(dmx(:,:))+abs(dmy(:,:)))/dx/dx) |
|---|
| 315 | |
|---|
| 316 | ! write(166,*) 'critere pour dt',time-dt,testdiag/aux,dt |
|---|
| 317 | |
|---|
| 318 | if (aux.gt.1.e-20) then |
|---|
| 319 | if (testdiag/aux.lt.dt) then |
|---|
| 320 | time=time-dt |
|---|
| 321 | dt=testdiag/aux*4. |
|---|
| 322 | call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) |
|---|
| 323 | end if |
|---|
| 324 | end if |
|---|
| 325 | |
|---|
| 326 | |
|---|
| 327 | ! print*,'Iteration =',nt,' temps = ',time-dt,' dt = ',dt |
|---|
| 328 | ! vient de step |
|---|
| 329 | !!$ if (time.lt.timemax-100) then |
|---|
| 330 | !!$ print*,'le modele remonte le temps de',timemax,'a',time |
|---|
| 331 | !!$ stop |
|---|
| 332 | !!$ endif |
|---|
| 333 | timemax=time |
|---|
| 334 | |
|---|
| 335 | if (time.lt.TGROUNDED) then |
|---|
| 336 | MARINE=.false. |
|---|
| 337 | else |
|---|
| 338 | MARINE=.true. |
|---|
| 339 | endif |
|---|
| 340 | |
|---|
| 341 | ! fin vient de step |
|---|
| 342 | |
|---|
| 343 | ! les variables dtdx et dtdx2 sont globales |
|---|
| 344 | Dtdx2=Dt/(Dx**2) |
|---|
| 345 | dtdx=dt/dx |
|---|
| 346 | |
|---|
| 347 | !!$write(166,*) 'time=',time |
|---|
| 348 | !!$ |
|---|
| 349 | !!$ |
|---|
| 350 | !!$write(166,*) 'advmx,uxbar' |
|---|
| 351 | !!$do j=43,45 |
|---|
| 352 | !!$ write(166,*) advmx(118:122,j) |
|---|
| 353 | !!$end do |
|---|
| 354 | !!$ |
|---|
| 355 | !!$write(166,*) 'advmy,uybar' |
|---|
| 356 | !!$do j=43,45 |
|---|
| 357 | !!$ write(166,*) advmy(118:122,j) |
|---|
| 358 | !!$end do |
|---|
| 359 | !!$ |
|---|
| 360 | !!$write(166,*) 'vieuxH,hdot' |
|---|
| 361 | !!$do j=43,45 |
|---|
| 362 | !!$ write(166,*) vieuxH(118:122,j) |
|---|
| 363 | !!$ write(166,*) Hdot(118:122,j) |
|---|
| 364 | !!$ |
|---|
| 365 | !!$end do |
|---|
| 366 | |
|---|
| 367 | |
|---|
| 368 | |
|---|
| 369 | call resol_adv_diff_2D (Dmx,Dmy,advmx,advmy,vieuxH) |
|---|
| 370 | |
|---|
| 371 | |
|---|
| 372 | |
|---|
| 373 | ! post traitements |
|---|
| 374 | !-------------------- |
|---|
| 375 | |
|---|
| 376 | ! epaisseur nulle sur les bords |
|---|
| 377 | H(1,:)=0. |
|---|
| 378 | H(nx,:)=0. |
|---|
| 379 | H(:,1)=0. |
|---|
| 380 | H(:,ny)=0. |
|---|
| 381 | |
|---|
| 382 | ! attention le bloc suivant sert a garder les ice-shelves stationnaires |
|---|
| 383 | |
|---|
| 384 | if (igrdline.eq.1) then |
|---|
| 385 | where (flot(:,:)) |
|---|
| 386 | Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt |
|---|
| 387 | H(:,:)=vieuxH(:,:) |
|---|
| 388 | end where |
|---|
| 389 | end if |
|---|
| 390 | |
|---|
| 391 | ! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord |
|---|
| 392 | where (H(:,:).lt.0.) |
|---|
| 393 | ablbord(:,:)=H(:,:)/dt |
|---|
| 394 | elsewhere |
|---|
| 395 | ablbord(:,:)=0. |
|---|
| 396 | endwhere |
|---|
| 397 | |
|---|
| 398 | H(:,:)=max(H(:,:),0.) ! pas d'epaisseur negative |
|---|
| 399 | |
|---|
| 400 | |
|---|
| 401 | ! calcul du masque "ice" |
|---|
| 402 | where (flot(:,:)) ! points flottants, sera éventuellement réévalué dans flottab |
|---|
| 403 | H(:,:)=max(H(:,:),1.) ! dans la partie marine l'épaisseur mini est 1 m |
|---|
| 404 | where(H(:,:).gt.1.) |
|---|
| 405 | ice(:,:)=1 |
|---|
| 406 | elsewhere |
|---|
| 407 | ice(:,:)=0 |
|---|
| 408 | endwhere |
|---|
| 409 | elsewhere ! points posés |
|---|
| 410 | where(H(:,:).gt.0.) |
|---|
| 411 | ice(:,:)=1 |
|---|
| 412 | elsewhere |
|---|
| 413 | ice(:,:)=0 |
|---|
| 414 | endwhere |
|---|
| 415 | endwhere |
|---|
| 416 | |
|---|
| 417 | ! calcul de hdot. Quand igrdline=1, on garde hdot pour l'estimation du bmelt d'equilibre |
|---|
| 418 | |
|---|
| 419 | if (igrdline.eq.1) then |
|---|
| 420 | where (.not.flot(:,:)) |
|---|
| 421 | Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt |
|---|
| 422 | endwhere |
|---|
| 423 | else |
|---|
| 424 | Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt |
|---|
| 425 | endif |
|---|
| 426 | |
|---|
| 427 | |
|---|
| 428 | ! si mk0=-1, on garde l'epaisseur precedente |
|---|
| 429 | ! en attendant un masque plus propre |
|---|
| 430 | |
|---|
| 431 | where(mk0(:,:).eq.-1) |
|---|
| 432 | H(:,:)=vieuxH(:,:) |
|---|
| 433 | Hdot(:,:)=0. |
|---|
| 434 | end where |
|---|
| 435 | |
|---|
| 436 | |
|---|
| 437 | |
|---|
| 438 | !calul de l'ablation sur les bords (pourrait n'être appelé que pour les sorties courtes) |
|---|
| 439 | if (isynchro.eq.1) call ablation_bord |
|---|
| 440 | |
|---|
| 441 | |
|---|
| 442 | |
|---|
| 443 | if (itracebug.eq.1) call tracebug(' Fin routine icethick')!, maxval(H) ',maxval(H(:,:)) |
|---|
| 444 | end subroutine icethick3 |
|---|
| 445 | |
|---|
| 446 | end module equat_adv_diff_2D |
|---|
![(please configure the [header_logo] section in trac.ini)](/grisli/chrome/site/Laurentide-16km-3D-Umag-water-bordure-GRISLI-small.png)
