module metdata 10,20
!-----------------------------------------------------------------------
!
! BOP
!
! !MODULE: metdata
!
! !DESCRIPTION
! Handles reading and interpolating offline meteorological data which
! is used to drive the dynamics.
!
! !USES
use shr_kind_mod
, only: r8 => shr_kind_r8, r4 => shr_kind_r4
use time_manager, only: get_curr_date, get_step_size, get_calendar
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, begchunk, endchunk
use time_manager, only: get_curr_calday, get_curr_date, get_step_size
use abortutils, only: endrun
use dynamics_vars, only: T_FVDYCORE_GRID
use phys_grid, only: scatter_field_to_chunk
#if ( defined SPMD )
use mpishorthand, only: mpicom, mpir8, mpiint,mpichar
use parutilitiesmodule, only: parcollective, BCSTOP, parscatter
#endif
use perf_mod
use cam_logfile, only: iulog
use pio, only: file_desc_t, pio_put_att, pio_global, pio_get_att, pio_inq_att, pio_inq_dimid, pio_inq_dimlen, &
pio_closefile, pio_get_var, pio_inq_varid
use cam_pio_utils, only: cam_pio_openfile
implicit none
private ! all unless made public
save
! !PUBLIC MEMBERS
public metdata_dyn_init ! subroutine to open files, allocate blocked arrays, etc
public metdata_phys_init ! subroutine to allocate chunked arrays
public advance_met ! subroutine to read more data and interpolate
public get_met_fields ! interface to set meteorology fields
public get_us_vs
public offline_met_defaultopts
public offline_met_setopts
public met_winds_on_walls
public write_met_restart
public read_met_restart
public met_rlx
public fix_mass
interface write_met_restart 1
Module procedure write_met_restart_bin
Module procedure write_met_restart_pio
end interface
interface read_met_restart 1
Module procedure read_met_restart_bin
Module procedure read_met_restart_pio
end interface
!------------------------------------------------------------------
! Interface to access the meteorology fields. Possible invocations
! are as follows:
! call get_met_fields( physics_state, us, vs , tend, dt )
! call get_met_fields( u, v )
! call get_met_fields( srfflx_state )
!------------------------------------------------------------------
Interface get_met_fields ! overload accessors 4
Module Procedure get_dyn_flds
Module Procedure get_uv_centered
Module Procedure get_srf_flds
Module Procedure get_ps
End Interface
real(r8), allocatable :: met_ps_next(:,:) ! PS interpolated to next timestep
real(r8), allocatable :: met_ps_curr(:,:) ! PS interpolated to next timestep
logical :: cell_wall_winds = .false. ! true => met data winds are cell centered
logical :: remove_met_file = .false. ! delete metdata file when finished with it
character(len=16) :: shflx_name = 'SHFLX'
character(len=16) :: qflx_name = 'QFLX'
real(r8) :: shflx_factor = 1._r8
real(r8) :: qflx_factor = 1._r8
! !REVISION HISTORY:
! 31 Oct 2003 Francis Vitt Creation
! 05 Feb 2004 F Vitt Removed reading/inperpolating PS for current timestep
! -- only met_ps_next is needed
! 10 Nov 2004 F Vitt Implemented ability to read from series of files
! 16 Dec 2004 F Vitt Added offline_met_defaultopts and offline_met_setopts
! 14 Jul 2005 W Sawyer Removed pmgrid, spmd_dyn dependencies
! 12 Apr 2006 W Sawyer Removed unneeded ghosting of met_us, met_vs
! 08 Apr 2010 J Edwards Replaced serial netcdf calls with pio interface
!
! EOP
!-----------------------------------------------------------------------
! $Id$
! $Author$
!-----------------------------------------------------------------------
type input2d
real(r8), dimension(:,:), pointer :: data
endtype input2d
type input3d
real(r8), dimension(:,:,:), pointer :: data
endtype input3d
real(r8), allocatable :: met_t(:,:,:) ! interpolated temperature
real(r8), allocatable :: met_u(:,:,:) ! interpolated zonal wind
real(r8), allocatable :: met_v(:,:,:) ! interpolated meridional wind
real(r8), allocatable :: met_us(:,:,:) ! interpolated zonal wind -staggered
real(r8), allocatable :: met_vs(:,:,:) ! interpolated meridional wind -staggered
real(r8), allocatable :: met_q(:,:,:) ! interpolated water vapor
real(r8), allocatable :: met_shflx(:,:)! interpolated surface pressure
real(r8), allocatable :: met_qflx(:,:) ! interpolated water vapor flux
real(r8), allocatable :: met_taux(:,:) ! interpolat
real(r8), allocatable :: met_tauy(:,:) ! interpolated
type(input3d) :: met_ti(2)
type(input3d) :: met_ui(2)
type(input3d) :: met_vi(2)
type(input3d) :: met_usi(2)
type(input3d) :: met_vsi(2)
type(input3d) :: met_qi(2)
type(input2d) :: met_psi_next(2)
type(input2d) :: met_psi_curr(2)
type(input2d) :: met_shflxi(2)
type(input2d) :: met_qflxi(2)
type(input2d) :: met_tauxi(2)
type(input2d) :: met_tauyi(2)
integer :: dateid ! var id of the date in the netCDF
integer :: secid ! var id of the sec data
real(r8) :: datatimem = -1.e36_r8 ! time of prv. values read in
real(r8) :: datatimep = -1.e36_r8 ! time of nxt. values read in
real(r8) :: datatimemn = -1.e36_r8 ! time of prv. values read in for next timestep
real(r8) :: datatimepn = -1.e36_r8 ! time of nxt. values read in for next timestep
integer, parameter :: nm=1 ! array index for previous (minus) data
integer, parameter :: np=2 ! array indes for next (plus) data
real(r8) :: curr_mod_time ! model time - calendar day
real(r8) :: next_mod_time ! model time - calendar day - next time step
character(len=256) :: curr_filename, next_filename, metdata_file
character(len=256) :: filenames_list = ''
type(file_desc_t) :: curr_fileid, next_fileid ! the id of the NetCDF file
real(r8), pointer, dimension(:) :: curr_data_times, next_data_times
real(r8) :: alpha = 1.0_r8 ! don't read in water vapor
! real(r8), private :: alpha = 0.0 ! read in water vaper each time step
real(r8), parameter :: D0_0 = 0.0_r8
real(r8), parameter :: D0_5 = 0.5_r8
real(r8), parameter :: D0_75 = 0.75_r8
real(r8), parameter :: D1_0 = 1.0_r8
real(r8), parameter :: days_per_month = 30.6_r8
real(r8), parameter :: days_per_non_leapyear = 365.0_r8
real(r8), parameter :: days_per_year = 365.25_r8
real(r8), parameter :: seconds_per_day = 86400.0_r8
real(r8), parameter :: fill_value = -9999.0_r8
logical :: online_test = .false.
logical :: debug = .false.
real(r8) :: met_rlx(pver)
integer :: met_levels
real(r8) :: rlx_top_km = 50._r8
real(r8) :: rlx_bot_km = 20._r8
real(r8) :: max_rlx = 1._r8
logical :: fix_mass = .true.
contains
!--------------------------------------------------------------------------
! Get the default runtime options
!--------------------------------------------------------------------------
subroutine offline_met_defaultopts( met_data_file_out, & 1
met_remove_file_out, &
met_cell_wall_winds_out, &
met_filenames_list_out, &
met_rlx_top_km_out, &
met_rlx_bot_km_out, &
met_max_rlx_out, &
met_fix_mass_out, &
met_shflx_name_out, &
met_shflx_factor_out, &
met_qflx_name_out, &
met_qflx_factor_out )
implicit none
character(len=256), intent(out), optional :: met_data_file_out
logical, intent(out), optional :: met_remove_file_out
logical, intent(out), optional :: met_cell_wall_winds_out
character(len=256), intent(out), optional :: met_filenames_list_out
real(r8), intent(out), optional :: met_rlx_top_km_out
real(r8), intent(out), optional :: met_rlx_bot_km_out
real(r8), intent(out), optional :: met_max_rlx_out
logical , intent(out), optional :: met_fix_mass_out
character(len=16), intent(out), optional :: met_shflx_name_out
real(r8), intent(out), optional :: met_shflx_factor_out
character(len=16), intent(out), optional :: met_qflx_name_out
real(r8), intent(out), optional :: met_qflx_factor_out
if ( present( met_data_file_out ) ) then
met_data_file_out = metdata_file
endif
if ( present( met_remove_file_out ) ) then
met_remove_file_out = remove_met_file
endif
if ( present( met_cell_wall_winds_out ) ) then
met_cell_wall_winds_out = cell_wall_winds
endif
if ( present( met_filenames_list_out ) ) then
met_filenames_list_out = filenames_list
endif
if ( present( met_rlx_top_km_out ) ) then
met_rlx_top_km_out = rlx_top_km
endif
if ( present( met_rlx_bot_km_out ) ) then
met_rlx_bot_km_out = rlx_bot_km
endif
if ( present( met_max_rlx_out ) ) then
met_max_rlx_out = max_rlx
endif
if ( present( met_fix_mass_out ) ) then
met_fix_mass_out = fix_mass
endif
if ( present( met_shflx_name_out ) ) then
met_shflx_name_out = shflx_name
endif
if ( present( met_shflx_factor_out ) ) then
met_shflx_factor_out = shflx_factor
endif
if ( present( met_qflx_name_out ) ) then
met_qflx_name_out = qflx_name
endif
if ( present( met_qflx_factor_out ) ) then
met_qflx_factor_out = qflx_factor
endif
end subroutine offline_met_defaultopts
!--------------------------------------------------------------------------
! Set runtime options
!--------------------------------------------------------------------------
subroutine offline_met_setopts( met_data_file_in, & 1
met_remove_file_in, &
met_cell_wall_winds_in, &
met_filenames_list_in, &
met_rlx_top_km_in, &
met_rlx_bot_km_in, &
met_max_rlx_in, &
met_fix_mass_in, &
met_shflx_name_in, &
met_shflx_factor_in, &
met_qflx_name_in, &
met_qflx_factor_in )
implicit none
character(len=256), intent(in), optional :: met_data_file_in
logical, intent(in), optional :: met_remove_file_in
logical, intent(in), optional :: met_cell_wall_winds_in
character(len=256), intent(in), optional :: met_filenames_list_in
real(r8), intent(in), optional :: met_rlx_top_km_in
real(r8), intent(in), optional :: met_rlx_bot_km_in
real(r8), intent(in), optional :: met_max_rlx_in
logical , intent(in), optional :: met_fix_mass_in
character(len=16), intent(in), optional :: met_shflx_name_in
real(r8), intent(in), optional :: met_shflx_factor_in
character(len=16), intent(in), optional :: met_qflx_name_in
real(r8), intent(in), optional :: met_qflx_factor_in
if ( present( met_data_file_in ) ) then
metdata_file = met_data_file_in
endif
if ( present( met_remove_file_in ) ) then
remove_met_file = met_remove_file_in
endif
if ( present( met_cell_wall_winds_in ) ) then
cell_wall_winds = met_cell_wall_winds_in
endif
if ( present( met_filenames_list_in ) ) then
filenames_list = met_filenames_list_in
endif
if ( present( met_rlx_top_km_in ) ) then
rlx_top_km = met_rlx_top_km_in
endif
if ( present( met_rlx_bot_km_in ) ) then
rlx_bot_km = met_rlx_bot_km_in
endif
if ( present( met_max_rlx_in ) ) then
max_rlx = met_max_rlx_in
endif
if ( present( met_fix_mass_in ) ) then
fix_mass = met_fix_mass_in
endif
if ( present( met_shflx_name_in ) ) then
shflx_name = met_shflx_name_in
endif
if ( present( met_shflx_factor_in ) ) then
shflx_factor = met_shflx_factor_in
endif
if ( present( met_qflx_name_in ) ) then
qflx_name = met_qflx_name_in
endif
if ( present( met_qflx_factor_in ) ) then
qflx_factor = met_qflx_factor_in
endif
if (masterproc) then
write(iulog,*)'Time-variant meteorological dataset (metdata_file) is: ', trim(metdata_file)
write(iulog,*)'Meteorological data file will be removed (remove_met_file): ', remove_met_file
write(iulog,*)'Meteorological winds are on cell walls (cell_wall_winds): ', cell_wall_winds
write(iulog,*)'Meteorological file names list file: ', trim(filenames_list)
write(iulog,*)'Meteorological relaxation top is (km): ', rlx_top_km
write(iulog,*)'Meteorological relaxation bottom is (km): ', rlx_bot_km
write(iulog,*)'Meteorological maximum relaxation is : ',max_rlx
write(iulog,*)'Offline driver mass fixer is trurned on (fix_mass): ',fix_mass
write(iulog,*)'Meteorological shflx field name : ', trim(shflx_name)
write(iulog,*)'Meteorological shflx multiplication factor : ', shflx_factor
write(iulog,*)'Meteorological qflx field name : ', trim(qflx_name)
write(iulog,*)'Meteorological qflx multiplication factor : ', qflx_factor
endif
end subroutine offline_met_setopts
!--------------------------------------------------------------------------
! Opens file, allocates arrays
!--------------------------------------------------------------------------
subroutine metdata_dyn_init(grid) 1,4
use infnan, only : nan
use cam_control_mod, only : nsrest
implicit none
! !INPUT PARAMETERS:
type (T_FVDYCORE_GRID), intent(in) :: grid
character(len=256) :: filen
integer :: im, jm, km, jfirst, jlast, kfirst, klast
integer :: ng_d, ng_s
integer :: ierr
im = grid%im
km = grid%km
jfirst = grid%jfirst
jlast = grid%jlast
kfirst = grid%kfirst
klast = grid%klast
ng_d = grid%ng_d
ng_s = grid%ng_s
if (nsrest/=1) then ! initial run or branch run
curr_filename = metdata_file
next_filename = ''
else
! restart run
! curr_filename & next_filename already set by restart_dynamics
endif
call open_met_datafile( grid, curr_filename, curr_fileid, curr_data_times, check_dims=.true. )
if ( len_trim(next_filename) > 0 ) &
call open_met_datafile( grid, next_filename, next_fileid, next_data_times )
!
! allocate space for data arrays ...
!
! dynamics grid
allocate( met_psi_next(nm)%data(im, jfirst:jlast) )
allocate( met_psi_next(np)%data(im, jfirst:jlast) )
allocate( met_psi_curr(nm)%data(im, jfirst:jlast) )
allocate( met_psi_curr(np)%data(im, jfirst:jlast) )
allocate( met_ps_next(im, jfirst:jlast) )
allocate( met_ps_curr(im, jfirst:jlast) )
allocate( met_us(im, jfirst-ng_d:jlast+ng_s, kfirst:klast) )
allocate( met_vs(im, jfirst-ng_s:jlast+ng_d, kfirst:klast) )
met_us = nan
met_vs = nan
if (cell_wall_winds) then
allocate( met_usi(nm)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_usi(np)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_vsi(nm)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_vsi(np)%data(im, jfirst:jlast, kfirst:klast) )
endif
if (.not. cell_wall_winds) then
allocate( met_u(im, jfirst-ng_d:jlast+ng_d, kfirst:klast) )
allocate( met_ui(nm)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_ui(np)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_v(im, jfirst-ng_s:jlast+ng_d, kfirst:klast) )
allocate( met_vi(nm)%data(im, jfirst:jlast, kfirst:klast) )
allocate( met_vi(np)%data(im, jfirst:jlast, kfirst:klast) )
endif
end subroutine metdata_dyn_init
!=================================================================================
subroutine metdata_phys_init 1,11
use cam_history, only : addfld,dyn_decomp,phys_decomp
call addfld ('MET_RLX ',' ',pver, 'A','Meteorology relax function',dyn_decomp)
call addfld ('MET_TAUX',' ',1, 'A','Meteorology taux',phys_decomp)
call addfld ('MET_TAUY',' ',1, 'A','Meteorology tauy',phys_decomp)
call addfld ('MET_SHFX',' ',1, 'A','Meteorology shflx',phys_decomp)
call addfld ('MET_QFLX',' ',1, 'A','Meteorology qflx',phys_decomp)
call addfld ('MET_PS',' ',1, 'A','Meteorology PS',dyn_decomp)
call addfld ('MET_T',' ',pver, 'A','Meteorology T',phys_decomp)
call addfld ('MET_U',' ',pver, 'A','Meteorology U',dyn_decomp)
call addfld ('MET_V',' ',pver, 'A','Meteorology V',dyn_decomp)
! allocate chunked arrays
allocate( met_ti(nm)%data(pcols,pver,begchunk:endchunk) )
allocate( met_ti(np)%data(pcols,pver,begchunk:endchunk) )
allocate( met_t(pcols,pver,begchunk:endchunk) )
allocate( met_qi(nm)%data(pcols,pver,begchunk:endchunk) )
allocate( met_qi(np)%data(pcols,pver,begchunk:endchunk) )
allocate( met_q(pcols,pver,begchunk:endchunk) )
allocate( met_shflxi(nm)%data(pcols,begchunk:endchunk) )
allocate( met_shflxi(np)%data(pcols,begchunk:endchunk) )
allocate( met_shflx(pcols,begchunk:endchunk) )
allocate( met_qflxi(nm)%data(pcols,begchunk:endchunk) )
allocate( met_qflxi(np)%data(pcols,begchunk:endchunk) )
allocate( met_qflx(pcols,begchunk:endchunk) )
allocate( met_tauxi(nm)%data(pcols,begchunk:endchunk) )
allocate( met_tauxi(np)%data(pcols,begchunk:endchunk) )
allocate( met_taux(pcols,begchunk:endchunk) )
allocate( met_tauyi(nm)%data(pcols,begchunk:endchunk) )
allocate( met_tauyi(np)%data(pcols,begchunk:endchunk) )
allocate( met_tauy(pcols,begchunk:endchunk) )
call set_met_rlx()
end subroutine metdata_phys_init
!-----------------------------------------------------------------------
! Reads more data if needed and interpolates data to current model time
!-----------------------------------------------------------------------
subroutine advance_met(grid) 1,7
use cam_history, only : outfld
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
real(r8) :: met_rlx_2d(grid%ifirstxy:grid%ilastxy,grid%km)
integer :: i,j,k, idim
call t_startf('MET__advance')
!
!
call get_model_time()
if ( ( curr_mod_time > datatimep ) .or. &
( next_mod_time > datatimepn ) ) then
call check_files(grid)
endif
if ( curr_mod_time > datatimep ) then
call read_next_metdata(grid)
end if
if ( next_mod_time > datatimepn ) then
call read_next_ps(grid)
end if
! need to inperpolate the data, regardless !
! each mpi tasks needs to interpolate
call interpolate_metdata(grid)
call t_stopf('MET__advance')
idim = grid%ilastxy - grid%ifirstxy + 1
do j = grid%jfirstxy, grid%jlastxy
do k = 1, grid%km
do i = grid%ifirstxy, grid%ilastxy
met_rlx_2d(i,k) = met_rlx(k)
enddo
enddo
call outfld('MET_RLX',met_rlx_2d, idim, j)
enddo
end subroutine advance_met
!-------------------------------------------------------------------
! Method to get some the meteorology data.
! Sets the following srfflx_state member fields to the
! meteorology data :
! qflx
! shflx
! taux
! tauy
!-------------------------------------------------------------------
subroutine get_srf_flds( srfflx_state2d ) 1,8
use camsrfexch_types, only: srfflx_state
use phys_grid, only: get_ncols_p
use cam_history, only: outfld
implicit none
type(srfflx_state), intent(inout), dimension(begchunk:endchunk) :: srfflx_state2d
integer :: c,ncol
call t_startf('MET__GET_SRF_FLDS')
do c=begchunk,endchunk
ncol = get_ncols_p(c)
srfflx_state2d(c)%wsx(:ncol) = met_taux(:ncol,c)
srfflx_state2d(c)%wsy(:ncol) = met_tauy(:ncol,c)
srfflx_state2d(c)%shf(:ncol) = met_shflx(:ncol,c)* shflx_factor
srfflx_state2d(c)%cflx(:ncol,1) = met_qflx(:ncol,c) * qflx_factor
end do ! Chunk loop
if (debug) then
write(iulog,*)'METDATA maxval(met_taux),minval(met_taux): ',maxval(met_taux),minval(met_taux)
write(iulog,*)'METDATA maxval(met_tauy),minval(met_tauy): ',maxval(met_tauy),minval(met_tauy)
write(iulog,*)'METDATA maxval(met_shflx),minval(met_shflx): ',maxval(met_shflx),minval(met_shflx)
write(iulog,*)'METDATA maxval(met_qflx),minval(met_qflx): ',maxval(met_qflx),minval(met_qflx)
endif
if ( debug ) then
do c = begchunk, endchunk
call outfld('MET_TAUX',srfflx_state2d(c)%wsx , pcols ,c )
call outfld('MET_TAUY',srfflx_state2d(c)%wsy , pcols ,c )
call outfld('MET_SHFX',srfflx_state2d(c)%shf , pcols ,c )
call outfld('MET_QFLX',srfflx_state2d(c)%cflx(:,1) , pcols ,c )
enddo
endif
call t_stopf('MET__GET_SRF_FLDS')
end subroutine get_srf_flds
!-------------------------------------------------------------------
! allows access to physics state fields
! q : water vapor
! ps : surface pressure
! t : temperature
!-------------------------------------------------------------------
subroutine get_dyn_flds( state, tend, dt ) 1,7
use physics_types, only: physics_state, physics_tend, physics_dme_adjust
use ppgrid, only: pcols, pver, begchunk, endchunk
use phys_grid, only: get_ncols_p
use cam_history, only: outfld
implicit none
type(physics_state), intent(inout), dimension(begchunk:endchunk) :: state
type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: tend
real(r8), intent(in ) :: dt ! model physics timestep
integer :: lats(pcols) ! array of latitude indices
integer :: lons(pcols) ! array of longitude indices
integer :: c, ncol, i,j,k
real(r8):: qini(pcols,pver) ! initial specific humidity
real(r8) :: tmp(pcols,pver)
call t_startf('MET__GET_DYN2')
do c = begchunk, endchunk
ncol = get_ncols_p(c)
do k=1,pver
do i=1,ncol
state(c)%t(i,k) = (1._r8-met_rlx(k))*state(c)%t(i,k) + met_rlx(k)*met_t(i,k,c)
qini(i,k) = state(c)%q(i,k,1)
! at this point tracer mixing ratios have already been
! converted from dry to moist
!!$ if ( moist_q_mmr .and. (.not. online_test)) then
state(c)%q(i,k,1) = alpha*state(c)%q(i,k,1) + &
(1-alpha)*met_q(i,k,c)
!!$ else
!!$ ! dry-to-moist conversion
!!$ state(c)%q(i,k,1) = alpha*state(c)%q(i,k,1) + &
!!$ (1.-alpha)*met_q(i,c,k) &
!!$ * state(c)%pdeldry(i,k)/state(c)%pdel(i,k)
!!$ endif
if ((state(c)%q(i,k,1) < D0_0).and. (alpha .ne. D1_0 )) state(c)%q(i,k,1) = D0_0
end do
end do
! now adjust mass of each layer now that water vapor has changed
if (( .not. online_test ) .and. (alpha .ne. D1_0 )) then
call physics_dme_adjust(state(c), tend(c), qini, dt)
endif
end do
if (debug) then
write(iulog,*)'METDATA maxval(met_t),minval(met_t): ', maxval(met_t),minval(met_t)
write(iulog,*)'METDATA maxval(met_ps_next),minval(met_ps_next): ', maxval(met_ps_next),minval(met_ps_next)
endif
do c = begchunk, endchunk
call outfld('MET_T ',state(c)%t , pcols ,c )
enddo
call t_stopf('MET__GET_DYN2')
end subroutine get_dyn_flds
!------------------------------------------------------------------------
! get the meteorological winds on the grid cell centers (A-grid)
!------------------------------------------------------------------------
subroutine get_uv_centered( grid, u, v ) 1,3
use cam_history, only: outfld
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
real(r8), intent(out) :: u(grid%im, grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d, &
grid%kfirst:grid%klast) ! u-wind on A-grid
real(r8), intent(out) :: v(grid%im, grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d, &
grid%kfirst:grid%klast) ! v-wind on A-grid
integer :: i,j,k
integer :: jm, jfirst, jlast, jfirstxy, jlastxy, kfirst, klast, ng_d, ng_s, ifirstxy, ilastxy
real(r8) :: u3s_tmp(grid%ifirstxy:grid%ilastxy,grid%km)
real(r8) :: v3s_tmp(grid%ifirstxy:grid%ilastxy,grid%km)
jm = grid%jm
jfirstxy= grid%jfirstxy
jlastxy = grid%jlastxy
jfirst = grid%jfirst
jlast = grid%jlast
kfirst = grid%kfirst
klast = grid%klast
ifirstxy= grid%ifirstxy
ilastxy = grid%ilastxy
ng_d = grid%ng_d
ng_s = grid%ng_s
u(:,:,:) = D0_0
v(:,:,:) = D0_0
u( :, max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast ) = &
met_u(:, max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )
v( :, max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast ) = &
met_v(:, max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast )
if (debug) write(iulog,*)'METDATA maxval(u),minval(u),maxval(v),minval(v) : ',&
maxval(u(:, max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )),&
minval(u(:, max(1,jfirst-ng_d):min(jm,jlast+ng_d), kfirst:klast )),&
maxval(v(:, max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast )),&
minval(v(:, max(1,jfirst-ng_s):min(jm,jlast+ng_d), kfirst:klast ))
if ( grid%twod_decomp .eq. 0 ) then
do j = jfirst, jlast
do k = kfirst, klast
do i = 1, grid%im
u3s_tmp(i,k) = u(i,j,k)
v3s_tmp(i,k) = v(i,j,k)
enddo
enddo
call outfld ('MET_U ', u3s_tmp, grid%im, j )
call outfld ('MET_V ', v3s_tmp, grid%im, j )
enddo
endif
end subroutine get_uv_centered
!------------------------------------------------------------------------
! get the meteorological surface pressure interp to dyn substep
!------------------------------------------------------------------------
subroutine get_ps( grid, ps, nsubsteps, n ) 1,2
use cam_history, only: outfld
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
real(r8), intent(out) :: ps(grid%im, grid%jfirst:grid%jlast)
integer, intent(in) :: nsubsteps
integer, intent(in) :: n
real(r8) :: num1, num2
integer :: j
num1 = n
num2 = nsubsteps
ps(:,:) = met_ps_curr(:,:) + num1*(met_ps_next(:,:)-met_ps_curr(:,:))/num2
if ( grid%twod_decomp .eq. 0 ) then
do j = grid%jfirst, grid%jlast
call outfld('MET_PS',ps(:,j), grid%im ,j )
enddo
endif
end subroutine get_ps
!------------------------------------------------------------------------
! get the meteorological winds on the grid cell walls (vorticity winds)
! us : staggered zonal wind
! vs : staggered meridional wind
!------------------------------------------------------------------------
subroutine get_us_vs( grid, us, vs ) 1
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
real(r8), intent(inout) :: us(grid%im, grid%jfirst-grid%ng_d:grid%jlast+grid%ng_s, &
grid%kfirst:grid%klast) ! u-wind on d-grid
real(r8), intent(inout) :: vs(grid%im, grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d, &
grid%kfirst:grid%klast) ! v-wind on d-grid
integer :: i,j,k
integer :: jm, jfirst, jlast, kfirst, klast, ng_d, ng_s
jm = grid%jm
jfirst = grid%jfirst
jlast = grid%jlast
kfirst = grid%kfirst
klast = grid%klast
ng_d = grid%ng_d
ng_s = grid%ng_s
call t_startf('MET__get_us_vs')
! vertical relaxation (blending) occurs in dyn_run (dyn_comp.F90)
us(:,:,:) = 1.e36_r8
vs(:,:,:) = 1.e36_r8
us( :, max(2,jfirst): min(jm,jlast), kfirst:klast) = &
met_us( :, max(2,jfirst): min(jm,jlast), kfirst:klast)
vs( :, max(1,jfirst): min(jm,jlast), kfirst:klast) = &
met_vs( :, max(1,jfirst): min(jm,jlast), kfirst:klast)
if (debug) write(iulog,*)grid%iam,': METDATA maxval(us),minval(us),maxval(vs),minval(vs) : ',&
maxval(us( :, max(2,jfirst): min(jm,jlast), kfirst:klast)),&
minval(us( :, max(2,jfirst): min(jm,jlast), kfirst:klast)),&
maxval(vs( :, max(1,jfirst): min(jm,jlast), kfirst:klast)),&
minval(vs( :, max(1,jfirst): min(jm,jlast), kfirst:klast))
!!$ if (debug) then
!!$ u3s_tmp = 1.e36
!!$ do j = jfirst, jlast
!!$ do k = kfirst, klast
!!$ do i = 1, im
!!$ if (j >= 2) u3s_tmp(i,k) = us(i,j,k)
!!$ v3s_tmp(i,k) = vs(i,j,k)
!!$ enddo
!!$ enddo
!!$ call outfld ('MET_US ', u3s_tmp, im, j )
!!$ call outfld ('MET_VS ', v3s_tmp, im, j )
!!$ enddo
!!$ endif
!!$
call t_stopf('MET__get_us_vs')
end subroutine get_us_vs
!-------------------------------------------------------------------------
! writes file names to restart file
!-------------------------------------------------------------------------
subroutine write_met_restart_pio(File) 1
type(file_desc_t), intent(inout) :: File
integer :: ierr
ierr = pio_put_att(File, PIO_GLOBAL, 'current_metdata_filename', curr_filename)
ierr = pio_put_att(File, PIO_GLOBAL, 'next_metdata_filename', next_filename)
end subroutine write_met_restart_pio
subroutine read_met_restart_pio(File) 1
type(file_desc_t), intent(inout) :: File
integer :: ierr, xtype, slen
ierr = pio_inq_att(File, PIO_GLOBAL, 'current_metdata_filename',xtype, slen)
ierr = pio_get_att(File, PIO_GLOBAL, 'current_metdata_filename', curr_filename)
curr_filename(slen+1:256) = ''
ierr = pio_inq_att(File, PIO_GLOBAL, 'next_metdata_filename',xtype, slen)
ierr = pio_get_att(File, PIO_GLOBAL, 'next_metdata_filename', next_filename)
next_filename(slen+1:256) = ''
end subroutine read_met_restart_pio
subroutine write_met_restart_bin( nrg ) 1,2
implicit none
integer,intent(in) :: nrg ! Unit number
integer :: ioerr ! error status
if (masterproc) then
write(nrg, iostat=ioerr) curr_filename
if (ioerr /= 0 ) then
write(iulog,*) 'WRITE ioerror ',ioerr,' on i/o unit = ',nrg
call endrun ('WRITE_RESTART_DYNAMICS')
end if
write(nrg, iostat=ioerr) next_filename
if (ioerr /= 0 ) then
write(iulog,*) 'WRITE ioerror ',ioerr,' on i/o unit = ',nrg
call endrun ('WRITE_RESTART_DYNAMICS')
end if
end if
end subroutine write_met_restart_bin
!-------------------------------------------------------------------------
! reads file names from restart file
!-------------------------------------------------------------------------
subroutine read_met_restart_bin( nrg ) 1,4
implicit none
integer,intent(in) :: nrg ! Unit number
integer :: ioerr ! error status
if (masterproc) then
read(nrg, iostat=ioerr) curr_filename
if (ioerr /= 0 ) then
write(iulog,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
call endrun ('READ_RESTART_DYNAMICS')
end if
read(nrg, iostat=ioerr) next_filename
if (ioerr /= 0 ) then
write(iulog,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
call endrun ('READ_RESTART_DYNAMICS')
end if
end if
#if ( defined SPMD )
call mpibcast ( curr_filename ,len(curr_filename) ,mpichar,0,mpicom)
call mpibcast ( next_filename ,len(next_filename) ,mpichar,0,mpicom)
#endif
end subroutine read_met_restart_bin
!-------------------------------------------------------------------------
! returns true if the met winds are defined on cell walls
!-------------------------------------------------------------------------
function met_winds_on_walls()
logical :: met_winds_on_walls
met_winds_on_walls = cell_wall_winds
end function met_winds_on_walls
! internal methods :
!-------------------------------------------------------------------------
! transfers cell-centered winds to cell walls
!-------------------------------------------------------------------------
subroutine transfer_windsToWalls(grid) 1
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
integer :: i,j,k
integer :: im, jfirst, jlast, kfirst, klast
im = grid%im
jfirst = grid%jfirst
jlast = grid%jlast
kfirst = grid%kfirst
klast = grid%klast
call t_startf('MET__transfer_windsToWalls')
!$omp parallel do private (i, j, k)
do k = kfirst, klast
do j = jfirst+1,jlast
do i = 1,im
met_us(i,j,k) = ( met_u(i,j,k) + met_u(i,j-1,k) )*D0_5
end do
end do
#if defined( SPMD )
if ( jfirst .gt. 1 ) then
do i = 1, im
! met_u is alread ghosted at this point
met_us(i,jfirst,k) = ( met_u(i,jfirst,k) + met_u(i,jfirst-1,k) )*D0_5
enddo
endif
#endif
do j = jfirst,jlast
met_vs(1,j,k) = ( met_v(1,j,k) + met_v(im,j,k) )*D0_5
do i = 2,im
met_vs(i,j,k) = ( met_v(i,j,k) + met_v(i-1,j,k) )*D0_5
end do
end do
end do
call t_stopf('MET__transfer_windsToWalls')
end subroutine transfer_windsToWalls
subroutine get_model_time() 5,9
implicit none
integer yr, mon, day, ncsec ! components of a date
call t_startf('MET__get_model_time')
call get_curr_date(yr, mon, day, ncsec)
curr_mod_time = get_time_float( yr, mon, day, ncsec )
next_mod_time = curr_mod_time + get_step_size()/seconds_per_day
call t_stopf('MET__get_model_time')
end subroutine get_model_time
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine check_files(grid) 2,15
use shr_sys_mod, only: shr_sys_system
use ioFileMod, only: getfil
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
!-----------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------
character(len=256) :: ctmp
character(len=256) :: loc_fname
integer :: istat
if (next_mod_time > curr_data_times(size(curr_data_times))) then
if ( .not. associated(next_data_times) ) then
! open next file...
next_filename = incr_filename( curr_filename )
call open_met_datafile( grid, next_filename, next_fileid, next_data_times )
endif
endif
if ( associated(next_data_times) ) then
if (curr_mod_time >= next_data_times(1)) then
! close current file ...
if (grid%iam == 0) then
call pio_closefile( curr_fileid )
! remove if requested
if( remove_met_file ) then
call getfil( curr_filename, loc_fname, 0 )
write(iulog,*) 'check_files: removing file = ',trim(loc_fname)
ctmp = 'rm -f ' // trim(loc_fname)
write(iulog,*) 'check_files: fsystem issuing command - '
write(iulog,*) trim(ctmp)
call shr_sys_system( ctmp, istat )
end if
endif
curr_filename = next_filename
curr_fileid = next_fileid
deallocate( curr_data_times )
allocate( curr_data_times( size( next_data_times ) ) )
curr_data_times(:) = next_data_times(:)
next_filename = ''
deallocate( next_data_times )
nullify( next_data_times )
endif
endif
end subroutine check_files
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
function incr_filename( filename ) 2,20
!-----------------------------------------------------------------------
! ... Increment or decrement a date string withing a filename
! the filename date section is assumed to be of the form
! yyyy-dd-mm
!-----------------------------------------------------------------------
use string_utils, only : incstr
use shr_file_mod, only : shr_file_getunit, shr_file_freeunit
implicit none
character(len=*), intent(in) :: filename ! present dynamical dataset filename
character(len=256) :: incr_filename ! next filename in the sequence
! set new next_filename ...
!-----------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------
integer :: pos, pos1, istat
character(len=256) :: fn_new, line
character(len=6) :: seconds
character(len=5) :: num
integer :: ios,unitnumber
if ( len_trim(filenames_list) .eq. 0) then
!-----------------------------------------------------------------------
! ... ccm type filename
!-----------------------------------------------------------------------
pos = len_trim( filename )
fn_new = filename(:pos)
write(iulog,*) 'incr_flnm: old filename = ',trim(fn_new)
if( fn_new(pos-2:) == '.nc' ) then
pos = pos - 3
end if
istat = incstr( fn_new(:pos), 1 )
if( istat /= 0 ) then
write(iulog,*) 'incr_flnm: incstr returned ', istat
write(iulog,*) ' while trying to decrement ',trim( fn_new )
call endrun
end if
else
! open filenames_list
write(iulog,*) 'incr_flnm: old filename = ',trim(filename)
write(iulog,*) 'incr_flnm: open filenames_list : ',filenames_list
unitnumber = shr_file_getUnit()
open( unit=unitnumber, file=filenames_list, iostat=ios, status="OLD")
if (ios /= 0) then
call endrun('not able to open filenames_list file: '//filenames_list)
endif
! read file names
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//filenames_list)
endif
do while( trim(line) /= trim(filename) )
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//filenames_list)
endif
enddo
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//filenames_list)
endif
fn_new = trim(line)
close(unit=unitnumber)
call shr_file_freeUnit(unitnumber)
endif
incr_filename = trim(fn_new)
write(iulog,*) 'incr_flnm: new filename = ',incr_filename
end function incr_filename
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine find_times( itms, fids, datatm, datatp, time ) 5,5
implicit none
integer, intent(out) :: itms(2) ! record numbers that bracket time
type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs
real(r8), intent(in) :: time ! time of interest
real(r8), intent(out):: datatm, datatp
integer np1 ! current forward time index of dataset
integer n,i !
integer :: curr_tsize, next_tsize, all_tsize
real(r8), allocatable, dimension(:):: all_data_times
curr_tsize = size(curr_data_times)
next_tsize = 0
if ( associated(next_data_times)) next_tsize = size(next_data_times)
all_tsize = curr_tsize + next_tsize
allocate( all_data_times( all_tsize ) )
all_data_times(:curr_tsize) = curr_data_times(:)
if (next_tsize > 0) all_data_times(curr_tsize+1:all_tsize) = next_data_times(:)
! find bracketing times
do n=1, all_tsize-1
np1 = n + 1
datatm = all_data_times(n)
datatp = all_data_times(np1)
if ( (time .ge. datatm) .and. (time .le. datatp) ) then
goto 20
endif
enddo
write(iulog,*)'FIND_TIMES: Failed to find dates bracketing desired time =', time
write(iulog,*)' datatm = ',datatm
write(iulog,*)' datatp = ',datatp
write(iulog,*)' all_data_times = ',all_data_times
call endrun
20 continue
deallocate( all_data_times )
itms(1) = n
itms(2) = np1
fids(:) = curr_fileid
do i=1,2
if ( itms(i) > curr_tsize ) then
itms(i) = itms(i) - curr_tsize
fids(i) = next_fileid
endif
enddo
end subroutine find_times
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine read_next_ps(grid) 1,4
use ncdio_atm, only: infld
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
integer :: recnos(2)
type(file_desc_t) :: fids(2)
character(len=8) :: varname
integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
logical :: readvar
if(online_test) then
varname='arch_PS'
else
varname='PS'
end if
jfirstxy= grid%jfirstxy
jlastxy = grid%jlastxy
ifirstxy= grid%ifirstxy
ilastxy = grid%ilastxy
call find_times( recnos, fids, datatimemn, datatimepn, next_mod_time )
call infld(varname, fids(1), 'lon', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
met_psi_next(nm)%data, readvar, grid_map='DYN', timelevel=recnos(1))
call infld(varname, fids(2), 'lon', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
met_psi_next(np)%data, readvar, grid_map='DYN', timelevel=recnos(2))
if(masterproc) write(iulog,*)'READ_NEXT_PS: Read meteorological data '
end subroutine read_next_ps
!------------------------------------------------------------------------
!------------------------------------------------------------------------
subroutine read_next_metdata(grid) 1,13
use ncdio_atm, only: infld
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
integer recnos(2), i !
type(file_desc_t) :: fids(2)
character(len=8) :: Uname, Vname, Tname, Qname, psname
integer :: im, jm, km
logical :: readvar
integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
call t_startf('MET__read_next_metdata')
jfirstxy= grid%jfirstxy
jlastxy = grid%jlastxy
ifirstxy= grid%ifirstxy
ilastxy = grid%ilastxy
im = grid%im
jm = grid%jm
km = grid%km
call find_times( recnos, fids, datatimem, datatimep, curr_mod_time )
!
! Set up hyperslab corners
!
if(online_test) then
Tname='arch_T'
Qname='arch_Q'
PSname='arch_PS'
if(cell_wall_winds) then
Uname='arch_US'
Vname='arch_VS'
else
Uname='arch_U'
Vname='arch_V'
end if
else
Tname='T'
Qname='Q'
PSname='PS'
if(cell_wall_winds) then
Uname='US'
Vname='VS'
else
Uname='U'
Vname='V'
end if
end if
do i=1,2
call infld(Tname, fids(i), 'lon', 'lev', 'lat', 1, pcols, 1, pver, &
begchunk, endchunk, met_ti(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
if (cell_wall_winds) then
call infld(Uname, fids(i), 'lon', 'lev', 'slat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
1,km, met_usi(i)%data, readvar, grid_map='DYN',timelevel=recnos(i))
call infld(Vname, fids(i), 'slon', 'lev', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
1,km,met_vsi(i)%data, readvar, grid_map='DYN',timelevel=recnos(i))
else
call infld(Uname, fids(i), 'lon', 'lev', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
1,km,met_ui(i)%data, readvar, grid_map='DYN',timelevel=recnos(i))
call infld(Vname, fids(i), 'lon', 'lev', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
1,km,met_vi(i)%data, readvar, grid_map='DYN',timelevel=recnos(i))
endif
call infld(Qname, fids(i), 'lon', 'lev', 'lat', 1, pcols, 1, pver, &
begchunk, endchunk, met_qi(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
call infld(shflx_name, fids(i), 'lon', 'lat', 1, pcols, begchunk, endchunk, &
met_shflxi(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
call infld(Qflx_name, fids(i), 'lon', 'lat', 1, pcols, begchunk, endchunk, &
met_qflxi(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
call infld('TAUX', fids(i), 'lon', 'lat', 1, pcols, begchunk, endchunk, &
met_tauxi(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
call infld('TAUY', fids(i), 'lon', 'lat', 1, pcols, begchunk, endchunk, &
met_tauyi(i)%data, readvar, grid_map='PHYS',timelevel=recnos(i))
call infld(PSname, fids(i), 'lon', 'lat', ifirstxy, ilastxy, jfirstxy, jlastxy, &
met_psi_curr(i)%data, readvar, grid_map='DYN',timelevel=recnos(i))
enddo
if(masterproc) write(iulog,*)'READ_NEXT_METDATA: Read meteorological data '
call t_stopf('MET__read_next_metdata')
end subroutine read_next_metdata
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine interpolate_metdata(grid) 1,6
#if defined( SPMD )
use mod_comm, only : mp_send4d_ns, mp_recv4d_ns
#endif
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
real(r4) fact1, fact2
real(r4) nfact1, nfact2
real(r8) deltat,deltatn
integer :: i,j,k
integer :: im, jm, km, jfirst, jlast, kfirst, klast, ng_d, ng_s
im = grid%im
jm = grid%jm
km = grid%km
jfirst = grid%jfirst
jlast = grid%jlast
kfirst = grid%kfirst
klast = grid%klast
ng_d = grid%ng_d
ng_s = grid%ng_s
call t_startf('MET__interpolate_metdata')
deltat = datatimep - datatimem
deltatn = datatimepn - datatimemn
fact1 = (datatimep - curr_mod_time)/deltat
! fact2 = (curr_mod_time - datatimem)/deltat
fact2 = D1_0-fact1
nfact1 = (datatimepn - next_mod_time)/deltatn
! nfact2 = (next_mod_time - datatimemn)/deltatn
nfact2 = D1_0-nfact1
met_t(:,:,:) = fact1*met_ti(nm)%data(:,:,:) + fact2*met_ti(np)%data(:,:,:)
met_q(:,:,:) = fact1*met_qi(nm)%data(:,:,:) + fact2*met_qi(np)%data(:,:,:)
if (.not. online_test) where (met_q .lt. D0_0) met_q = D0_0
met_ps_next(:,:) = nfact1*met_psi_next(nm)%data(:,:) + nfact2*met_psi_next(np)%data(:,:)
met_ps_curr(:,:) = fact1*met_psi_curr(nm)%data(:,:) + fact2*met_psi_curr(np)%data(:,:)
met_shflx(:,:) = fact1*met_shflxi(nm)%data(:,:) + fact2*met_shflxi(np)%data(:,:)
met_qflx(:,:) = fact1*met_qflxi(nm)%data(:,:) + fact2*met_qflxi(np)%data(:,:)
met_taux(:,:) = fact1*met_tauxi(nm)%data(:,:) + fact2*met_tauxi(np)%data(:,:)
met_tauy(:,:) = fact1*met_tauyi(nm)%data(:,:) + fact2*met_tauyi(np)%data(:,:)
if ( .not. cell_wall_winds ) then
met_u(1:im,jfirst:jlast,kfirst:klast) = fact1*met_ui(nm)%data(1:im,jfirst:jlast,kfirst:klast) &
+ fact2*met_ui(np)%data(1:im,jfirst:jlast,kfirst:klast)
met_v(1:im,jfirst:jlast,kfirst:klast) = fact1*met_vi(nm)%data(1:im,jfirst:jlast,kfirst:klast) &
+ fact2*met_vi(np)%data(1:im,jfirst:jlast,kfirst:klast)
! ghost u,v
#if defined( SPMD )
call mp_send4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast, &
kfirst, klast, ng_d, ng_d, met_u )
call mp_send4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast, &
kfirst, klast, ng_d, ng_s, met_v )
call mp_recv4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast, &
kfirst, klast, ng_d, ng_d, met_u )
call mp_recv4d_ns( grid%commxy, im, jm, km, 1, jfirst, jlast, &
kfirst, klast, ng_d, ng_s, met_v )
#endif
! average to cell walls (vorticity winds)
call transfer_windsToWalls(grid)
else
met_us(:,jfirst:jlast,kfirst:klast) = fact1*met_usi(nm)%data(:,jfirst:jlast,kfirst:klast) + &
fact2*met_usi(np)%data(:,jfirst:jlast,kfirst:klast)
met_vs(:,jfirst:jlast,kfirst:klast) = fact1*met_vsi(nm)%data(:,jfirst:jlast,kfirst:klast) + &
fact2*met_vsi(np)%data(:,jfirst:jlast,kfirst:klast)
endif
! ghost staggered u,v
! WS 2006.04.11: not necessary here since it will be done in cd_core
! write(iulog,*)'INTERPOLATE_METDATA: complete.'
call t_stopf('MET__interpolate_metdata')
end subroutine interpolate_metdata
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine get_dimension( fid, dname, dsize ) 11,2
implicit none
type(file_desc_t), intent(in) :: fid
character(*), intent(in) :: dname
integer, intent(out) :: dsize
integer :: dimid, ierr
ierr = pio_inq_dimid( fid, dname, dimid )
ierr = pio_inq_dimlen( fid, dimid, dsize )
end subroutine get_dimension
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine open_met_datafile( grid, fname, fileid, times, check_dims ) 3,11
use ioFileMod, only: getfil
implicit none
type (T_FVDYCORE_GRID), intent(in) :: grid
character(*), intent(in) :: fname
type(file_desc_t), intent(inout) :: fileid
real(r8), pointer :: times(:)
logical, optional, intent(in) :: check_dims
character(len=256) :: filen
integer :: year, month, day, dsize, i, timesize
integer :: dateid,secid
integer, allocatable , dimension(:) :: dates, datesecs
integer :: ierr
integer :: im, jm, km
im = grid%im
jm = grid%jm
km = grid%km
!
! open file and get fileid
!
call getfil( fname, filen, 0 )
call cam_pio_openfile( fileid, filen, 0 )
if (grid%iam == 0) write(iulog,*)'open_met_datafile: ',trim(filen)
call get_dimension( fileid, 'time', timesize )
if ( associated(times) ) deallocate(times)
allocate( times(timesize) )
allocate( dates(timesize) )
allocate( datesecs(timesize) )
ierr = pio_inq_varid( fileid, 'date', dateid )
ierr = pio_inq_varid( fileid, 'datesec', secid )
ierr = pio_get_var( fileid, dateid, dates )
ierr = pio_get_var( fileid, secid, datesecs )
do i=1,timesize
year = dates(i) / 10000
month = mod(dates(i),10000)/100
day = mod(dates(i),100)
times(i) = get_time_float( year, month, day, datesecs(i) )
enddo
deallocate( dates )
deallocate( datesecs )
!
! check that the data dim sizes match models dimensions
!
if (present(check_dims)) then
if (check_dims) then
call get_dimension( fileid, 'lon', dsize )
if (dsize /= im) then
write(iulog,*)'open_met_datafile: lonsiz=',dsize,' must = ',im
call endrun
endif
call get_dimension( fileid, 'lat', dsize )
if (dsize /= jm) then
write(iulog,*)'open_met_datafile: latsiz=',dsize,' must = ',jm
call endrun
endif
call get_dimension( fileid, 'lev', dsize )
#if ( defined WACCM_GHG || defined WACCM_MOZART )
met_levels = min( dsize, km )
#else
if (dsize /= km) then
write(iulog,*)'open_met_datafile: levsiz=',dsize,' must = ',km
call endrun
endif
met_levels = km
#endif
endif
endif
end subroutine open_met_datafile
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
function get_time_float( year, month, day, sec ) 2,3
! returns float representation of time -- number of days
! since 1 jan 0001 00:00:00.000
implicit none
integer, intent(in) :: year, month, day
integer, intent(in) :: sec
real(r8) :: get_time_float
! ref date is 1 jan 0001
integer :: refyr, refmn, refdy
real(r8) :: refsc, fltdy
integer :: doy(12)
! jan feb mar apr may jun jul aug sep oct nov dec
! 31 28 31 30 31 30 31 31 31 31 30 31
data doy / 1, 32, 60, 91,121,152,182,213,244,274,305,335 /
refyr = 1
refmn = 1
refdy = 1
refsc = D0_0
if ( get_calendar() == 'GREGORIAN' ) then
fltdy = greg2jday(year, month, day) - greg2jday(refyr,refmn,refdy)
else ! assume no_leap (all years are 365 days)
fltdy = (year - refyr)*days_per_non_leapyear + &
(doy(month)-doy(refmn)) + &
(day-refdy)
endif
get_time_float = fltdy + ((sec-refsc)/seconds_per_day)
endfunction get_time_float
!-----------------------------------------------------------------------
! ... Return Julian day number given Gregorian date.
!
! Algorithm from Hatcher,D.A., Simple Formulae for Julian Day Numbers
! and Calendar Dates, Q.Jl.R.astr.Soc. (1984) v25, pp 53-55.
!-----------------------------------------------------------------------
function greg2jday( year, month, day ) 2
implicit none
integer, intent(in) :: year, month, day
integer :: greg2jday
!-----------------------------------------------------------------------
! ... Local variables
!-----------------------------------------------------------------------
integer :: ap, mp
integer :: y, d, n, g
!-----------------------------------------------------------------------
! ... Modify year and month numbers
!-----------------------------------------------------------------------
ap = year - (12 - month)/10
mp = MOD( month-3,12 )
if( mp < 0 ) then
mp = mp + 12
end if
!-----------------------------------------------------------------------
! ... Julian day
!-----------------------------------------------------------------------
y = INT( days_per_year*( ap + 4712 ) )
d = INT( days_per_month*mp + D0_5 )
n = y + d + day + 59
g = INT( D0_75*INT( ap/100 + 49 ) ) - 38
greg2jday = n - g
end function greg2jday
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine set_met_rlx( ) 1,4
use pmgrid
use hycoef, only: hypm, ps0
integer :: k, k_cnt, k_top
real(r8), parameter :: h0 = 7._r8 ! scale height (km)
real(r8) :: p_top, p_bot
996 format( 'set_met_rlx: ',a15, I10 )
997 format( 'set_met_rlx: ',a15, E10.2 )
998 format( 'set_met_rlx: ',a15, PLEV(E10.2))
999 format( 'set_met_rlx: ',a15, PLEV(F10.5))
met_rlx(:) = 999._r8
#if ( defined WACCM_GHG || defined WACCM_MOZART )
p_top = ps0 * exp( - rlx_top_km/h0 )
p_bot = ps0 * exp( - rlx_bot_km/h0 )
if (masterproc) then
write(iulog,fmt=997) 'p_top = ',p_top
write(iulog,fmt=997) 'p_bot = ',p_bot
endif
where( hypm < p_top )
met_rlx = 0._r8
endwhere
where( hypm > p_bot )
met_rlx = max_rlx
endwhere
if ( any( met_rlx(:) /= max_rlx) ) then
k_top = max(plev - met_levels, 1)
do while ( met_rlx(k_top) /= 999._r8 )
k_top = k_top + 1
enddo
met_rlx(1:k_top) = 0._r8
k_cnt = count( met_rlx == 999._r8 )
if (masterproc) then
write(iulog,fmt=996) 'k_cnt = ',k_cnt
write(iulog,fmt=996) 'k_top = ',k_top
endif
do k = k_top,k_top+k_cnt
met_rlx(k) = max_rlx*real( k - k_top ) / real(k_cnt)
enddo
endif
if (masterproc) then
write(iulog,fmt=996) ' met_levels = ',met_levels
write(iulog,fmt=996) 'non-zero terms:',count( met_rlx /= 0._r8 )
endif
if ( met_levels < count( met_rlx /= 0._r8 ) ) then
call endrun('set_met_rlx: met_rlx_top is too high for the meteorology data')
endif
if (masterproc) then
write(iulog,fmt=998) 'press levels = ',hypm
write(iulog,fmt=999) ' met_rlx = ',met_rlx
endif
#else
met_rlx(:) = 1._r8
#endif
if ( any( (met_rlx > 1._r8) .or. (met_rlx < 0._r8) ) ) then
call endrun('Offline meteorology relaxation function not set correctly.')
endif
end subroutine set_met_rlx
end module metdata