!wrf:model_layer:physics
!
!ckay=KIRAN ALAPATY @ US EPA -- November 01, 2015
!
! Tim Glotfelty@CNSU; AJ Deng@PSU
!modified for use with FASDAS
!Flux Adjusting Surface Data Assimilation System to assimilate
!surface layer and soil layers temperature and moisture using
! surfance reanalsys
!Reference: Alapaty et al., 2008: Development of the flux-adjusting surface
! data assimilation system for mesoscale models. JAMC, 47, 2331-2350
!
!
MODULE module_fdda_psufddagd 2
CONTAINS
!
!-------------------------------------------------------------------
!
SUBROUTINE fddagd(itimestep,dx,dt,xtime, & 1,12
id,analysis_interval, end_fdda_hour, &
if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
guv, gt, gq, if_ramping, dtramp_min, &
grid_sfdda, &
analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
rinblw, &
u3d,v3d,th3d,t3d, &
qv3d, &
p3d,pi3d, &
u_ndg_old,v_ndg_old,t_ndg_old,q_ndg_old,mu_ndg_old, &
u_ndg_new,v_ndg_new,t_ndg_new,q_ndg_new,mu_ndg_new, &
u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
fasdas, SDA_HFX, SDA_QFX, & ! fasdas
HFX_FDDA,dz8w, & ! fasdas
pblh, ht, regime, znt, z, z_at_w, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!-------------------------------------------------------------------
implicit none
!-------------------------------------------------------------------
!
! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
! and released in December 2006. The surface analysis nudging capability was added and
! released in March 2009 with WRFV3.1.
!
!-- u3d 3d u-velocity staggered on u points
!-- v3d 3d v-velocity staggered on v points
!-- th3d 3d potential temperature (k)
!-- t3d temperature (k)
!-- qv3d 3d water vapor mixing ratio (kg/kg)
!-- p3d 3d pressure (pa)
!-- pi3d 3d exner function (dimensionless)
!-- rundgdten staggered u tendency due to
! fdda grid nudging (m/s/s)
!-- rvndgdten staggered v tendency due to
! fdda grid nudging (m/s/s)
!-- rthndgdten theta tendency due to
! fdda grid nudging (K/s)
!-- rqvndgdten qv tendency due to
! fdda grid nudging (kg/kg/s)
!-- rmundgdten mu tendency due to
! fdda grid nudging (Pa/s)
!-- ids start index for i in domain
!-- ide end index for i in domain
!-- jds start index for j in domain
!-- jde end index for j in domain
!-- kds start index for k in domain
!-- kde end index for k in domain
!-- ims start index for i in memory
!-- ime end index for i in memory
!-- jms start index for j in memory
!-- jme end index for j in memory
!-- kms start index for k in memory
!-- kme end index for k in memory
!-- its start index for i in tile
!-- ite end index for i in tile
!-- jts start index for j in tile
!-- jte end index for j in tile
!-- kts start index for k in tile
!-- kte end index for k in tile
!-------------------------------------------------------------------
!
INTEGER, INTENT(IN) :: itimestep, analysis_interval, end_fdda_hour
INTEGER, INTENT(IN) :: analysis_interval_sfc, end_fdda_hour_sfc
INTEGER, INTENT(IN) :: grid_sfdda
INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
if_no_pbl_nudging_q
INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
INTEGER, INTENT(IN) :: if_ramping
INTEGER , INTENT(IN) :: id
REAL, INTENT(IN) :: DT, dx, xtime, dtramp_min
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: qv3d, &
p3d, &
pi3d, &
th3d, &
t3d, &
z, &
z_at_w
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT) :: rundgdten, &
rvndgdten, &
rthndgdten, &
rqvndgdten
!
! FASDAS
!
INTEGER, INTENT(IN) :: fasdas
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: SDA_HFX, &
SDA_QFX
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT) :: HFX_FDDA
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: dz8w
!
! END FASDAS
!
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: rmundgdten
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: u_ndg_old, &
v_ndg_old, &
t_ndg_old, &
q_ndg_old, &
u_ndg_new, &
v_ndg_new, &
t_ndg_new, &
q_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: u10_ndg_old, &
v10_ndg_old, &
t2_ndg_old, &
th2_ndg_old, &
q2_ndg_old, &
rh_ndg_old, &
psl_ndg_old, &
ps_ndg_old, &
u10_ndg_new, &
v10_ndg_new, &
t2_ndg_new, &
th2_ndg_new, &
q2_ndg_new, &
rh_ndg_new, &
psl_ndg_new, &
ps_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: tob_ndg_old, &
tob_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: mu_ndg_old, &
mu_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: odis_ndg_old, odis_ndg_new
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: u3d, &
v3d
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
ht, &
regime, &
znt
REAL, INTENT(IN) :: guv, gt, gq
REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
REAL :: xtime_old, xtime_new, coef, val_analysis
INTEGER :: kpbl, dbg_level
REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
!BPR BEGIN
REAL :: tfac_sfc, actual_end_fdda_min_sfc
!BPR END
REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ) :: wpbl ! 1: u, 2: v, 3: t, 4: q
REAL, DIMENSION( kts:kte, 4 ) :: wzfac ! 1: u, 2: v, 3: t, 4: q
! TWG 2015 Pseudo Radiative Flux
REAL, DIMENSION( its:ite, kts:kte, jts:jte) :: rho_air
! END TWG 2015
LOGICAL , EXTERNAL :: wrf_dm_on_monitor
CHARACTER (LEN=256) :: message
INTEGER :: int4
int4 = 1 ! 1: temporal interpolation. else: target nudging toward *_ndg_new values
actual_end_fdda_min = end_fdda_hour*60.0
!BPR BEGIN
actual_end_fdda_min_sfc = end_fdda_hour*60.0
!BPR END
IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
!BPR BEGIN
actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0 + ABS(dtramp_min)
!BPR END
!BPR BEGIN
!IF( xtime > actual_end_fdda_min ) THEN
IF( ( xtime > actual_end_fdda_min ) .AND. ( xtime > actual_end_fdda_min_sfc ) ) THEN
!BPR END
! If xtime is greater than the end time, no need to calculate tendencies. Just set the tendencies
! to zero to turn off nudging and return.
DO j = jts, jte
DO k = kts, kte
DO i = its, ite
RUNDGDTEN(i,k,j) = 0.0
RVNDGDTEN(i,k,j) = 0.0
RTHNDGDTEN(i,k,j) = 0.0
RQVNDGDTEN(i,k,j) = 0.0
IF( k .EQ. kts ) RMUNDGDTEN(i,j) = 0.
ENDDO
ENDDO
ENDDO
RETURN
ENDIF
IF( analysis_interval <= 0 )CALL wrf_error_fatal
('In grid FDDA, gfdda_interval_m must be > 0')
xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
xtime_new = xtime_old + analysis_interval
IF( int4 == 1 ) THEN
coef = (xtime-xtime_old)/(xtime_new-xtime_old)
ELSE
coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
ENDIF
IF ( wrf_dm_on_monitor()) THEN
CALL get_wrf_debug_level( dbg_level )
IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
IF( xtime < end_fdda_hour*60.0 ) THEN
WRITE(message,'(a,i1,a,f10.3,a)') &
'D0',id,' 3-D analysis nudging reads new data at time = ', xtime, ' min.'
CALL wrf_message( TRIM(message) )
WRITE(message,'(a,i1,a,2f8.2,a)') &
'D0',id,' 3-D analysis nudging bracketing times = ', xtime_old, xtime_new, ' min.'
CALL wrf_message( TRIM(message) )
ENDIF
actual_end_fdda_min = end_fdda_hour*60.0
IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
! Find the mid point of the tile and print out the sample values
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
IF( guv > 0.0 ) THEN
DO k = kts, kte
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
CALL wrf_message( TRIM(message) )
ENDDO
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
' mu_ndg_old=', mu_ndg_old(i0,j0), ' mu_ndg_new=', mu_ndg_new(i0,j0)
CALL wrf_message( TRIM(message) )
DO k = kts, kte
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
CALL wrf_message( TRIM(message) )
ENDDO
ENDIF
IF( gt > 0.0 ) THEN
DO k = kts, kte
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
CALL wrf_message( TRIM(message) )
ENDDO
ENDIF
IF( gq > 0.0 ) THEN
DO k = kts, kte
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
' q_ndg_old=', q_ndg_old(i0,k,j0), ' q_ndg_new=', q_ndg_new(i0,k,j0)
CALL wrf_message( TRIM(message) )
ENDDO
ENDIF
IF( int4 == 1 ) then
WRITE(message,'(a,i1,a)') ' D0',id, &
' 3-D nudging towards the temporally interpolated analysis'
ELSE
WRITE(message,'(a,i1,a)') ' D0',id, &
' 3-D nudging towards the target analysis'
ENDIF
ENDIF
ENDIF
ENDIF
jtsv=MAX0(jts,jds+1)
itsu=MAX0(its,ids+1)
jtf=MIN0(jte,jde-1)
ktf=MIN0(kte,kde-1)
itf=MIN0(ite,ide-1)
!
! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t,
! if_no_pbl_nudging_q switches) are set to 1, compute the weighting function, wpbl(:,k,:,:),
! based on the PBL depth. wpbl = 1 above the PBL and wpbl = 0 in the PBL. If all
! the switches are set to zero, wpbl = 1 (default value).
!
wpbl(:,:,:,:) = 1.0
IF( if_no_pbl_nudging_uv == 1 .OR. grid_sfdda >= 1 ) THEN
DO j=jts,jtf
DO i=itsu,itf
kpbl = 1
zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
loop_ku: DO k=kts,ktf
! zagl = 0.5 * ( z(i-1,k,j)-ht(i-1,j) + z(i,k,j)-ht(i,j) )
zagl_bot = 0.5 * ( z_at_w(i-1,k, j)-ht(i-1,j) + z_at_w(i,k, j)-ht(i,j) )
zagl_top = 0.5 * ( z_at_w(i-1,k+1,j)-ht(i-1,j) + z_at_w(i,k+1,j)-ht(i,j) )
IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
kpbl = k
EXIT loop_ku
ENDIF
ENDDO loop_ku
DO k=kts,ktf
IF( k <= kpbl ) wpbl(i, k, j, 1) = 0.0
IF( k == kpbl+1 ) wpbl(i, k, j, 1) = 0.1
IF( k > kpbl+1 ) wpbl(i, k, j, 1) = 1.0
ENDDO
ENDDO
ENDDO
DO i=its,itf
DO j=jtsv,jtf
kpbl = 1
zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
loop_kv: DO k=kts,ktf
! zagl = 0.5 * ( z(i,k,j-1)-ht(i,j-1) + z(i,k,j)-ht(i,j) )
zagl_bot = 0.5 * ( z_at_w(i,k, j-1)-ht(i,j-1) + z_at_w(i,k, j)-ht(i,j) )
zagl_top = 0.5 * ( z_at_w(i,k+1,j-1)-ht(i,j-1) + z_at_w(i,k+1,j)-ht(i,j) )
IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
kpbl = k
EXIT loop_kv
ENDIF
ENDDO loop_kv
DO k=kts,ktf
IF( k <= kpbl ) wpbl(i, k, j, 2) = 0.0
IF( k == kpbl+1 ) wpbl(i, k, j, 2) = 0.1
IF( k > kpbl+1 ) wpbl(i, k, j, 2) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
IF( if_no_pbl_nudging_t == 1 .OR. grid_sfdda >= 1 ) THEN
DO j=jts,jtf
DO i=its,itf
kpbl = 1
zpbl = pblh(i,j)
loop_kt: DO k=kts,ktf
! zagl = z(i,k,j)-ht(i,j)
zagl_bot = z_at_w(i,k, j)-ht(i,j)
zagl_top = z_at_w(i,k+1,j)-ht(i,j)
IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
kpbl = k
EXIT loop_kt
ENDIF
ENDDO loop_kt
DO k=kts,ktf
IF( k <= kpbl ) wpbl(i, k, j, 3) = 0.0
IF( k == kpbl+1 ) wpbl(i, k, j, 3) = 0.1
IF( k > kpbl+1 ) wpbl(i, k, j, 3) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
IF( if_no_pbl_nudging_q == 1 .OR. grid_sfdda >= 1 ) THEN
DO j=jts,jtf
DO i=its,itf
kpbl = 1
zpbl = pblh(i,j)
loop_kq: DO k=kts,ktf
! zagl = z(i,k,j)-ht(i,j)
zagl_bot = z_at_w(i,k, j)-ht(i,j)
zagl_top = z_at_w(i,k+1,j)-ht(i,j)
IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
kpbl = k
EXIT loop_kq
ENDIF
ENDDO loop_kq
DO k=kts,ktf
IF( k <= kpbl ) wpbl(i, k, j, 4) = 0.0
IF( k == kpbl+1 ) wpbl(i, k, j, 4) = 0.1
IF( k > kpbl+1 ) wpbl(i, k, j, 4) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
!
! If the user-defined namelist switches (if_zfac_uv, if_zfac_t,
! if_zfac_q swithes) are set to 1, compute the weighting function, wzfac(k,:),
! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_q) below which analysis
! nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x). If all
! the switche are set to zero, wzfac = 1 (default value).
!
wzfac(:,:) = 1.0
IF( if_zfac_uv == 1 ) THEN
DO j=jts,jtf
DO i=itsu,itf
DO k=kts,ktf
IF( k <= k_zfac_uv ) wzfac(k, 1:2) = 0.0
IF( k == k_zfac_uv+1 ) wzfac(k, 1:2) = 0.1
IF( k > k_zfac_uv+1 ) wzfac(k, 1:2) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
IF( if_zfac_t == 1 ) THEN
DO j=jts,jtf
DO i=itsu,itf
DO k=kts,ktf
IF( k <= k_zfac_t ) wzfac(k, 3) = 0.0
IF( k == k_zfac_t+1 ) wzfac(k, 3) = 0.1
IF( k > k_zfac_t+1 ) wzfac(k, 3) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
IF( if_zfac_q == 1 ) THEN
DO j=jts,jtf
DO i=itsu,itf
DO k=kts,ktf
IF( k <= k_zfac_q ) wzfac(k, 4) = 0.0
IF( k == k_zfac_q+1 ) wzfac(k, 4) = 0.1
IF( k > k_zfac_q+1 ) wzfac(k, 4) = 1.0
ENDDO
ENDDO
ENDDO
ENDIF
!
! If if_ramping and dtramp_min are defined by user, compute a time weighting function, tfac,
! for analysis nudging so that at the end of the nudging period (which has to be at a
! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
!
! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at
! end_fdda_hour-ABS(dtramp_min).
!
! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at
! end_fdda_hour+ABS(dtramp_min).
! BPR BEGIN
! In this case, during the rampdown we nudge towards the most recent past analysis and ignore
! the future analysis (implemented by setting coef = 0.0)
! THE FOLLOWING COMMENT IS NO LONGER APPLICABLE
! In this case, the obs values are extrapolated using
! the obs tendency saved from the previous FDDA window. More specifically for extrapolation,
! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
! THE PRECEDING COMMENT IS NOT LONGER APPLICABLE
! BPR END
!
tfac = 1.0
!BPR BEGIN
tfac_sfc = 1.0
!BPR END
IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
IF( dtramp_min <= 0.0 ) THEN
actual_end_fdda_min = end_fdda_hour*60.0
ELSE
actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
ENDIF
IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN
tfac = 1.0
ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
!BPR BEGIN
!The original method assumed that to be here in the code with dtramp_min>0
!meant that we were after the valid time of *_ndg_new and so used the
!values *_ndg_old and *_ndg_new to extrapolate a current analysis value.
!HOWEVER, in practice once we get to the valid time of *_ndg_new WRF will
!read in the next pair of *_ndg_old and *_ndg_new values, even if we are
!at the beginning of the rampdown period. Since the *_ndg_new values have
!a valid time after the beginning of the rampdown period we do not want
!to use them as we would be using analyses valid after the supposed end
!time of the analysis nudging.
!IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/(analysis_interval*1.0)
!Use only the *_ndg_old values
IF( dtramp_min > 0.0 ) coef = 0.0
!BPR END
ELSE
tfac = 0.0
ENDIF
!BPR BEGIN
!Now calculate the same quantities for surface analysis nudging
!Add actual_end_fdda_min_sfc, tfac_sfc
IF( dtramp_min <= 0.0 ) THEN
actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0
ELSE
actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0 + dtramp_min
ENDIF
IF( xtime < actual_end_fdda_min_sfc-ABS(dtramp_min) )THEN
tfac_sfc = 1.0
ELSEIF( xtime >= actual_end_fdda_min_sfc-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min_sfc )THEN
tfac_sfc = ( actual_end_fdda_min_sfc - xtime ) / ABS(dtramp_min)
ELSE
tfac_sfc = 0.0
ENDIF
!BPR END
ENDIF
!
! Surface Analysis Nudging
!
IF( grid_sfdda >= 1 ) THEN
CALL SFDDAGD(itimestep,dx,dt,xtime, id, &
analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
rinblw, &
u3d,v3d,th3d,t3d, &
qv3d, &
p3d,pi3d, &
u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
pblh, ht, regime, znt, z, z_at_w, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
!BPR BEGIN
! actual_end_fdda_min, tfac &
actual_end_fdda_min_sfc, tfac_sfc &
!BPR END
!
! FASDAS
!
,fasdas, SDA_HFX, SDA_QFX &
!
! END FASDAS
!
)
ENDIF
!
! Compute 3-D nudging tendencies for u, v, t and q
!
DO j=jts,jtf
DO k=kts,ktf
DO i=itsu,itf
val_analysis = u_ndg_old(i,k,j) *( 1.0 - coef ) + u_ndg_new(i,k,j) * coef
RUNDGDTEN(i,k,j) = RUNDGDTEN(i,k,j) + guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
( val_analysis - u3d(i,k,j) )
ENDDO
ENDDO
ENDDO
DO j=jtsv,jtf
DO k=kts,ktf
DO i=its,itf
val_analysis = v_ndg_old(i,k,j) *( 1.0 - coef ) + v_ndg_new(i,k,j) * coef
RVNDGDTEN(i,k,j) = RVNDGDTEN(i,k,j) + guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
( val_analysis - v3d(i,k,j) )
ENDDO
ENDDO
ENDDO
DO j=jts,jtf
DO k=kts,ktf
DO i=its,itf
val_analysis = t_ndg_old(i,k,j) *( 1.0 - coef ) + t_ndg_new(i,k,j) * coef
RTHNDGDTEN(i,k,j) = RTHNDGDTEN(i,k,j) + gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
( val_analysis - th3d(i,k,j) + 300.0 )
!
!FASDAS
!
!TWG 2015 Pseudo Radiative Flux
rho_air(i,k,j) = p3d(i,k,j)/(287.0*th3d(i,k,j))
HFX_FDDA(i,k,j) = rho_air(i,k,j)*1004.0*RTHNDGDTEN(i,k,j)*dz8w(i,k,j)
!TWG 2015 END
!
!END FASDAS
!
val_analysis = q_ndg_old(i,k,j) *( 1.0 - coef ) + q_ndg_new(i,k,j) * coef
RQVNDGDTEN(i,k,j) = RQVNDGDTEN(i,k,j) + gq * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
( val_analysis - qv3d(i,k,j) )
ENDDO
ENDDO
ENDDO
!BPR BEGIN
!Diagnostic print
IF ( wrf_dm_on_monitor()) THEN
IF( dbg_level .GE. 10 ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
k0 = (kte-kts)/2+kts
WRITE(message,'(a,i1,a,f7.2,a,i4,a,i4,a,i4,a,f7.2,a,f4.2,a,f7.2,a,f4.2,a,f4.2 )') &
' D0',id,' At xtime=',xtime,' FDDA sample pot. temp. analysis (i=',i0,', j=',j0,' k=', &
k0,') = (t_ndg_old=', t_ndg_old(i0,k0,j0),') * ',1.0-coef,' + (t_ndg_new=', &
t_ndg_new(i0,k0,j0),') * ',coef, ' where tfac=',tfac
CALL wrf_message( TRIM(message) )
WRITE(message,'(a,i1,a,f7.2,a,i4,a,f7.2,a,f7.2)') &
' D0',id,' xtime_old=',xtime_old,' analysis_interval=',analysis_interval,' actual_end_fdda_min=', &
actual_end_fdda_min,' dtramp_min=',dtramp_min
CALL wrf_message( TRIM(message) )
END IF
END IF
!BPR END
END SUBROUTINE fddagd
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE sfddagd(itimestep,dx,dt,xtime, & 1,21
id, analysis_interval_sfc, end_fdda_hour_sfc, &
guv_sfc, gt_sfc, gq_sfc, rinblw, &
u3d,v3d,th3d,t3d, &
qv3d, &
p3d,pi3d, &
u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN, &
pblh, ht, regime, znt, z, z_at_w, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
!BPR BEGIN
! actual_end_fdda_min, tfac &
actual_end_fdda_min_sfc, tfac_sfc &
!BPR END
!
! FASDAS
!
,fasdas, SDA_HFX, SDA_QFX &
!
! END FASDAS
!
)
!-------------------------------------------------------------------
USE module_model_constants
implicit none
!-------------------------------------------------------------------
!
! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
! and released in December 2006. The surface analysis nudging capability was added and
! released in March 2009 with WRFV3.1.
!
!-- u3d 3d u-velocity staggered on u points
!-- v3d 3d v-velocity staggered on v points
!-- th3d 3d potential temperature (k)
!-- t3d temperature (k)
!-- qv3d 3d water vapor mixing ratio (kg/kg)
!-- p3d 3d pressure (pa)
!-- pi3d 3d exner function (dimensionless)
!-- rundgdten staggered u tendency due to
! fdda grid nudging (m/s/s)
!-- rvndgdten staggered v tendency due to
! fdda grid nudging (m/s/s)
!-- rthndgdten theta tendency due to
! fdda grid nudging (K/s)
!-- rqvndgdten qv tendency due to
! fdda grid nudging (kg/kg/s)
!-- rmundgdten mu tendency due to
! fdda grid nudging (Pa/s)
!-- ids start index for i in domain
!-- ide end index for i in domain
!-- jds start index for j in domain
!-- jde end index for j in domain
!-- kds start index for k in domain
!-- kde end index for k in domain
!-- ims start index for i in memory
!-- ime end index for i in memory
!-- jms start index for j in memory
!-- jme end index for j in memory
!-- kms start index for k in memory
!-- kme end index for k in memory
!-- its start index for i in tile
!-- ite end index for i in tile
!-- jts start index for j in tile
!-- jte end index for j in tile
!-- kts start index for k in tile
!-- kte end index for k in tile
!-------------------------------------------------------------------
!
INTEGER, INTENT(IN) :: itimestep, analysis_interval_sfc, end_fdda_hour_sfc
INTEGER , INTENT(IN) :: id
REAL, INTENT(IN) :: dx,DT, xtime
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: qv3d, &
p3d, &
pi3d, &
th3d, &
t3d, &
z, &
z_at_w
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT) :: rundgdten, &
rvndgdten, &
rthndgdten, &
rqvndgdten
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: rmundgdten
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: u10_ndg_old, &
v10_ndg_old, &
t2_ndg_old, &
th2_ndg_old, &
q2_ndg_old, &
rh_ndg_old, &
psl_ndg_old, &
ps_ndg_old, &
u10_ndg_new, &
v10_ndg_new, &
t2_ndg_new, &
th2_ndg_new, &
q2_ndg_new, &
rh_ndg_new, &
psl_ndg_new, &
ps_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: tob_ndg_old, &
tob_ndg_new
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: u3d, &
v3d
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: odis_ndg_old, odis_ndg_new
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
ht, &
regime, &
znt
REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, j0
REAL :: xtime_old_sfc, xtime_new_sfc, coef, val_analysis, es
INTEGER :: kpbl, dbg_level
REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac_sfc, actual_end_fdda_min_sfc
REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ), &
INTENT(INout) :: wpbl ! 1: u, 2: v, 3: t, 4: q
REAL, DIMENSION( kts:kte, 4 ), &
INTENT(IN) :: wzfac ! 1: u, 2: v, 3: t, 4: q
REAL, DIMENSION( its:ite, jts:jte) :: wndcor_u, wndcor_v
REAL, DIMENSION( its-2:ite+2, jts-2:jte+2) :: blw_old, blw_new
REAL, DIMENSION( its:ite, kts:kte, jts:jte) :: qsat
REAL :: m, b=1.8, blw, rindx, x
REAL :: difz, wr14, wr1z, wr24, wr2z, wndfac, reg, znt0
INTEGER, INTENT(IN) :: if_ramping
REAL, INTENT(IN) :: dtramp_min
LOGICAL , EXTERNAL :: wrf_dm_on_monitor
CHARACTER (LEN=256) :: message
INTEGER :: iwinds, idd, iqsat, int4
!
! FASDAS
!
INTEGER, INTENT(IN) :: fasdas
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: SDA_HFX, &
SDA_QFX
REAL :: stabFac, exf, DiffTN
REAL, DIMENSION( its:ite, kts:kte, jts:jte ) :: wpblfasdas
!
! END FASDAS
!
iwinds = 1 ! 1: Scale the surface wind analysis to the lowest model level,
! if the first model half-layer is greater than 10 meters
! and we are in the free convection regime (REGIME=4.0). else: no
idd = 1 ! 1: Obs data density correction is applied. else: no
iqsat = 1 ! 1: Remove super saturation. eles: no
int4 = 1 ! 1: temporal ionterpolation. else: target nudging toward *_ndg_new values
IF( analysis_interval_sfc <= 0 )CALL wrf_error_fatal('In grid sfc FDDA, sgfdda_interval_m must be > 0')
xtime_old_sfc = FLOOR(xtime/analysis_interval_sfc) * analysis_interval_sfc * 1.0
xtime_new_sfc = xtime_old_sfc + analysis_interval_sfc
IF( int4 == 1 ) THEN
coef = (xtime-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc) ! Temporal interpolation
ELSE
coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
ENDIF
IF ( wrf_dm_on_monitor()) THEN
CALL get_wrf_debug_level( dbg_level )
IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
IF( xtime < end_fdda_hour_sfc*60.0 ) THEN
WRITE(message,'(a,i1,a,f10.3,a)') &
'D0',id,' surface analysis nudging reads new data at time = ', xtime, ' min.'
CALL wrf_message( TRIM(message) )
WRITE(message,'(a,i1,a,2f8.2,a)') &
'D0',id,' surface analysis nudging bracketing times = ', xtime_old_sfc, xtime_new_sfc, ' min.'
CALL wrf_message( TRIM(message) )
ENDIF
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
! Find the mid point of the tile and print out the sample values
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
IF( guv_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
' D0',id,' sample surface analysis values at i,j=', i0, j0, &
' u10_ndg_old=', u10_ndg_old(i0,j0), ' u10_ndg_new=', u10_ndg_new(i0,j0)
CALL wrf_message( TRIM(message) )
WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
' D0',id,' sample surface analysis values at i,j=', i0, j0, &
' v10_ndg_old=', v10_ndg_old(i0,j0), ' v10_ndg_new=', v10_ndg_new(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
IF( gt_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
' D0',id,' sample surface analysis values at i,j=', i0, j0, &
' th2_ndg_old=', th2_ndg_old(i0,j0), ' th2_ndg_new=', th2_ndg_new(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
IF( gq_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
' D0',id,' sample surface analysis values at i,j=', i0, j0, &
' q2_ndg_old=', q2_ndg_old(i0,j0), ' q2_ndg_new=', q2_ndg_new(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
IF( iwinds == 1 ) &
WRITE(message,'(a,i1,a)') ' D0',id, &
' surface wind analysis s scaled to the lowest model level, if dz1 > 10m and REGIME=4.'
IF( idd == 1 ) &
WRITE(message,'(a,i1,a)') ' D0',id, &
' obs data density is used for additional weighting function'
IF( iqsat == 1 ) &
WRITE(message,'(a,i1,a)') ' D0',id, &
' super saturation is not allowed for q analysis'
IF( int4 == 1 ) then
WRITE(message,'(a,i1,a)') ' D0',id, &
' surface nudging towards the temporally interpolated analysis'
ELSE
WRITE(message,'(a,i1,a)') ' D0',id, &
' surface nudging towards the target analysis'
ENDIF
ENDIF
ENDIF
ENDIF
jtsv=MAX0(jts,jds+1)
itsu=MAX0(its,ids+1)
jtf=MIN0(jte,jde-1)
ktf=MIN0(kte,kde-1)
itf=MIN0(ite,ide-1)
!
! Compute the vertical weighting function to scale the surface wind analysis to
! the lowest model level, if the first model half-layer is greater
! than 10 meters and we are in the free convection regime (REGIME=4.0).
!
IF( iwinds == 1 ) THEN
wndcor_u(:,:) = 1.0
DO j=jts,jtf
DO i=itsu,itf
reg = 0.5 * ( regime(i-1, j) + regime(i, j) )
difz = 0.5 * ( z(i-1,1,j) - ht(i-1,j) &
+ z(i, 1,j) - ht(i, j) )
IF( reg > 3.5 .AND. difz > 10.0 ) THEN
znt0 = 0.5 * ( znt(i-1, j) + znt(i, j) )
IF( znt0 <= 0.2) THEN
wndcor_u(i,j) = 1.0+0.320*znt0**0.2
ELSE
wndcor_u(i,j) = 1.169+0.315*znt0
ENDIF
wr14 = log(40.0/0.05)
wr1z = log(difz/0.05)
wr24 = log(40.0/1.0)
wr2z = log(difz/1.0)
wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
wndcor_u(i,j) = wndfac*wndcor_u(i,j)
ENDIF
ENDDO
ENDDO
IF ( wrf_dm_on_monitor()) THEN
IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
' D0',id,' sample wndcor_u values at i,j=', i0, j0, &
' wndcor_u=', wndcor_u(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
ENDIF
ENDIF
ELSE
wndcor_u(:,:) = 1.0
ENDIF
IF( iwinds == 1 ) THEN
wndcor_v(:,:) = 1.0
DO j=jtsv,jtf
DO i=its,itf
reg = 0.5 * ( regime(i, j-1) + regime(i, j) )
difz = 0.5 * ( z(i,1,j-1) - ht(i,j-1) &
+ z(i,1,j ) - ht(i,j ) )
IF( reg > 3.5 .AND. difz > 10.0 ) THEN
znt0 = 0.5 * ( znt(i, j-1) + znt(i, j) )
IF( znt0 <= 0.2) THEN
wndcor_v(i,j) = 1.0+0.320*znt0**0.2
ELSE
wndcor_v(i,j) = 1.169+0.315*znt0
ENDIF
wr14 = log(40.0/0.05)
wr1z = log(difz/0.05)
wr24 = log(40.0/1.0)
wr2z = log(difz/1.0)
wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
wndcor_v(i,j) = wndfac*wndcor_v(i,j)
ENDIF
ENDDO
ENDDO
IF ( wrf_dm_on_monitor()) THEN
IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
' D0',id,' sample wndcor_v values at i,j=', i0, j0, &
' wndcor_v=', wndcor_v(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
ENDIF
ENDIF
ELSE
wndcor_v(:,:) = 1.0
ENDIF
!
! Compute saturation mixing ratio so that nudging to a super-saturated state
! is not allowed.
!
IF( iqsat == 1 ) THEN
DO j=jts,jtf
DO k=kts,ktf
DO i=its,itf
es = SVP1*EXP(SVP2*(t3d(i,k,j)-SVPT0)/(t3d(i,k,j)-SVP3)) * 10.0 ! mb
qsat(i,k,j) = EP_2*es/(p3d(i,k,j)/100.0-es)
ENDDO
ENDDO
ENDDO
IF ( wrf_dm_on_monitor()) THEN
IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
DO k = kts, kte
WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
' D0',id,' sample moisture values (g/kg) at i,k,j=', i0, k, j0, &
' qv3d=', qv3d(i0,k,j0)*1000.0, ' qsat=', qsat(i0,k,j0)*1000.0
CALL wrf_message( TRIM(message) )
ENDDO
ENDIF
ENDIF
ENDIF
ENDIF
!
! Obs data density weighting.
!
IF( idd == 1 ) THEN
IF( rinblw < 0.001 ) THEN
IF ( wrf_dm_on_monitor()) THEN
WRITE(message,'(a)') 'Error in rinblw, please specify a reasonable value ***'
CALL wrf_message( TRIM(message) )
ENDIF
CALL wrf_error_fatal('In grid FDDA')
ENDIF
rindx = rinblw*1000.0/dx
m = -0.8*2.0/rindx
DO j=MAX(jts-2,jds),MIN(jte+2,jde-1)
DO i=MAX(its-2,ids),MIN(ite+2,ide-1)
IF( odis_ndg_old(i,j) < 0.5*rinblw ) THEN
blw_old(i,j) = 1.0
ELSE
x = min( odis_ndg_old(i,j)*1000./dx, rindx )
blw_old(i,j) = m * x + b
ENDIF
IF( odis_ndg_new(i,j) < 0.5*rinblw ) THEN
blw_new(i,j) = 1.0
ELSE
x = min( odis_ndg_new(i,j)*1000./dx, rindx )
blw_new(i,j) = m * x + b
ENDIF
ENDDO
ENDDO
! Smoother applies one point outside the tile, but one point in from boundaries
CALL smther(blw_old, its-2,ite+2, jts-2,jte+2, 1, &
MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
CALL smther(blw_new, its-2,ite+2, jts-2,jte+2, 1, &
MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
WHERE ( blw_old > 1.0)
blw_old = 1.0
END WHERE
WHERE ( blw_new > 1.0)
blw_new = 1.0
END WHERE
WHERE ( blw_old < 0.0)
blw_old = 0.0
END WHERE
WHERE ( blw_new < 0.0)
blw_new = 0.0
END WHERE
IF ( wrf_dm_on_monitor()) THEN
IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
WRITE(message,'(a,i1,a,2i4,4(a,f10.4))') &
' D0',id,' sample blw values at i,j=', i0, j0, &
' odis_ndg_old=', odis_ndg_old(i0,j0), ' km odis_ndg_new=', odis_ndg_new(i0,j0), &
' km blw_old=', blw_old(i0,j0), ' blw_new=', blw_new(i0,j0)
CALL wrf_message( TRIM(message) )
ENDIF
ENDIF
ENDIF
ENDIF
!
! tfac_sfc for surface analysis nudging
!
IF( xtime >= actual_end_fdda_min_sfc-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min_sfc &
.AND. dtramp_min > 0.0 .AND. if_ramping == 1 ) &
!BPR BEGIN
!The original method assumed that to be here in the code with dtramp_min>0
!meant that we were after the valid time of *_ndg_new and so used the
!values *_ndg_old and *_ndg_new to extrapolate a current analysis value.
!HOWEVER, in practice once we get to the valid time of *_ndg_new WRF will
!read in the next pair of *_ndg_old and *_ndg_new values, even if we are
!at the beginning of the rampdown period. Since the *_ndg_new values have
!a valid time after the beginning of the rampdown period we do not want
!to use them as we would be using analyses valid after the supposed end
!time of the analysis nudging.
!coef = (xtime-xtime_old_sfc+analysis_interval_sfc)/(analysis_interval_sfc*1.0)
!Use only the *_ndg_old values
coef = 0.0
!BPR END
! print*, 'coef =', xtime_old_sfc, xtime, xtime_new_sfc, coef
!
! Compute surface analysis nudging tendencies for u, v, t and q
!
! FASDAS
!
IF( fasdas == 1 ) THEN
!Edit TWG2015 Add new variable wpblfasdas to avoid altering global wpbl
!field when using the FASDAS option
wpblfasdas = 1.0
DO j=jts,jtf
DO k=kts,ktf
DO i=its,itf
if( k == 1 ) then
wpblfasdas(i, k, j) = 0.0
else
wpblfasdas(i, k, j) = 1.0
endif
ENDDO
ENDDO
ENDDO
ENDIF
!
! END FASDAS
!
DO j=jts,jtf
DO k=kts,ktf
DO i=itsu,itf
IF( idd == 1 ) THEN
blw = 0.5* (blw_old(i-1,j)+blw_old(i,j)) * ( 1.0 - coef ) &
+ 0.5* (blw_new(i-1,j)+blw_new(i,j)) * coef
ELSE
blw = 1.0
ENDIF
val_analysis = u10_ndg_old(i,j) *( 1.0 - coef ) + u10_ndg_new(i,j) * coef
val_analysis = val_analysis * wndcor_u(i,j)
RUNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,1)) * wzfac(k,1) * tfac_sfc * blw * &
( val_analysis - u3d(i,1,j) )
ENDDO
ENDDO
ENDDO
DO j=jtsv,jtf
DO k=kts,ktf
DO i=its,itf
IF( idd == 1 ) THEN
blw = 0.5* (blw_old(i,j-1)+blw_old(i,j)) * ( 1.0 - coef ) &
+ 0.5* (blw_new(i,j-1)+blw_new(i,j)) * coef
ELSE
blw = 1.0
ENDIF
val_analysis = v10_ndg_old(i,j) *( 1.0 - coef ) + v10_ndg_new(i,j) * coef
val_analysis = val_analysis * wndcor_v(i,j)
RVNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,2)) * wzfac(k,2) * tfac_sfc * blw * &
( val_analysis - v3d(i,1,j) )
ENDDO
ENDDO
ENDDO
DO j=jts,jtf
DO k=kts,ktf
DO i=its,itf
IF( idd == 1 ) THEN
blw = blw_old(i,j) * ( 1.0 - coef ) + blw_new(i,j) * coef
ELSE
blw = 1.0
ENDIF
!ckay2015 and TWG2015 add stability factor for stable regime, ensure
!consistency between Direct and Indirect nudging, and Convert potential
!temperature tendency to temperature tendency
!
! FASDAS
!
IF( fasdas == 1 ) THEN
if(regime(i,j).le.1.1) then
stabFac = 1.0/3.0
else
stabFac = 1.0
end if
val_analysis = th2_ndg_old(i,j) *( 1.0 - coef ) + th2_ndg_new(i,j) * coef
!BPR BEGIN
!RTHNDGDTEN(i,k,j) = stabFac*gt_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,3) * tfac * blw * &
RTHNDGDTEN(i,k,j) = stabFac*gt_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,3) * tfac_sfc * blw * &
!BPR END
( val_analysis - th3d(i,1,j))
DiffTN = val_analysis - th3d(i,1,j)
if(k.eq.1) then
exf = (1.0E05/p3d(i,k,j))**(287./1004.)
SDA_HFX(i,j) = RTHNDGDTEN(i,k,j)/exf !TWG 2015
else
RTHNDGDTEN(i,k,j) = 0.0
end if
val_analysis = q2_ndg_old(i,j) *( 1.0 - coef ) + q2_ndg_new(i,j) * coef
IF( iqsat == 1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
!BPR BEGIN
!RQVNDGDTEN(i,k,j) = stabFac*gq_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,4) * tfac * blw * &
RQVNDGDTEN(i,k,j) = stabFac*gq_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,4) * tfac_sfc * blw * &
!BPR END
( val_analysis - qv3d(i,k,j) )
if(k.eq.1) then
SDA_QFX(i,j) = RQVNDGDTEN(i,k,j)
else
RQVNDGDTEN(i,k,j) = 0.0
end if
ELSE
val_analysis = th2_ndg_old(i,j) *( 1.0 - coef ) + th2_ndg_new(i,j) * coef
!BPR BEGIN
!RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac * blw * &
RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac_sfc * blw * &
!BPR END
( val_analysis - th3d(i,1,j))
val_analysis = q2_ndg_old(i,j) *( 1.0 - coef ) + q2_ndg_new(i,j) * coef
IF( iqsat == 1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
!BPR BEGIN
!RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac * blw * &
RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac_sfc * blw * &
!BPR END
( val_analysis - qv3d(i,k,j) )
ENDIF
!
! END FASDAS
!
ENDDO
ENDDO
ENDDO
!BPR BEGIN
!Diagnostic print
IF ( wrf_dm_on_monitor()) THEN
IF( dbg_level .GE. 10 ) THEN
i0 = (ite-its)/2+its
j0 = (jte-jts)/2+jts
WRITE(message,'(a,i1,a,f7.2,a,i4,a,i4,a,f7.2,a,f4.2,a,f7.2,a,f4.2,a,f4.2 )') &
' D0',id,' At xtime=',xtime,' SFC FDDA sample TH2 analysis (i=',i0,', j=',j0,') = (th2_ndg_old=', &
th2_ndg_old(i0,j0),') * ',1.0-coef,' + (th2_ndg_new=',th2_ndg_new(i0,j0),') * ',coef, &
' where tfac_sfc=',tfac_sfc
CALL wrf_message( TRIM(message) )
WRITE(message,'(a,i1,a,f7.2,a,i4,a,f7.2,a,f7.2)') &
' D0',id,' xtime_old_sfc=',xtime_old_sfc,' analysis_interval_sfc=',analysis_interval_sfc, &
' actual_end_fdda_min_sfc=', actual_end_fdda_min_sfc,' dtramp_min=',dtramp_min
CALL wrf_message( TRIM(message) )
END IF
END IF
!BPR END
END SUBROUTINE sfddagd
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE fddagdinit(id,rundgdten,rvndgdten,rthndgdten,rqvndgdten,rmundgdten,& 1,26
SDA_HFX, SDA_QFX, QNORM, HFX_BOTH, QFX_BOTH, fasdas,& !fasdas
HFX_FDDA, & !fasdas
run_hours, &
if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
guv, gt, gq, if_ramping, dtramp_min, end_fdda_hour, &
grid_sfdda, guv_sfc, gt_sfc, gq_sfc, &
restart, allowed_to_read, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!
INTEGER , INTENT(IN) :: id
LOGICAL, INTENT(IN) :: restart, allowed_to_read
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
rundgdten, &
rvndgdten, &
rthndgdten, &
rqvndgdten
!
! FASDAS
!
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: &
SDA_HFX, &
SDA_QFX, &
QNORM, HFX_BOTH, QFX_BOTH
INTEGER, INTENT(IN ) :: fasdas
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
HFX_FDDA
!
! END FASDAS
!
INTEGER, INTENT(IN) :: run_hours
INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
if_no_pbl_nudging_q, end_fdda_hour
INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
INTEGER, INTENT(IN) :: if_ramping, grid_sfdda
REAL, INTENT(IN) :: dtramp_min
REAL, INTENT(IN) :: guv, gt, gq
REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc
REAL :: actual_end_fdda_min
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: rmundgdten
INTEGER :: i, j, k
LOGICAL , EXTERNAL :: wrf_dm_on_monitor
CHARACTER (LEN=256) :: message
IF ( wrf_dm_on_monitor() ) THEN
IF( guv > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for wind is applied and Guv= ', guv
CALL wrf_message(TRIM(message))
ELSE IF( guv < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for wind is not applied and Guv= ', guv
CALL wrf_message(TRIM(message))
ENDIF
IF( gt > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for temperature is applied and Gt= ', gt
CALL wrf_message(TRIM(message))
ELSE IF( gt < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for temperature is not applied and Gt= ', gt
CALL wrf_message(TRIM(message))
ENDIF
IF( gq > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for water vapor mixing ratio is applied and Gq= ', gq
CALL wrf_message(TRIM(message))
ELSE IF( gq < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Gq must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' 3-D analysis nudging for water vapor mixing ratio is not applied and Gq= ', gq
CALL wrf_message(TRIM(message))
ENDIF
IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
WRITE(message,'(a,i1,a)') &
'D0',id,' 3-D analysis nudging for wind is turned off within the PBL.'
CALL wrf_message(TRIM(message))
ENDIF
IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
WRITE(message,'(a,i1,a)') &
'D0',id,' 3-D analysis nudging for temperature is turned off within the PBL.'
CALL wrf_message(TRIM(message))
ENDIF
IF( gq > 0.0 .AND. if_no_pbl_nudging_q == 1 ) THEN
WRITE(message,'(a,i1,a)') &
'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off within the PBL.'
CALL wrf_message(TRIM(message))
ENDIF
IF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
WRITE(message,'(a,i1,a,i3)') &
'D0',id,' 3-D analysis nudging for wind is turned off below layer', k_zfac_uv
CALL wrf_message(TRIM(message))
ENDIF
IF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
WRITE(message,'(a,i1,a,i3)') &
'D0',id,' 3-D analysis nudging for temperature is turned off below layer', k_zfac_t
CALL wrf_message(TRIM(message))
ENDIF
IF( gq > 0.0 .AND. if_zfac_q == 1 ) THEN
WRITE(message,'(a,i1,a,i3)') &
'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off below layer', &
k_zfac_q
CALL wrf_message(TRIM(message))
ENDIF
IF( grid_sfdda >=1 ) THEN
IF( guv_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surface analysis nudging for wind is applied and Guv_sfc= ', guv_sfc
CALL wrf_message(TRIM(message))
ELSE IF( guv_sfc < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Guv_sfc must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surface analysis nudging for wind is not applied and Guv_sfc= ', guv_sfc
CALL wrf_message(TRIM(message))
ENDIF
IF( gt_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surface analysis nudging for temperature is applied and Gt_sfc= ', gt_sfc
CALL wrf_message(TRIM(message))
ELSE IF( gt_sfc < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Gt_sfc must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surafce analysis nudging for temperature is not applied and Gt_sfc= ', gt_sfc
CALL wrf_message(TRIM(message))
ENDIF
IF( gq_sfc > 0.0 ) THEN
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surface analysis nudging for water vapor mixing ratio is applied and Gq_sfc= ', gq_sfc
CALL wrf_message(TRIM(message))
ELSE IF( gq_sfc < 0.0 ) THEN
CALL wrf_error_fatal('In grid FDDA, Gq_sfc must be positive.')
ELSE
WRITE(message,'(a,i1,a,e12.4)') &
'D0',id,' surface analysis nudging for water vapor mixing ratio is not applied and Gq_sfc= ', gq_sfc
CALL wrf_message(TRIM(message))
ENDIF
ENDIF
IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
IF( dtramp_min <= 0.0 ) THEN
actual_end_fdda_min = end_fdda_hour*60.0
ELSE
actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
ENDIF
IF( actual_end_fdda_min <= run_hours*60. ) THEN
WRITE(message,'(a,i1,a)') &
'D0',id,' analysis nudging is ramped down near the end of the nudging period,'
CALL wrf_message(TRIM(message))
WRITE(message,'(a,f6.2,a,f6.2,a)') &
' starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0, &
'h, ending at ', actual_end_fdda_min/60.0,'h.'
CALL wrf_message(TRIM(message))
ENDIF
ENDIF
ENDIF
IF(.not.restart) THEN
DO j = jts,jte
DO k = kts,kte
DO i = its,ite
rundgdten(i,k,j) = 0.
rvndgdten(i,k,j) = 0.
rthndgdten(i,k,j) = 0.
rqvndgdten(i,k,j) = 0.
! TWG 2015 Psuedo radiative flux
HFX_FDDA(i,k,j) = 0.
! TWG END
if(k.eq.kts) rmundgdten(i,j) = 0.
ENDDO
ENDDO
ENDDO
!
! FASDAS
!
IF( fasdas == 1 ) THEN
DO j = jts,jte
DO i = its,ite
SDA_HFX(I,J) = 0.0
SDA_QFX(I,J) = 0.0
QNORM(I,J) = 0.0
HFX_BOTH(I,J) = 0.0
QFX_BOTH(I,J) = 0.0
ENDDO
ENDDO
ENDIF
!
! END FASDAS
!
ENDIF
END SUBROUTINE fddagdinit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE smther(slab, idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd) 2
!
! PURPOSE: SPATIALLY SMOOTH DATA IN SLAB TO DAMPEN SHORT
! WAVELENGTH COMPONENTS.
!
! Implemented based on the same smoothing subroutine used in MM5, with modifications to
! remove the extra loop that causes unneccessary desmoothing. Aijun Deng (Penn State),
! December 2008
!
IMPLICIT NONE
INTEGER :: idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd
INTEGER :: i, j, k, kk
REAL :: avg
REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB
REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB_ORIG
REAL, DIMENSION(2) :: XNU
IF(NPASS.EQ.0)RETURN
XNU(1)=0.50
XNU(2)=-0.52
DO K=1,NPASS
KK = 2 - MOD(K,2)
DO J=JDIMST,JDIMND
DO I=IDIMST,IDIMND
SLAB_ORIG(I,J) = SLAB(I,J)
END DO
END DO
DO J=JST,JND
DO I=IST,IND
AVG = ( SLAB_ORIG(I+1,J ) + &
SLAB_ORIG(I-1,J ) + &
SLAB_ORIG(I ,J+1) + &
SLAB_ORIG(I ,J-1) ) * 0.25
SLAB(I,J)=SLAB_ORIG(I,J)+XNU(KK)*(AVG - SLAB_ORIG(I,J))
ENDDO
ENDDO
ENDDO
END SUBROUTINE smther
!-------------------------------------------------------------------
END MODULE module_fdda_psufddagd