diff --git a/CMakeLists.txt b/CMakeLists.txt index df1fc2344..4b0947f55 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -165,6 +165,7 @@ add_dependencies(fv3atm ccppdriver ccppdata ccppphys ccpp) target_link_libraries(fv3atm PUBLIC ccppdriver ccppdata ccppphys ccpp) target_include_directories(fv3atm PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/stochastic_physics) +target_include_directories(fv3atm PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/../stochastic_physics) target_compile_definitions(fv3atm PRIVATE "${_fv3atm_defs_private}") target_link_libraries(fv3atm PUBLIC fv3dycore @@ -192,7 +193,7 @@ endif() ### Install ############################################################################### install( - TARGETS fv3atm fv3dycore io ${CCPP_LIBRARIES} cpl stochastic_physics stochastic_physics_wrapper + TARGETS fv3atm fv3dycore io ${CCPP_LIBRARIES} cpl stochastic_physics_wrapper EXPORT fv3atm-config LIBRARY DESTINATION lib ARCHIVE DESTINATION lib) diff --git a/atmos_model.F90 b/atmos_model.F90 index 3730da692..04ecae56b 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -869,6 +869,7 @@ end subroutine update_atmos_model_state ! subroutine atmos_model_end (Atmos) + use get_stochy_pattern_mod, only: write_stoch_restart_atm type (atmos_data_type), intent(inout) :: Atmos !---local variables integer :: idx, seconds, ierr @@ -878,12 +879,12 @@ subroutine atmos_model_end (Atmos) call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst) - call stochastic_physics_wrapper_end(GFS_control) - if(restart_endfcst) then call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, & GFS_control, Atmos%domain) + call write_stoch_restart_atm('RESTART/atm_stoch.res.nc') endif + call stochastic_physics_wrapper_end(GFS_control) ! Fast physics (from dynamics) are finalized in atmosphere_end above; ! standard/slow physics (from IPD) are finalized in CCPP_step 'finalize'. diff --git a/ccpp/CMakeLists.txt b/ccpp/CMakeLists.txt index 760e09a8e..7f7575cd6 100644 --- a/ccpp/CMakeLists.txt +++ b/ccpp/CMakeLists.txt @@ -220,18 +220,6 @@ else (SIONLIB) message (STATUS "Disable SIONlib support") endif (SIONLIB) -#------------------------------------------------------------------------------ -# Set Intel MKL flags for preprocessor, compiler and linker (if defined) -if(MKL_DIR) - set (MKL_INC "-m64 -I${MKL_DIR}/include") - set (MKL_LIB "-L${MKL_DIR}/lib -Wl,-rpath,${MKL_DIR}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl") - set (CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${MKL_INC} ${MKL_LIB}") - ADD_DEFINITIONS(-DMKL) - message (STATUS "Enable Intel MKL support") -else(MKL_DIR) - message (STATUS "Disable Intel MKL support") -endif(MKL_DIR) - #------------------------------------------------------------------------------ # Set netCDF flags for preprocessor, compiler and linker (if defined) # Legacy settings for old make build diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90 index 913432c0e..44d7fe393 100644 --- a/ccpp/data/GFS_typedefs.F90 +++ b/ccpp/data/GFS_typedefs.F90 @@ -1065,11 +1065,15 @@ module GFS_typedefs !--- stochastic physics control parameters logical :: do_sppt + logical :: pert_clds + logical :: pert_radtend + logical :: pert_mp logical :: use_zmtnblck logical :: do_shum logical :: do_skeb integer :: skeb_npass integer :: lndp_type + real(kind=kind_phys) :: sppt_amp ! pjp cloud perturbations integer :: n_var_lndp logical :: lndp_each_step ! flag to indicate that land perturbations are applied at every time step, ! otherwise they are applied only after gcycle is run @@ -1282,7 +1286,7 @@ module GFS_typedefs real (kind=kind_phys), pointer :: acvt (:) => null() !< arrays used by cnvc90 top (cnvc90.f) !--- Stochastic physics properties calculated in physics_driver - real (kind=kind_phys), pointer :: dtdtr (:,:) => null() !< temperature change due to radiative heating per time step (K) + real (kind=kind_phys), pointer :: dtdtnp (:,:) => null() !< temperature change from physics that should not be perturbed with SPPT (k) real (kind=kind_phys), pointer :: dtotprcp (:) => null() !< change in totprcp (diag_type) real (kind=kind_phys), pointer :: dcnvprcp (:) => null() !< change in cnvprcp (diag_type) real (kind=kind_phys), pointer :: drain_cpl (:) => null() !< change in rain_cpl (coupling_type) @@ -1531,6 +1535,7 @@ module GFS_typedefs real (kind=kind_phys), pointer :: skebv_wts(:,:) => null() !< 10 meter v wind speed real (kind=kind_phys), pointer :: sppt_wts(:,:) => null() !< real (kind=kind_phys), pointer :: shum_wts(:,:) => null() !< + real (kind=kind_phys), pointer :: sfc_wts(:,:) => null() !< real (kind=kind_phys), pointer :: zmtnblck(:) => null() ! null() !< u momentum change due to physics real (kind=kind_phys), pointer :: dv3dt (:,:,:) => null() !< v momentum change due to physics @@ -1752,7 +1757,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: drain(:) => null() !< real (kind=kind_phys), pointer :: drain_in_m_sm1(:) => null() !< real (kind=kind_phys), pointer :: dtdt(:,:) => null() !< - real (kind=kind_phys), pointer :: dtdtc(:,:) => null() !< real (kind=kind_phys), pointer :: dtsfc1(:) => null() !< real (kind=kind_phys), pointer :: dtzm(:) => null() !< real (kind=kind_phys), pointer :: dt_mf(:,:) => null() !< @@ -3380,6 +3384,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- stochastic physics control parameters logical :: do_sppt = .false. + logical :: pert_mp = .false. + logical :: pert_clds = .false. + logical :: pert_radtend = .true. logical :: use_zmtnblck = .false. logical :: do_shum = .false. logical :: do_skeb = .false. @@ -3462,6 +3469,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & dlqf, rbcr, shoc_parm, psauras, prauras, wminras, & do_sppt, do_shum, do_skeb, & lndp_type, n_var_lndp, lndp_each_step, & + pert_mp,pert_clds,pert_radtend, & !--- Rayleigh friction prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, & ! --- Ferrier-Aligo @@ -4142,6 +4150,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! physics that are parsed here and then compared in init_stochastic_physics ! to the stochastic physics namelist parametersto ensure consistency. Model%do_sppt = do_sppt + Model%pert_mp = pert_mp + Model%pert_clds = pert_clds + Model%pert_radtend = pert_radtend Model%use_zmtnblck = use_zmtnblck Model%do_shum = do_shum Model%do_skeb = do_skeb @@ -5163,6 +5174,9 @@ subroutine control_print(Model) print *, ' ' print *, 'stochastic physics' print *, ' do_sppt : ', Model%do_sppt + print *, ' pert_mp : ', Model%pert_mp + print *, ' pert_clds : ', Model%pert_clds + print *, ' pert_radtend : ', Model%pert_radtend print *, ' do_shum : ', Model%do_shum print *, ' do_skeb : ', Model%do_skeb print *, ' lndp_type : ', Model%lndp_type @@ -5405,10 +5419,10 @@ subroutine tbd_create (Tbd, IM, Model) endif if (Model%do_sppt .or. Model%ca_global) then - allocate (Tbd%dtdtr (IM,Model%levs)) + allocate (Tbd%dtdtnp (IM,Model%levs)) allocate (Tbd%dtotprcp (IM)) allocate (Tbd%dcnvprcp (IM)) - Tbd%dtdtr = clear_val + Tbd%dtdtnp = clear_val Tbd%dtotprcp = clear_val Tbd%dcnvprcp = clear_val endif @@ -5680,7 +5694,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%skebv_wts(IM,Model%levs)) allocate (Diag%sppt_wts(IM,Model%levs)) allocate (Diag%shum_wts(IM,Model%levs)) - allocate (Diag%zmtnblck(IM)) + allocate (Diag%sfc_wts(IM,Model%n_var_lndp)) + allocate (Diag%zmtnblck(IM)) allocate (Diag%ca1 (IM)) allocate (Diag%ca2 (IM)) allocate (Diag%ca3 (IM)) @@ -6267,7 +6282,6 @@ subroutine interstitial_create (Interstitial, IM, Model) allocate (Interstitial%dqsfc1 (IM)) allocate (Interstitial%drain (IM)) allocate (Interstitial%dtdt (IM,Model%levs)) - allocate (Interstitial%dtdtc (IM,Model%levs)) allocate (Interstitial%dtsfc1 (IM)) allocate (Interstitial%dt_mf (IM,Model%levs)) allocate (Interstitial%dtzm (IM)) @@ -6961,7 +6975,6 @@ subroutine interstitial_phys_reset (Interstitial, Model) Interstitial%drain = clear_val Interstitial%dt_mf = clear_val Interstitial%dtdt = clear_val - Interstitial%dtdtc = clear_val Interstitial%dtsfc1 = clear_val Interstitial%dtzm = clear_val Interstitial%dudt = clear_val @@ -7313,7 +7326,6 @@ subroutine interstitial_print(Interstitial, Model, mpirank, omprank, blkno) write (0,*) 'sum(Interstitial%dqsfc1 ) = ', sum(Interstitial%dqsfc1 ) write (0,*) 'sum(Interstitial%drain ) = ', sum(Interstitial%drain ) write (0,*) 'sum(Interstitial%dtdt ) = ', sum(Interstitial%dtdt ) - write (0,*) 'sum(Interstitial%dtdtc ) = ', sum(Interstitial%dtdtc ) write (0,*) 'sum(Interstitial%dtsfc1 ) = ', sum(Interstitial%dtsfc1 ) write (0,*) 'sum(Interstitial%dtzm ) = ', sum(Interstitial%dtzm ) write (0,*) 'sum(Interstitial%dt_mf ) = ', sum(Interstitial%dt_mf ) diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta index 7dd7f040a..6643683e0 100644 --- a/ccpp/data/GFS_typedefs.meta +++ b/ccpp/data/GFS_typedefs.meta @@ -4267,6 +4267,31 @@ units = flag dimensions = () type = logical +[pert_mp] + standard_name = flag_for_stochastic_microphysics_perturbations + long_name = flag for stochastic microphysics perturbations + units = flag + dimensions = () + type = logical +[pert_clds] + standard_name = flag_for_stochastic_cloud_fraction_perturbations + long_name = flag for stochastic cloud fraction perturbations + units = flag + dimensions = () + type = logical +[sppt_amp] + standard_name = total_ampltiude_of_sppt_perturbation + long_name = toal ampltidue of stochastic sppt perturbation + units = none + dimensions = () + type = real + kind = kind_phys +[pert_radtend] + standard_name = flag_for_stochastic_radiative_heating_perturbations + long_name = flag for stochastic radiative heating perturbations + units = flag + dimensions = () + type = logical [use_zmtnblck] standard_name = flag_for_mountain_blocking long_name = flag for mountain blocking @@ -5386,10 +5411,10 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys -[dtdtr] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_on_physics_time_step - long_name = temp. change due to radiative heating per time step - units = K +[dtdtnp] + standard_name = tendency_of_air_temperature_to_withold_from_sppt + long_name = temp. change from physics that should not be perturbed by sppt + units = K s-1 dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys @@ -6504,6 +6529,13 @@ dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys +[sfc_wts] + standard_name = weights_for_stochastic_surface_physics_perturbation_flipped + long_name = weights for stochastic surface physics perturbation, flipped + units = none + dimensions = (horizontal_loop_extent,number_of_land_surface_variables_perturbed) + type = real + kind = kind_phys [zmtnblck] standard_name = level_of_dividing_streamline long_name = level of the dividing streamline @@ -8085,13 +8117,6 @@ dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys -[dtdtc] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_assuming_clear_sky - long_name = clear sky radiative (shortwave + longwave) heating rate at current time - units = K s-1 - dimensions = (horizontal_loop_extent,vertical_dimension) - type = real - kind = kind_phys [dtsfc1] standard_name = instantaneous_surface_upward_sensible_heat_flux long_name = surface upward sensible heat flux diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90 index f5a2d2105..fea122446 100644 --- a/ccpp/driver/GFS_diagnostics.F90 +++ b/ccpp/driver/GFS_diagnostics.F90 @@ -2033,6 +2033,28 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%shum_wts(:,:) enddo + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'sfc_wts1' + ExtDiag(idx)%desc = 'perturbation amplitude' + ExtDiag(idx)%unit = 'none' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sfc_wts(:,1) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'sfc_wts2' + ExtDiag(idx)%desc = 'perturbation amplitude' + ExtDiag(idx)%unit = 'none' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sfc_wts(:,2) + enddo + idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'ca1' diff --git a/ccpp/physics b/ccpp/physics index 7242a6de9..0f9b52d7c 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 7242a6de94466d49c9e9150c43b7c0e75a655bb0 +Subproject commit 0f9b52d7c7a832afd298ed37b1af766ced080a55 diff --git a/coupler_main.F90 b/coupler_main.F90 deleted file mode 100644 index 3f36fb550..000000000 --- a/coupler_main.F90 +++ /dev/null @@ -1,506 +0,0 @@ -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of fvGFS. * -!* * -!* fvGFS is free software; you can redistribute it and/or modify it * -!* and are expected to follow the terms of the GNU General Public * -!* License as published by the Free Software Foundation; either * -!* version 2 of the License, or (at your option) any later version. * -!* * -!* fvGFS is distributed in the hope that it will be useful, but * -!* WITHOUT ANY WARRANTY; without even the implied warranty of * -!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * -!* General Public License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** - -program coupler_main - -!----------------------------------------------------------------------- -! -! program that couples component models for the atmosphere, -! ocean (amip), land, and sea-ice using the exchange module -! -!----------------------------------------------------------------------- - -use time_manager_mod, only: time_type, set_calendar_type, set_time, & - set_date, days_in_month, month_name, & - operator(+), operator (<), operator (>), & - operator (/=), operator (/), operator (==),& - operator (*), THIRTY_DAY_MONTHS, JULIAN, & - NOLEAP, NO_CALENDAR, date_to_string, & - get_date - -use atmos_model_mod, only: atmos_model_init, atmos_model_end, & - update_atmos_model_dynamics, & - update_atmos_radiation_physics, & - update_atmos_model_state, & - atmos_data_type, atmos_model_restart - -use constants_mod, only: constants_init -#ifdef INTERNAL_FILE_NML -use mpp_mod, only: input_nml_file -#else -use fms_mod, only: open_namelist_file -#endif -use fms_mod, only: file_exist, check_nml_error, & - error_mesg, fms_init, fms_end, close_file, & - write_version_number, uppercase - -use mpp_mod, only: mpp_init, mpp_pe, mpp_root_pe, mpp_npes, mpp_get_current_pelist, & - mpp_set_current_pelist, stdlog, mpp_error, NOTE, FATAL, WARNING -use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, mpp_sync - -use mpp_io_mod, only: mpp_open, mpp_close, & - MPP_NATIVE, MPP_RDONLY, MPP_DELETE - -use mpp_domains_mod, only: mpp_get_global_domain, mpp_global_field, CORNER -use memutils_mod, only: print_memuse_stats -use sat_vapor_pres_mod,only: sat_vapor_pres_init - -use diag_manager_mod, only: diag_manager_init, diag_manager_end, & - get_base_date, diag_manager_set_time_end - -use data_override_mod, only: data_override_init - - -implicit none - -!----------------------------------------------------------------------- - -character(len=128) :: version = '$Id: coupler_main.F90,v 19.0.4.1.2.3 2014/09/09 23:51:59 Rusty.Benson Exp $' -character(len=128) :: tag = '$Name: ulm_201505 $' - -!----------------------------------------------------------------------- -!---- model defined-types ---- - - type (atmos_data_type) :: Atm - -!----------------------------------------------------------------------- -! ----- coupled model time ----- - - type (time_type) :: Time_atmos, Time_init, Time_end, & - Time_step_atmos, Time_step_ocean, & - Time_restart, Time_step_restart - integer :: num_cpld_calls, num_atmos_calls, nc, na, ret - -! ----- coupled model initial date ----- - - integer :: date_init(6) - integer :: calendar_type = -99 - -! ----- timing flags ----- - - integer :: initClock, mainClock, termClock - integer, parameter :: timing_level = 1 - -! ----- namelist ----- - integer, dimension(6) :: current_date = (/ 0, 0, 0, 0, 0, 0 /) - character(len=17) :: calendar = ' ' - logical :: force_date_from_namelist = .false. ! override restart values for date - integer :: months=0, days=0, hours=0, minutes=0, seconds=0 - integer :: dt_atmos = 0 - integer :: dt_ocean = 0 - integer :: restart_days = 0 - integer :: restart_secs = 0 - integer :: atmos_nthreads = 1 - logical :: memuse_verbose = .false. - logical :: use_hyper_thread = .false. - logical :: debug_affinity = .false. - integer :: ncores_per_node = 0 - - namelist /coupler_nml/ current_date, calendar, force_date_from_namelist, & - months, days, hours, minutes, seconds, & - dt_atmos, dt_ocean, atmos_nthreads, memuse_verbose, & - use_hyper_thread, ncores_per_node, debug_affinity, & - restart_secs, restart_days - -! ----- local variables ----- - character(len=32) :: timestamp - logical :: intrm_rst - -!####################################################################### - - call fms_init() - call mpp_init() - initClock = mpp_clock_id( 'Initialization' ) - call mpp_clock_begin (initClock) !nesting problem - - call fms_init - call constants_init - call sat_vapor_pres_init - - call coupler_init - call print_memuse_stats('after coupler init') - - call mpp_set_current_pelist() - call mpp_clock_end (initClock) !end initialization - mainClock = mpp_clock_id( 'Main loop' ) - termClock = mpp_clock_id( 'Termination' ) - call mpp_clock_begin(mainClock) !begin main loop - - do nc = 1, num_cpld_calls - - Time_atmos = Time_atmos + Time_step_atmos - - call update_atmos_model_dynamics (Atm) - - call update_atmos_radiation_physics (Atm) - - call update_atmos_model_state (Atm) - -!--- intermediate restart - if (intrm_rst) then - if ((nc /= num_cpld_calls) .and. (Time_atmos == Time_restart)) then - timestamp = date_to_string (Time_restart) - call atmos_model_restart(Atm, timestamp) - call coupler_res(timestamp) - Time_restart = Time_restart + Time_step_restart - endif - endif - - call print_memuse_stats('after full step') - - enddo - -!----------------------------------------------------------------------- - -#ifdef AVEC_TIMERS - call avec_timers_output -#endif - call mpp_set_current_pelist() - call mpp_clock_end(mainClock) - call mpp_clock_begin(termClock) - - call coupler_end - call mpp_set_current_pelist() - call mpp_clock_end(termClock) - - call fms_end - -!----------------------------------------------------------------------- - - stop - -contains - -!####################################################################### - - subroutine coupler_init - -!----------------------------------------------------------------------- -! initialize all defined exchange grids and all boundary maps -!----------------------------------------------------------------------- - integer :: total_days, total_seconds, unit, ierr, io - integer :: n, gnlon, gnlat - integer :: date(6), flags - type (time_type) :: Run_length - character(len=9) :: month - logical :: use_namelist - - logical, allocatable, dimension(:,:) :: mask - real, allocatable, dimension(:,:) :: glon_bnd, glat_bnd - integer :: omp_get_thread_num, get_cpu_affinity, base_cpu -!----------------------------------------------------------------------- -!----- initialization timing identifiers ---- - -!----- read namelist ------- -!----- for backwards compatibilty read from file coupler.nml ----- - -#ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=coupler_nml, iostat=io) - ierr = check_nml_error(io, 'coupler_nml') -#else - if (file_exist('input.nml')) then - unit = open_namelist_file () - else - call error_mesg ('program coupler', & - 'namelist file input.nml does not exist', FATAL) - endif - - ierr=1 - do while (ierr /= 0) - read (unit, nml=coupler_nml, iostat=io, end=10) - ierr = check_nml_error (io, 'coupler_nml') - enddo -10 call close_file (unit) -#endif - -!----- write namelist to logfile ----- - call write_version_number (version, tag) - if (mpp_pe() == mpp_root_pe()) write(stdlog(),nml=coupler_nml) - -!----- allocate and set the pelist (to the global pelist) ----- - allocate( Atm%pelist (mpp_npes()) ) - call mpp_get_current_pelist(Atm%pelist) - -!----- read restart file ----- - - if (file_exist('INPUT/coupler.res')) then - call mpp_open( unit, 'INPUT/coupler.res', action=MPP_RDONLY ) - read (unit,*,err=999) calendar_type - read (unit,*) date_init - read (unit,*) date - goto 998 !back to fortran-4 - ! read old-style coupler.res - 999 call mpp_close (unit) - call mpp_open (unit, 'INPUT/coupler.res', action=MPP_RDONLY, form=MPP_NATIVE) - read (unit) calendar_type - read (unit) date - 998 call mpp_close(unit) - else - force_date_from_namelist = .true. - endif - -!----- use namelist value (either no restart or override flag on) --- - - if ( force_date_from_namelist ) then - - if ( sum(current_date) <= 0 ) then - call error_mesg ('program coupler', & - 'no namelist value for current_date', FATAL) - else - date = current_date - endif - -!----- override calendar type with namelist value ----- - - select case( uppercase(trim(calendar)) ) - case( 'JULIAN' ) - calendar_type = JULIAN - case( 'NOLEAP' ) - calendar_type = NOLEAP - case( 'THIRTY_DAY' ) - calendar_type = THIRTY_DAY_MONTHS - case( 'NO_CALENDAR' ) - calendar_type = NO_CALENDAR - case default - call mpp_error ( FATAL, 'COUPLER_MAIN: coupler_nml entry calendar must '// & - 'be one of JULIAN|NOLEAP|THIRTY_DAY|NO_CALENDAR.' ) - end select - - endif - -!$ base_cpu = get_cpu_affinity() -!$ call omp_set_num_threads(atmos_nthreads) -!$OMP PARALLEL NUM_THREADS(atmos_nthreads) -!$ if(omp_get_thread_num() < atmos_nthreads/2 .OR. (.not. use_hyper_thread)) then -!$ call set_cpu_affinity(base_cpu + omp_get_thread_num()) -!$ else -!$ call set_cpu_affinity(base_cpu + omp_get_thread_num() + & -!$ ncores_per_node - atmos_nthreads/2) -!$ endif -!$ if (debug_affinity) then -!$ write(6,*) mpp_pe()," atmos ",get_cpu_affinity(), base_cpu, omp_get_thread_num() -!$ call flush(6) -!$ endif -!$OMP END PARALLEL - - call set_calendar_type (calendar_type) - -!----- write current/initial date actually used to logfile file ----- - - if ( mpp_pe() == mpp_root_pe() ) then - write (stdlog(),16) date(1),trim(month_name(date(2))),date(3:6) - endif - - 16 format (' current date used = ',i4,1x,a,2i3,2(':',i2.2),' gmt') - -!----------------------------------------------------------------------- -!------ initialize diagnostics manager ------ - - call diag_manager_init (TIME_INIT=date) - -!----- always override initial/base date with diag_manager value ----- - - call get_base_date ( date_init(1), date_init(2), date_init(3), & - date_init(4), date_init(5), date_init(6) ) - -!----- use current date if no base date ------ - - if ( date_init(1) == 0 ) date_init = date - -!----- set initial and current time types ------ - - Time_init = set_date (date_init(1), date_init(2), date_init(3), & - date_init(4), date_init(5), date_init(6)) - - Time_atmos = set_date (date(1), date(2), date(3), & - date(4), date(5), date(6)) - -!----------------------------------------------------------------------- -!----- compute the ending time (compute days in each month first) ----- -! -! (NOTE: if run length in months then starting day must be <= 28) - - if ( months > 0 .and. date(3) > 28 ) & - call error_mesg ('program coupler', & - 'if run length in months then starting day must be <= 28', FATAL) - - Time_end = Time_atmos - total_days = 0 - do n = 1, months - total_days = total_days + days_in_month(Time_end) - Time_end = Time_atmos + set_time (0,total_days) - enddo - - total_days = total_days + days - total_seconds = hours*3600 + minutes*60 + seconds - Run_length = set_time (total_seconds,total_days) - Time_end = Time_atmos + Run_length - - !Need to pass Time_end into diag_manager for multiple thread case. - call diag_manager_set_time_end(Time_end) - - -!----------------------------------------------------------------------- -!----- write time stamps (for start time and end time) ------ - - call mpp_open( unit, 'time_stamp.out', nohdrs=.TRUE. ) - - month = month_name(date(2)) - if ( mpp_pe() == mpp_root_pe() ) write (unit,20) date, month(1:3) - - call get_date (Time_end, date(1), date(2), date(3), & - date(4), date(5), date(6)) - month = month_name(date(2)) - if ( mpp_pe() == mpp_root_pe() ) write (unit,20) date, month(1:3) - - call mpp_close (unit) - - 20 format (6i4,2x,a3) - -!----------------------------------------------------------------------- -!----- compute the time steps ------ - -Time_step_atmos = set_time (dt_atmos,0) -Time_step_ocean = set_time (dt_ocean,0) -num_cpld_calls = Run_length / Time_step_ocean -num_atmos_calls = Time_step_ocean / Time_step_atmos -Time_step_restart = set_time (restart_secs, restart_days) -Time_restart = Time_atmos + Time_step_restart -intrm_rst = .false. -if (restart_days > 0 .or. restart_secs > 0) intrm_rst = .true. - -!----------------------------------------------------------------------- -!------------------- some error checks --------------------------------- - -!----- initial time cannot be greater than current time ------- - - if ( Time_init > Time_atmos ) call error_mesg ('program coupler', & - 'initial time is greater than current time', FATAL) - -!----- make sure run length is a multiple of ocean time step ------ - - if ( num_cpld_calls * Time_step_ocean /= Run_length ) & - call error_mesg ('program coupler', & - 'run length must be multiple of ocean time step', FATAL) - -! ---- make sure cpld time step is a multiple of atmos time step ---- - - if ( num_atmos_calls * Time_step_atmos /= Time_step_ocean ) & - call error_mesg ('program coupler', & - 'atmos time step is not a multiple of the ocean time step', FATAL) - -!------ initialize component models ------ - - call atmos_model_init (Atm, Time_init, Time_atmos, Time_step_atmos) - - call print_memuse_stats('after atmos model init') - - call mpp_get_global_domain(Atm%Domain, xsize=gnlon, ysize=gnlat) - allocate ( glon_bnd(gnlon+1,gnlat+1), glat_bnd(gnlon+1,gnlat+1) ) - call mpp_global_field(Atm%Domain, Atm%lon_bnd, glon_bnd, position=CORNER) - call mpp_global_field(Atm%Domain, Atm%lat_bnd, glat_bnd, position=CORNER) - - call data_override_init ( ) ! Atm_domain_in = Atm%domain, & - ! Ice_domain_in = Ice%domain, & - ! Land_domain_in = Land%domain ) - -!----------------------------------------------------------------------- -!---- open and close dummy file in restart dir to check if dir exists -- - - if (mpp_pe() == 0 ) then - call mpp_open( unit, 'RESTART/file' ) - call mpp_close(unit, MPP_DELETE) - endif - -!----------------------------------------------------------------------- - - end subroutine coupler_init - -!####################################################################### - subroutine coupler_res(timestamp) - character(len=32), intent(in) :: timestamp - - integer :: unit, date(6) - -!----- compute current date ------ - - call get_date (Time_atmos, date(1), date(2), date(3), & - date(4), date(5), date(6)) - -!----- write restart file ------ - - if (mpp_pe() == mpp_root_pe())then - call mpp_open( unit, 'RESTART/'//trim(timestamp)//'.coupler.res', nohdrs=.TRUE. ) - write( unit, '(i6,8x,a)' )calendar_type, & - '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - - write( unit, '(6i6,8x,a)' )date_init, & - 'Model start time: year, month, day, hour, minute, second' - write( unit, '(6i6,8x,a)' )date, & - 'Current model time: year, month, day, hour, minute, second' - call mpp_close(unit) - endif - end subroutine coupler_res - -!####################################################################### - - subroutine coupler_end - - integer :: unit, date(6) -!----------------------------------------------------------------------- - - call atmos_model_end (Atm) - -!----- compute current date ------ - - call get_date (Time_atmos, date(1), date(2), date(3), & - date(4), date(5), date(6)) - -!----- check time versus expected ending time ---- - - if (Time_atmos /= Time_end) call error_mesg ('program coupler', & - 'final time does not match expected ending time', WARNING) - -!----- write restart file ------ - - call mpp_open( unit, 'RESTART/coupler.res', nohdrs=.TRUE. ) - if (mpp_pe() == mpp_root_pe())then - write( unit, '(i6,8x,a)' )calendar_type, & - '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - - write( unit, '(6i6,8x,a)' )date_init, & - 'Model start time: year, month, day, hour, minute, second' - write( unit, '(6i6,8x,a)' )date, & - 'Current model time: year, month, day, hour, minute, second' - endif - call mpp_close(unit) - -!----- final output of diagnostic fields ---- - - call diag_manager_end (Time_atmos) - -!----------------------------------------------------------------------- - - end subroutine coupler_end - -!####################################################################### - -end program coupler_main - diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90 index 9c7b64bf1..776a89bf4 100644 --- a/module_fcst_grid_comp.F90 +++ b/module_fcst_grid_comp.F90 @@ -73,6 +73,7 @@ module module_fcst_grid_comp quilting, calendar_type, cpl, & cplprint_flag, force_date_from_configure, & num_restart_interval, frestart, restart_endfcst + use get_stochy_pattern_mod, only: write_stoch_restart_atm ! !----------------------------------------------------------------------- ! @@ -827,6 +828,7 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) atm_int_state%Time_restart = atm_int_state%Time_atstart + restart_inctime timestamp = date_to_string (atm_int_state%Time_restart) call atmos_model_restart(atm_int_state%Atm, timestamp) + call write_stoch_restart_atm('RESTART/'//trim(timestamp)//'.atm_stoch.res.nc') call wrt_atmres_timestamp(atm_int_state,timestamp) endif diff --git a/stochastic_physics/stochastic_physics_wrapper.F90 b/stochastic_physics/stochastic_physics_wrapper.F90 index 479270a9f..0d0b29164 100644 --- a/stochastic_physics/stochastic_physics_wrapper.F90 +++ b/stochastic_physics/stochastic_physics_wrapper.F90 @@ -95,9 +95,15 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) initalize_stochastic_physics: if (GFS_Control%kdt==0) then if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .GT. 0) ) then + allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz))) + allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz))) + do nb=1,Atm_block%nblks + xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:) + xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:) + end do ! Initialize stochastic physics - call init_stochastic_physics(GFS_Control%levs, GFS_Control%blksz, GFS_Control%dtp, & - GFS_Control%input_nml_file, GFS_Control%fn_nml, GFS_Control%nlunit, GFS_Control%do_sppt, GFS_Control%do_shum, & + call init_stochastic_physics(GFS_Control%levs, GFS_Control%blksz, GFS_Control%dtp, GFS_Control%sppt_amp, & + GFS_Control%input_nml_file, GFS_Control%fn_nml, GFS_Control%nlunit, xlon, xlat, GFS_Control%do_sppt, GFS_Control%do_shum, & GFS_Control%do_skeb, GFS_Control%lndp_type, GFS_Control%n_var_lndp, GFS_Control%use_zmtnblck, GFS_Control%skeb_npass, & GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, & GFS_Control%ak, GFS_Control%bk, nthreads, GFS_Control%master, GFS_Control%communicator, ierr) @@ -106,8 +112,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) return endif end if - allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz))) - allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz))) if (GFS_Control%do_sppt) then allocate(sppt_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs)) end if @@ -143,17 +147,13 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) allocate(zorll(1:Atm_block%nblks,maxval(GFS_Control%blksz))) endif - do nb=1,Atm_block%nblks - xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:) - xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:) - end do if ( GFS_Control%lndp_type .EQ. 1 ) then ! this scheme sets perts once allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%n_var_lndp)) - call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, & + call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%fhour, GFS_Control%blksz, & sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, & nthreads=nthreads) - ! Copy contiguous data back; no need to copy xlat/xlon, these are intent(in) - just deallocate + ! Copy contiguous data back do nb=1,Atm_block%nblks GFS_Data(nb)%Coupling%sfc_wts(:,:) = sfc_wts(nb,1:GFS_Control%blksz(nb),:) end do @@ -170,12 +170,11 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) endif else initalize_stochastic_physics - if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .EQ. 2) ) then - call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, & + call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%fhour, GFS_Control%blksz, & sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, & nthreads=nthreads) - ! Copy contiguous data back; no need to copy xlat/xlon, these are intent(in) - just deallocate + ! Copy contiguous data back if (GFS_Control%do_sppt) then do nb=1,Atm_block%nblks GFS_Data(nb)%Coupling%sppt_wts(:,:) = sppt_wts(nb,1:GFS_Control%blksz(nb),:) @@ -239,6 +238,7 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) GFS_Control%n_var_lndp, GFS_Control%lndp_var_list, GFS_Control%lndp_prt_list, & sfc_wts, xlon, xlat, stype, GFS_Control%pores, GFS_Control%resid,param_update_flag, & smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, snoalb, semis, zorll, ierr) + if (ierr/=0) then write(6,*) 'call to GFS_apply_lndp failed' return