!module diagno_L2_mod ! Nouvelle version, compatible remplimat 2008 Cat module diagno_mod ! nom pendant les tests !$ USE OMP_LIB use module3D_phy use module_choix implicit none real :: somint,test,delp,prec real, dimension(nx,ny) :: uxb1 real, dimension(nx,ny) :: uyb1 integer, dimension(nx,ny) :: imx_diag integer, dimension(nx,ny) :: imy_diag integer :: nxd1,nxd2 ! domaine selon x Dans l'appel rempli_L2 integer :: nyd1,nyd2 ! domaine selon y integer :: itour_pvi integer :: ifail_diagno ! pour recuperation d'erreur integer :: iplus1,jplus1 integer :: ctvisco,iumax,jumax real :: delumax,errmax real :: phiphi,bt2,d02,discr,ttau real :: sf3,sf1,epsxxm,epsyym,epsm,sf01,sf03 ! pour le calcul de la viscosite real :: viscm real :: sf_shelf ! coef mult enhancement factor pour shelves logical :: stopvisco,viscolin logical :: test_visc contains !------------------------------------------------------------------------------------ subroutine init_diagno namelist/diagno_rheol/sf01,sf03,pvimin ! attribution des coefficients de viscosite ! formats pour les ecritures dans 42 428 format(A) ! lecture des parametres du run block draghwat !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,diagno_rheol) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&diagno_rheol ! nom du bloc diagno_rheol' write(num_rep_42,*) write(num_rep_42,*) 'sf01 = ',sf01 write(num_rep_42,*) 'sf03 = ',sf03 write(num_rep_42,*) 'pvimin = ',pvimin write(num_rep_42,*)'/' write(num_rep_42,428) '! coefficients par rapport a la loi glace posee ' write(num_rep_42,428) '! sf01 : coefficient viscosite loi lineaire ' write(num_rep_42,428) '! sf03 : coefficient viscosite loi n=3 ' write(num_rep_42,428) '! pvimin : valeur de pvi pour les noeuds fictifs ~ 1.e3' write(num_rep_42,428) '! tres petit par rapport aux valeurs standards ~ 1.e10' write(num_rep_42,*) ! Precision utilisee dans de calcul prec = 1.e-2 itour_pvi=1 ! si prend les valeurs analytiques dans le shelf if (geoplace(1:5).eq.'mism3') then sf_shelf = 1. itour_pvi= 0 else sf_shelf = 0.4 end if return end subroutine init_diagno !------------------------------------------------------------------------------------ subroutine diagnoshelf ! Resolution numerique des equations diagnostiques if (itracebug.eq.1) call tracebug(' Entree dans diagnoshelf') itour_pvi=itour_pvi+1 ! boucle sur la viscosite (pour l'instant pas actif) ! pvi(:,:)=0. Taushelf(:,:)=0. ! attention le bloc suivant est pour debug !!$gzmx(:,:)=.false. !!$gzmy(:,:)=.false. !!$ilemx(:,:)=.false. !!$ilemy(:,:)=.false. !!$flgzmx(:,:)=flotmx(:,:) !!$flgzmy(:,:)=flotmy(:,:) call dragging ! doit etre appele avant imx_imy if (itour_pvi.le.1) then call calc_pvi ! calcule les viscosites integrees !$OMP PARALLEL !$OMP WORKSHARE where (flot(:,:).and.(H.gt.1)) ! valeur analytique pour les shelfs pvi(:,:) = (4./coef_Sflot/rog)**2/btt(:,:,1,1)/H(:,:) end where !$OMP END WORKSHARE !$OMP END PARALLEL ! avec couplage thermomecanique ! write(166,*) ' apres call calc_pvi',itour_pvi else call calc_pvi end if call imx_imy_nx_ny ! pour rempli_L2 : calcule les masques imx et imy qui ! donnent les cas de conditions aux limites ! ! version pour travailler sur tout le domaine nx ny if (geoplace(1:5).eq.'mism3') call mismip_boundary_cond ! appel a la routine rempl_L2 -------------------domaine nx x ny ------------ ! ! pour tout le domaine nxd1=1 nxd2=nx nyd1=1 nyd2=ny !call rempli_L2(1,nx,1,ny,uxbar,uybar,uxb1,uyb1,imx_diag,imy_diag,ifail_diagno) !nxd1=15 !nxd2=19 !nyd1=30 !nyd2=34 !nxd1=35 !nxd2=60 !nyd1=35 !nyd2=60 call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2), & imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno) ! Dans rempli_L2 ! uxprex(n1,n2) ! uyprec(n1,n2) vitesses de l'iteration precedente ! ! uxnew(n1,n2) ! uynew(n1,n2) uynew resultat de cette iteration ! ! imx(n1,n2) masque pour imposer les vitesses ou leur dérivee ! imy(n1,n2) masque pour imposer les vitesses ou leur dérivée ! eventuellement le domaine n1,n2 peut etre un sous-domaine de nx,ny ! attention il faudra alors appeler avec des sous-tableaux ! Dans l'appel uxbar -> uxprec et Uxb1 -> Uxnew !--------------------------------------------------------------------------- if (ifail_diagno.gt.0) then write(6,*) ' Probleme resolution systeme L2. ifail=',ifail_diagno STOP endif ! nouvelles vitesses !$OMP PARALLEL !$OMP WORKSHARE uxbar(:,:)=uxb1(:,:) uybar(:,:)=uyb1(:,:) !$OMP END WORKSHARE ! calcul de tobmx et tobmy (frottement basal) apres calcul des vitesses ! --------------------------------------------------------------------- !$OMP DO do j=1,ny do i=1,nx tobmx(i,j)=-betamx(i,j)*uxbar(i,j) tobmy(i,j)=-betamy(i,j)*uybar(i,j) enddo enddo !$OMP END DO ! Mise ne réserve des vitesses flgzmx et flgzmy !$OMP WORKSHARE where (flgzmx(:,:)) uxflgz(:,:)=uxbar(:,:) elsewhere uxflgz(:,:)=0. endwhere where (flgzmy(:,:)) uyflgz(:,:)=uybar(:,:) elsewhere uyflgz(:,:)=0. endwhere !$OMP END WORKSHARE !$OMP END PARALLEL !i=92 !j=152 !write(6,*) 'time',time, uxbar(92,152),gzmx(92,152),ilemx(92,152),flotmx(92,152), flgzmx(92,152) return end subroutine diagnoshelf !------------------------------------------------------------------- subroutine calc_pvi ! calcule les viscosites integrees pvi et pvm ! loi polynomiale + couplage thermomécanique ! ! Attention ne marche que si la loi est la loi en n=3 + n=1 ! y compris le pur glen (n=3) ou le pur Newtonien (n=1) ! -------------------------------------------------------------------- ! les deformations sont supposées indépendantes de la profondeur ! et sont calculées dans strain_rate (appelé par main) ! eps(i,j) -> eps ! ttau -> tau (2eme invariant du deviateur des contraintes) ! BT2 loi en n=3, phiphi loi en n=1 calculés dans flowlaw ! ! La viscosité est calculée partout y compris pour la glace posée (ou elle est moins ! précise qu'un calcul direct avec la loi de déformation.) ! Pour les noeuds posés mais ayant un voisin stream ou flottant, on calcule ! la viscosité avec stream/shelves ! le calcul se fait sur les noeuds majeurs !$ integer :: rang ,nb_taches !$ logical :: paral integer :: t1,t2,ir real :: temps, t_cpu_0, t_cpu_1, t_cpu, norme if (itracebug.eq.1) call tracebug(' Calc pvi') !$OMP PARALLEL PRIVATE(rang,iplus1,jplus1,sf3,sf1,BT2,phiphi,ttau,d02,discr) !$ paral = OMP_IN_PARALLEL() !$ rang=OMP_GET_THREAD_NUM() !$ nb_taches=OMP_GET_NUM_THREADS() !$OMP WORKSHARE pvi(:,:) = pvimin Abar(:,:) = 0. !$OMP END WORKSHARE !$OMP DO do j=1,ny do i=1,nx iplus1=min(i+1,nx) jplus1=min(j+1,ny) if (flot(i,j)) then ! noeuds flottants sf3=sf03*sf_shelf sf1=sf01*sf_shelf else if (gzmx(i,j).or.gzmx(iplus1,j).or.gzmy(i,j).or.gzmy(i,jplus1)) then sf1=sf01 sf3=max(sf03,0.01) ! pour les fleuves de glace, un peu de Glen else if (ilemx(i,j).or.ilemx(iplus1,j).or.ilemy(i,j).or.ilemy(i,jplus1)) then sf1=sf01 sf3=max(sf03,0.01) ! pour les iles aussi else ! sf1=1 ! sf3=1 sf1=sf01 ! pour la viscosite anisotrope (ici en longitudinal) sf3=sf03 endif do k=1,nz BT2=BTT(i,j,k,1)*sf3 ! changement du sf phiphi=BTT(i,j,k,2)*sf1 ! changement du sf if (BT2.lt.1.e-25) then ! pur newtonien visc(i,j,k)=1./phiphi ttau=2.*visc(i,j,k)*eps(i,j) else ! polynomial ! en mettant Bt2 en facteur d02=eps(i,j) discr=((phiphi/3.)**3.)/Bt2+d02**2 discr=discr**0.5 ttau=(d02+discr)**(1./3.)-(discr-d02)**(1./3.) ttau=ttau*Bt2**(-1./3.) visc(i,j,k)=Bt2*ttau*ttau+phiphi if (visc(i,j,k).gt.1.e-15) then visc(i,j,k)=1./visc(i,j,k) else visc(i,j,k)=1.e15 endif endif pvi(i,j)=pvi(i,j)+visc(i,j,k) Abar(i,j) =(Bt2/2.)**(-1./3.) + Abar(i,j) Taushelf(i,j)=Taushelf(i,j)+ttau end do pvi(i,j) = pvi(i,j)*H(i,j)/nz Abar(i,j) = (Abar(i,j) /nz)**(-3.) Taushelf(i,j)=Taushelf(i,j)/nz end do end do !$OMP END DO ! cas des noeuds fictifs, si l'épaisseur est très petite ! pvimin est très petit !$OMP WORKSHARE where (H(:,:).le.1.) pvi(:,:) = pvimin end where where (ramollo(:,:).ge..5) pvi(:,:) = pvimin end where debug_3D(:,:,27)=pvi(:,:) !$OMP END WORKSHARE ! attention run 35 !-------------------- !!$if (time.gt.10.) then !!$ where (flot(:,:)) !!$ pvi(:,:)=pvimin !!$ end where !!$end if ! calcul de la viscosite integree au milieu des mailles (pvm) !$OMP DO do i=2,nx do j=2,ny ! les lignes suivantes pour un pvm moyenne des pvi pvm(i,j)=0.25*((pvi(i,j)+pvi(i-1,j-1))+ & (pvi(i,j-1)+pvi(i-1,j))) end do end do !$OMP END DO !$OMP END PARALLEL end subroutine calc_pvi !------------------------------------------------------------------ subroutine imx_imy_nx_ny ! definition des masques ! pour rempli_L2 : calcule les masques imx et imy qui ! donnent les cas de conditions aux limites ! version pour travailler sur tout le domaine nx ny !---------------------------------------------------- ! -34 -3 Nord -23 ! !----------------------------------------! ! ! ! ! ! 1 (prescrite) ! ! -4 ! ou ! -2 ! West! 2 (L2) ! Est ! ! ! ! ! ! ! !----------------------------------------! ! -41 -1 Sud -12 !$OMP PARALLEL !$OMP WORKSHARE imx_diag(:,:)=0 imy_diag(:,:)=0 ! a l'interieur du domaine !------------------------- where (flgzmx(:,:)) imx_diag(:,:)=2 ! shelf ou stream elsewhere imx_diag(:,:)=1 ! vitesse imposee end where where (flgzmy(:,:)) imy_diag(:,:)=2 ! shelf ou stream elsewhere imy_diag(:,:)=1 ! vitesse imposee end where ! bord sud imx_diag(:,1)=-1 imy_diag(:,2)=-1 ! bord nord imx_diag(:,ny)=-3 imy_diag(:,ny)=-3 ! bord Est imx_diag(1,:)=0 ! hors domaine a cause des mailles alternees imx_diag(2,:)=-4 imy_diag(1,:)=-4 ! bord West imx_diag(nx,:)=-2 imy_diag(nx,:)=-2 ! Coins imx_diag(2,1)=-41 ! SW imy_diag(1,2)=-41 imx_diag(nx,1)=-12 ! SE imy_diag(nx,2)=-12 imx_diag(nx,ny)=-23 ! NE imy_diag(nx,ny)=-23 imx_diag(2,ny)=-34 ! NW imy_diag(1,ny)=-34 ! hors domaine imx_diag(1,:)=0 ! hors domaine a cause des mailles alternees imy_diag(:,1)=0 ! hors domaine a cause des mailles alternees !$OMP END WORKSHARE !$OMP END PARALLEL end subroutine imx_imy_nx_ny !___________________________________________________________________________ ! pour imposer les conditions de mismip sur les bords du fleuve ! a appeler apres imx_imy_nx_ny subroutine mismip_boundary_cond if (itracebug.eq.1) call tracebug(' Subroutine mismip_boundray_cond') ! Condition pas de flux sur les bords nord et sud imy_diag(:,2) = 1 imy_diag(:,3) = 1 imy_diag(:,ny-1) = 1 imy_diag(:,ny) = 1 Uybar(:,2) = 0. Uybar(:,3) = 0. Uybar(:,ny-1) = 0. Uybar(:,ny) = 0. ! condition pas de cisaillement sur les bords nord et sud imx_diag(:,2) = -1 imx_diag(:,ny-1) = -3 ! coins imx_diag(2,2) = -41 imx_diag(nx,2) = -12 imx_diag(2,ny-1) = -34 imx_diag(nx,ny-1) = -23 imx_diag(1,:) = 0 ! ces points sont hors grille end subroutine mismip_boundary_cond !end module diagno_L2_mod end module diagno_mod