source: trunk/SOURCES/Ant40_files/lect-anteis_mod.f90 @ 4

!> \file lect-anteis_mod.f90
!!Module pour la lecture de la topographie
!<

!> \namespace lect_topo_anteis
!! Module pour la lecture de la topographie
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<
module lect_topo_anteis

  use module3D_phy
  
character(len=80) :: filin
real,dimension(nx,ny) ::  xcc , ycc      !< coordeonnes en m
real, dimension(nx,ny,5) :: bidon          !< pour l'appel a courbure
contains
 
  subroutine input_topo

    !====================================== La reponse est 42 ===========
    write(num_rep_42,*)
    write(num_rep_42,*)' Fichiers en entree'
    write(num_rep_42,*)'----------------------'
    !==================================================================== 

    ! lecture adaptee aux fichiers intercomparaison EISMINT
    nxx=nx
    nyy=ny
    ! socle
    !       filin=TRIM(DIRNAMEINP)//'bedelev-2000-40km.dat'
    filin='bedelev-2000-40km.dat'
    call lect_eis(nx,ny,BSOC,filin,DIRNAMEINP)
    write(num_rep_42,*) 'fichier BSOC (socle) : ', filin 
    Bsoc0(:,:)=Bsoc(:,:)
    !  surface
    filin='surface-2000-40km.dat'
    call lect_eis(nx,ny,S,filin,DIRNAMEINP)
    write(num_rep_42,*) 'fichier surface : ', filin 

    !     epaisseur
    filin='icethic-2000-40km.dat'
    call lect_eis(nx,ny,H,filin,DIRNAMEINP)
    write(num_rep_42,*) 'fichier epaisseur : ', filin 
    H0(:,:)=H(:,:)

!!$    !     coupure du shelf
!!$    filin=TRIM(DIRNAMEINP)//'coupe.ijz'
!!$    open(11,file=filin)
!!$    !     distcent devient la condition de coupure
!!$    do k=1,nx*ny
!!$       read(11,*) i,j,distcent(i,j)
!!$    end do
!!$    close(11)
!!$    write(num_rep_42,*) 'coupure du shelf (pas activé ?) : ', filin 

    !     masque
    filin='maskHUY40km.dat'
    call lect_ieis(nx,ny,MK0,filin,DIRNAMEINP)
    write(num_rep_42,*) 'fichier masque : ', filin 

!!$
!!$    !------------------------------------------------------- juillet 2007 ----
!!$    ! corrections du socle pour favoriser les fleuves . A mettre aussi apres lect cptr
!!$
!!$    filin=DIRNAMEINP//'corrections-bsoc.dat'
!!$
!!$    write(num_rep_42,*) 'correction du socle dans',filin
!!$
!!$    open(777,file=filin)
!!$    read(777,*,end=778)
!!$    read(777,*,end=778)
!!$    do k=1,nx*ny
!!$       read(777,*,end=778) i,j,bsoc(i,j)
!!$       H(i,j)=S(i,j)-Bsoc(i,j)
!!$       B(i,j)=S(i,j)-H(i,j)
!!$       write(166,*) i,j,bsoc(i,j)
!!$    end do
!!$778 continue
!!$    write(num_rep_42,*) k,'points corriges'
!!$    close(777)
    !------------------------------------------------------------------------


    !     enlever les epaisseurs fictives de Philippe
    do I=1,NX
       do J=1,NY
          !     attention les valeurs de masque de Philippe sont differentes
          if ((MK0(I,J).gt.1).and.(H(i,j).le.218.00)) then ! mer libre
             H(i,j)=1. 
             S(I,J)=H(i,j)-(RO/ROW)*H(i,j)
          endif
       end do
    end do
    S0(:,:)=S(:,:)



    ! pour l'Antarctique masque mko vrai partout (version 2006)
    MK0(:,:)=1

    ! determination des flot
    do J=1,NY
       do I=1,NX
          if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then
             FLOT(I,J)=.TRUE.
          else
             FLOT(I,J)=.FALSE.
          endif
       enddo
    enddo


    ! calcul des courbures du socle

    call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), &
         bidon(:,:,4),socle_cry,bidon(:,:,5))

    socle_cry(:,:)=socle_cry(:,:)*dx*dx

    ! lecture des coordonnées geographiques

    filin=TRIM(DIRNAMEINP)//'coord-Ant-40km.dat'

    ! les coordonnees sont calculees en °dec avec GMT,
    ! les longitudes sont comprises entre -180 et +180 (negative a l'Ouest de
    ! Greenwich et positive a l'Est)
    open(unit=2004,file=filin,iostat=ios)
    do k=1,nx*ny
       read(2004,*) i,j,XCC(i,j),YCC(i,j),XLONG(i,j),YLAT(i,j)
    enddo
    close(2004)
    write(num_rep_42,*) 'fichier grille: ', filin 

    xmin=xcc(1,1)/1000.
    ymin=ycc(1,1)/1000.
    xmax=xcc(nx,ny)/1000.
    ymax=ycc(nx,ny)/1000.

    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel
    open(88,file=TRIM(DIRNAMEINP)//'grzone56-k000.pl')

    read(88,*)
    read(88,*)
    read(88,*)
    do j=1,ny
       do i=1,nx
          read(88,*) S_sealev(i,j),H_sealev(i,j),B_sealev(i,j),M_sealev(i,j)
       enddo
    enddo
    close(88)

    write(num_rep_42,*) 'fichier reference pour le niveau des mers : ', & 
         TRIM(DIRNAMEINP)//'grzone56-k000.pl'

    ! lecture du flux geothermique de Shapiro
    open(88,file=TRIM(DIRNAMEINP)//'ijphi-40km-ant.dat')

    write(num_rep_42,*) 'flux geothermique Shapiro : ', TRIM(DIRNAMEINP)//'ijphi-40km-ant.dat'

    do k=1,nx*ny
       read(88,*) i,j,ghf(i,j)
    end do

    ! pour passer les flux des mW/m2 au J/m2/an      
    ghf(:,:)=-SECYEAR/1000.*ghf(:,:)




  end subroutine input_topo


end module lect_topo_anteis