diff --git a/Makefile b/Makefile index 98a7d5ea75..8a82dd20dc 100644 --- a/Makefile +++ b/Makefile @@ -300,11 +300,11 @@ intel-mpi: # BUILDTARGET Intel compiler suite with Intel MPI library "CC_SERIAL = icc" \ "CXX_SERIAL = icpc" \ "FFLAGS_PROMOTION = -real-size 64" \ - "FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \ - "CFLAGS_OPT = -O3" \ - "CXXFLAGS_OPT = -O3" \ + "FFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback -no-wrap-margin -convert big_endian -free -align array64byte" \ + "CFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ + "CXXFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback" \ + "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback -no-wrap-margin" \ "CFLAGS_DEBUG = -g -traceback" \ "CXXFLAGS_DEBUG = -g -traceback" \ "LDFLAGS_DEBUG = -g -fpe0 -traceback" \ @@ -660,7 +660,7 @@ ifneq "$(LAPACK)" "" endif RM = rm -f -CPP = cpp -P -traditional +CPP = ${CXX_SERIAL} -E # Modified for use with the Intel C++ compiler RANLIB = ranlib ifdef CORE @@ -875,7 +875,7 @@ ifeq "$(findstring clean, $(MAKECMDGOALS))" "clean" # CHECK FOR CLEAN TARGET override AUTOCLEAN=false endif # END OF CLEAN TARGET CHECK -VER=$(shell git describe --dirty 2> /dev/null) +VER="v7.3.develop.EPA_FDDA" # Hard-coded specific version identifier #override CPPFLAGS += -DMPAS_GIT_VERSION=$(VER) ifeq "$(findstring v, $(VER))" "v" diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index b3127cce3b..96323553ee 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -801,6 +801,23 @@ +<<<<<<< atmosphere/develop-EPA_FDDA + + + + + + + + + + + + + + +======= +>>>>>>> develop @@ -1097,6 +1114,18 @@ #endif + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2066,6 +2123,86 @@ description="threshold above which the accumulated grid-scale precipitation is reset" possible_values="Positive real values"/> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + #ifdef DO_PHYSICS + + + + + + + + + + + + + + + diff --git a/src/core_atmosphere/physics/Makefile b/src/core_atmosphere/physics/Makefile index faf6b98d40..3fb9ce3436 100644 --- a/src/core_atmosphere/physics/Makefile +++ b/src/core_atmosphere/physics/Makefile @@ -21,6 +21,7 @@ OBJS = \ mpas_atmphys_driver.o \ mpas_atmphys_driver_cloudiness.o \ mpas_atmphys_driver_convection.o \ + mpas_atmphys_fdda.o \ mpas_atmphys_driver_gwdo.o \ mpas_atmphys_driver_lsm.o \ mpas_atmphys_driver_microphysics.o \ @@ -72,6 +73,7 @@ mpas_atmphys_control.o: \ mpas_atmphys_driver.o: \ mpas_atmphys_driver_cloudiness.o \ mpas_atmphys_driver_convection.o \ + mpas_atmphys_fdda.o \ mpas_atmphys_driver_gwdo.o \ mpas_atmphys_driver_lsm.o \ mpas_atmphys_driver_pbl.o \ @@ -93,6 +95,11 @@ mpas_atmphys_driver_convection.o: \ mpas_atmphys_utilities.o \ mpas_atmphys_vars.o +mpas_atmphys_fdda.o: \ + mpas_atmphys_constants.o \ + mpas_atmphys_utilities.o \ + mpas_atmphys_vars.o + mpas_atmphys_driver_gwdo.o: \ mpas_atmphys_vars.o diff --git a/src/core_atmosphere/physics/mpas_atmphys_control.F b/src/core_atmosphere/physics/mpas_atmphys_control.F index 9b7a08c5e0..7704756b3e 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_control.F +++ b/src/core_atmosphere/physics/mpas_atmphys_control.F @@ -49,6 +49,8 @@ module mpas_atmphys_control ! Laura D. Fowler (laura@ucar.edu) / 2014-09-18. ! * renamed "wsm6" to "mp_wsm6" and "kessler" to "mp_kessler". ! Laura D. Fowler (laura@ucar.edu) / 2016-03-09. +! * added four-dimensional data assimilation scheme (based on WRF FDDA scheme). +! O. Russell Bullock Jr. (bullock.russell@epa.gov) / 2016-03-21. ! * renamed "kain_fritsch" to "cu_kain_fritsch" and "tiedtke" to "cu_tiedtke". ! Laura D. Fowler (laura@ucar.edu) / 2016-03-22. ! * renamed "ysu" to "bl_ysu", "ysu_gwdo" to "bl_gwdo_ysu", and "monin_obukhov" to "sf_monin_obukhov". @@ -94,7 +96,8 @@ subroutine physics_namelist_check(configs) config_radt_cld_scheme, & config_radt_lw_scheme, & config_radt_sw_scheme, & - config_sfclayer_scheme + config_sfclayer_scheme, & + config_fdda_scheme !----------------------------------------------------------------------------------------------------------------- !call mpas_log_write('') @@ -110,6 +113,7 @@ subroutine physics_namelist_check(configs) call mpas_pool_get_config(configs,'config_radt_lw_scheme' ,config_radt_lw_scheme ) call mpas_pool_get_config(configs,'config_radt_sw_scheme' ,config_radt_sw_scheme ) call mpas_pool_get_config(configs,'config_sfclayer_scheme' ,config_sfclayer_scheme ) + call mpas_pool_get_config(configs,'config_fdda_scheme' ,config_fdda_scheme ) call mpas_log_write('') call mpas_log_write('----- Setting up physics suite '''//trim(config_physics_suite)//''' -----') @@ -128,6 +132,7 @@ subroutine physics_namelist_check(configs) if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction' if (trim(config_sfclayer_scheme) == 'suite') config_sfclayer_scheme = 'sf_monin_obukhov' if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'noah' + if (trim(config_fdda_scheme) == 'suite') config_fdda_scheme = 'off' else if (trim(config_physics_suite) == 'convection_permitting') then @@ -140,6 +145,7 @@ subroutine physics_namelist_check(configs) if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction' if (trim(config_sfclayer_scheme) == 'suite') config_sfclayer_scheme = 'sf_mynn' if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'noah' + if (trim(config_fdda_scheme) == 'suite') config_fdda_scheme = 'off' else if (trim(config_physics_suite) == 'none') then @@ -152,6 +158,7 @@ subroutine physics_namelist_check(configs) if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'off' if (trim(config_sfclayer_scheme) == 'suite') config_sfclayer_scheme = 'off' if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'off' + if (trim(config_fdda_scheme) == 'suite') config_fdda_scheme = 'off' else @@ -161,6 +168,17 @@ subroutine physics_namelist_check(configs) end if +!fdda scheme: + if(.not. (config_fdda_scheme .eq. 'off' .or. & + config_fdda_scheme .eq. 'analysis' .or. & + config_fdda_scheme .eq. 'scaled' )) then + + write(mpas_err_message,'(A,A10)') 'illegal value for config_fdda_scheme:', & + trim(config_fdda_scheme) + call physics_error_fatal(mpas_err_message) + + endif + !cloud microphysics scheme: if(.not. (config_microp_scheme .eq. 'off' .or. & config_microp_scheme .eq. 'mp_kessler' .or. & @@ -307,6 +325,7 @@ subroutine physics_namelist_check(configs) call mpas_log_write(' config_radt_sw_scheme = '//trim(config_radt_sw_scheme)) call mpas_log_write(' config_sfclayer_scheme = '//trim(config_sfclayer_scheme)) call mpas_log_write(' config_lsm_scheme = '//trim(config_lsm_scheme)) + call mpas_log_write(' config_fdda_scheme = '//trim(config_fdda_scheme)) call mpas_log_write('') end subroutine physics_namelist_check diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver.F b/src/core_atmosphere/physics/mpas_atmphys_driver.F index b120d0ccc4..75f7946e8e 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver.F @@ -19,6 +19,7 @@ module mpas_atmphys_driver use mpas_atmphys_driver_radiation_sw use mpas_atmphys_driver_sfclayer use mpas_atmphys_driver_oml + use mpas_atmphys_fdda use mpas_atmphys_constants use mpas_atmphys_interface use mpas_atmphys_update @@ -68,6 +69,7 @@ module mpas_atmphys_driver ! driver_radiation_sw : driver for short wave radiation schemes. ! driver_radiation_lw : driver for long wave radiation schemes. ! driver_sfclayer : driver for surface layer scheme. +! apply_fdda : applies four-dimensional data assimilation. ! update_convection_step1 : updates lifetime of deep convective clouds in Kain-Fritsch scheme. ! update_convection_step2 : updates accumulated precipitation output from convection schemes. ! update_radiation_diagnostics: updates accumualted radiation diagnostics from radiation schemes. @@ -85,13 +87,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 @@ -117,6 +119,7 @@ subroutine physics_driver(domain,itimestep,xtime_s) type(mpas_pool_type),pointer:: configs, & mesh, & state, & + fdda, & diag, & diag_physics, & tend_physics, & @@ -127,6 +130,7 @@ subroutine physics_driver(domain,itimestep,xtime_s) character(len=StrKIND),pointer:: config_bucket_update, & config_convection_scheme, & + config_fdda_scheme, & config_gwdo_scheme, & config_lsm_scheme, & config_pbl_scheme, & @@ -153,6 +157,7 @@ subroutine physics_driver(domain,itimestep,xtime_s) call mpas_timer_start('physics driver') call mpas_pool_get_config(domain%configs,'config_convection_scheme',config_convection_scheme) + call mpas_pool_get_config(domain%configs,'config_fdda_scheme' ,config_fdda_scheme ) call mpas_pool_get_config(domain%configs,'config_gwdo_scheme' ,config_gwdo_scheme ) call mpas_pool_get_config(domain%configs,'config_lsm_scheme' ,config_lsm_scheme ) call mpas_pool_get_config(domain%configs,'config_pbl_scheme' ,config_pbl_scheme ) @@ -169,13 +174,15 @@ 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' .or. & + config_fdda_scheme .ne. 'off') then block => domain % blocklist do while(associated(block)) 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 ) call mpas_pool_get_subpool(block%structs,'diag' ,diag ) call mpas_pool_get_subpool(block%structs,'diag_physics',diag_physics) call mpas_pool_get_subpool(block%structs,'atm_input' ,atm_input ) @@ -332,6 +339,17 @@ subroutine physics_driver(domain,itimestep,xtime_s) !$OMP END PARALLEL DO end if + !call to FDDA scheme within MPAS physics (not WRF physics) + if(config_fdda_scheme .ne. 'off') then + time_lev = 1 +!$OMP PARALLEL DO + do thread=1,nThreads + call apply_fdda(xtime_s,block%configs,mesh,state,time_lev,fdda,diag,diag_physics,tend_physics, & + cellSolveThreadStart(thread), cellSolveThreadEnd(thread)) + end do +!$OMP END PARALLEL DO + endif + !deallocate arrays shared by all physics parameterizations: call deallocate_forall_physics(block%configs) diff --git a/src/core_atmosphere/physics/mpas_atmphys_fdda.F b/src/core_atmosphere/physics/mpas_atmphys_fdda.F new file mode 100644 index 0000000000..ee7b966c9f --- /dev/null +++ b/src/core_atmosphere/physics/mpas_atmphys_fdda.F @@ -0,0 +1,254 @@ +!================================================================================================== + module mpas_atmphys_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 + + implicit none + private + public:: apply_fdda + + integer, private:: i,k + + +!MPAS grid nudging. +!O. Russell Bullock Jr. (send comments to bullock.russell@epa.gov). +!2016-03-30. +! +! subroutines in mpas_atmphys_fdda: +! ---------------------------------------------- +! apply_fdda : main fdda subroutine (called from subroutine physics_driver). +! + +! 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. + + contains + +!================================================================================================== + subroutine apply_fdda(xtime_s,configs,mesh,state,time_lev,fdda,diag,diag_physics,tend_physics, & + its,ite) +!================================================================================================== + +!input variables: + 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):: its,ite + integer,intent(in):: time_lev + real(kind=RKIND),intent(in):: xtime_s + +!inout variables: + type(mpas_pool_type),intent(inout):: tend_physics + +!local pointers: + integer,pointer:: index_qv,fdda_t_min_layer,fdda_q_min_layer,fdda_uv_min_layer + integer,dimension(:),pointer :: kpbl + 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 + real(kind=RKIND),dimension(:,:),pointer :: qv_fdda_old,th_fdda_old,u_fdda_old,v_fdda_old + real(kind=RKIND),dimension(:,:),pointer :: qv_fdda_new,th_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 + 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: + real(kind=RKIND),parameter:: rvrd = R_v/R_d + +!================================================================================================== + + 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_array(fdda,'qv_fdda_old',qv_fdda_old) + call mpas_pool_get_array(fdda,'th_fdda_old',th_fdda_old) + call mpas_pool_get_array(fdda,'u_fdda_old', u_fdda_old ) + call mpas_pool_get_array(fdda,'v_fdda_old', v_fdda_old ) + call mpas_pool_get_array(fdda,'qv_fdda_new',qv_fdda_new) + call mpas_pool_get_array(fdda,'th_fdda_new',th_fdda_new) + call mpas_pool_get_array(fdda,'u_fdda_new', u_fdda_new ) + call mpas_pool_get_array(fdda,'v_fdda_new', v_fdda_new ) + call mpas_pool_get_array(tend_physics,'rthfddaten',rthfddaten) + call mpas_pool_get_array(tend_physics,'rqvfddaten',rqvfddaten) + call mpas_pool_get_array(tend_physics,'rufddaten' ,rufddaten ) + call mpas_pool_get_array(tend_physics,'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_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) + +!check to see if 'old' FDDA target values need to be defined from initial conditions + if(first_call .and. .not. config_do_restart) then + qv_fdda_old = qv + th_fdda_old = theta_m / (1._RKIND + rvrd * qv) + u_fdda_old = u + v_fdda_old = v + first_call = .false. + endif + +!calculate time fraction within FDDA target data interval + tfrac = mod(xtime_s,fdda_int)/fdda_int + + fdda_select: select case(fdda_scheme) + + case ("scaled") + + srange = fdda_scale_max-fdda_scale_min + + DO i=its,ite + + dcEdge_m = len_disp / meshDensity(i)**0.25 + + if (dcEdge_m .lt. fdda_scale_max) then + if (dcEdge_m .lt. fdda_scale_min) then + sfac = 0.0 + else + sfac = (dcEdge_m-fdda_scale_min)/srange + endif + else + sfac = 1.0 + endif + + DO k=kts,kte + + theta_p = theta_m(k,i) / (1._RKIND + rvrd * 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-tfrac)*th_fdda_old(k,i)+tfrac*th_fdda_new(k,i) + rthfddaten(k,i) = fdda_t_coef * ( t_target - theta_p ) * sfac + else + rthfddaten(k,i) = 0. + 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-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. + 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-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-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. + rvfddaten(k,i) = 0. + endif + + ENDDO + ENDDO + + case("analysis") + + DO i=its,ite + DO k=kts,kte + + theta_p = theta_m(k,i) / (1._RKIND + rvrd * 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-tfrac)*th_fdda_old(k,i)+tfrac*th_fdda_new(k,i) + rthfddaten(k,i) = fdda_t_coef * ( t_target - theta_p ) + else + rthfddaten(k,i) = 0. + 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-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. + 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-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-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. + rvfddaten(k,i) = 0. + endif + + ENDDO + ENDDO + + case default + + end select fdda_select + +!check to see if this is the last time step within the current FDDA target data interval and +!if so, save the current "new" targets as the "old" targets for the next time step + if(mod(xtime_s+config_dt,fdda_int)/fdda_int .lt. tfrac) then + th_fdda_old = th_fdda_new + qv_fdda_old = qv_fdda_new + u_fdda_old = u_fdda_new + v_fdda_old = v_fdda_new + endif + + end subroutine apply_fdda + +!================================================================================================== + end module mpas_atmphys_fdda +!================================================================================================== diff --git a/src/core_atmosphere/physics/mpas_atmphys_manager.F b/src/core_atmosphere/physics/mpas_atmphys_manager.F index fe8ee5c27c..dc373e0d38 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_manager.F +++ b/src/core_atmosphere/physics/mpas_atmphys_manager.F @@ -386,7 +386,8 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) config_lsm_scheme, & config_microp_scheme, & config_radt_lw_scheme, & - config_radt_sw_scheme + config_radt_sw_scheme, & + config_fdda_scheme character(len=StrKIND),pointer:: config_conv_interval, & config_pbl_interval, & @@ -423,6 +424,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call mpas_pool_get_config(configs,'config_microp_scheme' ,config_microp_scheme ) call mpas_pool_get_config(configs,'config_radt_lw_scheme' ,config_radt_lw_scheme ) call mpas_pool_get_config(configs,'config_radt_sw_scheme' ,config_radt_sw_scheme ) + call mpas_pool_get_config(configs,'config_fdda_scheme' ,config_fdda_scheme ) call mpas_pool_get_config(configs,'config_conv_interval' ,config_conv_interval ) call mpas_pool_get_config(configs,'config_pbl_interval' ,config_pbl_interval ) diff --git a/src/core_atmosphere/physics/mpas_atmphys_todynamics.F b/src/core_atmosphere/physics/mpas_atmphys_todynamics.F index c7baeb213a..4474910a52 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_todynamics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_todynamics.F @@ -78,8 +78,9 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi real(kind=RKIND),dimension(:,:) :: tend_ru_physics, tend_rtheta_physics, tend_rho_physics !local variables: - character(len=StrKIND), pointer :: config_pbl_scheme, config_convection_scheme, & - config_radt_lw_scheme, config_radt_sw_scheme + character(len=StrKIND), pointer :: config_pbl_scheme, config_convection_scheme, & + config_radt_lw_scheme, config_radt_sw_scheme, & + config_fdda_scheme integer:: i,iCell,k,n integer,pointer:: index_qv, index_qc, index_qr, index_qi, index_qs, index_qg @@ -96,6 +97,7 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi real(kind=RKIND),dimension(:,:),pointer:: rthcuten,rqvcuten,rqccuten, & rqrcuten,rqicuten,rqscuten, & rucuten,rvcuten + real(kind=RKIND),dimension(:,:),pointer:: rthfddaten,rqvfddaten,rufddaten,rvfddaten real(kind=RKIND),dimension(:,:),pointer:: rthratenlw,rthratensw real(kind=RKIND),dimension(:,:),pointer:: rthdynten @@ -106,7 +108,7 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi real(kind=RKIND):: coeff real(kind=RKIND):: tem - real(kind=RKIND),dimension(:,:),pointer:: rublten_Edge,rucuten_Edge + real(kind=RKIND),dimension(:,:),pointer:: rublten_Edge,rucuten_Edge,rufddaten_Edge real(kind=RKIND),dimension(:,:),allocatable:: theta,tend_th @@ -121,6 +123,7 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi call mpas_pool_get_config(configs, 'config_convection_scheme', config_convection_scheme) call mpas_pool_get_config(configs, 'config_radt_lw_scheme', config_radt_lw_scheme) call mpas_pool_get_config(configs, 'config_radt_sw_scheme', config_radt_sw_scheme) + call mpas_pool_get_config(configs, 'config_fdda_scheme', config_fdda_scheme) call mpas_pool_get_array(state, 'theta_m', theta_m, 1) call mpas_pool_get_array(state, 'scalars', scalars, 1) @@ -159,6 +162,12 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi call mpas_pool_get_array(tend_physics, 'rqscuten', rqscuten) call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten) + call mpas_pool_get_array(tend_physics, 'rufddaten', rufddaten) + call mpas_pool_get_array(tend_physics, 'rvfddaten', rvfddaten) + call mpas_pool_get_array(tend_physics, 'rthfddaten', rthfddaten) + call mpas_pool_get_array(tend_physics, 'rqvfddaten', rqvfddaten) + call mpas_pool_get_array(tend_physics, 'rufddaten_Edge', rufddaten_Edge) + call mpas_pool_get_array(tend,'rt_diabatic_tend',tend_diabatic) call mpas_pool_get_array(tend_physics, 'rthratenlw', rthratenlw) @@ -203,16 +212,22 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi if (.not. associated(rqvcuten)) allocate(rqvcuten(0,0)) if (.not. associated(rqccuten)) allocate(rqccuten(0,0)) if (.not. associated(rqicuten)) allocate(rqicuten(0,0)) + if (.not. associated(rufddaten)) allocate(rufddaten(0,0)) + if (.not. associated(rvfddaten)) allocate(rvfddaten(0,0)) + if (.not. associated(rthfddaten)) allocate(rthfddaten(0,0)) + if (.not. associated(rqvfddaten)) allocate(rqvfddaten(0,0)) call physics_get_tend_work(block, mesh, nCells, nEdges, nCellsSolve, nEdgesSolve, & rk_step, dynamics_substep, & config_pbl_scheme, config_convection_scheme, config_radt_lw_scheme, config_radt_sw_scheme, & + config_fdda_scheme, & index_qv, index_qc, index_qr, index_qi, index_qs, index_ni, & rublten, rvblten, mass_edge, rublten_Edge, & tend_ru_physics, & rucuten, rvcuten, rucuten_Edge, & tend_th, tend_scalars, mass, rthblten, rqvblten, rqcblten, rqiblten, rniblten, & rthcuten, rqvcuten, rqccuten, rqrcuten, rqicuten, rqscuten, & + rufddaten,rvfddaten,rthfddaten,rqvfddaten,rufddaten_Edge, & rthratenlw, rthratensw, rthdynten, & tend_u_phys, tend_rtheta_adv, tend_diabatic, & theta_m, scalars, & @@ -240,6 +255,10 @@ subroutine physics_get_tend( block, mesh, state, diag, tend, tend_physics, confi if (size(rqvcuten) == 0) deallocate(rqvcuten) if (size(rqccuten) == 0) deallocate(rqccuten) if (size(rqicuten) == 0) deallocate(rqicuten) + if (size(rufddaten) == 0) deallocate(rufddaten) + if (size(rvfddaten) == 0) deallocate(rvfddaten) + if (size(rthfddaten) == 0) deallocate(rthfddaten) + if (size(rqvfddaten) == 0) deallocate(rqvfddaten) ! deallocate(theta) deallocate(tend_th) @@ -265,11 +284,13 @@ end subroutine physics_get_tend subroutine physics_get_tend_work(block, mesh, nCells, nEdges, nCellsSolve, nEdgesSolve, & rk_step, dynamics_substep, & config_pbl_scheme, config_convection_scheme, config_radt_lw_scheme, config_radt_sw_scheme, & + config_fdda_scheme, & index_qv, index_qc, index_qr, index_qi, index_qs, index_ni, & rublten, rvblten, mass_edge, rublten_Edge, tend_u, & rucuten, rvcuten, rucuten_Edge, & tend_th, tend_scalars, mass, rthblten, rqvblten, rqcblten, rqiblten, rniblten, & rthcuten, rqvcuten, rqccuten, rqrcuten, rqicuten, rqscuten, & + rufddaten,rvfddaten,rthfddaten,rqvfddaten,rufddaten_Edge, & rthratenlw, rthratensw, rthdynten, & tend_u_phys, tend_rtheta_adv, tend_diabatic, & theta_m, scalars, tend_theta, tend_theta_euler & @@ -288,6 +309,7 @@ subroutine physics_get_tend_work(block, mesh, nCells, nEdges, nCellsSolve, nEdge character(len=StrKIND), intent(in) :: config_convection_scheme character(len=StrKIND), intent(in) :: config_radt_lw_scheme character(len=StrKIND), intent(in) :: config_radt_sw_scheme + character(len=StrKIND), intent(in) :: config_fdda_scheme integer, intent(in) :: index_qv, index_qc, index_qr, index_qi, index_qs, index_ni real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rublten real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rvblten @@ -311,6 +333,11 @@ subroutine physics_get_tend_work(block, mesh, nCells, nEdges, nCellsSolve, nEdge real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rqrcuten real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rqicuten real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rqscuten + real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rufddaten + real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rvfddaten + real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rthfddaten + real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rqvfddaten + real (kind=RKIND), dimension(nVertLevels,nEdges+1), intent(inout) :: rufddaten_Edge real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rthratenlw real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(in) :: rthratensw real (kind=RKIND), dimension(nVertLevels,nCells+1), intent(inout) :: rthdynten @@ -403,6 +430,25 @@ subroutine physics_get_tend_work(block, mesh, nCells, nEdges, nCellsSolve, nEdge end select convection_select endif + !add coupled tendencies due to four-dimensional data assimilation (fdda): + if(config_fdda_scheme .ne. 'off') then + + rufddaten_Edge(:,:) = 0. + call tend_toEdges(block,mesh,rufddaten,rvfddaten,rufddaten_Edge) + do i = 1, nEdgesSolve + do k = 1, nVertLevels + tend_u(k,i)=tend_u(k,i)+rufddaten_Edge(k,i)*mass_edge(k,i) + enddo + enddo + + do i = 1, nCellsSolve + do k = 1, nVertLevels + tend_th(k,i) = tend_th(k,i) + rthfddaten(k,i)*mass(k,i) + tend_scalars(index_qv,k,i) = tend_scalars(index_qv,k,i) + rqvfddaten(k,i)*mass(k,i) + enddo + enddo + endif + !add coupled tendencies due to longwave radiation: if (config_radt_lw_scheme .ne. 'off') then do i = 1, nCellsSolve diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 7bc8d53146..4afc0c53af 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -551,6 +551,9 @@ + + +