source: trunk/SOURCES/GrIce2sea_files/lect_GrIce2sea_gen_nc.f90 @ 4

!> \file lect_GrIce2sa_gen_nc.f90
!!Module pour la lecture de la topographie
!<

!> \namespace lect_topo_green_gen
!! Module pour la lecture de la topographie
!! \author Cat
!! \date 31 oct 2011
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<

module lect_topo_green_gen

  use module3D_phy
  use param_phy_mod
  use interface_input
  use io_netcdf_grisli

  character(len=100) :: topo_surf      !< surface file name
  character(len=100) :: correc_surf    !< ice-real correction file name
! character(len=100) :: topo_thick    !< thickness file name
  character(len=100) :: topo_bed       !< bedrock file name
  character(len=100) :: mask_grounded  !< mask file name
  character(len=100) :: longitude      !< longitude file name
  character(len=100) :: latitude       !< latitude file name
  character(len=100) :: heatflux       !< geothermal heat flux
  logical            :: corr_firn      !< eventuellement correction firn (pour l'instant surface dejà corrigee)


  real,dimension(nx,ny) :: work_tab 

  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat

contains

  subroutine input_topo

    namelist/topo_groen_gen/topo_surf,correc_surf,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_groen_gen)

    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428)'!  module  lect_topo_ant_gen                                '
    write(num_rep_42,topo_groen_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'

! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd

    call lect_input(1,'Bsoc',1,Bsoc,topo_bed,file_ncdf)    ! socle
    Bsoc0(:,:) = Bsoc(:,:)

    call lect_input(1,'S',1,S,topo_surf,file_ncdf)         ! surface
    
!    call lect_input(1,'H',1,H,topo_thick,file_ncdf)       ! epaisseur calculee ensuite

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

    if (corr_firn) then                                    ! density correction     
       call lect_input(1,'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

    if (itracebug.eq.1)  call tracebug(' apres corr firn')
    H(:,:)=max(0.,H(:,:))
    H0(:,:)= H(:,:)
    S0(:,:)= S(:,:)


    mk0(:,:)=0.
    mk_init(:,:)=0.

    call lect_input(1,'work_tab2',2,work_tab,mask_grounded,file_ncdf)  ! masque


! definition du Mk_init dans ice2sea
!__________________________________________
!  outcrops         ->  eventuellement lu dans un autre fichier
!  ocean libre      ->  0
!  toundra          ->  1
!  pose avec glace  ->  2
!  epaisseur fixe   ->  3  masque pour ne pas traiter ces points
!  ice shelves      ->  4

    mk_init(:,:) = nint(work_tab(:,:))

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


! make consistency of input data files (flotation, thickness ...)
! Rappel : dans param_phy_mod
! coef_Sflot = (Row-Ro)/Row    S = coef_Sflot * H + sealevel pour les shelves
! coef_Bflot = -Ro/Row   B = coef_Bflot * H + sealevel pour les shelves


    do j=1,ny
       do i=1,nx
          if (mk_init(i,j).eq.0) then                     ! ocean
             H(i,j)    = 1.                               ! valeur bidon pour l'ocean
             S(i,j)    = coef_Sflot * H(i,j)
             B(i,j)    = coef_Bflot * H(i,j)
             Bsoc(i,j) = -100.                            ! on surcreuse le socle

             flot(i,j) = .true.

          else if  (mk_init(i,j).eq.1) then               ! toundra
             S(i,j)    = max(Bsoc(i,j),0.)
             H(i,j)    = S(i,j) - Bsoc(i,j) 
             B(i,j)    = Bsoc(i,j)

             flot(i,j) = .false.

          else if  (mk_init(i,j).eq.2) then               ! glace
             H(i,j)    = S(i,j) - Bsoc(i,j) 
             B(i,j)    = Bsoc(i,j)
             flot(i,j) = .false.   

          else if  (mk_init(i,j).eq.3) then               ! zones a ne pas toucher
             ! le traitement "fixe" se fait dans prescribe_mod 
             H(i,j)    = S(i,j) - Bsoc(i,j) 
             B(i,j)    = Bsoc(i,j)
             flot(i,j) = .false.   

          else if (mk_init(i,j).eq.4) then                ! ice shelf
             H(i,j)    = S(i,j) / coef_Sflot
             B(i,j)    = S(i,j) - H(i,j)
             Bsoc(i,j) = min (B(i,j)-20.,Bsoc(i,j))       ! on surcreuse le socle si besoin
             flot(i,j)=.true.
          end if

       end do
    end do

! garder les cartes debut de run
    Bsoc0(:,:) = Bsoc(:,:)
    H0(:,:)= H(:,:)
    S(:,:) = B(:,:) +H(:,:)
    S0(:,:)= S(:,:)


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


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


    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)        ! 

    if (itracebug.eq.1)  call tracebug(' apres lect lon/lat')


! pour recuperer le xmin xmax ymin ymax, on s'inspire de ce que Cat a fait pour l'antarc :
!   write(6,*) 'topo_surf',topo_surf
!   call Read_Ncdf_var('x',trim(topo_surf),Tab1D)    ! lit le tableau x
!   xmin = minval(Tab1D(:))
!   xmax = maxval(Tab1d(:))
!   call Read_Ncdf_var('y',trim(topo_surf),Tab1D)    ! lit le tableau y
!   ymin = minval(Tab1D(:))
!   ymax = maxval(Tab1d(:))

    open(22,file=trim(dirnameinp)//'long_ZBL_CISM_5km.dat')                         ! 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)


    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel
    !     pas de version correcte pour l'instant -> valeurs initiales
    S_sealev(:,:) = S0(:,:)
    H_sealev(:,:) = H0(:,:)
    B_sealev(:,:) = Bsoc(:,:)
    M_sealev(:,:) = Mk0(:,:)

    if (itracebug.eq.1)  call tracebug(' apres lect entete')
  
    call lect_input(1,'ghf',1,ghf,heatflux,file_ncdf)               ! flux geothermique


    ! 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(:,:),70.)           ! pour sea rise

    ! pour passer les flux des W/m2 au J/m2/an      
    ghf(:,:)=-SECYEAR*ghf(:,:)     !  les flux sont donnes en W/m2, correction pour passer a l'annee
    ghf0(:,:)=ghf(:,:)

    if (itracebug.eq.1)  call tracebug(' a la fin d input_topo')

  end subroutine input_topo


end module lect_topo_green_gen