Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add 'netcdf' and 'grib2' to input 'source' options in regional IC/LBC routines #96

Merged
merged 10 commits into from
May 7, 2021
33 changes: 22 additions & 11 deletions model/fv_regional_bc.F90
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,7 @@ module fv_regional_mod
integer :: a_step, p_step, k_step, n_step
!
character(len=80) :: data_source
character(len=26), parameter:: data_source_fv3gfs='FV3GFS NEMSIO/NETCDF/GRIB2'
contains

!-----------------------------------------------------------------------
Expand Down Expand Up @@ -1723,7 +1724,7 @@ subroutine regional_bc_data(Atm,bc_hour &
!*** Sensible temperature
!--------------------------
!
if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
if ( trim(data_source)==data_source_fv3gfs ) then
nlev=klev_in
var_name_root='t'
call read_regional_bc_file(is_input,ie_input,js_input,je_input &
Expand Down Expand Up @@ -3676,7 +3677,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &

! Compute true temperature using hydrostatic balance if not read from input.

if (data_source /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
if ( trim(data_source)/=data_source_fv3gfs ) then
do k=1,npz
BC_side%pt_BC(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*BC_side%q_BC(i,j,k,sphum)) )
enddo
Expand All @@ -3696,13 +3697,13 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
! From Jan-Huey Chen's HiRAM code
!-----------------------------------------------------------------------
!
! If the source is FV3GFS GAUSSIAN NEMSIO FILE then all the tracers are in the boundary files
! If the source is FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE then all the tracers are in the boundary files
! and will be read in.
! If the source is from old GFS or operational GSM then the tracers will be fixed in the boundaries
! and may not provide a very good result
!
if (cld_amt .gt. 0) BC_side%q_BC(:,:,:,cld_amt) = 0.
if (trim(data_source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
if ( trim(data_source)/=data_source_fv3gfs ) then
if ( Atm%flagstruct%nwat .eq. 6 ) then
do k=1,npz
do i=is,ie
Expand Down Expand Up @@ -3745,7 +3746,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
enddo
enddo
endif
endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
!
! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
Expand All @@ -3762,7 +3763,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &

call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, Atm%ptop)

if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
if ( trim(data_source)==data_source_fv3gfs ) then
do k=1,npz
do i=is,ie
BC_side%w_BC(i,j,k) = qn1(i,k)
Expand All @@ -3785,7 +3786,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
enddo
enddo

else !<-- datasource /= 'FV3GFS GAUSSIAN NEMSIO FILE'
else !<-- datasource /= 'FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE'
do k=1,npz
do i=is,ie
BC_side%w_BC(i,j,k) = qn1(i,k)/BC_side%delp_BC(i,j,k)*BC_side%delz_BC(i,j,k)
Expand Down Expand Up @@ -6636,21 +6637,31 @@ subroutine get_data_source(source,regional)
!
! This routine extracts the data source information if it is present in the datafile.
!
character (len = 80) :: source
integer :: ncids,sourceLength
character (len=80) :: source
logical :: lstatus,regional
!
! Use the fms call here so we can actually get the return code value.
! The term 'source' is specified by 'chgres_cube'
!
if (regional) then
lstatus = get_global_att_value('INPUT/gfs_data.nc',"source", source)
else
lstatus = get_global_att_value('INPUT/gfs_data.tile1.nc',"source", source)
endif
if (mpp_pe()==0) write(*,*) 'INPUT gfs_data source string=',source
if (.not. lstatus) then
if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute'
source='No Source Attribute'
endif

! data source string to global string --------
if ( trim(source)=='FV3GFS GAUSSIAN NEMSIO FILE' .or. &
trim(source)=='FV3GFS GAUSSIAN NETCDF FILE' .or. &
trim(source)=='FV3GFS GRIB2 FILE' ) then
source = 'FV3GFS NEMSIO/NETCDF/GRIB2'
if (mpp_pe()==0) write(*,*) 'New data source string=',source
endif

end subroutine get_data_source

!---------------------------------------------------------------------
Expand Down Expand Up @@ -6683,7 +6694,7 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
cld_amt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
!
source: if (trim(data_source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
source: if ( trim(data_source)==data_source_fv3gfs ) then
!
! if (cld_amt > 0) BC_side%q_BC(:,:,:,cld_amt) = 0.0 ! Moorthi
do k=1,npz
Expand All @@ -6705,7 +6716,7 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
enddo
enddo
!
else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO FILE
else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
!
! 20160928: Adjust the mixing ratios consistently...
do k=1,npz
Expand Down
22 changes: 11 additions & 11 deletions tools/external_ic.F90
Original file line number Diff line number Diff line change
Expand Up @@ -196,8 +196,8 @@ module external_ic_mod
real, parameter:: zvir = rvgas/rdgas - 1.
real(kind=R_GRID), parameter :: cnst_0p20=0.20d0
real, parameter :: deg2rad = pi/180.
character (len = 80),public :: source ! This tells what the input source was for the data
character(len=27), parameter :: source_fv3gfs = 'FV3GFS GAUSSIAN NEMSIO FILE'
character(len=80), public:: source
character(len=26), parameter:: data_source_fv3gfs='FV3GFS NEMSIO/NETCDF/GRIB2'
public get_external_ic, get_cubed_sphere_terrain
public remap_scalar, remap_dwinds

Expand Down Expand Up @@ -573,8 +573,8 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )

!
call get_data_source(source,Atm%flagstruct%regional)
if (trim(source) == source_fv3gfs) then
call mpp_error(NOTE, "READING FROM REGRIDDED FV3GFS NEMSIO FILE")
if ( trim(source)==data_source_fv3gfs ) then
call mpp_error(NOTE, "READING FROM REGRIDDED FV3GFS NEMSIO/NETCDF/GRIB2 FILE")
endif
!
!--- read in ak and bk from the gfs control file using fms_io read_data ---
Expand Down Expand Up @@ -842,7 +842,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )
snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat')
graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
ntclamt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
if (trim(source) == source_fv3gfs) then
if ( trim(source)==data_source_fv3gfs ) then
do k=1,npz
do j=js,je
do i=is,ie
Expand All @@ -861,7 +861,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos )
enddo
enddo
enddo
else
else ! not NEMSIO/NETCDF/GRIB2
!--- Add cloud condensate from GFS to total MASS
! 20160928: Adjust the mixing ratios consistently...
do k=1,npz
Expand Down Expand Up @@ -3231,7 +3231,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
endif

!$OMP parallel do default(none) &
!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,liq_aero,ice_aero,source, &
!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,liq_aero,ice_aero,source, &
!$OMP cld_amt,ncnst,npz,is,ie,js,je,km,k2,ak0,bk0,psc,t_in,zh,omga,qa,Atm,z500) &
!$OMP private(l,m,pst,pn,gz,pe0,pn0,pe1,pn1,dp2,qp,qn1,gz_fv)
do 5000 j=js,je
Expand Down Expand Up @@ -3375,7 +3375,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
!----------------------------------------------------
! Compute true temperature using hydrostatic balance
!----------------------------------------------------
if (trim(source) /= source_fv3gfs .or. .not. present(t_in)) then
if ( trim(source)/=data_source_fv3gfs .or. .not. present(t_in) ) then
do k=1,npz
#ifdef MULTI_GASES
Atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*virq(Atm%q(i,j,k,:)) )
Expand Down Expand Up @@ -3410,7 +3410,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
! seperate cloud water and cloud ice from Jan-Huey Chen's HiRAM code
! only use for NCEP IC and GFDL microphy
!-----------------------------------------------------------------------
if (trim(source) /= source_fv3gfs) then
if ( trim(source)/=data_source_fv3gfs ) then
if ((Atm%flagstruct%nwat .eq. 3 .or. Atm%flagstruct%nwat .eq. 6) .and. &
(Atm%flagstruct%ncep_ic .or. Atm%flagstruct%nggps_ic)) then
do k=1,npz
Expand Down Expand Up @@ -3459,7 +3459,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
enddo
enddo
endif
endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE

! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
Expand All @@ -3474,7 +3474,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
enddo
enddo
call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, Atm%ptop)
if (trim(source) == source_fv3gfs) then
if ( trim(source)==data_source_fv3gfs ) then
do k=1,npz
do i=is,ie
atm%w(i,j,k) = qn1(i,k)
Expand Down