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

!> \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
  use interface_input
  use io_netcdf_grisli

  implicit none

  character(len=100) :: topo_dep       ! Topo de départ
  character(len=100) :: topo_ref       ! Topo de référence
  character(len=100) :: grid_topo      ! fichier grille
  character(len=100) :: ghf_fich       ! fichier grille
  character(len=100) :: filin
  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

!  character(len=100) :: file1
!  character(len=100) :: file2

  real, dimension(nx,ny,5) :: bidon          ! pour l'appel a courbure
  real :: sealevel0

contains
 
  subroutine input_topo

    integer :: ios
    ! pour les lectures ncdf
    real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer

    namelist/topo_file/topo_ref,topo_dep,grid_topo,ghf_fich
    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
    read(num_param,topo_file)
    ! formats pour les ecritures dans 42
428 format(A)
    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428) '&topo_file                                  !  input_topo  '
    write(num_rep_42,'(A,A)') 'topo_ref = ', topo_ref
    write(num_rep_42,'(A,A)') 'topo_dep = ', topo_dep
    write(num_rep_42,'(A,A)') 'grid_topo =', grid_topo
    write(num_rep_42,'(A,A)') 'ghf_fich = ', ghf_fich
    write(num_rep_42,*)'/'                      
    write(num_rep_42,428) '! topo_ref= topo ref isostasie'
    write(num_rep_42,428) '! topo_dep= topo de depart'
    write(num_rep_42,428) '! grid_topo : fichier i,j,x,y,lon,lat'
    write(num_rep_42,428) '! ghf_fich  : fichier flux geothermique'  
    write(num_rep_42,*)

    topo_ref=trim(dirnameinp)//trim(topo_ref) 
    topo_dep=trim(dirnameinp)//trim(topo_dep) 
    grid_topo=trim(dirnameinp)//trim(grid_topo)
    ghf_fich=trim(dirnameinp)//trim(ghf_fich)


    ! lecture adaptee aux fichiers intercomparaison EISMINT
    nxx=nx
    nyy=ny

!!$!     lecture de la topo actuelle
!!$!     ---------------------------
!!$     open (20,file=TRIM(DIRNAMEINP)//file2,status='old')
!!$       
!!$     read(20,'(A80)') TITRE
!!$     read(20,*) NI,NJ,NXX,NYY,STEP
!!$     read(20,*)
!!$         do J=1,ny 
!!$          do I=1,nx
!!$             read (20,*)  S0(I,J),H0(I,J),BSOC0(I,J)
!!$             S0(i,j)=max(S0(i,j),0.)
!!$          end do
!!$        end do
!!$     close(20)
!!$
!!$     
!!$!     lecture de la topo de depart
!!$!     ---------------------------
!!$     open (20,file=TRIM(DIRNAMEINP)//file1,status='old')
!!$!         open (20,file='../INPUT-DATA/hemin.g50')
!!$     read(20,'(A80)') TITRE
!!$     read(20,*) NI,NJ,NXX,NYY,STEP
!!$     read(20,*)
!!$         do J=1,ny 
!!$          do I=1,nx
!!$             read (20,*) S(I,J),H(I,J),BSOC(I,J)
!!$          end do
!!$        end do
!!$     close(20)

    !     lecture de la topo de référence
    !     -------------------------------
    ! Cette topo sert a calculer le socle de reference pour l'isostasie
    ! voir init_iso et a avoir une surface de reference pour les temperatures
    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd
!    call lect_input(1,'Bsoc',1,Bsoc0,topo_ref,file_ncdf)    ! socle
!    call lect_input(1,'S',1,S0,topo_ref,file_ncdf)          ! surface
!    call lect_input(1,'H',1,H0,topo_ref,file_ncdf)          ! epaisseur

! lecture pour eviter plantage avec compile -O0
     call Read_Ncdf_var('Bsoc',topo_ref,tab)
     Bsoc0(:,:)  = tab(:,:)
     call Read_Ncdf_var('S',topo_ref,tab)
     S0(:,:)  = tab(:,:)
     call Read_Ncdf_var('H',topo_ref,tab)
     H0(:,:)  = tab(:,:)
!~ !cdc correction de 3 pts qui posent probleme (pente tres forte):
!~      S0(38,53)=1400.
!~      Bsoc0(38,53)=-1000.
!~      H0(38,53)=S0(38,53)-Bsoc0(38,53)
!~      
!~      S0(35,53)=1300.
!~      Bsoc0(35,53)=-1000.
!~      H0(35,53)=S0(35,53)-Bsoc0(35,53)
!~      
!~      S0(35,56)=1200.
!~      Bsoc0(35,56)=-1000.
!~      H0(35,56)=S0(35,56)-Bsoc0(35,56)

!cdc correction point pole sud :
!    S0(71,71)=(S0(71,70)+S0(71,72)+S0(70,71)+S0(72,71))/4.
!    Bsoc0(71,71)=(Bsoc0(71,70)+Bsoc0(71,72)+Bsoc0(70,71)+Bsoc0(72,71))/4.


!    where (S0(:,:).GT.0)
!       H0(:,:)=S0(:,:)-BSOC0(:,:)
!    elsewhere
!       H0(:,:)=1.
!    endwhere


    where (BSOC(:,:).LT.-9999.)
       BSOC(:,:)=-9999.
    endwhere

    sealevel0=0.  ! voir a passer dans le fichier parametre
    S0(:,:)=max(S0(:,:),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
    mk0(:,:) = 1                   ! mk0=0 pour les zones interdites

    !     lecture de la topo de depart
    !     ---------------------------
    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd
!    call lect_input(1,'Bsoc',1,Bsoc,topo_dep,file_ncdf)    ! socle
!    call lect_input(1,'S',1,S,topo_dep,file_ncdf)          ! surface
!    call lect_input(1,'H',1,H,topo_dep,file_ncdf)          ! epaisseur
    ! lecture pour eviter plantage avec compile -O0
     call Read_Ncdf_var('Bsoc',topo_dep,tab)
     Bsoc(:,:)  = tab(:,:)
     call Read_Ncdf_var('S',topo_dep,tab)
     S(:,:)  = tab(:,:)
     call Read_Ncdf_var('H',topo_dep,tab)
     H(:,:)  = tab(:,:)
     
!~ !cdc correction de 3 pts qui posent probleme (pente tres forte):
!~      S(39,54)=1400.
!~      Bsoc(39,54)=-1000.
!~      H(39,54)=S(39,54)-Bsoc(39,54)
!~      
!~      S(36,54)=1300.
!~      Bsoc(36,54)=-1000.
!~      H(36,54)=S(36,54)-Bsoc(36,54)
!~      
!~      S(36,57)=1200.
!~      Bsoc(36,57)=-1000.
!~      H(36,57)=S(36,57)-Bsoc(36,57)
    
!    S(71,71)=(S(71,70)+S(71,72)+S(70,71)+S(72,71))/4.
!    Bsoc(71,71)=(Bsoc(71,70)+Bsoc(71,72)+Bsoc(70,71)+Bsoc(72,71))/4.
!    where (S(:,:).GT.0)
!       H(:,:)=S(:,:)-BSOC(:,:)
!    elsewhere
!       H(:,:)=1.
!    endwhere
    S(:,:)=max(S(:,:),0.)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
    H(:,:)=max(H(:,:),1.)   ! pour avoir au moins 1 m


!!$
!!$    ! socle
!!$    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(:,:)
!!$
!!$
!!$    !     masque
!!$    filin='maskHUY40km.dat'
!!$    call lect_ieis(nx,ny,MK0,filin,DIRNAMEINP)
!!$    write(num_rep_42,*) 'fichier masque : ', filin 
!!$
!!$    !     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.).and.(H(I,J).gt.1.E-3)) 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=grid_topo,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)

    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=ghf_fich)

!    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
!    close(88)

!    call lect_input(1,'ghf',1,ghf,ghf_fich,file_ncdf)
    ! pour eviter plantage -O0
    call Read_Ncdf_var('ghf',ghf_fich,tab)
    ghf(:,:)  = tab(:,:)

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

!------------------------------------------------
! mko vrai partout (version 2006)
    MK0(:,:)=1

!------------------------------------------------ 


  end subroutine input_topo


end module lect_topo_anteis