read_climate_nudge_data_nlm.F90

Go to the documentation of this file.

module read_climate_nudge_data_nlm_mod

use fms_mod, only: open_namelist_file, check_nml_error, close_file, &
                   stdlog, mpp_pe, mpp_root_pe, write_version_number, &
                   string, error_mesg, fatal, note, file_exist
use mpp_mod, only: input_nml_file
use mpp_io_mod,    only: mpp_open, mpp_netcdf, mpp_rdonly,mpp_multi, mpp_single
use mpp_io_mod,    only: axistype, fieldtype, mpp_get_time_axis, mpp_get_atts
use mpp_io_mod,    only: mpp_get_fields, mpp_get_info, mpp_get_axes, mpp_get_times
use mpp_io_mod,    only: mpp_get_axis_data, mpp_read, mpp_close, mpp_get_default_calendar
use constants_mod, only: pi, grav, rdgas, rvgas
use fv_arrays_nlm_mod,  only: real4, real8, fvprc

implicit none
private

public :: read_climate_nudge_data_init, read_time, read_grid,  &
          read_climate_nudge_data, read_climate_nudge_data_end
public :: read_sub_domain_init

interface read_climate_nudge_data
   module procedure read_climate_nudge_data_2d
   module procedure read_climate_nudge_data_3d
end interface

  real(FVPRC), parameter :: p0 = 1.e5
  real(FVPRC), parameter :: d608 = rvgas/rdgas - 1.

  integer, parameter :: num_req_axes = 3
  integer, parameter :: index_lon = 1, index_lat = 2, index_lev = 3
  character(len=8), dimension(NUM_REQ_AXES) :: required_axis_names = &
                                     (/ 'lon', 'lat', 'lev' /)

  integer, parameter :: num_req_flds = 9
  integer, parameter :: index_p0 = 1, index_ak = 2, index_bk = 3, &
                        index_zs = 4, index_ps = 5,               &
                        index_t  = 6, index_q  = 7,               &
                        index_u  = 8, index_v  = 9
  character(len=8), dimension(NUM_REQ_FLDS) :: required_field_names = &
       (/ 'P0  ', 'hyai', 'hybi', 'PHI ', 'PS  ', 'T   ', 'Q   ', 'U   ', 'V   ' /)
 
  integer, parameter :: maxfiles = 53
  character(len=256) :: filenames(maxfiles)
  character(len=256) :: filename_tails(maxfiles)
  character(len=256) :: filename_head
  integer :: read_buffer_size
  integer :: nfiles = 0
  logical :: module_is_initialized = .false.

  namelist /read_climate_nudge_data_nml/ filename_tails, read_buffer_size, &
                                         filename_head

! dimensions for checking
  integer :: global_axis_size(num_req_axes), numtime, sub_domain_latitude_size
  integer, allocatable :: file_index(:)

type filedata_type
  integer :: ncid
  integer,  pointer :: length_axes(:)  ! length of all dimensions in file
  integer :: ndim, nvar, natt, ntim, varid_time
  integer :: time_offset
  integer, dimension(NUM_REQ_FLDS) :: field_index   ! varid for variables
  integer, dimension(NUM_REQ_AXES) :: axis_index    ! varid for dimensions
  type(axistype),  dimension(NUM_REQ_FLDS) :: axes
  type(fieldtype), dimension(NUM_REQ_FLDS) :: fields  
end type

  type(filedata_type), allocatable :: files(:)

CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine read_climate_nudge_data_init (nlon, nlat, nlev, ntime)
integer, intent(out) :: nlon, nlat, nlev, ntime
! returns dimension lengths of input data set
!   nlon, nlat  lat/lon grid size
!   nlev        number of levels
!   ntime       number of time levels

  integer :: iunit, ierr, io
  character(len=128) :: name
  integer :: istat, i, j, k, n, nd, siz(4), i1, i2
  type(axistype), allocatable :: axes(:)
  type(fieldtype), allocatable :: fields(:)

  if (module_is_initialized) return
  ! initial file names to blanks
  do n = 1, maxfiles
     do i = 1, len(filename_tails(n))
        filename_tails(n)(i:i) = ' '
        filenames(n)(i:i) = ' '
     enddo
  enddo

!----- read namelist -----
#ifdef INTERNAL_FILE_NML
  read (input_nml_file, nml=read_climate_nudge_data_nml, iostat=io)
  ierr = check_nml_error(io, 'read_climate_nudge_data_nml')
#else
  if (file_exist('input.nml') ) then
    iunit = open_namelist_file()
    ierr=1
    do while (ierr /= 0)
      read (iunit, nml=read_climate_nudge_data_nml, iostat=io, end=10)
      ierr = check_nml_error(io, 'read_climate_nudge_data_nml')
    enddo
10  call close_file (iunit)
  endif
#endif

!----- write version and namelist to log file -----

  iunit = stdlog()
  call write_version_number ( '0.0', 'fv3-jedi-lm' )
  if (mpp_pe() == mpp_root_pe()) write (iunit, nml=read_climate_nudge_data_nml)

  ! determine the number of files
  do n = 1, maxfiles
     if (filename_tails(n)(1:1) .eq. ' ') exit
     nfiles = n
  enddo
  do n=1,nfiles
    filenames(n) = trim(filename_head)//trim(filename_tails(n))
  end do

  allocate(files(nfiles))
  numtime = 0

! open input file(s)
  do n = 1, nfiles
     call mpp_open(files(n)%ncid, trim(filenames(n)), form=mpp_netcdf,action=mpp_rdonly,threading=mpp_multi, &
             fileset=mpp_single )
!     call error_mesg ('read_climate_nudge_data_nlm_mod', 'Buffer size for reading = '//trim(string(read_buffer_size)), NOTE)

     call mpp_get_info(files(n)%ncid, files(n)%ndim, files(n)%nvar, files(n)%natt, files(n)%ntim)

     allocate (files(n)%length_axes(files(n)%ndim))
     allocate (axes(files(n)%ndim))
     call mpp_get_axes(files(n)%ncid,axes)

     ! inquire dimension sizes
     do i = 1, files(n)%ndim
        call mpp_get_atts(axes(i), name=name, len=files(n)%length_axes(i))
        do j = 1, num_req_axes
           if (trim(name) .eq. trim(required_axis_names(j))) then
              call check_axis_size (j,files(n)%length_axes(i))
              files(n)%axes(j) = axes(i)
              files(n)%axis_index(j) = i
              exit
           endif
        enddo
     enddo
     deallocate(axes)
     ! time axis indexing
     files(n)%time_offset = numtime
     numtime = numtime + files(n)%ntim

     allocate(fields(files(n)%nvar))
     call mpp_get_fields(files(n)%ncid,fields)
     files(n)%field_index = 0
     do i = 1, files(n)%nvar
        call mpp_get_atts(fields(i), name=name, ndim=nd, siz=siz)
        do j = 1, num_req_flds
           if (trim(name) .eq. trim(required_field_names(j))) then
              files(n)%field_index(j) = i
              files(n)%fields(j) = fields(i)
              if (j .gt. 3) then
                 call check_resolution (siz(1:nd))
              endif
              exit
           endif
        enddo

        ! special case for surface geopotential (sometimes the name is PHIS)
        ! z1l: the following is not needed
!       ! rab: if not needed, should we remove?
!        if (trim(name) .eq. 'PHIS') then
!           Files(n)%field_index(INDEX_ZS) = i
!           call check_resolution (siz(1:nd))
!        endif
     enddo
     deallocate(fields)

  enddo ! "n" files loop

  ! setup file indexing
  allocate(file_index(numtime))
  i2 = 0
  do n = 1, nfiles
     i1 = i2+1
     i2 = i2+files(n)%ntim
     file_index(i1:i2) = n
  enddo

      sub_domain_latitude_size = global_axis_size(index_lat)

    ! output arguments
      nlon = global_axis_size(index_lon)
      nlat = global_axis_size(index_lat)
      nlev = global_axis_size(index_lev)
      ntime = numtime

      module_is_initialized = .true.

end subroutine read_climate_nudge_data_init

!###############################################################################

subroutine read_time ( times, units, calendar )
real(FVPRC),          intent(out) :: times(:)
character(len=*), intent(out) :: units, calendar
integer :: istat, i1, i2, n
type(axistype), save:: time_axis
character(len=32) :: default_calendar

   if (.not.module_is_initialized) then
     call error_mesg ('read_climate_nudge_data_nlm_mod/read_time',  &
                                        'module not initialized', fatal)
   endif

   if (size(times(:)) < numtime) then
      call error_mesg ('read_climate_nudge_data_nlm_mod', 'argument times too small in read_time', fatal)
   endif

 ! data
   i2 = 0
   do n = 1, nfiles
      i1 = i2+1
      i2 = i2+files(n)%ntim
      if( n == 1) then
         default_calendar = mpp_get_default_calendar()
         call mpp_get_time_axis( files(n)%ncid, time_axis)
         call mpp_get_atts(time_axis, units=units, calendar=calendar)
         if( trim(calendar) == trim(default_calendar)) calendar = 'gregorian  '
      endif
      call mpp_get_times(files(n)%ncid, times(i1:i2))
   enddo

! NOTE: need to do the conversion to time_type in this routine
!       this will allow different units and calendars for each file

end subroutine read_time

!###############################################################################

subroutine read_grid ( lon, lat, ak, bk )
real(FVPRC), intent(out), dimension(:) :: lon, lat, ak, bk

 real(FVPRC) :: pref
 integer :: istat

   if (.not.module_is_initialized) then
     call error_mesg ('read_climate_nudge_data_nlm_mod/read_grid',  &
                                        'module not initialized', fatal)
   endif


    ! static fields from first file only
      call mpp_get_axis_data(files(1)%axes(index_lon), lon)
      call mpp_get_axis_data(files(1)%axes(index_lat), lat)

    ! units are assumed to be degrees east and north
    ! convert to radians
      lon = lon * pi/180.
      lat = lat * pi/180.

    ! vertical coodinate
      if (files(1)%field_index(index_ak) .gt. 0) then
         call mpp_read(files(1)%ncid, files(1)%fields(index_ak), ak)
         if (files(1)%field_index(index_p0) .gt. 0) then
            call mpp_read(files(1)%ncid, files(1)%fields(index_p0), pref)
         else
            pref = p0
         endif
         ak = ak*pref
      else
         ak = 0.
      endif
 
      call mpp_read(files(1)%ncid, files(1)%fields(index_bk), bk)


end subroutine read_grid

!###############################################################################

subroutine read_sub_domain_init ( ylo, yhi, ydat, js, je )
 real(FVPRC),    intent(in)  :: ylo, yhi, ydat(:)
 integer, intent(out) :: js, je
 integer :: j

   if (.not.module_is_initialized) then
     call error_mesg ('read_climate_nudge_data_nlm_mod/read_sub_domain_init',  &
                                        'module not initialized', fatal)
   endif
   ! increasing data
   if (ydat(1) < ydat(2)) then
      js = 1
      do j = 1, size(ydat(:))-1
         if (ylo >= ydat(j) .and. ylo <= ydat(j+1)) then
            js = j
            exit
         endif
      enddo

      if (ylo < -1.5) then
         print *, 'ylo=',ylo
         print *, 'js,ydat=',js,ydat(js)
         print *, 'ydat=',ydat(:js+2)
      endif

      je = size(ydat(:))
      do j = js, size(ydat(:))-1
         if (yhi >= ydat(j) .and. yhi <= ydat(j+1)) then
            je = j+1
            exit
         endif
      enddo

      if (yhi > 1.5) then
         print *, 'yhi=',yhi
         print *, 'je,ydat=',je,ydat(je)
         print *, 'ydat=',ydat(je-2:)
      endif

   ! decreasing data (may not work)
   else
      call error_mesg ('read_climate_nudge_data_nlm_mod', 'latitude values for observational data decrease with increasing index', note)
      je = size(ydat(:))-1
      do j = 1, size(ydat(:))-1
         if (ylo >= ydat(j+1) .and. ylo <= ydat(j)) then
            je = j+1
            exit
         endif
      enddo

      js = 1
      do j = 1, je
         if (yhi >= ydat(j+1) .and. yhi <= ydat(j)) then
            js = j
            exit
         endif
      enddo

   endif

   sub_domain_latitude_size = je-js+1

 end subroutine read_sub_domain_init

!###############################################################################

subroutine read_climate_nudge_data_2d (itime, field, dat, is, js)
integer,          intent(in) :: itime
character(len=4), intent(in) :: field
real(FVPRC),             intent(out), dimension(:,:) :: dat
integer,          intent(in),  optional       :: is, js
integer :: istat, atime, n, this_index
integer :: nread(4), start(4)

   if (.not.module_is_initialized) then
     call error_mesg ('read_climate_nudge_data_nlm_mod',  &
                                        'module not initialized', fatal)
   endif
     ! time index check
      if (itime < 1 .or. itime > numtime) then
         call error_mesg ('read_climate_nudge_data_nlm_mod', 'itime out of range', fatal)
      endif

     ! check dimensions 
     if (present(js)) then
        if (size(dat,1) .ne. global_axis_size(index_lon) .or. &
            size(dat,2) .ne. sub_domain_latitude_size) then
            !write (*,'(a)') 'climate_nudge_data_mod: size dat2d = '//trim(string(size(dat,1)))//' x '//trim(string(size(dat,2)))// &
            !         '  <-vs->  '//trim(string(global_axis_size(INDEX_LON)))//' x '//trim(string(sub_domain_latitude_size))
            call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect 2d array dimensions', fatal)
        endif
     else
        if (size(dat,1) .ne. global_axis_size(index_lon) .or. &
            size(dat,2) .ne. global_axis_size(index_lat))     &
            call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect 2d array dimensions', fatal)
     endif

     ! check field
     if (field .eq. 'phis') then
        this_index = index_zs
     else if (field .eq. 'psrf') then
        this_index = index_ps
     else
         call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect field requested in read_climate_nudge_data_2d', fatal)
     endif
     
     ! file index and actual time index in file
     n = file_index(itime)
     atime = itime - files(n)%time_offset

     start = 1
     if (present(is)) start(1) = is
     if (present(js)) start(2) = js
     start(3) = atime

     nread = 1
     nread(1) = size(dat,1)
     nread(2) = size(dat,2)
     
     call mpp_read(files(n)%ncid, files(n)%fields(this_index), dat, start, nread)
  
      ! geopotential height (convert to m2/s2 if necessary)
     if (field .eq. 'phis') then
        if (maxval(dat) > 1000.*grav) then
          ! do nothing
        else
          dat = dat * grav
        endif
     endif

end subroutine read_climate_nudge_data_2d

!###############################################################################

subroutine read_climate_nudge_data_3d (itime, field, dat, is, js)
integer,          intent(in) :: itime
character(len=4), intent(in) :: field
real(FVPRC),             intent(out), dimension(:,:,:) :: dat
integer,          intent(in),  optional         :: is, js
integer :: istat, atime, n, this_index, start(4), nread(4)
!logical :: convert_virt_temp = .false.

   if (.not.module_is_initialized) then
     call error_mesg ('read_climate_nudge_data_nlm_mod',  &
                                        'module not initialized', fatal)
   endif

     ! time index check
     if (itime < 1 .or. itime > numtime) then
        call error_mesg ('read_climate_nudge_data_nlm_mod', 'itime out of range', fatal)
     endif

     ! check dimensions
     if (present(js)) then
        if (size(dat,1) .ne. global_axis_size(index_lon) .or. &
            size(dat,2) .ne. sub_domain_latitude_size    .or. &
            size(dat,3) .ne. global_axis_size(index_lev)) then
            !write (*,'(a)') 'climate_nudge_data_mod: size dat3d = '//trim(string(size(dat,1)))//' x '//trim(string(size(dat,2)))// &
            !                                        ' x '//trim(string(size(dat,3)))
            call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect 3d array dimensions', fatal)
        endif
     else
        if (size(dat,1) .ne. global_axis_size(index_lon) .or. &
            size(dat,2) .ne. global_axis_size(index_lat) .or. &
            size(dat,3) .ne. global_axis_size(index_lev))     &
            call error_mesg ('read_climate_nudge_mod', 'incorrect 3d array dimensions', fatal)
     endif

     ! check field
     if (field .eq. 'temp') then
        this_index = index_t
     else if (field .eq. 'qhum') then
        this_index = index_q
     else if (field .eq. 'uwnd') then
        this_index = index_u
     else if (field .eq. 'vwnd') then
        this_index = index_v
     else
        call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect field requested in read_climate_nudge_data_3d', fatal)
     endif
     

     ! file index and actual time index in file
     n = file_index(itime)
     atime = itime - files(n)%time_offset

     start = 1
     if (present(is)) start(1) = is
     if (present(js)) start(2) = js
     start(4) = atime

     nread = 1
     nread(1) = size(dat,1)
     nread(2) = size(dat,2)
     nread(3) = size(dat,3)

     call mpp_read(files(n)%ncid, files(n)%fields(this_index), dat, start, nread)

     ! convert virtual temp to temp
     ! necessary for some of the high resol AVN analyses
     !if (convert_virt_temp) then
     !   temp = temp/(1.+D608*qhum)
     !endif

end subroutine read_climate_nudge_data_3d

!###############################################################################

subroutine read_climate_nudge_data_end
integer :: istat, n

  if ( .not.module_is_initialized) return
  do n = 1, nfiles
     call mpp_close(files(n)%ncid)
  enddo
  deallocate (files)
  module_is_initialized = .false.

end subroutine read_climate_nudge_data_end

!###############################################################################

 subroutine check_axis_size (ind,lendim)
 integer, intent(in) :: ind,lendim

   ! once the axis size is set all subsuquent axes must be the same
   if (global_axis_size(ind) .gt. 0) then
      if (global_axis_size(ind) .ne. lendim) then
         call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect axis size for axis = '//trim(required_axis_names(ind)), fatal)
      endif
   else
      global_axis_size(ind) = lendim
   endif

 end subroutine check_axis_size

!------------------------------------

 subroutine check_resolution (axis_len)
 integer, intent(in) :: axis_len(:)

   if (size(axis_len(:)) .lt. 2) then
      call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect number of array dimensions', fatal)
   endif
   if (axis_len(1) .ne. global_axis_size(index_lon)) then
      call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect array dimension one', fatal)
   endif
   if (axis_len(2) .ne. global_axis_size(index_lat)) then
      call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect array dimension two', fatal)
   endif
   if (size(axis_len(:)) .gt. 3) then
      if (axis_len(3) .ne. global_axis_size(index_lev)) then
         call error_mesg ('read_climate_nudge_data_nlm_mod', 'incorrect array dimension three', fatal)
      endif
   endif

 end subroutine check_resolution

!###############################################################################

end module read_climate_nudge_data_nlm_mod