source: trunk/INPUT/MISMIP_3D/mismip3.f90 @ 80

program make_grisli_input

!-------------------------------------------------------!
!      Programm writen by V. Peyaud, nov 2011           !
!                                                       !
!      This is the squeleton of the programs            !
!        that generates NetCDF input file for Grisli    !
!                                                       !
!-------------------------------------------------------!

use netcdf
use  io_netcdf_grisli
!use sorties_ncdf_grisli
implicit none


! ---------------------------------------
!      5 km   10 km
! nx   321    161
! ny    21     11
! dx   5000.  10000.
integer,parameter :: nx = 321      ! length of the 1st dim
integer,parameter :: ny = 23       ! length of the 2nd dim
real,parameter    :: dx  =5000.    ! resolution (m, calcule a partir du fichier .nc)
real,parameter    :: dy  =5000.    ! resolution (m)
integer           :: i,j           ! indexes
! ---------------------------------------
character(len=30) :: fil_sortie  ! output file
integer                        :: status         ! utilise dans les appels netcdf
integer                        :: ncid
character(len=3), dimension(2) :: dimname        ! le nom de dimension dont on teste l'existence
real*8,dimension(:),  pointer  :: Tab1d          ! pointeur vers les variables dimensions (et 1D)
real*8,dimension(:,:),pointer  :: Tab            ! pointeur vers les variables 2D
! ---------------------------------------
character(len=20),dimension(3) :: dimnames2dxymaj       !< dimensions pour netcdf pour tableau 2d pour les noeud majeur
real,dimension(nx,ny) :: S, H, B, Bsoc,Ubar,Usurf,Uslid,betamx  ! ,allocatable
real,dimension(nx) :: xx
real,dimension(ny) :: yy
! ---------------------------------------
character(len=50) :: inputfile ! input file
real,dimension(nx) :: SW, HW, Bw, Bwsoc,Ubarw,Usurfw,Uslidw,betamxw


! coordinates en m

do i=1,nx
     xx(i) = -(nx-1)/2.*dx + (i-1)*dx
enddo
do j=1,ny
     yy(j) = -(ny-1)/2.*dy + (j-1)*dy
enddo 

!Use the standard output from Winnie

! m3_1Stnd05_class01_rec01.nc
!inputfile='profil_mismip3d_stnd5.nc' ! 5 km winnie standard run

!inputfile='profil_Winnie_StndA5_steady_float_1D.nc'        ! 5 km winnie standard run
inputfile='profil_Winnie_StndA5_steady_newbeta_1D.nc'       ! 5 km avec les valeurs de glissement et deformation de 
                                                            ! mismip

allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('S',trim(inputfile),Tab1d)

write(6,*) trim(inputfile),Tab1d(161)
    SW=Tab1d
deallocate(Tab1d)

allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('H',trim(inputfile),Tab1d)
    HW=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('B',trim(inputfile),Tab1d)
    BW=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('Bsoc',trim(inputfile),Tab1d)
    BWsoc=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('Ubar',trim(inputfile),Tab1d)
    Ubarw=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('Usurf',trim(inputfile),Tab1d)
    Usurfw=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('Uslid',trim(inputfile),Tab1d)
    Uslidw=Tab1d
deallocate(Tab1d)
allocate(Tab1d(nx))
    call Read_Ncdf_var1d_Real('Betamx',trim(inputfile),Tab1d)
    betamxw=Tab1d
deallocate(Tab1d)


do i=1,nx
    S(i,:)      = SW(i)
    H(i,:)      = HW(i)
    B(i,:)      = BW(i)
    Bsoc(i,:)   = BWsoc(i)
    Ubar(i,:)   = Ubarw(i)
    Usurf(i,:)  = Usurfw(i)
    Uslid(i,:)  = Uslidw(i)
    betamx(i,:) = betamxw(i)
 enddo
write(6,*) 'apres lecture sw(161)=',Sw(161),'     S(161,12)=',S(161,12)

!    Bsoc(i,:)=-800+(i-1)*dx
!enddo
!Bsoc(:,:) = 0.
!B(:,:)    = Bsoc(:,:)
!H(:,:)    = 100.
!S(:,:)    = B(:,:)+H(:,:)

!
! creation et fermeture du fichier output

write(*,*) ncid
fil_sortie='topo2grisli.nc'

status  = nf90_create(trim(fil_sortie),nf90_write,ncid)
if (status/=nf90_noerr) then    
       write(*,*)"unable to create netcdf file : ",trim(fil_sortie),ncid
endif

status  = nf90_close(ncid)
status  = nf90_open(trim(fil_sortie),nf90_write,ncid)
if (status/=nf90_noerr) then    
       write(*,*)"unable to close netcdf file : ",ncid
endif

dimname(1)='x'
dimname(2)='y'

write(6,*) "writing topo in NC file"
write(6,*) "n ",nx,ny


! ecriture dans le fichier output

call write_ncdf_dim(dimname(1),trim(fil_sortie),nx)
call write_ncdf_dim(dimname(2),trim(fil_sortie),ny)

! Allouer & Deallouer chaque fois pour eviter des conflits d'identifiants

allocate(Tab1d(nx))
    Tab1d(:)=xx(:)
    call write_ncdf_var(dimname(1),dimname(1),trim(fil_sortie), Tab1d ,'float')
deallocate(Tab1d)
allocate(Tab1d(ny))
    Tab1d(:)=yy(:)    
    call write_ncdf_var(dimname(2),dimname(2),trim(fil_sortie), Tab1d ,'float')

!call Copy_Ncdf_att_var(dimname(1),trim(filin),trim(fil_sortie))

allocate(Tab(nx,ny))
    Tab=S
    call write_ncdf_var('S',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=H
    call write_ncdf_var('H',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=B
    call write_ncdf_var('B',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=Bsoc
    call write_ncdf_var('Bsoc',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=Ubar
    call write_ncdf_var('Ubar',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=Uslid
    call write_ncdf_var('Uslid',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=Usurf
    call write_ncdf_var('Usurf',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)
allocate(Tab(nx,ny))
    Tab=betamx
    call write_ncdf_var('betamx',dimname(:),trim(fil_sortie), Tab ,'float')
deallocate(Tab)


!call write_ncdf_var(dimname(1),dimname(1),trim(fil_sortie),Tab1d,'float')
!call Copy_Ncdf_att_var(dimname(1),trim(filin),trim(fil_sortie))


!call write_ncdf_dim('x',trim(fil_sortie(j)),nx)
!call write_ncdf_dim('y',trim(fil_sortie(j)),ny)

!call write_ncdf_var('S',dimnames2dxymaj,trim(fil_sortie),   S  ,'double')
!call write_ncdf_var('H',dimnames2dxymaj,trim(fil_sortie(k)),   H   ,nrecs(k),idef(k,p),'double')
!call write_ncdf_var('Bsoc',dimnames2dxymaj,trim(fil_sortie(k)),Bsoc,nrecs(k),idef(k,p),'double')
!call write_ncdf_var('B',dimnames2dxymaj,trim(fil_sortie(k)),  ,B   ,nrecs(k),idef(k,p),'double')

!            read(72,*)
!            read(72,*)
!            read(72,*)
!            read(72,*)
!            read(72,*) long_name
!            read(72,*) standard_name
!            read(72,*) unit
!            read(72,*) descriptions
!            read(72,*)


end program make_grisli_input