mod_realinput_nicam.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_realinput_nicam
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_tracer
  use scale_prc, only: &
     myrank => prc_myrank,  &
     prc_abort
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: parentatmossetupnicam
  public :: parentatmosopennicam
  public :: parentatmosinputnicam
  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 :: read3D (:,:,:)
 
  real(SP), allocatable :: read3DL(:,:,:)
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine parentatmossetupnicam( &
      dims,    &
      timelen, &
      basename_org )
    use scale_file, only: &
         file_get_shape
    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)
 
    !---------------------------------------------------------------------------
 
    log_info("ParentAtmosSetupNICAM",*) 'Real Case/Atmos Input File Type: NICAM-NETCDF'
    basename = "ms_pres"//trim(basename_org)
    call file_get_shape( basename, "ms_pres", & ! (in)
                         dims_ncm(:),         & ! (out)
                         single=.true.        ) ! (in)
    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( read3d( dims(1), dims(2), dims(3) ) )
 
    return
  end subroutine parentatmossetupnicam
 
  !-----------------------------------------------------------------------------
  subroutine parentatmosopennicam( &
       lon_org,      &
       lat_org,      &
       cz_org,       &
       basename_num, &
       dims          )
    use scale_const, only: &
         d2r => const_d2r
    use scale_file, only: &
         file_open, &
         file_read
    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 :: fid
    integer :: k, i, j
 
    basename = "ms_pres"//trim(basename_num)
    call file_open( basename, fid, single=.true., postfix="" )
 
    call file_read( fid, "lon", read1dx(:) )
    do j = 1, dims(3)
       lon_org(:,j)  = read1dx(:) * d2r
    enddo
 
    call file_read( fid, "lat", read1dy(:) )
    do i = 1, dims(2)
       lat_org(i,:)  = read1dy(:) * d2r
    enddo
 
    call file_read( fid, "lev", read1dz(:) )
    do j = 1, dims(3)
    do i = 1, dims(2)
       cz_org(3:,i,j) = read1dz(:)
    end do
    end do
 
    call file_read( fid, "ms_pres", read3d(:,:,:) )
    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 ( read3d(i,j,k-2) <= 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
          else
             cz_org(k,i,j) = read3d(i,j,k)
          endif
       enddo
    enddo
    enddo
 
    return
  end subroutine parentatmosopennicam
 
  !-----------------------------------------------------------------------------
  subroutine parentatmosinputnicam( &
       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_file, only: &
       file_read
    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
    !---------------------------------------------------------------------------
 
    basename = "ms_u"//trim(basename_num)
    call file_read( basename, "ms_u", read3d(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = 1, dims(1)
          if ( abs(abs(read3d(i,j,k))-300.0_sp) < sqrt(eps) ) then ! missing value
             velx_org(k+2,i,j) = 0.0_rp
          else
             velx_org(k+2,i,j) = read3d(i,j,k)
          end if
       end do
       velx_org(1:2,i,j) = 0.0_rp
    end do
    end do
 
    basename = "ms_v"//trim(basename_num)
    call file_read( basename, "ms_v", read3d(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = 1, dims(1)
          if ( abs(abs(read3d(i,j,k))-300.0_sp) < sqrt(eps) ) then ! missing value
             vely_org(k+2,i,j) = 0.0_rp
          else
             vely_org(k+2,i,j) = read3d(i,j,k)
          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 FILE_read( basename, "ss_t2m", tsfc_org(:,:), step=it, rankid=myrank, single=.true., postfix="" )
 
    basename = "ms_tem"//trim(basename_num)
    call file_read( basename, "ms_tem", read3d(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read3d(i,j,k) <= 50.0_rp ) then ! missing value
             temp_org(1:k+2,i,j) = temp_org(k+3,i,j)
             exit
          else
             temp_org(k+2,i,j) = read3d(i,j,k)
          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 file_read( basename, "ss_slp", slp_org(:,:), step=it, rankid=myrank, single=.true., postfix="" )
 
    basename = "ms_pres"//trim(basename_num)
    call file_read( basename, "ms_pres", read3d(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read3d(i,j,k) < 0.0_rp ) then ! missing value
             exit
          else
             pres_org(k+2,i,j) = read3d(i,j,k)
          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 file_read( basename, "ss_q2m", qvsfc_org(:,:), step=it, rankid=myrank, single=.true., postfix="" )
 
    qtrc_org(:,:,:,:) = 0.0_rp
    basename = "ms_qv"//trim(basename_num)
    call file_read( basename, "ms_qv", read3d(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, dims(3)
    do i = 1, dims(2)
       do k = dims(1), 1, -1
          if ( read3d(i,j,k) < 0.0_rp ) then ! missing value
             exit
          else
             qtrc_org(k+2,i,j,i_qv) = read3d(i,j,k)
          end if
       end do
       qtrc_org(1:k+2,i,j,i_qv) = qvsfc_org(i,j)
    end do
    end do
 
 
    return
  end subroutine parentatmosinputnicam
 
  !-----------------------------------------------------------------------------
  subroutine parentlandsetupnicam( &
       ldims,    &
       basename_org )
    use scale_file, only: &
         file_get_shape
    implicit none
 
    integer,          intent(out) :: ldims(3)
    character(len=*), intent(in)  :: basename_org
 
    character(len=H_LONG) :: basename
    integer :: dims_ncm(4)
 
    !---------------------------------------------------------------------------
 
    log_info("ParentLandSetupNICAM",*) 'Real Case/Land Input File Type: NICAM-NETCDF'
    basename = "la_tg"//trim(basename_org)
    call file_get_shape( basename, "la_tg", & ! (in)
                         dims_ncm(:),       & ! (out)
                         single=.true.      ) ! (in)
    ! 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)                       ) )
    end if
 
    allocate( read3dl( ldims(2), ldims(3), ldims(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 scale_file, only: &
         file_read
    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
    !---------------------------------------------------------------------------
 
    basename = "la_tg"//trim(basename_num)
    call file_read( basename, "lev", lz_org(:), single=.true., postfix="" )
 
    call file_read( basename, "lon", read1dx(:), single=.true., postfix="" )
    do j = 1, ldims(3)
       llon_org(:,j)  = read1dx(:) * d2r
    enddo
 
    call file_read( basename, "lat", read1dy(:), single=.true., postfix="" )
    do i = 1, ldims(2)
       llat_org(i,:)  = read1dy(:) * d2r
    enddo
 
    basename = "lsmask"//trim(basename_num)
    call file_read( basename, "lsmask", lmask_org(:,:), step=it, rankid=myrank, single=.true., allow_missing=.true., missing_value=undef, postfix="" )
 
    basename = "la_tg"//trim(basename_num)
    call file_read( basename, "la_tg", read3dl(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
    do j = 1, ldims(3)
    do i = 1, ldims(2)
    do k = 1, ldims(1)
       tg_org(k,i,j) = read3dl(i,j,k)
    end do
    end do
    end do
 
    if( use_file_landwater ) then
       basename = "la_wg"//trim(basename_num)
       call file_read( basename, "la_wg", read3dl(:,:,:), step=it, rankid=myrank, single=.true., postfix="" )
       do j = 1, ldims(3)
       do i = 1, ldims(2)
       do k = 1, ldims(1)
          strg_org(k,i,j) = read3dl(i,j,k)
       end do
       end do
       end do
    end if
 
    basename = "ss_tem_sfc"//trim(basename_num)
    call file_read( basename, "ss_tem_sfc", lst_org(:,:), step=it, rankid=myrank, single=.true., postfix="" )
 
    ! 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 scale_file, only: &
         file_get_shape
    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)
 
    !---------------------------------------------------------------------------
 
    log_info("ParentOceanSetupNICAM",*) 'Real Case/Ocean Input File Type: NICAM-NETCDF'
 
    basename = "oa_sst"//trim(basename_org)
    call file_get_shape( basename, "oa_sst", &
                         dims_ncm(:),        &
                         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)                       ) )
    end if
 
    return
  end subroutine parentoceansetupnicam
 
  !-----------------------------------------------------------------------------
  subroutine parentoceanopennicam( &
       olon_org,     &
       olat_org,     &
       omask_org,    &
       basename_num, &
       odims         )
    use scale_const, only: &
         d2r   => const_d2r, &
         undef => const_undef
    use scale_file, only: &
         file_read
    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
 
    basename = "oa_sst"//trim(basename_num)
    call file_read( basename, "lon", read1dx(:), single=.true., postfix="" )
    do j = 1, odims(2)
       olon_org(:,j)  = read1dx(:) * d2r
    enddo
 
    call file_read( basename, "lat", read1dy(:), single=.true., postfix="" )
    do i = 1, odims(1)
       olat_org(i,:)  = read1dy(:) * d2r
    enddo
 
    basename = "lsmask"//trim(basename_num)
    call file_read( basename, "lsmask", omask_org(:,:), step=1, rankid=myrank, single=.true., allow_missing=.true., missing_value=undef, postfix="" )
 
    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 scale_file, only: &
         file_read
    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
    !---------------------------------------------------------------------------
 
    basename = "oa_sst"//trim(basename_num)
    call file_read( basename, "oa_sst", sst_org(:,:), step=1, rankid=myrank, single=.true., postfix="" )
 
    !basename = "oa_ice"//trim(basename_num)
    !call FILE_read( basename, "oa_ice", ice_org(:,:), step=1, rankid=myrank, single=.true., postfix="" )
 
 
    ! 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