PrepareNetcdfOutput.f90

Go to the documentation of this file.

!==============================================================================
!  contains subroutines  PrepareNetcdfOutput, InitNCVar, m_check 
!==============================================================================
module netcdftmmodule
    
contains
 
! ==============================================================================
 
!==============================================================================
SUBROUTINE preparenetcdfoutput(nspc)
    use filelistmodule, only: nc_file_name,nc_file_id, biodiv_nc,biomass_nc, &
                              number_nc, seed_nc, ingrowth_nc, npp_nc, basalarea_nc,&
                              antagonist_nc, pollen_nc, lai_nc, light_nc, heightstruct_nc 
 
    use netcdf, only: nf90_char, nf90_global, nf90_int, nf90_netcdf4, nf90_real,&
                      nf90_unlimited, nf90_def_dim, nf90_put_att, nf90_put_var, &
                      nf90_create, nf90_def_var, nf90_enddef
    use all_par, only: unitofspatialdata, lonrealstart, latrealstart,&
                       cellsidelength, simustartyear, maxlat, maxlon, &
                       numyrs, writeoutput, maxhc, maxlc, reportintervals, &
                       nreportintervals, spec, experimentid
    IMPLICIT NONE
    
    INTEGER, INTENT(in) :: nspc
    
    INTEGER, dimension(:) :: dimYeSpXY(4), dimLimits(3)
    INTEGER :: dimId_year, dimId_speciesLength,dimId_species, dimId_x, dimId_y,  dimId_hc, dimId_lc
    INTEGER :: varId_year,varId_species, varId_x, varId_y, varId_crs, varId_hc, varId_lc
    
    
    REAL :: lats(maxlat)
    REAL :: lons(maxlon)
    INTEGER :: years(numyrs)
    INTEGER :: species(nspc+1)
    INTEGER :: hc(maxhc+1)
    INTEGER :: lc(maxlc)
    
    INTEGER :: i, numReportYears, actualYear, rpti 
 
    character*6:: yearChar
 
    ! variables for date/time
    INTEGER date_time(8)
    character*5:: zone
    character*25 :: dateTimeString
 
 
    ! if netcdf output is not set, we do nothing
    if(.not. writeoutput%netcdf) return
 
    ! create netcdf file; the NC_File_ID is provided by the nf90_create function
    call m_check(nf90_create(nc_file_name, nf90_netcdf4, nc_file_id), __line__, __file__)  
  
    ! get date/time and store them to a string in a format Netcdf likes
    call date_and_time(zone=zone, values= date_time)
    write(datetimestring, '(I4.4,"-",2(I2.2,"-"),"T",2(I2.2,":"),I2.2,A5)')  date_time(1),date_time(2),&
        date_time(3), date_time(5), date_time(6), date_time(7), zone
    ! define global attributes for file
    call m_check(nf90_put_att(nc_file_id, nf90_global , "title", "TreeMig model output"), __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "institution", &
            "Swiss Federal Institute for Forest, Snow and Landscape Research WSL"), __line__,__file__)
 
    call m_check(nf90_put_att(nc_file_id, nf90_global , "institute_id", "WSL"),          __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "experiment_id", experimentid),  __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "model_id", "TreeMig"),          __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "contact", "maintainer_treemig@wsl.ch"),&
                __line__,__file__)
 
    call m_check(nf90_put_att(nc_file_id, nf90_global , "reference", "TreeMig model described by Lischke et al., "&
        "2006 (https://www.dora.lib4ri.ch/wsl/islandora/object/wsl%3A6049)"),            __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "treeMig_version", "0.0.1"),     __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "conventions", "CF-1.10"),       __line__,__file__)
    call m_check(nf90_put_att(nc_file_id, nf90_global , "creation_date", datetimestring),&
                __line__,__file__)
 
    ! define the dimensions. NetCDF will hand back an ID for each. 
    
    ! setup time dimension
    call m_check(nf90_def_dim(nc_file_id, "time",    nf90_unlimited,  dimid_year), __line__,__file__)
 
    ! setup species dimension
    call m_check(nf90_def_dim(nc_file_id, "species", nspc+1,    dimid_species),    __line__,__file__)
 
    ! setup x,y coordinate dimension
    call m_check(nf90_def_dim(nc_file_id, "x",       maxlon,  dimid_x),            __line__,__file__)
    call m_check(nf90_def_dim(nc_file_id, "y",       maxlat,  dimid_y),            __line__,__file__)
    call m_check(nf90_def_dim(nc_file_id, "string80", 80,    dimid_specieslength), __line__,__file__) ! ? maybe number of species
 
    ! define time variable and its attributes
    call m_check(nf90_def_var(nc_file_id, "time", nf90_int, dimid_year,varid_year), __line__,__file__)
 
! axis: T unit year standard_name: time    
        write(yearchar,"(I6)")  simustartyear
        call m_check(nf90_put_att(nc_file_id, varid_year, "axis", "T"),             __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_year, "long_name", "time"),     __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_year, "standard_name", "time"), __line__,__file__)
 
        call m_check(nf90_put_att(nc_file_id, varid_year, "units", "years since "//trim(yearchar)//"-1-1 0:0:0"),&
                __line__,__file__)
 
! define species variable  
    call m_check(nf90_def_var(nc_file_id, "species", nf90_char, &
            (/ dimid_specieslength, dimid_species /), varid_species),&
            __line__,__file__)
! define x variable  and its attributes
    call m_check(nf90_def_var(nc_file_id, "x",     nf90_real, dimid_x,  varid_x), __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_x, "standard_name", "projection_x_coordinate"),&
                __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_x, "long_name", "x coordinate of projection"),&
                __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_x, "units", "m"),__line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_x, "axis", "X"),__line__,__file__)
! define y variable  and its attributes
    call m_check(nf90_def_var(nc_file_id, "y",     nf90_real, dimid_y,      varid_y),__line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_y, "standard_name", "projection_y_coordinate"),&
                __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_y, "long_name", "y coordinate of projection"),&
                __line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_y, "units", "m"),__line__,__file__)
        call m_check(nf90_put_att(nc_file_id, varid_y, "axis", "Y"),__line__,__file__)
! define coordinates
    call m_check(nf90_def_var(nc_file_id, "crs",     nf90_int, [INTEGER::],   varId_crs),__LINE__,__FILE__)
! from where?
    
  
    dimyespxy(1:4) =  (/ dimid_year,dimid_species, dimid_x, dimid_y /)   
! here finally the different output variables are defined, allocated and initialized
    if(writeoutput%biodiv)&
        call initncvar(biodiv_nc, 2, (/ dimid_year,dimid_x, dimid_y/), (/maxlon,maxlat/))
 
    
    dimlimits(1:3) = (/nspc+1,maxlon, maxlat/)
    if(writeoutput%biomass)&
        call initncvar(biomass_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%number)&
        call initncvar(number_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%seeds)&
        call initncvar(seed_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%ingrowth)&
        call initncvar(ingrowth_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%NPP)&
        call initncvar(npp_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%basalArea)&
        call initncvar(basalarea_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%antagonists)&
        call initncvar(antagonist_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%pollen)&        
        call initncvar(pollen_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%lai)&        
        call initncvar(lai_nc, 3, dimyespxy, dimlimits)
    if(writeoutput%hstruct) then
        do i=1,size(heightstruct_nc)
            call  initncvar(heightstruct_nc(i), 3, dimyespxy, dimlimits)
        end do
    end if
! initializing the light output
    if(writeoutput%light) then        
        call m_check(nf90_def_dim(nc_file_id, "hc",       maxhc+1,  dimid_hc),__line__,__file__)
        call m_check(nf90_def_dim(nc_file_id, "lc",       maxlc,  dimid_lc),  __line__,__file__)
 
        call m_check(nf90_def_var(nc_file_id, "hc",    nf90_int,  dimid_hc,    varid_hc),&
         __line__,__file__)
            call m_check(nf90_put_att(nc_file_id, varid_hc, "long_name", "height class"),&
            __line__,__file__)
 
 
        call m_check(nf90_def_var(nc_file_id, "lc",    nf90_int,  dimid_lc,    varid_lc),&
         __line__,__file__)
            call m_check(nf90_put_att(nc_file_id, varid_lc, "long_name", "light class"),&
            __line__,__file__)
 
          
        call initncvar(light_nc, 4, (/ dimid_year,dimid_hc, dimid_lc, dimid_x, dimid_y /),&
                                    (/maxhc+1,maxlc,maxlon, maxlat/))
        do i=0,(maxhc)
            hc(i+1) = i
        end do
        do i=1,maxlc
            lc(i) = i
        end do
    end if  
! the values for the x and y variable are collected
    do i=0,(maxlon-1)
        lons(i+1) = (lonrealstart+i*cellsidelength/unitofspatialdata)
    end do
    do i=0,(maxlat-1)
        lats(i+1) = (latrealstart+i*cellsidelength/unitofspatialdata)
    end do    
! the years for the year variable are collected 
    numreportyears = 0
    do actualyear=simustartyear, simustartyear+numyrs-1
        rpti = 1
        do i=1,nreportintervals  ! determine current report interval
            if (actualyear >= reportintervals(i,2) .and.  actualyear <= reportintervals(i,3) )  then
                rpti=reportintervals(i,1)
                exit
            end if 
        end do
        if(mod(actualyear-simustartyear, rpti) == 0)then
            years(numreportyears+1) = actualyear
            numreportyears = numreportyears+1
        end if
    end do
 
    do i=0,(nspc)
        species(i+1) = (i+1)
    end do
    ! ready with the definitions
    call m_check( nf90_enddef(nc_file_id),__line__,__file__)
    ! now the x, y and year values are added to the corresponding variables
    call m_check( nf90_put_var(nc_file_id, varid_y, lats), __line__,__file__)
 
    call m_check( nf90_put_var(nc_file_id, varid_x, lons), __line__,__file__)
 
    call m_check( nf90_put_var(nc_file_id, varid_year, years(1:numreportyears)), __line__,__file__)
 
    call m_check( nf90_put_var(nc_file_id, varid_species, (/(spec(i)%nl, i=1, nspc),"all_species"/)),&
            __line__,__file__)
 
    if(writeoutput%light)then        
        call m_check( nf90_put_var(nc_file_id, varid_hc, hc), __line__,__file__)
 
        call m_check( nf90_put_var(nc_file_id, varid_lc, lc), __line__,__file__)
 
    end if 
 
end SUBROUTINE preparenetcdfoutput
 
!==============================================================================
! @brief InitNCVar  
! ***************************************** 
!==============================================================================
SUBROUTINE initncvar(var, numDims, NCDims, limits)
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: ieee_arithmetic
    use filelistmodule
    use netcdf
    use loggermodule
    IMPLICIT NONE
    
    TYPE(ncvar), INTENT(inout) :: var
    INTEGER, dimension(:), INTENT(in) :: NCDims
    INTEGER, INTENT(in) :: numDims
    INTEGER, dimension(:), INTENT(in) :: limits
    
    REAL :: fillval
    INTEGER :: alloc_st
    character (len=500) :: errormsg
    fillval = ieee_value(fillval, ieee_quiet_nan) ! ?????
    ! defines the variable and its attributes 
    call m_check(nf90_def_var(nc_file_id, var%name, nf90_real, ncdims, var%varId, deflate_level = 4),&
                __line__,__file__)
    ! sets a fillval ????
    call m_check(nf90_put_att(nc_file_id,  var%varId, "_FillValue", fillval), __line__,__file__)
 
    ! sets a unit
    call m_check(nf90_put_att(nc_file_id,  var%varId, "units", var%attr_unit), __line__,__file__)
 
 
    ! depending on the number of dimensions, the arrays for this variable is allocated and filled with fillval
    select case (numdims)
        case (4)
            allocate (var%values4D(1:limits(1), 1:limits(2), 1:limits(3), 1:limits(4)), stat=alloc_st, errmsg=errormsg)
            var%values4D = fillval
        case (3)
            allocate (var%values3D(1:limits(1), 1:limits(2), 1:limits(3)), stat=alloc_st, errmsg=errormsg)
            var%values3D = fillval
        case (2)
            allocate (var%values2D(1:limits(1), 1:limits(2)), stat=alloc_st, errmsg=errormsg)
            var%values2D = fillval
        case default
            call logerror("Invalid number of NETCDF dimensions")
            
            error stop
    end select
    call checkalloc(alloc_st, errormsg, var%name//" Output")
end SUBROUTINE initncvar
end module netcdftmmodule
 
!==============================================================================
!==============================================================================
     SUBROUTINE m_check(status, lineNumber, fileName)
         use loggermodule
         use netcdf
         INTEGER, INTENT (in) :: status
         INTEGER, INTENT (in) :: lineNumber
         character(len=*), INTENT (in) :: fileName
         if(status /= nf90_noerr) then 
             write(logmessage,"(A,I3.3)") "Error in "//filename//" on line ", linenumber
            call logerror(logmessage)
             call logerror(trim(nf90_strerror(status)))
             error stop "Execution has been stopped due to a netcdf error"
         end if
    end SUBROUTINE m_check
!==============================================================================