From 947477ae459552b67c7452daf200d6bed7a1c4d4 Mon Sep 17 00:00:00 2001 From: David Wright Date: Tue, 29 Sep 2020 13:22:35 +0000 Subject: [PATCH] Add FVCOM lake surface conditions to UFS (TRY 2) Attempt #2 to resolve out of sync issues with repository. This feature takes lake surface data generated by FVCOM and replaces variables in sfc_data.nc (created by chgres_cube) with it. Input data into fvcom_to_FV3.exe needs to be in .nc format and already horizontally interpolated to the destination grid used in sfc_data.nc. Some QC is done within the code to create a minimum lake ice value of 15%. This code is needed to run PR #293 in the regional_workflow. Fortran code is strongly based upon work done by Eric James (GSL) and Ming Hu (GSL) to get FVCOM surface conditions into HRRRv4. Testing has been conducted on Jet using Intel compiler (18.0.5.274) --- sorc/CMakeLists.txt | 1 + sorc/fvcom_tools.fd/CMakeLists.txt | 22 + sorc/fvcom_tools.fd/kinds.f90 | 48 + sorc/fvcom_tools.fd/module_ncio.f90 | 2360 +++++++++++++++++++++++ sorc/fvcom_tools.fd/module_nwp.f90 | 280 +++ sorc/fvcom_tools.fd/module_nwp_base.f90 | 126 ++ sorc/fvcom_tools.fd/process_FVCOM.f90 | 217 +++ sorc/fvcom_tools.fd/readme.md | 39 + 8 files changed, 3093 insertions(+) create mode 100644 sorc/fvcom_tools.fd/CMakeLists.txt create mode 100644 sorc/fvcom_tools.fd/kinds.f90 create mode 100644 sorc/fvcom_tools.fd/module_ncio.f90 create mode 100644 sorc/fvcom_tools.fd/module_nwp.f90 create mode 100644 sorc/fvcom_tools.fd/module_nwp_base.f90 create mode 100755 sorc/fvcom_tools.fd/process_FVCOM.f90 create mode 100644 sorc/fvcom_tools.fd/readme.md 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