diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 425e0349d1..0d1e08ae7d 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -420,7 +420,90 @@
possible_values=".true. or .false."/>
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
@@ -1170,6 +1261,18 @@
#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -2116,6 +2235,28 @@
#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/core_atmosphere/dynamics/Makefile b/src/core_atmosphere/dynamics/Makefile
index 9b0741568d..16f40d0cef 100644
--- a/src/core_atmosphere/dynamics/Makefile
+++ b/src/core_atmosphere/dynamics/Makefile
@@ -6,7 +6,7 @@ OBJS = mpas_atm_time_integration.o \
all: $(OBJS)
-mpas_atm_time_integration.o: mpas_atm_boundaries.o mpas_atm_iau.o mpas_atm_dissipation_models.o
+mpas_atm_time_integration.o: mpas_atm_boundaries.o mpas_atm_iau.o mpas_atm_dissipation_models.o mpas_atm_fdda.o
mpas_atm_boundaries.o:
diff --git a/src/core_atmosphere/dynamics/mpas_atm_fdda.F b/src/core_atmosphere/dynamics/mpas_atm_fdda.F
new file mode 100644
index 0000000000..1a61eba5f5
--- /dev/null
+++ b/src/core_atmosphere/dynamics/mpas_atm_fdda.F
@@ -0,0 +1,434 @@
+!==================================================================================================
+! MPAS four-dimensional data assimilation (FDDA) grid nudging.
+!
+! O. Russell Bullock Jr. (send comments to bullock.russell@epa.gov).
+! 2016-03-30.
+!
+! subroutines in mpas_atm_fdda:
+! ----------------------------------------------
+! apply_fdda : main fdda subroutine (called from subroutine physics_driver).
+! mpas_atm_update_fdda_state : subroutine to read in fdda input stream and save/shift to correct time-level
+!
+!
+
+! add-ons and modifications to sourcecode:
+! ----------------------------------------
+! Original sourcecode by O. Russell Bullock Jr. (bullock.russell@epa.gov) / 2016-03-30.
+!
+! * Added various controls for length-scaled FDDA application
+! O. Russell Bullock Jr. (bullock.russell@epa.gov) / 2016-08-11.
+! * Modified to only require reads of "fdda_new" data arrays. "fdda_old" is initially
+! set to initial conditions. At end of FDDA data interval, "fdda_old" is set to
+! "fdda_new" and updated "fdda_new" data comes from FDDA input stream.
+! O. Russell Bullock Jr. (bullock.russell@epa.gov) / 2016-08-30.
+! * Modified to add variable transformation and adding nudging tendencies to tend_ru, tend_rtheta
+! May Wong (mwong@ucar.edu) / 2025-07-30
+! * Modified to add subroutine mpas_atm_update_fdda_state to update the fdda target data arrays
+! from the input stream.
+! May Wong (mwong@ucar.edu) / 2025-07-30
+
+ module mpas_atm_fdda
+
+ use mpas_kind_types
+ use mpas_derived_types
+ use mpas_pool_routines
+
+ use mpas_atmphys_constants
+ use mpas_atmphys_utilities
+ use mpas_atmphys_vars
+ use mpas_atmphys_todynamics, only : tend_toEdges
+ use mpas_log, only : mpas_log_write
+
+ implicit none
+ private
+ public:: apply_fdda, mpas_atm_update_fdda_state
+
+ integer, private:: i,k
+
+
+ contains
+
+ !***********************************************************************
+ !
+ ! routine apply_fdda
+ !
+ !> \brief Applies FDDA nudging tendencies to the RHS tendencies for dynamics
+ !> \author O. Russell Bullock Jr.
+ !> \date 30 March 2016
+ !> \details
+ !> This routine computes the FDDA weighted nudging tendencies and adds them
+ !> to the 'physics' tendencies for the dynamics.
+ !>
+ !> Variable transformation and adding to the physics tendencies were added.
+ !> May Wong (mwong@ucar.edu) / 2025-07-30
+ !
+ !-----------------------------------------------------------------------
+ subroutine apply_fdda(block,xtime_s,configs,mesh,state,time_lev,fdda,diag,diag_physics,tend,tend_physics,tend_fdda, &
+ tend_ru, tend_rtheta)
+
+ !input variables:
+ type(block_type), intent(in) :: block
+ type(mpas_pool_type),intent(in):: configs
+ type(mpas_pool_type),intent(in):: mesh
+ type(mpas_pool_type),intent(in):: state
+ type(mpas_pool_type),intent(in):: fdda
+ type(mpas_pool_type),intent(in):: diag
+ type(mpas_pool_type),intent(in):: diag_physics
+ integer,intent(in):: time_lev
+ real(kind=RKIND),intent(in):: xtime_s
+
+ !inout variables:
+ type(mpas_pool_type),intent(inout):: tend_physics
+ type(mpas_pool_type),intent(inout):: tend
+ type(mpas_pool_type),intent(inout):: tend_fdda
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: tend_ru
+ real (kind=RKIND), dimension(:,:), intent(inout) :: tend_rtheta
+
+ !local pointers:
+ integer,pointer:: index_qv,fdda_t_min_layer,fdda_q_min_layer,fdda_uv_min_layer
+ integer,dimension(:),pointer :: kpbl
+ integer, pointer :: nCellsSolve, nVertLevels, nEdgesSolve
+ real(kind=RKIND),pointer:: fdda_t_coef,fdda_q_coef,fdda_uv_coef
+ real(kind=RKIND),pointer:: fdda_int,fdda_scale_min,fdda_scale_max
+ real(kind=RKIND),pointer:: config_dt
+ real(kind=RKIND),pointer:: len_disp
+ real(kind=RKIND),dimension(:),pointer :: meshDensity
+ real(kind=RKIND),dimension(:,:),pointer :: rthfddaten,rqvfddaten
+ real(kind=RKIND),dimension(:,:),pointer :: rufddaten,rvfddaten,rufddaten_Edge
+ real(kind=RKIND),dimension(:,:),pointer :: qv_fdda_old,theta_fdda_old,u_fdda_old,v_fdda_old
+ real(kind=RKIND),dimension(:,:),pointer :: qv_fdda_new,theta_fdda_new,u_fdda_new,v_fdda_new
+ real(kind=RKIND),dimension(:,:),pointer :: u,v
+ real(kind=RKIND),dimension(:,:),pointer :: qv
+ real(kind=RKIND),dimension(:,:),pointer :: theta_m
+ real(kind=RKIND),dimension(:,:,:),pointer:: scalars, tend_scalars
+ real(kind=RKIND),dimension(:,:),pointer:: rho_zz
+ real(kind=RKIND),dimension(:,:),pointer:: rho_edge
+ real (kind=RKIND), dimension(:,:), pointer :: tend_th ! tendency for potential temperature
+ character(len=StrKIND),pointer:: fdda_scheme
+ logical,pointer:: fdda_t, fdda_t_in_pbl
+ logical,pointer:: fdda_q, fdda_q_in_pbl
+ logical,pointer:: fdda_uv,fdda_uv_in_pbl
+ logical,pointer:: config_do_restart
+
+ !local variables:
+ real(kind=RKIND) :: tfrac
+ real(kind=RKIND) :: t_target,q_target,u_target,v_target
+ real(kind=RKIND) :: sfac
+ real(kind=RKIND) :: theta_p
+ real(kind=RKIND) :: srange,dcEdge_m
+
+ integer :: i,k
+
+ logical :: first_call = .true.
+
+ !local parameters:
+
+ !==================================================================================================
+
+ call mpas_pool_get_config(configs,'config_do_restart' ,config_do_restart )
+ call mpas_pool_get_config(configs,'config_dt' ,config_dt )
+ call mpas_pool_get_config(configs,'config_fdda_scheme', fdda_scheme)
+ call mpas_pool_get_config(configs,'config_fdda_int',fdda_int)
+ call mpas_pool_get_config(configs,'config_fdda_t',fdda_t)
+ call mpas_pool_get_config(configs,'config_fdda_q',fdda_q)
+ call mpas_pool_get_config(configs,'config_fdda_uv',fdda_uv)
+ call mpas_pool_get_config(configs,'config_fdda_t_in_pbl',fdda_t_in_pbl)
+ call mpas_pool_get_config(configs,'config_fdda_q_in_pbl',fdda_q_in_pbl)
+ call mpas_pool_get_config(configs,'config_fdda_uv_in_pbl',fdda_uv_in_pbl)
+ call mpas_pool_get_config(configs,'config_fdda_t_min_layer',fdda_t_min_layer)
+ call mpas_pool_get_config(configs,'config_fdda_q_min_layer',fdda_q_min_layer)
+ call mpas_pool_get_config(configs,'config_fdda_uv_min_layer',fdda_uv_min_layer)
+ call mpas_pool_get_config(configs,'config_fdda_t_coef',fdda_t_coef)
+ call mpas_pool_get_config(configs,'config_fdda_q_coef',fdda_q_coef)
+ call mpas_pool_get_config(configs,'config_fdda_uv_coef',fdda_uv_coef)
+ call mpas_pool_get_config(configs,'config_fdda_scale_min',fdda_scale_min)
+ call mpas_pool_get_config(configs,'config_fdda_scale_max',fdda_scale_max)
+ call mpas_pool_get_config(configs,'config_len_disp' ,len_disp)
+
+ call mpas_pool_get_array(mesh,'meshDensity',meshDensity)
+ call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve)
+ call mpas_pool_get_dimension(mesh, 'nEdgesSolve', nEdgesSolve)
+ call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
+
+ call mpas_pool_get_array(fdda,'qv_fdda',qv_fdda_old , 1)
+ call mpas_pool_get_array(fdda,'theta_fdda',theta_fdda_old, 1)
+ call mpas_pool_get_array(fdda,'u_fdda', u_fdda_old , 1)
+ call mpas_pool_get_array(fdda,'v_fdda', v_fdda_old , 1)
+ call mpas_pool_get_array(fdda,'qv_fdda',qv_fdda_new , 2)
+ call mpas_pool_get_array(fdda,'theta_fdda',theta_fdda_new, 2)
+ call mpas_pool_get_array(fdda,'u_fdda', u_fdda_new , 2)
+ call mpas_pool_get_array(fdda,'v_fdda', v_fdda_new , 2)
+ call mpas_pool_get_array(tend_fdda,'rthfddaten',rthfddaten)
+ call mpas_pool_get_array(tend_fdda,'rqvfddaten',rqvfddaten)
+ call mpas_pool_get_array(tend_fdda,'rufddaten_Edge' ,rufddaten_Edge )
+ call mpas_pool_get_array(tend_fdda,'rufddaten' ,rufddaten )
+ call mpas_pool_get_array(tend_fdda,'rvfddaten' ,rvfddaten )
+ call mpas_pool_get_array(diag_physics,'kpbl',kpbl)
+ call mpas_pool_get_array(diag,'uReconstructZonal' ,u)
+ call mpas_pool_get_array(diag,'uReconstructMeridional',v)
+ call mpas_pool_get_array(tend,'scalars_tend' , tend_scalars )
+
+ call mpas_pool_get_dimension(state,'index_qv',index_qv)
+ call mpas_pool_get_array(state,'scalars',scalars,time_lev)
+ qv => scalars(index_qv,:,:)
+ call mpas_pool_get_array(state,'theta_m',theta_m,time_lev)
+
+ call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2)
+ call mpas_pool_get_array(diag , 'rho_edge', rho_edge)
+
+ allocate(tend_th(nVertLevels,nCellsSolve))
+
+
+ !calculate time fraction within FDDA target data interval
+ tfrac = mod(xtime_s,fdda_int)/fdda_int
+ !call mpas_log_write('inside FDDA: tfrac = $r', realArgs=(/tfrac/))
+
+ fdda_select: select case(fdda_scheme)
+
+ case ("scaled")
+
+ srange = fdda_scale_max-fdda_scale_min
+
+ DO i=1, nCellsSolve
+
+ dcEdge_m = len_disp / meshDensity(i)**0.25_RKIND
+
+ if (dcEdge_m .lt. fdda_scale_max) then
+ if (dcEdge_m .lt. fdda_scale_min) then
+ sfac = 0.0_RKIND
+ else
+ sfac = (dcEdge_m-fdda_scale_min)/srange
+ endif
+ else
+ sfac = 1.0_RKIND
+ endif
+
+ DO k=1,nVertLevels
+
+ theta_p = theta_m(k,i) / (1._RKIND + rvord * qv(k,i))
+
+ if (fdda_t .and. k.gt.fdda_t_min_layer .and. (fdda_t_in_pbl .or. k.gt.kpbl(i))) then
+ t_target = (1.0_RKIND-tfrac)*theta_fdda_old(k,i)+tfrac*theta_fdda_new(k,i)
+ rthfddaten(k,i) = fdda_t_coef * ( t_target - theta_p ) * sfac
+
+ else
+ rthfddaten(k,i) = 0._RKIND
+ endif
+
+ if (fdda_q .and. k.gt.fdda_q_min_layer .and. (fdda_q_in_pbl .or. k.gt.kpbl(i))) then
+ q_target = (1.0_RKIND-tfrac)*qv_fdda_old(k,i)+tfrac*qv_fdda_new(k,i)
+ rqvfddaten(k,i) = fdda_q_coef * ( q_target - qv(k,i) ) * sfac
+ else
+ rqvfddaten(k,i) = 0._RKIND
+ endif
+
+ if (fdda_uv .and. k.gt.fdda_uv_min_layer .and. (fdda_uv_in_pbl .or. k.gt.kpbl(i))) then
+ u_target = (1.0_RKIND-tfrac)*u_fdda_old(k,i)+tfrac*u_fdda_new(k,i)
+ rufddaten(k,i) = fdda_uv_coef * ( u_target - u(k,i) ) * sfac
+ v_target = (1.0_RKIND-tfrac)*v_fdda_old(k,i)+tfrac*v_fdda_new(k,i)
+ rvfddaten(k,i) = fdda_uv_coef * ( v_target - v(k,i) ) * sfac
+ else
+ rufddaten(k,i) = 0._RKIND
+ rvfddaten(k,i) = 0._RKIND
+ endif
+
+ ENDDO
+ ENDDO
+
+ case("analysis")
+
+ DO i = 1, nCellsSolve
+ DO k = 1, nVertLevels
+
+ theta_p = theta_m(k,i) / (1._RKIND + rvord * qv(k,i))
+
+ if (fdda_t .and. k.gt.fdda_t_min_layer .and. (fdda_t_in_pbl .or. k.gt.kpbl(i))) then
+ t_target = (1.0_RKIND-tfrac)*theta_fdda_old(k,i)+tfrac*theta_fdda_new(k,i)
+ rthfddaten(k,i) = fdda_t_coef * ( t_target - theta_p )
+ !if (i.eq.1 .and. k.eq.1 ) call mpas_log_write('t_target, theta_p, theta_fdda_old, theta_fdda_new = $r, $r, $r, $r', &
+ ! realArgs=(/t_target, theta_p, theta_fdda_old(k,i), theta_fdda_new(k,i)/) )
+ else
+ rthfddaten(k,i) = 0._RKIND
+ endif
+
+ if (fdda_q .and. k.gt.fdda_q_min_layer .and. (fdda_q_in_pbl .or. k.gt.kpbl(i))) then
+ q_target = (1.0_RKIND-tfrac)*qv_fdda_old(k,i)+tfrac*qv_fdda_new(k,i)
+ rqvfddaten(k,i) = fdda_q_coef * ( q_target - qv(k,i) )
+ else
+ rqvfddaten(k,i) = 0._RKIND
+ endif
+
+ if (fdda_uv .and. k.gt.fdda_uv_min_layer .and. (fdda_uv_in_pbl .or. k.gt.kpbl(i))) then
+ u_target = (1.0_RKIND-tfrac)*u_fdda_old(k,i)+tfrac*u_fdda_new(k,i)
+ rufddaten(k,i) = fdda_uv_coef * ( u_target - u(k,i) )
+ v_target = (1.0_RKIND-tfrac)*v_fdda_old(k,i)+tfrac*v_fdda_new(k,i)
+ rvfddaten(k,i) = fdda_uv_coef * ( v_target - v(k,i) )
+ else
+ rufddaten(k,i) = 0._RKIND
+ rvfddaten(k,i) = 0._RKIND
+ endif
+
+ ENDDO
+ ENDDO
+ !call mpas_log_write('in apply_fdda: min,max rthfddaten = $r, $r', realArgs=(/minval(rthfddaten), maxval(rthfddaten)/) )
+ !call mpas_log_write('in apply_fdda: min,max rqvfddaten = $r, $r', realArgs=(/minval(rqvfddaten), maxval(rqvfddaten)/) )
+ !call mpas_log_write('in apply_fdda: min,max rufddaten = $r, $r', realArgs=(/minval(rufddaten), maxval(rufddaten)/) )
+ !call mpas_log_write('in apply_fdda: min,max rvfddaten = $r, $r', realArgs=(/minval(rvfddaten), maxval(rvfddaten)/) )
+
+ case default
+
+ end select fdda_select
+
+
+ ! MW: add transformation of the u,v wind nudging tendencies to those for horizontal normal velocity at the edge
+ rufddaten_Edge(:,:) = 0.0_RKIND
+ call tend_toEdges(block, mesh, rufddaten, rvfddaten, rufddaten_Edge)
+
+ ! MW: add coupled tendencies for u on edges
+ do i = 1, nEdgesSolve
+ do k = 1, nVertLevels
+ tend_ru(k,i)=tend_ru(k,i)+rufddaten_Edge(k,i)*rho_edge(k,i)
+ enddo
+ enddo
+
+ ! MW: add coupled tendencies for other variables at cell centers
+ tend_th = 0.0_RKIND
+ do i = 1, nCellsSolve
+ do k = 1, nVertLevels
+ tend_th(k,i) = tend_th(k,i) + rthfddaten(k,i)*rho_zz(k,i)
+ tend_scalars(index_qv,k,i) = tend_scalars(index_qv,k,i) + rqvfddaten(k,i)*rho_zz(k,i)
+ enddo
+ enddo
+
+
+ ! if non-hydrostatic core, convert the tendency for the potential temperature to a
+ ! tendency for the modified potential temperature
+ do i = 1, nCellsSolve
+ do k = 1,nVertLevels
+ tend_th(k,i) = (1. + rvord * scalars(index_qv,k,i)) * tend_th(k,i) &
+ + rvord * theta_m(k,i) * rqvfddaten(k,i)*rho_zz(k,i)
+ tend_rtheta(k,i) = tend_rtheta(k,i) + tend_th(k,i)
+ end do
+ end do
+
+
+ deallocate(tend_th)
+
+ end subroutine apply_fdda
+
+
+ !***********************************************************************
+ !
+ ! routine mpas_atm_update_fdda_state
+ !
+ !> \brief Reads new fdda input data and updates the FDDA tendencies
+ !> \author Michael Duda, May Wong
+ !> \date 30 July 2025
+ !> \details
+ !> This routine reads from the 'fdda' stream all variables in the 'fdda'
+ !> pool. When called with firstCall=.true., the latest time before the
+ !> present is read into time level 2 of the fdda pool; otherwise, the
+ !> contents of time level 2 are shifted to time level 1, the earliest
+ !> time strictly later than the present is read into time level 2, and
+ !> the tendencies for all fields in the fdda pool are computed and stored
+ !> in time level 1.
+ !>
+ !> This subroutine for FDDA tendencies is based on mpas_atm_update_bdy_tend.
+ !
+ !-----------------------------------------------------------------------
+ subroutine mpas_atm_update_fdda_state(clock, streamManager, block, firstCall, ierr)
+
+ use mpas_constants, only : rvord
+ use mpas_stream_manager, only : mpas_stream_mgr_read
+ use mpas_log, only : mpas_log_write
+ use mpas_derived_types, only : MPAS_STREAM_MGR_NOERR, MPAS_LOG_ERR
+ use mpas_timekeeping, only : mpas_get_time
+
+ implicit none
+
+ type (mpas_clock_type), intent(in) :: clock
+ type (MPAS_streamManager_type), intent(inout) :: streamManager
+ type (block_type), intent(inout) :: block
+ logical, intent(in) :: firstCall
+ integer, intent(out) :: ierr
+
+ character(len=StrKIND) :: lbc_intv_start_string
+ character(len=StrKIND) :: lbc_intv_end_string
+
+ type (mpas_pool_type), pointer :: mesh
+ type (mpas_pool_type), pointer :: state
+ type (mpas_pool_type), pointer :: fdda
+ real (kind=RKIND) :: dt
+
+ integer, pointer :: nCells
+ integer, pointer :: nEdges
+ integer, pointer :: index_qv
+
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real (kind=RKIND), dimension(:,:), pointer :: ru
+ real (kind=RKIND), dimension(:,:), pointer :: rho_edge
+ real (kind=RKIND), dimension(:,:), pointer :: w
+ real (kind=RKIND), dimension(:,:), pointer :: theta
+ real (kind=RKIND), dimension(:,:), pointer :: rtheta_m
+ real (kind=RKIND), dimension(:,:), pointer :: rho_zz
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+
+ real (kind=RKIND), dimension(:,:), pointer :: qv_fdda
+ real (kind=RKIND), dimension(:,:), pointer :: theta_fdda
+ real (kind=RKIND), dimension(:,:), pointer :: u_fdda
+ real (kind=RKIND), dimension(:,:), pointer :: v_fdda
+
+
+ integer :: dd_intv, s_intv, sn_intv, sd_intv
+ type (MPAS_Time_Type) :: currTime
+ type (MPAS_TimeInterval_Type) :: lbc_interval
+ character(len=StrKIND) :: read_time
+ integer :: iEdge
+ integer :: cell1, cell2
+
+
+ ierr = 0
+
+ call mpas_pool_get_subpool(block % structs, 'mesh', mesh)
+ call mpas_pool_get_subpool(block % structs, 'state', state)
+ call mpas_pool_get_subpool(block % structs, 'fdda', fdda)
+
+ if (firstCall) then
+ call MPAS_stream_mgr_read(streamManager, streamID='fdda', timeLevel=1, whence=MPAS_STREAM_LATEST_BEFORE, &
+ actualWhen=read_time, ierr=ierr)
+ if (ierr /= MPAS_STREAM_MGR_NOERR) then
+ call mpas_log_write('Could not read from ''fdda'' stream on or before the current date '// &
+ 'to update target FDDA fields', messageType=MPAS_LOG_ERR)
+ ierr = 1
+ end if
+
+ call MPAS_stream_mgr_read(streamManager, streamID='fdda', timeLevel=2, whence=MPAS_STREAM_EARLIEST_STRICTLY_AFTER, &
+ actualWhen=read_time, ierr=ierr)
+ if (ierr /= MPAS_STREAM_MGR_NOERR) then
+ call mpas_log_write('Could not read from ''fdda'' stream after the current date '// &
+ 'to update target FDDA fields', messageType=MPAS_LOG_ERR)
+
+ ierr = 1
+ end if
+
+ else
+ call mpas_pool_shift_time_levels(fdda)
+ call MPAS_stream_mgr_read(streamManager, streamID='fdda', timeLevel=2, whence=MPAS_STREAM_EARLIEST_STRICTLY_AFTER, &
+ actualWhen=read_time, ierr=ierr)
+ if (ierr /= MPAS_STREAM_MGR_NOERR) then
+ call mpas_log_write('Could not read from ''fdda'' stream after the current date '// &
+ 'to update target FDDA fields', messageType=MPAS_LOG_ERR)
+
+ ierr = 1
+ end if
+ end if
+ if (ierr /= 0) then
+ return
+ end if
+
+ end subroutine mpas_atm_update_fdda_state
+
+
+ end module mpas_atm_fdda
diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
index 227fbde862..71d8390014 100644
--- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
+++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
@@ -36,6 +36,7 @@ module atm_time_integration
use mpas_atm_boundaries, only : nSpecZone, nRelaxZone, nBdyZone, mpas_atm_get_bdy_state, mpas_atm_get_bdy_tend ! regional_MPAS addition
use mpas_atm_iau
+ use mpas_atm_fdda
use mpas_atm_dissipation_models
! Provides definition of halo_exchange_routine
@@ -1786,7 +1787,7 @@ subroutine mpas_atm_dynamics_finalize(domain)
end subroutine mpas_atm_dynamics_finalize
- subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group)
+ subroutine atm_timestep(domain, dt, nowTime, itimestep, xtime_s, exchange_halo_group)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
!
@@ -1802,6 +1803,7 @@ subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group)
real (kind=RKIND), intent(in) :: dt
type (MPAS_Time_type), intent(in) :: nowTime
integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: xtime_s
procedure (halo_exchange_routine) :: exchange_halo_group
@@ -1826,7 +1828,7 @@ subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group)
call mpas_atm_pre_dynamics(domain)
if (trim(config_time_integration) == 'SRK3') then
- call atm_srk3(domain, dt, itimestep, exchange_halo_group)
+ call atm_srk3(domain, dt, itimestep, xtime_s, exchange_halo_group)
else
call mpas_log_write('Unknown time integration option '//trim(config_time_integration), messageType=MPAS_LOG_ERR)
call mpas_log_write('Currently, only ''SRK3'' is supported.', messageType=MPAS_LOG_CRIT)
@@ -1852,7 +1854,7 @@ subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group)
end subroutine atm_timestep
- subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
+ subroutine atm_srk3(domain, dt, itimestep, xtime_s, exchange_halo_group)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! time-split RK3 scheme
@@ -1871,6 +1873,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: xtime_s
procedure (halo_exchange_routine) :: exchange_halo_group
integer :: thread
@@ -1912,14 +1915,17 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
integer :: nCells, nEdges, nVertices, nEdgesSolve, nVertLevels, num_scalars
character(len=StrKIND), pointer :: config_IAU_option
+ character(len=StrKIND), pointer :: config_fdda_scheme
type (mpas_pool_type), pointer :: state
type (mpas_pool_type), pointer :: diag
type (mpas_pool_type), pointer :: diag_physics
type (mpas_pool_type), pointer :: mesh
type (mpas_pool_type), pointer :: tend
+ type (mpas_pool_type), pointer :: tend_fdda => null()
type (mpas_pool_type), pointer :: tend_physics
type (mpas_pool_type), pointer :: lbc ! regional_MPAS addition
+ type (mpas_pool_type), pointer :: fdda
real (kind=RKIND), dimension(:,:), pointer :: w
real (kind=RKIND), dimension(:,:), pointer :: u, uReconstructZonal, uReconstructMeridional, uReconstructX, uReconstructY, uReconstructZ
@@ -1957,6 +1963,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
call mpas_pool_get_config(block % configs, 'config_microp_scheme', config_microp_scheme)
call mpas_pool_get_config(block % configs, 'config_convection_scheme', config_convection_scheme)
#endif
+ call mpas_pool_get_config(block % configs, 'config_fdda_scheme', config_fdda_scheme)
!
! Retrieve field structures
@@ -1966,9 +1973,11 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
call mpas_pool_get_subpool(block % structs, 'diag', diag)
call mpas_pool_get_subpool(block % structs, 'tend', tend)
call mpas_pool_get_subpool(block % structs, 'tend_physics', tend_physics)
+ call mpas_pool_get_subpool(block % structs, 'tend_fdda', tend_fdda)
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics)
#endif
+ call mpas_pool_get_subpool(block % structs, 'fdda', fdda)
!
! Retrieve dimensions
@@ -2173,6 +2182,17 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
tend_rtheta_physics, tend_rho_physics)
end if
+ !
+ ! FDDA - Four-dimensional data assimilation (Grid point nudging)
+ !
+ ! MW: Why are tend_ru_physics and tend_rtheta_physics global variables, but not tend_scalars?
+ if (trim(config_fdda_scheme) /= 'off') then
+ call apply_fdda( block, xtime_s, block % configs, mesh, state, 1, fdda, diag, diag_physics, tend, tend_physics, tend_fdda, &
+ tend_ru_physics, tend_rtheta_physics)
+ end if
+
+
+
DYNAMICS_SUBSTEPS : do dynamics_substep = 1, dynamics_split
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index eef1951e36..6d49af72a5 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -626,6 +626,7 @@ function atm_core_run(domain) result(ierr)
use mpas_timer, only : mpas_timer_start, mpas_timer_stop
use mpas_atm_boundaries, only : mpas_atm_update_bdy_tend
use mpas_atm_diagnostics_manager, only : mpas_atm_diag_update, mpas_atm_diag_compute, mpas_atm_diag_reset
+ use mpas_atm_fdda, only : mpas_atm_update_fdda_state
implicit none
@@ -635,6 +636,7 @@ function atm_core_run(domain) result(ierr)
real (kind=RKIND) :: dt
logical, pointer :: config_do_restart
logical, pointer :: config_apply_lbcs
+ character(len=StrKIND),pointer:: config_fdda_scheme
type (block_type), pointer :: block_ptr
type (MPAS_Time_Type) :: currTime
@@ -728,6 +730,7 @@ function atm_core_run(domain) result(ierr)
#endif
call mpas_pool_get_config(domain % blocklist % configs, 'config_apply_lbcs', config_apply_lbcs)
+ call mpas_pool_get_config(domain % blocklist % configs, 'config_fdda_scheme', config_fdda_scheme)
!
! Read initial boundary state
@@ -748,6 +751,27 @@ function atm_core_run(domain) result(ierr)
end do
end if
+ !
+ ! Read initial FDDA target state
+ !
+ if ( trim(config_fdda_scheme) /= 'off' .and. &
+ MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='fdda', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call mpas_atm_update_fdda_state(clock, domain % streamManager, block_ptr, .true., ierr)
+ if (ierr /= 0) then
+ currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
+ call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp)
+ call mpas_log_write('Failed to process FDDA data on or before '//trim(timeStamp), messageType=MPAS_LOG_ERR)
+ return
+ end if
+
+ block_ptr => block_ptr % next
+ end do
+ end if
+
+
+
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
itimestep = 1
@@ -776,10 +800,33 @@ function atm_core_run(domain) result(ierr)
end do
end if
+ !
+ ! Read future FDDA target state and shift time levels
+ !
+ if ( itimestep.ne.1 .and. trim(config_fdda_scheme) /= 'off' .and. &
+ MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='fdda', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call mpas_atm_update_fdda_state(clock, domain % streamManager, block_ptr, .false., ierr)
+ if (ierr /= 0) then
+ call mpas_log_write('Failed to process FDDA data at next time after '//trim(timeStamp), messageType=MPAS_LOG_ERR)
+ return
+ end if
+
+ block_ptr => block_ptr % next
+ end do
+ end if
+
+
+
+
! Regardless of whether boundary tendencies were updated, above, we do not want to read the 'lbc_in' stream
! as a general input stream, below.
call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)
+ ! MW: Similar to the 'lbc_in' stream, we want to reset the alarm for the fdda stream
+ call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='fdda', direction=MPAS_STREAM_INPUT, ierr=ierr)
+
!
! Read external field updates
@@ -818,6 +865,8 @@ function atm_core_run(domain) result(ierr)
call mpas_timer_stop("time integration")
call mpas_log_write(' Timing for integration step: $r s', realArgs=(/real(integ_stop_time - integ_start_time, kind=RKIND)/))
+
+
! Move time level 2 fields back into time level 1 for next time step
call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state)
call mpas_pool_shift_time_levels(state)
@@ -827,7 +876,6 @@ function atm_core_run(domain) result(ierr)
call mpas_advance_clock(clock)
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
-
call mpas_dmpar_get_time(diag_start_time)
!
@@ -1041,7 +1089,7 @@ subroutine atm_do_timestep(domain, dt, itimestep)
call chemistry_step()
#endif
- call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group)
+ call atm_timestep(domain, dt, currTime, itimestep, xtime_s, exchange_halo_group)
end subroutine atm_do_timestep
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver.F b/src/core_atmosphere/physics/mpas_atmphys_driver.F
index 8e31672657..6d7bbf8582 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver.F
@@ -88,13 +88,13 @@ module mpas_atmphys_driver
! Laura D. Fowler (laura@ucar.edu) / 2014-05-15.
! * renamed config_conv_deep_scheme to config_convection_scheme.
! Laura D. Fowler (laura@ucar.edu) / 2014-09-18.
-! * in the call to driver_convection, added block%configs needed for the implementation of the
-! Grell-Freitas convection scheme.
-! Laura D. Fowler (laura@ucar.edu) / 2016-03-30.
! * modified the call to the subroutines driver_sfclayer and driver_pbl for the implementation
! of the MYNN surface layer scheme and PBL schemes. itimestep and block%configs are added to
! the argument list.
! Laura D. Fowler (laura@ucar.edu) / 2015-01-06.
+! * in the call to driver_convection, added block%configs needed for the implementation of the
+! Grell-Freitas convection scheme.
+! Laura D. Fowler (laura@ucar.edu) / 2016-03-30.
! * now only call subroutine update_convection_step2 when config_convection_scheme is not off.
! Laura D. Fowler (laura@ucar.edu) / 2016-04-13.
! * modified call to driver_cloudiness to accomodate the calculation of the cloud fraction with the Thompson
@@ -177,7 +177,7 @@ subroutine physics_driver(domain,itimestep,xtime_s)
config_pbl_scheme .ne. 'off' .or. &
config_radt_lw_scheme .ne. 'off' .or. &
config_radt_sw_scheme .ne. 'off' .or. &
- config_sfclayer_scheme .ne. 'off') then
+ config_sfclayer_scheme .ne. 'off') then
block => domain % blocklist
do while(associated(block))
diff --git a/src/core_atmosphere/physics/mpas_atmphys_todynamics.F b/src/core_atmosphere/physics/mpas_atmphys_todynamics.F
index 39f3aea841..8ac6dc760d 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_todynamics.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_todynamics.F
@@ -26,7 +26,11 @@ module mpas_atmphys_todynamics
implicit none
private
+<<<<<<< HEAD
+ public:: physics_get_tend, pre_physics_get_tend, post_physics_get_tend, tend_toEdges
+=======
public:: physics_get_tend, pre_physics_get_tend, post_physics_get_tend
+>>>>>>> 91c5eac175eebeaf4206bacd5cb50c39dff3c152
!Interface between the physics parameterizations and the non-hydrostatic dynamical core.
@@ -215,7 +219,6 @@ subroutine physics_get_tend(block,mesh,state,diag,tend,tend_physics,configs,rk_s
tend_rtheta_physics(:,:) = 0._RKIND
tend_rho_physics(:,:) = 0._RKIND
-
!in case some variables are not allocated due to their associated packages. We need to make their pointers
!associated here to avoid triggering run-time. checks when calling physics_get_tend_work:
if(.not. associated(rucuten) ) allocate(rucuten(0,0) )
@@ -274,6 +277,7 @@ subroutine physics_get_tend(block,mesh,state,diag,tend,tend_physics,configs,rk_s
if(size(rnifablten) == 0) deallocate(rnifablten)
if(size(rnwfablten) == 0) deallocate(rnwfablten)
+
deallocate(tend_th)
end subroutine physics_get_tend
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index 831901565a..48697c7177 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -377,6 +377,7 @@
+
@@ -622,6 +623,9 @@
+
+
+
@@ -690,6 +694,23 @@
+
+
+
+
+
+
+
+
+
+
+
@@ -1219,6 +1240,53 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index fb67a5164b..a074404d07 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -63,6 +63,7 @@ subroutine init_atm_setup_case(domain, stream_manager)
type (mpas_pool_type), pointer :: diag
type (mpas_pool_type), pointer :: diag_physics
type (mpas_pool_type), pointer :: lbc_state
+ type (mpas_pool_type), pointer :: fdda_state
integer, pointer :: config_init_case
logical, pointer :: config_static_interp
@@ -80,7 +81,7 @@ subroutine init_atm_setup_case(domain, stream_manager)
real (kind=RKIND), pointer :: Time
type (MPAS_Time_type) :: curr_time, stop_time, start_time
- type (MPAS_TimeInterval_type) :: clock_interval, lbc_stream_interval, surface_stream_interval
+ type (MPAS_TimeInterval_type) :: clock_interval, lbc_stream_interval, surface_stream_interval, fdda_stream_interval
type (MPAS_TimeInterval_type) :: time_since_start
character(len=StrKIND) :: timeStart,timeString
@@ -389,6 +390,69 @@ subroutine init_atm_setup_case(domain, stream_manager)
block_ptr => block_ptr % next
end do
+ else if (config_init_case == 11 ) then
+
+ call mpas_log_write('FDDA analysis nudging case')
+
+ !
+ ! Check that the first-guess interval (which is the same as the clock timestep)
+ ! matches the output interval of the 'fdda' stream
+ !
+ clock_interval = mpas_get_clock_timestep(domain % clock, ierr=ierr)
+ fdda_stream_interval = MPAS_stream_mgr_get_stream_interval(stream_manager, 'fdda', MPAS_STREAM_OUTPUT, ierr)
+ if (clock_interval /= fdda_stream_interval) then
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('The intermediate file interval specified by ''config_fg_interval''', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('does not match the output_interval for the ''fdda'' stream.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Please correct the namelist.init_atmosphere and/or', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('streams.init_atmosphere files.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+
+ curr_time = mpas_get_clock_time(domain % clock, MPAS_NOW)
+ stop_time = mpas_get_clock_time(domain % clock, MPAS_STOP_TIME)
+
+ do while (curr_time <= stop_time)
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
+ call mpas_pool_get_subpool(block_ptr % structs, 'fg', fg)
+ call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
+ call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
+ call mpas_pool_get_subpool(block_ptr % structs, 'fdda_state', fdda_state)
+
+ call mpas_pool_get_array(state, 'xtime', xtime)
+
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nCells', nCells)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nEdges', nEdges)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nVertLevels', nVertLevels)
+
+ call mpas_get_time(curr_time, dateTimeString=timeString)
+ xtime = timeString ! Set field valid time, xtime, to the current time in the time loop
+
+ call init_atm_case_fdda(timeString, block_ptr, mesh, nCells, nEdges, nVertLevels, fg, state, &
+ diag, fdda_state, block_ptr % dimensions, block_ptr % configs)
+
+
+ block_ptr => block_ptr % next
+ end do
+
+ call mpas_stream_mgr_write(stream_manager, streamID='fdda', ierr=ierr)
+ call mpas_stream_mgr_reset_alarms(stream_manager, streamID='fdda', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
+
+ call mpas_advance_clock(domain % clock)
+ curr_time = mpas_get_clock_time(domain % clock, MPAS_NOW)
+
+ end do
+
+ !
+ ! Ensure that no output alarms are still ringing for the 'fdda' stream after
+ ! we exit the time loop above; the main run routine may write out all other
+ ! output streams with ringing alarms.
+ !
+ call mpas_stream_mgr_reset_alarms(stream_manager, streamID='fdda', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
+
else if (config_init_case == 13 ) then
call mpas_log_write(' CAM-MPAS grid ')
@@ -6903,6 +6967,841 @@ real (kind=RKIND) function atm_get_sounding_2( variable, height )
end function atm_get_sounding_2
!-----------
+<<<<<<< HEAD
+
+ !-----------------------------------------------------------------------
+ ! routine init_atm_case_fdda
+ !
+ !> \brief Computes theta_fdda, qv_fdda, u_fdda, v_fdda fields for FDDA nudging
+ !> \author Michael Duda, May Wong
+ !> \date 2025-07-30
+ !> \details
+ !> This routine is similar to the init_atm_case_lbc routine in that it reads
+ !> atmospheric fields from "intermediate" files and horizontally and vertically
+ !> interpolates them to an MPAS mesh.However, it also includes the reconstructed
+ !> wind fields as originally in the FDDA nudging code. This routine is responsible for
+ !> producing only those fields that are needed for FDDA nudging.
+ !
+ !-----------------------------------------------------------------------
+ subroutine init_atm_case_fdda(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, fdda_state, dims, configs)
+
+ use mpas_dmpar, only : mpas_dmpar_min_real, mpas_dmpar_max_real
+ use init_atm_read_met, only : met_data, read_met_init, read_met_close, read_next_met_field
+ use init_atm_llxy, only : proj_info, map_init, map_set, latlon_to_ij, PROJ_LATLON, PROJ_GAUSS, DEG_PER_RAD
+ use init_atm_hinterp, only : interp_sequence, FOUR_POINT, SIXTEEN_POINT, W_AVERAGE4, SEARCH
+ use mpas_hash, only : hashtable, mpas_hash_init, mpas_hash_destroy, mpas_hash_search, mpas_hash_size, mpas_hash_insert
+
+ implicit none
+
+ character(len=*), intent(in) :: timestamp
+ type (block_type), intent(inout), target :: block
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertLevels
+ type (mpas_pool_type), intent(inout) :: fg
+ type (mpas_pool_type), intent(inout) :: state
+ type (mpas_pool_type), intent(inout) :: diag
+ type (mpas_pool_type), intent(inout):: fdda_state
+ type (mpas_pool_type), intent(inout):: dims
+ type (mpas_pool_type), intent(inout):: configs
+
+ type (parallel_info), pointer :: parinfo
+ type (dm_info), pointer :: dminfo
+
+ real (kind=RKIND), parameter :: t0b = 250.0
+
+ type (met_data) :: field
+ type (proj_info) :: proj
+
+ real (kind=RKIND), dimension(:), pointer :: dzu, fzm, fzp
+ real (kind=RKIND), dimension(:), pointer :: vert_level, latPoints, lonPoints
+ real (kind=RKIND), dimension(:,:), pointer :: zgrid, zz
+ real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, t, rt
+ real (kind=RKIND), dimension(:), pointer :: destField1d
+ real (kind=RKIND), dimension(:,:), pointer :: destField2d
+ real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+
+ real (kind=RKIND) :: target_z
+ integer :: iCell, iEdge, i, k, nVertLevelsP1
+ integer, pointer :: nCellsSolve
+ integer :: nInterpPoints, ndims
+
+ integer :: nfglevels_actual
+ integer, pointer :: index_qv
+
+ integer, dimension(5) :: interp_list
+ real (kind=RKIND) :: msgval
+
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+ real (kind=RKIND), dimension(:,:), pointer :: sorted_arr
+
+ integer :: sfc_k
+
+ integer :: it
+ real (kind=RKIND) :: p_check
+
+ integer :: istatus
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: rslab
+
+ real (kind=RKIND) :: flux
+ real (kind=RKIND) :: lat, lon, x, y
+
+ real (kind=RKIND) :: p0
+
+ real (kind=RKIND) :: etavs, ztemp
+
+ real (kind=RKIND) :: rs, rcv
+
+ ! calculation of the water vapor mixing ratio:
+ real (kind=RKIND) :: sh_max,sh_min,global_sh_max,global_sh_min
+
+ character (len=StrKIND), pointer :: config_met_prefix
+ logical, pointer :: config_use_spechumd
+ integer, pointer :: config_nfglevels
+ integer, pointer :: config_theta_adv_order
+ real (kind=RKIND), pointer :: config_coef_3rd_order
+
+ character (len=StrKIND), pointer :: config_extrap_airtemp
+ integer :: extrap_airtemp
+
+ real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
+ real (kind=RKIND), dimension(:), pointer :: latEdge, lonEdge
+ real (kind=RKIND), dimension(:), pointer :: angleEdge
+
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real (kind=RKIND), dimension(:,:), pointer :: w
+ real (kind=RKIND), dimension(:,:), pointer :: theta
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:), pointer :: relhum
+ real (kind=RKIND), dimension(:,:), pointer :: spechum
+ real (kind=RKIND), dimension(:,:), pointer :: ru
+ real (kind=RKIND), dimension(:,:), pointer :: rw
+
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructY
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructZ
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructZonal
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructMeridional
+ real (kind=RKIND), dimension(:,:), pointer :: u_fg
+ real (kind=RKIND), dimension(:,:), pointer :: v_fg
+ real (kind=RKIND), dimension(:,:), pointer :: z_fg
+ real (kind=RKIND), dimension(:,:), pointer :: t_fg
+ real (kind=RKIND), dimension(:,:), pointer :: rh_fg
+ real (kind=RKIND), dimension(:,:), pointer :: sh_fg
+ real (kind=RKIND), dimension(:,:), pointer :: p_fg
+ real (kind=RKIND), dimension(:), pointer :: soilz
+
+ type (hashtable), allocatable :: level_hash
+ logical :: too_many_fg_levs
+ integer :: level_value
+
+ character (len=StrKIND) :: errstring
+
+ real (kind=RKIND) :: max_zgrid_local, max_zgrid_global
+
+
+ call mpas_log_write('Interpolating FDDA input fields at time '//trim(timestamp))
+
+ call mpas_pool_get_config(configs, 'config_met_prefix', config_met_prefix)
+ call mpas_pool_get_config(configs, 'config_use_spechumd', config_use_spechumd)
+ call mpas_pool_get_config(configs, 'config_nfglevels', config_nfglevels)
+ call mpas_pool_get_config(configs, 'config_theta_adv_order', config_theta_adv_order)
+ call mpas_pool_get_config(configs, 'config_coef_3rd_order', config_coef_3rd_order)
+
+ call mpas_pool_get_config(configs, 'config_extrap_airtemp', config_extrap_airtemp)
+ if (trim(config_extrap_airtemp) == 'constant') then
+ extrap_airtemp = 0
+ else if (trim(config_extrap_airtemp) == 'linear') then
+ extrap_airtemp = 1
+ else if (trim(config_extrap_airtemp) == 'lapse-rate') then
+ extrap_airtemp = 2
+ else
+ call mpas_log_write('*************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('* Invalid value for namelist variable config_extrap_airtemp *', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('*************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+ call mpas_log_write("Using option '" // trim(config_extrap_airtemp) // "' for vertical extrapolation of temperature")
+
+ dminfo => block % domain % dminfo
+
+ call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
+ call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
+ call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
+ call mpas_pool_get_array(mesh, 'angleEdge', angleEdge)
+
+ call mpas_pool_get_array(mesh, 'zb', zb)
+ call mpas_pool_get_array(mesh, 'zb3', zb3)
+
+ call mpas_pool_get_array(mesh, 'zgrid', zgrid)
+ call mpas_pool_get_array(mesh, 'dzu', dzu)
+ call mpas_pool_get_array(mesh, 'fzm', fzm)
+ call mpas_pool_get_array(mesh, 'fzp', fzp)
+ call mpas_pool_get_array(mesh, 'zz', zz)
+
+ call mpas_pool_get_array(diag, 'exner_base', pb)
+ call mpas_pool_get_array(diag, 'rho_base', rb)
+ call mpas_pool_get_array(diag, 'theta_base', tb)
+ call mpas_pool_get_array(diag, 'rtheta_base', rtb)
+ call mpas_pool_get_array(diag, 'exner', p)
+ call mpas_pool_get_array(diag, 'pressure_base', ppb)
+ call mpas_pool_get_array(diag, 'pressure_p', pp)
+ call mpas_pool_get_array(diag, 'pressure', pressure)
+ call mpas_pool_get_array(diag, 'relhum', relhum)
+ call mpas_pool_get_array(diag, 'spechum', spechum)
+ call mpas_pool_get_array(diag, 'ru', ru)
+ call mpas_pool_get_array(diag, 'rw', rw)
+
+ call mpas_pool_get_array(state, 'rho_zz', rho_zz)
+ call mpas_pool_get_array(diag, 'rho_p', rr)
+ call mpas_pool_get_array(state, 'theta_m', t)
+ call mpas_pool_get_array(diag, 'rtheta_p', rt)
+ call mpas_pool_get_array(fdda_state, 'scalars_fdda', scalars)
+ call mpas_pool_get_array(fdda_state, 'uEdge_fdda', u)
+ call mpas_pool_get_array(fdda_state, 'theta_fdda', theta)
+ call mpas_pool_get_array(fdda_state, 'rho_fdda', rho)
+ call mpas_pool_get_array(fdda_state, 'uReconstructX_fdda', uReconstructX)
+ call mpas_pool_get_array(fdda_state, 'uReconstructY_fdda', uReconstructY)
+ call mpas_pool_get_array(fdda_state, 'uReconstructZ_fdda', uReconstructZ)
+ call mpas_pool_get_array(fdda_state, 'u_fdda', uReconstructZonal)
+ call mpas_pool_get_array(fdda_state, 'v_fdda', uReconstructMeridional)
+
+ call mpas_pool_get_array(mesh, 'latCell', latCell)
+ call mpas_pool_get_array(mesh, 'lonCell', lonCell)
+ call mpas_pool_get_array(mesh, 'latEdge', latEdge)
+ call mpas_pool_get_array(mesh, 'lonEdge', lonEdge)
+
+ call mpas_pool_get_array(fg, 'u', u_fg)
+ call mpas_pool_get_array(fg, 'v', v_fg)
+ call mpas_pool_get_array(fg, 'z', z_fg)
+ call mpas_pool_get_array(fg, 't', t_fg)
+ call mpas_pool_get_array(fg, 'rh', rh_fg)
+ call mpas_pool_get_array(fg, 'sh', sh_fg)
+ call mpas_pool_get_array(fg, 'p', p_fg)
+
+ call mpas_pool_get_dimension(dims, 'nCellsSolve', nCellsSolve)
+ nVertLevelsP1 = nVertLevels + 1
+
+ call mpas_pool_get_dimension(state, 'index_qv', index_qv)
+
+ etavs = (1.0_RKIND - 0.252_RKIND) * pii / 2.0_RKIND
+ rcv = rgas / (cp - rgas)
+ p0 = 1.0e+05_RKIND
+
+ scalars(:,:,:) = 0.0_RKIND
+
+ !
+ ! Check that we have what looks like a valid zgrid field. If the max value for zgrid is zero,
+ ! the input file likely does not contain vertical grid information.
+ !
+ max_zgrid_local = maxval(zgrid(:,1:nCellsSolve))
+ call mpas_dmpar_max_real(dminfo, max_zgrid_local, max_zgrid_global)
+ if (max_zgrid_global == 0.0_RKIND) then
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('The maximum value of the zgrid field is 0. Please ensure that the ''input'' stream ', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('contains valid vertical grid information.', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_CRIT)
+ end if
+
+
+ !
+ ! Horizontally interpolate meteorological data
+ !
+ allocate(vert_level(config_nfglevels))
+ vert_level(:) = -1.0
+
+ ! adding a check that the timestamp is actually of length >= 13
+ if ( len(trim(timestamp)) < 13 ) then
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('init_atm_case_fdda: Error in reading timestamp for processing meteorological data files', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_CRIT)
+ end if
+
+ call read_met_init(trim(config_met_prefix), .false., timestamp(1:13), istatus)
+
+ if (istatus /= 0) then
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Error opening initial meteorological data file '//trim(config_met_prefix)//':'//timestamp(1:13), &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('********************************************************************************', &
+ messageType=MPAS_LOG_CRIT)
+ end if
+
+ allocate(level_hash)
+ call mpas_hash_init(level_hash)
+ too_many_fg_levs = .false.
+
+ call read_next_met_field(field, istatus)
+
+ do while (istatus == 0)
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = SEARCH
+ interp_list(3) = 0
+
+ msgval = -1.e30
+
+ if (trim(field % field) == 'UU' .or. &
+ trim(field % field) == 'VV' .or. &
+ trim(field % field) == 'TT' .or. &
+ trim(field % field) == 'RH' .or. &
+ trim(field % field) == 'SPECHUMD' .or. &
+ trim(field % field) == 'GHT' .or. &
+ trim(field % field) == 'SOILHGT' .or. &
+ trim(field % field) == 'PRES' .or. &
+ trim(field % field) == 'PRESSURE') then
+
+ if (trim(field % field) /= 'SOILHGT') then
+
+ ! Since the hash table can only store integers, transfer the bit pattern from
+ ! the real-valued xlvl into an integer; that the result is not an integer version
+ ! of the level is not important, since we only want to test uniqueness of levels
+ level_value = transfer(field % xlvl, level_value)
+ if (.not. mpas_hash_search(level_hash, level_value)) then
+ call mpas_hash_insert(level_hash, level_value)
+ if (mpas_hash_size(level_hash) > config_nfglevels) then
+ too_many_fg_levs = .true.
+ end if
+ end if
+
+ !
+ ! In case we have more than config_nfglevels levels, just keep cycling through
+ ! the remaining fields in the intermediate file for the purpose of counting how
+ ! many unique levels are found using the code above
+ !
+ if (too_many_fg_levs) then
+ call read_next_met_field(field, istatus)
+ cycle
+ end if
+
+ do k=1,config_nfglevels
+ if (vert_level(k) == field % xlvl .or. vert_level(k) == -1.0) exit
+ end do
+ if (vert_level(k) == -1.0) vert_level(k) = field % xlvl
+ else
+ k = 1
+ end if
+
+ !
+ ! Set up projection
+ !
+ call map_init(proj)
+
+ if (field % iproj == PROJ_LATLON) then
+ call map_set(PROJ_LATLON, proj, &
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
+ else if (field % iproj == PROJ_GAUSS) then
+ call map_set(PROJ_GAUSS, proj, &
+ nlat = nint(field % deltalat), &
+ loninc = 360.0_RKIND / real(field % nx,RKIND), &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
+ end if
+
+
+ !
+ ! Horizontally interpolate the field at level k
+ !
+ if (trim(field % field) == 'UU') then
+ call mpas_log_write('Interpolating U at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nEdges
+ latPoints => latEdge
+ lonPoints => lonEdge
+ call mpas_pool_get_array(fg, 'u', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'VV') then
+ call mpas_log_write('Interpolating V at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nEdges
+ latPoints => latEdge
+ lonPoints => lonEdge
+ call mpas_pool_get_array(fg, 'v', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'TT') then
+ call mpas_log_write('Interpolating TT at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 't', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'RH') then
+ call mpas_log_write('Interpolating RH at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'rh', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'SPECHUMD') then
+ call mpas_log_write('Interpolating SPECHUMD at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'sh', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'GHT') then
+ call mpas_log_write('Interpolating GHT at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'z', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'PRES') then
+ call mpas_log_write('Interpolating PRES at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'p', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'PRESSURE') then
+ call mpas_log_write('Interpolating PRESSURE at $i $r', intArgs=(/k/), realArgs=(/vert_level(k)/))
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'p', destField2d)
+ ndims = 2
+ else if (trim(field % field) == 'SOILHGT') then
+ call mpas_log_write('Interpolating SOILHGT')
+ nInterpPoints = nCells
+ latPoints => latCell
+ lonPoints => lonCell
+ call mpas_pool_get_array(fg, 'soilz', destField1d)
+ ndims = 1
+ end if
+
+ allocate(rslab(-2:field % nx+3, field % ny))
+ rslab(1:field % nx, 1:field % ny) = field % slab(1:field % nx, 1:field % ny)
+ rslab(0, 1:field % ny) = field % slab(field % nx, 1:field % ny)
+ rslab(-1, 1:field % ny) = field % slab(field % nx-1, 1:field % ny)
+ rslab(-2, 1:field % ny) = field % slab(field % nx-2, 1:field % ny)
+ rslab(field % nx+1, 1:field % ny) = field % slab(1, 1:field % ny)
+ rslab(field % nx+2, 1:field % ny) = field % slab(2, 1:field % ny)
+ rslab(field % nx+3, 1:field % ny) = field % slab(3, 1:field % ny)
+
+ do i=1,nInterpPoints
+ lat = latPoints(i)*DEG_PER_RAD
+ lon = lonPoints(i)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ else if (x >= real(field%nx)+0.5) then
+ lon = lon - 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (y < 0.5) then
+ y = 1.0
+ else if (y >= real(field%ny)+0.5) then
+ y = real(field%ny)
+ end if
+ if (ndims == 1) then
+ destField1d(i) = interp_sequence(x, y, 1, rslab, -2, field%nx + 3, 1, field%ny, 1, 1, msgval, interp_list, 1)
+ else if (ndims == 2) then
+ destField2d(k,i) = interp_sequence(x, y, 1, rslab, -2, field%nx + 3, 1, field%ny, 1, 1, msgval, interp_list, 1)
+ end if
+ end do
+
+ deallocate(rslab)
+
+ end if
+
+ deallocate(field % slab)
+ call read_next_met_field(field, istatus)
+ end do
+
+ call read_met_close()
+ level_value = mpas_hash_size(level_hash)
+ call mpas_hash_destroy(level_hash)
+ deallocate(level_hash)
+
+ if (too_many_fg_levs) then
+ write(errstring,'(a,i4)') ' Please increase config_nfglevels to at least ', level_value
+ call mpas_log_write('*******************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Error: The meteorological data file has more than config_nfglevels.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write(trim(errstring), messageType=MPAS_LOG_ERR)
+ call mpas_log_write(' in the namelist and re-run.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('*******************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+
+ !
+ ! Check how many distinct levels we actually found in the meteorological data
+ !
+ do k=1,config_nfglevels
+ if (vert_level(k) == -1.0) exit
+ end do
+ nfglevels_actual = k-1
+ call mpas_log_write('*************************************************')
+ call mpas_log_write('Found $i levels in the first-guess data', intArgs=(/nfglevels_actual/))
+ call mpas_log_write('*************************************************')
+
+
+ !
+ ! For isobaric data, fill in the 3-d pressure field; otherwise, ensure
+ ! that the surface pressure and height fields are filled in
+ !
+ if (minval(p_fg(1:nfglevels_actual,1:nCellsSolve)) == 0.0 .and. &
+ maxval(p_fg(1:nfglevels_actual,1:nCellsSolve)) == 0.0) then
+ call mpas_log_write('Setting pressure field for isobaric data')
+ do k=1,config_nfglevels
+ if (vert_level(k) /= 200100.0) then
+ p_fg(k,:) = vert_level(k)
+ end if
+ end do
+ else
+ call mpas_pool_get_array(fg, 'z', z_fg)
+ call mpas_pool_get_array(fg, 'soilz', soilz)
+ call mpas_log_write('Assuming model-level input data')
+ do k=1,config_nfglevels
+ if (vert_level(k) == 200100.0) then
+ z_fg(k,:) = soilz(:)
+ end if
+ end do
+ end if
+
+
+ !
+ ! Compute normal wind component and store in fg % u
+ !
+ do iEdge=1,nEdges
+ do k=1,nfglevels_actual
+ u_fg(k,iEdge) = cos(angleEdge(iEdge)) * u_fg(k,iEdge) &
+ + sin(angleEdge(iEdge)) * v_fg(k,iEdge)
+ end do
+ end do
+
+
+ !
+ ! Vertically interpolate meteorological data
+ !
+ allocate(sorted_arr(2,nfglevels_actual))
+
+ do iCell=1,nCells
+
+ ! T
+ sorted_arr(:,:) = -999.0
+ do k = 1, nfglevels_actual
+ sorted_arr(1,k) = z_fg(k,iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = t_fg(k,iCell)
+ end do
+ call mpas_quicksort(nfglevels_actual, sorted_arr)
+ do k = 1, nVertLevels
+ target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
+ t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, &
+ extrap=extrap_airtemp, ierr=istatus)
+ if (istatus /= 0) then
+ write(errstring,'(a,i4,a,i10)') 'Error in interpolation of t(k,iCell) for k=', k, ', iCell=', iCell
+ call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write(trim(errstring), messageType=MPAS_LOG_ERR)
+ call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+ end do
+
+
+ ! RH
+ sorted_arr(:,:) = -999.0
+ relhum(:,iCell) = 0._RKIND
+ do k = 1, nfglevels_actual
+ sorted_arr(1,k) = z_fg(k,iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = rh_fg(k,iCell)
+ end do
+ call mpas_quicksort(nfglevels_actual, sorted_arr)
+ do k = nVertLevels, 1, -1
+ target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
+ relhum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ !bugfix-from-v8.3.1 if (target_z < z_fg(1,iCell) .and. k < nVertLevels) relhum(k,iCell) = relhum(k+1,iCell)
+ end do
+
+
+ ! SPECHUM: if first-guess values are negative, set those values to zero before
+ ! vertical interpolation.
+ sorted_arr(:,:) = -999.0
+ spechum(:,iCell) = 0._RKIND
+ do k = 1, nfglevels_actual
+ sorted_arr(1,k) = z_fg(k,iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = max(0._RKIND,sh_fg(k,iCell))
+ end do
+ call mpas_quicksort(nfglevels_actual, sorted_arr)
+ do k = nVertLevels, 1, -1
+ target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
+ spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ !bugfix-from-v8.3.1 if (target_z < z_fg(1,iCell) .and. k < nVertLevels) spechum(k,iCell) = spechum(k+1,iCell)
+ end do
+
+
+ ! PRESSURE
+ sorted_arr(:,:) = -999.0
+ do k = 1, nfglevels_actual
+ sorted_arr(1,k) = z_fg(k,iCell)
+ if (vert_level(k) == 200100.0) then
+ sorted_arr(1,k) = 99999.0
+ sfc_k = k
+ p_fg(k,iCell) = 1.0 ! Any value that has valid log is fine...
+ end if
+ sorted_arr(2,k) = log(p_fg(k,iCell))
+ end do
+ call mpas_quicksort(nfglevels_actual, sorted_arr)
+ do k = 1, nVertLevels
+ target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
+ pressure(k,iCell) = exp(vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1))
+ end do
+
+ end do
+
+
+ do iEdge=1,nEdges
+
+ ! U
+ sorted_arr(:,:) = -999.0
+ do k=1,nfglevels_actual
+ sorted_arr(1,k) = 0.5 * (z_fg(k,cellsOnEdge(1,iEdge)) + z_fg(k,cellsOnEdge(2,iEdge)))
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = u_fg(k,iEdge)
+ end do
+ call mpas_quicksort(nfglevels_actual, sorted_arr)
+ do k=1,nVertLevels
+ target_z = 0.25 * (zgrid(k,cellsOnEdge(1,iEdge)) + zgrid(k+1,cellsOnEdge(1,iEdge)) &
+ + zgrid(k,cellsOnEdge(2,iEdge)) + zgrid(k+1,cellsOnEdge(2,iEdge)))
+ u(k,iEdge) = vertical_interp(target_z, nfglevels_actual-1, sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1)
+ end do
+
+ end do
+
+
+ !
+ ! Reconstruct zonal and meridional winds for diagnostic puposes:
+ !
+ call mpas_rbf_interp_initialize(mesh)
+ call mpas_init_reconstruct(mesh)
+ call mpas_reconstruct(mesh, u, &
+ uReconstructX, &
+ uReconstructY, &
+ uReconstructZ, &
+ uReconstructZonal, &
+ uReconstructMeridional &
+ )
+
+
+ deallocate(sorted_arr)
+
+ ! Diagnose the water vapor mixing ratios:
+ global_sh_min = 0._RKIND
+ global_sh_max = 0._RKIND
+ if(config_use_spechumd) then
+ sh_min = minval(spechum(:,1:nCellsSolve))
+ sh_max = maxval(spechum(:,1:nCellsSolve))
+ call mpas_dmpar_min_real(dminfo,sh_min,global_sh_min)
+ call mpas_dmpar_max_real(dminfo,sh_max,global_sh_max)
+ endif
+ call mpas_log_write('')
+ call mpas_log_write('--- global_sh_min = $r', realArgs=(/global_sh_min/))
+ call mpas_log_write('--- global_sh_max = $r', realArgs=(/global_sh_max/))
+ call mpas_log_write('')
+
+ call mpas_log_write('--- config_use_spechumd = $l', logicArgs=(/config_use_spechumd/))
+ if(.not. config_use_spechumd .or. (global_sh_min==0._RKIND .and. global_sh_max==0._RKIND)) then
+ !--- calculate the saturation mixing ratio and interpolated first-guess relative humidity:
+ if (config_use_spechumd) then
+ call mpas_log_write('config_use_spechumd=T, but specific humidity was not found in ' &
+ //trim(config_met_prefix)//':'//timestamp(1:13), messageType=MPAS_LOG_WARN)
+ end if
+ call mpas_log_write(' *** initializing water vapor mixing ratio using first-guess relative humidity')
+ call mpas_log_write('')
+
+ do k = 1, nVertLevels
+ do iCell = 1, nCells
+ !
+ ! Note: the RH field provided by ungrib should always be with respect to liquid water,
+ ! hence, we can always call rslf; see the routine fix_gfs_rh in WPS/ungrib/src/rrpr.F .
+ !
+ rs = rslf(pressure(k,iCell),t(k,iCell))
+ scalars(index_qv,k,iCell) = 0.01_RKIND*rs*relhum(k,iCell)
+ enddo
+ enddo
+ else
+ !--- use the interpolated first-guess specific humidity:
+ call mpas_log_write(' *** initializing water vapor mixing ratio using first-guess specific humidity')
+ call mpas_log_write('')
+ do k = 1, nVertLevels
+ do iCell = 1, nCells
+ scalars(index_qv,k,iCell) = spechum(k,iCell)/(1._RKIND-spechum(k,iCell))
+ enddo
+ enddo
+ endif
+
+
+ !
+ ! Diagnose fields needed in initial conditions file (u, w, rho, theta)
+ ! NB: At this point, "rho_zz" is simple dry density, and "theta_m" is regular potential temperature
+ !
+ do iCell=1,nCells
+
+ do k=1,nVertLevels
+ ! PI
+ p(k,iCell) = (pressure(k,iCell) / p0) ** (rgas / cp)
+
+ ! THETA - can compute this using PI instead
+ t(k,iCell) = t(k,iCell) * (p0 / pressure(k,iCell)) ** (rgas / cp)
+
+ ! RHO_ZZ
+ ! MW note. this is from _gfs case and is different from _lbc case.
+ rho_zz(k,iCell) = pressure(k,iCell) / rgas / (p(k,iCell) * t(k,iCell) &
+ * (1.0 + (rvord - 1.0) * scalars(index_qv,k,iCell)))
+ rho_zz(k,iCell) = rho_zz(k,iCell) / (1.0 + scalars(index_qv,k,iCell))
+ end do
+ end do
+
+
+
+
+ !
+ ! Reference state based on a dry isothermal atmosphere
+ !
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ ztemp = 0.5*(zgrid(k+1,iCell)+zgrid(k,iCell))
+ ppb(k,iCell) = p0*exp(-gravity*ztemp/(rgas*t0b)) ! pressure_base
+ pb (k,iCell) = (ppb(k,iCell)/p0)**(rgas/cp) ! exner_base
+ rb (k,iCell) = ppb(k,iCell)/(rgas*t0b) ! rho_base
+ tb (k,iCell) = t0b/pb(k,iCell) ! theta_base
+ rtb(k,iCell) = rb(k,iCell)*tb(k,iCell) ! rtheta_base
+ p (k,iCell) = pb(k,iCell) ! exner
+ pp (k,iCell) = 0. ! pressure_p
+ rr (k,iCell) = 0. ! rho_p
+ end do
+ end do
+
+
+
+ do iCell=1,nCells
+ k = 1
+
+! WCS 20130821 - couple with vertical metric, note: rr is coupled here
+ rho_zz(k,iCell) = ((pressure(k,iCell) / p0)**(cv / cp)) * (p0 / rgas) &
+ / (t(k,iCell)*(1.0 + 1.61*scalars(index_qv,k,iCell))) / zz(k,iCell)
+ rr(k,iCell) = rho_zz(k,iCell) - rb(k,iCell)
+
+ do k=2,nVertLevels
+ it = 0
+ p_check = 2.0 * 0.0001
+ do while ( (it < 30) .and. (p_check > 0.0001) )
+
+ p_check = pp(k,iCell)
+
+! WCS 20130821 - MPAS hydrostatic relation
+ pp(k,iCell) = pp(k-1,iCell) - (fzm(k)*rr(k,iCell) + fzp(k)*rr(k-1,iCell))*gravity*dzu(k) &
+ - (fzm(k)*rho_zz(k,iCell)*scalars(index_qv,k,iCell) &
+ + fzp(k)*rho_zz(k-1,iCell)*scalars(index_qv,k-1,iCell))*gravity*dzu(k)
+ pressure(k,iCell) = pp(k,iCell) + ppb(k,iCell)
+ p(k,iCell) = (pressure(k,iCell) / p0) ** (rgas / cp)
+
+! WCS 20130821 - couple with vertical metric
+ rho_zz(k,iCell) = pressure(k,iCell) / rgas &
+ / (p(k,iCell)*t(k,iCell)*(1.0 + 1.61*scalars(index_qv,k,iCell)))/zz(k,iCell)
+ rr(k,iCell) = rho_zz(k,iCell) - rb(k,iCell)
+
+ p_check = abs(p_check - pp(k,iCell))
+
+ it = it + 1
+ end do
+ end do
+ end do
+
+
+
+ ! Compute theta_m and rho-tilde
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ t(k,iCell) = t(k,iCell) * (1.0 + 1.61*scalars(index_qv,k,iCell))
+! WCS 20130821 - decouple rr from vertical metric
+ rr(k,iCell) = rr(k,iCell)*zz(k,iCell)
+ end do
+ end do
+
+
+
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ ru(k,iEdge) = u(k,iEdge) * 0.5*(rho_zz(k,cellsOnEdge(1,iEdge)) + rho_zz(k,cellsOnEdge(2,iEdge)))
+ end do
+ end do
+
+
+
+
+ rw(:,:) = 0.0
+ do iCell=1,nCellsSolve
+
+ do i=1,nEdgesOnCell(iCell)
+ iEdge=edgesOnCell(i,iCell)
+
+ do k = 2, nVertLevels
+ flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ rw(k,iCell) = rw(k,iCell) - (fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))*zb(k,1,iEdge)*flux
+ else
+ rw(k,iCell) = rw(k,iCell) + (fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))*zb(k,2,iEdge)*flux
+ end if
+
+ if (config_theta_adv_order ==3) then
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ rw(k,iCell) = rw(k,iCell) &
+ + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order* &
+ (fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))*zb3(k,1,iEdge)*flux
+ else
+ rw(k,iCell) = rw(k,iCell) &
+ - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order* &
+ (fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))*zb3(k,2,iEdge)*flux
+ end if
+ end if
+
+ end do
+
+ end do
+
+ end do
+
+
+
+ deallocate(vert_level)
+
+
+
+
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell)
+ theta(k,iCell) = t(k,iCell) / (1.0_RKIND + 1.61_RKIND * scalars(index_qv,k,iCell))
+ end do
+ end do
+
+
+ end subroutine init_atm_case_fdda
+
+
+=======
+>>>>>>> 91c5eac175eebeaf4206bacd5cb50c39dff3c152
!-----------------------------------------------------------------------
! routine init_atm_case_cam_mpas
diff --git a/src/core_init_atmosphere/mpas_init_atm_core_interface.F b/src/core_init_atmosphere/mpas_init_atm_core_interface.F
index f277a4a72f..f66a3ac2e8 100644
--- a/src/core_init_atmosphere/mpas_init_atm_core_interface.F
+++ b/src/core_init_atmosphere/mpas_init_atm_core_interface.F
@@ -114,7 +114,7 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
integer :: ierr
logical :: lexist
- logical, pointer :: initial_conds, sfc_update, lbcs
+ logical, pointer :: initial_conds, sfc_update, lbcs, fdda
logical, pointer :: gwd_stage_in, gwd_stage_out, vertical_stage_in, vertical_stage_out, met_stage_in, met_stage_out
logical, pointer :: gwd_gsl_stage_out
logical, pointer :: config_native_gwd_static, config_static_interp, config_vertical_grid, config_met_interp
@@ -146,6 +146,9 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
nullify(lbcs)
call mpas_pool_get_package(packages, 'lbcsActive', lbcs)
+ nullify(fdda)
+ call mpas_pool_get_package(packages, 'fddaActive', fdda)
+
nullify(gwd_stage_in)
call mpas_pool_get_package(packages, 'gwd_stage_inActive', gwd_stage_in)
@@ -205,6 +208,13 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
mp_thompson_aers_in = .false.
end if
+ if (config_init_case == 11) then
+ fdda = .true.
+ else
+ fdda = .false.
+ end if
+
+
if (config_init_case == 7) then
!
@@ -256,6 +266,19 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
inquire(file="QNWFA_QNIFA_SIGMA_MONTHLY.dat",exist=lexist)
if((lexist .and. met_stage_out) .or. (lexist .and. met_stage_in)) mp_thompson_aers_in = .true.
+ initial_conds = .false. ! Also, turn off the initial_conds package to avoid writing the IC "output" stream
+ !
+ ! When interpolating FDDA input fields, we need all inputs that would be needed for the interpolation
+ ! of ICs (similar to LBC), so met_stage_in = .true.
+ !
+ else if (config_init_case == 11) then
+ gwd_stage_in = .false.
+ gwd_stage_out = .false.
+ vertical_stage_in = .false.
+ vertical_stage_out = .false.
+ met_stage_in = .true.
+ met_stage_out = .true.
+
initial_conds = .false. ! Also, turn off the initial_conds package to avoid writing the IC "output" stream
else if (config_init_case == 13) then
@@ -285,7 +308,7 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
call mpas_pool_get_package(packages, 'first_guess_fieldActive', first_guess_field)
first_guess_field = .false.
- if ((config_init_case == 7 .and. config_met_interp) .or. config_init_case == 9) then
+ if ((config_init_case == 7 .and. config_met_interp) .or. config_init_case == 9 .or. config_init_case == 11) then
first_guess_field = .true.
end if