source: trunk/SOURCES/Hudson_files/output_hudson_mod.f90 @ 334

module  output_hudson

USE module3D_phy
use sedim_declar


implicit none

! 1_initialization
!------------------
INTEGER,parameter :: dtt_gliss=20 ! nbr de dtt pour calcul moy glissante
INTEGER,parameter :: mem=2000 ! nbr de dtt pour garder en memoire pour apres parametrer la fusion basale
REAL,dimension(dtt_gliss) ::  VOLtot, VOLshelf, oldVOLtot, oldVOLshelf 
REAL :: VOLtotgliss2 ! moyenne glissante au temps t avec des poids sur les 5 derniers pas de temps avant l'actuel 
REAL :: VOLtotgliss1 ! moyenne glissante au temps t-1 avec des poids sur les 5 derniers pas de temps avant l'actuel 
REAL :: VOLshelfgliss2 !  moyenne glissante au temps tavec des poids sur les 5 derniers pas de temps avant l'actuel 
REAL :: VOLshelfgliss1 !  moyenne glissante au temps t-1avec des poids sur les 5 derniers pas de temps avant l'actuel 
REAL :: deltaVOLtotdt ! variation de volume total sur moy glissante 
REAL :: deltaVOLshelfdt ! variation de volume flottant sur moy glissante 
REAL,DIMENSION(mem) ::  memdeltaVOLtot, memdeltaVOLshelf
integer ::  NPtot
integer ::  NPshelf
integer :: nptemp_tot
integer ::  NPtemp_sed
integer :: nbsed
REAL ::  Hmean_sed
REAL ::  tbmean_sed

! definition des points particuliers
integer,parameter  :: nbpoints=7
integer,parameter :: n1=6 ! nombre de sorties avant les points
integer,dimension(nbpoints) :: ipos,jpos
integer,parameter :: nbsortie=nbpoints*4+n1


! positions en km
real,dimension(nbpoints) :: xpos,ypos

! sorties pour ces points 
real,dimension(nbsortie) :: sortie


CONTAINS

subroutine init_outshort

!ndisp sorite courte tous les ndisp
NDISP=100
VOLtot(:)=0.
VOLshelf(:)=0.

end subroutine init_outshort
!VOLtot(:)=0.
!VOLshelf(:)=0.
!_________________________________________________________________________
subroutine shortoutput

oldVOLtot(:)=VOLtot(:)
oldVOLshelf(:)=VOLshelf(:)

!cdc if (itracebug.eq.1)  call tracebug(' Heino: entree dans routine output_heino')

! 1_initialization
!------------------
!VOLtot=0.
!VOLtot(:)=0.
!VOLshelf(:)=0.
!VOLtotgliss2=0.
!VOLtotgliss1=0.
!VOLshelfgliss2=0.
!VOLshelfgliss1=0.
NPtot=0
NPshelf=0
nptemp_tot=0
NPtemp_sed=0
Hmean_sed=0.
tbmean_sed=0.
nbsed=0.
! definition des positions des points particuliers en km
ypos(:)=2000.

xpos(1)=3900.
xpos(2)=3800.
xpos(3)=3700.
xpos(4)=3500.
xpos(5)=3200.
xpos(6)=2900.
xpos(7)=2600.


do k=1,nbpoints
   ipos(k)=nint(xpos(k)/dx*1000.)+1
   jpos(k)=nint(ypos(k)/dx*1000.)+1
end do



!oldVOLtot=eoshift(VOLtot,shift=1,boundary=0.)
!oldVOLshelf=eoshift(VOLshelf,shift=1,boundary=0.)

! 2_preparing outputs
!--------------------      
  do I=1,NX 
     do J=1,NY
        NPtot=NPtot+1   ! nombre de points
        VOLtot(dtt_gliss)=VOLtot(dtt_gliss)+H(I,J)  ! volume total
        if (flot(i,j))then
           NPshelf=NPshelf+1
           VOLshelf(dtt_gliss)=VOLshelf(dtt_gliss)+H(I,J)
        endif
        if (t(i,j,nz).ge.tpmp(i,j,nz)) then ! points base temperee
           nptemp_tot=nptemp_tot+1 ! 
           if (mksedim(i,j).ge.2) nptemp_sed=nptemp_sed+1 ! seulement ceux du sediment
        endif
        if (mksedim(i,j).ge.2) then
           Hmean_sed=Hmean_sed+H(i,j) ! epaisseur moyenne zone sediment
           tbmean_sed=tbmean_sed+(t(i,j,nz)-tpmp(i,j,nz))+273.15
           nbsed=nbsed+1
        endif
    end do
  end do

!VOLtot=eoshift(VOLtot,shift=1,boundary=0.)
!VOLshelf=eoshift(VOLshelf,shift=1,boundary=0.)



VOLtotgliss1=VOLtotgliss2

VOLtotgliss2=(0.1*VOLtot(1)+0.15*VOLtot(2)+0.2*VOLtot(3)+0.25*VOLtot(4)+0.3*VOLtot(5))/12+(0.1*VOLtot(6)+0.15*VOLtot(7)+0.2*VOLtot(8)+0.25*VOLtot(9)+0.3*VOLtot(10))/6+(0.1*VOLtot(11)+0.15*VOLtot(12)+0.2*VOLtot(13)+0.25*VOLtot(14)+VOLtot(15))/4+(0.1*VOLtot(16)+0.15*VOLtot(17)+0.2*VOLtot(18)+0.25*VOLtot(19)+0.3*VOLtot(20))/2

!VOLtotgliss2=0.1*VOLtot(1)+0.15*VOLtot(2)+0.2*VOLtot(3)+0.25*VOLtot(4)+0.3*VOLtot(5)
deltaVOLtotdt=(VOLtotgliss2-VOLtotgliss1)/dtt

memdeltaVOLtot(mem)=deltaVOLtotdt


VOLshelfgliss1=VOLshelfgliss2

VOLshelfgliss2=(0.1*VOLshelf(1)+0.15*VOLshelf(2)+0.2*VOLshelf(3)+0.25*VOLshelf(4)+0.3*VOLshelf(5))/12+(0.1*VOLshelf(6)+0.15*VOLshelf(7)+0.2*VOLshelf(8)+0.25*VOLshelf(9)+0.3*VOLshelf(10))/6+(0.1*VOLshelf(11)+0.15*VOLshelf(12)+0.2*VOLshelf(13)+0.25*VOLshelf(14)+VOLshelf(15))/4+(0.1*VOLshelf(16)+0.15*VOLshelf(17)+0.2*VOLshelf(18)+0.25*VOLshelf(19)+0.3*VOLshelf(20))/2

!VOLshelfgliss2=0.1*VOLshelf(1)+0.15*VOLshelf(2)+0.2*VOLshelf(3)+0.25*VOLshelf(4)+0.3*VOLshelf(5)
deltaVOLshelfdt=(VOLshelfgliss2-VOLshelfgliss1)/dtt

memdeltaVOLshelf(mem)=deltaVOLshelfdt

!debug
!print*
!print*,'output_hudson'
!print*,time,dtt,dt
!print*,VOLtot(:)
!print*,VOLshelf(:)
!print*,deltaVOLshelfdt, 'deltavolshelfdt'
!print*,memdeltaVOLshelf(:), 'memdeltavolshelf(:)'
!print*,memdeltaVOLshelf(1000), 'memdeltavolshelf(1000)'
!print*,VOLtotgliss1,VOLtotgliss2,VOLshelfgliss1,VOLshelfgliss2
!print*


nbsed=max(nbsed,1)  ! pour eviter une division par 0.

sortie(1)=time
sortie(2)=VOLtot(19)*dx*dx/1.e15
!sortie(3)=VOLtot(4)*dx*dx/1.e15
sortie(3)=VOLshelf(19)*dx*dx/1.e15
!sortie(4)=oldVOLtot(5)*dx*dx/1.e15
!sortie(5)=oldVOLtot(4)*dx*dx/1.e15
!sortie(5)=oldVOLshelf(5)*dx*dx/1.e15
!sortie(5)=VOLshelf(6)*dx*dx/1.e15
sortie(4)=deltaVOLtotdt*dx*dx/1.e15
sortie(5)=deltaVOLshelfdt*dx*dx/1.e15
!jalv sortie(3)=nptemp_tot*dx*dx/1.e12
!jalv sortie(4)=hmean_sed/nbsed/1.e3
!sortie(5)=nptemp_sed*dx*dx/1.e12
!sortie(6)=tbmean_sed/nbsed
sortie(6)=VOLtotgliss2*dx*dx/1.e15
l=n1+1
do k=1,7
   i1=ipos(k)
   j1=jpos(k)
!   sortie(l)=VOLtotgliss1*dx*dx/1.e15    !sortie 7
   sortie(l)=memdeltaVOLshelf(1000)    !sortie 7
   sortie(l+1)=VOLshelfgliss2*dx*dx/1.e15   !sortie 8
   sortie(l+2)=(memdeltaVOLshelf(1)+memdeltaVOLshelf(500)+memdeltaVOLshelf(1000)+memdeltaVOLshelf(1500))*dx*dx/4.e15    
!   sortie(l+2)=VOLshelfgliss1*dx*dx/1.e15   !sortie 9
!   sortie(l)=H(i1,j1)/1000.
!   sortie(l+1)=T(i1,j1,nz)-tpmp(i1,j1,nz)
!   sortie(l+2)=chalgliss_maj(i1,j1)/secyear
!   sortie(l+2)=rog*Hmx(i1,j1)*sdx(i1,j1)*uxbar(i1,j1)
!   sortie(l+3)=rog*Hmx(i1+1,j1)*sdx(i1+1,j1)*uxbar(i1+1,j1) 
   sortie(l+3)=0.5*(uxbar(i1,j1)+uxbar(i1+1,j1))
   l=l+4
end do

!decalage circulaire des variables contenant le volume de glace sur les derniers pas de temps
VOLtot=eoshift(VOLtot,shift=1,boundary=0.)
VOLshelf=eoshift(VOLshelf,shift=1,boundary=0.)

memdeltaVOLshelf=eoshift(memdeltaVOLshelf,shift=1,boundary=0.)


! 2_writing outputs
!------------------    
!     **** short display **** 


  write(num_ritz,fmt='(34(es14.6,1x))') (sortie(k),k=1,nbsortie)
  

end subroutine shortoutput

end module output_hudson