diff --git a/sorc/CMakeLists.txt b/sorc/CMakeLists.txt index 92838fa68..0bc698c8d 100644 --- a/sorc/CMakeLists.txt +++ b/sorc/CMakeLists.txt @@ -18,3 +18,4 @@ add_subdirectory(chgres_cube.fd) add_subdirectory(orog_mask_tools.fd) add_subdirectory(sfc_climo_gen.fd) add_subdirectory(vcoord_gen.fd) +add_subdirectory(fvcom_tools.fd) diff --git a/sorc/fvcom_tools.fd/CMakeLists.txt b/sorc/fvcom_tools.fd/CMakeLists.txt new file mode 100644 index 000000000..7f0bec898 --- /dev/null +++ b/sorc/fvcom_tools.fd/CMakeLists.txt @@ -0,0 +1,22 @@ +set(fortran_src + kinds.f90 + module_ncio.f90 + module_nwp_base.f90 + module_nwp.f90 + process_FVCOM.f90) + + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -convert big_endian") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8 -fconvert=big-endian") +endif() + +set(exe_name fvcom_to_FV3) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries( + ${exe_name} + MPI::MPI_Fortran + NetCDF::NetCDF_Fortran) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/fvcom_tools.fd/kinds.f90 b/sorc/fvcom_tools.fd/kinds.f90 new file mode 100644 index 000000000..d10a7d797 --- /dev/null +++ b/sorc/fvcom_tools.fd/kinds.f90 @@ -0,0 +1,48 @@ +module kinds +!$$$ module documentation block +! . . . . +! module: kinds +! abstract: Module to hold specification kinds for variable declaration. +! This module is based on (copied from) Paul vanDelst's +! type_kinds module found in the community radiative transfer +! model +! +! module history log: +! +! Subroutines Included: +! +! Functions Included: +! +! remarks: +! The numerical data types defined in this module are: +! r_single - specification kind for single precision (4-byte) real variable +! i_kind - generic specification kind for default integer +! r_kind - generic specification kind for default floating point +! +! +! attributes: +! language: f90 +! +!$$$ end documentation block + implicit none + private +! +! for name string + integer, parameter, public :: len_sta_name = 8 + +! Integer type definitions below + +! Integer types + integer, parameter, public :: i_kind = 4 + integer, parameter, public :: i_short = 2 + integer, parameter, public :: i_byte = 1 +! Real types + integer, parameter, public :: r_single = 4 ! single precision + integer, parameter, public :: r_kind = 8 + +! + real(r_single),parameter,public :: rmissing=-99999.0 + real(i_kind),parameter,public :: imissing=-99999 + real(r_kind),parameter,public :: drmissing=-99999.0 + +end module kinds diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90 new file mode 100644 index 000000000..5f118ad93 --- /dev/null +++ b/sorc/fvcom_tools.fd/module_ncio.f90 @@ -0,0 +1,2360 @@ +module module_ncio +! +! module: functions to read and write netcdf files +! +! Ming Hu +! +! program history log: +! 2017-11-01 Hu initial build +! +! Subroutines Included: +! + + use netcdf + implicit none + + public :: ncio +! set default to private + private +! + type :: ncio + character(len=256) :: filename + integer :: ncid, status + integer :: debug_level + + integer :: nDims + integer :: ends(4) + integer :: xtype + character(len=40) :: dimname(4) + contains + procedure :: open => open_nc + procedure :: close => close_nc + +! read in dimension from the nc file + procedure :: get_dim => get_dim_nc + +! read in attribute from the nc file + generic :: get_att => get_att_nc_int,get_att_nc_real,get_att_nc_string + procedure :: get_att_nc_int + procedure :: get_att_nc_real + procedure :: get_att_nc_string + +! read in a 1d, 2d, 3d, or 4d field from the nc file + generic :: get_var => get_var_nc_double_1d, get_var_nc_double_2d, & + get_var_nc_double_3d, & + get_var_nc_real_1d,get_var_nc_real_2d, & + get_var_nc_real_3d, & + get_var_nc_short_1d,get_var_nc_short_2d, & + get_var_nc_int_1d,get_var_nc_int_2d, & + get_var_nc_int_3d, & + get_var_nc_char_1d,get_var_nc_char_2d, & + get_var_nc_char_3d + procedure :: get_var_nc_short + procedure :: get_var_nc_short_1d + procedure :: get_var_nc_short_2d + procedure :: get_var_nc_int + procedure :: get_var_nc_int_1d + procedure :: get_var_nc_int_2d + procedure :: get_var_nc_int_3d + procedure :: get_var_nc_real + procedure :: get_var_nc_real_1d + procedure :: get_var_nc_real_2d + procedure :: get_var_nc_real_3d + procedure :: get_var_nc_double + procedure :: get_var_nc_double_1d + procedure :: get_var_nc_double_2d + procedure :: get_var_nc_double_3d + procedure :: get_var_nc_char + procedure :: get_var_nc_char_1d + procedure :: get_var_nc_char_2d + procedure :: get_var_nc_char_3d + +! replace 1d, 2d, 3d, or 4d field from the nc file + generic :: replace_var => replace_var_nc_double_1d, replace_var_nc_double_2d, & + replace_var_nc_double_3d, & + replace_var_nc_real_1d,replace_var_nc_real_2d, & + replace_var_nc_real_3d, & + replace_var_nc_int_1d,replace_var_nc_int_2d, & + replace_var_nc_int_3d, & + replace_var_nc_char_1d,replace_var_nc_char_2d, & + replace_var_nc_char_3d + procedure :: replace_var_nc_int + procedure :: replace_var_nc_int_1d + procedure :: replace_var_nc_int_2d + procedure :: replace_var_nc_int_3d + procedure :: replace_var_nc_real + procedure :: replace_var_nc_real_1d + procedure :: replace_var_nc_real_2d + procedure :: replace_var_nc_real_3d + procedure :: replace_var_nc_double + procedure :: replace_var_nc_double_1d + procedure :: replace_var_nc_double_2d + procedure :: replace_var_nc_double_3d + procedure :: replace_var_nc_char + procedure :: replace_var_nc_char_1d + procedure :: replace_var_nc_char_2d + procedure :: replace_var_nc_char_3d + + procedure :: handle_err + + procedure :: convert_theta2t_2dgrid +!Add a new 3d variable to output file (David.M.Wright) + procedure :: add_new_var => add_new_var_3d + end type ncio + +contains + +subroutine open_nc(this,filename,action,debug_level) +! +! open a netcdf file, set initial debug level +! +! prgmmr: Ming Hu org: GSD date: 2017-11-01 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: filename + character(len=1),intent(in) :: action + integer,intent(in),optional :: debug_level + + integer :: ncid, status + + this%debug_level=20 + if(present(debug_level)) this%debug_level=debug_level + + this%filename=trim(filename) +! open existing netCDF dataset + if(action=="r" .or. action=="R") then + status = nf90_open(path = trim(filename), mode = nf90_nowrite, ncid = ncid) + elseif(action=="w" .or. action=="W") then + status = nf90_open(path = trim(filename), mode = nf90_write, ncid = ncid) + else + write(*,*) 'unknow action :', action + stop 123 + endif + if (status /= nf90_noerr) call this%handle_err(status) + this%ncid=ncid + + if(this%debug_level>0) then + write(*,*) '>>> open file: ',trim(this%filename) + endif + +end subroutine open_nc + +subroutine close_nc(this) +! +! initial netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-04-10 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + + integer :: ncid, status + + ncid=this%ncid +! +! close netCDF dataset + status = nf90_close(ncid) + if (status /= nf90_noerr) call this%handle_err(status) + +end subroutine close_nc + +subroutine get_att_nc_real(this,attname,rval) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + real, intent(out) :: rval + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), rval) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_real + +subroutine get_att_nc_int(this,attname,ival) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + integer, intent(out) :: ival + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), ival) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_int + +subroutine get_att_nc_string(this,attname,string) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + character(len=*), intent(out) :: string + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), string) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_string + + +subroutine get_dim_nc(this,dimname,dimvalue) +! +! get dimensions in netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-11-01 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*), intent(in) :: dimname + integer,intent(out) :: dimvalue + + integer :: ncid, status + integer :: DimId + +! open existing netCDF dataset + ncid=this%ncid + +! get dimension from exisiting NC file + status = nf90_inq_dimid(ncid,trim(dimname), DimId) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_Inquire_Dimension(ncid, DimId, len = dimvalue) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_dim_nc + +!==========================begin replace_var ========================== +subroutine replace_var_nc_char_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + character, intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_char_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + + call this%replace_var_nc_char(varname,ilength,field) +! +end subroutine replace_var_nc_char_1d + +subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + character, intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_char_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1) + endif +! + call this%replace_var_nc_char(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_char_2d + +subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + character, intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_char_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1,1) + endif + + call this%replace_var_nc_char(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_char_3d +! +subroutine replace_var_nc_char(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + character, intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_char' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_CHAR) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_char +!--- replace_var_nc_char + +!---- replace real +subroutine replace_var_nc_real_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(4), intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_real_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif +! + call this%replace_var_nc_real(varname,ilength,field) +! +end subroutine replace_var_nc_real_1d + +subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(4), intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_real_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_real(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_real_2d + +subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(4), intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_real_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_real(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_real_3d + +subroutine replace_var_nc_real(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(4), intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_real' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_FLOAT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_real + +!---- repalce double +subroutine replace_var_nc_double_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(8), intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_double_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif +! + call this%replace_var_nc_double(varname,ilength,field) +! +end subroutine replace_var_nc_double_1d + +subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(8), intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_double_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_double(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_double_2d + +subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(8), intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_double_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_double(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_double_3d +! + +subroutine replace_var_nc_double(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(8), intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_double' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_DOUBLE) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_double +! +!---- replace int +subroutine replace_var_nc_int_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer, intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_int_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + + call this%replace_var_nc_int(varname,ilength,field) +! +end subroutine replace_var_nc_int_1d + +subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer, intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_int_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_int(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_int_2d + +subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + integer, intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_int_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_int(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_int_3d +! +subroutine replace_var_nc_int(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer, intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_int' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_INT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_int +! +!==========================end of replace_var ========================== + +!==========================begin get_var ========================== +subroutine get_var_nc_double_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(8), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_double_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_double(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_double_1d + +subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(8), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_double_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_double(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_double_2d + +subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(8), intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_double_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_double(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_double_3d +! +subroutine get_var_nc_double(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(8), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_double' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_DOUBLE) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_DOUBLE,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(a,I10)') ' data type : ',this%xtype + write(*,'(a,I10)')' dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(a,I5,I10,2x,a)') ' rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_double + +subroutine get_var_nc_real_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(4), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_real_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_real(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_real_1d + +subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(4), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_real_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_real(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_real_2d + +subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(4), intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_real_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_real(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_real_3d +! +subroutine get_var_nc_real(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(4), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_real' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_FLOAT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_FLOAT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_real + +subroutine get_var_nc_int_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer, intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_int_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_int(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_int_1d + +subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer, intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_int_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_int(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_int_2d + +subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + integer, intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_int_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_int(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_int_3d +! +subroutine get_var_nc_int(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer, intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_int' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_INT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_int +! +subroutine get_var_nc_short_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer(2), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_short_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_short(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_short_1d +! +subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer(2), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer(2),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_short_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_short(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_short_2d +! +subroutine get_var_nc_short(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer(2), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_short' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_SHORT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_SHORT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_short + +subroutine get_var_nc_char_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + character, intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_char_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_char(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_char_1d + +subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + character, intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_char_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_char(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_char_2d + +subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + character, intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_char_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_char(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1,1) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_char_3d +! +subroutine get_var_nc_char(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + character, intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_char' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_CHAR) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_CHAR,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_char +!==========================end of get_var ========================== + +subroutine handle_err(this,status) + use netcdf + implicit none + class(ncio) :: this +! + integer, intent ( in) :: status + if(status /= nf90_noerr) then + print *, trim(nf90_strerror(status)) + stop "Stopped" + end if +end subroutine handle_err + +subroutine convert_theta2t_2dgrid(this,nx,ny,ps,t2) +! convertt theta T to T + implicit none + class(ncio) :: this + + integer :: nx,ny + real, intent(in ) :: ps(nx,ny) + real, intent(inout) :: t2(nx,ny) + + integer :: i,j + real(8) :: rd,cp,rd_over_cp + + + rd = 2.8705e+2_8 + cp = 1.0046e+3_8 ! specific heat of air @pressure (J/kg/K) + rd_over_cp = rd/cp + + do j=1,ny + do i=1,nx + t2(i,j)=t2(i,j)*(ps(i,j)/1000.0)**rd_over_cp - 273.15 + enddo + enddo + +end subroutine convert_theta2t_2dgrid + +subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units) +! +! prgmmr: David.M.Wright org: UM/GLERL date: 2020-09-01 +! +! abstract: Add a new variable to sfc_data.nc with dimensions +! (Time, yaxis_1, xaxis_1) +! +! Input: varname = Name of variable to be created in netcdf file +! dname1,dname2,dname3 = 1st,2nd, and 3rd dimension names +! lname = long name output for netcdf variable +! units = units to use in netcdf variable +! + use netcdf + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname,dname1,dname2,dname3 & + ,lname,units + integer :: status, ncid, dim1id, dim2id, dim3id, varid + + status = nf90_inq_dimid(this%ncid, dname1, dim1id) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_inq_dimid(this%ncid, dname2, dim2id) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_inq_dimid(this%ncid, dname3, dim3id) + if (status /= nf90_noerr) call this%handle_err(status) + + status = nf90_def_var(this%ncid, varname, nf90_double, & + (/ dim1id, dim2id, dim3id /), varid) + if (status /= nf90_noerr) call this%handle_err(status) + + status = nf90_put_att(this%ncid, varid, 'long_name', lname) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_put_att(this%ncid, varid, 'units', units) + if (status /= nf90_noerr) call this%handle_err(status) + +end subroutine add_new_var_3d + +end module module_ncio diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90 new file mode 100644 index 000000000..40965cdcf --- /dev/null +++ b/sorc/fvcom_tools.fd/module_nwp.f90 @@ -0,0 +1,280 @@ +! module_nwp.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This module defines FV3LAM and FVCOM forecast data structure and the method to +! read and write observations from and to those data structures. It is used by +! ingest_FVCOM.f90. +! +! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF +! to get FVCOM data into the model. + +module module_nwp + + use kinds, only: r_kind, r_single, i_short, rmissing + use module_nwp_base, only: nwpbase +! use module_map_utils, only: map_util + use module_ncio, only: ncio + + implicit none + + public :: fcst_nwp + public :: nwp_type + + private + type :: nwp_type + character(len=6) :: datatype + integer :: numvar, xlat, xlon, xtime + integer :: i_mask, i_sst, i_ice, i_sfcT, i_iceT + character(len=20), allocatable :: varnames(:) + character(len=20), allocatable :: latname + character(len=20), allocatable :: lonname + character(len=20), allocatable :: dimnameEW + character(len=20), allocatable :: dimnameNS + character(len=20), allocatable :: dimnameTIME + + real(r_kind), allocatable :: nwp_mask(:,:,:) + real(r_kind), allocatable :: nwp_sst(:,:,:) + real(r_kind), allocatable :: nwp_ice(:,:,:) + real(r_kind), allocatable :: nwp_sfcT(:,:,:) + real(r_kind), allocatable :: nwp_iceT(:,:,:) + end type nwp_type + + type, extends(nwp_type) :: fcst_nwp + type(nwpbase), pointer :: head => NULL() + type(nwpbase), pointer :: tail => NULL() + contains + procedure :: initial => initial_nwp + procedure :: list_initial => list_initial_nwp + procedure :: read_n => read_nwp + procedure :: finish => finish_nwp + end type fcst_nwp + + type(ncio) :: ncdata +! type(map_util) :: map + + contains + + subroutine initial_nwp(this,itype) + +! This subroutine defines the number of variables and their names for +! each NWP data type. The indices of the variables are +! also defined for later reference. + + class(fcst_nwp) :: this + + character(len=6), intent(in) :: itype + +! FVCOM grid + if (itype==' FVCOM') then + this%datatype = itype + this%numvar = 4 + + this%i_mask = 1 + this%i_sst = 2 + this%i_ice = 3 + this%i_iceT = 4 + this%i_sfcT = 0 + + allocate(this%varnames(this%numvar)) + this%varnames(1) = 'glmask' + this%varnames(2) = 'tsfc' + this%varnames(3) = 'aice' + this%varnames(4) = 'tisfc' + + allocate(this%latname) + allocate(this%lonname) + this%latname = 'lat' + this%lonname = 'lon' + + allocate(this%dimnameEW) + allocate(this%dimnameNS) + allocate(this%dimnameTIME) + this%dimnameEW = 'lon' + this%dimnameNS = 'lat' + this%dimnameTIME = 'Time' + +! FV3LAM grid + + else if (trim(itype)=='FV3LAM') then + this%datatype = itype + this%numvar = 4 + + this%i_mask = 1 + this%i_sst = 2 + this%i_ice = 3 + this%i_iceT = 4 + this%i_sfcT = 0 + + allocate(this%varnames(this%numvar)) + this%varnames(1) = 'slmsk' + this%varnames(2) = 'tsea' + this%varnames(3) = 'fice' + this%varnames(4) = 'tisfc' + + allocate(this%latname) + allocate(this%lonname) + this%latname = 'yaxis_1' + this%lonname = 'xaxis_1' + + allocate(this%dimnameEW) + allocate(this%dimnameNS) + allocate(this%dimnameTIME) + this%dimnameEW = 'xaxis_1' + this%dimnameNS = 'yaxis_1' + this%dimnameTIME = 'Time' + +! If the data type does not match one of the known types, exit. + + else + write(*,*) 'Unknown data type:', itype + stop 1234 + end if + + this%head => NULL() + this%tail => NULL() + + write(*,*) 'Finished initial_nwp' + write(*,*) ' ' + + end subroutine initial_nwp + + subroutine list_initial_nwp(this) + +! This subroutine lists the setup for NWP data that was done by +! the initial_nwp subroutine. + + class(fcst_nwp) :: this + + integer :: k + + write(*,*) 'List initial setup for ', this%datatype + write(*,*) 'number of variables ', this%numvar + write(*,*) 'variable index: mask, sst, ice, sfcT' + write(*,'(15x,10I3)') this%i_mask, this%i_sst, this%i_ice, & + & this%i_sfcT + write(*,*) 'variable name:' + do k=1,this%numvar + write(*,*) k,trim(this%varnames(k)) + enddo + + write(*,*) 'Finished list_initial_nwp' + write(*,*) ' ' + + end subroutine list_initial_nwp + + subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT) + +! This subroutine initializes arrays to receive the NWP data, +! and opens the file and gets the data. + + class(fcst_nwp) :: this + + character(len=5), intent(in) :: itype + character(len=*), intent(in) :: filename + + integer, intent(in) :: time_to_get + integer, intent(inout) :: numlon, numlat, numtimes +! real(r_single), intent(inout) :: mask(:,:), sst(:,:), ice(:,:), sfcT(:,:) + real(r_kind), intent(inout) :: mask(:,:),sst(:,:),ice(:,:),sfcT(:,:) & + ,iceT(:,:) + +! Open the file using module_ncio.f90 code, and find the number of +! lat/lon points + + call ncdata%open(trim(filename),'r',200) + call ncdata%get_dim(this%dimnameEW,this%xlon) + call ncdata%get_dim(this%dimnameNS,this%xlat) + call ncdata%get_dim(this%dimnameTIME,this%xtime) + + write(*,*) 'number of longitudes for file ', filename, this%xlon + numlon = this%xlon + write(*,*) 'number of latitudes for file ', filename, this%xlat + numlat = this%xlat + write(*,*) 'number of times for file ', filename, this%xtime + numtimes = this%xtime + +! Allocate all the arrays to receive data + + allocate(this%nwp_mask(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_sst(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_ice(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_sfcT(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_iceT(this%xlon,this%xlat,this%xtime)) + +! Get variables from the data file, but only if the variable is +! defined for that data type. + + if (this%i_mask .gt. 0) then + call ncdata%get_var(this%varnames(this%i_mask),this%xlon, & + this%xlat,this%xtime,this%nwp_mask) + mask = this%nwp_mask(:,:,1) + end if + if (this%i_sst .gt. 0) then + call ncdata%get_var(this%varnames(this%i_sst),this%xlon, & + this%xlat,this%xtime,this%nwp_sst) + sst = this%nwp_sst(:,:,time_to_get) + end if + if (this%i_ice .gt. 0) then + call ncdata%get_var(this%varnames(this%i_ice),this%xlon, & + this%xlat,this%xtime,this%nwp_ice) + ice = this%nwp_ice(:,:,time_to_get) + end if + if (this%i_sfcT .gt. 0) then + call ncdata%get_var(this%varnames(this%i_sfcT),this%xlon, & + this%xlat,this%xtime,this%nwp_sfcT) + sfcT = this%nwp_sfcT(:,:,time_to_get) + end if + if (this%i_iceT .gt. 0) then + call ncdata%get_var(this%varnames(this%i_iceT),this%xlon, & + this%xlat,this%xtime,this%nwp_iceT) + iceT = this%nwp_iceT(:,:,time_to_get) + end if + +! Close the netCDF file. + + call ncdata%close + + write(*,*) 'Finished read_nwp' + write(*,*) ' ' + + end subroutine read_nwp + + subroutine finish_nwp(this) + + class(fcst_nwp) :: this + + type(nwpbase), pointer :: thisobs,thisobsnext + + deallocate(this%varnames) + deallocate(this%latname) + deallocate(this%lonname) + deallocate(this%dimnameEW) + deallocate(this%dimnameNS) + deallocate(this%dimnameTIME) + deallocate(this%nwp_mask) + deallocate(this%nwp_sst) + deallocate(this%nwp_ice) + deallocate(this%nwp_sfcT) + deallocate(this%nwp_iceT) + + thisobs => this%head + if(.NOT.associated(thisobs)) then + write(*,*) 'No memory to release' + return + endif + do while(associated(thisobs)) +! write(*,*) 'destroy ==',thisobs%name + + thisobsnext => thisobs%next + call thisobs%destroy() + thisobs => thisobsnext + enddo + + write(*,*) 'Finished finish_nwp' + write(*,*) ' ' + + end subroutine finish_nwp + +end module module_nwp diff --git a/sorc/fvcom_tools.fd/module_nwp_base.f90 b/sorc/fvcom_tools.fd/module_nwp_base.f90 new file mode 100644 index 000000000..d32e841e4 --- /dev/null +++ b/sorc/fvcom_tools.fd/module_nwp_base.f90 @@ -0,0 +1,126 @@ +! module_nwp_base.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This module defines nwp observation data structure and the method to +! read and write observations from and to those data structures. It is used by +! ingest_FVCOM.f90. +! +! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF. +! + +module module_nwp_base + + use kinds, only: r_kind, r_single, rmissing + + implicit none + + public :: nwpbase + public :: nwplocation + + private + +! Define a nwp observation type. + + type nwplocation + real(r_single) :: lon ! stroke longitude + real(r_single) :: lat ! stroke latitiude + end type nwplocation + +! Define a nwp observation type to contain actual data. + + type, extends(nwplocation) :: nwpbase +! HOW DOES THIS POINTER THING WORK? + type(nwpbase), pointer :: next => NULL() + real(r_single) :: time ! observation time + integer :: numvar ! number of variables in this obs type +! real(r_single), allocatable :: obs(:) ! observation value (# numvar) + real(r_kind), allocatable :: obs(:) + logical :: ifquality ! do these obs include quality info? +! GLM has flash_quality_flag + integer, allocatable :: quality(:) ! if so, quality flags + contains + procedure :: list => list_obsbase + procedure :: alloc => alloc_obsbase + procedure :: destroy => destroy_obsbase + end type nwpbase + + contains + + subroutine list_obsbase(this) + +! This subroutine lists the contents of a base nwp observation + + class(nwpbase) :: this + + integer :: i, numvar + +! Write out the lon, lat, and time of the ob + + write(*,'(a,3f10.3)') 'LIGHTNING OB: longitude, latitude, time =', & + this%lon, this%lat, this%time + +! Loop through all variables and print out obs and quality + + numvar = this%numvar + if (numvar >= 1) then +! MULTI-DIMENSIONAL EXAMPLE IN module_obs_base.f90 + write(*,'(a4,10F12.2)') 'obs=', (this%obs(i),i=1,numvar) + if(this%ifquality) & + write(*,'(a4,10I12)') 'qul=', (this%quality(i),i=1,numvar) + else + write(*,*) 'No obs for this location' + endif + + end subroutine list_obsbase + + subroutine alloc_obsbase(this,numvar,ifquality) + +! This subroutine allocates memory for base nwp observation variables +! Input variables: +! numvar : number of variables in this ob type +! itquality: does this observation include quality information? + + class(nwpbase) :: this + + integer, intent(in) :: numvar + + logical, intent(in), optional :: ifquality + + if (numvar >= 1) then + this%numvar = numvar + + if(allocated(this%obs)) deallocate(this%obs) + allocate(this%obs(numvar)) + + this%ifquality=.false. + if(present(ifquality)) this%ifquality = ifquality + if(this%ifquality) allocate(this%quality(numvar)) + + else + write(*,*) 'alloc_obsbase Error: dimension must be larger than 0:', numvar + stop 1234 + endif + + end subroutine alloc_obsbase + + subroutine destroy_obsbase(this) + +! This subroutine releases memory associated with nwp observations + + class(nwpbase) :: this + + this%numvar = 0 + this%time = 0 + + if(allocated(this%obs)) deallocate(this%obs) + + this%ifquality=.false. + if(allocated(this%quality)) deallocate(this%quality) + + this%next => NULL() + + end subroutine destroy_obsbase + +end module module_nwp_base diff --git a/sorc/fvcom_tools.fd/process_FVCOM.f90 b/sorc/fvcom_tools.fd/process_FVCOM.f90 new file mode 100755 index 000000000..da3be6f9b --- /dev/null +++ b/sorc/fvcom_tools.fd/process_FVCOM.f90 @@ -0,0 +1,217 @@ +! process_FVCOM.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This is the code to put lake surface temperature and aerial ice +! concentration from GLERL-provided FVCOM forecast files (which have +! already been mapped to the FV3-LAM grid) into sfc_data.nc. +! +! This script will take two variables from the command line: +! 1. Name of FV3 sfc data file (e.g. sfc_data.tile7.halo0.nc) +! 2. Name of FVCOM data file (e.g. fvcom.nc) +! +! To run the script, use the following example, modifying file +! names as needed: +! ./fvcom_to_FV3 sfc_data.tile7.halo0.nc fvcom.nc +! +! Code is strongly based upon Eric James' (ESRL/GSL) work +! to update HRRR/WRF Great Lakes' temperature data with FVCOM. +! Code also relies heavily on Ming Hu's ncio module. + +program process_FVCOM + + use mpi + use kinds, only: r_kind, i_kind, r_single + use module_ncio, only: ncio + use module_nwp, only: fcst_nwp + + implicit none + +! MPI variables + integer :: npe, mype, mypeLocal,ierror +! + +! New object-oriented declarations + + type(ncio) :: geo + type(fcst_nwp) :: fcst + +! Grid variables + + character*180 :: geofile + character*2 :: workPath + character*1 :: char1 + + integer :: MAP_PROJ, NLON, NLAT + integer :: fv3lon, fv3lat, fv3times + integer :: lbclon, lbclat, lbctimes + integer :: i, j, t1, t2 + integer :: num_args, ix + + real :: rad2deg = 180.0/3.1415926 + real :: userDX, userDY, CEN_LAT, CEN_LON + real :: userTRUELAT1, userTRUELAT2, MOAD_CEN_LAT, STAND_LON + real :: truelat1, truelat2, stdlon, lat1, lon1, r_earth + real :: knowni, knownj, dx + real :: one, pi, deg2rad + + character(len=180) :: fv3file + character(len=180) :: fvcomfile + character(len=180), dimension(:), allocatable :: args + + real(r_kind), allocatable :: fv3ice(:,:), fv3sst(:,:) + real(r_kind), allocatable :: fv3sfcT(:,:), fv3mask(:,:) + real(r_kind), allocatable :: fv3iceT(:,:) + real(r_kind), allocatable :: lbcice(:,:), lbcsst(:,:) + real(r_kind), allocatable :: lbcsfcT(:,:), lbcmask(:,:) + real(r_kind), allocatable :: lbciceT(:,:) + +! Declare namelists +! SETUP (general control namelist) : + + integer :: update_type + +! namelist/setup/update_type, t2 + +! MPI setup + call MPI_INIT(ierror) + call MPI_COMM_SIZE(mpi_comm_world,npe,ierror) + call MPI_COMM_RANK(mpi_comm_world,mype,ierror) +! +! NCEP LSF has to use all cores allocated to run this application +! but this if check can make sure only one core run through the real code. +! +if(mype==0) then + +! Get file names from command line arguements + num_args = command_argument_count() + allocate(args(num_args)) + + do ix = 1, num_args + call get_command_argument(ix,args(ix)) + end do + + fv3file=trim(args(1)) + write(*,*) trim(fv3file) + fvcomfile=trim(args(2)) + write(*,*) trim(fvcomfile) + + t2 = 1 +! Obtain grid parameters + + workPath='./' + write(geofile,'(a,a)') trim(workPath), trim(fv3file) + write(*,*) 'sfc data file', trim(geofile) + call geo%open(trim(geofile),'r',200) + call geo%get_dim("xaxis_1",NLON) + call geo%get_dim("yaxis_1",NLAT) + write(*,*) 'NLON,NLAT:', NLON, NLAT + call geo%close + + write(*,*) 'Finished reading sfc_data grid information.' + write(*,*) ' ' + +! Allocate variables for I/O + + allocate(fv3ice(nlon,nlat)) + allocate(fv3sfcT(nlon,nlat)) + allocate(fv3sst(nlon,nlat)) + allocate(fv3mask(nlon,nlat)) + allocate(fv3iceT(nlon,nlat)) + + allocate(lbcice(nlon,nlat)) + allocate(lbcsfcT(nlon,nlat)) + allocate(lbcsst(nlon,nlat)) + allocate(lbcmask(nlon,nlat)) + allocate(lbciceT(nlon,nlat)) + +! Read fv3 sfc_data.nc before update + +! fv3file='sfc_data.nc' + t1=1 + + call fcst%initial('FV3LAM') + call fcst%list_initial + call fcst%read_n(trim(fv3file),'FV3LAM',fv3lon,fv3lat,fv3times,t1,fv3mask,fv3sst,fv3ice,fv3sfcT,fv3iceT) + call fcst%finish + + + write(*,*) 'fv3times: ', fv3times + write(*,*) 'time to use: ', t1 + +! Read FVCOM input datasets + +! fvcomfile='fvcom.nc' +! Space infront of ' FVCOM' below is important!! + call fcst%initial(' FVCOM') + call fcst%list_initial + call fcst%read_n(trim(fvcomfile),' FVCOM',lbclon,lbclat,lbctimes,t2,lbcmask,lbcsst,lbcice,lbcsfcT,lbciceT) + call fcst%finish + +! Check that the dimensions match + + if (lbclon .ne. nlon .or. lbclat .ne. nlat) then + write(*,*) 'ERROR: FVCOM/FV3 dimensions do not match:' + write(*,*) 'lbclon: ', lbclon + write(*,*) 'nlon: ', nlon + write(*,*) 'lbclat: ', lbclat + write(*,*) 'nlat: ', nlat + stop 135 + endif + + write(*,*) 'lbctimes: ', lbctimes + write(*,*) 'time to use: ', t2 + +! Update with FVCOM fields and process +! ice cover data. Ice fraction is set +! to a minimum of 15% due to FV3-LAM +! raising any value below 15% to 15%. + + do j=1,nlat + do i=1,nlon + if (lbcmask(i,j) > 0. .and. lbcsst(i,j) .ge. -90.0) then + !If ice fraction below 15%, set to 0 + if (lbcice(i,j) < 0.15) then + lbcice(i,j) = 0.0 + endif + fv3ice(i,j) = lbcice(i,j) + !If ice in FVCOM, but not in FV3-LAM, change to ice + if (lbcice(i,j) > 0. .and. fv3mask(i,j) == 0.) then + fv3mask(i,j) = 2. + endif + !If ice in FV3-LAM and not FVCOM, remove it from FV3-LAM + if (fv3mask(i,j) == 2. .and. lbcice(i,j) == 0.) then + fv3mask(i,j) = 0. + endif + fv3sst(i,j) = lbcsst(i,j) + 273.15 + fv3sfcT(i,j) = lbcsst(i,j) + 273.15 + fv3iceT(i,j) = lbcsst(i,j) + 273.15 + !If ice exists in FVCOM, change ice surface temp + if (lbcice(i,j) > 0.) then + fv3iceT(i,j) = lbciceT(i,j) + 273.15 + end if + endif + enddo + enddo + +! Write out sfc file again + + call geo%open(trim(fv3file),'w',300) + call geo%replace_var("tsea",NLON,NLAT,fv3sst) + call geo%replace_var("fice",NLON,NLAT,fv3ice) + call geo%replace_var("slmsk",NLON,NLAT,fv3mask) + call geo%replace_var("tisfc",NLON,NLAT,fv3iceT) +! Add_New_Var takes names of (Variable,Dim1,Dim2,Dim3,Long_Name,Units) + call geo%add_new_var('glmsk','xaxis_1','yaxis_1','Time','glmsk','none') + call geo%replace_var('glmsk',NLON,NLAT,lbcmask) + call geo%close + + write(6,*) "=== LOWBC UPDATE SUCCESS ===" + +endif ! mype==0 + +call MPI_FINALIZE(ierror) + + +end program process_FVCOM diff --git a/sorc/fvcom_tools.fd/readme.md b/sorc/fvcom_tools.fd/readme.md new file mode 100644 index 000000000..252c7ddf0 --- /dev/null +++ b/sorc/fvcom_tools.fd/readme.md @@ -0,0 +1,39 @@ +**fvcom_to_FV3.exe** + +**Introduction:** + This code replaces lake surface and lake ice temperature along + with aerial ice concentration generated from Great Lakes + Operational Forecast System (GLOFS), an FVCOM-based model, into + sfc_data.nc. + **NOTE** that the variables in the input files must reside on + the same grid. This means data from FVCOM must be horizontally + interpolated to the FV3 grid. This routine will also force a + minimum ice concentration of 15%. If ice concentration is less + than 15% in FVCOM, it will be set to 0% to avoid FV3 from + changing values less than 15% to 15% and generating unrealistic + lake ice temperatures. + +**Library Dependencies:** + Installation depends on the netCDF library and cmake. + +**Running:** + This routine will take two variables from the command line: + 1. Name of FV3 sfc data file (e.g. sfc_data.tile7.halo0.nc) + which is generated from chgres_cube.exe. + 2. Name of FVCOM data file in netcdf format (e.g. fvcom.nc) + + To run the script, use the following example, modifying file + names as needed: + ./fvcom_to_FV3 sfc_data.tile7.halo0.nc fvcom.nc + Output will be to the sfc data file and include lake surface + and lake ice temperature, and lake ice concentration from FVCOM. + + +This routine is *strongly* based upon Eric James' (ESRL/GSL) work + to update HRRR/WRF Great Lakes' temperature data with FVCOM. + It also relies heavily on Ming Hu's (ESRL/GSL) ncio module. + +**For more information, please contact:** + David Wright + University of Michigan and GLERL + dmwright@umich.edu