!> \file remplimat-ant-0.5-40km.f90 !! TOOOOOOO DOOOOOO !< !====================================================================! SUBROUTINE REMPLIDOM(N1,N2,uxprec,uyprec,uxnew,uynew,ifail) ! ! ! ! Vince 17 Jan ! ! Calcul de la vitesse ! ! des coins -> 25 fev ! ! ! ! Modification du schema numerique 3 Mars 95 (stabilite) ! ! Aout 95 (optimisation) ! ! Modifications Catherine Mars 97 ! ! gradient conjugue utilisant blas ! ! ! ! Modification catherine Avril 97 ! ! DECOMPOSITION DE DOMAINE ! ! Catherine Mai 97 application a l\'Antarctique ! ! ! ! Reecriture Vincent novembre 2003 ! ! On remanie les conditions aux limites ! ! ------------------------------------------------------ ! ! numdom numero du domaine traite ! ! id1,id2 bornes selon x ! ! jd1,jd2 bornes selon y ! ! ! ! uxprex(n1,n2) ! ! uyprec(n1,n2) vitesses de l\'iteration precedente ! ! ! ! uxnew(n1,n2) ! ! uynew(n1,n2) uynew resultat de cette iteration ! ! ! !====================================================================! USE module3D_phy USE EQ_ELLIPTIQUE_MOD IMPLICIT NONE ! dimensionnement momentane de slv real, dimension(nx,ny) :: slv ! 1_Declaration des variables ! ---------------------------- integer nblig integer aaai,aaaj real aaa1 integer iinf,isup,jinf,jsup integer iii,jjj ! position du noeud de l'echec de sgbsv real aaa INTEGER LDB,NRHS,LDAB,N1,N2 INTEGER kzero,knonul,ITER LOGICAL CONJGR!,IREORD INTEGER nb_err ! nombre de passage dans sgbsv ! apres le 3eme passage infructueux on stoppoe !!!!!! INTEGER IPIV(NPTMAX) INTEGER CTEURLIG,SNR,INEARX INTEGER IFLOTMX(NX,NY),IFLOTMY(NX,NY),ifail ! dimensions liees au domaine ! INTEGER NUMDOM,ID1,ID2,JD1,JD2 REAL UXPREC(NX,NY),UYPREC(NX,NY),UXNEW(NX,NY),UYNEW(NX,NY) ! !Pour le remplissage de la matice REAL SCAL INTEGER IIMAT !INTEGER,DIMENSION(NX,NY) :: INEARX ! la simple variable est devenu ! tableau !!!!!!!! voir si on peut la remplacer par un tableau deja !!!!!!!!!existant car elle est utile que temporairement !!!!!!!!! dans §2_Initialisation des variables REAL,DIMENSION(NX,NY) :: TUMIMJ, TUMIJ, TUMIPJ, & TUIMJ , TUIJ , TUIPJ , & TUPIMJ, TUPIJ, TUPIPJ, & TVMIJ,TVMIPJ,TVIJ,TVIPJ !Pour le calcul des U REAL,DIMENSION(NX,NY) :: SVMIMJ, SVMIJ, SVMIPJ, & SVIMJ , SVIJ , SVIPJ , & SVPIMJ, SVPIJ, SVPIPJ, & SUIMJ,SUIJ,SUPIMJ,SUPIJ !Pour le calcul des V REAL,DIMENSION(NX,NY) :: OPPOSX,OPPOSY!Pression laterale contrebalancant le mouvement pour u et v logical,DIMENSION(NX,NY) ::okumatest,okvmatest ! test de controle pour ! ligu, ligv transporte depuis 3D integer, dimension(NX,NY) :: LIGU ! numero de ligne de U dans remplidom integer, dimension(NX,NY) :: LIGV ! numero de ligne de V dans remplidom real upr,unw, emmax real moteur,beta!,unorm INTEGER,DIMENSION(:,:),allocatable:: posligu, posligv !transpose de ligu,ligv INTEGER,DIMENSION(:),allocatable:: ipiv !7eme parametre de sgbsv (pour permutations) CHARACTER(len=8),DIMENSION(NX,NY,2) :: METHOD INTEGER,DIMENSION(NX,NY,2) :: METHOD2 integer :: num_err = 266 ! sorties debug open(num_err,file='erreurs-remplimat') write(num_err,*)' time=',time uxnew(:,:)=0. uynew(:,:)=0. ! attention, pour l'instant slv est local et juste attribue ici slv(:,:)=sealevel kzero =0 nb_err=0 ! 2_Initialisation des variables ! ---------------------------------- METHOD(:,:,:)='000000' METHOD2(:,:,:)=-2 NRHS=1 LDAB=2*KL+KU+1 LDB=NPTMAX INEARX=1 CTEURLIG = 0 iii=0; jjj=0 LIGU(:,:)=0 LIGV(:,:)=0 OKUMAT(:,:)=.FALSE. OKVMAT(:,:)=.FALSE. OKUMATEST(:,:)=.FALSE.! TEST OKVMATEST(:,:)=.FALSE. ! attribution de iflotmx et iflotmy pour remplacer les operations ! sur les logiques, iflot est vrai quand soit flotmx ou gzmx est vrai (soit flgzmx) do J=1,NY do I=1,NX if (FLGZMX(I,J)) then IFLOTMX(I,J)=1 FROTMX(I,J)=min(abs(betamx(i,j)),betamax) else IFLOTMX(I,J)=0. FROTMX(I,J)=abs(betamx(i,j)) endif if (FLGZMY(I,J)) then IFLOTMY(I,J)=1 FROTMY(I,J)= min(abs(betamy(i,j)),betamax) else IFLOTMY(I,J)=0 FROTMY(I,J)=abs(betamy(i,j)) endif end do end do !testvincent ! print*,'emplimat',time,sealevel,tafor !si le point ou un des voisins et stream/shelf inearx.ge.1 !!!!!DO I=ID1,ID2 !!!!!!!!DO J=JD1,JD2 DO j=2,ny !J=1,NY DO i=2,nx! I=1,NX !!! Faire des tests sur les shifts pour amelioration code !!! plutot somme de sous-tableaux IF (ICE(I,J).EQ.1) THEN INEARX=IFLOTMX(I-1,J) +IFLOTMX(I,J) +IFLOTMX(I+1,J) & +IFLOTMX(I-1,J-1)+IFLOTMX(I,J-1)+IFLOTMX(I+1,J-1) & +IFLOTMX(I-1,J+1)+IFLOTMX(I,J+1)+IFLOTMX(I+1,J+1) & +IFLOTMY(I-1,J) +IFLOTMY(I,J) +IFLOTMY(I+1,J) & +IFLOTMY(I-1,J-1)+IFLOTMY(I,J-1)+IFLOTMY(I+1,J-1) & +IFLOTMY(I-1,J+1)+IFLOTMY(I,J+1)+IFLOTMY(I+1,J+1) ELSE INEARX=0 END IF !print*,('i=',i,'j=',j,'INEARX=',INEARX) IF (INEARX.ge.1) THEN OKUMAT(I,J) = .TRUE. OKUMAT(I+1,J) = .TRUE. ! redondant ??? OKVMAT(I,J) = .TRUE. OKVMAT(I,J+1) = .TRUE. ! attention a la methode de calcul, les points Ui+1j sont ejectes aux ! calandes ENDIF END DO END DO DO I=1,NX DO J=1,NY IF (OKUMAT(I,J)) THEN cteurlig=cteurlig+1 LIGU(i,j)=cteurlig ENDIF IF (OKVMAT(I,J)) THEN cteurlig=cteurlig+1 LIGV(i,j)=cteurlig ENDIF END DO END DO ! nblig est le nombre d'equation lineaire a resoudre dans mmat NBLIG = CTEURLIG !!!!!!Coeff pour mmat et sgbsv !!!!!! ! print*,'avant largeur de bande, ku=',ku,'kl=',kl call larg_bande(ku,kl) ! print*,'apres largeur de bande, ku=',ku,'kl=',kl if (ku.lt.0.or.kl.lt.0) then print*,'ku=',ku,'kl=',kl print*,'pas de shelves a traiter : sortie de remplimat' ifail=0 return endif kdc2=ku+kl+1 ! print*,'kdc2=',kdc2 ldab=2*kl+ku+1 ! print*,'ldab=',ldab ldb=nblig !test sur l'ancienne taille de mmat ! si mmat est trop petite, on l'agrandie a la taille nessessaire if (nblig.gt.nptmax.or.ldab.gt.LDBMAX) then !print*,'redefinition des matrices',nblig,nptmax,ldab,LDBMAX nptmax=max(nptmax,nblig) LDBMAX=max(LDBMAX,ldab) call DEFINITION_MATRICE(.FALSE.,nptmax,ldbmax) endif allocate(POSLIGU(nblig,2),stat=err) if (err/=0) then print *,"Erreur à l'allocation de POSLIGU",err stop 4 end if allocate(POSLIGV(nblig,2),stat=err) if (err/=0) then print *,"Erreur à l'allocation de POSLIGV",err stop 4 end if allocate(IPIV(nblig),stat=err) if (err/=0) then print *,"Erreur à l'allocation de IPIV",err stop 4 end if cteurlig=0 POSLIGU(:,:)=0 POSLIGV(:,:)=0 DO I=1,NX DO J=1,NY IF (OKUMAT(I,J)) THEN cteurlig=cteurlig+1 LIGU(i,j)=cteurlig POSLIGU(cteurlig,1)=i POSLIGU(cteurlig,2)=j ENDIF IF (OKVMAT(I,J)) THEN cteurlig=cteurlig+1 LIGV(i,j)=cteurlig POSLIGV(cteurlig,1)=i POSLIGV(cteurlig,2)=j ENDIF END DO END DO ! do i=1,cteurlig ! print *,'u',POSLIGU(i,1),POSLIGU(i,2),'v',POSLIGV(i,1),POSLIGV(i,2) ! enddo ! nblig est le nombre d'equation lineaire a resoudre dans mmat if (NBLIG.ne.CTEURLIG) then print *,"Erreur lors du second comptage" stop endif 942 continue MMAT(:,:)=0 BDR(:,:)=0 BDRO(:)=0 !initialisation des tuij,suij... a 0. TUMIMJ=0.; TUMIJ=0.; TUMIPJ=0. TUIMJ=0. ; TUIJ=0. ; TUIPJ=0. TUPIMJ=0.; TUPIJ=0.; TUPIPJ=0. TVMIJ=0.;TVMIPJ=0.;TVIJ=0.;TVIPJ=0. SVMIMJ=0.; SVMIJ=0.; SVMIPJ=0. SVIMJ=0. ; SVIJ=0. ; SVIPJ=0. SVPIMJ=0.; SVPIJ=0.; SVPIPJ=0. SUIMJ=0.;SUIJ=0.;SUPIMJ=0.;SUPIJ=0. !print*,'NBLIG =',NBLIG ! Initialisations et reperages ! ---------------------------- ! IREORD=.FALSE. ! T pour faire le rearrangement de la matrice ! 0 pour aller le lire dans un fichier ! F pour traiter sans rearrangement ! On defini une valeur par defaut de kl et ku, ! 3_Boucle de determination des coefficients de la matrice ! ------------------------------------------------------ !!!!!!!!!! ATTENTION, le sens de l'enroulement pourrait etre modifié DO I=1,NX DO J=1,NY !!!!!!!!!!!!!!!!!!!!! !!!!!TEST POUR Uij!!! !!!!!!!!!!!!!!!!!!!!! IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY)) THEN TUIJ(I,J)=1 SVIJ(I,J)=1 METHOD(i,j,1)='bordx' METHOD2(i,j,1)=-1 METHOD(i,j,2)='bordy' METHOD2(i,j,2)=-1 ELSE IF (.NOT.FLGZMX(I,J)) THEN !!Test sur FLGZMX ! points poses en x et n'etant pas sur le bord du domaine TUIJ(I,J)=1 OPPOSX(i,j)=UXBAR(I,J) if (okumat(i,j)) then METHOD(i,j,1)='poseX' METHOD2(i,j,1)=0 else METHOD(i,j,1)='rienX' METHOD2(i,j,1)=-1 endif ELSE !!!Test sur FLGZMX, ici les points sont stream/shelf !!on resout l'equation elliptique !--------------------------------------------------- IF (ICE(i,j)==1) THEN!!!!!TEST SUR ICE, ici en glace SELECT CASE (FRONT(i,j)) CASE (4) !le point a tous ses voisins => Calcul general ! CALL U_GENERAL METHOD(i,j,1)='4U_GEN' METHOD2(i,j,1)=1 ! CASE (3) !le point a 3 voisins IF (frontfacex(i,j)==0) THEN !!!!!!!!!!!!!!!!!!!!!!!front //Y if (ICE(i-1,j+frontfacey(i,j))==1) then CALL U_GENERAL METHOD(i,j,1)='3U_GEN' METHOD2(i,j,1)=1 else CALL U_CISAILL(frontfacey(i,j)) METHOD(i,j,1)='3U_CIS' METHOD2(i,j,1)=3 endif ELSE !!!!!!!!!!!!!!!!!!!!!!!front //X if (frontfacex(i,j)==-1) then !!front i-1 | i (a droite) CALL U_PRESSIO(-1) METHOD(i,j,1)='3U_P-1' METHOD2(i,j,1)=2 else ! CALL U_PRESSIO(+1) CALL U_GENERAL METHOD(i,j,1)='3U_GEN' METHOD2(i,j,1)=1 endif ENDIF ! ! CASE (2) !le point a 2 voisins IF (ISOLX(i,j)) THEN !langue de glace non confinée avec bords //x || CALL U_TONGUE_LE(frontfacey(i,j)) METHOD(i,j,1)='2U_T_LE' METHOD2(i,j,1)=6 ELSEIF (ISOLY(i,j)) THEN !langue de glace non confinée avec bords //y = CALL U_TONGUE METHOD(i,j,1)='2U_T' METHOD2(i,j,1)=5 ELSE ! ON A 2 FRONTS NON-OPPOSES if (frontfacex(i,j)==-1) then CALL U_COIN(-1) METHOD(i,j,1)='2U_C-1' METHOD2(i,j,1)=4 else !!! CALL U_PRESSIO(+1)!pour Ui+1j!!!!!CALL U_COIN IF((FRONT(i-1,j)==4).or.((FRONT(i-1,j)==3).and. & (frontfacex(i-1,j)==-1)))then!!test pour Uij CALL U_GENERAL METHOD(i,j,1)='2U_GEN' METHOD2(i,j,1)=1 ELSE CALL U_CISAILL(frontfacey(i,j)) METHOD(i,j,1)='2U_CIS' METHOD2(i,j,1)=3 ENDIF endif ENDIF ! ! CASE (1) IF (ISOLX(i,j)) THEN !langue de glace non confinée avec bords //x || CALL U_TONGUE_LE(frontfacey(i,j)) METHOD(i,j,1)='1U_T_LE' METHOD2(i,j,1)=6 ELSE !langue de glace non confinée avec bords //y = if (.not.ISOLY(i,j)) then!test write(num_err,*) 'test sur u case=1 pb'! STOP !test endif !test if (frontfacex(i,j)==1) then CALL U_TONGUE METHOD(i,j,1)='1U_T' METHOD2(i,j,1)=5 else CALL U_TONGUEFRONT(-1) METHOD(i,j,1)='1U_TF-1' METHOD2(i,j,1)=5 endif ENDIF CASE (0) METHOD(i,j,1)='0U_' METHOD2(i,j,1)=9 TUIJ(i,j)=1 OPPOSX(i,j)=5! ! PRINT*,'front =0','i=',i,'j=',j CASE DEFAULT write(num_err,*) 'ERREUR front est pas 0,1, 2, 3, 4' STOP END SELECT ELSE!!!!!TEST SUR ICE, ij pas glace IF (ISOLX(i-1,j)) THEN CALL U_TONGUE_RI(frontfacey(i-1,j)) METHOD(i,j,1)='U_T_RI' METHOD2(i,j,1)=6 ELSEIF (ISOLY(i-1,j)) THEN CALL U_TONGUEFRONT(+1) METHOD(i,j,1)='U_TF+1' METHOD2(i,j,1)=5 ELSEIF (FRONT(i-1,j)==3) THEN CALL U_PRESSIO(+1) METHOD(i,j,1)='U_P+1' METHOD2(i,j,1)=2 ! TUIJ(i,j)=1 ! opposx(i,j)=2 ! print*, 'U_PRESSIO(+1) evite',i,j ELSEIF (FRONT(i-1,j)==2) THEN ! if (ifail.gt.0.and.i==iii.and.j==jjj) then ! CALL U_PRESSIO(+1) ! METHOD(i,j,1)='U_P+1-' ! print*,'modif',i,j,METHOD(i,j,1) if (ifail.gt.0.and.METHOD(i,j,1).eq.'U-chang') then CALL U_PRESSIO(+1) METHOD2(i,j,1)=10 ! METHOD(i,j,1)='U_P+1-' write(num_err,*)'modif',i,j,METHOD(i,j,1) else CALL U_COIN(+1) METHOD2(i,j,1)=4 METHOD(i,j,1)='U_C+1' endif ELSE METHOD(i,j,1)='U_ORFAN' METHOD2(i,j,1)=999 TUIJ(i,j)=1 OPPOSX(i,j)=0!777666 ! print*, 'point flgzmX sans glace', i,j ENDIF ENDIF!!!!!TEST SUR ICE, fin du test ENDIF !!Test sur flgzmx !!!!!!!!!!!!!!!!!!!!! !!!!!TEST POUR Vij!!! !!!!!!!!!!!!!!!!!!!!! IF (.NOT.FLGZMY(I,J)) THEN !!Test sur FLGZMX/Y SVIJ(I,J)=1 OPPOSY(i,j)=UYBAR(I,J) !testvincent if (okvmat(i,j)) then METHOD(i,j,2)='poseY' METHOD2(i,j,2)=0 else METHOD(i,j,2)='rienY' METHOD2(i,j,2)=-1 endif ELSE !!!Test sur FLGZMY, ici les points sont stream/shelf !! on resout l'equation elliptique ! ---------------------------------------------------- IF (ICE(i,j)==1) THEN !test sur ICE, ici point en glace SELECT CASE (FRONT(i,j)) CASE (4) !le point a tous ses voisins => calcul general ! CALL V_GENERAL METHOD(i,j,2)='4V_GE' METHOD2(i,j,2)=1 ! CASE (3) !le point a 3 voisins ! IF (frontfacex(i,j)==0) THEN !!!!!!!!!!!!!!!!!!!!!!!front //Y if (frontfacey(i,j)==-1) then !!front j-1 | j en aval CALL V_PRESSIO(-1) METHOD(i,j,2)='3V_P-1' METHOD2(i,j,2)=2 else !!front j | j+1 en amont ! CALL V_PRESSIO(+1) CALL V_GENERAL METHOD(i,j,2)='3V_GE' METHOD2(i,j,2)=1 endif ELSE !!!!!!!!!!!!!!!!!!!!!!!front //X !test pour le calcul de Vij if(ICE(i+frontfacex(i,j),j-1)==1) then CALL V_GENERAL METHOD(i,j,2)='3V_GEN' METHOD2(i,j,2)=1 else CALL V_CISAILL(frontfacex(i,j)) METHOD(i,j,2)='3V_CIS' METHOD2(i,j,2)=3 endif ENDIF ! ! CASE (2) ! IF (ISOLX(i,j)) THEN !langue de glace non confinée avec bords //x || CALL V_TONGUE METHOD(i,j,2)='2V_T' METHOD2(i,j,2)=5 ELSEIF (ISOLY(i,j)) THEN !langue de glace non confinée avec bords //y = CALL V_TONGUE_DO(frontfacex(i,j)) METHOD(i,j,2)='2V_T_DO' METHOD2(i,j,2)=6 ELSE ! ON A 2 FRONTS NON-OPPOSES if (frontfacey(i,j)==-1) then CALL V_COIN(-1) METHOD(i,j,2)='2V_C-1' METHOD2(i,j,2)=4 else !!!!!!!!!!!!! CALL V_PRESSIO(+1)!!!!!!!!!!!!!!!!CALL V_COIN !TESTS POUR VOIR DE QUEL TYPE EST VIJ IF((FRONT(i,j-1)==4).or.((FRONT(i,j-1)==3).and. & (frontfacey(i,j-1)==-1)))then CALL V_GENERAL METHOD(i,j,2)='2V_GE' METHOD2(i,j,2)=1 ELSE CALL V_CISAILL(frontfacex(i,j)) METHOD(i,j,2)='2V_CIS' METHOD2(i,j,2)=3 ENDIF endif ENDIF ! CASE (1) IF (ISOLX(i,j)) THEN !!!!!!!!!!!!!!!!!!! if (frontfacex(i,j).eq.1) then if (frontfacey(i,j).eq.1) then CALL V_TONGUE METHOD(i,j,2)='1V_T' METHOD2(i,j,2)=5 else CALL V_TONGUEFRONT(-1) METHOD(i,j,2)='1V_TF-1' METHOD2(i,j,2)=5 endif ELSEIF (ISOLY(i,j)) THEN CALL V_TONGUE_DO(frontfacex(i,j)) METHOD(i,j,2)='1V_T_DO' METHOD2(i,j,2)=6 ELSE write(num_err,*) 'ERREUR test sur case 1 ds remplimat' STOP ENDIF CASE (0) METHOD(i,j,2)='0V_' METHOD2(i,j,2)=9 SVIJ(i,j)=1 OPPOSY(i,j)=5! CASE DEFAULT write(num_err,*) 'ERREUR front est pas 0,1, 2, 3, 4' STOP END SELECT ELSE !!!!test sur ice(i,j) point non en glace if (isolx(i,j-1)) then CALL V_TONGUEFRONT(+1) METHOD(i,j,2)='V_TF+1' METHOD2(i,j,2)=5 elseif (isoly(i,j-1)) then CALL V_TONGUE_UP(frontfacex(i,j-1)) METHOD(i,j,2)='V_T_UP' METHOD2(i,j,2)=6 elseif (front(i,j-1)==3) then if (ifail.gt.0.and.METHOD(i,j,2).eq.'V-chang') then CALL U_COIN(+1) METHOD2(i,j,1)=10 write(num_err,*) 'modif',i,j,METHOD(i,j,1) else CALL V_PRESSIO(+1) METHOD(i,j,2)='V_P+1' METHOD2(i,j,2)=2 endif elseif (front(i,j-1)==2) then CALL V_COIN(+1) METHOD(i,j,2)='V_C+1' METHOD2(i,j,2)=4 else METHOD(i,j,2)='V_ORFAN' METHOD2(i,j,2)=995 SVIJ(i,j)=1 OPPOSY(i,j)=0!888666 ! print*, 'point flgzmY sans ice',i,j endif ENDIF !!test sur ice(i,j), fin du test ENDIF !!Test sur flgzmy ENDIF ENDDO ENDDO write(num_err,*) 'remplimat he' ! 4_Boucle de remplissage de la matrice. ! ------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! on determine la position de chaque vitesse Uij et Vij !!! dec 2003 pour l'instant on considere tous les points ! DO J=1,NY ! DO I=1,NX ! LIGU(i,j)=CTEURLIG ! LIGV(i,j)=CTEURLIG+1 ! cteurlig=cteurlig+2 ! ENDDO ! ENDDO ! nblig=cteurlig-1 !stop iinf=51 isup=53 jinf=53 jsup=55 CTEURLIG = 0 aaa=0. ! print*,'remplissage de mmat' ! print*,('shape=',shape(mmat)) ! print*,('size =',size(mmat)) ! print*,'remplissage de bdr' ! print*,('shape=',shape(bdr)) ! print*,('size =',size(bdr)) DO I=1,NX DO J=1,NY test_oku : IF (OKUMAT(I,J)) THEN CTEURLIG = CTEURLIG + 1 ! print*,i,j,'okumat=',OKUMAT(I,J),'CTEUR=',CTEURLIG,'LIGU',LIGU(I,J) test_shelfx : IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY) & .or.(.not.FLGZMX(I,J))) THEN mmat(kdc2,ligu(i,j)) = 1.0 ! BDR(ligu(i,j),1) =UXBAR(I,J) BDR(ligu(i,j),1) =OPPOSX(I,J) if (OPPOSX(I,J).ne.UXBAR(I,J)) then write(num_err,*) 'pb OPPOSX pose',i,j write(num_err,*) METHOD(i,j,1),OPPOSX(I,J),UXBAR(I,J) endif ELSE ! Classement de la matrice pour Uij ! --------------------------------- SCAL = TUIJ(i,j) mmat(kdc2,ligu(i,j)) = 1.0 iimat=kdc2+ligu(i,j) !print*,'uij mmat(kdc2,ligu)=mmat(',kdc2,ligu(i,j),')' BDR(ligu(i,j),1) = OPPOSX(i,j)/SCAL opposx(i,j) = OPPOSX(i,j)/SCAL if (LIGU(I-1,J-1).gt.0) then SNR = LIGU(I-1,J-1) ! les Ui-1 !print*, 'mmat(',iimat-snr,snr,') TUMIMJ' !print*,'iimat=',iimat,'snr=',snr mmat(iimat-snr,snr) = TUMIMJ(i,j)/SCAL endif if (LIGU(I-1,J).gt.0) then SNR = LIGU(I-1,J) mmat(iimat-snr,snr) = TUMIJ(i,j)/SCAL endif if (LIGU(I-1,J+1).gt.0) then SNR = LIGU(I-1,J+1) mmat(iimat-snr,snr) = TUMIPJ(i,j)/SCAL endif if (LIGU(I,J-1).gt.0) then SNR = LIGU(I,J-1) ! les Ui mmat(iimat-snr,snr) = TUIMJ(i,j)/SCAL endif if (LIGU(I,J+1).gt.0) then SNR = LIGU(I,J+1) mmat(iimat-snr,snr) = TUIPJ(i,j)/SCAL endif if (LIGU(I+1,J-1).gt.0) then SNR = LIGU(I+1,J-1) ! les Ui+1 mmat(iimat-snr,snr) = TUPIMJ(i,j)/SCAL endif if (LIGU(I+1,J).gt.0) then SNR = LIGU(I+1,J) mmat(iimat-snr,snr) = TUPIJ(i,j)/SCAL endif if (LIGU(I+1,J+1).gt.0) then SNR = LIGU(I+1,J+1) mmat(iimat-snr,snr) = TUPIPJ(i,j)/SCAL endif if (LIGV(I-1,J).gt.0) then SNR = LIGV(I-1,J) ! les Vi-1 mmat(iimat-snr,snr) = TVMIJ(i,j)/SCAL endif if (LIGV(I-1,J+1).gt.0) then SNR = LIGV(I-1,J+1) mmat(iimat-snr,snr) = TVMIPJ(i,j)/SCAL endif if (LIGV(I,J).gt.0) then SNR = LIGV(I,J) ! les Vi mmat(iimat-snr,snr) = TVIJ(i,j)/SCAL endif if (LIGV(I,J+1).gt.0) then SNR = LIGV(I,J+1) !print*, 'mmat(',iimat-snr,snr,') TVIPJ' !print*,'iimat=',iimat,'snr=',snr mmat(iimat-snr,snr) = TVIPJ(i,j)/SCAL endif ENDIF test_shelfx ENDIF test_oku test_okv : IF (OKVMAT(I,J)) THEN CTEURLIG = CTEURLIG + 1 ! print*,i,j,'okvmat=',OKUMAT(I,J),'CTEUR=',CTEURLIG,'LIGV',LIGV(I,J) test_shelfy : IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY) & .or.(.not.FLGZMY(I,J))) THEN mmat(kdc2,ligv(i,j)) = 1.0 ! BDR(ligv(i,j),1) =UYBAR(I,J) BDR(ligv(i,j),1) =OPPOSY(I,J) ELSE !test_shelfy ! Classement de la matrice pour Vij ! --------------------------------- SCAL = SVIJ(i,j) iimat=kdc2+ligv(i,j) !print*,'vij mmat(kdc2,ligv)=mmat(',kdc2,ligv(i,j),')' mmat(kdc2,ligv(i,j)) = 1.0 BDR(ligv(i,j),1) = OPPOSY(i,j)/SCAL if ((i.ge.67).and.(i.le.71)) then write(6,*)'i,j,opposy,scal',i,j,OPPOSY(i,j),scal,OPPOSY(i,j)/scal,method(i,j,2) end if opposy(i,j)=opposy(i,j)/scal if (LIGV(I-1,J-1).gt.0) then SNR = LIGV(I-1,J-1) ! les Vi-1 mmat(iimat-snr,snr) = SVMIMJ(i,j)/SCAL endif if (LIGV(I-1,J).gt.0) then SNR = LIGV(I-1,J) mmat(iimat-snr,snr) = SVMIJ(i,j)/SCAL endif if (LIGV(I-1,J+1).gt.0) then SNR = LIGV(I-1,J+1) !print*, 'mmat(',iimat-snr,snr,') SVMIPJ' !print*,'iimat=',iimat,'snr=',snr mmat(iimat-snr,snr) = SVMIPJ(i,j)/SCAL endif if (LIGV(I,J-1).gt.0) then SNR = LIGV(I,J-1) ! les Vi mmat(iimat-snr,snr) = SVIMJ(i,j)/SCAL endif if (LIGV(I,J+1).gt.0) then SNR = LIGV(I,J+1) mmat(iimat-snr,snr) = SVIPJ(i,j)/SCAL endif if (LIGV(I+1,J-1).gt.0) then SNR = LIGV(I+1,J-1) ! les Vi+1 mmat(iimat-snr,snr) = SVPIMJ(i,j)/SCAL endif if (LIGV(I+1,J).gt.0) then SNR = LIGV(I+1,J) mmat(iimat-snr,snr) = SVPIJ(i,j)/SCAL endif if (LIGV(I+1,J+1).gt.0) then SNR = LIGV(I+1,J+1) mmat(iimat-snr,snr) = SVPIPJ(i,j)/SCAL endif if (LIGU(I,J-1).gt.0) then SNR = LIGU(I,J-1) ! les Ui mmat(iimat-snr,snr) = SUIMJ(i,j)/SCAL endif if (LIGU(I,J).gt.0) then SNR = LIGU(I,J) mmat(iimat-snr,snr) = SUIJ(i,j)/SCAL endif if (LIGU(I+1,J-1).gt.0) then SNR = LIGU(I+1,J-1) !les Ui+1 mmat(iimat-snr,snr) = SUPIMJ(i,j)/SCAL endif if (LIGU(I+1,J).gt.0) then SNR = LIGU(I+1,J) mmat(iimat-snr,snr) = SUPIJ(i,j)/SCAL endif ENDIF test_shelfy ENDIF test_okv !print*,('classement des V fini') ENDDO ENDDO !print*,'CTEURLIG',CTEURLIG !print*,'NBLIG =',NBLIG NBLIG = CTEURLIG - 1 !print*,'fin de boucle de remplissage' ! stop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 5_Resolution du systeme lineaire ! ------------------------------ ifail=0 conjgr=.false. if (CONJGR) then !!!!!!!!!!!!debut test sur CONJGR DO K=1,NBLIG BDRO(K) = BDR(K,1) END DO ! nouveau conjgrad (utilisant blas) ! ------------------------- ! ----------------a enlever apres test----------------------- CONJGR = .true. DO I=1,NX DO J=1,NY ! DO I=ID1,ID2 ! DO J=JD1,JD2 IF (OKUMAT(I,J)) THEN SAV(LIGU(I,J))=UXPREC(I,J) END IF IF (OKVMAT(I,J)) THEN SAV(LIGV(I,J))=UYPREC(I,J) END IF END DO END DO ! -----------------fin de la partie a enlever--------------- ! write(6,*)'avant conjgrad time=',time,' ldab=',ldab ! call conjgrad(nblig,kl,ku,mmat,ldab,BDR,SAV,rsq,iter,ifail) ! write(6,*)'time',time,' numdom',numdom,' iter',iter ! CALL CONJGRAD(BDRO,NBLIG,SAV,RSQ,ifail,ITER) if (iter.gt.1) write(77,*) 'time=',time,' iter=',iter IF (ifail.EQ.1) THEN write(*,*) 'gradients conjugues insuffisants' CONJGR=.FALSE. TIMECG=TIME goto 911 ENDIF ELSE !!!!!!!!!!!!continue test sur CONJGR ! write(6,*) ! write(6,*) 'avant entree sgbsv, domaine',numdom ! write(6,*) 'nblig',nblig,' kl',kl,' ku',ku,' nrhs',nrhs ! write(6,*) 'ldab',ldab,' ldb',ldb knonul=0 emmax=0. !----pour mettre mmat=identite----! ! do i=1,ldbmax ! do j=1,nptmax ! if (i.eq.ku+kl+1) then ! mmat(i,j)=1. ! else ! mmat(i,j)=0. ! endif ! end do ! end do !---------------------------------! do i=1,nblig bdro(i)=bdr(i,1) end do ! pour verifier qu'on a pas oublie de point do i=1,nblig if (mmat(kl+ku+1,i).eq.0.) then print*,(mmat(kl+ku+1,i)) print*,'kl+ku+1',kl+ku+1 write(6,*) '******* ATTENTION point oublie numero' & ,i,'****************************' endif end do !----chercher le maximum de la matrice-----------------! ! emmax=0. ! do i=1,nblig ! do j=1,ldab ! emmax=max(abs(mmat(j,i)),emmax)!-! ! end do ! end do ! write(6,*) emmax !------------------------------------------------------! !open(UNIT=994,file='indexes') ! do i=1,nx ! do j=1,ny ! write(994,*) 'i ,j ,LIGU,LIGV' ! write(994,*) i,j,LIGU(i,j),LIGV(i,j) ! enddo ! enddo ! close (994) !open(UNIT=995,file='bdr') ! do cteurlig=1,nlig ! write(995,*) bdr(cteurlig),cteurlig ! enddo ! close (995) !open(UNIT=996,file='mmat_columns_before') !write(996,*)'bdr(u(10,10))=',bdr(LIGU(10,10)) !write(996,*)'bdr(v(10,10))=',bdr(LIGV(10,10)) !do i=1,2*ku+ku+1 ! write(996,*) mmat(i,LIGU(10,10)),mmat(i,LIGV(10,10)) !enddo !close (996) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! write(6,*)'avant sgbsv time=',time,' ldab=',ldab ! write(6,*)'nblig',nblig,'ku=',ku,'kl=',kl,'nrhs',nrhs ! write(6,*)'ipiv ?','ldb ',ldb,'ifail',ifail !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*, 'min value, location',minval(mmat),minloc(mmat) ! print*, 'max value, location',maxval(mmat),maxloc(mmat) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! sorties des Tuij debug_3D(:,:,26)=TUIJ(:,:) ! write(6,*)'avant sgbsv', nblig,kl,ku,nrhs,ldbmax ! quelques sorties netcdf where (okumat(:,:)) debug_3D(:,:,39)=1 elsewhere debug_3D(:,:,39)=0 end where where (okvmat(:,:)) debug_3D(:,:,40)=1 elsewhere debug_3D(:,:,40)=0 end where debug_3d(:,:,35) = opposx(:,:) debug_3d(:,:,36) = opposy(:,:) CALL SGBSV(nblig,kl,ku,nrhs,mmat,ldbmax,ipiv,bdr,ldb,ifail) ! write(6,*)'apres sgbsv' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! open(UNIT=998,file='sgbsv_vit_u.dat') ! open(UNIT=999,file='sgbsv_vit_v.dat') ! do i=1,nblig ! if (bdro(i).ne.bdr(i,1)) then !write(6,*) i, bdro(i)-bdr(i,1) ! endif ! end do ! call DEBUGAGE ! write(num_ic_by,*) mmat ! write(6,*) 'fin de l ecriture de mmat' IF (ifail.NE.0) THEN nb_err=nb_err+1 WRITE(*,*) ' ======================================' WRITE(*,*) ' ERREUR DANS L''ELIMINATION GAUSSIENNE ' WRITE(*,*) ' ERREUR N ',ifail WRITE(*,*) ' ======================================' if (posligu(ifail,1).gt.0) then iii=posligu(ifail,1) jjj=posligu(ifail,2) print*,'Ux',iii,jjj,method(iii,jjj,1) METHOD(iii,jjj,1)='U-chang' elseif (posligv(ifail,1).gt.0) then iii=posligv(ifail,1) jjj=posligv(ifail,2) print*,'Vy',iii,jjj,method(iii,jjj,2) METHOD(iii,jjj,2)='V-chang' endif if (geoplace.eq.'anteis1'.and.H(71,71).lt.10.) then print*,time,ifail print*,'stop apres glace disparut a Pole sud h=',h(71,71) stop endif call DEBUGAGE !!! if (nb_err.lt.2) then !!! goto 942 !!! else print*,'=apres',nb_err,'passage infructueux : on stop' ! pause return ! STOP !!! endif END IF ! call DEBUGAGE DO I=1,NX DO J=1,NY ! if ((I.ge.ID1).and.(I.le.ID2).and. & ! (J.ge.JD1).and.(J.le.JD2).and.(OKUMAT(I,J))) then if (OKUMAT(I,J)) then UXNEW(I,J)=BDR(LIGU(I,J),1) SAV(LIGU(I,J))=UXNEW(I,J) else UXNEW(I,J)=UXPREC(I,J) endif ! if ((I.ge.ID1).and.(I.le.ID2).and. & ! (J.ge.JD1).and.(J.le.JD2).and.(OKVMAT(I,J))) then if (OKVMAT(I,J)) then UYNEW(I,J)=BDR(LIGV(I,J),1) SAV(LIGV(I,J))=UYNEW(I,J) else UYNEW(I,J)=UYPREC(I,J) endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !open(UNIT=998,file='sgbsv_vit_u.dat') !open(UNIT=999,file='sgbsv_vit_v.dat') ! write (998,*) i,j,UXNEW(I,J) ! write (999,*) i,j,UYNEW(I,J) if (ICE(i,j)==0) then if (.not.flgzmx(i,j)) then UXNEW(I,J)=3.3333333 endif if (.not.flgzmy(i,j)) then UYNEW(I,J)=3.3333333 endif endif !close (998) !close (999) !!!!!!!!!!!!!test vincent : limiter vitesses apres sgbsv UXNEW(I,J)=max(-3900.,UXNEW(I,J)) UXNEW(I,J)=min( 3900.,UXNEW(I,J)) UYNEW(I,J)=max(-3900.,UYNEW(I,J)) UYNEW(I,J)=min( 3900.,UYNEW(I,J)) END DO END DO where ((Hmx(:,:).le.1.).and.(flgzmx(:,:))) uxnew(:,:)=0. end where where ((Hmy(:,:).le.1.).and.(flgzmy(:,:))) uynew(:,:)=0. end where debug_3D(:,:,27)=uxnew(:,:) ! print*,'???',UXNEW(57,40),UXNEW(58,40),UXNEW(59,40) ! pause ! close (998) !!!!!!! close (999) ENDIF !!!!!!!!!!!!fin test sur CONJGR ! impression des vitesses ! open (57,file='Dom'//domname//domtot2(1:jn)) ! write(57,*) 'time=',time ! write(57,*)'i j vitesse old-vitesse delta-vitesse ' ! do i=2,nx-1 ! do j=2,ny-1 ! UNW = SQRT(((UXnew(I,J)+UXnew(I+1,J))/2)**2 ! & + ((UYnew(I,J)+UYnew(I,J+1))/2.)**2) ! UPR = SQRT(((UXPREC(I,J)+UXPREC(I+1,J))/2)**2 ! & + ((UYPREC(I,J)+UYPREC(I,J+1))/2.)**2) ! write(57,958) i,j,unw, upr,unw-upr ! end do ! end do ! close(57) !958 format(i3,1x,i3,2x,3(e10.4,2x)) ! print*,'fin de remplimat' 911 continue close(num_err) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 6_routines pour referencer les coefficients utiles pour remplir la matrice ! ------------------------------------------------------------------------ CONTAINS ! Les 10 procedures suivantes sont internes a remplidom ! U_* V_* determinent les procedures de calcul du shelf ! Cas general U_GENERAL V_GENERAL ! face a un front action (normale au front) de la pression U_PERSSION V_PERSSION ! en bordure du font condition de non cisaillzment U_CISAILL V_CISAILL ! langue non confinee ISOLXBACKW ISOLYBACKW calcul des vitesses en backward ! languee non confinee ISOLXFOREW ISOLYFOREW calcul des ! vitesses en foreward ! LARGEUR_DE_BANDE determine les ku et kl minimums ! DEBUG ecrit une serie de fichiers TABLE_qch qui donne les coeff de ! mmat, la methode utilisees SUBROUTINE U_GENERAL !Rempli la matice pour le calcul de Uij IMPLICIT NONE ! Calcul des coefficients de l'equation en U: ! ------------------------------------------- beta=FROTMX(I,J)*dx*dx ! Terme en u(i,j): ! _______________ TUIJ(I,J) = -4.*PVI(I,J) - 4.*PVI(I-1,J) & - PVM(i,j+1)-PVM(I,J)-beta ! Terme en u(i-1,j): ! _________________ TUMIJ(I,J) = 4.*PVI(I-1,J) ! Terme en u(i+1,j): ! _________________ TUPIJ(I,J) = 4.*PVI(I,J) ! Terme en u(i,j-1): ! _________________ TUIMJ(I,J) = PVM(I,J) ! Terme en u(i,j+1): ! _________________ TUIPJ(I,J) = PVM(I,J+1) ! Terme en v(i,j): ! _______________ TVIJ(I,J) = -2.*PVI(I,J)-PVM(I,J) ! Terme en v(i-1,j): ! _________________ TVMIJ(I,J) = 2.*PVI(I-1,J)+PVM(I,J) ! Terme en v(i-1,j+1): ! ___________________ TVMIPJ(I,J) = -2.*PVI(I-1,J)-PVM(I,J+1) ! Terme en v(i,j+1): ! _________________ TVIPJ(I,J) = 2.*PVI(I,J)+PVM(I,J+1) moteur= RO*G*HMX(I,J)* & (S(I,J)-S(I-1,J))/dx moteur=min(moteur,moteurmax) moteur=max(moteur,-moteurmax) OPPOSX(I,J) = moteur*dx*dx !print*,'u_gen OPPOS',OPPOSX(I,J),i,j !print*,'motmax=',moteurmax,'motmax*dx2=',moteurmax*dx*dx END SUBROUTINE U_GENERAL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_GENERAL !Rempli la matice pour le calcul de Vij IMPLICIT NONE ! print*, 'V_GENERAL',i,j ! Calcul des coefficients de l'equation en V: ! ------------------------------------------- beta=FROTMY(I,J)*dx*dx ! Terme en v(i,j): ! ________________ SVIJ(I,J) = -4.*PVI(I,J)-4.*PVI(I,J-1) & - PVM(I+1,J)-PVM(I,J)-beta ! Terme en v(i,j-1): ! __________________ SVIMJ(I,J) = 4.*PVI(I,J-1) ! Terme en v(i,j+1): ! __________________ SVIPJ(I,J) = 4.*PVI(I,J) ! Terme en v(i-1,j): ! __________________ SVMIJ(I,J) = PVM(I,J) ! Terme en v(i+1,j): ! __________________ SVPIJ(I,J) = PVM(I+1,J) ! Terme en u(i,j): ! ________________ SUIJ(I,J) = -2.*PVI(I,J)-PVM(I,J) ! Terme en u(i,j-1): ! __________________ SUIMJ(I,J) = 2.*PVI(I,J-1)+PVM(I,J) ! Terme en u(i+1,j-1): ! ___________________ SUPIMJ(I,J) = -2.*PVI(I,J-1)-PVM(I+1,J) ! Terme en u(i+1,j): ! ___________________ SUPIJ(I,J) = 2.*PVI(I,J)+PVM(I+1,J) moteur= RO*G*HMY(I,J)* & (S(I,J)-S(I,J-1))/dx moteur=min(moteur,moteurmax) moteur=max(moteur,-moteurmax) OPPOSY(i,j) = moteur*dx*dx END SUBROUTINE V_GENERAL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE U_PRESSIO(DIR) ! sert quand front=3 et frontfacex=+-1 !Rempli la matice pour le calcul de Uij IMPLICIT NONE INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1 ! Equivaut frontfacey if (DIR.EQ.1) then !front a droite (case vide a droite de Uij) beta=FROTMX(I,J)*dx beta=min(beta,0.3*4*pvi(I-1,J)) TUIJ(I,J) = 4*PVI(I-1,J) - beta TUMIJ(I,J) =-4*PVI(I-1,J) TVMIPJ(I,J) = 2*PVI(I-1,J) TVMIJ(I,J) =-2*PVI(I-1,J) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I-1,J)**2)*dx/2. OPPOSX(I,J) = ROG*(H(I-1,J)**2)*dx/2. - ROWG*((min(slv(i-1,j)-B(i-1,j),0.))**2)*dx/2. else !front a gauche (case vide a gauche de Uij) beta=FROTMX(I,J)*dx beta=min(beta,0.3*4*pvi(I,J)) TUIJ(I,J) =-4*PVI(I,J) + beta TUPIJ(I,J) = 4*PVI(I,J) TVIPJ(I,J) = 2*PVI(I,J) TVIJ(I,J) =-2*PVI(I,J) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSX(I,J) = ROG*(H(I,J)**2)*dx/2. - ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. endif END SUBROUTINE U_PRESSIO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_PRESSIO(DIR) ! sert quand front=3 et frontfacey=+-1 !Rempli la matice pour le calcul de Vij IMPLICIT NONE INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1 ! Equivaut frontfacey if (DIR.EQ.1) then!front en haut (case vide EN HAUT de Vij) beta=FROTMY(I,J)*dx beta=min(beta,0.3*4*pvi(I,J-1)) SVIJ(I,J) = 4*pvi(I,J-1) - beta SVIMJ(I,J) =-4*pvi(I,J-1) SUPIMJ(I,J) = 2*pvi(I,J-1) SUIMJ(I,J) =-2*pvi(I,J-1) OPPOSY(I,J) = ROG*(H(I,J-1)**2)*dx/2. - ROWG*((min(slv(i,j-1)-B(i,j-1),0.))**2)*dx/2. else!front en bas (case vide en bas de Vij) ! print*, 'V_PRESSIO(-1)',i,j beta=FROTMY(I,J)*dx beta=min(beta,0.3*4*pvi(I,J)) SVIJ(I,J) =-4*pvi(I,J) + beta SVIPJ(I,J) = 4*pvi(I,J) SUPIJ(I,J) = 2*pvi(I,J) SUIJ(I,J) =-2*pvi(I,J) !!!! OPPOSY(I,J) = ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSY(I,J) = ROG*(H(I,J)**2)*dx/2. - ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. end if END SUBROUTINE V_PRESSIO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE U_COIN(CAS) ! sert quand front=2 et frontfacex=+-1,frontfacey=+-1 !Rempli la matice pour le calcul de Uij IMPLICIT NONE integer CAS !valeur de frontfacex :si le front est a-droite ou !a-gauche du point en glace if (cas==1) then beta=FROTMX(i,j)*dx beta=min(beta,0.3*PVI(i-1,j)) ! TUIJ(I,J) = 6*PVI(i-1,j) - beta ! TUMIJ(I,J) =-6*PVI(i-1,j) TUIJ(I,J) = PVI(i-1,j) - beta TUMIJ(I,J) =-PVI(i-1,j) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I-1,J)**2)*dx/2. OPPOSX(i,j)= ROG*(H(I-1,J)**2)*dx/2. - ROWG*((min(slv(i-1,j)-B(i-1,j),0.))**2)*dx/2. else beta=FROTMX(i,j)*dx beta=min(beta,0.3*PVI(i,j)) ! TUIJ(I,J) =-6*PVI(i,j) + beta ! TUPIJ(I,J) = 6*PVI(i,j) TUIJ(I,J) =-PVI(i,j) + beta TUPIJ(I,J) = PVI(i,j) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSX(i,j)= ROG*(H(I,J)**2)*dx/2.- ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. endif END SUBROUTINE U_COIN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_COIN(CAS) ! sert quand front=2 et frontfacex=+-1,frontfacey=+-1 !Rempli la matice pour le calcul de Vij IMPLICIT NONE integer CAS !valeur de frontfacey : si le front est au-dessus ou !au-dessous du point en glace if (cas==1) then beta=FROTMY(I,J)*dx beta=min(beta,0.3*PVI(i,j-1)) ! SVIJ(I,J) = 6*PVI(i,j-1) - beta ! SVIMJ(I,J) =-6*PVI(i,j-1) SVIJ(I,J) = PVI(i,j-1) - beta SVIMJ(I,J) =-PVI(i,j-1) !!!! OPPOSY(I,J) = ROG*(1.-RO/ROW)*(H(I,J-1)**2)*dx/2. OPPOSY(i,j)= ROG*(H(I,J-1)**2)*dx/2. - ROWG*((min(slv(i,j-1)-B(i,j-1),0.))**2)*dx/2. else beta=FROTMY(I,J)*dx beta=min(beta,0.3*PVI(i,j)) ! SVIJ(i,j) =-6*PVI(i,j) + beta ! SVIPJ(i,j) = 6*PVI(i,j) SVIJ(i,j) =-PVI(i,j) + beta SVIPJ(i,j) = PVI(i,j) !!!! OPPOSY(i,j)= ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSY(i,j)= ROG*(H(I,J)**2)*dx/2. - ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. endif END SUBROUTINE V_COIN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE U_CISAILL(DIR) !Rempli la matice pour le calcul de Uij IMPLICIT NONE INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1 ! Equivaut frontfacey IF (DIR==1) THEN if (.not.(flotmy(i,j).and.flotmy(i-1,j))) then ! if (.not.(flgzmy(i,j).and.flgzmy(i-1,j))) then !On ne calcule pas le cisaillement sur le bord des streams if (ice(i+1,j)==0.and.ice(i-2,j)==0) then TUIJ(I,J) = 1 TUIMJ(I,J) =-1 else TUIJ(I,J) = 1 TUPIJ(I,J) =-0.5 TUMIJ(I,J) =+0.5 endif else !print*, 'U_CISAILL(+1)',i,j !Le front est en haut, on calcule EXY=0 vers le bas (en arriere) !du/dy TUIJ(I,J) = 1 TUIMJ(I,J) =-1 !dv/dx TVIPJ(I,J) = 1 TVMIPJ(I,J)=-1 endif ELSE if (.not.(flotmy(i,j+1).and.flotmy(i-1,j+1))) then ! if (.not.(flgzmy(i,j+1).and.flgzmy(i-1,j+1))) then !On ne calcule pas le cisaillement sur le bord des streams if (ice(i+1,j)==0.and.ice(i-2,j)==0) then TUIJ(I,J) = 1 TUIPJ(I,J) =-1 else TUIJ(I,J) = 1 TUPIJ(I,J) =-0.5 TUMIJ(I,J) =+0.5 endif else !Le front est en bas, on calcule EXY=0 vers le haut(en arriere) !du/dy TUIJ(I,J) =-1 TUIPJ(I,J) = 1 !dv/dx TVIJ(I,J) = 1 TVMIJ(I,J) =-1 endif ENDIF OPPOSX(I,J)= 0 END SUBROUTINE U_CISAILL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_CISAILL(DIR) !Rempli la matice pour le calcul de Vij IMPLICIT NONE INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1 ! Equivaut frontfacex IF (DIR==1) THEN ! if (i.eq.54) print*,i,j,flotmx(i,j),flotmx(i,j-1) if ((.not.flotmy(i,j)).or..not.(flotmx(i,j).and.flotmx(i,j-1))) then ! if (.not.(flgzmx(i,j).and.flgzmx(i,j-1))) then !On ne calcule pas le cisaillement sur le bord des streams if (ice(i,j+1)==0.and.ice(i,j-2)==0) then SVIJ(I,J) = 1 SVMIJ(I,J) =-1 else SVIJ(I,J) = 1 SVIPJ(I,J) =-0.5 SVIMJ(I,J) =+0.5 endif else !Le front est a droite, on calcule EXY=0 vers la gauche (en arriere) !dv/dx SVIJ(I,J) = 1 SVMIJ(I,J) =-1 !du/dy SUPIJ(I,J) = 1 SUPIMJ(I,J) =-1 endif ELSE if ((.not.flotmy(i,j)).or..not.(flotmx(i+1,j).and.flotmx(i+1,j-1))) then ! if (.not.(flgzmx(i+1,j).and.flgzmx(i+1,j-1))) then !On ne calcule pas le cisaillement sur le bord des streams if (ice(i,j+1)==0.and.ice(i,j-2)==0) then SVIJ(I,J) = 1 SVPIJ(I,J) =-1 else SVIJ(I,J) = 1 SVIPJ(I,J) =-0.5 SVIMJ(I,J) =+0.5 endif else !Le front est a gauche, on calcule EXY=0 vers la droite (en arriere) !dv/dx SVIJ(I,J) =-1 SVPIJ(I,J) = 1 !du/dy SUIJ(I,J) = 1 SUIMJ(I,J) =-1 endif ENDIF OPPOSY(I,J)= 0 END SUBROUTINE V_CISAILL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_TONGUE !la zone a 2 voisins alignes selon l'axe Y : on traite Vij beta=FROTMY(I,J)*dx*dx SVIJ(I,J) =-6*PVI(i,j-1) - 6*PVI(i,j) - beta SVIMJ(I,J) = 6*PVI(i,j-1) SVIPJ(I,J) = 6*PVI(i,j) moteur= ROG*HMX(I,J)* & (S(I,J)-S(I,J-1))/dx moteur=min(moteur,moteurmax) moteur=max(moteur,-moteurmax) OPPOSY(I,J)=moteur*dx*dx END SUBROUTINE V_TONGUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_TONGUEFRONT(cas) !la langue de glace est verticale (selon les Y): calcul Vij sur le front integer CAS !valeur de frontfacey : si le front est au-dessus ou !au-dessous du point en glace if (cas==1) then!le point en glace est dessous: ICE(i,j-1)=1 beta=FROTMY(I,J)*dx SVIJ(I,J) = 6*PVI(i,j-1) - beta SVIMJ(I,J) =-6*PVI(i,j-1) ! SVIJ(I,J) = PVI(i,j-1) - beta ! SVIMJ(I,J) =-PVI(i,j-1) !!!! OPPOSY(I,J) = ROG*(1.-RO/ROW)*(H(I,J-1)**2)*dx/2. OPPOSY(i,j)= ROG*(H(I,J-1)**2)*dx/2. - ROWG*((min(slv(i,j-1)-B(i,j-1),0.))**2)*dx/2. else !cas=-1 le point en glace est dessus : ICE(i,j)=1 beta=FROTMY(I,J)*dx SVIJ(i,j) =-6*PVI(i,j) - beta SVIPJ(i,j) = 6*PVI(i,j) ! SVIJ(i,j) =-PVI(i,j) - beta ! SVIPJ(i,j) = PVI(i,j) !!!! OPPOSY(i,j)= ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSY(i,j)= ROG*(H(I,J)**2)*dx/2.- ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. endif END SUBROUTINE V_TONGUEFRONT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE U_TONGUE !la zone a 2 voisins alignes selon l'axe X : on traite Uij beta=FROTMX(i,j)*dx*dx TUIJ(I,J) =-6*PVI(i-1,j) - 6*PVI(i,j) - beta TUMIJ(I,J) = 6*PVI(i-1,j) TUPIJ(I,J) = 6*PVI(i,j) moteur= RO*G*HMX(I,J)* & (S(I,J)-S(I-1,J))/dx moteur=min(moteur,moteurmax) moteur=max(moteur,-moteurmax) OPPOSX(I,J)=moteur*dx*dx END SUBROUTINE U_TONGUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE U_TONGUEFRONT(cas) !langue de glace horizontale (selon les X): calcul Uij sur le front integer CAS !valeur de frontfacex :si le front est a-droite ou !a-gauche du point en glace if (cas==1) then!le point en glace est a gauche: ICE(i-1,j)=1 beta=FROTMX(i,j)*dx TUIJ(I,J) = 6*PVI(i-1,j) - beta TUMIJ(I,J) =-6*PVI(i-1,j) ! TUIJ(I,J) = PVI(i-1,j) - beta ! TUMIJ(I,J) =-PVI(i-1,j) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I-1,J)**2)*dx/2. OPPOSX(I,J) = ROG*(H(I-1,J)**2)*dx/2. - ROWG*((min(slv(i-1,j)-B(i-1,j),0.))**2)*dx/2. else !cas=-1 le point en glace est a droite: ICE(i,j)=1 beta=FROTMX(i,j)*dx TUIJ(I,J) =-6*PVI(i,j) - beta TUPIJ(I,J) = 6*PVI(i,j) ! TUIJ(I,J) =-PVI(i,j) - beta ! TUPIJ(I,J) = PVI(i,j) !!!! OPPOSX(I,J) = ROG*(1.-RO/ROW)*(H(I,J)**2)*dx/2. OPPOSX(I,J) = ROG*(H(I,J)**2)*dx/2. - ROWG*((min(slv(i,j)-B(i,j),0.))**2)*dx/2. endif END SUBROUTINE U_TONGUEFRONT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE V_TONGUE_UP(front) integer front ! valeur de frontfacey(i,j-1) OPPOSY(I,J)= 0 if (front==0) then!voisin en glace (en bas) a gauche et a doite SVIJ(I,J) = 1 SVMIJ(I,J) =-.25 SVMIMJ(I,J)=-.25 SVPIJ(I,J) =-.25 SVPIMJ(I,J)=-.25 SUPIMJ(I,J)=-.5 SUIMJ(I,J) = .5 elseif (front==1) then!voisin en glace (en bas) a gauche SVIJ(I,J) = 1 SVMIJ(I,J) =-.5 SVMIMJ(I,J)=-.5 SUPIMJ(I,J)=-.5 SUIMJ(I,J) = .5 else !!voisin en glace (en bas) a doite SVIJ(I,J) = 1 SVPIJ(I,J) =-.5 SVPIMJ(I,J)=-.5 SUPIMJ(I,J)=-.5 SUIMJ(I,J) = .5 endif END SUBROUTINE V_TONGUE_UP SUBROUTINE V_TONGUE_DO(front)!DO=down integer front ! valeur de frontfacey(i,j) OPPOSY(I,J)= 0 if (front==0) then!voisin en glace a gauche et a doite SVIJ(I,J) = 1 SVMIPJ(I,J)=-.25 SVMIJ(I,J) =-.25 SVPIPJ(I,J)=-.25 SVPIJ(I,J) =-.25 SUPIJ(I,J) = .5 SUIJ(I,J) =-.5 elseif (front==1) then!voisin en glace a gauche SVIJ(I,J) = 1 SVMIPJ(I,J)=-.5 SVMIJ(I,J) =-.5 SUPIJ(I,J) = .5 SUIJ(I,J) =-.5 else!voisin en glace a doite SVIJ(I,J) = 1 SVPIPJ(I,J)=-.5 SVPIJ(I,J) =-.5 SUPIJ(I,J) = .5 SUIJ(I,J) =-.5 endif END SUBROUTINE V_TONGUE_DO SUBROUTINE U_TONGUE_LE(front)!LE=left integer front ! valeur de frontfacex(i,j) OPPOSX(I,J)= 0 if (front==0) then!voisin en glace en haut et en bas TUIJ(I,J) = 1 TUPIMJ(I,J)=-.25 TUIMJ(I,J) =-.25 TUPIPJ(i,j)=-.25 TUIPJ(i,j) =-.25 TVIPJ(I,J) = .5 TVIJ(I,J) =-.5 elseif(front==1) then!voisin en glace en bas TUIJ(I,J) = 1 TUPIMJ(I,J)=-.5 TUIMJ(I,J) =-.5 TVIPJ(I,J) = .5 TVIJ(I,J) =-.5 else!voisin en glace en haut TUIJ(i,j) =1 TUIPJ(i,j) =-0.5 TUPIPJ(i,j)=-0.5 TVIPJ(i,j) = 0.5 TVIJ(i,j) =-0.5 endif END SUBROUTINE U_TONGUE_LE SUBROUTINE U_TONGUE_RI(front) integer front ! valeur de frontfacex(i-1,j) OPPOSX(I,J)= 0 if (front==0) then!voisin en glace (a gauche) en haut et en bas TUIJ(I,J) = 1 TUIMJ(I,J) =-.25 TUMIMJ(I,J)=-.25 TUMIPJ(i,j)=-.25 TUIPJ(i,j) =-.25 TVMIPJ(I,J)=-.5 TVMIJ(I,J) = .5 elseif (front==1) then!voisin en glace (a gauche) en bas TUIJ(I,J) = 1 TUIMJ(I,J) =-.5 TUMIMJ(I,J)=-.5 TVMIPJ(I,J)=-.5 TVMIJ(I,J) = .5 else!!voisin en glace (a gauche) en haut TUIJ(i,j) =1 TUMIPJ(i,j)=-0.5 TUIPJ(i,j) =-0.5 TVMIPJ(i,j)=-0.5 TVMIJ(i,j) = 0.5 endif END SUBROUTINE U_TONGUE_RI !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! -------------------------------------------------------------- ! subroutine de calcul de la largeur de bande ! version sous-domaine ! -------------------------------------------------------------- subroutine larg_bande(ksur,ksous) implicit none INTEGER KSUR,KSOUS,KMIN,KMAX,IKSUR,IKSOUS,JKSUR,JKSOUS,ktest KSOUS=-1 KSUR=-1 ! noeuds U ! do I=1,NX ! do J=1,NY do I=2,NX-1 do J=2,NY-1 if (OKUMAT(I,J)) then kmin=ligu(i,j) kmax=ligu(i,j) ktest=ligu(i-1,j-1) !Ui-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i-1,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i-1,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i,j-1) !Ui if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i+1,j-1) !Ui-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i+1,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i+1,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i-1,j) !Vi-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i-1,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i,j) !Vi if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ! ------------------------------- fin de recherche min max de la ligne if (LIGU(I,J)-KMIN.gt.KSOUS) then KSOUS=LIGU(I,J)-KMIN IKSOUS=I JKSOUS=J endif if (KMAX-LIGU(I,J).gt.KSUR) then KSUR=KMAX-LIGU(I,J) IKSUR=I JKSUR=J endif endif end do end do ! ----------------- impression de largeur de bande ------------------- ! write(6,*)'largeur de bande pour les noeuds internes en U' ! write(6,*)'ksous=',ksous,' pour i=',iksous,' pour j=',jksous ! write(6,*)'ksur=',ksur,' pour i=',iksur,' pour j=',jksur ! -------------------------------------------------------------------- ! noeuds V ! KSOUS=-1 ! KSUR=-1 ! do I=1,NX ! do J=1,NY do I=2,NX-1 do J=2,NY-1 if (OKVMAT(I,J)) then kmin=ligv(i,j) kmax=ligv(i,j) ktest=ligv(i-1,j-1) !Vi-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i-1,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i-1,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i,j-1) !Vi if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i+1,j-1) !Vi-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i+1,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligv(i+1,j+1) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i,j-1) !Ui if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i+1,j-1) !Ui-1 if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ktest=ligu(i+1,j) if (ktest.gt.0) then kmin=min(ktest,kmin) kmax=max(ktest,kmax) endif ! ------------------------------- fin de recherche min max de la ligne if (LIGV(I,J)-KMIN.gt.KSOUS) then KSOUS=LIGV(I,J)-KMIN IKSOUS=I JKSOUS=J endif if (KMAX-LIGV(I,J).gt.KSUR) then KSUR=KMAX-LIGV(I,J) IKSUR=I JKSUR=J endif endif end do end do ! ----------------- impression de largeur de bande ------------------- ! write(6,*)'largeur de bande pour les noeuds internes en V' ! write(6,*)'ksous=',ksous,' pour i=',iksous,' pour j=',jksous ! write(6,*)'ksur=',ksur,' pour i=',iksur,' pour j=',jksur ! --------------------------------------------------------------- end subroutine larg_bande !---------------------------------------------------------------------- ! FIN DU CALCUL DE LARGEUR DE BANDE !---------------------------------------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE DEBUGAGE implicit none logical yesno ! character aaa integer iii,jjj character(len=80) :: filin real,dimension(-1:1,-1:1) :: Au, Bv real,dimension(-1:0,0:1) :: Av real,dimension(0:1,-1:0) :: Bu yesno=.false. if (yesno) then do i=2,nx-1 do j=2,ny-1 Au(-1,-1)=TUMIMJ(I,J) Au(-1,0)=TUMIJ(I,J) Au(-1,1)=TUMIPJ(I,J) Au(0,-1)=TUIMJ(I,J) Au(0,0)=TUIJ(I,J) Au(0,1)=TUIPJ(I,J) Au(1,-1)=TUPIMJ(I,J) Au(1,0)=TUPIJ(I,J) Au(1,1)=TUPIPJ(I,J) !!!!!!!! Av(-1,0)=TVMIJ(I,J) Av(-1,1)=TVMIPJ(I,J) Av(0,0)=TVIJ(I,J) Av(0,1)=TVIPJ(I,J) !!!!!!!! Bv(-1,-1)=SVMIMJ(I,J) Bv(-1,0)=SVMIJ(I,J) Bv(-1,1)=SVMIPJ(I,J) Bv(0,-1)=SVIMJ(I,J) Bv(0,0)=SVIJ(I,J) Bv(0,1)=SVIPJ(I,J) Bv(1,-1)=SVPIMJ(I,J) Bv(1,0)=SVPIJ(I,J) Bv(1,1)=SVPIPJ(I,J) !!!!!!!! Bu(0,0)=SUIJ(I,J) Bu(0,-1)=SUIMJ(I,J) Bu(1,0)=SUPIJ(I,J) Bu(1,-1)=SUPIMJ(I,J) do iii=-1,1 do jjj=-1,1 if ((Au(iii,jjj).ne.0).and. & (abs(OPPOSX(I+iii,J+jjj)/Au(0,0)).gt.100)) then write(6,*) 'erreur i,j=',i,j,'Au(',iii,jjj,')' write(6,*) 'OPPOSX=',OPPOSX(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2) ! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop write(6,*)'Au(0,0)=',Au(0,0),'OPPOSX/Au(0,0)',OPPOSX(I+iii,J+jjj)/Au(0,0) endif if ((Bv(iii,jjj).ne.0).and. & (abs(OPPOSY(I+iii,J+jjj)/Bv(0,0)).gt.100)) then write(6,*) 'erreuri,j=',i,j,'Bv(',iii,jjj,')' write(6,*) 'OPPOSY=',OPPOSY(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2) write(6,*)'Bv(0,0)=',Bv(0,0),'OPPOSY/Bv(0,0)',OPPOSY(I+iii,J+jjj)/Bv(0,0) ! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop endif enddo enddo do iii=-1,0 do jjj=0,1 if ((Av(iii,jjj).ne.0).and. & (abs(OPPOSY(I+iii,J+jjj)/Au(0,0)).gt.100)) then write(6,*) 'erreur i,j=',i,j,'Av(',iii,jjj,')' write(6,*) 'OPPOSY=',OPPOSY(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2) write(6,*) 'AV(0,0)=',Av(0,0),'OPPOSY/Av(0,0)',OPPOSY(I+iii,J+jjj)/Av(0,0) ! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop endif enddo enddo do jjj=-1,0 do iii=0,1 if ((Bu(iii,jjj).ne.0).and. & (abs(OPPOSX(I+iii,J+jjj)/Bv(0,0)).gt.100)) then write(6,*) 'erreur i,j=',i,j,'Bu(',iii,jjj,')' write(6,*) 'OPPOSX=',OPPOSX(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2) write(6,*) 'Bu(0,0)=',Bu(0,0),'OPPOSX/Bu(0,0)',OPPOSX(I+iii,J+jjj)/Bu(0,0) ! if(iii.ne.0.or.jjj.ne.0) read*, aaa!stop endif enddo enddo write(6,*)'===================================================' enddo enddo endif !!!!!!!!!!!impression de la methode!!!!!!!!!!!!!!!!! filin = 'METHOD_uv' open (UNIT=num_file3,file=filin) write(num_file3,*) geoplace,'time=',time write(num_file3,*) NX*NY,DX,NX,NY do j=1,ny do i=1,nx ! do i=35,60 ! do j=35,60 write(num_file3,*) i,j, METHOD(i,j,1),METHOD(i,j,2) enddo enddo close (num_file3) filin = 'TABLE_met_U' open (UNIT=num_file3,file=filin) write(num_file3,*) geoplace,'time=',time write(num_file3,*) NX*NY,DX,NX,NY do j=1,ny do i=1,nx write(num_file3,*) i,j, METHOD2(i,j,1) enddo enddo close (num_file3) filin = 'TABLE_met_V' open (UNIT=num_file3,file=filin) write(num_file3,*) geoplace,'time=',time write(num_file3,*) NX*NY,DX,NX,NY do j=1,ny do i=1,nx write(num_file3,*) i,j, METHOD2(i,j,2) enddo enddo close (num_file3) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!impression tuij!!!!!!!!!!!!!!!!! filin = 'TABLE_tuij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tuij(i,j)!/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tumimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tumimj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tumij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tumij(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tuimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tuimj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tumipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tumipj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tuimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tuimj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tuipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tuipj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tupimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tupimj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tupij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tupij(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tupipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tupipj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tvij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tvij(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tvmij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tvmij(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tvmipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tvmipj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_tvipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,tvipj(i,j)/tuij(i,j) enddo enddo close (num_file4) filin = 'TABLE_opposx' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do j=1,ny do i=1,nx write(num_file4,*) i,j,opposx(i,j)/tuij(i,j) enddo enddo close (num_file4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!impression svij!!!!!!!!!!!!!!!!! filin = 'TABLE_svij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svij(i,j)!/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svmimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svmimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svmipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svmipj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svipj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svpimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svpimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svpij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svpij(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_svpipj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,svpipj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_suij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,suij(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_suimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,suimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_supimj' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,supimj(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_supij' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,supij(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_opposy' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do j=1,ny do i=1,nx write(num_file4,*) i,j,opposy(i,j)/svij(i,j) enddo enddo close (num_file4) filin = 'TABLE_FRONT' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny IF (front(I,J).EQ.4) THEN INEARX=IFLOTMX(I-1,J) +IFLOTMX(I,J) +IFLOTMX(I+1,J) & +IFLOTMX(I-1,J-1)+IFLOTMX(I,J-1)+IFLOTMX(I+1,J-1) & +IFLOTMX(I-1,J+1)+IFLOTMX(I,J+1)+IFLOTMX(I+1,J+1) & +IFLOTMY(I-1,J) +IFLOTMY(I,J) +IFLOTMY(I+1,J) & +IFLOTMY(I-1,J-1)+IFLOTMY(I,J-1)+IFLOTMY(I+1,J-1) & +IFLOTMY(I-1,J+1)+IFLOTMY(I,J+1)+IFLOTMY(I+1,J+1) if (INEARX==0) then write(num_file4,*) i,j,' 6' else write(num_file4,*) i,j,front(i,j) endif ELSE write(num_file4,*) i,j,front(i,j) INEARX=1000 END IF enddo enddo close (num_file4) filin = 'sgbsv_u.dat' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,uxbar(i,j) enddo enddo close (num_file4) filin = 'sgbsv_v.dat' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,uybar(i,j) enddo enddo close (num_file4) filin = 'TABLE_h' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,h(i,j) enddo enddo close (num_file4) filin = 'TABLE_ICE' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,ICE(i,j) enddo enddo close (num_file4) filin = 'TABLE_index' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny if (okumat(i,j).or.okvmat(i,j)) then write(num_file4,*) i,j,ligu(i,j),ligv(i,j) endif enddo enddo close (num_file4) filin = 'TABLE_pvi' open(UNIT=num_file4,file=filin) write(num_file4,*) geoplace,'time=',time write(num_file4,*) NX*NY,DX,NX,NY do i=1,nx do j=1,ny write(num_file4,*) i,j,pvi(i,j) enddo enddo close (num_file4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END SUBROUTINE DEBUGAGE ! RETURN END SUBROUTINE REMPLIDOM