module_nwp.f90

Go to the documentation of this file.

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
      integer :: xlat
      integer :: xlon
      integer :: xtime
      integer :: datelen
      integer :: i_mask
      integer :: i_sst
      integer :: i_ice
      integer :: i_sfcT
      integer :: i_iceT
      integer :: i_sfcTl
      integer :: i_zorl
      integer :: i_hice
      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
      character(len=20), allocatable :: dimnameDATE
      character(len=1), allocatable :: times(:,:)

      real(r_kind), allocatable :: nwp_mask_c(:,:)
      real(r_kind), allocatable :: nwp_sst_c(:,:,:)
      real(r_kind), allocatable :: nwp_ice_c(:,:,:)
      real(r_kind), allocatable :: nwp_sfct_c(:,:,:)
      real(r_kind), allocatable :: nwp_icet_c(:,:,:)
      real(r_kind), allocatable :: nwp_zorl_c(:,:,:)
      real(r_kind), allocatable :: nwp_hice_c(:,:,:)

      real(r_kind), allocatable :: nwp_mask_w(:,:)
      real(r_kind), allocatable :: nwp_sst_w(:,:)
      real(r_kind), allocatable :: nwp_ice_w(:,:)
      real(r_kind), allocatable :: nwp_sfct_w(:,:)
      real(r_kind), allocatable :: nwp_icet_w(:,:)
      real(r_kind), allocatable :: nwp_sfctl_w(:,:)
      real(r_kind), allocatable :: nwp_zorl_w(:,:)
      real(r_kind), allocatable :: nwp_hice_w(:,:)

   end type nwp_type

   type, extends(nwp_type) :: fcst_nwp
      ! The pointers are carryover from when I inherited the code from
      ! GSL's work with HRRR for a similar use. I am not sure with
      ! object based coding in Fortran if it needs to have parts
      ! initialized to gain access to the procedures within it. - D. Wright.
      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 :: get_time_ind => get_time_ind_nwp 
         procedure :: finish => finish_nwp 
   end type fcst_nwp

   type(ncio) :: ncdata
!   type(map_util) :: map

   contains

      subroutine initial_nwp(this,itype,wcstart)
         class(fcst_nwp) :: this

         character(len=6), intent(in) :: itype
         character(len=4), intent(in) :: wcstart

!        FVCOM grid
         if (itype==' FVCOM') then
            this%datatype = itype
            this%numvar = 5

            this%i_mask = 1
            this%i_sst = 2
            this%i_ice = 3
            this%i_iceT = 4
            this%i_hice = 5
            this%i_sfcT = 0
            this%i_zorl = 0

            allocate(this%varnames(this%numvar))
            this%varnames(1) = 'glmask'
            this%varnames(2) = 'tsfc'
            this%varnames(3) = 'aice'
            this%varnames(4) = 'tisfc'
            this%varnames(5) = 'vice'

            allocate(this%latname)
            allocate(this%lonname)
            this%latname = 'lat'
            this%lonname = 'lon'

            allocate(this%dimnameEW)
            allocate(this%dimnameNS)
            allocate(this%dimnameTIME)
            allocate(this%dimnameDATE)
            this%dimnameEW = 'lon'
            this%dimnameNS = 'lat'
            this%dimnameTIME = 'Time'
            this%dimnameDATE = 'DateStrLen'

!        FV3LAM grid

         else if (trim(itype)=='FV3LAM' .AND. wcstart=='warm') then
            this%datatype = itype
            this%numvar = 8

            this%i_mask = 1
            this%i_sst = 2
            this%i_ice = 3
            this%i_iceT = 4
            this%i_sfcT = 5
            this%i_sfcTl= 6
            this%i_zorl = 7
            this%i_hice = 8

            allocate(this%varnames(this%numvar))
            this%varnames(1) = 'slmsk'
            this%varnames(2) = 'tsea'
            this%varnames(3) = 'fice'
            this%varnames(4) = 'tisfc'
            this%varnames(5) = 'tsfc'
            this%varnames(6) = 'tsfcl'
            this%varnames(7) = 'zorli'
            this%varnames(8) = 'hice'

            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'

         else if (trim(itype)=='FV3LAM' .AND. wcstart=='cold') then
            this%datatype = itype
            this%numvar = 6

            this%i_mask = 1
            this%i_sst = 2
            this%i_ice = 3
            this%i_iceT = 4
            this%i_zorl = 5
            this%i_hice = 6
            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'
            this%varnames(5) = 'zorl'
            this%varnames(6) = 'hice'

            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(6,*) 'Unknown data type:', itype
            stop 1234
         end if

         this%head => null()
         this%tail => null()

         write(6,*) 'Finished initial_nwp'
         write(6,*) ' '

      end subroutine initial_nwp

      subroutine list_initial_nwp(this)

         class(fcst_nwp) :: this

         integer :: k

         write(6,*) 'List initial setup for ', this%datatype
         write(6,*) 'number of variables ', this%numvar
         write(6,*) 'variable index: mask, sst, ice, sfcT, sfcTl'
         write(6,'(15x,10I3)') this%i_mask, this%i_sst, this%i_ice, &
      &      this%i_sfcT, this%i_sfcTl
         write(6,*) 'variable name:'
         do k=1,this%numvar
            write(6,*) k,trim(this%varnames(k))
         enddo

         write(6,*) 'Finished list_initial_nwp'
         write(6,*) ' '

      end subroutine list_initial_nwp

      subroutine read_nwp(this,filename,itype,wcstart,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT,sfcTl,zorl,hice,ybegin,yend)

         class(fcst_nwp) :: this

         character(len=6), intent(in) :: itype
         character(len=*), intent(in) :: filename
         character(len=4), intent(in) :: wcstart

         integer, intent(in) :: time_to_get
         integer, intent(in) :: ybegin,yend
         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(:,:),sfctl(:,:),zorl(:,:),hice(:,:)

!
!        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(6,*) 'number of longitudes for file ', filename, this%xlon
         numlon = this%xlon
         write(6,*) 'number of latitudes for file ', filename, this%xlat
         !numlat = this%xlat
         numlat = yend-ybegin+1
         write(6,*) 'number of times for file ', filename, this%xtime
         numtimes = this%xtime
         write(6,*) 'the range of Y for this domain is=',ybegin,yend

!        Allocate all the arrays to receive data
         if (wcstart == 'cold' .OR. itype == ' FVCOM') then
            allocate(this%nwp_mask_c(this%xlon,this%xlat))
            allocate(this%nwp_sst_c(this%xlon,this%xlat,this%xtime))
            allocate(this%nwp_ice_c(this%xlon,this%xlat,this%xtime))
            allocate(this%nwp_sfcT_c(this%xlon,this%xlat,this%xtime))
            allocate(this%nwp_iceT_c(this%xlon,this%xlat,this%xtime))
            allocate(this%nwp_zorl_c(this%xlon,this%xlat,this%xtime))
            allocate(this%nwp_hice_c(this%xlon,this%xlat,this%xtime))

!        Get variables from the data file, but only if the variable is
!        defined for that data type.

            write(6,*) 'itype = ', itype
            write(6,*) 'wcstart = ', wcstart
            write(6,*) 'xlat = ', this%xlat
            write(6,*) 'xlon = ', this%xlon
            write(6,*) 'xtime = ', this%xtime

            if (this%i_mask .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_mask),this%xlon,  &
                                   this%xlat,this%nwp_mask_c)
               mask = this%nwp_mask_c(:,ybegin:yend)
            end if
            if (this%i_sst .gt. 0) then
               write(6,*) 'get sst for cold or FVCOM'
               call ncdata%get_var(this%varnames(this%i_sst),this%xlon,  &
                                   this%xlat,this%xtime,this%nwp_sst_c)
               sst = this%nwp_sst_c(:,ybegin:yend,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_c)
               ice = this%nwp_ice_c(:,ybegin:yend,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_c)
               sfct = this%nwp_sfcT_c(:,ybegin:yend,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_c)
                icet = this%nwp_iceT_c(:,ybegin:yend,time_to_get)
            end if
            if (this%i_zorl .gt. 0) then
                call ncdata%get_var(this%varnames(this%i_zorl),this%xlon,  &
                                    this%xlat,this%xtime,this%nwp_zorl_c)
                zorl = this%nwp_zorl_c(:,ybegin:yend,time_to_get)
            end if 
            if (this%i_hice .gt. 0) then
                call ncdata%get_var(this%varnames(this%i_hice),this%xlon,  &
                                    this%xlat,this%xtime,this%nwp_hice_c)
                hice = this%nwp_hice_c(:,ybegin:yend,time_to_get)
            end if 
 
         else if (wcstart == 'warm') then
            allocate(this%nwp_mask_w(this%xlon,this%xlat))
            allocate(this%nwp_sst_w(this%xlon,this%xlat))
            allocate(this%nwp_ice_w(this%xlon,this%xlat))
            allocate(this%nwp_sfcT_w(this%xlon,this%xlat))
            allocate(this%nwp_iceT_w(this%xlon,this%xlat))
            allocate(this%nwp_sfcTl_w(this%xlon,this%xlat))
            allocate(this%nwp_zorl_w(this%xlon,this%xlat))
            allocate(this%nwp_hice_w(this%xlon,this%xlat))
!        Get variables from the data file, but only if the variable is
!        defined for that data type.
      
            write(6,*) 'itype = ', itype
            write(6,*) 'wcstart =', wcstart
            write(6,*) 'xlat = ', this%xlat
            write(6,*) 'xlon = ', this%xlon
            write(6,*) 'xtime = ', this%xtime

            if (this%i_mask .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_mask),this%xlon,  &
                                   this%xlat,this%nwp_mask_w)
               mask = this%nwp_mask_w(:,ybegin:yend)
            end if
            if (this%i_sst .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_sst),this%xlon,  &
                                   this%xlat,this%nwp_sst_w)
               sst = this%nwp_sst_w(:,ybegin:yend)
            end if
            if (this%i_ice .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_ice),this%xlon,  &
                                   this%xlat,this%nwp_ice_w)
               ice = this%nwp_ice_w(:,ybegin:yend)
            end if
            if (this%i_sfcT .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_sfcT),this%xlon,  &
                                   this%xlat,this%nwp_sfcT_w)
               sfct = this%nwp_sfcT_w(:,ybegin:yend)
            end if
            if (this%i_iceT .gt. 0) then
                call ncdata%get_var(this%varnames(this%i_iceT),this%xlon,  &
                                    this%xlat,this%nwp_iceT_w)
                icet = this%nwp_iceT_w(:,ybegin:yend)
            end if 
            if (this%i_sfcTl .gt. 0) then
               call ncdata%get_var(this%varnames(this%i_sfcTl),this%xlon,  &
                                   this%xlat,this%nwp_sfcTl_w)
               sfctl = this%nwp_sfcTl_w(:,ybegin:yend)
            end if
            if (this%i_zorl .gt. 0) then
                call ncdata%get_var(this%varnames(this%i_zorl),this%xlon,  &
                                    this%xlat,this%nwp_zorl_w)
                zorl = this%nwp_zorl_w(:,ybegin:yend)
            end if 
            if (this%i_hice .gt. 0) then
                call ncdata%get_var(this%varnames(this%i_hice),this%xlon,  &
                                    this%xlat,this%nwp_hice_w)
                hice = this%nwp_hice_w(:,ybegin:yend)
            end if 

         else
            write(6,*) 'Choose either "warm" or "cold" for file'
            stop 'Error in wcstart. Check spelling or if variable was assigned'
         end if
!        Close the netCDF file.

         call ncdata%close

         write(6,*) 'Finished read_nwp'
         write(6,*) ' '

      end subroutine read_nwp

      subroutine finish_nwp(this,itype,wcstart)

         class(fcst_nwp) :: this
         character(len=6), intent(in) :: itype
         character(len=4), intent(in) :: wcstart

         type(nwpbase), pointer :: thisobs,thisobsnext

         deallocate(this%varnames)
         deallocate(this%latname)
         deallocate(this%lonname)
         deallocate(this%dimnameEW)
         deallocate(this%dimnameNS)
         deallocate(this%dimnameTIME)
         if (wcstart == 'cold' .OR. itype==' FVCOM') then
            deallocate(this%nwp_mask_c)
            deallocate(this%nwp_sst_c)
            deallocate(this%nwp_ice_c)
            deallocate(this%nwp_sfcT_c)
            deallocate(this%nwp_iceT_c)
            deallocate(this%nwp_zorl_c)
            deallocate(this%nwp_hice_c)
            if (itype==' FVCOM') deallocate(this%dimnameDATE)
         else if (wcstart == 'warm') then
            deallocate(this%nwp_mask_w)
            deallocate(this%nwp_sst_w)
            deallocate(this%nwp_ice_w)
            deallocate(this%nwp_sfcT_w)
            deallocate(this%nwp_iceT_w)
            deallocate(this%nwp_sfcTl_w)
            deallocate(this%nwp_zorl_w)
            deallocate(this%nwp_hice_w)
         else
            write(6,*) 'no deallocation'
         end if

         thisobs => this%head
         if(.NOT.associated(thisobs)) then
            write(6,*) 'No memory to release'
            return
         endif
         do while(associated(thisobs))

            thisobsnext => thisobs%next
            call thisobs%destroy()
            thisobs => thisobsnext
         enddo

         write(6,*) 'Finished finish_nwp'
         write(6,*) ' '

      end subroutine finish_nwp

      subroutine get_time_ind_nwp(this,filename,instr,outindex)

         class(fcst_nwp) :: this

         character(len=*), intent(in) :: filename
         character(len=*), intent(in) :: instr
         integer, intent(out) :: outindex

         character(len=26) :: temp
         integer :: foundind
         integer :: k,i
 
!        Open the file using module_ncio.f90 code, and find the length of
!        time in the file
         call ncdata%open(trim(filename),'r',200)
         call ncdata%get_dim(this%dimnameTIME,this%xtime)
         call ncdata%get_dim(this%dimnameDATE,this%datelen)
         write(6,*) 'xtime = ', this%xtime
         write(6,*) 'datelen = ', this%datelen
         allocate(this%times(this%datelen,this%xtime))
         call ncdata%get_var('Times',this%datelen,this%xtime,this%times)

         foundind = 0

         do k=1,this%xtime,1
            do i = 1,len(temp),1
               temp(i:i) = this%times(i,k)
            end do
            if (trim(temp) == trim(instr)) then !If times are equal return k
              outindex = k
              foundind = 1
            end if
         end do
         if (foundind == 0) then
            outindex = -999
            deallocate(this%times)
            call ncdata%close
            write(6,*) 'WARNING: Supplied time not found in file: ', trim(instr)
            write(6,*) 'Stoppping fvcom_to_FV3 and proceeding without using FVCOM data'
            stop
         end if

         deallocate(this%times)
         call ncdata%close

      end subroutine get_time_ind_nwp

end module module_nwp