source: trunk/SOURCES/Old-sources/conserv-mass-adv-diff_mod-EGU-2009.f90 @ 334

!> \file conserv-mass-adv-diff_mod-EGU-2009.f90
!! Module qui calcule l evolution de l epaisseur en resolvant 
!! une equation d'advection diffusion
!<

!> \namespace equat_adv_diff_2D  
!! Module qui calcule l evolution de l epaisseur en resolvant 
!! une equation d'advection diffusion
!! \author Cat
!! \date Fevrier 2007
!! @note Used modules
!! @note   - use module3D_phy
!! @note   - use reso_adv_diff_2D
!<
module equat_adv_diff_2D                             ! Cat Fevrier 2007

use module3D_phy
use reso_adv_diff_2D


implicit none

real, dimension(nx,ny) :: advmx           !< partie advective en x
real, dimension(nx,ny) :: advmy           !< partie advective en y
real, dimension(nx,ny) :: dmx             !< partie advective en x
real, dimension(nx,ny) :: dmy             !< partie advective en y
real, dimension(nx,ny) :: VieuxH          !< ancienne valeur de H
logical, dimension(nx,ny) :: zonemx       !< pour separer advection-diffusion
logical,  dimension(nx,ny) :: zonemy      !< pour separer advection-diffusion
real :: adv_frac                          !< fraction du flux traitee en advection
integer :: itesti
integer :: itour

! test d'epaisseur prescrite
real,dimension(nx,ny) :: H_p      !< H value if prescribed
integer,dimension(nx,ny) :: i_Hp  !< 1 if H is prescribed on this node, else 0
real ::itest                      !< variable travail
real :: dh_gael                   !< variation d'epaisseur imposee
contains
!> SUBROUTINE: init_icethick  
!! Routine qui permet d'initialiser les parametres icethick
!! @note  definition du parametre qui gerent la repartition advection diffusion
!>
subroutine init_icethick


write(num_rep_42,*)'Conservation de la masse avec equation advection-diffusion '
write(num_rep_42,*)'la repartition depend de adv_frac'
write(num_rep_42,*)'>1  -> advection seule'
write(num_rep_42,*)'0   -> diffusion seule'
write(num_rep_42,*)'0<*<1   -> fraction de l advection'
write(num_rep_42,*) '-1 -> zones diffusion + zones advecttion'

!varname='adv_frac'

adv_frac=2. ! 0.8

!adv_frac=0
write(num_rep_42,*)'adv_frac=',adv_frac

 call init_reso_adv_diff_2D

 Dx1=1./Dx

itour=0


! lecture des points prescrits (à mettre ensuite dans les lectures standard)  16 avril 2009

!Version avec la ligne actuelle
goto 78
open(123,file="../INPUT/ANTEIS1/grounding-line-gael-40km.ijg")


i_HP(:,:)=0
H_p(:,:)=-9999.


do k=1,nx*ny
   read(123,*,end=77) i,j
   i_HP(i,j)=1
   H_p(i,j)=H0(i,j)
   if (H0(i,j) .lt.-Bsoc(i,j)*row/ro) then 
      write(6,*) i,j,H0(i,j),Bsoc(i,j),mk(i,j)
   endif
   H_p(i,j)=max(H0(i,j),-Bsoc(i,j)*row/ro+50.) ! pour etre sur que le point ne flotte pas.
end do
77 continue

close(123)

78 continue ! version avec la perturbation de Gael
open(123,file="../INPUT/ANTEIS1/delta-H-Gael-ij.ijz")

i_HP(:,:)=0
H_p(:,:)=-9999.

do k=1,nx*ny
   read(123,*,end=79) i,j,dH_gael   ! attention il faut enlever dh_gael
   i_HP(i,j)=1
   H_p(i,j)=H(i,j)-dH_gael          ! on enleve la perturbation Gael de la topo reprise dans le cptr
                                    ! il peut être devenu flottant ou pas
   H(i,j)=H_p(i,j)
end do
79 continue
close(123)

end subroutine init_icethick

!------------------------------------------------------------------
!> SUBROUTINE: icethick3 
!! Routine de calcul de l'epaisseur de glace 
!>
subroutine icethick3

  implicit none
  real,dimension(nx,ny) :: Dminx,Dminy
  real,dimension(nx,ny) :: Uxdiff,Uydiff           ! vitesse due a la diffusion
  integer :: it1,it2,it3

  real aux                                         ! pour le calcul du critere
  real maxdia                                      ! sur le pas de temps
  integer imaxdia,jmaxdia





  if (itracebug.eq.1)  call tracebug(' Entree dans routine icethick')


  ! Copie de H sur vieuxH
  ! --------------------
  where (mk0(:,:).eq.0)
     vieuxH(:,:)=0.
  elsewhere 
     vieuxH(:,:)=H(:,:)
  end where

  ! mk0 est le masque à ne jamais dépasser
  ! mk0=0  ->  pas de glace 
  ! mk0=-1 ->  on garde la même epaisseur
  ! pour l'Antarctique masque mko=1 partout


  Maxdia = -1.0
  imaxdia=0
  jmaxdia=0

  ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante
  do i=2,nx-1
     do j=2,ny-1
        aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) &
             + (abs(uybar(i,j)) + abs(uybar(i,j+1)))
        aux=aux*dx1*0.5
        if (aux.gt.maxdia) then
           maxdia=aux
           imaxdia=i
           jmaxdia=j
        endif
     end do
  end do


  if (maxdia.ge.(testdiag/dtmax)) then
     dt = testdiag/Maxdia
     dt=max(dt,dtmin)
  else
     dt = dtmax
  end if

  ! next_time ajuste le pas de temps pour faciliter la synchronisation 
  ! avec le pas de temps long (dtt)

  call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)


!!$write(166,*)'--------------------------------------------------------'
!!$write(166,*)'apres old calcul de dt',time,dt

  ! calcul diffusivite - vitesse 
  !----------------------------------------------------------------- 

  Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif
  Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar
  ! pour eviter des problemes dans le cas 0< adv_frac <1 
  dmx(:,:)=diffmx(:,:)
  dmy(:,:)=diffmy(:,:)

  ! schema amont pour la diffusion en x
  where (Uxbar(:,:).ge.0.)
     dmx(:,:)=diffmx(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=1)
     uxdiff(:,:)=min(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot
  elsewhere
     dmx(:,:)=diffmx(:,:)*H(:,:)
     uxdiff(:,:)=max(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot
  end where



  ! schema amont pour la diffusion en y
  where (uybar(:,:).ge.0.)
     dmy(:,:)=diffmy(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=2)
     uydiff(:,:)=min(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot
  elsewhere
     dmy(:,:)=dmy(:,:)*H(:,:)
     uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot
  end where



  ! tests pour la répartition advection - diffusion
  !------------------------------------------------

  if (adv_frac.gt.1.) then ! advection seulement
     advmx(:,:)=Uxbar(:,:)
     advmy(:,:)=Uybar(:,:)
     Dmx(:,:)=0.
     Dmy(:,:)=0.

  else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement
     advmx(:,:)=0.
     advmy(:,:)=0.

     ! D reste la valeur au dessus)

  else if ((adv_frac.ge.1.e-8).and.(adv_frac.le.1.)) then        ! advection-diffusion)

     ! selon x --------------------------------

     ! advection = adv_frac* tout ce qui n'est pas diffusion
     ! diffusion = ce qui peut être exprime en diffusion 
     !             + une partie (1-adv_frac) de l'advection

     where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general
        advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion
        advmx(:,:)=advmx(:,:)*adv_frac                         ! la fraction adv_frac du precedent
        Dminx(:,:)=-(Uxbar(:,:)-advmx(:,:))/sdx(:,:)          ! ce qui reste exprime en diffusivite
        ! a multiplier par H

        ! schema amont pour la diffusivite
        where (uxbar(:,:).ge.0.)
           dmx(:,:)=dminx(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=1)
        elsewhere
           dmx(:,:)=dminx(:,:)*H(:,:)
        end where


     elsewhere   ! zones trop plates ou trop fines
        Dmx(:,:)=0.
        advmx(:,:)=Uxbar(:,:)
     end where


     ! selon y --------------------------------

     ! advection = adv_frac* tout ce qui n'est pas diffusion
     ! diffusion = ce qui peut être exprime en diffusion 
     !             + une partie (1-adv_frac) de l'advection

     where ((abs(sdy(:,:)).gt.1.e-8).and.(Hmy(:,:).gt.2.))      ! cas general
        advmy(:,:)=(Uybar(:,:)-Uydiff(:,:))                   ! tout ce qui n'est pas diffusion
        advmy(:,:)=advmy(:,:)*adv_frac                         ! la fraction adv_frac du precedent
        Dminy(:,:)=-(Uybar(:,:)-advmy(:,:))/sdy(:,:)          ! ce qui reste exprime en diffusivite
        ! a multiplier par H

        ! schema amont pour la diffusivite
        where (uybar(:,:).ge.0.)
           dmy(:,:)=dminy(:,:)*eoshift(H(:,:),shift=-1,boundary=0.,dim=2)
        elsewhere
           dmy(:,:)=dminy(:,:)*H(:,:)
        end where



     elsewhere   ! zones trop plates ou trop fines
        Dmy(:,:)=0.
        advmy(:,:)=Uybar(:,:)
     end where


  else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs

     zonemx(:,:)=flgzmx(:,:)
     zonemy(:,:)=flgzmy(:,:)

     ! selon x --------------
     where (zonemx(:,:))
        advmx(:,:)=Uxbar(:,:)
        Dmx(:,:)=0.
     elsewhere
        advmx(:,:)=0.
     end where

     ! selon y --------------
     where (zonemy(:,:))
        advmy(:,:)=Uybar(:,:)
        Dmy(:,:)=0.
     elsewhere
        advmy(:,:)=0.
     end where

  end if


  where (flot(:,:) )
     tabtest(:,:)=1.
  elsewhere
     tabtest(:,:)=0.
  end where

  !------------------------------------------------------------------detect
!!$ tabtest(:,:)=dmx(:,:)
!!$
!!$call detect_assym(nx,ny,0,41,1,0,1,0,tabtest,itesti)
!!$
!!$if (itesti.gt.0) then
!!$   write(6,*) 'asymetrie sur dmx avant resol pour time=',time,itesti
!!$   stop
!!$else
!!$   write(6,*) ' pas d asymetrie sur dmx avant resol pour time=',time,itesti
!!$end if
  !----------------------------------------------------------- fin detect  


  !------------------------------------------------------------------detect
!!$ tabtest(:,:)=dmy(:,:)
!!$
!!$call detect_assym(nx,ny,0,41,1,0,1,1,tabtest,itesti)
!!$
!!$if (itesti.gt.0) then
!!$   write(6,*) 'asymetrie sur dmy avant resol pour time=',time,itesti
!!$   stop
!!$else
!!$   write(6,*) ' pas d asymetrie sur dmy avant resol pour time=',time,itesti
!!$end if
  !----------------------------------------------------------- fin detect  

  !------------------------------------------------------------------detect
!!$ tabtest(:,:)=advmx(:,:)
!!$
!!$call detect_assym(nx,ny,0,41,1,0,1,0,tabtest,itesti)
!!$
!!$if (itesti.gt.0) then
!!$   write(6,*) 'asymetrie sur advmx avant resol pour time=',time,itesti
!!$   stop
!!$else
!!$   write(6,*) ' pas d asymetrie sur advdmx avant resol pour time=',time,itesti
!!$end if
  !----------------------------------------------------------- fin detect  

  !------------------------------------------------------------------detect
!!$ tabtest(:,:)=advmy(:,:)
!!$
!!$call detect_assym(nx,ny,0,41,1,0,-1,1,tabtest,itesti)
!!$
!!$if (itesti.gt.0) then
!!$   write(6,*) 'asymetrie sur advmy avant resol pour time=',time,itesti
!!$   stop
!!$else
!!$   write(6,*) ' pas d asymetrie sur advmy avant resol pour time=',time,itesti
!!$end if
  !----------------------------------------------------------- fin detect  

  ! nouveau calcul de dt
  aux=maxval( (abs(advmx(:,:))+abs(advmy(:,:)))/dx+(abs(dmx(:,:))+abs(dmy(:,:)))/dx/dx)

  ! write(166,*) 'critere pour dt',time-dt,testdiag/aux,dt

  if (aux.gt.1.e-20) then
     if (testdiag/aux.lt.dt) then
        time=time-dt
        dt=testdiag/aux*4.
        call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
     end if
  end if


  !      print*,'Iteration =',nt,' temps = ',time-dt,' dt = ',dt

  ! vient de step
!!$     if (time.lt.timemax-100) then
!!$        print*,'le modele remonte le temps de',timemax,'a',time
!!$        stop
!!$     endif
  timemax=time

  if (time.lt.TGROUNDED) then
     MARINE=.false.
  else
     MARINE=.true.
  endif

  ! fin vient de step

  ! les variables dtdx et dtdx2 sont globales
  Dtdx2=Dt/(Dx**2)
  dtdx=dt/dx

!!$write(166,*) 'time=',time
!!$
!!$
!!$write(166,*) 'advmx,uxbar'
!!$do j=43,45
!!$   write(166,*) advmx(118:122,j)
!!$end do
!!$
!!$write(166,*) 'advmy,uybar'
!!$do j=43,45
!!$   write(166,*) advmy(118:122,j)
!!$end do
!!$
!!$write(166,*) 'vieuxH,hdot'
!!$do j=43,45
!!$   write(166,*) vieuxH(118:122,j)
!!$   write(166,*) Hdot(118:122,j)
!!$
!!$end do

  debug_3D(:,:,45)=dmx(:,:)
  debug_3D(:,:,46)=dmy(:,:)
  debug_3D(:,:,47)=advmx(:,:)
  debug_3D(:,:,48)=advmy(:,:)



  ! call resol_adv_diff_2D (Dmx,Dmy,advmx,advmy,vieuxH)

  ! test epaisseur prescrite : se rajoute au i_hp et h_p déja donnés. déclaré dans init
  where (flot(:,:))
     H_p(:,:)   = H0(:,:)
     i_hp(:,:) = 1
     !   elsewhere   
     !      i_hp(:,:) = 0
  endwhere



  call resol_adv_diff_2D_Hpresc (Dmx,Dmy,advmx,advmy,vieuxH,H_p,i_Hp)

  ! post traitements
  !--------------------

  ! epaisseur nulle sur les bords
  H(1,:)=0.
  H(nx,:)=0.
  H(:,1)=0.
  H(:,ny)=0.

  ! attention le bloc suivant sert a garder les ice-shelves stationnaires

  if (igrdline.eq.1) then
     where (flot(:,:))
        Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt
        H(:,:)=vieuxH(:,:)
     end where
  end if

  ! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord 
  where (H(:,:).lt.0.)
     ablbord(:,:)=H(:,:)/dt
  elsewhere
     ablbord(:,:)=0.
  endwhere

  H(:,:)=max(H(:,:),0.)       ! pas d'epaisseur negative


  ! calcul du masque "ice"
  where (flot(:,:))        ! points flottants, sera éventuellement réévalué dans flottab
     H(:,:)=max(H(:,:),1.) ! dans la partie marine l'épaisseur mini est 1 m
     where(H(:,:).gt.1.)
        ice(:,:)=1
     elsewhere
        ice(:,:)=0
     endwhere
  elsewhere            ! points posés
     where(H(:,:).gt.0.)
        ice(:,:)=1
     elsewhere
        ice(:,:)=0
     endwhere
  endwhere

  ! calcul de hdot. Quand igrdline=1, on garde hdot pour l'estimation du bmelt d'equilibre

  if (igrdline.eq.1) then
     where (.not.flot(:,:))
        Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt
     endwhere
  else
     Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt
  endif


  ! si mk0=-1, on garde l'epaisseur precedente
  ! en attendant un masque plus propre

  where(mk0(:,:).eq.-1)
     H(:,:)=vieuxH(:,:)
     Hdot(:,:)=0.
  end where



  !calul de l'ablation sur les bords (pourrait n'être appelé que pour les sorties courtes)
  if (isynchro.eq.1) call ablation_bord



  if (itracebug.eq.1)  call tracebug(' Fin routine icethick')! maxval(H)',maxval(H(:,:))
end subroutine icethick3

end module equat_adv_diff_2D