!> \file flottab2-0.7.f90 !! Module pour determiner les endroits ou la glace !! flotte , les iles, et la position du front de glace !< !> \namespace flottab_mod !! Determine les endroits ou la glace !! flotte , les iles, et la position du front de glace !! \author Vincent & Cat !! \date 10 juillet 2005 !! @note Used module !! @note - use module3D_phy !! @note - use module_choix !! @todo "nasty Island". If the bedrock is above sealevel force grounded (mk_init) !< module flottab_mod USE module3D_phy use module_choix IMPLICIT NONE real :: surnet !< surnet hauteur de glace au dessus de la mer real :: archim !< test de flottaison integer:: itestf logical,dimension(nx,ny) :: gz1mx,gz1my logical,dimension(nx,ny) :: fl1mx,fl1my real,dimension(nx,ny) :: uxs1 !< uxbar a l'entree de flottab real,dimension(nx,ny) :: uys1 !< uybar a l'entree de flottab integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y ! Variables pour la determination des differents shelfs/stream ! (representés comme des taches ou l'on resoud l'eq elliptique) !________________________________________________________________ integer,parameter :: n_ta_max=2000!< nombre de tache max integer,dimension(nx,ny) :: table_out !< pour les numeros des taches integer,dimension(nx,ny) :: tablebis !< pour les numeros des taches integer,dimension(0:n_ta_max) :: compt !< contient les equivalence entre les taches integer,dimension(0:n_ta_max) :: nb_pts_tache !< indique le nombre de points par tache logical,dimension(0:n_ta_max) :: iceberg !< T si iceberg, F si calotte posee logical,dimension(nx,ny) :: mask_tache_ij !< masque de travail !< vrai pour toute la tache de i,j integer,dimension(2) :: smax_coord !< pour le maxloc des iles ! Variables pour determiner le point le plus haut (surf) ! d'une ile completement stream !_________________________________________________________ ! icetrim : T si ice stream, F si calotte posee(vertical shear) logical,dimension(n_ta_max) :: icetrim integer :: ii,jj integer :: smax_i integer :: smax_j real :: smax_ integer :: numtache integer :: nb_pt contains ! ----------------------------------------------------------------------------------- !> SUBROUTINE: flottab() !! Cette routine determine les endroits ou la glace !! flotte , les iles, et la position du front de glace !! @note Il y a 4 sortes de zone !! @note - Pose ! @note - Grounding zone et streams gzmx et gzmy ! @note - Iles ilemx, ilemy ! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur !> subroutine flottab() ! ! Vince 5 Jan 95 ! Modifie 20 Jan 95 ! Modifie le 30 Novembre 98 ! Passage f90 + determination des fronts Vincent dec 2003 ! nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005 ! Re-nettoyage par Cat en aout 2006. ! Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging ! pour les points cotiers, toujours fait dans flottab ! ! ----------------- ! ! Cette routine determine les endroits ou la glace ! flotte , les iles, et la position du front de glace ! ! Il y a 4 sortes de zone ! Pose ! Grounding zone et streams gzmx et gzmy ! Iles ilemx, ilemy ! flottant flot sur le noeud majeur, flotmx sur le noeud mineur ! passage dans flottab tous les pas de temps dt ! ! ! _________________________________________________________ if (itracebug.eq.1) call tracebug(' Entree dans routine flottab') SHELFY = .FALSE. ! cas particulier des runs paleo ou on impose un masque grounded if (igrdline.eq.2) then where ( mk_init(:,:).eq.1) ! pose flot(:,:) = .False. H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro) ! pour avoir archim=10 S(:,:) = Bsoc(:,:) + H(:,:) elsewhere flot(:,:) = .True. end where end if ! 1-INITIALISATION ! ---------------- ! initialisation des variables pour detecter les points qui se mettent ! a flotter entre 2 dtt appel_new_flot=.false. do j=1,ny do i=1,nx new_flot_point(i,j)=.false. new_flotmx(i,j)=.false. new_flotmy(i,j)=.false. enddo enddo ! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m ICE(:,:)=0 front(:,:)=0 frontfacex(:,:)=0 frontfacey(:,:)=0 isolx(:,:)=.FALSE. isoly(:,:)=.FALSE. cotemx(:,:)=.false. cotemy(:,:)=.false. boost=.false. ! fin de l'initialisation !_____________________________________________________________________ ! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS ! ------------------------------------- do j=1,ny do i=1,nx uxs1(i,j)=uxbar(i,j) uys1(i,j)=uybar(i,j) archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then ! le point flotte mk(i,j)=1 ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then ! il ne flottait pas avant FLOT(I,J)=.true. BOOST=.false. ! Attention : avec le bloc dessous il faut faire un calcul de flottaison a la lecture if (igrdline.eq.1) then ! en cas de grounding line prescrite flot(i,j)=.false. H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10 new_flot_point(i,j)=.false. endif else ! isynchro=0 ou il flottait déja if (.not.FLOT(I,J)) then ! il ne flottait pas (isynchro=0) new_flot_point(i,j)=.true. ! signale un point qui ne flottait pas ! au pas de temps precedent flot(i,j)=.true. if (igrdline.eq.1) then ! en cas de grounding line prescrite flot(i,j)=.false. H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10 new_flot_point(i,j)=.false. endif endif endif ex_pose else ! le point ne flotte pas mk(i,j)=0 if(FLOT(I,J)) then ! mais il flottait avant FLOT(I,J)=.false. BOOST=.false. endif endif arch ! Si la glace flotte -> PRUDENCE !!! ! S et B sont alors determines avec la condition de flottabilite if (flot(i,j)) then shelfy = .true. surnet=H(i,j)*(1.-ro/row) S(i,j)=surnet+sealevel B(i,j)=S(i,j)-H(i,j) end if end do end do !!$ do i=1,nx !!$ do j=1,ny !!$ if (flot(i,j)) then !!$ mk(i,j)=1 !!$ else !!$ mk(i,j)=0 !!$ endif !!$ end do !!$ end do !----------------------------------------------------------------------- domain_x: do j=1,ny do i=2,nx ! 3_x A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS ! AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX ! ------------------------------------------- ! fl1 est l'ancienne valeur de flotmx gz1mx(i,j)=gzmx(i,j) ! gz1 est l'ancienne valeur de gzmx fl1mx(i,j)=flotmx(i,j) ! fl1 est l'ancienne valeur de flotmx flotmx(i,j)=flot(i,j).and.flot(i-1,j) ! test pour detecter les nouveaux flotmx entre 2 dtt : if (flotmx(i,j).and.(new_flot_point(i,j).or. & new_flot_point(i-1,j))) then appel_new_flot=.true. new_flotmx(i,j)=.true. endif ! premiere determination de gzmx !__________________________________________________________________________ ! gzmx si un des deux voisins est flottant et l'autre posé ! i-1 i gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & ! F P .or.(.not.flot(i,j).and.flot(i-1,j))) ! P F ! A condition d'etre assez proche de la flottaison ! sur le demi noeud condition archim < 100 m archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j) gzmx(i,j)=gzmx(i,j).and.(archim.le.100.) cotemx(i,j)=gzmx(i,j) end do end do domain_x if (itracebug.eq.1) call tracebug(' routine flottab apres domain_x') ! 3_y B- NOUVELLE DEFINITION DE FLOTMY ! -------------------------------- domain_y: do j=2,ny do i=1,nx gz1my(i,j)=gzmy(i,j) ! gz1 est l'ancienne valeur de gzmy fl1my(i,j)=flotmy(i,j) ! fl1 est l'ancienne valeur de flotmy flotmy(i,j)=flot(i,j).and.flot(i,j-1) ! test pour detecter les nouveaux flotmy entre 2 dtt : if (flotmy(i,j).and.(new_flot_point(i,j).or. & new_flot_point(i,j-1))) then appel_new_flot=.true. new_flotmy(i,j)=.true. endif ! premiere determination de gzmy !__________________________________________________________________________ ! gzmy si un des deux voisins est flottant et l'autre posé gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1)) & .or.(.not.flot(i,j).and.flot(i,j-1))) ! A condition d'etre assez proche de la flottaison ! sur le demi noeud condition archim > 100 m archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j) gzmy(i,j)=gzmy(i,j).and.(archim.le.100.) cotemy(i,j)=gzmy(i,j) end do end do domain_y !------------------------------------------------------------------------------------- ! attention : pour expériences Heino ! gzmy(i,j)=gzmy_heino(i,j) ! shelfy=.true. ! shelfy=.false. !____________________________________________________________________________________ !!$ !!$ !!$! 4- Condition sur les bords !!$ !!$ !!$ do i=2,nx !!$ flotmx(i,1) = (flot(i,1).or.flot(i-1,1)) !!$ flotmy(i,1) = .false. !!$ end do !!$ !!$ do j=2,ny !!$ flotmy(1,j) = (flot(1,j).or.flot(1,j-1)) !!$ flotmx(1,j) = .false. !!$ end do !!$ !!$ flotmx(1,1) = .false. !!$ flotmy(1,1) = .false. !!$ ! 4- determination des iles ! ------------------------- ilemx(:,:)=.false. ilemy(:,:)=.false. ! selon x ilesx: do j=2,ny-1 do i=3,nx-2 ! F G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. & (sdx(i,j).LT.1.E-02)) then ilemx(i,j)=.true. ilemx(i+1,j)=.true. ! F G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((flot(i-1,j).and..not.flot(i,j) & .and..not.flot(i+1,j)).and.flot(i+2,j).and. & (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then ilemx(i,j)=.true. ilemx(i+1,j)=.true. ilemx(i+2,j)=.true. ! F G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((flot(i-2,j).and..not.flot(i-1,j) & .and..not.flot(i,j)).and.flot(i+1,j).and. & (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then ilemx(i-1,j)=.true. ilemx(i,j)=.true. ilemx(i+1,j)=.true. ! F G G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((i.lt.nx-2) & .and.(flot(i-2,j).and..not.flot(i-1,j) & .and..not.flot(i,j)).and..not.flot(i+1,j) & .and.flot(i+2,j).and. & (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 & .and.sdx(i+1,j).LT.1.E-02)) then ilemx(i-1,j)=.true. ilemx(i,j)=.true. ilemx(i+1,j)=.true. ilemx(i+2,j)=.true. endif end do end do ilesx ! selon y ilesy: do j=3,ny-2 do i=2,nx-1 ! F G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. & (sdy(i,j).LT.1.E-02)) then ilemy(i,j)=.true. ilemy(i,j+1)=.true. ! F G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((flot(i,j-1).and..not.flot(i,j) & .and..not.flot(i,j+1)).and.flot(i,j+2).and. & (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then ilemy(i,j)=.true. ilemy(i,j+1)=.true. ilemy(i,j+2)=.true. ! F G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((flot(i,j-2).and..not.flot(i,j-1) & .and..not.flot(i,j)).and.flot(i,j+1).and. & (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then ilemy(i,j-1)=.true. ilemy(i,j)=.true. ilemy(i,j+1)=.true. ! F G G G F ! x ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m) else if ((j.lt.ny-2) & .and.(flot(i,j-2).and..not.flot(i,j-1) & .and..not.flot(i,j)).and..not.flot(i,j+1) & .and.flot(i,j+2).and. & (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 & .and.sdy(i,j+1).LT.1.E-02)) then ilemy(i,j-1)=.true. ilemy(i,j)=.true. ilemy(i,j+1)=.true. ilemy(i,j+2)=.true. endif end do end do ilesy ! fin des iles !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans flottab avant dragging asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans flottab avant dragging pas d asymetrie sur H pour time=',time !!$ !!$end if ! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy !---------------------------------------------------------------------------------- if (iter_beta.eq.0) then Call dragging endif if (itracebug.eq.1) call tracebug(' routine flottab apres call dragging') !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans flottab apres dragging asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans flottab aapres dragging pas d asymetrie sur H pour time=',time !!$ !!$end if ! 6- calcule les vitesses des points qui sont devenus gzm do j=1,ny do i=2,nx-1 ! si le point etait posé (non gz) et devient gzmx ! definir la direction de la vitesse (moyenne des points) if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. & (i.gt.2).and.(i.lt.nx)) then uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2. endif end do end do do j=2,ny-1 do i=1,nx ! si le point etait posé (non gz) et devient gzmy ! definir la direction de la vitesse (moyenne des points) if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. & (j.gt.2).and.(j.lt.ny)) then uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2. endif end do end do ! 7-On determine finalement la position des noeuds stream ou shelf ! ------------------------------------------------------------- if (nt.ge.2) then ! pour ne pas faire ce calcul lors du premier passage uxbar(:,:)=uxs1(:,:) uybar(:,:)=uys1(:,:) endif flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:))) & .or.(.not.marine.and.flotmx(:,:)) flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:))) & .or.(.not.marine.and.flotmy(:,:)) ! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line !--------------------------------------------------------------------------------- ! fbm est vrai si le point est flottant mais un des voisins est pose !_________________________________________________________________________ do j=2,ny-1 do i=2,nx-1 ! if (i.gt.2.AND.i.lt.nx) then fbm(i,j)=flot(i,j).and. & ((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) & .or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1))) ! endif end do end do ! 9-On determine maintenant la position du front de glace ! ------------------------------------------------------- ! C'est ici que l'on determine la position et le type de front ! print*,'on est dans flottab pour definir les fronts' !!$do i=3,nx-2 !!$ do j=3,ny-2 !!$ if (h(i,j).gt.1.1) ice(i,j)=1 !!$ end do !!$end do where (flot(:,:)) where (H(:,:).gt.(1.1)) ice(:,:)=1 elsewhere ice(:,:)=0 end where elsewhere where (H(:,:).gt.(.1)) ice(:,:)=1 elsewhere ice(:,:)=0 end where end where call DETERMIN_TACHE synchro : if (isynchro.eq.1) then !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans flottab avant DETERMIN_TACHE asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H pour time=',time !!$ !!$end if !----------------------------------------------! !On determine les differents ice strean/shelf ! call DETERMIN_TACHE ! !----------------------------------------------! !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans flottab apres DETERMIN_TACHE asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H pour time=',time !!$ !!$end if ! / test aurel debug_3d(:,:,43)=table_out(:,:) iceberg(:)=.true. do j=1,ny do i=1,nx if(.not.flot(i,j)) then iceberg(table_out(i,j))=.false. end if end do end do do j=1,ny do i=1,nx debug_3d(i,j,56)=nb_pts_tache(table_out(i,j)) if (iceberg(table_out(i,j))) then debug_3d(i,j,57)=1 else debug_3d(i,j,57)=0 end if end do end do ! fin test aurel / !ice(:,:)=0 Attention, voir si ca marche toujours pour l'Antarctique et heminord ! !On compte comme englacé uniquement les calottes dont une partie est posée do i=3,nx-2 do j=3,ny-2 test1: if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg if (nb_pts_tache(table_out(i,j)).ge.1) then ice(i,j)=1 if (nb_pts_tache(table_out(i,j)).le.10) then ! les petites iles sont en sia ! write(6,*) 'petite ile ',i,j flgzmx(i,j)=.false. flgzmx(i+1,j)=.false. flgzmy(i,j)=.false. flgzmy(i,j+1)=.false. gzmx(i:i+1,j)=.false. gzmy(i,j:j+1)=.false. endif ! ici on est sur une tache non iceberg >= 5 points ! on teste si la tache n'est pas completement ice stream test2: if (icetrim(table_out(i,j))) then ! on a une ile d'ice stream mask_tache_ij(:,:)=.false. mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j)) ! pour toute la tache smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) smax_i=smax_coord(1) smax_j=smax_coord(2) !!$ smax_i=0 ; smax_j=0 ; smax_=sealevel !!$ do ii=3,nx-2 !!$ do jj=3,ny-2 !!$ if (table_out(ii,jj).eq.table_out(i,j)) then !!$ if (s(ii,jj).gt.smax_) then !!$ smax_ =s(ii,jj) !!$ smax_i=ii !!$ smax_j=jj !!$ endif !!$ endif !!$ end do !!$ end do gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. if (Smax_.le.sealevel) then write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' write(num_tracebug,*)'time=',time,' coord:',smax_i,smax_j end if endif test2 end if ! endif deplace else ! on est sur un iceberg ! test1 ice(i,j)=0 h(i,j)=1. surnet=H(i,j)*(1.-ro/row) S(i,j)=surnet+sealevel B(i,j)=S(i,j)-H(i,j) endif test1 end do end do !---------------------------------------------- ! On caracterise le front des ice shelfs/streams ! call DETERMIN_FRONT !---------------------------------------------- !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans flottab apres DETERMIN_front asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H pour time=',time !!$ !!$end if endif synchro ! correction momentanée pour symetrie Heino !where ((.not.flot(:,:)).and.(ice(:,:).eq.0)) H(:,:)=0. !fin de routine flottab2 !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j) end subroutine flottab !-------------------------------------------------------------------- !> SUBROUTINE: determin_tache !! Routine pour la dtermination du numero de tache a effectuer !> subroutine determin_tache implicit none integer :: indice integer :: label ! no des taches rencontrées dans le mask integer :: label_max ! no temporaire maxi de tache rencontrées ! integer :: mask_nb = 4 integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales integer :: vartemp ! variable temporaire pour reorganiser compt ! integer,dimension(mask_nb) :: mask integer,dimension(mask_nb) :: mask ! 1-initialisation !----------------- label_max=1 ! numero de la tache, la premiere tache est notée 1 label=1 do i=1,n_ta_max compt(i)=i enddo ! table_in = .false. table_out(:,:) = 0 iceberg(:) = .true. icetrim (:) = .true. nb_pts_tache(:) = 0 ! open(unit=100,file="tache.data",status='replace') ! 2-reperage des taches !---------------------- do i=2,nx-1 do j=2,ny-1 IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------! if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes ! un des voisins est deja en glace mask(1) = table_out(i-1,j) mask(2) = table_out(i,j-1) label = label_max !on determine la valeur de la tache minimun (>0) presente ds le masque do indice=1,mask_nb if (mask(indice).gt.0) label=min(label,mask(indice)) enddo !on fixe la valeur de la tache voisine minimun au point etudie (via label) table_out(i,j)=label !si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F. if (.not.FLOT(I,J)) then iceberg(label)=.false. endif !si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F. if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then icetrim(label)=.false. endif ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt ! i.e. les plus grands numeros correspondent au plus petit ! on lui affecte le numero de la tache fondamentale avec un signe - ! pour indiquer le changement do indice=1,mask_nb if(mask(indice).gt.label) then compt(mask(indice))=-label endif enddo else !aucun des voisins est une tache table_out(i,j)= label_max compt(label_max)=label_max !si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F. if (.not.FLOT(I,J)) then iceberg(label_max)=.false. endif !si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F. if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. & (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then icetrim(label)=.false. endif label_max = label_max+1 if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max endif else !on est pas sur une tache---------------------------------------------- table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) endif !--------------------------------------------------------------------- enddo enddo ! On reorganise compt en ecrivant le numero de la tache fondamentale ! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité) ! On indique aussi le nb de point que contient chaque taches (nb_pts_tache) do indice=1,label_max vartemp = compt(indice) if (compt(indice).lt.0) then compt(indice)= compt(-vartemp) if (.not.iceberg(indice)) iceberg(-vartemp)=.false. if (.not.icetrim(indice)) icetrim(-vartemp)=.false. endif enddo do i=1,nx do j=1,ny if (table_out(i,j).ne.0) then table_out(i,j)=compt(table_out(i,j)) nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1 endif enddo enddo !!$tablebis(:,:)=table_out(:,:) !!$do j=1,ny !!$ do i=1,nx !!$ if (tablebis(i,j).ne.0) then ! tache de glace !!$ numtache=table_out(i,j) !!$ nb_pt=count(table_out(:,:).eq.numtache) ! compte tous les points de la tache !!$ nb_pts_tache(table_out(i,j))=nb_pt ! !!$ !!$ where (table_out(:,:).eq.numtache) !!$ tablebis(:,:)=0 ! la table de tache est remise a 0 pour eviter de repasser !!$ end where !!$ write(6,*) i,j,nb_pt,table_out(i,j) !!$ endif !!$ enddo !!$enddo !!$do j=1,ny !!$ do i=1,nx !!$ debug_3d(i,j,56)=nb_pts_tache(table_out(i,j)) !!$ end do !!$end do !!$ end subroutine determin_tache !---------------------------------------------------------------------- !> SUBROUTINE: determin_front !!Routine pour la determination du front !> subroutine determin_front integer :: i_moins1,i_plus1,i_plus2 integer :: j_moins1,j_plus1,j_plus2 do i=3,nx-2 do j=3,ny-2 surice:if (ice(i,j).eq.0) then do pmx=-1,1,2 do pmy=-1,1,2 diagice : if (ice(i+pmx,j+pmy).eq.1) then if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues ! de glaces diagonales en coin(26dec04) if ((ice(i+2*pmx,j).eq.1.and.ice(i+2*pmx,j+pmy).eq.0).or.& (ice(i,j+2*pmy).eq.1.and.ice(i+pmx,j+2*pmy).eq.0)) then ice(i,j)=1 h(i,j)=max(1.,h(i,j)) endif ! test (i) pour eviter les langues de glaces diagonales : ! mouvement du cheval aux echecs if ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.1)) then if (ice(i+2*pmx,j+pmy).eq.1.and. & (ice(i+2*pmx,j+2*pmy)+ice(i,j+2*pmy)).ge.1) then ice(i,j)=1 h(i,j)=max(1.,h(i,j)) endif if (ice(i+pmx,j+2*pmy).eq.1.and. & (ice(i+2*pmx,j+2*pmy)+ice(i+2*pmx,j)).ge.1) then ice(i,j)=1 h(i,j)=max(1.,h(i,j)) endif ! test (ii) pour eviter les langues de glaces diagonales : ! le point glace ice(i+pmx,j+pmy) a : ! - ses 4 voisins frontaux en glace ! - mais 2 voisins vides diagonalement opposes elseif ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.2) & .and.ice(i+2*pmx,j+2*pmy).eq.0) then ! test (iii) pour faire les tests (i) et (ii) ice(i,j)=1 h(i,j)=max(1.,h(i,j)) ice(i+2*pmx,j+2*pmy)=1 h(i+2*pmx,j+2*pmy)=max(1.,h(i+2*pmx,j+2*pmy)) endif endif endif diagice enddo enddo endif surice end do end do !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans front avant remplissage baies asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H pour time=',time !!$ !!$end if ! print*,'dans remplissage baies',time baies: do k=1,2 do j=1,ny do i=1,nx surice_xy: if (ice(i,j).eq.0) then i_moins1=max(i-1,1) j_moins1=max(j-1,1) i_plus1=min(i+1,nx) j_plus1=min(j+1,ny) i_plus2=min(i+2,nx) j_plus2=min(j+2,ny) ! test (iii) pour trouver les baies de largeur 1 ou 2 cases ! et combler les trous si ce sont des baies ! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles ! baies horizontales if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1 ice(i,j)=1 H(i,j)=max(1.,H(i,j)) endif endif if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1 ice(i,j)=1 H(i,j)=max(1.,H(i,j)) endif endif endif surice_xy end do end do end do baies !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans front apres remplissage baies asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H pour time=',time !!$ !!$end if do i=2,nx-1 do j=2,ny-1 if (ice(i,j).eq.1) then ! test si ice=1 ! if ice, on determine front... ! ainsi, front=0 sur les zones = 0 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) !front= le nb de faces en contact avec un voisin englacé endif end do end do ! traitement des bords. On considere que l'exterieur n'a pas de glace ! attention ce n'est vrai que sur la grande grille do j=2,ny-1 i=1 front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1)) i=nx front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1)) end do do i=2,nx-1 j=1 front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)) j=ny front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1)) end do ! traitement des coins front(1,1)=ice(2,1)+ice(2,1) front(1,ny)=ice(2,ny)+ice(1,ny-1) front(nx,1)=ice(nx,2)+ice(nx-1,1) front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny) !!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf) !!$if (itestf.gt.0) then !!$ write(6,*) 'dans front apres front asymetrie sur H pour time=',time !!$ stop !!$else !!$ write(6,*) 'dans front apres front pas d asymetrie sur H pour time=',time !!$ !!$end if ! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales) ! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs ! de n'importe quelle taille). ! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites do j=1,ny do i=1,nx-1 if (front(i,j).eq.1) then if (front(i+1,j).eq.1) then ice(i,j)=0 ice(i+1,j)=0 front(i,j)=0 front(i+1,j)=0 endif endif end do end do do j=1,ny-1 do i=1,nx if (front(i,j).eq.1) then if (front(i,j+1).eq.1) then ice(i,j)=0 ice(i,j+1)=0 front(i,j)=0 front(i,j+1)=0 endif end if end do end do !isolx signifie pas de voisins en x !isoly signifie pas de voisins en y !remarque : !si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0) ! calcul de frontfacex et isolx do j=1,ny do i=2,nx-1 if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x if ((ice(i-1,j)+ice(i+1,j)).eq.0) then isolx(i,j)=.true. elseif (ice(i-1,j).eq.0) then frontfacex(i,j)=-1 ! front i-1 |i i+1 else frontfacex(i,j)=+1 ! front i-1 i| i+1 endif endif end if !fin du test il y a un front end do end do ! calcul de frontfacey et isoly do j=2,ny-1 do i=1,nx if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3) if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y if ((ice(i,j-1)+ice(i,j+1)).eq.0) then isoly(i,j)=.true. !front j-1 |j| j+1 elseif (ice(i,j-1).eq.0) then frontfacey(i,j)=-1 !front j-1 |j j+1 else frontfacey(i,j)=+1 !front j-1 j| j+1 endif endif end if !fin du test il y a un front end do end do ! traitement des bords. On considere que l'exterieur n'a pas de glace ! attention ce n'est vrai que sur la grande grille do j=2,ny-1 i=1 if (front(i,j).ge.1) then if (ice(i+1,j).eq.0) then isolx(i,j)=.true. else frontfacex(i,j)=-1 endif end if i=nx if (front(i,j).ge.1) then if (ice(i-1,j).eq.0) then isolx(i,j)=.true. else frontfacex(i,j)=1 endif end if end do do i=2,nx-1 j=1 if (front(i,j).ge.1) then if (ice(i,j+1).eq.0) then isoly(i,j)=.true. else frontfacey(i,j)=-1 endif end if j=ny if (front(i,j).ge.1) then if (ice(i,j-1).eq.0) then isoly(i,j)=.true. else frontfacey(i,j)=1 endif end if end do return end subroutine determin_front !------------------------------------------------------------------------------ end module flottab_mod