source: trunk/SOURCES/Old-sources/lect-Ant_CISM_gen_dat.f90 @ 334

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

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

module lect_topo_ant_CISM_gen

use module3D_phy
use param_phy_mod
use interface_input
character(len=100) :: topo_surf      !< surface
character(len=100) :: correc_surf    !< ice-real correction
character(len=100) :: topo_thick     !< thickness
character(len=100) :: topo_bed       !< bedrock
character(len=100) :: mask_grounded  !< mask
character(len=100) :: longitude      !< longitude
character(len=100) :: latitude       !< latitude
character(len=100) :: heatflux       !< geothermal heat flux
logical            :: corr_firn      !< eventuellement correction firn (deja fait dans Le Brocq)
real,dimension(nx,ny) :: work_tab          
Character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

contains

  subroutine input_topo

    namelist/topo_CISM_gen/topo_surf,correc_surf,topo_thick,topo_bed,mask_grounded,longitude,latitude,heatflux
    if (itracebug.eq.1)  call tracebug(' Entree dans routine input_topo')

428 format(A)

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

    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428)'!  module  lect_topo_ant_CISM                               '
    write(num_rep_42,topo_CISM_gen)
    write(num_rep_42,428)'!___________________________________________________________' 


    if ((trim(correc_surf).eq.'no').or.(trim(correc_surf).eq.'NO')) then 
       corr_firn=.false.
    else
       corr_firn=.true.
    end if

    topo_surf     = trim(dirnameinp)//trim(topo_surf)
    correc_surf   = trim(dirnameinp)//trim(correc_surf)
    topo_thick    = trim(dirnameinp)//trim(topo_thick)
    topo_bed      = trim(dirnameinp)//trim(topo_bed)
    mask_grounded = trim(dirnameinp)//trim(mask_grounded)
    longitude     = trim(dirnameinp)//trim(longitude)
    latitude      = trim(dirnameinp)//trim(latitude)
    heatflux      = trim(dirnameinp)//trim(heatflux)

    ! Transformation des fichiers .dat en .nc

    file_ncdf= trim(dirnameinp)//trim(runname)//'.nc'

    call lect_input(3,'Bsoc',1,Bsoc,topo_bed,file_ncdf)
    ! lecture adaptee aux fichiers ZBL.dat
    !call lect_datfile(nx,ny,Bsoc,1,topo_bed)               ! bedrock

    Bsoc0(:,:) = Bsoc(:,:)

    call lect_input(3,'S',1,S,topo_surf,file_ncdf)
    !call lect_datfile(nx,ny,S,1,topo_surf)                 ! surface 


    call lect_input(3,'H',1,H,topo_thick,file_ncdf)
    !call lect_datfile(nx,ny,H,1,topo_thick)                ! thickness

if (itracebug.eq.1)  call tracebug(' avant corr firn')

    if (corr_firn) then                                    ! density correction
       !call lect_datfile(nx,ny,work_tab,1,correc_surf)    ! for CISM files but not 
       ! Cat treatment of Le Brocq
       ! already done
if (itracebug.eq.1)  call tracebug(' apres corr firn')
       call lect_input(3,'work_tab1',1,work_tab,correc_surf,file_ncdf)

       where (H(:,:).gt.0.)
          work_tab(:,:)=min(H(:,:),work_tab(:,:))   
          S(:,:) = S(:,:) - work_tab(:,:)                  ! corrected surface
          H(:,:) = H(:,:) - work_tab(:,:)                  ! corrected thickness
       end where
    end if
    H(:,:)=max(0.,H(:,:))
    H0(:,:)= H(:,:)
    S0(:,:)= S(:,:)

!!$call lect_datfile(nx,ny,work_tab,1,mask_grounded)   ! masque
!!$mk0(:,:) = int(work_tab)
    mk0(:,:)=0.

    call lect_input(3,'work_tab2',2,work_tab,mask_grounded,file_ncdf)
    !call lect_datfile(nx,ny,work_tab,2,mask_grounded)    ! masque

    ! On lit la colonne 2 : 0 pour la glace 1 ailleurs y compris les iles qu'on veut enlever
    ! ij_grisli2mask.sh points-iles-a-enlever.dat

    ! Le masque iceberg est 0 quand glace et 1 quand mer (evite des incoherences sur les iles)
    ! ce masque est fait en bricolant le *berg.nc
if (itracebug.eq.1)  call tracebug(' avant worktab')


    mk0(:,:)=work_tab(:,:)

! bloc pour enlever les iles : merde sur les outcrop
where ((mk0(:,:).eq.1).and.(H(:,:).gt.1.))
   H(:,:)=0.
   S(:,:)=0.
   Bsoc(:,:)=min(Bsoc(:,:),-100.)
   B(:,:)=0.
end where
   Bsoc0(:,:) = Bsoc(:,:)
   H0(:,:)= H(:,:)
   S0(:,:)= S(:,:)

  if (itracebug.eq.1)  call tracebug(' apres worktab')


    write(266,*) 'petite ile qui merde'
    do j=227,224,-1
       write(266,*) (H(i,j),i=43,46)
    end do
    write(266,*) 'socle'
    do j=227,224,-1
       write(266,*) (Bsoc(i,j),i=43,46)
    end do


! pour l'Antarctique masque mko vrai partout (version 2006)
 MK0(:,:)=1  ! mis apres la lecture cptr


    ! 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

    ! le calcul des courbures du socle devrait etre dans le dragging

    call lect_input(1,'Xlong',1,Xlong,longitude,file_ncdf)
    !call lect_datfile(nx,ny,Xlong,1,longitude)      ! lecture des coordonnées geographiques
    
    call lect_input(1,'Ylat',1,Ylat,latitude,file_ncdf)
    !call lect_datfile(nx,ny,Ylat,1,latitude)        !

    open(22,file=longitude)                         ! juste pour avoir l'en tete
    read(22,*) i,i,xmin,xmax,ymin,ymax
    xmin=xmin/1000.
    ymin=ymin/1000.
    xmax=xmax/1000.
    ymax=ymax/1000.
    close(22)


    call lect_input(1,'ghf',1,ghf,heatflux,file_ncdf)
    !call lect_datfile(nx,ny,ghf,1,heatflux)
 
    ! for the purposes of SeaRISE, GTHF should be capped at a value of -0.07 w/m^2,
    ! and the field provided in the netCDF file does not include these updated values

    ghf(:,:)=min(ghf(:,:),0.07)

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

  end subroutine input_topo


end module lect_topo_ant_CISM_gen