source: trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90 @ 10

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

!> \namespace equat_adv_diff_2D_vect
!!  Calcule l evolution de l epaisseur en resolvant 
!!  une equation d'advection diffusion
!! @note Version avec :  call resol_adv_diff_2D_vect (Dmx,Dmy,advmx,advmy,H_p,i_Hp,bilmass,vieuxH,H)
!! @note - le terme vecteur est passe en dummy parametres
!! l epaisseur peut etre imposee   
!! \author CatRitz
!! \date june 2009
!! @note  Used modules:
!! @note  - use module3D_phy
!! @note  - use reso_adv_diff_2D
!! @todo essayer dans le futur la methode de Richard dUX/dx = H dU/dx + U dH/dx
!<
module equat_adv_diff_2D_vect                          ! Cat nouvelle mouture juin 2009
use module3D_phy
use reso_adv_diff_2D_vect
use prescribe_H

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
real, dimension(nx,ny) :: bilmass         !< bilan de masse pour la colonne
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


contains

!> SUBROUTINE: init_icethick
!! Definition et initialisation des parametres qui gerent la repartition advection diffusion
!! @return adv_frac

subroutine init_icethick

implicit none

namelist/mass_conserv/adv_frac,V_limit  


rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,mass_conserv)


write(num_rep_42,*)'! Conservation de la masse avec equation advection-diffusion '
write(num_rep_42,mass_conserv)
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 advection'
write(num_rep_42,*)'! V_limit depend de la calotte : '
write(num_rep_42,*)'! typiquement 3000 m/an en Antarctique, 10000 m/an au Groenland'
write(num_rep_42,*)'!___________________________________________________________' 


call init_reso_adv_diff_2D     
call init_prescribe_H        ! initialize grounding line mask

Dx1=1./Dx
itour=0

end subroutine init_icethick

!------------------------------------------------------------------
!> SUBROUTINE: icethick3
!! Routine pour le calcule 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 (i_Hp(:,:).eq.1)                 ! previous prescribed thickness
   vieuxH(:,:)=Hp(:,:)                     ! H may have been changed by calving
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

! attention limitation ! 
uxbar(:,:)=max(uxbar(:,:),-V_limit)
uxbar(:,:)=min(uxbar(:,:),V_limit)
uybar(:,:)=max(uybar(:,:),-V_limit)
uybar(:,:)=min(uybar(:,:),V_limit)

! 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)


! 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.
if (itracebug.eq.1) call tracebug("aux tres petit,deuxieme appel a next_time dt*4")
      call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
   end if
end if


timemax=time

! etait avant dans step
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


if (geoplace(1:5).eq.'mism3') then       ! pas de flux sur les limites laterales
   dmy(:,2)       = 0.
   dmy(:,3)       = 0.
   dmy(:,ny-1)    = 0.
   dmy(:,ny)      = 0.
   advmy(:,2)     = 0.
   advmy(:,3)     = 0.
   advmy(:,ny-1)  = 0.
   advmy(:,ny)    = 0.
   uybar(:,2)     = 0
   uybar(:,3)     = 0
   uybar(:,ny-1)  = 0
   uybar(:,ny)    = 0   
end if


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


bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance

! diverses precription de l'epaisseur en fonction de l'objectif
!-------------------------------------------------------------------


  call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions)


   


 if ((igrdline.eq.1)) then               ! present grounding line

!    if ((time.lt.time_gl(1)).or.(nt.lt.2)) then      ATTENTION test seulement si version prescribe-H_mod.f90
       call prescribe_present_H_gl       ! prescribe present thickness at the grounding line
!    else
!       call prescribe_retreat
!    endif
 end if

 if ((igrdline.eq.2)) then               ! paleo grounding line
    call prescribe_paleo_gl_shelf
 end if

 if ((igrdline.eq.3)) then  
!    where (flot(:,:))
!       mk_gr(:,:) = 0
!    elsewhere
!       mk_gr(:,:) = 1
!    end where
    if (itracebug.eq.1) call tracebug(" avant time_step_recul")
    call time_step_recul                 ! version prescribe-H-i2s_mod.f90 avec use proto_recul_mod 
    if (itracebug.eq.1) call tracebug(" apres time_step_recul")
!    call prescr_ice2sea_retreat         ! version prescribe-H_mod.f90
 endif


!if (time.ge.t_break) then                ! action sur les ice shelves
!   call melt_ice_shelves                 ! ATTENTION version prescribe-H_mod.f90
!   call break_all_ice_shelves
!end if


!!$! impose une variation d'épaisseur d'apres une carte a un temps donne (a affiner plus tard pour des runs transient)
!!$if (time.eq.t_break) then
!!$   call lect_prescribed_deltaH
!!$else if (time.gt.t_break) then
!!$   call prescribe_deltaH 
!!$end if
!!$
!!$if (time.eq.t_break) then                            ! si appele apres lect_prescribed_deltaH, cumule les effets 
!!$   call break_all_ice_shelves
!!$end if

debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)

debug_3D(:,:,86)=Hp(:,:)
debug_3D(:,:,85)=i_Hp(:,:)

!!$where (i_hp(:,:).eq.1)
!!$   vieuxh(:,:)=Hp(:,:)
!!$endwhere


! Appel a la routine de resolution resol_adv_diff_2D_vect
!------------------------------------------------------------
! On ecrit les vitesses sous la forme
! Ux = -Dmx * sdx + advmx

! Dmx, Dmy sont les termes diffusifs de la vitesse
! advmx, advmy sont les termes advectifs.
! la repartition entre les deux depend de adv_frac.

! i_Hp est un tableau des points ou l'epaisseur est imposee a la valeur Hp
! bilmass est le bilan de masse (surface et fond)
! vieuxH est la precedente valeur de H
! H est le retour de la nouvelle valeur.


call resol_adv_diff_2D_vect (Dmx,Dmy,advmx,advmy,Hp,i_Hp,bilmass,vieuxH,H)

! 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


! eventuellement retirer apres spinup ou avoir un cas serac
Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt

if (igrdline.ne.3) then
   Hdot(:,:)=min(Hdot(:,:),10.)
   Hdot(:,:)=max(Hdot(:,:),-10.)
endif

where (i_hp(:,:).ne.1)
   H(:,:)=vieuxH(:,:)+Hdot(:,:)*dt
end where



! 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_vect