![(please configure the [header_logo] section in trac.ini)](/grisli/chrome/site/Laurentide-16km-3D-Umag-water-bordure-GRISLI-small.png)
Changeset 396
- Timestamp:
- 03/24/23 15:10:41 (13 months ago)
- File:
-
- 1 edited
-
branches/GRISLIv3/SOURCES/flottab2-0.7.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/GRISLIv3/SOURCES/flottab2-0.7.f90
r332 r396 15 15 !< 16 16 module flottab_mod 17 17 18 18 !$ USE OMP_LIB 19 USE module3D_phy 20 use module_choix 21 22 23 24 IMPLICIT NONE 25 26 real :: surnet !< surnet hauteur de glace au dessus de la mer 27 real :: archim !< test de flottaison 28 ! real, parameter :: Hmin=1.001 !< Hmin pour être considere comme point ice 29 30 integer:: itestf 31 32 logical,dimension(nx,ny) :: gz1mx,gz1my 33 logical,dimension(nx,ny) :: fl1mx,fl1my 34 35 36 real,dimension(nx,ny) :: uxs1 !< uxbar a l'entree de flottab 37 real,dimension(nx,ny) :: uys1 !< uybar a l'entree de flottab 38 39 40 integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y 41 42 43 ! Variables pour la determination des differents shelfs/stream 44 ! (representés comme des taches ou l'on resoud l'eq elliptique) 45 !________________________________________________________________ 46 integer,parameter :: n_ta_max=2000!< nombre de tache max 47 integer,dimension(nx,ny) :: table_out !< pour les numeros des taches 48 integer,dimension(nx,ny) :: tablebis !< pour les numeros des taches 49 integer,dimension(0:n_ta_max) :: compt !< contient les equivalence entre les taches 50 integer,dimension(0:n_ta_max) :: nb_pts_tache !< indique le nombre de points par tache 51 logical,dimension(0:n_ta_max) :: iceberg1D !< T si iceberg, F si calotte posee 52 53 logical,dimension(nx,ny) :: mask_tache_ij !< masque de travail 54 !< vrai pour toute la tache de i,j 55 integer,dimension(2) :: smax_coord !< pour le maxloc des iles 56 57 ! Variables pour determiner le point le plus haut (surf) 58 ! d'une ile completement stream 59 !_________________________________________________________ 60 61 ! icetrim : T si ice stream, F si calotte posee(vertical shear) 62 logical,dimension(n_ta_max) :: icetrim 63 64 integer :: ii,jj 65 integer :: smax_i 66 integer :: smax_j 67 real :: smax_ 68 integer :: numtache 69 integer :: nb_pt 70 real :: petit_H=0.001 ! pour test ice sur zone flottante 19 USE module3D_phy, only:nx,ny 20 21 implicit none 22 23 real :: surnet !< surnet hauteur de glace au dessus de la mer 24 real :: archim !< test de flottaison 25 ! real, parameter :: Hmin=1.001 !< Hmin pour être considere comme point ice 26 27 integer:: itestf 28 29 logical,dimension(nx,ny) :: gz1mx,gz1my 30 logical,dimension(nx,ny) :: fl1mx,fl1my 31 32 33 real,dimension(nx,ny) :: uxs1 !< uxbar a l'entree de flottab 34 real,dimension(nx,ny) :: uys1 !< uybar a l'entree de flottab 35 36 37 integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y 38 39 40 ! Variables pour la determination des differents shelfs/stream 41 ! (representés comme des taches ou l'on resoud l'eq elliptique) 42 !________________________________________________________________ 43 integer,parameter :: n_ta_max=2000!< nombre de tache max 44 integer,dimension(nx,ny) :: table_out !< pour les numeros des taches 45 integer,dimension(nx,ny) :: tablebis !< pour les numeros des taches 46 integer,dimension(0:n_ta_max) :: compt !< contient les equivalence entre les taches 47 integer,dimension(0:n_ta_max) :: nb_pts_tache !< indique le nombre de points par tache 48 logical,dimension(0:n_ta_max) :: iceberg1D !< T si iceberg, F si calotte posee 49 50 logical,dimension(nx,ny) :: mask_tache_ij !< masque de travail 51 !< vrai pour toute la tache de i,j 52 integer,dimension(2) :: smax_coord !< pour le maxloc des iles 53 54 ! Variables pour determiner le point le plus haut (surf) 55 ! d'une ile completement stream 56 !_________________________________________________________ 57 58 ! icetrim : T si ice stream, F si calotte posee(vertical shear) 59 logical,dimension(n_ta_max) :: icetrim 60 61 integer :: ii,jj 62 integer :: smax_i 63 integer :: smax_j 64 real :: smax_ 65 integer :: numtache 66 integer :: nb_pt 67 real :: petit_H=0.001 ! pour test ice sur zone flottante 71 68 contains 72 ! -----------------------------------------------------------------------------------73 !> SUBROUTINE: flottab()74 !! Cette routine determine les endroits ou la glace75 !! flotte , les iles, et la position du front de glace76 !! @note Il y a 4 sortes de zone77 !! @note - Pose78 ! @note - Grounding zone et streams gzmx et gzmy79 ! @note - Iles ilemx, ilemy80 ! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur81 !>82 subroutine flottab ()69 ! ----------------------------------------------------------------------------------- 70 !> SUBROUTINE: flottab() 71 !! Cette routine determine les endroits ou la glace 72 !! flotte , les iles, et la position du front de glace 73 !! @note Il y a 4 sortes de zone 74 !! @note - Pose 75 ! @note - Grounding zone et streams gzmx et gzmy 76 ! @note - Iles ilemx, ilemy 77 ! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur 78 !> 79 subroutine flottab 83 80 ! 84 81 ! Vince 5 Jan 95 … … 105 102 ! 106 103 ! _________________________________________________________ 107 108 !~ print*,'debut flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel 109 !~ print*,'debut flottab',flot(132,183),ice(132,183) 110 !print*,'H(90,179) flottab 1',H(90,179),ice(90,179), flot(90,179) 111 104 105 use runparam, only:itracebug,nt 106 use module3D_phy, only:shelfy,igrdline,mk_init,flot,H,sealevel_2d,Bsoc,S,H,B,appel_new_flot,& 107 new_flot_point,new_flotmx,new_flotmy,ice,front,frontfacex,frontfacey,isolx,isoly,cotemx,& 108 cotemy,boost,iceberg,uxbar,uybar,mk,gzmx,gzmy,flotmx,flotmy,hmx,hmy,isynchro,ilemx,ilemy,& 109 flgzmx,flgzmy,marine,fbm,bm,bmelt,debug_3D,dt 110 use param_phy_mod, only:row,ro 111 use module_choix, only:mstream_dragging ! subroutine qui depend du dragging 112 113 implicit none 114 115 integer :: i,j 116 112 117 if (itracebug.eq.1) call tracebug(' Entree dans routine flottab') 113 118 … … 223 228 224 229 else if ((H(i,j).LE.0.).and.(archim.LT.0.)) then ! terre deglace qui devient ocean 225 !cdc ice(i,j)=0226 !cdc H(i,j)=1.227 !cdc 1m H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))230 !cdc ice(i,j)=0 231 !cdc H(i,j)=1. 232 !cdc 1m H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01)) 228 233 flot(i,j)=.true. !cdc points ocean sont flot meme sans glace 229 234 H(i,j)=0. … … 246 251 end do 247 252 !$OMP END DO 248 253 249 254 !!$ do i=1,nx 250 255 !!$ do j=1,ny … … 377 382 !$OMP END WORKSHARE 378 383 379 ! afq -- 17/01/19: on supprime les iles.380 ! ! selon x381 ! !$OMP DO382 ! ilesx: do j=2,ny-1383 ! do i=3,nx-2384 ! ! F G F385 ! ! x386 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)387 ! if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &388 ! (sdx(i,j).LT.1.E-02)) then389 ! ilemx(i,j)=.true.390 ! ilemx(i+1,j)=.true.391 392 ! ! F G G F393 ! ! x394 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)395 ! else if ((flot(i-1,j).and..not.flot(i,j) &396 ! .and..not.flot(i+1,j)).and.flot(i+2,j).and. &397 ! (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then398 ! ilemx(i,j)=.true.399 ! ilemx(i+1,j)=.true.400 ! ilemx(i+2,j)=.true.401 402 ! ! F G G F403 ! ! x404 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)405 ! else if ((flot(i-2,j).and..not.flot(i-1,j) &406 ! .and..not.flot(i,j)).and.flot(i+1,j).and. &407 ! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then408 ! ilemx(i-1,j)=.true.409 ! ilemx(i,j)=.true.410 ! ilemx(i+1,j)=.true.411 412 ! ! F G G G F413 ! ! x414 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)415 ! else if ((i.lt.nx-2) &416 ! .and.(flot(i-2,j).and..not.flot(i-1,j) &417 ! .and..not.flot(i,j)).and..not.flot(i+1,j) &418 ! .and.flot(i+2,j).and. &419 ! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 &420 ! .and.sdx(i+1,j).LT.1.E-02)) then421 ! ilemx(i-1,j)=.true.422 ! ilemx(i,j)=.true.423 ! ilemx(i+1,j)=.true.424 ! ilemx(i+2,j)=.true.425 426 ! endif427 428 ! end do429 ! end do ilesx430 ! !$OMP END DO431 432 ! ! selon y433 ! !$OMP DO434 ! ilesy: do j=3,ny-2435 ! do i=2,nx-1436 ! ! F G F437 ! ! x438 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)439 ! if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. &440 ! (sdy(i,j).LT.1.E-02)) then441 ! ilemy(i,j)=.true.442 ! ilemy(i,j+1)=.true.443 444 ! ! F G G F445 ! ! x446 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)447 ! else if ((flot(i,j-1).and..not.flot(i,j) &448 ! .and..not.flot(i,j+1)).and.flot(i,j+2).and. &449 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then450 ! ilemy(i,j)=.true.451 ! ilemy(i,j+1)=.true.452 ! ilemy(i,j+2)=.true.453 454 ! ! F G G F455 ! ! x456 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)457 ! else if ((flot(i,j-2).and..not.flot(i,j-1) &458 ! .and..not.flot(i,j)).and.flot(i,j+1).and. &459 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then460 ! ilemy(i,j-1)=.true.461 ! ilemy(i,j)=.true.462 ! ilemy(i,j+1)=.true.463 464 ! ! F G G G F465 ! ! x466 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)467 ! else if ((j.lt.ny-2) &468 ! .and.(flot(i,j-2).and..not.flot(i,j-1) &469 ! .and..not.flot(i,j)).and..not.flot(i,j+1) &470 ! .and.flot(i,j+2).and. &471 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 &472 ! .and.sdy(i,j+1).LT.1.E-02)) then473 ! ilemy(i,j-1)=.true.474 ! ilemy(i,j)=.true.475 ! ilemy(i,j+1)=.true.476 ! ilemy(i,j+2)=.true.477 ! endif478 ! end do479 ! end do ilesy480 ! !$OMP END DO481 ! ! fin des iles384 ! afq -- 17/01/19: on supprime les iles. 385 ! ! selon x 386 ! !$OMP DO 387 ! ilesx: do j=2,ny-1 388 ! do i=3,nx-2 389 ! ! F G F 390 ! ! x 391 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 392 ! if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. & 393 ! (sdx(i,j).LT.1.E-02)) then 394 ! ilemx(i,j)=.true. 395 ! ilemx(i+1,j)=.true. 396 397 ! ! F G G F 398 ! ! x 399 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 400 ! else if ((flot(i-1,j).and..not.flot(i,j) & 401 ! .and..not.flot(i+1,j)).and.flot(i+2,j).and. & 402 ! (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then 403 ! ilemx(i,j)=.true. 404 ! ilemx(i+1,j)=.true. 405 ! ilemx(i+2,j)=.true. 406 407 ! ! F G G F 408 ! ! x 409 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 410 ! else if ((flot(i-2,j).and..not.flot(i-1,j) & 411 ! .and..not.flot(i,j)).and.flot(i+1,j).and. & 412 ! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then 413 ! ilemx(i-1,j)=.true. 414 ! ilemx(i,j)=.true. 415 ! ilemx(i+1,j)=.true. 416 417 ! ! F G G G F 418 ! ! x 419 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 420 ! else if ((i.lt.nx-2) & 421 ! .and.(flot(i-2,j).and..not.flot(i-1,j) & 422 ! .and..not.flot(i,j)).and..not.flot(i+1,j) & 423 ! .and.flot(i+2,j).and. & 424 ! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 & 425 ! .and.sdx(i+1,j).LT.1.E-02)) then 426 ! ilemx(i-1,j)=.true. 427 ! ilemx(i,j)=.true. 428 ! ilemx(i+1,j)=.true. 429 ! ilemx(i+2,j)=.true. 430 431 ! endif 432 433 ! end do 434 ! end do ilesx 435 ! !$OMP END DO 436 437 ! ! selon y 438 ! !$OMP DO 439 ! ilesy: do j=3,ny-2 440 ! do i=2,nx-1 441 ! ! F G F 442 ! ! x 443 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 444 ! if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. & 445 ! (sdy(i,j).LT.1.E-02)) then 446 ! ilemy(i,j)=.true. 447 ! ilemy(i,j+1)=.true. 448 449 ! ! F G G F 450 ! ! x 451 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 452 ! else if ((flot(i,j-1).and..not.flot(i,j) & 453 ! .and..not.flot(i,j+1)).and.flot(i,j+2).and. & 454 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then 455 ! ilemy(i,j)=.true. 456 ! ilemy(i,j+1)=.true. 457 ! ilemy(i,j+2)=.true. 458 459 ! ! F G G F 460 ! ! x 461 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 462 ! else if ((flot(i,j-2).and..not.flot(i,j-1) & 463 ! .and..not.flot(i,j)).and.flot(i,j+1).and. & 464 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then 465 ! ilemy(i,j-1)=.true. 466 ! ilemy(i,j)=.true. 467 ! ilemy(i,j+1)=.true. 468 469 ! ! F G G G F 470 ! ! x 471 ! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) 472 ! else if ((j.lt.ny-2) & 473 ! .and.(flot(i,j-2).and..not.flot(i,j-1) & 474 ! .and..not.flot(i,j)).and..not.flot(i,j+1) & 475 ! .and.flot(i,j+2).and. & 476 ! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 & 477 ! .and.sdy(i,j+1).LT.1.E-02)) then 478 ! ilemy(i,j-1)=.true. 479 ! ilemy(i,j)=.true. 480 ! ilemy(i,j+1)=.true. 481 ! ilemy(i,j+2)=.true. 482 ! endif 483 ! end do 484 ! end do ilesy 485 ! !$OMP END DO 486 ! ! fin des iles 482 487 483 488 !$OMP END PARALLEL … … 487 492 !---------------------------------------------------------------------------------- 488 493 call mstream_dragging 489 494 490 495 if (itracebug.eq.1) call tracebug(' routine flottab apres call dragging') 491 496 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) … … 544 549 .or.(.not.marine.and.flotmx(:,:)) 545 550 where (hmx(:,:).eq.0.) 546 flgzmx(:,:) = .false.551 flgzmx(:,:) = .false. 547 552 endwhere 548 553 flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:))) & 549 554 .or.(.not.marine.and.flotmy(:,:)) 550 555 where (hmy(:,:).eq.0.) 551 flgzmy(:,:) = .false.556 flgzmy(:,:) = .false. 552 557 endwhere 553 558 !$OMP END WORKSHARE … … 584 589 !!$ end do 585 590 !!$end do 586 ! print*, 'flolottab debug', H(71,25),flot(71,25),bm(71,25),ice(71,25),time587 !print*,'H(90,179) flottab 2',H(90,179),ice(90,179),flot(90,179)591 ! print*, 'flolottab debug', H(71,25),flot(71,25),bm(71,25),ice(71,25),time 592 !print*,'H(90,179) flottab 2',H(90,179),ice(90,179),flot(90,179) 588 593 !$OMP WORKSHARE 589 594 where (flot(:,:)) 590 !cdc 1m where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:)))591 where (H(:,:).gt.max(BM(:,:)-Bmelt(:,:)+petit_H,0.)*dt)592 !cdc where (H(:,:).gt.0.)595 !cdc 1m where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:))) 596 where (H(:,:).gt.max(BM(:,:)-Bmelt(:,:)+petit_H,0.)*dt) 597 !cdc where (H(:,:).gt.0.) 593 598 ice(:,:)=1 594 599 elsewhere … … 607 612 !$OMP END WORKSHARE 608 613 !$OMP END PARALLEL 609 ! print*,'flottab',time,H(191,81),bm(191,81),bmelt(191,81),(BM(191,81)-Bmelt(191,81)+petit_H)*dt,ice(191,81)610 !print*,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)*ro/row -sealevel611 614 ! print*,'flottab',time,H(191,81),bm(191,81),bmelt(191,81),(BM(191,81)-Bmelt(191,81)+petit_H)*dt,ice(191,81) 615 !print*,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)*ro/row -sealevel 616 612 617 !call determin_front ! cette version ne conserve pas la masse !!! 613 618 call determin_front_tof ! version simplifiee … … 617 622 debug_3D(:,:,119) = ice(:,:)*mk(:,:) 618 623 624 625 end subroutine flottab 626 !-------------------------------------------------------------------- 627 628 !> SUBROUTINE: determin_tache 629 !! Routine pour la dtermination du numero de tache a effectuer 630 !> 631 subroutine determin_tache 619 632 620 end subroutine flottab 621 !-------------------------------------------------------------------- 622 623 !> SUBROUTINE: determin_tache 624 !! Routine pour la dtermination du numero de tache a effectuer 625 !> 626 subroutine determin_tache 627 628 !$ USE OMP_LIB 629 630 implicit none 631 integer :: indice 632 integer :: label ! no des taches rencontrées dans le mask 633 integer :: label_max ! no temporaire maxi de tache rencontrées 634 ! integer :: mask_nb = 4 635 integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales 636 ! integer,dimension(mask_nb) :: mask 637 integer,dimension(mask_nb) :: mask 638 639 640 ! 1-initialisation 641 !----------------- 642 label_max=1 ! numero de la tache, la premiere tache est notée 1 643 label=1 644 do i=1,n_ta_max 645 compt(i)=i 646 enddo 647 ! table_in = .false. 648 !$OMP PARALLEL 649 !$OMP WORKSHARE 650 table_out(:,:) = 0 651 iceberg1D(:) = .true. 652 icetrim (:) = .true. 653 nb_pts_tache(:) = 0 654 !$OMP END WORKSHARE 655 !$OMP END PARALLEL 656 ! open(unit=100,file="tache.data",status='replace') 657 658 ! 2-reperage des taches 659 !---------------------- 660 661 do j=2,ny-1 662 do i=2,nx-1 663 664 IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------! 665 666 if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes 667 ! un des voisins est deja en glace 668 mask(1) = table_out(i-1,j) 669 mask(2) = table_out(i,j-1) 670 label = label_max 671 672 !on determine la valeur de la tache minimun (>0) presente ds le masque 673 do indice=1,mask_nb 674 if (mask(indice).gt.0) label=min(label,mask(indice)) 675 enddo 676 !cdc label=min(label,minval(mask(:), mask=mask > 0)) 677 678 !on fixe la valeur de la tache voisine minimun au point etudie (via label) 679 table_out(i,j)=label 680 !si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F. 681 if (.not.FLOT(I,J)) then 682 iceberg1D(label)=.false. 683 endif 684 685 !si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F. 686 if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & 687 (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then 688 icetrim(label)=.false. 689 endif 690 691 ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt 692 ! on lui affecte le numero de la tache fondamentale 693 694 do indice=1,mask_nb 695 if(mask(indice).gt.label) then 696 compt(mask(indice))=-label 697 endif 698 enddo 699 ! exemple on est sur le point X : 5 X 700 do indice=1,mask_nb ! 20 701 if(mask(indice).gt.label) then ! mask(2)=20 > 5 702 compt(mask(indice))=label ! compt(20)=5 703 if (.not.iceberg1D(mask(indice))) iceberg1D(label)=.false. ! si la tache n'etais pas un iceberg => iceberg =.false. iceberg(5)=.false. 704 if (.not.icetrim(mask(indice))) icetrim(label)=.false. 705 where (table_out(:,:).eq.mask(indice)) ! where table_out(:,:)=mask(2)=20 706 table_out(:,:)=label ! table_out(:,:)=label=5 707 endwhere 708 endif 709 enddo 710 711 else !aucun des voisins est une tache 712 table_out(i,j)= label_max 713 compt(label_max)=label_max 714 !si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F. 715 if (.not.FLOT(I,J)) then 716 iceberg1D(label_max)=.false. 717 endif 718 719 !si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F. 720 if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & 721 (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then 722 icetrim(label)=.false. 723 endif 724 725 label_max = label_max+1 726 if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches icebergs=',label_max 727 endif 728 729 730 else !on est pas sur une tache---------------------------------------------- 731 table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) 732 endif !--------------------------------------------------------------------- 733 734 735 enddo 736 enddo 737 738 739 740 ! On reorganise compt en ecrivant le numero de la tache fondamentale 741 ! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité) 742 ! On indique aussi le nb de point que contient chaque taches (nb_pts_tache) 743 744 !$OMP PARALLEL 745 !$OMP DO 746 do j=1,ny 747 do i=1,nx 748 if (table_out(i,j).ne.0) then 749 nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1 750 endif 633 use module3D_phy, only:ice,flot,gzmx,gzmy,S,flgzmx,flgzmy,sealevel,num_tracebug,time,iceberg,& 634 debug_3D 635 636 !$ USE OMP_LIB 637 638 implicit none 639 integer :: i,j 640 integer :: indice 641 integer :: label ! no des taches rencontrées dans le mask 642 integer :: label_max ! no temporaire maxi de tache rencontrées 643 ! integer :: mask_nb = 4 644 integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales 645 ! integer,dimension(mask_nb) :: mask 646 integer,dimension(mask_nb) :: mask 647 648 649 ! 1-initialisation 650 !----------------- 651 label_max=1 ! numero de la tache, la premiere tache est notée 1 652 label=1 653 do i=1,n_ta_max 654 compt(i)=i 751 655 enddo 752 enddo 753 !$OMP END DO 754 !$OMP END PARALLEL 755 756 !On compte comme englacé uniquement les calottes dont une partie est posée 757 !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 758 !$OMP DO 759 do j=3,ny-2 760 do i=3,nx-2 761 test1: if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg 762 if (nb_pts_tache(table_out(i,j)).ge.1) then 763 ice(i,j)=1 764 ! ici on est sur une tache non iceberg >= 5 points 765 ! on teste si la tache n'est pas completement ice stream 656 ! table_in = .false. 657 !$OMP PARALLEL 658 !$OMP WORKSHARE 659 table_out(:,:) = 0 660 iceberg1D(:) = .true. 661 icetrim (:) = .true. 662 nb_pts_tache(:) = 0 663 !$OMP END WORKSHARE 664 !$OMP END PARALLEL 665 ! open(unit=100,file="tache.data",status='replace') 666 667 ! 2-reperage des taches 668 !---------------------- 669 670 do j=2,ny-1 671 do i=2,nx-1 672 673 IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------! 674 675 if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes 676 ! un des voisins est deja en glace 677 mask(1) = table_out(i-1,j) 678 mask(2) = table_out(i,j-1) 679 label = label_max 680 681 !on determine la valeur de la tache minimun (>0) presente ds le masque 682 do indice=1,mask_nb 683 if (mask(indice).gt.0) label=min(label,mask(indice)) 684 enddo 685 !cdc label=min(label,minval(mask(:), mask=mask > 0)) 686 687 !on fixe la valeur de la tache voisine minimun au point etudie (via label) 688 table_out(i,j)=label 689 !si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F. 690 if (.not.FLOT(I,J)) then 691 iceberg1D(label)=.false. 692 endif 693 694 !si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F. 695 if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & 696 (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then 697 icetrim(label)=.false. 698 endif 699 700 ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt 701 ! on lui affecte le numero de la tache fondamentale 702 703 do indice=1,mask_nb 704 if(mask(indice).gt.label) then 705 compt(mask(indice))=-label 706 endif 707 enddo 708 ! exemple on est sur le point X : 5 X 709 do indice=1,mask_nb ! 20 710 if(mask(indice).gt.label) then ! mask(2)=20 > 5 711 compt(mask(indice))=label ! compt(20)=5 712 if (.not.iceberg1D(mask(indice))) iceberg1D(label)=.false. ! si la tache n'etais pas un iceberg => iceberg =.false. iceberg(5)=.false. 713 if (.not.icetrim(mask(indice))) icetrim(label)=.false. 714 where (table_out(:,:).eq.mask(indice)) ! where table_out(:,:)=mask(2)=20 715 table_out(:,:)=label ! table_out(:,:)=label=5 716 endwhere 717 endif 718 enddo 719 720 else !aucun des voisins est une tache 721 table_out(i,j)= label_max 722 compt(label_max)=label_max 723 !si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F. 724 if (.not.FLOT(I,J)) then 725 iceberg1D(label_max)=.false. 726 endif 727 728 !si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F. 729 if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & 730 (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then 731 icetrim(label)=.false. 732 endif 733 734 label_max = label_max+1 735 if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches icebergs=',label_max 736 endif 737 738 739 else !on est pas sur une tache---------------------------------------------- 740 table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) 741 endif !--------------------------------------------------------------------- 742 743 744 enddo 745 enddo 746 747 748 749 ! On reorganise compt en ecrivant le numero de la tache fondamentale 750 ! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité) 751 ! On indique aussi le nb de point que contient chaque taches (nb_pts_tache) 752 753 !$OMP PARALLEL 754 !$OMP DO 755 do j=1,ny 756 do i=1,nx 757 if (table_out(i,j).ne.0) then 758 nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1 759 endif 760 enddo 761 enddo 762 !$OMP END DO 763 !$OMP END PARALLEL 764 765 !On compte comme englacé uniquement les calottes dont une partie est posée 766 !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 767 !$OMP DO 768 do j=3,ny-2 769 do i=3,nx-2 770 test1: if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg 771 if (nb_pts_tache(table_out(i,j)).ge.1) then 772 ice(i,j)=1 773 ! ici on est sur une tache non iceberg >= 5 points 774 ! on teste si la tache n'est pas completement ice stream 775 776 test2: if (icetrim(table_out(i,j))) then ! on a une ile d'ice stream 777 778 mask_tache_ij(:,:)=.false. 779 mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j)) ! pour toute la tache 780 781 smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 782 smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 783 smax_i=smax_coord(1) 784 smax_j=smax_coord(2) 785 786 gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 787 gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 788 flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 789 flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 790 791 if (Smax_.le.sealevel) then ! afq -- WARNING: en fait avec le niveau marin local cest plus complique que ca! 792 !if (Smax_.le.sealevel_2d(i,j)) then ! afq --WARNING: il faudrait regarder point par point sur la tache... 793 write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 794 write(num_tracebug,*)'time=',time,' coord:',smax_i,smax_j 795 end if 796 endif test2 797 end if ! endif deplace 798 !cdc transfere dans calving : 799 else ! on est sur un iceberg ! test1 800 iceberg(i,j)=iceberg1D(table_out(i,j)) 801 !~ ice(i,j)=0 802 !~ h(i,j)=0. !1. afq, we should put everything in calving! 803 !~ surnet=H(i,j)*(1.-ro/row) 804 !~ S(i,j)=surnet+sealevel 805 !~ B(i,j)=S(i,j)-H(i,j) 806 endif test1 807 end do 808 end do 809 !$OMP END DO 810 !$OMP END PARALLEL 811 812 debug_3D(:,:,124)=real(table_out(:,:)) 813 814 end subroutine determin_tache 815 !---------------------------------------------------------------------- 816 !> SUBROUTINE: determin_front 817 !!Routine pour la determination du front 818 !> 819 subroutine determin_front 820 821 use module3D_phy, only:ice,H,front,frontfacex,frontfacey,isolx,isoly 822 766 823 767 test2: if (icetrim(table_out(i,j))) then ! on a une ile d'ice stream 824 !!$ USE OMP_LIB 825 implicit none 768 826 769 mask_tache_ij(:,:)=.false. 770 mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j)) ! pour toute la tache 771 772 smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 773 smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 774 smax_i=smax_coord(1) 775 smax_j=smax_coord(2) 776 777 gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 778 gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 779 flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 780 flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 781 782 if (Smax_.le.sealevel) then ! afq -- WARNING: en fait avec le niveau marin local cest plus complique que ca! 783 !if (Smax_.le.sealevel_2d(i,j)) then ! afq --WARNING: il faudrait regarder point par point sur la tache... 784 write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 785 write(num_tracebug,*)'time=',time,' coord:',smax_i,smax_j 786 end if 787 endif test2 788 end if ! endif deplace 789 !cdc transfere dans calving : 790 else ! on est sur un iceberg ! test1 791 iceberg(i,j)=iceberg1D(table_out(i,j)) 792 !~ ice(i,j)=0 793 !~ h(i,j)=0. !1. afq, we should put everything in calving! 794 !~ surnet=H(i,j)*(1.-ro/row) 795 !~ S(i,j)=surnet+sealevel 796 !~ B(i,j)=S(i,j)-H(i,j) 797 endif test1 798 end do 799 end do 800 !$OMP END DO 801 !$OMP END PARALLEL 802 803 debug_3D(:,:,124)=real(table_out(:,:)) 804 805 end subroutine determin_tache 806 !---------------------------------------------------------------------- 807 !> SUBROUTINE: determin_front 808 !!Routine pour la determination du front 809 !> 810 subroutine determin_front 811 !!$ USE OMP_LIB 812 integer :: i_moins1,i_plus1,i_plus2 813 integer :: j_moins1,j_plus1,j_plus2 814 815 !$OMP PARALLEL 816 !$OMP DO 817 do j=3,ny-2 818 do i=3,nx-2 819 820 surice:if (ice(i,j).eq.0) then 821 do pmx=-1,1,2 822 do pmy=-1,1,2 823 824 diagice : if (ice(i+pmx,j+pmy).eq.1) then 825 826 if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues 827 ! de glaces diagonales en coin(26dec04) 828 if ((ice(i+2*pmx,j).eq.1.and.ice(i+2*pmx,j+pmy).eq.0).or.& 829 (ice(i,j+2*pmy).eq.1.and.ice(i+pmx,j+2*pmy).eq.0)) then 830 ice(i,j)=1 831 h(i,j)=max(1.,h(i,j)) 832 endif 833 834 ! test (i) pour eviter les langues de glaces diagonales : 835 ! mouvement du cheval aux echecs 836 837 if ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.1)) then 838 if (ice(i+2*pmx,j+pmy).eq.1.and. & 839 (ice(i+2*pmx,j+2*pmy)+ice(i,j+2*pmy)).ge.1) then 840 ice(i,j)=1 841 h(i,j)=max(1.,h(i,j)) 842 endif 843 if (ice(i+pmx,j+2*pmy).eq.1.and. & 844 (ice(i+2*pmx,j+2*pmy)+ice(i+2*pmx,j)).ge.1) then 845 ice(i,j)=1 846 h(i,j)=max(1.,h(i,j)) 847 endif 848 849 ! test (ii) pour eviter les langues de glaces diagonales : 850 ! le point glace ice(i+pmx,j+pmy) a : 851 ! - ses 4 voisins frontaux en glace 852 ! - mais 2 voisins vides diagonalement opposes 853 854 elseif ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.2) & 855 .and.ice(i+2*pmx,j+2*pmy).eq.0) then 856 857 ! test (iii) pour faire les tests (i) et (ii) 858 ice(i,j)=1 859 h(i,j)=max(1.,h(i,j)) 860 ice(i+2*pmx,j+2*pmy)=1 861 h(i+2*pmx,j+2*pmy)=max(1.,h(i+2*pmx,j+2*pmy)) 862 endif 863 endif 864 endif diagice 865 enddo 866 enddo 867 endif surice 868 end do 869 end do 870 !$OMP END DO 871 !$OMP ENd PARALLEL 827 integer :: i,j,k 828 integer :: i_moins1,i_plus1,i_plus2 829 integer :: j_moins1,j_plus1,j_plus2 830 831 !$OMP PARALLEL 832 !$OMP DO 833 do j=3,ny-2 834 do i=3,nx-2 835 836 surice:if (ice(i,j).eq.0) then 837 do pmx=-1,1,2 838 do pmy=-1,1,2 839 840 diagice : if (ice(i+pmx,j+pmy).eq.1) then 841 842 if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues 843 ! de glaces diagonales en coin(26dec04) 844 if ((ice(i+2*pmx,j).eq.1.and.ice(i+2*pmx,j+pmy).eq.0).or.& 845 (ice(i,j+2*pmy).eq.1.and.ice(i+pmx,j+2*pmy).eq.0)) then 846 ice(i,j)=1 847 h(i,j)=max(1.,h(i,j)) 848 endif 849 850 ! test (i) pour eviter les langues de glaces diagonales : 851 ! mouvement du cheval aux echecs 852 853 if ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.1)) then 854 if (ice(i+2*pmx,j+pmy).eq.1.and. & 855 (ice(i+2*pmx,j+2*pmy)+ice(i,j+2*pmy)).ge.1) then 856 ice(i,j)=1 857 h(i,j)=max(1.,h(i,j)) 858 endif 859 if (ice(i+pmx,j+2*pmy).eq.1.and. & 860 (ice(i+2*pmx,j+2*pmy)+ice(i+2*pmx,j)).ge.1) then 861 ice(i,j)=1 862 h(i,j)=max(1.,h(i,j)) 863 endif 864 865 ! test (ii) pour eviter les langues de glaces diagonales : 866 ! le point glace ice(i+pmx,j+pmy) a : 867 ! - ses 4 voisins frontaux en glace 868 ! - mais 2 voisins vides diagonalement opposes 869 870 elseif ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.2) & 871 .and.ice(i+2*pmx,j+2*pmy).eq.0) then 872 873 ! test (iii) pour faire les tests (i) et (ii) 874 ice(i,j)=1 875 h(i,j)=max(1.,h(i,j)) 876 ice(i+2*pmx,j+2*pmy)=1 877 h(i+2*pmx,j+2*pmy)=max(1.,h(i+2*pmx,j+2*pmy)) 878 endif 879 endif 880 endif diagice 881 enddo 882 enddo 883 endif surice 884 end do 885 end do 886 !$OMP END DO 887 !$OMP ENd PARALLEL 872 888 873 889 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) … … 881 897 882 898 883 ! print*,'dans remplissage baies',time884 885 baies: do k=1,2886 !$OMP PARALLEL887 !$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2)888 do j=1,ny889 do i=1,nx890 891 surice_xy: if (ice(i,j).eq.0) then892 i_moins1=max(i-1,1)893 j_moins1=max(j-1,1)894 i_plus1=min(i+1,nx)895 j_plus1=min(j+1,ny)896 i_plus2=min(i+2,nx)897 j_plus2=min(j+2,ny)898 899 ! test (iii) pour trouver les baies de largeur 1 ou 2 cases900 ! et combler les trous si ce sont des baies901 ! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles902 ! baies horizontales903 904 if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then905 if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1906 ice(i,j)=1907 H(i,j)=max(1.,H(i,j))908 endif909 endif910 911 912 if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then913 if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1914 ice(i,j)=1915 H(i,j)=max(1.,H(i,j))916 endif917 endif918 919 endif surice_xy920 end do921 end do922 !$OMP END DO923 !$OMP END PARALLEL924 end do baies899 ! print*,'dans remplissage baies',time 900 901 baies: do k=1,2 902 !$OMP PARALLEL 903 !$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2) 904 do j=1,ny 905 do i=1,nx 906 907 surice_xy: if (ice(i,j).eq.0) then 908 i_moins1=max(i-1,1) 909 j_moins1=max(j-1,1) 910 i_plus1=min(i+1,nx) 911 j_plus1=min(j+1,ny) 912 i_plus2=min(i+2,nx) 913 j_plus2=min(j+2,ny) 914 915 ! test (iii) pour trouver les baies de largeur 1 ou 2 cases 916 ! et combler les trous si ce sont des baies 917 ! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles 918 ! baies horizontales 919 920 if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then 921 if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1 922 ice(i,j)=1 923 H(i,j)=max(1.,H(i,j)) 924 endif 925 endif 926 927 928 if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then 929 if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1 930 ice(i,j)=1 931 H(i,j)=max(1.,H(i,j)) 932 endif 933 endif 934 935 endif surice_xy 936 end do 937 end do 938 !$OMP END DO 939 !$OMP END PARALLEL 940 end do baies 925 941 926 942 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) … … 933 949 !!$end if 934 950 935 !$OMP PARALLEL936 !$OMP DO937 do j=2,ny-1938 do i=2,nx-1939 940 if (ice(i,j).eq.1) then ! test si ice=1941 942 ! if ice, on determine front...943 ! ainsi, front=0 sur les zones = 0944 945 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))946 !front= le nb de faces en contact avec un voisin englacé947 endif948 end do949 end do950 !$OMP END DO951 952 ! traitement des bords. On considere que l'exterieur n'a pas de glace953 ! attention ce n'est vrai que sur la grande grille954 955 !$OMP DO PRIVATE(i)956 do j=2,ny-1957 i=1958 front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))959 i=nx960 front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))961 end do962 !$OMP END DO963 964 !$OMP DO PRIVATE(j)965 do i=2,nx-1966 j=1967 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))968 j=ny969 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))970 end do971 !$OMP END DO972 973 ! traitement des coins974 975 front(1,1)=ice(2,1)+ice(2,1)976 front(1,ny)=ice(2,ny)+ice(1,ny-1)977 front(nx,1)=ice(nx,2)+ice(nx-1,1)978 front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)951 !$OMP PARALLEL 952 !$OMP DO 953 do j=2,ny-1 954 do i=2,nx-1 955 956 if (ice(i,j).eq.1) then ! test si ice=1 957 958 ! if ice, on determine front... 959 ! ainsi, front=0 sur les zones = 0 960 961 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 962 !front= le nb de faces en contact avec un voisin englacé 963 endif 964 end do 965 end do 966 !$OMP END DO 967 968 ! traitement des bords. On considere que l'exterieur n'a pas de glace 969 ! attention ce n'est vrai que sur la grande grille 970 971 !$OMP DO PRIVATE(i) 972 do j=2,ny-1 973 i=1 974 front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 975 i=nx 976 front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1)) 977 end do 978 !$OMP END DO 979 980 !$OMP DO PRIVATE(j) 981 do i=2,nx-1 982 j=1 983 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)) 984 j=ny 985 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1)) 986 end do 987 !$OMP END DO 988 989 ! traitement des coins 990 991 front(1,1)=ice(2,1)+ice(2,1) 992 front(1,ny)=ice(2,ny)+ice(1,ny-1) 993 front(nx,1)=ice(nx,2)+ice(nx-1,1) 994 front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny) 979 995 980 996 !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) … … 987 1003 !!$end if 988 1004 989 ! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales) 990 ! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs 991 ! de n'importe quelle taille). 992 ! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites 993 994 !$OMP DO 995 do j=1,ny 996 do i=1,nx-1 997 if (front(i,j).eq.1) then 998 if (front(i+1,j).eq.1) then 999 ice(i,j)=0 1000 ice(i+1,j)=0 1001 front(i,j)=0 1002 front(i+1,j)=0 1003 endif 1004 endif 1005 end do 1006 end do 1007 !$OMP END DO 1008 1009 !$OMP DO 1010 do j=1,ny-1 1011 do i=1,nx 1012 if (front(i,j).eq.1) then 1013 if (front(i,j+1).eq.1) then 1014 ice(i,j)=0 1015 ice(i,j+1)=0 1016 front(i,j)=0 1017 front(i,j+1)=0 1018 endif 1019 end if 1020 end do 1021 end do 1022 !$OMP END DO 1023 1024 !isolx signifie pas de voisins en x 1025 !isoly signifie pas de voisins en y 1026 !remarque : 1027 !si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0) 1028 1029 ! calcul de frontfacex et isolx 1030 !$OMP DO 1031 do j=1,ny 1032 do i=2,nx-1 1033 1034 if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) 1035 1036 if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x 1037 1038 if ((ice(i-1,j)+ice(i+1,j)).eq.0) then 1039 isolx(i,j)=.true. 1040 elseif (ice(i-1,j).eq.0) then 1041 frontfacex(i,j)=-1 ! front i-1 |i i+1 1042 else 1043 frontfacex(i,j)=+1 ! front i-1 i| i+1 1044 endif 1045 endif 1046 end if !fin du test il y a un front 1047 1048 end do 1049 end do 1050 !$OMP END DO 1051 1052 ! calcul de frontfacey et isoly 1053 !$OMP DO 1054 do j=2,ny-1 1055 do i=1,nx 1056 1057 if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) 1058 1059 if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y 1060 1061 if ((ice(i,j-1)+ice(i,j+1)).eq.0) then 1062 isoly(i,j)=.true. !front j-1 |j| j+1 1063 elseif (ice(i,j-1).eq.0) then 1064 frontfacey(i,j)=-1 !front j-1 |j j+1 1065 else 1066 frontfacey(i,j)=+1 !front j-1 j| j+1 1067 endif 1068 endif 1069 end if !fin du test il y a un front 1070 1071 end do 1072 end do 1073 !$OMP END DO 1074 1075 1076 ! traitement des bords. On considere que l'exterieur n'a pas de glace 1077 ! attention ce n'est vrai que sur la grande grille 1078 1079 !$OMP DO PRIVATE(i) 1080 do j=2,ny-1 1081 i=1 1082 if (front(i,j).ge.1) then 1083 if (ice(i+1,j).eq.0) then 1084 isolx(i,j)=.true. 1085 else 1086 frontfacex(i,j)=-1 1087 endif 1088 end if 1089 i=nx 1090 if (front(i,j).ge.1) then 1091 if (ice(i-1,j).eq.0) then 1092 isolx(i,j)=.true. 1093 else 1094 frontfacex(i,j)=1 1095 endif 1096 end if 1097 end do 1098 !$OMP END DO 1099 1100 !$OMP DO PRIVATE(j) 1101 do i=2,nx-1 1102 j=1 1103 if (front(i,j).ge.1) then 1104 if (ice(i,j+1).eq.0) then 1105 isoly(i,j)=.true. 1106 else 1107 frontfacey(i,j)=-1 1108 endif 1109 end if 1110 j=ny 1111 if (front(i,j).ge.1) then 1112 if (ice(i,j-1).eq.0) then 1113 isoly(i,j)=.true. 1114 else 1115 frontfacey(i,j)=1 1116 endif 1117 end if 1118 end do 1119 !$OMP END DO 1120 !$OMP END PARALLEL 1121 1122 return 1123 end subroutine determin_front 1124 !------------------------------------------------------------------------------ 1125 1126 subroutine determin_front_tof 1127 1128 integer,dimension(nx,ny) :: nofront ! tableau de travail (points au dela du front) 1129 !$OMP PARALLEL 1130 !$OMP DO 1131 do j=2,ny-1 1132 do i=2,nx-1 1133 1134 if (ice(i,j).eq.1) then ! test si ice=1 1135 1136 ! if ice, on determine front... 1137 ! ainsi, front=0 sur les zones = 0 1138 1139 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 1140 !front= le nb de faces en contact avec un voisin englacé 1141 endif 1142 end do 1143 end do 1144 !$OMP END DO 1145 1146 ! traitement des bords. On considere que l'exterieur n'a pas de glace 1147 ! attention ce n'est vrai que sur la grande grille 1148 1149 !$OMP DO PRIVATE(i) 1150 do j=2,ny-1 1151 i=1 1152 front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 1153 i=nx 1154 front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1)) 1155 end do 1156 !$OMP END DO 1157 1158 !$OMP DO PRIVATE(j) 1159 do i=2,nx-1 1160 j=1 1161 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)) 1162 j=ny 1163 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1)) 1164 end do 1165 !$OMP END DO 1166 !$OMP BARRIER 1167 ! traitement des coins 1168 1169 front(1,1)=ice(2,1)+ice(2,1) 1170 front(1,ny)=ice(2,ny)+ice(1,ny-1) 1171 front(nx,1)=ice(nx,2)+ice(nx-1,1) 1172 front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny) 1173 1174 ! Les points à plus d'un point de grille du bord sont front=-1 1175 !$OMP WORKSHARE 1176 nofront(:,:)=0 1177 !$OMP END WORKSHARE 1178 !$OMP BARRIER 1179 1180 !$OMP DO 1181 do j=2,ny-1 1182 do i=2,nx-1 1183 if ((ice(i,j).eq.0).and.(front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).eq.0) then 1184 nofront(i,j)=-1 1185 endif 1186 enddo 1187 enddo 1188 !$OMP END DO 1189 !$OMP BARRIER 1190 1191 !$OMP WORKSHARE 1192 where (nofront(:,:).eq.-1) 1193 front(:,:)=-1 1194 endwhere 1195 !$OMP END WORKSHARE 1196 !$OMP END PARALLEL 1197 1198 end subroutine determin_front_tof 1199 1200 1201 !> SUBROUTINE: determin_marais 1202 !! afq -- Routine pour l'identification des "marais" 1203 !! Un marais est une tache "shelf" entouré de points grounded 1204 !! Copie sauvage de determin_tache, adapte au probleme du marais 1205 !> 1206 subroutine determin_marais 1207 1208 !$ USE OMP_LIB 1209 1210 implicit none 1211 1212 integer :: indice 1213 integer :: label ! no des taches rencontrées dans le mask 1214 integer :: label_max ! no temporaire maxi de tache rencontrées 1215 integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales 1216 integer,dimension(mask_nb) :: mask ! numero de tache des points adjacents 1217 1218 integer,dimension(nx,ny) :: table_out_marais !< numeros de tache d'un point ij 1219 integer,dimension(0:n_ta_max) :: compt_marais !< contient les equivalence entre les taches 1220 integer,dimension(0:n_ta_max) :: nb_pts_marais !< indique le nombre de points par tache 1221 logical,dimension(0:n_ta_max) :: marais !< T si flottants entoure de poses, F sinon 1222 1223 1224 ! 1-initialisation 1225 !----------------- 1226 label_max=1 ! numero de la tache, la premiere tache est notée 1 1227 label=1 1228 do i=1,n_ta_max 1229 compt_marais(i)=i 1230 enddo 1231 !$OMP PARALLEL 1232 !$OMP WORKSHARE 1233 table_out_marais(:,:) = 0 1234 marais(:) = .true. 1235 nb_pts_marais(:) = 0 1236 !$OMP END WORKSHARE 1237 !$OMP END PARALLEL 1238 1239 ! 2-reperage des taches 1240 !---------------------- 1241 1242 do j=2,ny-1 1005 ! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales) 1006 ! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs 1007 ! de n'importe quelle taille). 1008 ! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites 1009 1010 !$OMP DO 1011 do j=1,ny 1012 do i=1,nx-1 1013 if (front(i,j).eq.1) then 1014 if (front(i+1,j).eq.1) then 1015 ice(i,j)=0 1016 ice(i+1,j)=0 1017 front(i,j)=0 1018 front(i+1,j)=0 1019 endif 1020 endif 1021 end do 1022 end do 1023 !$OMP END DO 1024 1025 !$OMP DO 1026 do j=1,ny-1 1027 do i=1,nx 1028 if (front(i,j).eq.1) then 1029 if (front(i,j+1).eq.1) then 1030 ice(i,j)=0 1031 ice(i,j+1)=0 1032 front(i,j)=0 1033 front(i,j+1)=0 1034 endif 1035 end if 1036 end do 1037 end do 1038 !$OMP END DO 1039 1040 !isolx signifie pas de voisins en x 1041 !isoly signifie pas de voisins en y 1042 !remarque : 1043 !si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0) 1044 1045 ! calcul de frontfacex et isolx 1046 !$OMP DO 1047 do j=1,ny 1048 do i=2,nx-1 1049 1050 if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) 1051 1052 if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x 1053 1054 if ((ice(i-1,j)+ice(i+1,j)).eq.0) then 1055 isolx(i,j)=.true. 1056 elseif (ice(i-1,j).eq.0) then 1057 frontfacex(i,j)=-1 ! front i-1 |i i+1 1058 else 1059 frontfacex(i,j)=+1 ! front i-1 i| i+1 1060 endif 1061 endif 1062 end if !fin du test il y a un front 1063 1064 end do 1065 end do 1066 !$OMP END DO 1067 1068 ! calcul de frontfacey et isoly 1069 !$OMP DO 1070 do j=2,ny-1 1071 do i=1,nx 1072 1073 if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) 1074 1075 if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y 1076 1077 if ((ice(i,j-1)+ice(i,j+1)).eq.0) then 1078 isoly(i,j)=.true. !front j-1 |j| j+1 1079 elseif (ice(i,j-1).eq.0) then 1080 frontfacey(i,j)=-1 !front j-1 |j j+1 1081 else 1082 frontfacey(i,j)=+1 !front j-1 j| j+1 1083 endif 1084 endif 1085 end if !fin du test il y a un front 1086 1087 end do 1088 end do 1089 !$OMP END DO 1090 1091 1092 ! traitement des bords. On considere que l'exterieur n'a pas de glace 1093 ! attention ce n'est vrai que sur la grande grille 1094 1095 !$OMP DO PRIVATE(i) 1096 do j=2,ny-1 1097 i=1 1098 if (front(i,j).ge.1) then 1099 if (ice(i+1,j).eq.0) then 1100 isolx(i,j)=.true. 1101 else 1102 frontfacex(i,j)=-1 1103 endif 1104 end if 1105 i=nx 1106 if (front(i,j).ge.1) then 1107 if (ice(i-1,j).eq.0) then 1108 isolx(i,j)=.true. 1109 else 1110 frontfacex(i,j)=1 1111 endif 1112 end if 1113 end do 1114 !$OMP END DO 1115 1116 !$OMP DO PRIVATE(j) 1243 1117 do i=2,nx-1 1244 if ((ice(i,j).ge.1).and.flot(i,j)) then ! on est sur la glace qui flotte-------------------! 1245 1246 if (((ice(i-1,j).ge.1).and.flot(i-1,j)).or.((ice(i,j-1).ge.1).and.flot(i,j-1))) then !masque de 2 cases adjacentes 1247 ! un des voisins est deja en glace 1248 mask(1) = table_out_marais(i-1,j) 1249 mask(2) = table_out_marais(i,j-1) 1250 label = label_max 1251 1252 ! on determine la valeur de la tache minimun (>0) presente ds le masque 1253 do indice=1,mask_nb 1254 if (mask(indice).gt.0) label=min(label,mask(indice)) 1255 enddo 1256 1257 ! on fixe la valeur de la tache voisine minimun au point etudie (via label) 1258 table_out_marais(i,j)=label 1259 1260 !si un des voisins n'est pas glace alors la tache n'est pas un marais 1261 if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then 1262 marais(label)=.false. 1118 j=1 1119 if (front(i,j).ge.1) then 1120 if (ice(i,j+1).eq.0) then 1121 isoly(i,j)=.true. 1122 else 1123 frontfacey(i,j)=-1 1263 1124 endif 1264 1265 ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt_marais 1266 ! on lui affecte le numero de la tache fondamentale 1267 1268 ! exemple on est sur le point X : 5 X 1269 do indice=1,mask_nb ! 20 1270 if(mask(indice).gt.label) then ! mask(2)=20 > 5 1271 compt_marais(mask(indice))=label ! compt_marais(20)=5 1272 if (.not.marais(mask(indice))) marais(label)=.false. ! si la tache n'etais pas un marais => marais =.false. marais(-(-5))=.false. 1273 where (table_out_marais(:,:).eq.mask(indice)) ! where table_out_marais(:,:)=mask(2)=20 1274 table_out_marais(:,:)=label ! table_out_marais(:,:)=label=5 1275 endwhere 1276 endif 1277 enddo 1278 1279 else !aucun des voisins est une tache 1280 table_out_marais(i,j)= label_max 1281 compt_marais(label_max)=label_max 1282 1283 ! si un des voisins n'est pas glace alors la tache n'est pas un marais 1284 if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then 1285 marais(label_max)=.false. 1125 end if 1126 j=ny 1127 if (front(i,j).ge.1) then 1128 if (ice(i,j-1).eq.0) then 1129 isoly(i,j)=.true. 1130 else 1131 frontfacey(i,j)=1 1286 1132 endif 1287 label_max = label_max+1 1288 if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches=',label_max 1289 endif 1290 else ! on est pas sur une tache-------------------------------------------- 1291 table_out_marais(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) 1292 endif !--------------------------------------------------------------------- 1133 end if 1134 end do 1135 !$OMP END DO 1136 !$OMP END PARALLEL 1137 1138 return 1139 end subroutine determin_front 1140 !------------------------------------------------------------------------------ 1141 1142 subroutine determin_front_tof 1143 1144 use module3D_phy, only:front,ice 1145 1146 implicit none 1147 1148 integer :: i,j 1149 integer,dimension(nx,ny) :: nofront ! tableau de travail (points au dela du front) 1150 !$OMP PARALLEL 1151 !$OMP DO 1152 do j=2,ny-1 1153 do i=2,nx-1 1154 1155 if (ice(i,j).eq.1) then ! test si ice=1 1156 1157 ! if ice, on determine front... 1158 ! ainsi, front=0 sur les zones = 0 1159 1160 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 1161 !front= le nb de faces en contact avec un voisin englacé 1162 endif 1163 end do 1164 end do 1165 !$OMP END DO 1166 1167 ! traitement des bords. On considere que l'exterieur n'a pas de glace 1168 ! attention ce n'est vrai que sur la grande grille 1169 1170 !$OMP DO PRIVATE(i) 1171 do j=2,ny-1 1172 i=1 1173 front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) 1174 i=nx 1175 front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1)) 1176 end do 1177 !$OMP END DO 1178 1179 !$OMP DO PRIVATE(j) 1180 do i=2,nx-1 1181 j=1 1182 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)) 1183 j=ny 1184 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1)) 1185 end do 1186 !$OMP END DO 1187 !$OMP BARRIER 1188 ! traitement des coins 1189 1190 front(1,1)=ice(2,1)+ice(2,1) 1191 front(1,ny)=ice(2,ny)+ice(1,ny-1) 1192 front(nx,1)=ice(nx,2)+ice(nx-1,1) 1193 front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny) 1194 1195 ! Les points à plus d'un point de grille du bord sont front=-1 1196 !$OMP WORKSHARE 1197 nofront(:,:)=0 1198 !$OMP END WORKSHARE 1199 !$OMP BARRIER 1200 1201 !$OMP DO 1202 do j=2,ny-1 1203 do i=2,nx-1 1204 if ((ice(i,j).eq.0).and.(front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).eq.0) then 1205 nofront(i,j)=-1 1206 endif 1207 enddo 1293 1208 enddo 1294 enddo 1295 1296 marais(0)=.false. 1297 1298 !$OMP PARALLEL 1299 !$OMP DO 1300 do j=1,ny 1301 do i=1,nx 1302 if (table_out_marais(i,j).ne.0) then 1303 nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1 1304 endif 1305 flot_marais(i,j) = marais(table_out_marais(i,j)) 1209 !$OMP END DO 1210 !$OMP BARRIER 1211 1212 !$OMP WORKSHARE 1213 where (nofront(:,:).eq.-1) 1214 front(:,:)=-1 1215 endwhere 1216 !$OMP END WORKSHARE 1217 !$OMP END PARALLEL 1218 1219 end subroutine determin_front_tof 1220 1221 1222 !> SUBROUTINE: determin_marais 1223 !! afq -- Routine pour l'identification des "marais" 1224 !! Un marais est une tache "shelf" entouré de points grounded 1225 !! Copie sauvage de determin_tache, adapte au probleme du marais 1226 !> 1227 subroutine determin_marais 1228 1229 use module3D_phy, only:ice,flot,flot_marais,debug_3D 1230 !$ USE OMP_LIB 1231 1232 implicit none 1233 1234 integer :: i,j 1235 integer :: indice 1236 integer :: label ! no des taches rencontrées dans le mask 1237 integer :: label_max ! no temporaire maxi de tache rencontrées 1238 integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales 1239 integer,dimension(mask_nb) :: mask ! numero de tache des points adjacents 1240 1241 integer,dimension(nx,ny) :: table_out_marais !< numeros de tache d'un point ij 1242 integer,dimension(0:n_ta_max) :: compt_marais !< contient les equivalence entre les taches 1243 integer,dimension(0:n_ta_max) :: nb_pts_marais !< indique le nombre de points par tache 1244 logical,dimension(0:n_ta_max) :: marais !< T si flottants entoure de poses, F sinon 1245 1246 1247 ! 1-initialisation 1248 !----------------- 1249 label_max=1 ! numero de la tache, la premiere tache est notée 1 1250 label=1 1251 do i=1,n_ta_max 1252 compt_marais(i)=i 1306 1253 enddo 1307 enddo 1308 !$OMP END DO 1309 !$OMP END PARALLEL 1310 1311 debug_3D(:,:,125)=real(table_out_marais(:,:)) 1312 1313 end subroutine determin_marais 1254 !$OMP PARALLEL 1255 !$OMP WORKSHARE 1256 table_out_marais(:,:) = 0 1257 marais(:) = .true. 1258 nb_pts_marais(:) = 0 1259 !$OMP END WORKSHARE 1260 !$OMP END PARALLEL 1261 1262 ! 2-reperage des taches 1263 !---------------------- 1264 1265 do j=2,ny-1 1266 do i=2,nx-1 1267 if ((ice(i,j).ge.1).and.flot(i,j)) then ! on est sur la glace qui flotte-------------------! 1268 1269 if (((ice(i-1,j).ge.1).and.flot(i-1,j)).or.((ice(i,j-1).ge.1).and.flot(i,j-1))) then !masque de 2 cases adjacentes 1270 ! un des voisins est deja en glace 1271 mask(1) = table_out_marais(i-1,j) 1272 mask(2) = table_out_marais(i,j-1) 1273 label = label_max 1274 1275 ! on determine la valeur de la tache minimun (>0) presente ds le masque 1276 do indice=1,mask_nb 1277 if (mask(indice).gt.0) label=min(label,mask(indice)) 1278 enddo 1279 1280 ! on fixe la valeur de la tache voisine minimun au point etudie (via label) 1281 table_out_marais(i,j)=label 1282 1283 !si un des voisins n'est pas glace alors la tache n'est pas un marais 1284 if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then 1285 marais(label)=.false. 1286 endif 1287 1288 ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt_marais 1289 ! on lui affecte le numero de la tache fondamentale 1290 1291 ! exemple on est sur le point X : 5 X 1292 do indice=1,mask_nb ! 20 1293 if(mask(indice).gt.label) then ! mask(2)=20 > 5 1294 compt_marais(mask(indice))=label ! compt_marais(20)=5 1295 if (.not.marais(mask(indice))) marais(label)=.false. ! si la tache n'etais pas un marais => marais =.false. marais(-(-5))=.false. 1296 where (table_out_marais(:,:).eq.mask(indice)) ! where table_out_marais(:,:)=mask(2)=20 1297 table_out_marais(:,:)=label ! table_out_marais(:,:)=label=5 1298 endwhere 1299 endif 1300 enddo 1301 1302 else !aucun des voisins est une tache 1303 table_out_marais(i,j)= label_max 1304 compt_marais(label_max)=label_max 1305 1306 ! si un des voisins n'est pas glace alors la tache n'est pas un marais 1307 if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then 1308 marais(label_max)=.false. 1309 endif 1310 label_max = label_max+1 1311 if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches=',label_max 1312 endif 1313 else ! on est pas sur une tache-------------------------------------------- 1314 table_out_marais(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) 1315 endif !--------------------------------------------------------------------- 1316 enddo 1317 enddo 1318 1319 marais(0)=.false. 1320 1321 !$OMP PARALLEL 1322 !$OMP DO 1323 do j=1,ny 1324 do i=1,nx 1325 if (table_out_marais(i,j).ne.0) then 1326 nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1 1327 endif 1328 flot_marais(i,j) = marais(table_out_marais(i,j)) 1329 enddo 1330 enddo 1331 !$OMP END DO 1332 !$OMP END PARALLEL 1333 1334 debug_3D(:,:,125)=real(table_out_marais(:,:)) 1335 1336 end subroutine determin_marais 1314 1337 1315 1338 end module flottab_mod
Note: See TracChangeset
for help on using the changeset viewer.
