mod_realinput_nicam.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
module mod_realinput_nicam
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_tracer
  use scale_process, only: &
     myrank => prc_myrank,  &
     prc_mpistop
  use scale_external_io, only: &
     inicam

  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: parentatomsetupnicam
  public :: parentatomopennicam
  public :: parentatominputnicam
  public :: parentlandsetupnicam
  public :: parentlandinputnicam
  public :: parentoceansetupnicam
  public :: parentoceanopennicam
  public :: parentoceaninputnicam

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), parameter :: missval_tg   = 50.0_rp
  real(RP), parameter :: missval_strg = 0.0_rp
  real(RP), parameter :: missval_sst  = 273.1506_rp


  real(SP), allocatable :: read1DX(:)
  real(SP), allocatable :: read1DY(:)
  real(SP), allocatable :: read1DZ(:)
  real(SP), allocatable :: read3DS(:,:,:,:)
  real(SP), allocatable :: read4D (:,:,:,:)

  real(SP), allocatable :: read1DLZ(:)
  real(SP), allocatable :: read4DL(:,:,:,:)

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine parentatomsetupnicam( &
      dims,    &
      timelen, &
      basename_org )
    use gtool_file, only: &
         filegetshape
    implicit none

    integer,          intent(out) :: dims(6)
    integer,          intent(out) :: timelen
    character(len=*), intent(in)  :: basename_org

    character(len=H_LONG) :: basename
    integer :: dims_ncm(4)

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

    if( io_l ) write(io_fid_log,*) '+++ Real Case/Atom Input File Type: NICAM-NETCDF'
    basename = "ms_pres"//trim(basename_org)
    call filegetshape( dims_ncm(:), trim(basename), "ms_pres", 1, single=.true. )
    timelen = dims_ncm(4)

    ! full level
    dims(1) = dims_ncm(3)
    dims(2) = dims_ncm(1)
    dims(3) = dims_ncm(2)
    ! half level
    ! nicam lat-lon data doesn't have staggered grid system
    dims(4) = dims(1)
    dims(5) = dims(2)
    dims(6) = dims(3)

    allocate( read1dx( dims(2)                      ) )
    allocate( read1dy( dims(3)                      ) )
    allocate( read1dz( dims(1)                      ) )
    allocate( read3ds( 1,       dims(2), dims(3), 1 ) )
    allocate( read4d( dims(1), dims(2), dims(3), 1 ) )

    return
  end subroutine parentatomsetupnicam

  !-----------------------------------------------------------------------------
  subroutine parentatomopennicam( &
       lon_org,      &
       lat_org,      &
       cz_org,       &
       basename_num, &
       dims          )
    use scale_const, only: &
         d2r => const_d2r
    use gtool_file, only: &
         fileread
    use scale_external_io, only: &
         externalfileread
    implicit none

    real(RP),         intent(out) :: lon_org(:,:)
    real(RP),         intent(out) :: lat_org(:,:)
    real(RP),         intent(out) :: cz_org(:,:,:)
    character(len=*), intent(in)  :: basename_num
    integer,          intent(in)  :: dims(6)

    character(len=H_LONG) :: basename

    integer :: k, i, j

    if( io_l ) write(io_fid_log,*) '+++ ScaleLib/IO[realinput]/Categ[AtomOpenNICAM]'

    basename = "ms_pres"//trim(basename_num)
    call fileread( read1dx(:), trim(basename), "lon", 1, 1, single=.true. )
    do j = 1, dims(3)
       lon_org(:,j)  = read1dx(:) * d2r
    enddo

    call fileread( read1dy(:), trim(basename), "lat", 1, 1, single=.true. )
    do i = 1, dims(2)
       lat_org(i,:)  = read1dy(:) * d2r
    enddo

    call fileread( read1dz(:), trim(basename), "lev", 1, 1, single=.true. )
    do j = 1, dims(3)
    do i = 1, dims(2)
       cz_org(3:,i,j) = read1dz(:)
    end do
    end do

    call externalfileread( read4d(:,:,:,:), trim(basename), "ms_pres", 1, 1, myrank, inicam, single=.true. )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1)+2, 3, -1
          ! missing data implies under ground
          ! For more accurate surface height, we need topograph data
          ! So surface data is not used at this moment
          if ( read4d(k-2,i,j,1) <= 0.0_rp ) then
!             cz_org(k,i,j) = max( ( cz_org(k+1,i,j) + cz_org(k,i,j) ) * 0.5_RP, 0.0_RP )
             cz_org(1:k,i,j) = 0.0_rp
             exit
          endif
       enddo
    enddo
    enddo

    return
  end subroutine parentatomopennicam

  !-----------------------------------------------------------------------------
  subroutine parentatominputnicam( &
       velz_org,     &
       velx_org,     &
       vely_org,     &
       pres_org,     &
       temp_org,     &
       qtrc_org,     &
       basename_num, &
       dims,         &
       it            )
    use scale_const, only: &
       rdry => const_rdry, &
       cpdry => const_cpdry, &
       p00 => const_pre00, &
       d2r => const_d2r, &
       eps => const_eps
    use scale_external_io, only: &
         externalfileread, &
         externalfilereadoffset
    use scale_atmos_hydrometeor, only: &
       i_qv
    implicit none

    real(RP),         intent(out) :: velz_org(:,:,:)
    real(RP),         intent(out) :: velx_org(:,:,:)
    real(RP),         intent(out) :: vely_org(:,:,:)
    real(RP),         intent(out) :: pres_org(:,:,:)
    real(RP),         intent(out) :: temp_org(:,:,:)
    real(RP),         intent(out) :: qtrc_org(:,:,:,:)
    character(len=*), intent(in)  :: basename_num
    integer,          intent(in)  :: dims(6)
    integer,          intent(in)  :: it

    real(RP) :: tsfc_org (dims(2),dims(3))
    real(RP) :: slp_org  (dims(2),dims(3))
    real(RP) :: qvsfc_org(dims(2),dims(3))

    real(RP) :: pott
    real(RP) :: rovcp
    real(RP) :: cpovr

    integer :: k, i, j

    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '+++ ScaleLib/IO[realinput]/Categ[AtomInputNICAM]'

    basename = "ms_u"//trim(basename_num)
    call externalfilereadoffset( read4d(:,:,:,:), trim(basename), "ms_u", it, it, myrank, inicam, single=.true. )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = 1, dims(1)
          if ( abs(abs(read4d(k,i,j,1))-300.0_sp) < sqrt(eps) ) then ! missing value
             velx_org(k+2,i,j) = 0.0_rp
          else
             velx_org(k+2,i,j) = real( read4D(k,i,j,1), kind=rp )
          end if
       end do
       velx_org(1:2,i,j) = 0.0_rp
    end do
    end do

    basename = "ms_v"//trim(basename_num)
    call externalfilereadoffset( read4d(:,:,:,:),  &
                                 trim(basename),   &
                                 "ms_v",           &
                                 it, it,      &
                                 myrank,           &
                                 inicam,           &
                                 single=.true.     )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = 1, dims(1)
          if ( abs(abs(read4d(k,i,j,1))-300.0_sp) < sqrt(eps) ) then ! missing value
             vely_org(k+2,i,j) = 0.0_rp
          else
             vely_org(k+2,i,j) = real( read4D(k,i,j,1), kind=rp )
          end if
       end do
       vely_org(1:2,i,j) = 0.0_rp
    end do
    end do

    velz_org(:,:,:) = 0.0_rp


    ! The surface height is not available, so t2 is not used, at this moment.
!    basename = "ss_t2m"//trim(basename_num)
!    call ExternalFileRead( read3DS(:,:,:,:), trim(basename), "ss_t2m", it, it, myrank, iNICAM, single=.true. )
!    tsfc_org(:,:) = real( read3DS(1,:,:,1), kind=RP )
    basename = "ms_tem"//trim(basename_num)
    call externalfilereadoffset( read4d(:,:,:,:), trim(basename), "ms_tem", it, it, myrank, inicam, single=.true. )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read4d(k,i,j,1) <= 50.0_sp ) then ! missing value
!             temp_org(k+2,i,j) = tsfc_org(i,j)
             exit
          else
             temp_org(1:k+2,i,j) = real( read4D(k,i,j,1), kind=rp )
          end if
!          if( k==1 ) temp_org(2,i,j) = tsfc_org(i,j) ! no missing value case
       end do
    end do
    end do

    rovcp = rdry/cpdry
    cpovr = cpdry/rdry

    basename = "ss_slp"//trim(basename_num)
    call externalfilereadoffset( read3ds(:,:,:,:), trim(basename), "ss_slp", it, it, myrank, inicam, single=.true. )
    slp_org(:,:) = real( read3DS(1,:,:,1), kind=rp )
    basename = "ms_pres"//trim(basename_num)
    call externalfileread( read4d(:,:,:,:), trim(basename), "ms_pres", it, it, myrank, inicam, single=.true. )
    pres_org(3:,:,:) = read4d(:,:,:,1)
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read4d(k,i,j,1) < 0.0_sp ) then ! missing value
             exit
          else
             pres_org(k+2,i,j) = real( read4D(k,i,j,1), kind=rp )
          end if
       end do
       ! If data has missing value, k is the top layer having missing value,
       ! otherwise k is zero.
       pott = temp_org(k+3,i,j) * (p00/pres_org(k+3,i,j))**rovcp ! lowest level
!       pres_org(k+2,i,j) = P00 * (temp_org(k+2,i,j)/pott)**CPovR ! surface
       pres_org(1:k+2,i,j) = slp_org(i,j)                        ! sea level
       temp_org(1:k+2,i,j) = pott * (slp_org(i,j)/p00)**rovcp    ! sea level
    end do
    end do

    basename = "ss_q2m"//trim(basename_num)
    call externalfileread( read3ds(:,:,:,:), trim(basename), "ss_q2m", it, it, myrank, inicam, single=.true. )
    qvsfc_org(:,:) = real( read3DS(1,:,:,1), kind=rp )
    qtrc_org(:,:,:,:) = 0.0_rp
    basename = "ms_qv"//trim(basename_num)
    call externalfileread( read4d(:,:,:,:), trim(basename), "ms_qv", it, it, myrank, inicam, single=.true. )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read4d(k,i,j,1) < 0.0_sp ) then ! missing value
             qtrc_org(1:k+2,i,j,i_qv) = qvsfc_org(i,j) ! surface and sealevel
             exit
          else
             qtrc_org(k+2,i,j,i_qv) = real( read4D(k,i,j,1), kind=rp )
          end if
       end do
       qtrc_org(1:2,i,j,i_qv) = qvsfc_org(i,j)
    end do
    end do


    return
  end subroutine parentatominputnicam

  !-----------------------------------------------------------------------------
  subroutine parentlandsetupnicam( &
       ldims,    &
       basename_org )
    use gtool_file, only: &
         filegetshape
    implicit none

    integer,          intent(out) :: ldims(3)
    character(len=*), intent(in)  :: basename_org

    character(len=H_LONG) :: basename
    integer :: dims_ncm(4)

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

    if( io_l ) write(io_fid_log,*) '+++ Real Case/Land Input File Type: NICAM-NETCDF'
    basename = "la_tg"//trim(basename_org)
    call filegetshape( dims_ncm(:), trim(basename), "la_tg", 1, single=.true. )
    ! land
    ldims(1) = dims_ncm(3) ! vertical grid of land model
    ldims(2) = dims_ncm(1)
    ldims(3) = dims_ncm(2)

    if ( .not. allocated(read1dx) ) then
       allocate( read1dx( ldims(2)                       ) )
       allocate( read1dy( ldims(3)                       ) )
       allocate( read3ds( 1,       ldims(2), ldims(3), 1 ) )
    end if

    allocate( read1dlz( ldims(1)                      ) )
    allocate( read4dl( ldims(1), ldims(2), ldims(3), 1 ) )

    return
  end subroutine parentlandsetupnicam

  !-----------------------------------------------------------------------------
  subroutine parentlandinputnicam( &
      tg_org,             &
      strg_org,           &
      lst_org,            &
      llon_org,           &
      llat_org,           &
      lz_org,             &
      topo_org,           &
      lmask_org,          &
      basename_num,       &
      ldims,              &
      use_file_landwater, &
      it                  )
    use scale_const, only: &
         undef => const_undef, &
         d2r   => const_d2r,   &
         tem00 => const_tem00, &
         eps   => const_eps
    use gtool_file, only: &
         fileread
    use scale_external_io, only: &
         externalfileread, &
         externalfilereadoffset
    implicit none

    real(RP), intent(out) :: tg_org(:,:,:)
    real(RP), intent(out) :: strg_org(:,:,:)
    real(RP), intent(out) :: lst_org(:,:)
    real(RP), intent(out) :: llon_org(:,:)
    real(RP), intent(out) :: llat_org(:,:)
    real(RP), intent(out) :: lz_org(:)
    real(RP), intent(out) :: topo_org(:,:)
    real(RP), intent(out) :: lmask_org(:,:)
    character(len=*), intent(in) :: basename_num
    integer,          intent(in) :: ldims(3)
    logical,          intent(in) :: use_file_landwater   ! use land water data from files
    integer,          intent(in) :: it

    ! work
    integer :: k, i, j

    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '+++ ScaleLib/IO[realinput]/Categ[LandInputNICAM]'

    basename = "la_tg"//trim(basename_num)
    call fileread( read1dlz(:), trim(basename), "lev", 1, 1, single=.true. )
    lz_org(:) = read1dlz(:)

    call fileread( read1dx(:), trim(basename), "lon", 1, 1, single=.true. )
    do j = 1, ldims(3)
       llon_org(:,j)  = read1dx(:) * d2r
    enddo

    call fileread( read1dy(:), trim(basename), "lat", 1, 1, single=.true. )
    do i = 1, ldims(2)
       llat_org(i,:)  = read1dy(:) * d2r
    enddo

    basename = "lsmask"//trim(basename_num)
    call externalfileread( read3ds(:,:,:,:), &
                           trim(basename),   &
                           "lsmask",         &
                           it, it, &
                           myrank,           &
                           inicam,           &
                           single=.true.,    &
                           option=.true.     )
    lmask_org(:,:) = real( read3DS(1,:,:,1), kind=rp )

    basename = "la_tg"//trim(basename_num)
    call externalfilereadoffset( read4dl(:,:,:,:), &
                                 trim(basename),   &
                                 "la_tg",          &
                                 it, it,      &
                                 myrank,           &
                                 inicam,           &
                                 single=.true.     )
    tg_org(:,:,:) = real( read4DL(:,:,:,1), kind=rp )

    if( use_file_landwater ) then
       basename = "la_wg"//trim(basename_num)
       call externalfilereadoffset( read4dl(:,:,:,:), &
                                    trim(basename),   &
                                    "la_wg",          &
                                    it, it,      &
                                    myrank,           &
                                    inicam,           &
                                    single=.true.     )
       strg_org(:,:,:) = real( read4DL(:,:,:,1), kind=rp )
    end if

    basename = "ss_tem_sfc"//trim(basename_num)
    call externalfileread( read3ds(:,:,:,:), &
                           trim(basename),   &
                           "ss_tem_sfc",     &
                           it, it,       &
                           myrank,           &
                           inicam,           &
                           single=.true.     )
    lst_org(:,:) = real( read3DS(1,:,:,1), kind=rp )


    ! replace missing value
    do j = 1, ldims(3)
    do i = 1, ldims(2)
    do k = 1, ldims(1)
       if ( abs(tg_org(k,i,j)  -missval_tg  ) < eps ) tg_org(k,i,j)   = undef
       if ( abs(strg_org(k,i,j)-missval_strg) < eps ) strg_org(k,i,j) = undef
    end do
    end do
    end do

    topo_org = undef

    return
  end subroutine parentlandinputnicam

  !-----------------------------------------------------------------------------
  subroutine parentoceansetupnicam( &
       odims,    &
       timelen, &
       basename_org )
    use gtool_file, only: &
         filegetshape
    implicit none

    integer,          intent(out) :: odims(2)
    integer,          intent(out) :: timelen
    character(len=*), intent(in)  :: basename_org

    character(len=H_LONG) :: basename
    integer :: dims_ncm(4)

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

    if( io_l ) write(io_fid_log,*) '+++ Real Case/Ocean Input File Type: NICAM-NETCDF'

    basename = "oa_sst"//trim(basename_org)
    call filegetshape( dims_ncm(:), trim(basename), "oa_sst", 1, single=.true. )
    odims(1) = dims_ncm(1)
    odims(2) = dims_ncm(2)

    timelen = dims_ncm(4)

    if ( .not. allocated(read1dx) ) then
       allocate( read1dx( odims(1)                       ) )
       allocate( read1dy( odims(2)                       ) )
       allocate( read3ds( 1,       odims(1), odims(2), 1 ) )
    end if

    return
  end subroutine parentoceansetupnicam

  !-----------------------------------------------------------------------------
  subroutine parentoceanopennicam( &
       olon_org,     &
       olat_org,     &
       omask_org,    &
       basename_num, &
       odims         )
    use scale_const, only: &
         d2r => const_d2r
    use gtool_file, only: &
         fileread
    use scale_external_io, only: &
         externalfileread
    implicit none

    real(RP),         intent(out) :: olon_org(:,:)
    real(RP),         intent(out) :: olat_org(:,:)
    real(RP),         intent(out) :: omask_org(:,:)
    character(len=*), intent(in)  :: basename_num
    integer,          intent(in)  :: odims(2)

    character(len=H_LONG) :: basename

    integer :: k, i, j

    if( io_l ) write(io_fid_log,*) '+++ ScaleLib/IO[realinput]/Categ[OceanOpenNICAM]'

    basename = "oa_sst"//trim(basename_num)
    call fileread( read1dx(:), trim(basename), "lon", 1, 1, single=.true. )
    do j = 1, odims(2)
       olon_org(:,j)  = read1dx(:) * d2r
    enddo

    call fileread( read1dy(:), trim(basename), "lat", 1, 1, single=.true. )
    do i = 1, odims(1)
       olat_org(i,:)  = read1dy(:) * d2r
    enddo

    basename = "lsmask"//trim(basename_num)
    call externalfileread( read3ds(:,:,:,:), &
                           trim(basename),   &
                           "lsmask",         &
                           1, 1, &
                           myrank,           &
                           inicam,           &
                           single=.true.,    &
                           option=.true.     )
    omask_org(:,:) = real( read3DS(1,:,:,1), kind=rp )

    return
  end subroutine parentoceanopennicam

  !-----------------------------------------------------------------------------
  subroutine parentoceaninputnicam( &
       tw_org,       &
       sst_org,      &
       basename_num, &
       odims,        &
       omask_org,    &
       it            )
    use scale_const, only: &
         undef => const_undef, &
         d2r   => const_d2r,   &
         tem00 => const_tem00, &
         eps   => const_eps
    use gtool_file, only: &
         fileread
    use scale_external_io, only: &
         externalfileread, &
         externalfilereadoffset
    implicit none

    real(RP), intent(out) :: tw_org(:,:)
    real(RP), intent(out) :: sst_org(:,:)
    character(len=*), intent(in) :: basename_num
    integer,          intent(in) :: odims(2)
    real(RP),         intent(in) :: omask_org(:,:)
    integer,          intent(in) :: it

    ! work
    integer :: i, j

    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '+++ ScaleLib/IO[realinput]/Categ[OceanInputNICAM]'


    basename = "oa_sst"//trim(basename_num)
    call externalfilereadoffset( read3ds(:,:,:,:), &
                                 trim(basename),   &
                                 "oa_sst",         &
                                 1, 1,      &
                                 myrank,           &
                                 inicam,           &
                                 single=.true.     )
    sst_org(:,:) = real( read3DS(1,:,:,1), kind=rp )

    !basename = "oa_ice"//trim(basename_num)
    !call ExternalFileRead( read3DS(:,:,:,:), trim(basename), &
    !                       "oa_ice", 1, 1, myrank, iNICAM, single=.true. )
    !ice_org(:,:) = real( read3DS(1,:,:,1), kind=RP )


    ! SST: retrieve SST data around coast
    do j = 1, odims(2)
    do i = 1, odims(1)
       if ( abs(omask_org(i,j)-1.0_rp) < eps ) then ! land data
          cycle
       else
          sst_org(i,j) = ( sst_org(i,j) - missval_sst*omask_org(i,j) ) &
                       / ( 1.0_rp - omask_org(i,j) )
       end if
    end do
    end do

    tw_org = sst_org

    return
  end subroutine parentoceaninputnicam

end module mod_realinput_nicam