source: trunk/SOURCES/relaxation_mod-0.3.f90 @ 23

!> \file relaxation_mod-0.3.f90
!! Module pour la resolution de l'equation de relaxation.
!<

!> \namespace RELAXATION_MOD
!! Module pour la resolution de l'equation de relaxation.
!! \author Christophe DUMAS
!! \date fevrier 2000
!! @note Cette routine est appellee par ICETHICK3
!! @note Modele couple "ice sheet - ice shelves"
!<

MODULE RELAXATION_MOD

CONTAINS
!> SUBROUTINE: relaxation
!! Routine pour la methode de relaxation
!! \param NXX     [in]         defini la taille des tableaux
!! \param NYY     [in]         defini la taille des tableaux 
!! \param Mk      [in]         masks (ice sheet, max, above water, below water, 1)
!! \param MARINE
!! \param FLOT
!! \param DT      [in]         pas de temps court 
!! \param DX      [in]         pas en x
!! \param vieuxH  [in]         H au pas de temps precedent
!! \param BM      [in]         mass balance   'o'
!! \param BMELT   [in]         basal melting  'o'
!! \param UXBAR   [in]         vertically integrated velocity '>'
!! \param UYBAR   [in]         vertically integrated velocity '^'
!! \param H       [out]        ice thickness  'o'
!! \param HDOT    [out]        ice thickness derivee / t  'o'
!>
  subroutine RELAXATION(NXX,NYY,MK,MARINE,FLOT,DT,DX,vieuxH,BM,BMELT,UXBAR,UYBAR,H,HDOT)

    IMPLICIT NONE


! declaration des variables en entree
    INTEGER, intent(in) :: NXX, NYY !< defini la taille des tableaux
    INTEGER,dimension(NXX,NYY), intent(in) :: MK     !< masks (ice sheet, max, above water, below water, 1)
    LOGICAL, intent(in) :: MARINE
    LOGICAL,dimension(NXX,NYY), intent(in) :: FLOT
    REAL, intent(in) ::  DT                          !< pas de temps court
    REAL, intent(in) ::  DX                          !< pas en x
    REAL,dimension(NXX,NYY), intent(in) :: vieuxH    !< H au pas de temps precedent
    REAL,dimension(NXX,NYY), intent(in) :: BM        !< mass balance   'o'
    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
    REAL,dimension(NXX,NYY), intent(in) :: UXBAR     !< vertically integrated velocity '>'
    REAL,dimension(NXX,NYY), intent(in) :: UYBAR     !< vertically integrated velocity '^'


! declaration des variables en sortie
    REAL,dimension(NXX,NYY), intent(out) :: H      ! ice thickness  'o'
    REAL,dimension(NXX,NYY), intent(out) :: HDOT     ! ice thickness derivee / t  'o'


! declaration des variables locales
    INTEGER :: I,J
    REAL  :: TESTH 
    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
    REAL,dimension(NXX,NYY) :: DELTAH
    REAL :: RESTE,DELH,VH
    INTEGER :: ntour
    integer :: mbord
    REAL  :: DTSDX 
    LOGICAL :: STOPP


!   print*,'dans relaxation mod'
    H(:,:)= vieuxH(:,:)

!   attribution des coefficients  arelax ....
!   ----------------------------------------

    dtsdx=dt/DX
  
    do J=2,NYY-1
       do I=2,NXX-1


             FRELAX(I,J)= VIEUXH(I,J)+(BM(I,J)-BMELT(I,J))*DT
!          selon x
!                                     ecoulement vers la droite
             if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).ge.0.)) then
                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
                BRELAX(I,J)=0.
                CRELAX(I,J)=UXBAR(I+1,J)*dtsdx
              
!                                     ecoulement vers la gauche
             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).le.0.)) then
                ARELAX(I,J)=0.
                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
                CRELAX(I,J)=-UXBAR(I,J)*dtsdx
              
!                                     ecoulement divergent
             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).ge.0.)) then
                ARELAX(I,J)=0.
                BRELAX(I,J)=0.
                CRELAX(I,J)=(UXBAR(I+1,J)-UXBAR(I,J))*dtsdx
              
!                                     ecoulement convergent
             else if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).le.0.)) then
                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
                CRELAX(I,J)=0.
             endif

!          selon y 
!                                 ecoulement vers le haut (sens axe y)
             if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).ge.0.)) then
                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
                ERELAX(I,J)=0.
                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)*dtsdx))
                
!                                     ecoulement vers le bas
             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).le.0.)) then
                DRELAX(I,J)=0.
                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
                CRELAX(I,J)=1.+(CRELAX(I,J)-(UYBAR(I,J)*dtsdx))
           
!                                     ecoulement divergent
             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).ge.0.)) then
                DRELAX(I,J)=0.
                ERELAX(I,J)=0.
                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)-UYBAR(I,J))*dtsdx)
           
!                                     ecoulement convergent
             else if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).le.0.)) then
                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
                CRELAX(I,J)=1.+CRELAX(I,J)
           
             endif
!          endif
       end do
    end do 

! Boucle de relaxation :
! ----------------------

    STOPP = .false.
    dtsdx=dt/dx
    ntour=0

    Do 362 WHILE(.NOT.STOPP)
       ntour=ntour+1
     
       do j=2,NYY-1
          do i=2,NXX-1
            
             RESTE = (((ARELAX(I,J)*H(I-1,J) + BRELAX(I,J)*H(I+1,J)) &
                       + (DRELAX(I,J)*H(I,J-1) + ERELAX(I,J)*H(I,J+1))) &
                       + CRELAX(I,J)*H(I,J))- FRELAX(I,J)
     
             DELTAH(I,J) = RESTE/CRELAX(I,J)
!           H(I,J) = H(I,J) - DELTAH(I,J)
          end do
       end do

       do j=2,NYY-1
          do i=2,NXX-1
             H(I,J) = H(I,J) - DELTAH(I,J)
          end do
       end do
         
! critere d'arret:
! ----------------         

       Delh=0
       Vh=0

       DO j=2,NYY-1
          DO i=2,NXX-1
!   write(6,*) I,J,delh,deltah(i,j)
             Delh=Delh+deltah(i,j)**2
!            Vh=Vh+h(i,j)**2.
          END DO
       END DO

!      testh=SQRT(Delh/Vh)
       if (delh.gt.0.) then
          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
       else
          testh=0.
       endif
       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
     
     
362  Continue

! Conditions aux limites:
! -----------------------

!   ************ valeurs de H pour les coins ********

       DO i=2,NXX-1
          H(i,1) = 2.*H(i,2) - H(i,3)
          H(i,NYY) = 2.*H(i,NYY-1) - H(i,NYY-2)
       END DO
       
       DO J=2,NYY-1
          H(1,J) = 2.*H(2,j) - H(3,j)
          H(NXX,J) = 2.*H(NXX-1,j) - H(NXX-2,j)
       END DO

     
       h(1,1)=(h(1,2)+h(2,1))/2.
       Hdot(1,1)=0.
       h(1,NYY)=(h(1,NYY-1)+h(2,NYY))/2.
       hdot(1,NYY)=0.
       h(NXX,1)=(h(NXX,2)+h(NXX-1,1))/2.
       hdot(nxx,1)=0.
       h(NXX,NYY)=(h(NXX,NYY-1)+h(NXX-1,NYY))/2.
       hdot(nxx,nyy)=0.
! =========================================================================


!     open(87,file='test-H')
!     write(87,*)'time=',time

       do J=2,NYY-1
          do I=2,NXX-1
             HDOT(I,J)=(H(i,j)-vieuxH(i,j))/DT 



!         if (abs(HDOT(I,J)).gt.hdotmax) then
!  hdotmax=hdot(i,j)
!  idotmax=i
! jdotmax=j
! end if

!         write(87,*)i,j,h(i,j),hdot(i,j)
          end do
       end do


!         close(87) 


  end subroutine relaxation

end module relaxation_mod