source: trunk/LMDZ.COMMON/libf/evolution/xios_data.F90 @ 3989

MODULE xios_data

use netcdf, only: nf90_open, nf90_close, nf90_inquire_dimension, nf90_inq_dimid, nf90_noerr, nf90_nowrite, nf90_get_var, nf90_inq_varid

implicit none

character(19), parameter :: file1_daily  = "Xoutdaily4pem_Y1.nc"
character(19), parameter :: file2_daily  = "Xoutdaily4pem_Y2.nc"
character(20), parameter :: file1_yearly = "Xoutyearly4pem_Y1.nc"
character(20), parameter :: file2_yearly = "Xoutyearly4pem_Y2.nc"
character(256)           :: msg      ! For reading
integer                  :: fID, vID ! For reading

interface get_var
    module procedure get_var_1d, get_var_2d, get_var_3d, get_var_4d
end interface get_var

!=======================================================================
contains
!=======================================================================

!=======================================================================
SUBROUTINE load_xios_data(ngrid,nslope,nsoil_PCM,nsol,h2ofrost_PCM,co2frost_PCM,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
                         ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts,min_h2oice,min_co2ice)

use grid_conversion, only: lonlat2vect
use soil,            only: do_soil
use tendencies,      only: compute_tend
use metamorphism,    only: compute_frost

implicit none

! Arguments
!----------
integer,                       intent(in) :: ngrid, nslope, nsoil_PCM, nsol
real, dimension(ngrid,nslope), intent(in) :: h2ofrost_PCM, co2frost_PCM
real, dimension(ngrid),                       intent(out) :: ps_avg
real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1, watersurf_density_avg, d_h2oice, d_co2ice, min_h2oice, min_co2ice
real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg
real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts
real, dimension(ngrid,nsoil_PCM,nslope,nsol), intent(out) :: tsoil_ts, watersoil_density_ts

! Local variables
!----------------
integer                               :: islope, isoil, isol, nlon, nlat
real, dimension(:,:),     allocatable :: var_read_2d
real, dimension(:,:,:),   allocatable :: var_read_3d
real, dimension(:,:,:,:), allocatable :: var_read_4d
character(:),             allocatable :: num ! For reading slope variables
real, dimension(ngrid,nslope,2)       :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost

! Code
!-----
! Initialization
min_h2operice = 0.
min_co2perice = 0.
min_h2ofrost = 0.
min_co2frost = 0.
if (nslope == 1) then ! No slope
    allocate(character(0) :: num)
else ! Using slopes
    allocate(character(8) :: num)
endif

!~~~~~~~~~~~~~~~~~~~~~~~~ Year 1 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
! Open the NetCDF file of XIOS outputs
write(*,*) "Opening "//file1_yearly//"..."
call error_msg(nf90_open(file1_yearly,nf90_nowrite,fID),'open',file1_yearly)

! Get the dimensions
call error_msg(nf90_inq_dimid(fID,'lon',vID),'inq',file1_yearly)
call error_msg(nf90_inquire_dimension(fID,vID,len = nlon),'inq',file1_yearly)
call error_msg(nf90_inq_dimid(fID,'lat',vID),'inq',file1_yearly)
call error_msg(nf90_inquire_dimension(fID,vID,len = nlat),'inq',file1_yearly)

! Allocate and read the variables
allocate(var_read_2d(nlon,nlat),var_read_3d(nlon,nlat,nsoil_PCM))
do islope = 1,nslope
    if (nslope /= 1) then
        num='  '
        write(num,'(i2.2)') islope
        num = '_slope'//num
    endif
    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,1))
    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,1))
    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,1))
    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,1))
    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg_y1(:,islope))
enddo

! Close the NetCDF file of XIOS outputs
call error_msg(nf90_close(fID),"close",file1_yearly)
write(*,*) '> Data from "'//file1_yearly//'" downloaded!'

!~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
! Open the NetCDF file of XIOS outputs
write(*,*) "Opening "//file2_yearly//"..."
call error_msg(nf90_open(file2_yearly,nf90_nowrite,fID),'open',file2_yearly)

! Allocate and read the variables
call get_var('ps',var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,ps_avg)
do islope = 1,nslope
    if (nslope /= 1) then
        num='  '
        write(num,'(i2.2)') islope
        num = '_slope'//num
    endif
    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg(:,islope))
    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,2))
    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,2))
    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,2))
    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,2))
    if (do_soil) then
        call get_var('soiltemp'//num,var_read_3d)
        do isoil = 1,nsoil_PCM
            call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isoil),tsoil_avg(:,isoil,islope))
        enddo
        call get_var('waterdensity_surface'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,watersurf_density_avg(:,islope))
    endif
enddo
! Close the NetCDF file of XIOS outputs
call error_msg(nf90_close(fID),"close",file1_yearly)
write(*,*) '> Data from "'//file1_yearly//'" downloaded!'

!~~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Daily data ~~~~~~~~~~~~~~~~~~~~~~~~~
! Open the NetCDF file of XIOS outputs
write(*,*) "Opening "//file2_daily//"..."
call error_msg(nf90_open(file2_daily,nf90_nowrite,fID),'open',file2_daily)

! Allocate and read the variables
deallocate(var_read_2d,var_read_3d)
allocate(var_read_3d(nlon,nlat,nsol),var_read_4d(nlon,nlat,nsoil_PCM,nsol))
call get_var('ps',var_read_3d)
do isol = 1,nsoil_PCM
    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),ps_ts(:,isol))
enddo
call get_var('h2o_layer1',var_read_3d)
do isol = 1,nsoil_PCM
    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_h2o_ts(:,isol))
enddo
call get_var('co2_layer1',var_read_3d)
do isol = 1,nsoil_PCM
    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_co2_ts(:,isol))
enddo
if (do_soil) then
    do islope = 1,nslope
        if (nslope /= 1) then
            num='  '
            write(num,'(i2.2)') islope
            num = '_slope'//num
        endif
        call get_var('soiltemp'//num,var_read_4d)
        do isol = 1,nsol
            do isoil = 1,nsoil_PCM
                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),tsoil_ts(:,isoil,islope,isol))
            enddo
        enddo
        call get_var('waterdensity_soil'//num,var_read_4d)
        do isol = 1,nsol
            do isoil = 1,nsoil_PCM
                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),watersoil_density_ts(:,isoil,islope,isol))
            enddo
        enddo
    enddo
endif
deallocate(var_read_3d,var_read_4d,num)

! Close the NetCDF file of XIOS outputs
call error_msg(nf90_close(fID),"close",file1_daily)
write(*,*) '> Data from "'//file2_daily//'" downloaded!'

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Compute frost from yearly minima
call compute_frost(ngrid,nslope,h2ofrost_PCM,min_h2ofrost,co2frost_PCM,min_co2frost)

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Compute ice tendencies from yearly minima
write(*,*) '> Computing surface ice tendencies'
call compute_tend(ngrid,nslope,min_h2operice + min_h2ofrost,d_h2oice)
write(*,*) 'H2O ice tendencies (min/max):', minval(d_h2oice), maxval(d_h2oice)
call compute_tend(ngrid,nslope,min_co2perice + min_co2frost,d_co2ice)
write(*,*) 'CO2 ice tendencies (min/max):', minval(d_co2ice), maxval(d_co2ice)

END SUBROUTINE load_xios_data
!=======================================================================

!=======================================================================
SUBROUTINE get_timelen(filename,timelen)

use netcdf

implicit none

! Arguments
! ---------
character(*), intent(in)  :: filename
integer,      intent(out) :: timelen

! Local variables
! ---------------
integer :: ncid  ! File ID
integer :: dimid ! Dimension ID
integer :: ierr  ! Return codes

! Code
! ----
! Open the NetCDF file
ierr = nf90_open(trim(filename),NF90_NOWRITE,ncid)
if (ierr /= nf90_noerr) then
    write(*,*) "Error opening file:", trim(nf90_strerror(ierr))
    error stop
endif

! Get the dimension ID for 'time_counter'
ierr = nf90_inq_dimid(ncid,"time_counter",dimid)
if (ierr /= nf90_noerr) then
    write(*,*) "Error getting dimid 'time_counter':", trim(nf90_strerror(ierr))
    error stop
endif

! Get the size of the dimension 'time_counter'
ierr = nf90_inquire_dimension(ncid,dimid,len = timelen)
if (ierr /= nf90_noerr) then
    write(*,*) "Error getting dimension length:", trim(nf90_strerror(ierr))
    error stop
endif

! Close the file
ierr = nf90_close(ncid)
if (ierr /= nf90_noerr) then
    write(*,*) "Error closing file:", trim(nf90_strerror(ierr))
    error stop
endif

END SUBROUTINE get_timelen
!=======================================================================

!=======================================================================
SUBROUTINE error_msg(ierr,typ,nam) 

implicit none

integer,      intent(in) :: ierr !--- NetCDF error code
character(*), intent(in) :: typ  !--- Type of operation
character(*), intent(in) :: nam  !--- Field/File name

if (ierr == nf90_noerr) return
select case(typ)
    case('inq')  ; msg = "Dim/Field <"//trim(nam)//"> is missing"
    case('get')  ; msg = "Reading failed for <"//trim(nam)//">"
    case('put')  ; msg = "Writing failed for <"//trim(nam)//">"
    case('open') ; msg = "File opening failed for <"//trim(nam)//">"
    case('close'); msg = "File closing failed for <"//trim(nam)//">"
    case default
        write(*,*) 'There is no message for this error.'
        error stop
end select
call abort_gcm(__FILE__,trim(msg),ierr)

END SUBROUTINE error_msg
!=======================================================================

!=======================================================================
SUBROUTINE get_var_1d(var,v)

implicit none

character(*), intent(in) :: var
real, dimension(:), intent(out) :: v

call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
call error_msg(nf90_get_var(fID,vID,v),"get",var)

END SUBROUTINE get_var_1d
!=======================================================================

!=======================================================================
SUBROUTINE get_var_2d(var,v)

implicit none

character(*), intent(in) :: var
real, dimension(:,:), intent(out) :: v

call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
call error_msg(nf90_get_var(fID,vID,v),"get",var)

END SUBROUTINE get_var_2d
!=======================================================================

!=======================================================================
SUBROUTINE get_var_3d(var,v)

implicit none

character(*), intent(in) :: var
real, dimension(:,:,:), intent(out) :: v

call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
call error_msg(nf90_get_var(fID,vID,v),"get",var)

END SUBROUTINE get_var_3d
!=======================================================================

!=======================================================================
SUBROUTINE get_var_4d(var,v)

implicit none

character(*), intent(in) :: var
real, dimension(:,:,:,:), intent(out) :: v

call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
call error_msg(nf90_get_var(fID,vID,v),"get",var)

END SUBROUTINE get_var_4d
!=======================================================================

END MODULE xios_data