From a231b93f943f1820203d888d197bf1fd70f017c8 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 1 Aug 2018 11:47:21 -0600 Subject: [PATCH 1/4] Remove legacy physics/GFS_suite_setup_scm.F90 --- physics/GFS_suite_setup_scm.F90 | 263 -------------------------------- 1 file changed, 263 deletions(-) delete mode 100644 physics/GFS_suite_setup_scm.F90 diff --git a/physics/GFS_suite_setup_scm.F90 b/physics/GFS_suite_setup_scm.F90 deleted file mode 100644 index 7dc7f15c2..000000000 --- a/physics/GFS_suite_setup_scm.F90 +++ /dev/null @@ -1,263 +0,0 @@ -module GFS_suite_setup_scm - - implicit none - - private - -!---------------- -! Public entities -!---------------- - public GFS_suite_setup_scm_init, GFS_suite_setup_scm_run, GFS_suite_setup_scm_finalize - - CONTAINS -!******************************************************************************************* - -!-------------- -! GFS initialze -!-------------- - -!> \section arg_table_GFS_suite_setup_scm_init Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------------|-------------------------------------------------------------|-------------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| -!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_control_type | | inout | F | -!! | Statein | FV3-GFS_Statein_type | Fortran DDT containing FV3-GFS prognostic state data in from dycore | DDT | 0 | GFS_statein_type | | inout | F | -!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_stateout_type | | inout | F | -!! | Sfcprop | FV3-GFS_Sfcprop_type | Fortran DDT containing FV3-GFS surface fields | DDT | 0 | GFS_sfcprop_type | | inout | F | -!! | Coupling | FV3-GFS_Coupling_type | derived type GFS_coupling_type in FV3 | DDT | 0 | GFS_coupling_type | | inout | F | -!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_grid_type | | inout | F | -!! | Tbd | FV3-GFS_Tbd_type | derived type GFS_tbd_type in FV3 | DDT | 0 | GFS_tbd_type | | inout | F | -!! | Cldprop | FV3-GFS_Cldprop_type | derived type GFS_cldprop_type in FV3 | DDT | 0 | GFS_cldprop_type | | inout | F | -!! | Radtend | FV3-GFS_Radtend_type | derived type GFS_radtend_type in FV3 | DDT | 0 | GFS_radtend_type | | inout | F | -!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_diag_type | | inout | F | -!! | Interstitial | FV3-GFS_Interstitial_type | derived type GFS_interstitial_type in FV3 | DDT | 0 | GFS_interstitial_type | | inout | F | -!! | Init_parm | FV3-GFS_Init_type | dervied type GFS_init_type in FV3 | DDT | 0 | GFS_init_type | | in | F | -!! | n_ozone_layers | vertical_dimension_of_ozone_forcing_data_from_host | number of vertical layers in ozone forcing data coming from host | count | 0 | integer | | in | F | -!! | n_ozone_lats | number_of_latitutde_points_in_ozone_forcing_data_from_host | number of latitude points in ozone forcing data coming from host | count | 0 | integer | | in | F | -!! | n_ozone_times | number_of_time_levels_in_ozone_forcing_data_from_host | number of time levels in ozone forcing data coming from host | count | 0 | integer | | in | F | -!! | n_ozone_coefficients | number_of_coefficients_in_ozone_forcing_data_from_host | number of coeffcients in ozone forcing data coming from host | count | 0 | integer | | in | F | -!! | ozone_lat | latitude_of_ozone_forcing_data_from_host | latitude value of the ozone forcing data coming from host | degree | 1 | real | kind_phys | in | F | -!! | ozone_pres | natural_log_of_ozone_forcing_data_pressure_levels_from_host | natural logarithm of the pressure levels of the ozone forcing data | Pa | 1 | real | kind_phys | in | F | -!! | ozone_time | time_levels_in_ozone_forcing_data_from_host | time values of the ozone forcing data coming from host | day | 1 | real | kind_phys | in | F | -!! | ozone_forcing_in | ozone_forcing_from_host | ozone forcing data from host | various | 4 | real | kind_phys | in | F | -!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | -!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | -!! - subroutine GFS_suite_setup_scm_init (Model, Statein, Stateout, Sfcprop, & - Coupling, Grid, Tbd, Cldprop, Radtend, Diag, & - Interstitial, Init_parm, n_ozone_lats, & - n_ozone_layers, n_ozone_times, n_ozone_coefficients, & - ozone_lat, ozone_pres, ozone_time, ozone_forcing_in, & - errmsg, errflg) - - use machine, only: kind_phys - use GFS_typedefs, only: GFS_init_type, & - GFS_statein_type, GFS_stateout_type, & - GFS_sfcprop_type, GFS_coupling_type, & - GFS_control_type, GFS_grid_type, & - GFS_tbd_type, GFS_cldprop_type, & - GFS_radtend_type, GFS_diag_type, & - GFS_interstitial_type - use funcphys, only: gfuncphys - use module_microphysics, only: gsmconst - use cldwat2m_micro, only: ini_micro - use aer_cloud, only: aer_cloud_init - use module_ras, only: ras_init - use ozne_def, only: latsozp, levozp, timeoz, oz_coeff, oz_lat, oz_pres, oz_time, ozplin - use GFS_rrtmg_setup, only: GFS_rrtmg_setup_init - - !--- interface variables - type(GFS_control_type), intent(inout) :: Model - type(GFS_statein_type), intent(inout) :: Statein - type(GFS_stateout_type), intent(inout) :: Stateout - type(GFS_sfcprop_type), intent(inout) :: Sfcprop - type(GFS_coupling_type), intent(inout) :: Coupling - type(GFS_grid_type), intent(inout) :: Grid - type(GFS_tbd_type), intent(inout) :: Tbd - type(GFS_cldprop_type), intent(inout) :: Cldprop - type(GFS_radtend_type), intent(inout) :: Radtend - type(GFS_diag_type), intent(inout) :: Diag - type(GFS_interstitial_type), intent(inout) :: Interstitial - type(GFS_init_type), intent(in) :: Init_parm - - integer, intent(in) :: n_ozone_lats, n_ozone_layers, n_ozone_coefficients, n_ozone_times - real(kind=kind_phys), intent(in) :: ozone_lat(:), ozone_pres(:), ozone_time(:), ozone_forcing_in(:,:,:,:) - - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - ! - ! !--- local variables - real(kind=kind_phys), allocatable :: si(:) - real(kind=kind_phys), parameter :: p_ref = 101325.0d0 - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - -! !--- set control properties (including namelist read) - call Model%init (Init_parm%nlunit, Init_parm%fn_nml, & - Init_parm%me, Init_parm%master, & - Init_parm%logunit, Init_parm%isc, & - Init_parm%jsc, Init_parm%nx, Init_parm%ny, & - Init_parm%levs, Init_parm%cnx, Init_parm%cny, & - Init_parm%gnx, Init_parm%gny, & - Init_parm%dt_dycore, Init_parm%dt_phys, & - Init_parm%bdat, Init_parm%cdat, & - Init_parm%tracer_names, Init_parm%blksz) - - !allocate memory for the variables stored in ozne_def and set them - allocate(oz_lat(n_ozone_lats), oz_pres(n_ozone_layers), oz_time(n_ozone_times+1)) - allocate(ozplin(n_ozone_lats, n_ozone_layers, n_ozone_coefficients, n_ozone_times)) - latsozp = n_ozone_lats - levozp = n_ozone_layers - timeoz = n_ozone_times - oz_coeff = n_ozone_coefficients - oz_lat = ozone_lat - oz_pres = ozone_pres - oz_time = ozone_time - ozplin = ozone_forcing_in - - call Statein%create(1, Model) - call Stateout%create(1, Model) - call Sfcprop%create(1, Model) - call Coupling%create(1, Model) - call Grid%create(1, Model) - call Tbd%create(1, 1, Model) - call Cldprop%create(1, Model) - call Radtend%create(1, Model) - !--- internal representation of diagnostics - call Diag%create(1, Model) - !--- internal representation of interstitials for CCPP physics - call Interstitial%create(1, Model) - -! !--- populate the grid components - call GFS_grid_populate (Grid, Init_parm%xlon, Init_parm%xlat, Init_parm%area) - - !--- read in and initialize ozone and water - if (Model%ntoz > 0) then - call setindxoz (Init_parm%blksz, Grid%xlat_d, Grid%jindx1_o3, & - Grid%jindx2_o3, Grid%ddy_o3) - endif - - if (Model%h2o_phys) then - call setindxh2o (Init_parm%blksz, Grid%xlat_d, Grid%jindx1_h, & - Grid%jindx2_h, Grid%ddy_h) - endif - -! !--- Call gfuncphys (funcphys.f) to compute all physics function tables. - call gfuncphys () -! - call gsmconst (Model%dtp, Model%me, .TRUE.) -! -! !--- define sigma level for radiation initialization -! !--- The formula converting hybrid sigma pressure coefficients to sigma coefficients follows Eckermann (2009, MWR) -! !--- ps is replaced with p0. The value of p0 uses that in http://www.emc.ncep.noaa.gov/officenotes/newernotes/on461.pdf -! !--- ak/bk have been flipped from their original FV3 orientation and are defined sfc -> toa - allocate(si(Model%levr+1)) - si = (Init_parm%ak + Init_parm%bk * p_ref - Init_parm%ak(Model%levr+1)) & - / (p_ref - Init_parm%ak(Model%levr+1)) - call GFS_rrtmg_setup_init (si, Model%levr, Model%ictm, Model%isol, & - Model%ico2, Model%iaer, Model%ialb, Model%iems, & - Model%ntcw, Model%num_p2d, Model%num_p3d, Model%npdf3d, & - Model%ntoz, Model%iovr_sw, Model%iovr_lw, Model%isubc_sw, & - Model%isubc_lw, Model%crick_proof, Model%ccnorm, & - Model%imp_physics, Model%norad_precip, Model%idate, & - Model%iflip, Interstitial%im, Interstitial%faerlw, & - Interstitial%faersw, Interstitial%aerodp, & - Model%me, errmsg, errflg) - if (errflg/=0) then - print *, "An error occured in GFS_rrtmg_setup_init: " // trim(errmsg) - stop - endif - deallocate (si) -! -! !--- initialize Morrison-Gettleman microphysics - if (Model%ncld == 2) then - call ini_micro (Model%mg_dcs, Model%mg_qcvar, Model%mg_ts_auto_ice) - call aer_cloud_init () - endif -! -! !--- initialize ras - if (Model%ras) call ras_init (Model%levs, Model%me) -! -! !--- initialize soil vegetation - call set_soilveg(Model%me, Model%isot, Model%ivegsrc, Model%nlunit) -! -! !--- lsidea initialization - if (Model%lsidea) then - print *,' LSIDEA is active but needs to be reworked for FV3 - shutting down' - stop - !--- NEED TO get the logic from the old phys/gloopb.f initialization area - endif -! -! !--- sncovr may not exist in ICs from chgres. -! !--- FV3GFS handles this as part of the IC ingest -! !--- this not is placed here to alert users to the need to study -! !--- the FV3GFS_io.F90 module - - end subroutine GFS_suite_setup_scm_init - -!> \section arg_table_GFS_suite_setup_scm_run Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------------|-------------------------------------------------------------|-------------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| -!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | -!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | -!! - subroutine GFS_suite_setup_scm_run (errmsg, errflg) - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - end subroutine GFS_suite_setup_scm_run - -!> \section arg_table_GFS_suite_setup_scm_finalize Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------------|-------------------------------------------------------------|-------------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| -!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | -!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | -!! - subroutine GFS_suite_setup_scm_finalize (errmsg, errflg) - - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - - end subroutine GFS_suite_setup_scm_finalize - - !------------------ - ! GFS_grid_populate - !------------------ - subroutine GFS_grid_populate (Grid, xlon, xlat, area) - use machine, only: kind_phys - use physcons, only: pi => con_pi - use GFS_typedefs, only: GFS_grid_type - - implicit none - - type(GFS_grid_type) :: Grid - real(kind=kind_phys), intent(in) :: xlon(:,:) - real(kind=kind_phys), intent(in) :: xlat(:,:) - real(kind=kind_phys), intent(in) :: area(:,:) - - !--- local variables - integer :: n_columns, i - - n_columns = size(Grid%xlon) - - do i=1, n_columns - Grid%xlon = xlon(i,1) - Grid%xlat = xlat(i,1) - Grid%xlat_d(i) = xlat(i,1) * 180.0_kind_phys/pi - Grid%sinlat(i) = sin(Grid%xlat(i)) - Grid%coslat(i) = sqrt(1.0_kind_phys - Grid%sinlat(i)*Grid%sinlat(i)) - Grid%area(i) = area(i,1) - Grid%dx(i) = sqrt(area(i,1)) - end do - - end subroutine GFS_grid_populate - -end module GFS_suite_setup_scm From e9ace0a1fa4f95a2a5ea03df42feac7f72b9fc6c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 2 Aug 2018 11:39:31 -0600 Subject: [PATCH 2/4] physics/GFS_PBL_generic.f90: filter metadata tables through preprocessor --- physics/GFS_PBL_generic.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/physics/GFS_PBL_generic.f90 b/physics/GFS_PBL_generic.f90 index caaae694e..3c8a79b67 100644 --- a/physics/GFS_PBL_generic.f90 +++ b/physics/GFS_PBL_generic.f90 @@ -12,6 +12,7 @@ subroutine GFS_PBL_generic_pre_finalize() end subroutine GFS_PBL_generic_pre_finalize !> \brief This scheme sets up the vertically diffused tracer array for any PBL scheme based on the microphysics scheme chosen +#if 0 !! \section arg_table_GFS_PBL_generic_pre_run Argument Table !! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | !! |------------------------------|--------------------------------------------------------|-------------------------------------------------------------------------------------|---------------|------|-----------|-----------|--------|----------| @@ -40,6 +41,7 @@ end subroutine GFS_PBL_generic_pre_finalize !! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! +#endif subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_wsm6, ltaerosol, qgrs, qgrs_water_vapor, qgrs_liquid_cloud, qgrs_ice_cloud, qgrs_ozone, & qgrs_cloud_droplet_num_conc, qgrs_cloud_ice_num_conc, qgrs_water_aer_num_conc, qgrs_ice_aer_num_conc, qgrs_rain, & @@ -139,7 +141,7 @@ subroutine GFS_PBL_generic_post_finalize () end subroutine GFS_PBL_generic_post_finalize - +#if 0 !> \section arg_table_GFS_PBL_generic_post_run Argument Table !! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | !! |------------------------------|-----------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------|---------------|------|-----------|-----------|--------|----------| @@ -208,6 +210,7 @@ end subroutine GFS_PBL_generic_post_finalize !! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! +#endif subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, ntoz, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_wsm6, ltaerosol, cplflx, lssav, ldiag3d, lsidea, hybedmf, dvdftra, dusfc1, dvsfc1, dtsfc1, dqsfc1, dtf, & dudt, dvdt, dtdt, htrsw, htrlw, xmu,& From 7ca4e99992f2e829f692b28057cd8eecd757e83e Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 2 Aug 2018 11:45:14 -0600 Subject: [PATCH 3/4] Add stochastic_physics directory containing FV3/stochastic_physics with bugfixes and CCPP modifications --- stochastic_physics/compns_stochy.F90 | 209 +++++++ stochastic_physics/dezouv_stochy.f | 269 +++++++++ stochastic_physics/dozeuv_stochy.f | 267 +++++++++ stochastic_physics/epslon_stochy.f | 93 ++++ stochastic_physics/four_to_grid_stochy.F | 271 +++++++++ stochastic_physics/function2 | 5 + stochastic_physics/function_indlsev | 3 + stochastic_physics/function_indlsod | 3 + stochastic_physics/get_lats_node_a_stochy.f | 92 ++++ stochastic_physics/get_ls_node_stochy.f | 81 +++ stochastic_physics/get_stochy_pattern.F90 | 518 ++++++++++++++++++ stochastic_physics/getcon_lag_stochy.f | 89 +++ stochastic_physics/getcon_spectral.F90 | 277 ++++++++++ stochastic_physics/glats_stochy.f | 109 ++++ stochastic_physics/gozrineo_stochy.f | 180 ++++++ .../initialize_spectral_mod.F90 | 268 +++++++++ stochastic_physics/pln2eo_stochy.f | 287 ++++++++++ stochastic_physics/setlats_a_stochy.f | 195 +++++++ stochastic_physics/setlats_lag_stochy.f | 127 +++++ stochastic_physics/spectral_layout.f | 31 ++ stochastic_physics/stochastic_physics.F90 | 410 ++++++++++++++ stochastic_physics/stochy_ccpp.F90 | 438 +++++++++++++++ stochastic_physics/stochy_data_mod.F90 | 393 +++++++++++++ stochastic_physics/stochy_gg_def.f | 9 + .../stochy_internal_state_mod.F90 | 136 +++++ stochastic_physics/stochy_layout_lag.f | 13 + stochastic_physics/stochy_namelist_def.F90 | 37 ++ .../stochy_patterngenerator.F90 | 357 ++++++++++++ stochastic_physics/stochy_resol_def.f | 44 ++ stochastic_physics/sumfln_stochy.f | 294 ++++++++++ 30 files changed, 5505 insertions(+) create mode 100644 stochastic_physics/compns_stochy.F90 create mode 100644 stochastic_physics/dezouv_stochy.f create mode 100644 stochastic_physics/dozeuv_stochy.f create mode 100644 stochastic_physics/epslon_stochy.f create mode 100644 stochastic_physics/four_to_grid_stochy.F create mode 100644 stochastic_physics/function2 create mode 100644 stochastic_physics/function_indlsev create mode 100644 stochastic_physics/function_indlsod create mode 100644 stochastic_physics/get_lats_node_a_stochy.f create mode 100644 stochastic_physics/get_ls_node_stochy.f create mode 100644 stochastic_physics/get_stochy_pattern.F90 create mode 100644 stochastic_physics/getcon_lag_stochy.f create mode 100644 stochastic_physics/getcon_spectral.F90 create mode 100644 stochastic_physics/glats_stochy.f create mode 100644 stochastic_physics/gozrineo_stochy.f create mode 100644 stochastic_physics/initialize_spectral_mod.F90 create mode 100644 stochastic_physics/pln2eo_stochy.f create mode 100644 stochastic_physics/setlats_a_stochy.f create mode 100644 stochastic_physics/setlats_lag_stochy.f create mode 100644 stochastic_physics/spectral_layout.f create mode 100644 stochastic_physics/stochastic_physics.F90 create mode 100644 stochastic_physics/stochy_ccpp.F90 create mode 100644 stochastic_physics/stochy_data_mod.F90 create mode 100644 stochastic_physics/stochy_gg_def.f create mode 100644 stochastic_physics/stochy_internal_state_mod.F90 create mode 100644 stochastic_physics/stochy_layout_lag.f create mode 100644 stochastic_physics/stochy_namelist_def.F90 create mode 100644 stochastic_physics/stochy_patterngenerator.F90 create mode 100644 stochastic_physics/stochy_resol_def.f create mode 100644 stochastic_physics/sumfln_stochy.f diff --git a/stochastic_physics/compns_stochy.F90 b/stochastic_physics/compns_stochy.F90 new file mode 100644 index 000000000..29bd79b27 --- /dev/null +++ b/stochastic_physics/compns_stochy.F90 @@ -0,0 +1,209 @@ +module compns_stochy_mod + + implicit none + + contains + +!----------------------------------------------------------------------- + subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) +!$$$ Subprogram Documentation Block +! +! Subprogram: compns Check and compute namelist frequencies +! Prgmmr: Iredell Org: NP23 Date: 1999-01-26 +! +! Abstract: This subprogram checks global spectral model namelist +! frequencies in hour units for validity. If they are valid, +! then the frequencies are computed in timestep units. +! The following rules are applied: +! 1. the timestep must be positive; +! +! Program History Log: +! 2016-10-11 Phil Pegion make the stochastic physics stand alone +! +! Usage: call compns_stochy (me,deltim,nlunit, stochy_namelist,iret) +! Input Arguments: +! deltim - real timestep in seconds +! Output Arguments: +! iret - integer return code (0 if successful or +! between 1 and 8 for which rule above was broken) +! stochy_namelist +! +! Attributes: +! Language: Fortran 90 +! +!$$$ + + + use stochy_namelist_def + + implicit none + + + integer, intent(out) :: iret + integer, intent(in) :: nlunit,me,sz_nml + character(len=*), intent(in) :: input_nml_file(sz_nml) + character(len=64), intent(in) :: fn_nml + real, intent(in) :: deltim + real tol + integer k,ios + +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! + namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, & + iseed_shum,iseed_sppt,shum,shum_tau,& + shum_lscale,fhstoch,stochini,skeb_varspect_opt,sppt_sfclimit, & + skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & + skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& + shum_sigefold,skebint,skeb_npass,use_zmtnblck + namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts + pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land + + tol=0.01 ! tolerance for calculations +! spectral resolution defintion + ntrunc=-999 + lon_s=-999 + lat_s=-999 + ! can specify up to 5 values for the stochastic physics parameters + ! (each is an array of length 5) + sppt = -999. ! stochastic physics tendency amplitude + shum = -999. ! stochastic boundary layer spf hum amp + skeb = -999. ! stochastic KE backscatter amplitude + ! mg, sfcperts + pertz0 = -999. ! momentum roughness length amplitude + pertshc = -999. ! soil hydraulic conductivity amp + pertzt = -999. ! mom/heat roughness length amplitude + pertlai = -999. ! leaf area index amplitude + pertvegf = -999. ! vegetation fraction amplitude + pertalb = -999. ! albedo perturbations amplitude +! logicals + do_sppt = .false. + use_zmtnblck = .false. + do_shum = .false. + do_skeb = .false. + ! mg, sfcperts + do_sfcperts = .false. + sppt_land = .false. + nsfcpert = 0 +! for sfcperts random patterns + sfc_lscale = -999. ! length scales + sfc_tau = -999. ! time scales + iseed_sfc = 0 ! random seeds (if 0 use system clock) +! for SKEB random patterns. + skeb_vfilt = 0 + skebint = 0 + skeb_npass = 11 ! number of passes of smoother for dissipation estiamte + sppt_tau = -999. ! time scales + shum_tau = -999. + skeb_tau = -999. + skeb_vdof = 5 ! proxy for vertical correlation, 5 is close to 40 passes of the 1-2-1 filter in the GFS + skebnorm = 0 ! 0 - random pattern is stream function, 1- pattern is kenorm, 2- pattern is vorticity + sppt_lscale = -999. ! length scales + shum_lscale = -999. + skeb_lscale = -999. + iseed_sppt = 0 ! random seeds (if 0 use system clock) + iseed_shum = 0 + iseed_skeb = 0 +! parameters to control vertical tapering of stochastic physics with +! height + sppt_sigtop1 = 0.1 + sppt_sigtop2 = 0.025 + skeb_sigtop1 = 0.1 + skeb_sigtop2 = 0.025 + shum_sigefold = 0.2 +! reduce amplitude of sppt near surface (lowest 2 levels) + sppt_sfclimit = .false. +! gaussian or power law variance spectrum for skeb (0: gaussian, 1: +! power law). If power law, skeb_lscale interpreted as a power not a +! length scale. + skeb_varspect_opt = 0 + sppt_logit = .false. ! logit transform for sppt to bounded interval [-1,+1] + fhstoch = -999.0 ! forecast hour to dump random patterns + stochini = .false. ! true= read in pattern, false=initialize from seed + +#ifdef INTERNAL_FILE_NML + read(input_nml_file, nml=nam_stochy) +#else + rewind (nlunit) + open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios) + read(nlunit,nam_stochy) +#endif +#ifdef INTERNAL_FILE_NML + read(input_nml_file, nml=nam_sfcperts) +#else + rewind (nlunit) + open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios) + read(nlunit,nam_sfcperts) +#endif + + if (me == 0) then + print *,' in compns_stochy' + print*,'skeb=',skeb + endif + +! PJP stochastic physics additions + IF (sppt(1) > 0 ) THEN + do_sppt=.true. + ENDIF + IF (shum(1) > 0 ) THEN + do_shum=.true. +! shum parameter has units of 1/hour, to remove time step +! dependence. +! change shum parameter units from per hour to per timestep + DO k=1,5 + IF (shum(k) .gt. 0.0) shum(k)=shum(k)*deltim/3600.0 + ENDDO + ENDIF + IF (skeb(1) > 0 ) THEN + do_skeb=.true. + if (skebnorm==0) then ! stream function norm + skeb=skeb*1.111e3*sqrt(deltim) + !skeb=skeb*5.0e5/sqrt(deltim) + endif + if (skebnorm==1) then ! stream function norm + skeb=skeb*0.00222*sqrt(deltim) + !skeb=skeb*1/sqrt(deltim) + endif + if (skebnorm==2) then ! vorticty function norm + skeb=skeb*1.111e-9*sqrt(deltim) + !skeb=skeb*5.0e-7/sqrt(deltim) + endif +! adjust skeb values for resolution. +! scaling is such that a value of 1.0 at T574 with a 900 second +! timestep produces well-calibrated values of forecast spread. +! DO k=1,5 +! IF (skeb(k) .gt. 0.0) THEN +! skeb(k)=skeb(k)*deltim/(ntrunc*(ntrunc+1))*365765.0 ! 365765 is a scale factor so the base SKEB value in the namelist is 1.0 +! skeb(k)=skeb(k)*deltim/(ntrunc*(ntrunc+1))*2000.0 ! 2000 is new scale factor so the base SKEB value in the namelist is 1.0 +! ENDIF +! ENDDO + ENDIF +! compute frequencty to estimate dissipation timescale + IF (skebint == 0.) skebint=deltim + nsskeb=nint(skebint/deltim) ! skebint in seconds + IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN + WRITE(0,*) "SKEB interval is invalid",skebint + iret=9 + return + ENDIF +! mg, sfcperts + IF (pertz0(1) > 0 .OR. pertshc(1) > 0 .OR. pertzt(1) > 0 .OR. & + pertlai(1) > 0 .OR. pertvegf(1) > 0 .OR. pertalb(1) > 0) THEN + do_sfcperts=.true. + ENDIF +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! +! All checks are successful. +! + if (me == 0) then + print *, 'stochastic physics' + print *, ' do_sppt : ', do_sppt + print *, ' do_shum : ', do_shum + print *, ' do_skeb : ', do_skeb + print *, ' do_sfcperts : ', do_sfcperts + endif + iret = 0 +! + return + end subroutine compns_stochy + +end module compns_stochy_mod diff --git a/stochastic_physics/dezouv_stochy.f b/stochastic_physics/dezouv_stochy.f new file mode 100644 index 000000000..32f4e7dba --- /dev/null +++ b/stochastic_physics/dezouv_stochy.f @@ -0,0 +1,269 @@ + module dezouv_stochy_mod + + implicit none + + contains + + subroutine dezouv_stochy(dev,zod,uev,vod,epsedn,epsodn, + & snnp1ev,snnp1od,ls_node) +cc + +cc + use stochy_resol_def + use spectral_layout_mod + use machine + implicit none +cc + real(kind_dbl_prec) dev(len_trie_ls,2) + real(kind_dbl_prec) zod(len_trio_ls,2) + real(kind_dbl_prec) uev(len_trie_ls,2) + real(kind_dbl_prec) vod(len_trio_ls,2) +cc + real(kind_dbl_prec) epsedn(len_trie_ls) + real(kind_dbl_prec) epsodn(len_trio_ls) +cc + real(kind_dbl_prec) snnp1ev(len_trie_ls) + real(kind_dbl_prec) snnp1od(len_trio_ls) +cc + integer ls_node(ls_dim,3) +cc +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +cc + integer l,locl,n +cc + integer indev,indev1,indev2 + integer indod,indod1,indod2 + integer inddif +cc + real(kind_dbl_prec) rl +cc + real(kind_dbl_prec) cons0 !constant +cc + integer indlsev,jbasev + integer indlsod,jbasod + real(kind_evod) rerth +cc + include 'function2' +cc +cc +cc...................................................................... +cc +cc + cons0 = 0.d0 !constant + rerth =6.3712e+6 ! radius of earth (m) +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) +cc + uev(indlsev(l,l),1) = cons0 !constant + uev(indlsev(l,l),2) = cons0 !constant +cc +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + 1 + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + uev(indev,1) = -epsedn(indev) + x * zod(indev-inddif,1) +cc + uev(indev,2) = -epsedn(indev) + x * zod(indev-inddif,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + vod(indev-inddif,1) = epsodn(indev-inddif) + x * dev(indev,1) +cc + vod(indev-inddif,2) = epsodn(indev-inddif) + x * dev(indev,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + if ( l .ge. 1 ) then + rl = l + do indev = indev1 , indev2 +cc u(l,n)=-i*l*d(l,n)/(n*(n+1)) +cc + uev(indev,1) = uev(indev,1) + 1 + rl * dev(indev,2) + 2 / snnp1ev(indev) +cc + uev(indev,2) = uev(indev,2) + 1 - rl * dev(indev,1) + 2 / snnp1ev(indev) +cc + enddo + endif +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasod=ls_node(locl,3) + indod1 = indlsod(L+1,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indod2 = indlsod(jcap ,L) + else + indod2 = indlsod(jcap+1,L) - 1 + endif + if ( l .ge. 1 ) then + rl = l + do indod = indod1 , indod2 +cc u(l,n)=-i*l*d(l,n)/(n*(n+1)) +cc + vod(indod,1) = vod(indod,1) + 1 + rl * zod(indod,2) + 2 / snnp1od(indod) +cc + vod(indod,2) = vod(indod,2) + 1 - rl * zod(indod,1) + 2 / snnp1od(indod) +cc + enddo + endif +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) - 1 + else + indev2 = indlsev(jcap ,L) - 1 + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + uev(indev,1) = uev(indev ,1) + 1 + epsodn(indev-inddif) * zod(indev-inddif,1) +cc + uev(indev,2) = uev(indev ,2) + 1 + epsodn(indev-inddif) * zod(indev-inddif,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + 1 + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + vod(indev-inddif,1) = vod(indev-inddif,1) + 1 - epsedn(indev) * dev(indev ,1) +cc + vod(indev-inddif,2) = vod(indev-inddif,2) + 1 - epsedn(indev) * dev(indev ,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + indod1 = indlsod(L+1,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) + indod2 = indlsod(jcap ,L) + else + indev2 = indlsev(jcap ,L) + indod2 = indlsod(jcap+1,L) + endif + do indev = indev1 , indev2 +cc + uev(indev,1) = uev(indev,1) * rerth + uev(indev,2) = uev(indev,2) * rerth +cc + enddo +cc + do indod = indod1 , indod2 +cc + vod(indod,1) = vod(indod,1) * rerth + vod(indod,2) = vod(indod,2) * rerth +cc + enddo +cc + enddo +cc + return + end + + end module dezouv_stochy_mod diff --git a/stochastic_physics/dozeuv_stochy.f b/stochastic_physics/dozeuv_stochy.f new file mode 100644 index 000000000..4ff5ad8f2 --- /dev/null +++ b/stochastic_physics/dozeuv_stochy.f @@ -0,0 +1,267 @@ + module dozeuv_stochy_mod + + implicit none + + contains + + subroutine dozeuv_stochy(dod,zev,uod,vev,epsedn,epsodn, + & snnp1ev,snnp1od,ls_node) +cc + use stochy_resol_def + use spectral_layout_mod + use machine + implicit none +cc + real(kind_dbl_prec) dod(len_trio_ls,2) + real(kind_dbl_prec) zev(len_trie_ls,2) + real(kind_dbl_prec) uod(len_trio_ls,2) + real(kind_dbl_prec) vev(len_trie_ls,2) +cc + real(kind_dbl_prec) epsedn(len_trie_ls) + real(kind_dbl_prec) epsodn(len_trio_ls) +cc + real(kind_dbl_prec) snnp1ev(len_trie_ls) + real(kind_dbl_prec) snnp1od(len_trio_ls) +cc + integer ls_node(ls_dim,3) +cc +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +cc + integer l,locl,n +cc + integer indev,indev1,indev2 + integer indod,indod1,indod2 + integer inddif +cc + real(kind_dbl_prec) rl +cc + real(kind_dbl_prec) cons0 !constant +cc + integer indlsev,jbasev + integer indlsod,jbasod + real(kind_evod) rerth +cc + include 'function2' +cc +cc +cc...................................................................... +cc +cc + cons0 = 0.d0 !constant + rerth =6.3712e+6 ! radius of earth (m) +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) +cc + vev(indlsev(l,l),1) = cons0 !constant + vev(indlsev(l,l),2) = cons0 !constant +cc +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + uod(indev-inddif,1) = -epsodn(indev-inddif) + x * zev(indev,1) +cc + uod(indev-inddif,2) = -epsodn(indev-inddif) + x * zev(indev,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + 1 + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + vev(indev,1) = epsedn(indev) + x * dod(indev-inddif,1) +cc + vev(indev,2) = epsedn(indev) + x * dod(indev-inddif,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasod=ls_node(locl,3) + indod1 = indlsod(L+1,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indod2 = indlsod(jcap ,L) + else + indod2 = indlsod(jcap+1,L) - 1 + endif + if ( l .ge. 1 ) then + rl = l + do indod = indod1 , indod2 +cc u(l,n)=-i*l*d(l,n)/(n*(n+1)) +cc + uod(indod,1) = uod(indod,1) + 1 + rl * dod(indod,2) + 2 / snnp1od(indod) +cc + uod(indod,2) = uod(indod,2) + 1 - rl * dod(indod,1) + 2 / snnp1od(indod) +cc + enddo + endif +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + if ( l .ge. 1 ) then + rl = l + do indev = indev1 , indev2 +cc u(l,n)=-i*l*d(l,n)/(n*(n+1)) +cc + vev(indev,1) = vev(indev,1) + 1 + rl * zev(indev,2) + 2 / snnp1ev(indev) +cc + vev(indev,2) = vev(indev,2) + 1 - rl * zev(indev,1) + 2 / snnp1ev(indev) +cc + enddo + endif +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + 1 + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + uod(indev-inddif,1) = uod(indev-inddif,1) + 1 + epsedn(indev) * zev(indev ,1) +cc + uod(indev-inddif,2) = uod(indev-inddif,2) + 1 + epsedn(indev) * zev(indev ,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) - 1 + else + indev2 = indlsev(jcap ,L) - 1 + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + do indev = indev1 , indev2 +cc + vev(indev,1) = vev(indev ,1) + 1 - epsodn(indev-inddif) * dod(indev-inddif,1) +cc + vev(indev,2) = vev(indev ,2) + 1 - epsodn(indev-inddif) * dod(indev-inddif,2) +cc + enddo +cc + enddo +cc +cc...................................................................... +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + indod1 = indlsod(L+1,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) + indod2 = indlsod(jcap ,L) + else + indev2 = indlsev(jcap ,L) + indod2 = indlsod(jcap+1,L) + endif + do indod = indod1 , indod2 +cc + uod(indod,1) = uod(indod,1) * rerth + uod(indod,2) = uod(indod,2) * rerth +cc + enddo +cc + do indev = indev1 , indev2 +cc + vev(indev,1) = vev(indev,1) * rerth + vev(indev,2) = vev(indev,2) * rerth +cc + enddo +cc + enddo +cc + return + end + + end module dozeuv_stochy_mod diff --git a/stochastic_physics/epslon_stochy.f b/stochastic_physics/epslon_stochy.f new file mode 100644 index 000000000..c7aace515 --- /dev/null +++ b/stochastic_physics/epslon_stochy.f @@ -0,0 +1,93 @@ + module epslon_stochy_mod + + implicit none + + contains + + subroutine epslon_stochy(epse,epso,epsedn,epsodn, + & ls_node) +cc + use stochy_resol_def + use spectral_layout_mod + use machine + implicit none +cc + real(kind_dbl_prec) epse(len_trie_ls) + real(kind_dbl_prec) epso(len_trio_ls) +cc + real(kind_dbl_prec) epsedn(len_trie_ls) + real(kind_dbl_prec) epsodn(len_trio_ls) +cc + integer ls_node(ls_dim,3) +cc +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +cc + integer l,locl,n +cc + integer indev + integer indod +cc + real(kind_dbl_prec) f1,f2,rn,val +cc + real(kind_dbl_prec) cons0 !constant +cc + integer indlsev,jbasev + integer indlsod,jbasod +cc + include 'function2' +cc +cc + cons0=0.0d0 !constant +cc +cc +cc...................................................................... +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + indev=indlsev(l,l) + epse (indev)=cons0 !constant + epsedn(indev)=cons0 !constant + indev=indev+1 +cc + + do n=l+2,jcap+1,2 + rn=n + f1=n*n-l*l + f2=4*n*n-1 + val=sqrt(f1/f2) + epse (indev)=val + epsedn(indev)=val/rn + indev=indev+1 + enddo +cc + enddo +cc +cc +cc...................................................................... +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasod=ls_node(locl,3) + indod=indlsod(l+1,l) +cc + do n=l+1,jcap+1,2 + rn=n + f1=n*n-l*l + f2=4*n*n-1 + val=sqrt(f1/f2) + epso (indod)=val + epsodn(indod)=val/rn + indod=indod+1 + enddo +cc + enddo +cc + return + end + + end module epslon_stochy_mod diff --git a/stochastic_physics/four_to_grid_stochy.F b/stochastic_physics/four_to_grid_stochy.F new file mode 100644 index 000000000..5f26a0a7e --- /dev/null +++ b/stochastic_physics/four_to_grid_stochy.F @@ -0,0 +1,271 @@ + module four_to_grid_mod + + use stochy_ccpp, only: num_parthds_stochy => ompthreads + + implicit none + + contains + + subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, + & lon_dim_coef,lon_dim_grid,lons_lat,lot) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + use machine + implicit none +!! + real(kind=kind_dbl_prec) syn_gr_a_1(lon_dim_coef,lot) + real(kind=kind_dbl_prec) syn_gr_a_2(lon_dim_grid,lot) + integer lon_dim_coef + integer lon_dim_grid + integer lons_lat + integer lot +!________________________________________________________ +#ifdef MKL + integer*8 plan +#else + real(kind=kind_dbl_prec) aux1crs(42002) + real(kind=kind_dbl_prec) scale_ibm + integer ibmsign + integer init +#endif + integer lot_thread + integer num_threads + integer nvar_thread_max + integer nvar_1 + integer nvar_2 + integer thread +#ifdef MKL + include "fftw/fftw3.f" + integer NULL +#else + external dcrft + external scrft +#endif +!________________________________________________________ + num_threads = min(num_parthds_stochy,lot) + + nvar_thread_max = (lot+num_threads-1)/num_threads + + if ( kind_dbl_prec == 8 ) then !------------------------------------ +#ifdef MKL +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,plan) +#else +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+shared(ibmsign,scale_ibm) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) +#endif + do thread=1,num_threads ! start of thread loop .............. + nvar_1=(thread-1)*nvar_thread_max + 1 + nvar_2=min(nvar_1+nvar_thread_max-1,lot) + + lot_thread=nvar_2 - nvar_1 + 1 + + if (nvar_2 >= nvar_1) then +#ifdef MKL + !call dfftw_plan_many_dft_c2r( + ! plan, 1, N, m, & + ! X, NULL, 1, dimx, & + ! Y, NULL, 1, dimy, & + ! fftw_flag) + call dfftw_plan_many_dft_c2r( & + & plan, 1, lons_lat, lot_thread, & + & syn_gr_a_1, NULL, 1, size(syn_gr_a_1,dim=1), & + & syn_gr_a_2, NULL, 1, size(syn_gr_a_2,dim=1), & + & FFTW_ESTIMATE) + call dfftw_execute(plan) + call dfftw_destroy_plan(plan) +#else + init = 1 + ibmsign = -1 + scale_ibm = 1.0d0 + + call dcrft(init, + & syn_gr_a_1(1,nvar_1) ,lon_dim_coef/2, + & syn_gr_a_2(1,nvar_1) ,lon_dim_grid, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000) + + init = 0 + call dcrft(init, + & syn_gr_a_1(1,nvar_1) ,lon_dim_coef/2, + & syn_gr_a_2(1,nvar_1) ,lon_dim_grid, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000) +#endif + endif + + enddo ! fin thread loop ...................................... + else !------------------------------------------------------------ +#ifdef MKL +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,plan) +#else +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+shared(ibmsign,scale_ibm) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) +#endif + do thread=1,num_threads ! start of thread loop .............. + nvar_1 = (thread-1)*nvar_thread_max + 1 + nvar_2 = min(nvar_1+nvar_thread_max-1,lot) + + lot_thread = nvar_2 - nvar_1 + 1 + + if (nvar_2 >= nvar_1) then +#ifdef MKL + !call sfftw_plan_many_dft_c2r( + ! plan, 1, N, m, & + ! X, NULL, 1, dimx, & + ! Y, NULL, 1, dimy, & + ! fftw_flag) + call sfftw_plan_many_dft_c2r( & + & plan, 1, lons_lat, lot_thread, & + & syn_gr_a_1, NULL, 1, size(syn_gr_a_1,dim=1), & + & syn_gr_a_2, NULL, 1, size(syn_gr_a_2,dim=1), & + & FFTW_ESTIMATE) + call sfftw_execute(plan) + call sfftw_destroy_plan(plan) +#else + init = 1 + ibmsign = -1 + scale_ibm = 1.0d0 + call scrft(init, + & syn_gr_a_1(1,nvar_1) ,lon_dim_coef/2, + & syn_gr_a_2(1,nvar_1) ,lon_dim_grid, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000, + & aux1crs(22001),0) + init = 0 + call scrft(init, + & syn_gr_a_1(1,nvar_1) ,lon_dim_coef/2, + & syn_gr_a_2(1,nvar_1) ,lon_dim_grid, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000, + & aux1crs(22001),0) +#endif + endif + enddo ! fin thread loop ...................................... + endif !----------------------------------------------------------- +!! + return + end + subroutine grid_to_four(anl_gr_a_2,anl_gr_a_1, + & lon_dim_grid,lon_dim_coef,lons_lat,lot) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + use machine + implicit none +!! + real(kind=kind_dbl_prec) anl_gr_a_2(lon_dim_grid,lot) + real(kind=kind_dbl_prec) anl_gr_a_1(lon_dim_coef,lot) + integer lon_dim_grid + integer lon_dim_coef + integer lons_lat + integer lot +!________________________________________________________ + real(kind=kind_dbl_prec) aux1crs(42002) + real(kind=kind_dbl_prec) scale_ibm,rone + integer ibmsign + integer init + integer lot_thread + integer num_threads + integer nvar_thread_max + integer nvar_1,nvar_2 + integer thread +!________________________________________________________ +#ifdef MKL + write(0,*) "ERROR in grid_to_four: srcft and drcft ", + & " must be replaced with MKL's FFTW calls. ABORT." + call sleep(5) + stop +#endif + num_threads=min(num_parthds_stochy,lot) + + nvar_thread_max=(lot+num_threads-1)/num_threads + + if ( kind_dbl_prec == 8 ) then !------------------------------------ +!$omp parallel do shared(anl_gr_a_1,anl_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+shared(ibmsign,scale_ibm,rone) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) + + do thread=1,num_threads ! start of thread loop .............. + nvar_1 = (thread-1)*nvar_thread_max + 1 + nvar_2 = min(nvar_1+nvar_thread_max-1,lot) + + if (nvar_2 >= nvar_1) then + lot_thread = nvar_2 - nvar_1 + 1 + + init = 1 + ibmsign = 1 + rone = 1.0d0 + scale_ibm = rone/lons_lat + call drcft(init, + & anl_gr_a_2(1,nvar_1), lon_dim_grid, + & anl_gr_a_1(1,nvar_1), lon_dim_coef/2, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000) + init = 0 + call drcft(init, + & anl_gr_a_2(1,nvar_1), lon_dim_grid, + & anl_gr_a_1(1,nvar_1), lon_dim_coef/2, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000) + + endif + enddo ! fin thread loop ...................................... + else !------------------------------------------------------------ +!$omp parallel do shared(anl_gr_a_1,anl_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+shared(ibmsign,scale_ibm,rone) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) + + do thread=1,num_threads ! start of thread loop .............. + nvar_1 = (thread-1)*nvar_thread_max + 1 + nvar_2 = min(nvar_1+nvar_thread_max-1,lot) + + if (nvar_2 >= nvar_1) then + lot_thread=nvar_2 - nvar_1 + 1 + + init = 1 + ibmsign = 1 + rone = 1.0d0 + scale_ibm = rone/lons_lat + call srcft(init, + & anl_gr_a_2(1,nvar_1), lon_dim_grid, + & anl_gr_a_1(1,nvar_1), lon_dim_coef/2, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000, + & aux1crs(22001),0) + init = 0 + call srcft(init, + & anl_gr_a_2(1,nvar_1), lon_dim_grid, + & anl_gr_a_1(1,nvar_1), lon_dim_coef/2, + & lons_lat,lot_thread,ibmsign,scale_ibm, + & aux1crs,22000, + & aux1crs(22001),20000, + & aux1crs(22001),0) + + endif + enddo ! fin thread loop ...................................... + endif !----------------------------------------------------------- +!! + return + end + + end module four_to_grid_mod diff --git a/stochastic_physics/function2 b/stochastic_physics/function2 new file mode 100644 index 000000000..f3328235b --- /dev/null +++ b/stochastic_physics/function2 @@ -0,0 +1,5 @@ +!cc + indlsev(n,l) = jbasev + (n-l)/2 + 1 +!cc + indlsod(n,l) = jbasod + (n-l)/2 + 1 +!cc diff --git a/stochastic_physics/function_indlsev b/stochastic_physics/function_indlsev new file mode 100644 index 000000000..94e605cfb --- /dev/null +++ b/stochastic_physics/function_indlsev @@ -0,0 +1,3 @@ +!cc + indlsev(n,l) = jbasev + (n-l)/2 + 1 +!cc diff --git a/stochastic_physics/function_indlsod b/stochastic_physics/function_indlsod new file mode 100644 index 000000000..16c4e7996 --- /dev/null +++ b/stochastic_physics/function_indlsod @@ -0,0 +1,3 @@ +!cc + indlsod(n,l) = jbasod + (n-l)/2 + 1 +!cc diff --git a/stochastic_physics/get_lats_node_a_stochy.f b/stochastic_physics/get_lats_node_a_stochy.f new file mode 100644 index 000000000..62bc7b293 --- /dev/null +++ b/stochastic_physics/get_lats_node_a_stochy.f @@ -0,0 +1,92 @@ + module get_lats_node_a_stochy_mod + + implicit none + + contains + + subroutine get_lats_node_a_stochy(me_fake,global_lats_a, + & lats_nodes_a_fake,gl_lats_index, + & global_time_sort_index,iprint) +cc + use stochy_resol_def + use spectral_layout_mod + implicit none +cc + integer gl_lats_index,gl_start + integer me_fake + integer global_lats_a(latg) + integer lats_nodes_a_fake + integer iprint +cc + integer ijk + integer jptlats + integer lat + integer node,nodesio + integer global_time_sort_index(latg) + integer nodes_tmp +cc +c +!jw if (liope) then +!jw if (icolor.eq.2) then +!jw nodesio=1 +!jw else + nodesio=nodes +!jw endif +!jw else +!jw nodesio=nodes +!jw endif +!! +cc + lat = 1 + nodes_tmp = nodes +!jw if (liope .and. icolor .eq. 2) nodes_tmp = nodes -1 + + gl_start = gl_lats_index +cc............................................. + do ijk=1,latg +cc + do node=1,nodes_tmp + if (node.eq.me_fake+1) then + gl_lats_index=gl_lats_index+1 + global_lats_a(gl_lats_index) = global_time_sort_index(lat) + endif + lat = lat + 1 + if (lat .gt. latg) go to 200 + enddo +cc + do node=nodes_tmp,1,-1 + if (node.eq.me_fake+1) then + gl_lats_index=gl_lats_index+1 + global_lats_a(gl_lats_index) = global_time_sort_index(lat) + endif + lat = lat + 1 + if (lat .gt. latg) go to 200 + enddo +cc + enddo +cc............................................. +cc + 200 continue +cc +cc............................................. +cc +!jw if (liope .and. icolor .eq. 2) gl_start = 0 + do node=1,nodes_tmp + if (node.eq.me_fake+1) then + lats_nodes_a_fake=gl_lats_index-gl_start +c$$$ print*,' setting lats_nodes_a_fake = ', +c$$$ . lats_nodes_a_fake + endif + enddo + + if(iprint.eq.1) print 220 + 220 format ('completed loop 200 in get_lats_a ') +c + if(iprint.eq.1) + & print*,'completed get_lats_node, lats_nodes_a_fake=', + & lats_nodes_a_fake +cc + return + end + + end module get_lats_node_a_stochy_mod diff --git a/stochastic_physics/get_ls_node_stochy.f b/stochastic_physics/get_ls_node_stochy.f new file mode 100644 index 000000000..51d9f85c3 --- /dev/null +++ b/stochastic_physics/get_ls_node_stochy.f @@ -0,0 +1,81 @@ + module get_ls_node_stochy_mod + + implicit none + + contains + + subroutine get_ls_node_stochy(me_fake,ls_node,ls_max_node_fake, + c iprint) +! + use stochy_resol_def + use spectral_layout_mod + implicit none +! + integer me_fake, ls_max_node_fake, iprint + integer ls_node(ls_dim) + + integer ijk, jptls, l, node, nodesio +! +!jw if (liope) then +!jw if (icolor.eq.2) then +!jw nodesio=1 +!jw else + + nodesio = nodes + +!jw endif +!jw else +!jw nodesio=nodes +!jw endif +!! + ls_node = -1 +! + jptls = 0 + l = 0 +!............................................. + do ijk=1,jcap1 +! + do node=1,nodesio + if (node == me_fake+1) then + jptls = jptls + 1 + ls_node(jptls) = l + endif + l = l + 1 + if (l > jcap) go to 200 + enddo +! + do node=nodesio,1,-1 + if (node == me_fake+1) then + jptls = jptls + 1 + ls_node(jptls) = l + endif + l = l + 1 + if (l > jcap) go to 200 + enddo +! + enddo +!............................................. +! + 200 continue +! +!............................................. +! + if(iprint == 1) print *, 'completed loop 200 in get_ls_node' + ls_max_node_fake = 0 + do ijk=1,ls_dim + if(ls_node(ijk) >= 0) then + ls_max_node_fake = ijk + if(iprint == 1) + & print 230, me_fake, ijk, ls_node(ijk) + endif + 230 format ('me_fake=',i5,' get_ls_node ls_node(', i5, ')=',i5) + enddo +! + if(iprint == 1) + & print*,'completed get_ls_node, ls_max_node_fake=', + & ls_max_node_fake +! + return + end + + end module get_ls_node_stochy_mod diff --git a/stochastic_physics/get_stochy_pattern.F90 b/stochastic_physics/get_stochy_pattern.F90 new file mode 100644 index 000000000..acbd70e32 --- /dev/null +++ b/stochastic_physics/get_stochy_pattern.F90 @@ -0,0 +1,518 @@ +module get_stochy_pattern_mod + use machine, only : kind_dbl_prec, kind_evod + use stochy_ccpp, only : nodes => mpisize, stochy_la2ga + use stochy_resol_def, only : latg, latg2, levs, lonf, skeblevs + use spectral_layout_mod, only : ipt_lats_node_a, lat1s_a, lats_dim_a, & + lats_node_a, lon_dim_a, len_trie_ls, & + len_trio_ls, ls_dim + use stochy_namelist_def, only : nsfcpert, ntrunc, stochini + use stochy_data_mod, only : gg_lats, gg_lons, inttyp, nskeb, nshum, nsppt, & + rad2deg, rnlat, rpattern_sfc, rpattern_skeb, & + rpattern_shum, rpattern_sppt, skebu_save, & + skebv_save, skeb_vwts, skeb_vpts, wlon + use stochy_gg_def, only : coslat_a + use stochy_patterngenerator_mod, only: random_pattern, ndimspec, & + patterngenerator_advance + use stochy_internal_state_mod, only: stochy_internal_state + use stochy_ccpp, only : is_master, mp_reduce_sum, mpicomm + use GFS_typedefs, only: GFS_control_type, GFS_grid_type + use mersenne_twister, only: random_seed + use dezouv_stochy_mod, only: dezouv_stochy + use dozeuv_stochy_mod, only: dozeuv_stochy + use four_to_grid_mod, only: four_to_grid + use sumfln_stochy_mod, only: sumfln_stochy + implicit none + private + + public get_random_pattern_fv3,get_random_pattern_fv3_vect + public get_random_pattern_sfc_fv3 + public dump_patterns + logical :: first_call=.true. + contains + +subroutine get_random_pattern_fv3(rpattern,npatterns,& + gis_stochy,Model,Grid,nblks,maxlen,pattern_2d) + +! generate a random pattern for stochastic physics + implicit none + type(random_pattern), intent(inout) :: rpattern(npatterns) + type(stochy_internal_state) :: gis_stochy + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid(nblks) + integer,intent(in):: npatterns,nblks,maxlen + + integer i,j,l,lat,ierr,n,nn,k,nt + real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d + + integer :: num2d +! logical lprint + + real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg + real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a) + integer kmsk0(lonf,gis_stochy%lats_node_a),len + real(kind=kind_dbl_prec) :: globalvar,globalvar0 + real(kind=kind_dbl_prec) :: pattern_2d(nblks,maxlen) + real(kind=kind_dbl_prec) :: pattern_1d(maxlen) + real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk + integer :: blk + + kmsk0 = 0 + glolal = 0. + do n=1,npatterns + call patterngenerator_advance(rpattern(n),1,.false.) + call scalarspect_to_gaugrid( & + rpattern(n)%spec_e,rpattern(n)%spec_o,wrk2d,& + gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& + gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,gis_stochy%lonsperlat,& + gis_stochy%plnev_a,gis_stochy%plnod_a,1) + glolal = glolal + wrk2d(:,:,1) + enddo + + allocate(workg(lonf,latg)) + workg = 0. + do j=1,gis_stochy%lats_node_a + lat=gis_stochy%global_lats_a(ipt_lats_node_a-1+j) + do i=1,lonf + workg(i,lat) = glolal(i,j) + enddo + enddo + + call mp_reduce_sum(workg,lonf,latg) + +! interpolate to cube grid + + allocate(rslmsk(lonf,latg)) + do blk=1,nblks + len=size(Grid(blk)%xlat,1) + pattern_1d = 0 + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + pattern_2d(blk,:)=pattern_1d(:) + end associate + enddo + deallocate(rslmsk) + deallocate(workg) + +end subroutine get_random_pattern_fv3 + + +subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& + gis_stochy,Model,Grid,nblks,maxlen,pattern_3d) + +! generate a random pattern for stochastic physics + implicit none + type(random_pattern), intent(inout) :: rpattern(npatterns) + type(stochy_internal_state), target :: gis_stochy + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid(nblks) + integer,intent(in):: npatterns,nblks,maxlen + + integer i,j,l,lat,ierr,n,nn,k,nt + real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d + + integer :: num2d +! logical lprint + + real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg + real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a) + integer kmsk0(lonf,gis_stochy%lats_node_a),len + real(kind=kind_dbl_prec) :: globalvar,globalvar0 + real(kind=kind_dbl_prec) :: pattern_3d(nblks,maxlen,nsfcpert) + real(kind=kind_dbl_prec) :: pattern_1d(maxlen) + real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk + integer :: blk + + do k=1,nsfcpert + kmsk0 = 0 + glolal = 0. + do n=1,npatterns + if (is_master()) print *, 'Random pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k)) + call scalarspect_to_gaugrid( & + rpattern(n)%spec_e(:,:,k),rpattern(n)%spec_o(:,:,k),wrk2d,& + gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& + gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,gis_stochy%lonsperlat,& + gis_stochy%plnev_a,gis_stochy%plnod_a,1) + glolal = glolal + wrk2d(:,:,1) + enddo + + allocate(workg(lonf,latg)) + workg = 0. + do j=1,gis_stochy%lats_node_a + lat=gis_stochy%global_lats_a(ipt_lats_node_a-1+j) + do i=1,lonf + workg(i,lat) = glolal(i,j) + enddo + enddo + + call mp_reduce_sum(workg,lonf,latg) + if (is_master()) print *, 'workg after mp_reduce_sum for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(workg), maxval(workg) + +! interpolate to cube grid + + allocate(rslmsk(lonf,latg)) + do blk=1,nblks + len=size(Grid(blk)%xlat,1) + pattern_1d = 0 + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + pattern_3d(blk,:,k)=pattern_1d(:) + end associate + enddo + if (is_master()) print *, '3D pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k)) + deallocate(rslmsk) + deallocate(workg) + + enddo ! loop over k, nsfcpert + +end subroutine get_random_pattern_sfc_fv3 + + +subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& + gis_stochy,Model,Grid,nblks,maxlen,upattern_3d,vpattern_3d) + +! generate a random pattern for stochastic physics + implicit none + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid(nblks) + type(stochy_internal_state), target :: gis_stochy + type(random_pattern), intent(inout) :: rpattern(npatterns) + + real(kind=kind_evod), dimension(len_trie_ls,2,1) :: vrtspec_e,divspec_e + real(kind=kind_evod), dimension(len_trio_ls,2,1) :: vrtspec_o,divspec_o + integer:: npatterns,nblks,blk,len,maxlen + + real(kind=kind_dbl_prec) :: upattern_3d(nblks,maxlen,levs) + real(kind=kind_dbl_prec) :: vpattern_3d(nblks,maxlen,levs) + real(kind=kind_dbl_prec) :: pattern_1d(maxlen) + real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk + integer i,j,l,lat,ierr,n,nn,k,nt + real(kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2du,wrk2dv + + integer :: num2d +! logical lprint + + real, allocatable, dimension(:,:) :: workgu,workgv + integer kmsk0(lonf,gis_stochy%lats_node_a),i1,i2,j1 + real(kind=kind_dbl_prec) :: globalvar,globalvar0 + kmsk0 = 0 + allocate(workgu(lonf,latg)) + allocate(workgv(lonf,latg)) + allocate(rslmsk(lonf,latg)) + if (first_call) then + allocate(skebu_save(nblks,maxlen,skeblevs)) + allocate(skebv_save(nblks,maxlen,skeblevs)) + do k=2,skeblevs + workgu = 0. + workgv = 0. + do n=1,npatterns + if (.not. stochini) call patterngenerator_advance(rpattern(n),k,first_call) + ! ke norm (convert streamfunction forcing to vorticity forcing) + divspec_e = 0; divspec_o = 0. + do nn=1,2 + vrtspec_e(:,nn,1) = gis_stochy%kenorm_e*rpattern(n)%spec_e(:,nn,k) + vrtspec_o(:,nn,1) = gis_stochy%kenorm_o*rpattern(n)%spec_o(:,nn,k) + enddo + ! convert to winds + call vrtdivspect_to_uvgrid(& + divspec_e,divspec_o,vrtspec_e,vrtspec_o,& + wrk2du,wrk2dv,& + gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& + gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,gis_stochy%lonsperlat,& + gis_stochy%epsedn,gis_stochy%epsodn,gis_stochy%snnp1ev,gis_stochy%snnp1od,& + gis_stochy%plnev_a,gis_stochy%plnod_a,1) + do i=1,lonf + do j=1,gis_stochy%lats_node_a + lat=gis_stochy%global_lats_a(ipt_lats_node_a-1+j) + workgu(i,lat) = workgu(i,lat) + wrk2du(i,j,1) + workgv(i,lat) = workgv(i,lat) + wrk2dv(i,j,1) + enddo + enddo + enddo + call mp_reduce_sum(workgu,lonf,latg) + call mp_reduce_sum(workgv,lonf,latg) +! interpolate to cube grid + do blk=1,nblks + len=size(Grid(blk)%xlat,1) + pattern_1d = 0 + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + skebu_save(blk,:,k)=pattern_1d(:) + call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + skebv_save(blk,:,k)=-1*pattern_1d(:) + end associate + enddo + enddo + endif + do k=1,skeblevs-1 + skebu_save(:,:,k)=skebu_save(:,:,k+1) + skebv_save(:,:,k)=skebv_save(:,:,k+1) + do n=1,npatterns + rpattern(n)%spec_e(:,:,k)=rpattern(n)%spec_e(:,:,k+1) + rpattern(n)%spec_o(:,:,k)=rpattern(n)%spec_o(:,:,k+1) + enddo + enddo + +! get pattern for last level + workgu = 0. + workgv = 0. + do n=1,npatterns +! if (stochini.AND. first_call) then +! print*,'skipping advance' +! else + call patterngenerator_advance(rpattern(n),skeblevs,first_call) +! endif +! ke norm (convert streamfunction forcing to vorticity forcing) + divspec_e = 0; divspec_o = 0. + do nn=1,2 + vrtspec_e(:,nn,1) = gis_stochy%kenorm_e*rpattern(n)%spec_e(:,nn,skeblevs) + vrtspec_o(:,nn,1) = gis_stochy%kenorm_o*rpattern(n)%spec_o(:,nn,skeblevs) + enddo + ! convert to winds + call vrtdivspect_to_uvgrid(& + divspec_e,divspec_o,vrtspec_e,vrtspec_o,& + wrk2du,wrk2dv,& + gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& + gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,gis_stochy%lonsperlat,& + gis_stochy%epsedn,gis_stochy%epsodn,gis_stochy%snnp1ev,gis_stochy%snnp1od,& + gis_stochy%plnev_a,gis_stochy%plnod_a,1) + do i=1,lonf + do j=1,gis_stochy%lats_node_a + lat=gis_stochy%global_lats_a(ipt_lats_node_a-1+j) + workgu(i,lat) = workgu(i,lat) + wrk2du(i,j,1) + workgv(i,lat) = workgv(i,lat) + wrk2dv(i,j,1) + enddo + enddo + enddo + call mp_reduce_sum(workgu,lonf,latg) + call mp_reduce_sum(workgv,lonf,latg) +! interpolate to cube grid + do blk=1,nblks + len=size(Grid(blk)%xlat,1) + pattern_1d = 0 + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + skebu_save(blk,:,skeblevs)=pattern_1d(:) + call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) + skebv_save(blk,:,skeblevs)=-1*pattern_1d(:) + end associate + enddo + deallocate(rslmsk) + deallocate(workgu) + deallocate(workgv) +! interpolate in the vertical ! consider moving to cubed sphere side, more memory, but less interpolations + do k=1,Model%levs + do blk=1,nblks + upattern_3d(blk,:,k) = skeb_vwts(k,1)*skebu_save(blk,:,skeb_vpts(k,1))+skeb_vwts(k,2)*skebu_save(blk,:,skeb_vpts(k,2)) + vpattern_3d(blk,:,k) = skeb_vwts(k,1)*skebv_save(blk,:,skeb_vpts(k,1))+skeb_vwts(k,2)*skebv_save(blk,:,skeb_vpts(k,2)) + enddo + enddo + first_call=.false. + +end subroutine get_random_pattern_fv3_vect + +subroutine scalarspect_to_gaugrid(& + trie_ls,trio_ls,datag,& + ls_node,ls_nodes,max_ls_nodes,& + lats_nodes_a,global_lats_a,lonsperlat,& + plnev_a,plnod_a,nlevs) + + + implicit none + real(kind=kind_dbl_prec), intent(in) :: trie_ls(len_trie_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(in) :: trio_ls(len_trio_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(out) :: datag(lonf,lats_node_a,nlevs) + integer, intent(in) :: ls_node(ls_dim,3),ls_nodes(ls_dim,nodes),& + nlevs,max_ls_nodes(nodes),lats_nodes_a(nodes),global_lats_a(latg),lonsperlat(latg) + real(kind=kind_dbl_prec),intent(in) :: plnev_a(len_trie_ls,latg2),plnod_a(len_trio_ls,latg2) +! local vars + real(kind=kind_dbl_prec) for_gr_a_1(lon_dim_a,nlevs,lats_dim_a) + real(kind=kind_dbl_prec) for_gr_a_2(lonf,nlevs,lats_dim_a) + integer i,j,k + integer l,lan,lat + integer lons_lat + + call sumfln_stochy(trie_ls,& + trio_ls,& + lat1s_a,& + plnev_a,plnod_a,& + nlevs,ls_node,latg2,& + lats_dim_a,nlevs,for_gr_a_1,& + ls_nodes,max_ls_nodes,& + lats_nodes_a,global_lats_a,& + lats_node_a,ipt_lats_node_a,& + lonsperlat,lon_dim_a,latg,0) + + do lan=1,lats_node_a + lat = global_lats_a(ipt_lats_node_a-1+lan) + lons_lat = lonsperlat(lat) + CALL FOUR_TO_GRID(for_gr_a_1(1,1,lan),for_gr_a_2(1,1,lan),& + lon_dim_a,lonf,lons_lat,nlevs) + enddo + + datag = 0. + do lan=1,lats_node_a + lat = global_lats_a(ipt_lats_node_a-1+lan) + lons_lat = lonsperlat(lat) + do k=1,nlevs + do i=1,lons_lat + datag(i,lan,k) = for_gr_a_2(i,k,lan) + enddo + enddo + enddo + + return + end subroutine scalarspect_to_gaugrid + +subroutine dump_patterns(sfile) + implicit none + character*120 :: sfile + integer :: stochlun,k,n + stochlun=99 + if (is_master()) then + if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0) then + OPEN(stochlun,file=sfile,form='unformatted') + print*,'open ',sfile,' for output' + endif + endif + if (nsppt > 0) then + do n=1,nsppt + call write_pattern(rpattern_sppt(n),1,stochlun) + enddo + endif + if (nshum > 0) then + do n=1,nshum + call write_pattern(rpattern_shum(n),1,stochlun) + enddo + endif + if (nskeb > 0) then + do n=1,nskeb + do k=1,skeblevs + call write_pattern(rpattern_skeb(n),k,stochlun) + enddo + enddo + endif + close(stochlun) + end subroutine dump_patterns + subroutine write_pattern(rpattern,lev,lunptn) + implicit none + type(random_pattern), intent(inout) :: rpattern + integer, intent(in) :: lunptn,lev + real(kind_dbl_prec), allocatable :: pattern2d(:) + integer nm,nn,ierr,arrlen,isize + integer,allocatable :: isave(:) + arrlen=2*ndimspec + + allocate(pattern2d(arrlen)) + pattern2d=0.0 + ! fill in apprpriate pieces of array + !print*,'before collection...',me,maxval(rpattern%spec_e),maxval(rpattern%spec_o) & + ! ,minval(rpattern%spec_e),minval(rpattern%spec_o) + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + pattern2d(nm) = rpattern%spec_e(nn,1,lev) + pattern2d(ndimspec+nm) = rpattern%spec_e(nn,2,lev) + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + pattern2d(nm) = rpattern%spec_o(nn,1,lev) + pattern2d(ndimspec+nm) = rpattern%spec_o(nn,2,lev) + enddo + call mp_reduce_sum(pattern2d,arrlen) + ! write only on root process + if (is_master()) then + print*,'writing out random pattern (min/max/size)',& + minval(pattern2d),maxval(pattern2d),size(pattern2d) + !print*,'max/min pattern=',maxval(pattern2d),minval(pattern2d) + write(lunptn) ntrunc + call random_seed(size=isize) ! get seed size + allocate(isave(isize)) ! get seed + call random_seed(get=isave,stat=rpattern%rstate) ! write seed + write(lunptn) isave + write(lunptn) pattern2d + endif + deallocate(pattern2d) + end subroutine write_pattern + subroutine vrtdivspect_to_uvgrid(& + trie_di,trio_di,trie_ze,trio_ze,& + uug,vvg,& + ls_node,ls_nodes,max_ls_nodes,& + lats_nodes_a,global_lats_a,lonsperlar,& + epsedn,epsodn,snnp1ev,snnp1od,plnev_a,plnod_a,nlevs) + + implicit none + real(kind=kind_dbl_prec), intent(in) :: trie_di(len_trie_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(in) :: trio_di(len_trio_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(in) :: trie_ze(len_trie_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(in) :: trio_ze(len_trio_ls,2,nlevs) + real(kind=kind_dbl_prec), intent(out) :: uug(lonf,lats_node_a,nlevs) + real(kind=kind_dbl_prec), intent(out) :: vvg(lonf,lats_node_a,nlevs) + integer, intent(in) :: ls_node(ls_dim,3),ls_nodes(ls_dim,nodes),& + nlevs,max_ls_nodes(nodes),lats_nodes_a(nodes),global_lats_a(latg),lonsperlar(latg) + real(kind=kind_dbl_prec),intent(in) :: epsedn(len_trie_ls),& + epsodn(len_trio_ls),snnp1ev(len_trie_ls),snnp1od(len_trio_ls),& + plnev_a(len_trie_ls,latg2),plnod_a(len_trio_ls,latg2) +! local vars + real(kind=kind_dbl_prec) trie_ls(len_trie_ls,2,2*nlevs) + real(kind=kind_dbl_prec) trio_ls(len_trio_ls,2,2*nlevs) + real(kind=kind_dbl_prec) for_gr_a_1(lon_dim_a,2*nlevs,lats_dim_a) + real(kind=kind_dbl_prec) for_gr_a_2(lonf,2*nlevs,lats_dim_a) + integer i,j,k + integer l,lan,lat + integer lons_lat + real (kind=kind_dbl_prec) tx1 + + do k=1,nlevs + call dezouv_stochy(trie_di(1,1,k), trio_ze(1,1,k),& + trie_ls(1,1,k), trio_ls(1,1,nlevs+k),& + epsedn,epsodn,snnp1ev,snnp1od,ls_node) + call dozeuv_stochy(trio_di(1,1,k), trie_ze(1,1,k),& + trio_ls(1,1,k), trie_ls(1,1,nlevs+k),& + epsedn,epsodn,snnp1ev,snnp1od,ls_node) + enddo + + call sumfln_stochy(trie_ls,& + trio_ls,& + lat1s_a,& + plnev_a,plnod_a,& + 2*nlevs,ls_node,latg2,& + lats_dim_a,2*nlevs,for_gr_a_1,& + ls_nodes,max_ls_nodes,& + lats_nodes_a,global_lats_a,& + lats_node_a,ipt_lats_node_a,& + lonsperlar,lon_dim_a,latg,0) + + do lan=1,lats_node_a + lat = global_lats_a(ipt_lats_node_a-1+lan) + lons_lat = lonsperlar(lat) + CALL FOUR_TO_GRID(for_gr_a_1(1,1,lan),for_gr_a_2(1,1,lan),& + lon_dim_a,lonf,lons_lat,2*nlevs) + enddo + + uug = 0.; vvg = 0. + do lan=1,lats_node_a + lat = global_lats_a(ipt_lats_node_a-1+lan) + lons_lat = lonsperlar(lat) + tx1 = 1. / coslat_a(lat) + do k=1,nlevs + do i=1,lons_lat + uug(i,lan,k) = for_gr_a_2(i,k,lan) * tx1 + vvg(i,lan,k) = for_gr_a_2(i,nlevs+k,lan) * tx1 + enddo + enddo + enddo + + return + end subroutine vrtdivspect_to_uvgrid +end module get_stochy_pattern_mod diff --git a/stochastic_physics/getcon_lag_stochy.f b/stochastic_physics/getcon_lag_stochy.f new file mode 100644 index 000000000..e90f51279 --- /dev/null +++ b/stochastic_physics/getcon_lag_stochy.f @@ -0,0 +1,89 @@ + module getcon_lag_stochy_mod + + implicit none + + contains + + subroutine getcon_lag_stochy(lats_nodes_a,global_lats_a, + & lats_nodes_h,global_lats_h_sn, + & lonsperlat,xhalo,yhalo) + use stochy_resol_def, only : jcap,latg,latg2,lonf + use spectral_layout_mod, only : me,nodes + + use stochy_gg_def, only : colrad_a,sinlat_a + use stochy_layout_lag, only : + & ipt_lats_node_h,lat1s_h,lats_dim_h, + & lats_node_h,lats_node_h_max,lon_dim_h + use setlats_lag_stochy_mod, only: setlats_lag_stochy + implicit none +! + integer yhalo,xhalo +! + integer, dimension(nodes) :: lats_nodes_a, lats_nodes_h + integer, dimension(latg) :: lonsperlat, global_lats_a + + integer, dimension(latg+2*yhalo*nodes) :: global_lats_h_sn +! + integer i,j,l,n,lat,i1,i2,node,nodesio + integer, dimension(latg+2*yhalo*nodes) :: global_lats_h_ns +! + if (me == 0) print 100, jcap, me +100 format ('getcon_h jcap= ',i4,2x,'me=',i3) + + do lat = 1, latg2 + lonsperlat(latg+1-lat) = lonsperlat(lat) + end do + nodesio = nodes + +! print*,'con_h me,nodes,nodesio = ',me,nodes,nodesio + + call setlats_lag_stochy(lats_nodes_a,global_lats_a, + & lats_nodes_h,global_lats_h_ns,yhalo) + +! reverse order for use in set_halos + + i1 = 1 + i2 = 0 + do n=1,nodes + j = 0 + i2 = i2 + lats_nodes_h(n) + do i=i1,i2 + j = j + 1 + global_lats_h_sn(i) = global_lats_h_ns(i2+1-j) + enddo + i1 = i2 + 1 + enddo + + 830 format(10(i4,1x)) + lats_dim_h = 0 + do node=1,nodes + lats_dim_h = max(lats_dim_h, lats_nodes_h(node)) + enddo + lats_node_h = lats_nodes_h(me+1) + lats_node_h_max = 0 + do i=1,nodes + lats_node_h_max = max(lats_node_h_max, lats_nodes_h(i)) + enddo + ipt_lats_node_h = 1 + if ( me > 0 ) then + do node=1,me + ipt_lats_node_h = ipt_lats_node_h + lats_nodes_h(node) + enddo + endif + do j=1,latg2 + sinlat_a(j) = cos(colrad_a(j)) + enddo + do l=0,jcap + do lat = 1, latg2 + if ( l <= min(jcap,lonsperlat(lat)/2) ) then + lat1s_h(l) = lat + go to 200 + endif + end do + 200 continue + end do + lon_dim_h = lonf + 1 + xhalo + xhalo !even/odd + return + end + + end module getcon_lag_stochy_mod diff --git a/stochastic_physics/getcon_spectral.F90 b/stochastic_physics/getcon_spectral.F90 new file mode 100644 index 000000000..7eaa48ae8 --- /dev/null +++ b/stochastic_physics/getcon_spectral.F90 @@ -0,0 +1,277 @@ +module getcon_spectral_mod + + implicit none + + contains + + subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & + lats_nodes_a,global_lats_a, & + lonsperlat,latsmax, & + lats_nodes_ext,global_lats_ext, & + epse,epso,epsedn,epsodn, & + snnp1ev,snnp1od, & + plnev_a,plnod_a,pddev_a,pddod_a, & + plnew_a,plnow_a,colat1) + +! program log: +! 20110220 henry juang update code to fit mass_dp and ndslfv +! + use epslon_stochy_mod, only: epslon_stochy + use get_lats_node_a_stochy_mod, only: get_lats_node_a_stochy + use get_ls_node_stochy_mod, only: get_ls_node_stochy + use glats_stochy_mod, only: glats_stochy + use gozrineo_a_stochy_mod, only: gozrineo_a_stochy + use pln2eo_a_stochy_mod, only: pln2eo_a_stochy + use setlats_a_stochy_mod, only: setlats_a_stochy + use stochy_resol_def + use spectral_layout_mod + use stochy_gg_def + use stochy_internal_state_mod + + implicit none +! + integer i,j,k,l,lat,lan,lons_lat,n + integer ls_node(ls_dim,3),ierr +! +! ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +! ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +! ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +! + integer ls_nodes(ls_dim,nodes) + integer, dimension(nodes) :: max_ls_nodes, lats_nodes_a + integer, dimension(latg) :: global_lats_a, lonsperlat +! + integer lats_nodes_ext(nodes) + integer global_lats_ext(latg+2*jintmx+2*nypt*(nodes-1)) +! + real(kind=kind_dbl_prec), dimension(len_trie_ls) :: epse, epsedn, snnp1ev + real(kind=kind_dbl_prec), dimension(len_trio_ls) :: epso, epsodn, snnp1od +! + real(kind=kind_dbl_prec), dimension(len_trie_ls,latg2) :: plnev_a, pddev_a, plnew_a + real(kind=kind_dbl_prec), dimension(len_trio_ls,latg2) :: plnod_a, pddod_a, plnow_a +! + real(kind=kind_dbl_prec), allocatable:: colrad_dp(:), wgt_dp(:),& + wgtcs_dp(:), rcs2_dp(:), epse_dp(:), epso_dp(:),& + epsedn_dp(:), epsodn_dp(:),plnev_dp(:), plnod_dp(:),& + pddev_dp(:), pddod_dp(:),plnew_dp(:), plnow_dp(:) +! + integer iprint,locl,node,& + len_trie_ls_nod, len_trio_ls_nod,& + indev, indod, indlsev,jbasev,indlsod,jbasod +! + integer gl_lats_index, latsmax + integer global_time_sort_index_a(latg) +! + real fd2 +!! + include 'function2' +! + real(kind=kind_dbl_prec) global_time_a(latg) +! + real(kind=kind_dbl_prec), parameter :: cons0 = 0.d0, cons0p5 = 0.5d0,& + cons1 = 1.d0, cons0p92 = 0.92d0 + real(kind=kind_dbl_prec) colat1 +! + gl_lats_index = 0 + global_lats_a = -1 + do lat = 1,latg !my intialize global_time_a to lonsperlat + global_time_a(lat) = lonsperlat(lat) + enddo + + do lat = 1, latg2 + lonsperlat(latg+1-lat) = lonsperlat(lat) + end do + do node=1,nodes + call get_lats_node_a_stochy( node-1, global_lats_a,lats_nodes_a(node),& + gl_lats_index,global_time_sort_index_a, iprint) + enddo + call setlats_a_stochy(lats_nodes_a,global_lats_a,iprint, lonsperlat) + + iprint = 0 + do node=1,nodes + call get_ls_node_stochy( node-1, ls_nodes(1,node),max_ls_nodes(node), iprint ) + enddo +! + len_trie_ls_max = 0 + len_trio_ls_max = 0 + do node=1,nodes +! + len_trie_ls_nod = 0 + len_trio_ls_nod = 0 + do locl=1,max_ls_nodes(node) + l=ls_nodes(locl,node) + len_trie_ls_nod = len_trie_ls_nod+(jcap+3-l)/2 + len_trio_ls_nod = len_trio_ls_nod+(jcap+2-l)/2 + enddo + len_trie_ls_max = max(len_trie_ls_max,len_trie_ls_nod) + len_trio_ls_max = max(len_trio_ls_max,len_trio_ls_nod) +! + enddo +! + iprint = 0 +! + lats_dim_a = 0 + do node=1,nodes + lats_dim_a = max(lats_dim_a,lats_nodes_a(node)) + enddo + lats_node_a = lats_nodes_a(me+1) + + lats_node_a_max = 0 + do i=1,nodes + lats_node_a_max = max(lats_node_a_max, lats_nodes_a(i)) + enddo + latsmax = lats_node_a_max + +! + ipt_lats_node_ext = 1 +! + ipt_lats_node_a = 1 + if ( me > 0 ) then + do node=1,me + ipt_lats_node_a = ipt_lats_node_a + lats_nodes_a(node) + enddo + endif + +! + iprint = 0 +! + if ( kind_dbl_prec == 8 ) then !------------------------------------ + call glats_stochy(latg2,colrad_a,wgt_a,wgtcs_a,rcs2_a,iprint) + call epslon_stochy(epse,epso,epsedn,epsodn,ls_node) + call pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a,ls_node,latg2) + call gozrineo_a_stochy(plnev_a,plnod_a,pddev_a,pddod_a, & + plnew_a,plnow_a,epse,epso,rcs2_a,wgt_a,ls_node,latg2) +! + else !------------------------------------------------------------ + allocate ( colrad_dp(latg2) ) + allocate ( wgt_dp(latg2) ) + allocate ( wgtcs_dp(latg2) ) + allocate ( rcs2_dp(latg2) ) +! + allocate ( epse_dp(len_trie_ls) ) + allocate ( epso_dp(len_trio_ls) ) + allocate ( epsedn_dp(len_trie_ls) ) + allocate ( epsodn_dp(len_trio_ls) ) +! + allocate ( plnev_dp(len_trie_ls) ) + allocate ( plnod_dp(len_trio_ls) ) + allocate ( pddev_dp(len_trie_ls) ) + allocate ( pddod_dp(len_trio_ls) ) + allocate ( plnew_dp(len_trie_ls) ) + allocate ( plnow_dp(len_trio_ls) ) + + call glats_stochy(latg2,colrad_dp,wgt_dp,wgtcs_dp,rcs2_dp,iprint) +! + do i=1,latg2 + colrad_a(i) = colrad_dp(i) + wgt_a(i) = wgt_dp(i) + wgtcs_a(i) = wgtcs_dp(i) + rcs2_a(i) = rcs2_dp(i) + enddo +! + call epslon_stochy(epse_dp,epso_dp,epsedn_dp,epsodn_dp,ls_node) +! + do i=1,len_trie_ls + epse(i) = epse_dp(i) + epsedn(i) = epsedn_dp(i) + enddo +! + do i=1,len_trio_ls + epso(i) = epso_dp(i) + epsodn(i) = epsodn_dp(i) + enddo +! + do lat=1,latg2 +! + call pln2eo_a_stochy(plnev_dp,plnod_dp,epse_dp,epso_dp,colrad_dp(lat),ls_node,1) +! + call gozrineo_a_stochy(plnev_dp,plnod_dp,pddev_dp,pddod_dp, plnew_dp,plnow_dp,& + epse_dp,epso_dp,rcs2_dp(lat),wgt_dp(lat),ls_node,1) +! + do i=1,len_trie_ls + plnev_a(i,lat) = plnev_dp(i) + pddev_a(i,lat) = pddev_dp(i) + plnew_a(i,lat) = plnew_dp(i) + enddo + do i=1,len_trio_ls + plnod_a(i,lat) = plnod_dp(i) + pddod_a(i,lat) = pddod_dp(i) + plnow_a(i,lat) = plnow_dp(i) + enddo + enddo +! + deallocate ( colrad_dp, wgt_dp, wgtcs_dp, rcs2_dp , & + epse_dp, epso_dp, epsedn_dp, epsodn_dp, & + plnev_dp, plnod_dp, pddev_dp, pddod_dp , & + plnew_dp, plnow_dp ) + endif !----------------------------------------------------------- +! +! + do locl=1,ls_max_node + l = ls_node(locl,1) + jbasev = ls_node(locl,2) + indev = indlsev(l,l) + do n = l, jcap, 2 + snnp1ev(indev) = n*(n+1) + indev = indev+1 + end do + end do +! +! + do locl=1,ls_max_node + l = ls_node(locl,1) + jbasod = ls_node(locl,3) + if ( l <= jcap-1 ) then + indod = indlsod(l+1,l) + do n = l+1, jcap, 2 + snnp1od(indod) = n*(n+1) + indod = indod+1 + end do + end if + end do +! +! + do locl=1,ls_max_node + l = ls_node(locl,1) + jbasev = ls_node(locl,2) + jbasod = ls_node(locl,3) + if (mod(L,2) == mod(jcap+1,2)) then ! set even (n-l) terms of top row to zero + snnp1ev(indlsev(jcap+1,l)) = cons0 + else ! set odd (n-l) terms of top row to zero + snnp1od(indlsod(jcap+1,l)) = cons0 + endif + enddo +! + do j=1,latg + if( j <= latg2 ) then + sinlat_a(j) = cos(colrad_a(j)) + else + sinlat_a(j) = -cos(colrad_a(latg+1-j)) + endif + coslat_a(j) = sqrt(1.-sinlat_a(j)*sinlat_a(j)) + enddo +! + do L=0,jcap + do lat = 1, latg2 + if ( L <= min(jcap,lonsperlat(lat)/2) ) then + lat1s_a(L) = lat + go to 200 + endif + end do + 200 continue + end do +! + + do j=1,lats_node_a + lat = global_lats_a(ipt_lats_node_a-1+j) + if ( lonsperlat(lat) == lonf ) then + lon_dims_a(j) = lonfx + else + lon_dims_a(j) = lonsperlat(lat) + 2 + endif + enddo +! + return + end + +end module getcon_spectral_mod diff --git a/stochastic_physics/glats_stochy.f b/stochastic_physics/glats_stochy.f new file mode 100644 index 000000000..4ed9d38f4 --- /dev/null +++ b/stochastic_physics/glats_stochy.f @@ -0,0 +1,109 @@ + module glats_stochy_mod + + implicit none + + contains + + subroutine glats_stochy(lgghaf,colrad,wgt,wgtcs,rcs2,iprint) +! +! Jan 2013 Henry Juang increase precision by kind_qdt_prec=16 +! to help wgt (Gaussian weighting) + use machine + implicit none + integer iter,k,k1,l2,lgghaf,iprint +! +! increase precision for more significant digit to help wgt + real(kind=kind_qdt_prec) drad,dradz,p1,p2,phi,pi,rad,rc +! real(kind=kind_qdt_prec) drad,dradz,eps,p1,p2,phi,pi,rad,rc + real(kind=kind_qdt_prec) rl2,scale,si,sn,w,x +! + real(kind=kind_dbl_prec), dimension(lgghaf) :: colrad, wgt, + & wgtcs, rcs2 +! + real(kind=kind_dbl_prec), parameter :: cons0 = 0.d0, cons1 = 1.d0, + & cons2 = 2.d0, cons4 = 4.d0, + & cons180 = 180.d0, + & cons360 = 360.d0, + & cons0p25 = 0.25d0 + real(kind=kind_qdt_prec), parameter :: eps = 1.d-20 +! +! for better accuracy to select smaller number +! eps = 1.d-12 +! eps = 1.d-20 +! + if(iprint == 1) print 101 + 101 format (' i colat colrad wgt', 12x, 'wgtcs', + & 10x, 'iter res') + si = cons1 + l2 = 2*lgghaf + rl2 = l2 + scale = cons2/(rl2*rl2) + k1 = l2-1 + pi = atan(si)*cons4 +! dradz = pi / cons360 / 10.0 +! for better accuracy to start iteration + dradz = pi / float(lgghaf) / 200.0 + rad = cons0 + do k=1,lgghaf + iter = 0 + drad = dradz +1 call poly(l2,rad,p2) +2 p1 = p2 + iter = iter + 1 + rad = rad + drad + call poly(l2,rad,p2) + if(sign(si,p1) == sign(si,p2)) go to 2 + if(drad < eps)go to 3 + rad = rad-drad + drad = drad * cons0p25 + go to 1 +3 continue + colrad(k) = rad + phi = rad * cons180 / pi + call poly(k1,rad,p1) + x = cos(rad) + w = scale * (cons1 - x*x)/ (p1*p1) + wgt(k) = w + sn = sin(rad) + w = w/(sn*sn) + wgtcs(k) = w + rc = cons1/(sn*sn) + rcs2(k) = rc + call poly(l2,rad,p1) + if(iprint == 1) + & print 102,k,phi,colrad(k),wgt(k),wgtcs(k),iter,p1 + 102 format(1x,i3,2x,f6.2,2x,f10.7,2x,e14.7,2x,e14.7,2x,i4,2x,e14.7) + enddo + if(iprint == 1) print 100,lgghaf +100 format(1h ,'shalom from 0.0e0 glats for ',i3) +! + return + end + + subroutine poly(n,rad,p) + use machine +! + implicit none +! + integer i,n +! +! increase precision for more significant digit to help wgt + real(kind=kind_qdt_prec) floati,g,p,rad,x,y1,y2,y3 +! + real(kind=kind_dbl_prec), parameter :: cons1 = 1.d0 +! + x = cos(rad) + y1 = cons1 + y2 = x + do i=2,n + g = x*y2 + floati = i + y3 = g - y1 + g - (g-y1)/floati + y1 = y2 + y2 = y3 + enddo + p = y3 + return + end + + end module glats_stochy_mod diff --git a/stochastic_physics/gozrineo_stochy.f b/stochastic_physics/gozrineo_stochy.f new file mode 100644 index 000000000..01210ff92 --- /dev/null +++ b/stochastic_physics/gozrineo_stochy.f @@ -0,0 +1,180 @@ + module gozrineo_a_stochy_mod + + implicit none + + contains + + subroutine gozrineo_a_stochy(plnev_a,plnod_a, + & pddev_a,pddod_a, + & plnew_a,plnow_a, + & epse,epso,rcs2_a,wgt_a,ls_node,num_lat) +cc + use stochy_resol_def + use spectral_layout_mod + use machine + implicit none +cc + real(kind=kind_dbl_prec) plnev_a(len_trie_ls,latg2) + real(kind=kind_dbl_prec) plnod_a(len_trio_ls,latg2) + real(kind=kind_dbl_prec) pddev_a(len_trie_ls,latg2) + real(kind=kind_dbl_prec) pddod_a(len_trio_ls,latg2) + real(kind=kind_dbl_prec) plnew_a(len_trie_ls,latg2) + real(kind=kind_dbl_prec) plnow_a(len_trio_ls,latg2) +cc + real(kind=kind_dbl_prec) epse(len_trie_ls) + real(kind=kind_dbl_prec) epso(len_trio_ls) +cc + real(kind=kind_dbl_prec) rcs2_a(latg2) + real(kind=kind_dbl_prec) wgt_a(latg2) +cc + integer ls_node(ls_dim,3) +cc + integer num_lat +cc +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +cc + integer l,lat,locl,n +cc + integer indev,indev1,indev2 + integer indod,indod1,indod2 + integer inddif +cc + real(kind=kind_dbl_prec) rn,rnp1,wcsa +cc + real(kind=kind_dbl_prec) cons0 !constant + real(kind=kind_dbl_prec) cons2 !constant + real rerth +cc + integer indlsev,jbasev + integer indlsod,jbasod +cc + include 'function2' +cc +cc + cons0 = 0.d0 !constant + cons2 = 2.d0 !constant + rerth =6.3712e+6 ! radius of earth (m) +cc +cc + do lat=1,num_lat +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) +cc + rn=l +cc +cc + pddev_a(indlsev(l,l),lat) = -epso(indlsod(l+1,l)) + & * plnod_a(indlsod(l+1,l),lat) * rn + indev1 = indlsev(L,L) + 1 + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap-1,L) + else + indev2 = indlsev(jcap ,L) + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + rn =l+2 + rnp1=l+2+1 + do indev = indev1 , indev2 +cc + pddev_a(indev,lat) = epse(indev) + & * plnod_a(indev-inddif ,lat) * rnp1 + & - epso(indev-inddif+1) + & * plnod_a(indev-inddif+1,lat) * rn +cc + rn = rn + cons2 !constant + rnp1 = rnp1 + cons2 !constant + enddo +cc +cc...................................................................... + indev1 = indlsev(L,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) - 1 + else + indev2 = indlsev(jcap ,L) - 1 + endif + indod1 = indlsod(l+1,l) + inddif = indev1 - indod1 +cc + rn =l+1 + rnp1=l+1+1 + do indev = indev1 , indev2 +cc + pddod_a(indev-inddif,lat) = epso(indev-inddif) + & * plnev_a(indev ,lat) * rnp1 + & - epse(indev+1) + & * plnev_a(indev+1,lat) * rn +cc + rn = rn + cons2 !constant + rnp1 = rnp1 + cons2 !constant + enddo +cc + enddo +cc +cc...................................................................... +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) +cc + if (mod(L,2).eq.mod(jcap+1,2)) then +cc +cc set the even (n-l) terms of the top row to zero + pddev_a(indlsev(jcap+1,l),lat) = cons0 !constant +cc + else +cc +cc set the odd (n-l) terms of the top row to zero + pddod_a(indlsod(jcap+1,l),lat) = cons0 !constant +cc + endif +cc + enddo +cc +cc...................................................................... +cc + wcsa=rcs2_a(lat)/rerth +cc +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + indev1 = indlsev(L,L) + indod1 = indlsod(L+1,L) + if (mod(L,2).eq.mod(jcap+1,2)) then + indev2 = indlsev(jcap+1,L) + indod2 = indlsod(jcap ,L) + else + indev2 = indlsev(jcap ,L) + indod2 = indlsod(jcap+1,L) + endif + do indev = indev1 , indev2 +cc + pddev_a(indev,lat) = pddev_a(indev,lat) * wcsa + plnew_a(indev,lat) = plnev_a(indev,lat) * wgt_a(lat) +cc + enddo +cc + do indod = indod1 , indod2 +cc + pddod_a(indod,lat) = pddod_a(indod,lat) * wcsa + plnow_a(indod,lat) = plnod_a(indod,lat) * wgt_a(lat) +cc + enddo +cc + enddo +cc + enddo +cc + return + end + + end module gozrineo_a_stochy_mod diff --git a/stochastic_physics/initialize_spectral_mod.F90 b/stochastic_physics/initialize_spectral_mod.F90 new file mode 100644 index 000000000..13d1dd2f0 --- /dev/null +++ b/stochastic_physics/initialize_spectral_mod.F90 @@ -0,0 +1,268 @@ +! !module: stochy_initialize_spectral +! --- initialize module of the +! gridded component of the stochastic physics patteern +! generator, which is in spectral space +! +! !description: gfs dynamics gridded component initialize module. +! +! !revision history: +! +! oct 11 2016 P.Pegion copy of gsm/dynamics to create stand alone version +! +! !interface: +! + module initialize_spectral_mod +! +!!uses: +! + use machine + use spectral_layout_mod, only : ipt_lats_node_a, lats_node_a_max,lon_dim_a,len_trie_ls,len_trio_ls & + ,nodes,ls_max_node,lats_dim_a,ls_dim,nodes_comp,lat1s_a + use stochy_layout_lag, only : lat1s_h + use stochy_internal_state_mod + use spectral_layout_mod,only:lon_dims_a + use stochy_resol_def + use stochy_namelist_def + use stochy_ccpp, only : is_master, num_parthds_stochy => ompthreads + use stochy_gg_def, only : wgt_a,sinlat_a,coslat_a,colrad_a,wgtcs_a,rcs2_a,lats_nodes_h,global_lats_h + use getcon_spectral_mod, only: getcon_spectral + use get_ls_node_stochy_mod, only: get_ls_node_stochy + use getcon_lag_stochy_mod, only: getcon_lag_stochy +#ifndef IBM + USE omp_lib +#endif + + implicit none + + contains + + subroutine initialize_spectral(gis_stochy, rc) + +! this subroutine set up the internal state variables, +! allocate internal state arrays for initializing the gfs system. +!---------------------------------------------------------------- +! + implicit none +! +! type(stochy_internal_state), pointer, intent(inout) :: gis_stochy + type(stochy_internal_state), intent(inout) :: gis_stochy + integer, intent(out) :: rc + integer :: ierr, npe_single_member, iret,latghf + integer :: i, j, k, l, n, locl + logical :: file_exists=.false. + integer, parameter :: iunit=101 + +!------------------------------------------------------------------- + +! set up gfs internal state dimension and values for dynamics etc +!------------------------------------------------------------------- +! print*,'before allocate lonsperlat,',& +! allocated(gis_stochy%lonsperlat),'latg=',latg +! +! gis_stochy%nodes=mpp_npes() +! print*,'mpp_npes=',mpp_npes() + nodes = gis_stochy%nodes + npe_single_member = gis_stochy%npe_single_member + + lon_dim_a = lon_s + 2 + jcap=ntrunc + jcap1 = jcap+1 + jcap2 = jcap+2 + latg = lat_s + latg2 = latg/2 + lonf = lon_s + lnt = jcap2*jcap1/2 + lnuv = jcap2*jcap1 + lnt2 = lnt + lnt + lnt22 = lnt2 + 1 + lnte = (jcap2/2)*((jcap2/2)+1)-1 + lnto = (jcap2/2)*((jcap2/2)+1)-(jcap2/2) + lnted = lnte + lntod = lnto + + gis_stochy%lnt2 = lnt2 + + allocate(lat1s_a(0:jcap)) + allocate(lon_dims_a(latg)) + + allocate(wgt_a(latg2)) + allocate(wgtcs_a(latg2)) + allocate(rcs2_a(latg2)) + +!! create io communicator and comp communicator +!! + nodes_comp=nodes +! +! if (is_master()) then +! print*,'number of threads is',num_parthds_stochy +! print*,'number of mpi procs is',nodes +! endif +! + ls_dim = (jcap1-1)/nodes+1 +! print*,'allocating lonsperlat',latg + allocate(gis_stochy%lonsperlat(latg)) +! print*,'size=',size(gis_stochy%lonsperlat) + + + inquire (file="lonsperlat.dat", exist=file_exists) + if ( .not. file_exists ) then + !call mpp_error(FATAL,'Requested lonsperlat.dat data file does not exist') + gis_stochy%lonsperlat(:)=lonf + else + open (iunit,file='lonsperlat.dat',status='old',form='formatted', & + action='read',iostat=iret) + if (iret /= 0) then + write(0,*) 'error while reading lonsperlat.dat' + rc = 1 + return + end if + rewind iunit + read (iunit,*,iostat=iret) latghf,(gis_stochy%lonsperlat(i),i=1,latghf) + if (latghf+latghf /= latg) then + write(0,*)' latghf=',latghf,' not equal to latg/2=',latg/2 + if (iret /= 0) then + write(0,*) 'lonsperlat file has wrong size' + rc = 1 + return + end if + endif + do i=1,latghf + gis_stochy%lonsperlat(latg-i+1) = gis_stochy%lonsperlat(i) + enddo + close(iunit) + endif +!! +!cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx +! +! write(0,*)'before allocate ls_nodes,',allocated(gis_stochy%ls_nodes),& +! 'ls_dim=', ls_dim,'nodes=',nodes + allocate ( gis_stochy%ls_node (ls_dim*3) ) + allocate ( gis_stochy%ls_nodes(ls_dim,nodes) ) + allocate ( gis_stochy%max_ls_nodes(nodes) ) +! + allocate ( gis_stochy%lats_nodes_a_fix(nodes)) ! added for mGrid +! + allocate ( gis_stochy%lats_nodes_a(nodes) ) + allocate ( gis_stochy%global_lats_a(latg) ) +! + allocate ( gis_stochy%lats_nodes_ext(nodes) ) + allocate ( gis_stochy%global_lats_ext(latg+2*jintmx+2*nypt*(nodes-1)) ) + +! internal parallel structure. Weiyu. +!--------------------------------------------------- + ALLOCATE(gis_stochy%TRIE_LS_SIZE (npe_single_member)) + ALLOCATE(gis_stochy%TRIO_LS_SIZE (npe_single_member)) + ALLOCATE(gis_stochy%TRIEO_LS_SIZE (npe_single_member)) + ALLOCATE(gis_stochy%LS_MAX_NODE_GLOBAL(npe_single_member)) + ALLOCATE(gis_stochy%LS_NODE_GLOBAL (LS_DIM*3, npe_single_member)) + + gis_stochy%LS_NODE_GLOBAL = 0 + gis_stochy%LS_MAX_NODE_GLOBAL = 0 + gis_stochy%TRIEO_TOTAL_SIZE = 0 + + DO i = 1, npe_single_member + CALL GET_LS_NODE_STOCHY(i-1, gis_stochy%LS_NODE_GLOBAL(1, i), & + gis_stochy%LS_MAX_NODE_GLOBAL(i), gis_stochy%IPRINT) + gis_stochy%TRIE_LS_SIZE(i) = 0 + gis_stochy%TRIO_LS_SIZE(i) = 0 + DO LOCL = 1, gis_stochy%LS_MAX_NODE_GLOBAL(i) + gis_stochy%LS_NODE_GLOBAL(LOCL+ LS_DIM, i) = gis_stochy%TRIE_LS_SIZE(i) + gis_stochy%LS_NODE_GLOBAL(LOCL+ 2*LS_DIM, i) = gis_stochy%TRIO_LS_SIZE(i) + + L = gis_stochy%LS_NODE_GLOBAL(LOCL, i) + + gis_stochy%TRIE_LS_SIZE(i) = gis_stochy%TRIE_LS_SIZE(i) + (JCAP+3-L)/2 + gis_stochy%TRIO_LS_SIZE(i) = gis_stochy%TRIO_LS_SIZE(i) + (JCAP+2-L)/2 + END DO + gis_stochy%TRIEO_LS_SIZE(i) = gis_stochy%TRIE_LS_SIZE(i) + gis_stochy%TRIO_LS_SIZE(i) + 3 + gis_stochy%TRIEO_TOTAL_SIZE = gis_stochy%TRIEO_TOTAL_SIZE + gis_stochy%TRIEO_LS_SIZE(i) + END DO + + +!--------------------------------------------------- +! + gis_stochy%iprint = 0 + call get_ls_node_stochy( gis_stochy%me, gis_stochy%ls_node, ls_max_node, gis_stochy%iprint ) +! +! + len_trie_ls = 0 + len_trio_ls = 0 + do locl=1,ls_max_node + gis_stochy%ls_node(locl+ ls_dim) = len_trie_ls + gis_stochy%ls_node(locl+2*ls_dim) = len_trio_ls + l = gis_stochy%ls_node(locl) + len_trie_ls = len_trie_ls+(jcap+3-l)/2 + len_trio_ls = len_trio_ls+(jcap+2-l)/2 + enddo +! if (gis_stochy%me == 0) print *,'ls_node=',gis_stochy%ls_node(1:ls_dim),'2dim=', & +! gis_stochy%ls_node(ls_dim+1:2*ls_dim),'3dim=', & +! gis_stochy%ls_node(2*ls_dim+1:3*ls_dim) +! +! + allocate ( gis_stochy%epse (len_trie_ls) ) + allocate ( gis_stochy%epso (len_trio_ls) ) + allocate ( gis_stochy%epsedn(len_trie_ls) ) + allocate ( gis_stochy%epsodn(len_trio_ls) ) + allocate ( gis_stochy%kenorm_e(len_trie_ls) ) + allocate ( gis_stochy%kenorm_o(len_trio_ls) ) +! + allocate ( gis_stochy%snnp1ev(len_trie_ls) ) + allocate ( gis_stochy%snnp1od(len_trio_ls) ) +! + allocate ( gis_stochy%plnev_a(len_trie_ls,latg2) ) + allocate ( gis_stochy%plnod_a(len_trio_ls,latg2) ) + allocate ( gis_stochy%pddev_a(len_trie_ls,latg2) ) + allocate ( gis_stochy%pddod_a(len_trio_ls,latg2) ) + allocate ( gis_stochy%plnew_a(len_trie_ls,latg2) ) + allocate ( gis_stochy%plnow_a(len_trio_ls,latg2) ) + + allocate(colrad_a(latg2)) + allocate(sinlat_a(latg)) + allocate(coslat_a(latg)) + allocate(lat1s_h(0:jcap)) +! + if(gis_stochy%iret/=0) then + write(0,*) 'incompatible namelist - aborted in stochy' + rc = 1 + return + end if +!! + gis_stochy%lats_nodes_ext = 0 + call getcon_spectral(gis_stochy%ls_node, gis_stochy%ls_nodes, & + gis_stochy%max_ls_nodes, gis_stochy%lats_nodes_a, & + gis_stochy%global_lats_a, gis_stochy%lonsperlat, & + gis_stochy%lats_node_a_max, gis_stochy%lats_nodes_ext, & + gis_stochy%global_lats_ext, gis_stochy%epse, & + gis_stochy%epso, gis_stochy%epsedn, & + gis_stochy%epsodn, gis_stochy%snnp1ev, & + gis_stochy%snnp1od, gis_stochy%plnev_a, & + gis_stochy%plnod_a, gis_stochy%pddev_a, & + gis_stochy%pddod_a, gis_stochy%plnew_a, & + gis_stochy%plnow_a, gis_stochy%colat1) +! + gis_stochy%lats_node_a = gis_stochy%lats_nodes_a(gis_stochy%me+1) + gis_stochy%ipt_lats_node_a = ipt_lats_node_a + +! if (gis_stochy%me == 0) & +! write(0,*)'after getcon_spectral lats_node_a=',gis_stochy%lats_node_a & +! ,'ipt_lats_node_a=',gis_stochy%ipt_lats_node_a +! + if (.not. allocated(lats_nodes_h)) allocate (lats_nodes_h(nodes)) + if (.not. allocated(global_lats_h)) allocate (global_lats_h(latg+2*gis_stochy%yhalo*nodes)) + call getcon_lag_stochy(gis_stochy%lats_nodes_a,gis_stochy%global_lats_a, & + lats_nodes_h, global_lats_h, & + gis_stochy%lonsperlat,gis_stochy%xhalo,gis_stochy%yhalo) + +! +! + allocate ( gis_stochy%trie_ls (len_trie_ls,2,lotls) ) + allocate ( gis_stochy%trio_ls (len_trio_ls,2,lotls) ) + +! if (gis_stochy%me == 0) then +! print*, ' lats_dim_a=', lats_dim_a, ' lats_node_a=', gis_stochy%lats_node_a +! endif + rc=0 + + end subroutine initialize_spectral + + end module initialize_spectral_mod diff --git a/stochastic_physics/pln2eo_stochy.f b/stochastic_physics/pln2eo_stochy.f new file mode 100644 index 000000000..4c6576ba6 --- /dev/null +++ b/stochastic_physics/pln2eo_stochy.f @@ -0,0 +1,287 @@ + module pln2eo_a_stochy_mod + + implicit none + + contains + + subroutine pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a, + & ls_node,num_lat) +! +! use x-number method to archieve accuracy due to recursive to avoid +! underflow and overflow if necessary by henry juang 2012 july +! + use stochy_resol_def + use spectral_layout_mod + use machine + implicit none +! +! define x number constant for real8 start + integer, parameter :: in_f = 960 , in_h = in_f/2 + real(kind=kind_dbl_prec), parameter :: bb_f = 2.d0 ** ( in_f ) + real(kind=kind_dbl_prec), parameter :: bs_f = 2.d0 ** (-in_f ) + real(kind=kind_dbl_prec), parameter :: bb_h = 2.d0 ** ( in_h ) + real(kind=kind_dbl_prec), parameter :: bs_h = 2.d0 ** (-in_h ) +! define x number constant end + +cc + real(kind=kind_dbl_prec) plnev_a(len_trie_ls,latg2) + real(kind=kind_dbl_prec) plnod_a(len_trio_ls,latg2) +cc + real(kind=kind_dbl_prec) epse(len_trie_ls) + real(kind=kind_dbl_prec) epso(len_trio_ls) +cc + real(kind=kind_dbl_prec) colrad_a(latg2) +cc + integer ls_node(ls_dim,3) +cc + integer num_lat +cc +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +cc + integer l,lat,locl,max_l,n +cc + integer indev + integer indod +cc +! need index for alp to be x-number + integer id, ialp1, ialp2, ialp3, iprod + integer ialp10(0:jcap) + real(kind=kind_dbl_prec) aa, bb, w + + real(kind=kind_dbl_prec) a,alp1,alp2,alp3,b + real(kind=kind_dbl_prec) cos2,fl,prod,sinlat,coslat +cc + real(kind=kind_dbl_prec) alp10(0:jcap) +cc + real(kind=kind_dbl_prec) cons0,cons0p5,cons1,cons2,cons3 !constant +cc +cc + integer indlsev,jbasev + integer indlsod,jbasod +cc + include 'function2' +cc +cc + cons0=0.0d0 !constant + cons0p5=0.5d0 !constant + cons1=1.0d0 !constant + cons2=2.0d0 !constant + cons3=3.0d0 !constant +cc +cc + max_l=-1 + do locl=1,ls_max_node + max_l = max ( max_l, ls_node(locl,1) ) + enddo +cc +cc + do lat=1,num_lat +cc + sinlat = cos(colrad_a(lat)) + cos2=cons1-sinlat*sinlat !constant + coslat = sqrt(cos2) + +! use x number for alp10 + alp10(0) = sqrt(0.5) + ialp10(0) = 0 + + do l=1,max_l + fl = l + prod=coslat*sqrt(cons1+cons1/(cons2*fl)) + iprod=0 + w = abs(prod) + if( w.ge.bb_h ) then + prod = prod * bs_f + iprod = iprod + 1 + elseif( w.lt.bs_h ) then + prod = prod * bb_f + iprod = iprod - 1 + endif + alp10(l)=alp10(l-1)*prod + ialp10(l)=ialp10(l-1)+iprod + w = abs(alp10(l)) + if( w.ge.bb_h ) then + alp10(l) = alp10(l) * bs_f + ialp10(l) = ialp10(l) + 1 + elseif( w.lt.bs_h ) then + alp10(l) = alp10(l) * bb_f + ialp10(l) = ialp10(l) - 1 + endif + enddo +cc + do locl=1,ls_max_node + l=ls_node(locl,1) + jbasev=ls_node(locl,2) + jbasod=ls_node(locl,3) + n=l + fl=l +! get m=normalized x number for alp1 start + alp1=alp10(l) + ialp1=ialp10(l) + + indev=indlsev(n ,l) + indod=indlsod(n+1,l) +! x2f plnev_a(indev ,lat)=alp1 +! x2f start + if( ialp1.eq.0 ) then + plnev_a(indev ,lat)=alp1 + elseif( ialp1.eq.-1 ) then + plnev_a(indev ,lat)=alp1 * bs_f + elseif( ialp1.lt.-1 ) then + plnev_a(indev ,lat)=0.0 +!! plnev_a(indev ,lat)=alp1 * bs_f * bs_f + else + plnev_a(indev ,lat)=alp1 * bb_f + endif +! x2f end + +! xltime alp2=sqrt(cons2*fl+cons3)*sinlat*alp1 !constant +! xltime start + prod=sqrt(cons2*fl+cons3)*sinlat + iprod=0 + w = abs(prod) + if( w.ge.bb_h ) then + prod = prod * bs_f + iprod = iprod + 1 + elseif( w.lt.bs_h ) then + prod = prod * bb_f + iprod = iprod - 1 + endif + alp2=alp1*prod + ialp2 = ialp1 + iprod +! xltime end +! norm alp2 start + w = abs(alp2) + if( w.ge.bb_h ) then + alp2 = alp2*bs_f + ialp2 = ialp2 + 1 + elseif( w.lt.bs_h ) then + alp2 = alp2*bb_f + ialp2 = ialp2 - 1 + endif +! norm alp2 end + +! x2f plnod_a(indod ,lat)=alp2 +! x2f start + if( ialp2.eq.0 ) then + plnod_a(indod ,lat)=alp2 + elseif( ialp2.eq.-1 ) then + plnod_a(indod ,lat)=alp2 * bs_f + elseif( ialp2.lt.-1 ) then + plnod_a(indod ,lat)=0.0 +!! plnod_a(indod ,lat)=alp2 * bs_f * bs_f + else + plnod_a(indod ,lat)=alp2 * bb_f + endif +! x2f end +cc + do n=l+2,jcap+1 + if(mod(n+l,2).eq.0) then + indev=indev+1 +! xlsum2 start + aa = sinlat / epse(indev) + bb = epso(indod) / epse(indev) + id = ialp2 - ialp1 + if( id.eq.0 ) then + alp3 = aa*alp2 - bb*alp1 + ialp3 = ialp1 + elseif( id.eq.1 ) then + alp3 = aa*alp2 - bb*alp1*bs_f + ialp3 = ialp2 + elseif( id.eq.-1 ) then + alp3 = aa*alp2*bs_f - bb*alp1 + ialp3 = ialp1 + elseif( id.gt.1 ) then + alp3 = aa*alp2 + ialp3 = ialp2 + else + alp3 = - bb*alp1 + ialp3 = ialp1 + endif +! xlsum2 end +! xnorm alp3 start + w = abs(alp3) + if( w.ge.bb_h ) then + alp3 = alp3*bs_f + ialp3 = ialp3 + 1 + elseif( w.lt.bs_h ) then + alp3 = alp3*bb_f + ialp3 = ialp3 - 1 + endif +! xnorm alp3 end + +! x2f alp3 start + if( ialp3.eq.0 ) then + plnev_a(indev,lat)=alp3 + elseif( ialp3.eq.-1 ) then + plnev_a(indev,lat)=alp3 * bs_f + elseif( ialp3.lt.-1 ) then + plnev_a(indev,lat)=0.0 + else + plnev_a(indev,lat)=alp3 * bb_f + endif +! x2f alp3 end + + else + indod=indod+1 + +! xlsum2 start + aa = sinlat / epso(indod) + bb = epse(indev) / epso(indod) + id = ialp2 - ialp1 + if( id.eq.0 ) then + alp3 = aa*alp2 - bb*alp1 + ialp3 = ialp1 + elseif( id.eq.1 ) then + alp3 = aa*alp2 - bb*alp1*bs_f + ialp3 = ialp2 + elseif( id.eq.-1 ) then + alp3 = aa*alp2*bs_f - bb*alp1 + ialp3 = ialp1 + elseif( id.gt.1 ) then + alp3 = aa*alp2 + ialp3 = ialp2 + else + alp3 = - bb*alp1 + ialp3 = ialp1 + endif +! xlsum2 end +! xnorm alp3 start + w = abs(alp3) + if( w.ge.bb_h ) then + alp3 = alp3*bs_f + ialp3 = ialp3 + 1 + elseif( w.lt.bs_h ) then + alp3 = alp3*bb_f + ialp3 = ialp3 - 1 + endif +! xnorm alp3 end + +! x2f alp3 start + if( ialp3.eq.0 ) then + plnod_a(indod,lat)=alp3 + elseif( ialp3.eq.-1 ) then + plnod_a(indod,lat)=alp3 * bs_f + elseif( ialp3.lt.-1 ) then + plnod_a(indod,lat)=0.0 + else + plnod_a(indod,lat)=alp3 * bb_f + endif +! x2f alp3 end + endif + alp1=alp2 + alp2=alp3 + ialp1 = ialp2 + ialp2 = ialp3 + enddo +cc + enddo +cc + enddo +cc + return + end + + end module pln2eo_a_stochy_mod diff --git a/stochastic_physics/setlats_a_stochy.f b/stochastic_physics/setlats_a_stochy.f new file mode 100644 index 000000000..4eb4b876f --- /dev/null +++ b/stochastic_physics/setlats_a_stochy.f @@ -0,0 +1,195 @@ + module setlats_a_stochy_mod + + implicit none + + contains + + subroutine setlats_a_stochy(lats_nodes_a,global_lats_a, + & iprint,lonsperlat) +! + use stochy_resol_def , only : latg,lonf + use spectral_layout_mod , only : nodes,me +! + implicit none +! + integer, dimension(latg) :: global_lats_a, lonsperlat + integer lats_nodes_a(nodes) + + integer iprint,opt,ifin,nodesio + &, jcount,jpt,lat,lats_sum,node,i,ii + &, ngrptg,ngrptl,ipe,irest,idp + &, ngrptgh,nodesioh +! &, ilatpe,ngrptg,ngrptl,ipe,irest,idp +! + integer,allocatable :: lats_hold(:,:) +! + allocate ( lats_hold(latg,nodes) ) +! +! iprint = 1 + iprint = 0 + opt = 1 ! reduced grid + if (opt == 2) lonsperlat = lonf ! full grid + lats_nodes_a = 0 +! if (liope .and. icolor == 2) then +! nodesio = 1 +! else + nodesio = nodes +! endif +! + ngrptg = 0 + do lat=1,latg + do i=1,lonsperlat(lat) + ngrptg = ngrptg + 1 + enddo + enddo + +! +! ngrptg contains total number of grid points. +! +! distribution of the grid + nodesioh = nodesio / 2 + + if (nodesioh*2 /= nodesio) then +! ilatpe = ngrptg / nodesio + ngrptl = 0 + ipe = 0 + irest = 0 + idp = 1 + + do lat=1,latg + ifin = lonsperlat(lat) + ngrptl = ngrptl + ifin + +! if (me == 0) +! &write(2000+me,*)'in setlats lat=',lat,' latg=',latg,' ifin=',ifin +! &,' ngrptl=',ngrptl,' nodesio=',nodesio,' ngrptg=',ngrptg +! &,' irest=',irest + + if (ngrptl*nodesio <= ngrptg+irest) then + lats_nodes_a(ipe+1) = lats_nodes_a(ipe+1) + 1 + lats_hold(idp,ipe+1) = lat + idp = idp + 1 +! if (me == 0) +! & write(2000+me,*)' nodesio1=',nodesio,' idp=',idp,' ipe=',ipe + else + ipe = ipe + 1 + if (ipe <= nodesio) lats_hold(1,ipe+1) = lat + idp = 2 + irest = irest + ngrptg - (ngrptl-ifin)*nodesio + ngrptl = ifin + lats_nodes_a(ipe+1) = lats_nodes_a(ipe+1) + 1 +! if (me == 0) +! & write(2000+me,*)' nodesio1=',nodesio,' idp=',idp,' ipe=',ipe + endif +! if (me == 0) +! & write(2000+me,*)' lat=',lat,' lats_nodes_a=',lats_nodes_a(ipe+1) +! &,' ipe+1=',ipe+1 + enddo + else + nodesioh = nodesio/2 + ngrptgh = ngrptg/2 + ngrptl = 0 + ipe = 0 + irest = 0 + idp = 1 + + do lat=1,latg/2 + ifin = lonsperlat(lat) + ngrptl = ngrptl + ifin + +! if (me == 0) +! &write(0,*)'in setlats lat=',lat,' latg=',latg,' ifin=',ifin +! &,' ngrptl=',ngrptl,' nodesio=',nodesio,' ngrptg=',ngrptg +! &,' irest=',irest,' ngrptgh=',ngrptgh,' nodesioh=',nodesioh + + if (ngrptl*nodesioh <= ngrptgh+irest .or. lat == latg/2) then + lats_nodes_a(ipe+1) = lats_nodes_a(ipe+1) + 1 + lats_hold(idp,ipe+1) = lat +! lats_nodes_a(nodesio-ipe) = lats_nodes_a(nodesio-ipe) + 1 +! lats_hold(idp,nodesio-ipe) = latg+1-lat + idp = idp + 1 +! if (me == 0) +! & write(0,*)' nodesio1=',nodesioh,' idp=',idp,' ipe=',ipe + else + ipe = ipe + 1 + if (ipe <= nodesioh) then + lats_hold(1,ipe+1) = lat +! lats_hold(1,nodesio-ipe) = latg+1-lat + endif + idp = 2 + irest = irest + ngrptgh - (ngrptl-ifin)*nodesioh + ngrptl = ifin + lats_nodes_a(ipe+1) = lats_nodes_a(ipe+1) + 1 +! lats_nodes_a(nodesio-ipe) = lats_nodes_a(nodesio-ipe) + 1 +! if (me == 0) +! & write(0,*)' nodesio1h=',nodesioh,'idp=',idp,' ipe=',ipe + endif +! if (me == 0) +! & write(0,*)' lat=',lat,' lats_nodes_a=',lats_nodes_a(ipe+1) +! &,' ipe+1=',ipe+1 + enddo + do node=1, nodesioh + ii = nodesio-node+1 + jpt = lats_nodes_a(node) + lats_nodes_a(ii) = jpt + do i=1,jpt + lats_hold(jpt+1-i,ii) = latg+1-lats_hold(i,node) + enddo + enddo + + + endif +!! +!!........................................................ +!! + jpt = 0 + do node=1,nodesio +! write(2000+me,*)'node=',node,' lats_nodes_a=',lats_nodes_a(node) +! &, ' jpt=',jpt,' nodesio=',nodesio + if ( lats_nodes_a(node) > 0 ) then + do jcount=1,lats_nodes_a(node) + global_lats_a(jpt+jcount) = lats_hold(jcount,node) +! write(2000+me,*)' jpt+jcount=',jpt+jcount +! &, 'global_lats_a=',global_lats_a(jpt+jcount) + enddo + endif + jpt = jpt + lats_nodes_a(node) + enddo +!! + deallocate (lats_hold) + if ( iprint /= 1 ) return +!! + if (me == 0) then + jpt=0 + do node=1,nodesio + if ( lats_nodes_a(node) > 0 ) then + print 600 + lats_sum=0 + do jcount=1,lats_nodes_a(node) + lats_sum=lats_sum + lonsperlat(global_lats_a(jpt+jcount)) + print 700, node-1, + x node, lats_nodes_a(node), + x jpt+jcount, global_lats_a(jpt+jcount), + x lonsperlat(global_lats_a(jpt+jcount)), + x lats_sum + enddo + endif + jpt=jpt+lats_nodes_a(node) + enddo +! + print 600 +! + 600 format ( ' ' ) +! + 700 format ( 'setlats me=', i4, + x ' lats_nodes_a(', i4, ' )=', i4, + x ' global_lats_a(', i4, ' )=', i4, + x ' lonsperlat=', i5, + x ' lats_sum=', i6 ) +! + endif + + return + end + + end module setlats_a_stochy_mod diff --git a/stochastic_physics/setlats_lag_stochy.f b/stochastic_physics/setlats_lag_stochy.f new file mode 100644 index 000000000..ff4e4f0a6 --- /dev/null +++ b/stochastic_physics/setlats_lag_stochy.f @@ -0,0 +1,127 @@ + module setlats_lag_stochy_mod + + implicit none + + contains + + subroutine setlats_lag_stochy(lats_nodes_a, global_lats_a, + & lats_nodes_h, global_lats_h, yhalo) +! + use stochy_resol_def, only : latg + use spectral_layout_mod, only : me,nodes + implicit none +! + integer yhalo +! + integer lats_nodes_a(nodes), lats_nodes_h(nodes) + &, global_lats_a(latg) + &, global_lats_h(latg+2*yhalo*nodes) +! + integer jj,jpt_a,jpt_h,lat_val,nn,nodes_lats + &, j1, j2, iprint +! + lats_nodes_h = 0 +! + nodes_lats = 0 + do nn=1,nodes + if (lats_nodes_a(nn) > 0) then + lats_nodes_h(nn) = lats_nodes_a(nn) + yhalo + yhalo + nodes_lats = nodes_lats + 1 + endif + enddo +! + global_lats_h = 0 +! +! set non-yhalo latitudes +! + jpt_a = 0 + jpt_h = yhalo + do nn=1,nodes + if (lats_nodes_a(nn) > 0) then + do jj=1,lats_nodes_a(nn) + jpt_a = jpt_a + 1 + jpt_h = jpt_h + 1 + global_lats_h(jpt_h) = global_lats_a(jpt_a) + enddo + jpt_h = jpt_h + yhalo + yhalo + endif + enddo +! + j1 = latg + (yhalo+yhalo) * nodes_lats + do jj=1,yhalo + j2 = yhalo - jj + global_lats_h(jj) = global_lats_a(1) + j2 ! set north pole yhalo + global_lats_h(j1-j2) = global_lats_a(latg) + 1 - jj ! set south pole yhalo + enddo +! + if (lats_nodes_a(1) /= latg) then +! +! set non-polar south yhalos + jpt_h = 0 + do nn=1,nodes-1 + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h-yhalo) + do jj=1,yhalo + global_lats_h(jpt_h-yhalo+jj) = min(lat_val+jj,latg) + enddo + enddo +! +! set non-polar north yhalos + jpt_h = 0 + do nn=1,nodes-1 + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h+yhalo+1) + do jj=1,yhalo + global_lats_h(jpt_h+yhalo-(jj-1)) = max(lat_val-jj,1) + enddo + enddo +! + endif +! + + iprint = 0 +! iprint = 1 + if (iprint == 1 .and. me == 0) then +! + write(me+6000,'("setlats_h yhalo=",i3," nodes=",i3/)') + & yhalo,nodes +! + do nn=1,nodes + write(me+6000,'("lats_nodes_a(",i4,")=",i4," ", + & " lats_nodes_h(",i4,")=",i4)') + & nn, lats_nodes_a(nn), + & nn, lats_nodes_h(nn) + enddo +! + jpt_a = 0 + do nn=1,nodes + if (lats_nodes_a(nn) > 0) then + write(me+6000,'(" ")') + do jj=1,lats_nodes_a(nn) + jpt_a=jpt_a+1 + write(me+6000,'(2i4," global_lats_a(",i4,")=",i4)') + & nn, jj, jpt_a, global_lats_a(jpt_a) + enddo + endif + enddo +! + jpt_h=0 + do nn=1,nodes + if (lats_nodes_h(nn).gt.0) then + write(me+6000,'(" ")') + do jj=1,lats_nodes_h(nn) + jpt_h=jpt_h+1 + write(me+6000,'(2i4," global_lats_h(",i4,")=",i4)') + & nn, jj, jpt_h, global_lats_h(jpt_h) + enddo + endif + enddo +! + close(6000+me) + endif +! close(6000+me) +! + return + end + + end module setlats_lag_stochy_mod diff --git a/stochastic_physics/spectral_layout.f b/stochastic_physics/spectral_layout.f new file mode 100644 index 000000000..687e2da82 --- /dev/null +++ b/stochastic_physics/spectral_layout.f @@ -0,0 +1,31 @@ + module spectral_layout_mod + + implicit none +! +! program log: +! 20161011 philip pegion : make stochastic pattern generator standalone +! +! 20180731 dom heinzeller : todo: cleanup nodes, me, ... (defined multiple times in confusing ways in several files) +! + + integer nodes, nodes_comp,nodes_io, + & me,lon_dim_a, + & ls_dim, + & ls_max_node, + & lats_dim_a, + & lats_node_a, + & lats_node_a_max, + & ipt_lats_node_a, + & len_trie_ls, + & len_trio_ls, + & len_trie_ls_max, + & len_trio_ls_max, + & me_l_0, + + & lats_dim_ext, + & lats_node_ext, + & ipt_lats_node_ext +! + INTEGER ,ALLOCATABLE :: lat1s_a(:), lon_dims_a(:),lon_dims_ext(:) + + end module spectral_layout_mod diff --git a/stochastic_physics/stochastic_physics.F90 b/stochastic_physics/stochastic_physics.F90 new file mode 100644 index 000000000..a9338d0bb --- /dev/null +++ b/stochastic_physics/stochastic_physics.F90 @@ -0,0 +1,410 @@ +module stochastic_physics + +use stochy_ccpp, only : is_initialized, & + is_master, & + mpicomm, & + mpirank, & + mpiroot, & + mpisize, & + ompthreads + +implicit none + +private + +public :: stochastic_physics_init, stochastic_physics_run, stochastic_physics_finalize + +contains + +!> \section arg_table_stochastic_physics_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_control_type | | inout | F | +!! | nthreads | omp_threads | number of OpenMP threads available for physics schemes | count | 0 | integer | | in | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +subroutine stochastic_physics_init(Model, nthreads, errmsg, errflg) +use stochy_internal_state_mod +use stochy_data_mod, only : nshum,rpattern_shum,init_stochdata,rpattern_sppt,nsppt,rpattern_skeb,nskeb,gg_lats,gg_lons,& + rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_skeb,vfact_sppt,vfact_shum,skeb_vpts,skeb_vwts,sl +use stochy_resol_def, only : latg,lonf,skeblevs +use stochy_gg_def, only : colrad_a +use stochy_namelist_def +use physcons, only: con_pi +use spectral_layout_mod, only:me +use GFS_typedefs, only: GFS_control_type + +implicit none +type(GFS_control_type), intent(inout) :: Model +integer, intent(in) :: nthreads +character(len=*), intent(out) :: errmsg +integer, intent(out) :: errflg + +integer :: nblks +integer :: iret +real*8 :: PRSI(Model%levs),PRSL(Model%levs),dx +real, allocatable :: skeb_vloc(:) +integer :: k,kflip,latghf,nodes,blk,k2 +character*2::proc + +! Initialize CCPP error handling variables +errmsg = '' +errflg = 0 + +! Set/update shared variables in stochy_ccpp +mpicomm = Model%communicator +mpirank = Model%me +mpiroot = Model%master +mpisize = Model%ntasks +ompthreads = nthreads +is_initialized = .true. + +! ------------------------------------------ + +nblks = size(Model%blksz) + +! replace +rad2deg=180.0/con_pi +INTTYP=0 ! bilinear interpolation +me=Model%me +nodes=Model%ntasks +gis_stochy%me=me +gis_stochy%nodes=nodes +call init_stochdata(Model%levs,Model%dtp,Model%input_nml_file,Model%fn_nml,Model%nlunit,iret) +! check namelist entries for consistency +if (Model%do_sppt.neqv.do_sppt) then + write(errmsg,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_sppt and sppt' + errflg = 1 +else if (Model%do_shum.neqv.do_shum) then + write(errmsg,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_shum and shum' + errflg = 1 +else if (Model%do_skeb.neqv.do_skeb) then + write(errmsg,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_skeb and skeb' + errflg = 1 +else if (Model%do_sfcperts.neqv.do_sfcperts) then ! mg, sfc-perts + write(errmsg,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_sfcperts and pertz0 / pertshc / pertzt / pertlai / pertvegf / pertalb' + errflg = 1 +end if +! update remaining model configuration parameters from namelist +Model%use_zmtnblck=use_zmtnblck +Model%skeb_npass=skeb_npass +Model%nsfcpert=nsfcpert ! mg, sfc-perts +Model%pertz0=pertz0 ! mg, sfc-perts +Model%pertzt=pertzt ! mg, sfc-perts +Model%pertshc=pertshc ! mg, sfc-perts +Model%pertlai=pertlai ! mg, sfc-perts +Model%pertalb=pertalb ! mg, sfc-perts +Model%pertvegf=pertvegf ! mg, sfc-perts +if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return +allocate(sl(Model%levs)) +do k=1,Model%levs + sl(k)= 0.5*(Model%ak(k)/101300.+Model%bk(k)+Model%ak(k+1)/101300.0+Model%bk(k+1)) ! si are now sigmas +! if(is_master())print*,'sl(k)',k,sl(k),Model%ak(k),Model%bk(k) +enddo +if (do_sppt) then + allocate(vfact_sppt(Model%levs)) + do k=1,Model%levs + if (sl(k) .lt. sppt_sigtop1 .and. sl(k) .gt. sppt_sigtop2) then + vfact_sppt(k) = (sl(k)-sppt_sigtop2)/(sppt_sigtop1-sppt_sigtop2) + else if (sl(k) .lt. sppt_sigtop2) then + vfact_sppt(k) = 0.0 + else + vfact_sppt(k) = 1.0 + endif + enddo + if (sppt_sfclimit) then + vfact_sppt(2)=vfact_sppt(3)*0.5 + vfact_sppt(1)=0.0 + endif + if (is_master()) then + do k=1,MOdel%levs + print *,'sppt vert profile',k,sl(k),vfact_sppt(k) + enddo + endif +endif +if (do_skeb) then + !print*,'allocating skeb stuff',skeblevs + allocate(vfact_skeb(Model%levs)) + allocate(skeb_vloc(skeblevs)) ! local + allocate(skeb_vwts(Model%levs,2)) ! save for later + allocate(skeb_vpts(Model%levs,2)) ! save for later + do k=1,Model%levs + if (sl(k) .lt. skeb_sigtop1 .and. sl(k) .gt. skeb_sigtop2) then + vfact_skeb(k) = (sl(k)-skeb_sigtop2)/(skeb_sigtop1-skeb_sigtop2) + else if (sl(k) .lt. skeb_sigtop2) then + vfact_skeb(k) = 0.0 + else + vfact_skeb(k) = 1.0 + endif + if (is_master()) print *,'skeb vert profile',k,sl(k),vfact_skeb(k) + enddo +! calculate vertical interpolation weights + do k=1,skeblevs + skeb_vloc(k)=sl(1)-real(k-1)/real(skeblevs-1.0)*(sl(1)-sl(Model%levs)) + enddo +! surface +skeb_vwts(1,2)=0 +skeb_vpts(1,1)=1 +! top +skeb_vwts(Model%levs,2)=1 +skeb_vpts(Model%levs,1)=skeblevs-2 +! internal +DO k=2,Model%levs-1 + DO k2=1,skeblevs-1 + IF (sl(k) .LE. skeb_vloc(k2) .AND. sl(k) .GT. skeb_vloc(k2+1)) THEN + skeb_vpts(k,1)=k2 + skeb_vwts(k,2)=(skeb_vloc(k2)-sl(k))/(skeb_vloc(k2)-skeb_vloc(k2+1)) + ENDIF + ENDDO +ENDDO +deallocate(skeb_vloc) +if (is_master()) then +DO k=1,Model%levs + print*,'skeb vpts ',skeb_vpts(k,1),skeb_vwts(k,2) +ENDDO +endif +skeb_vwts(:,1)=1.0-skeb_vwts(:,2) +skeb_vpts(:,2)=skeb_vpts(:,1)+1.0 +endif + +if (do_shum) then + allocate(vfact_shum(Model%levs)) + do k=1,Model%levs + vfact_shum(k) = exp((sl(k)-1.)/shum_sigefold) + if (sl(k).LT. 2*shum_sigefold) then + vfact_shum(k)=0.0 + endif + if (is_master()) print *,'shum vert profile',k,sl(k),vfact_shum(k) + enddo +endif +! get interpolation weights +! define gaussian grid lats and lons +latghf=latg/2 +!print *,'define interp weights',latghf,lonf +!print *,allocated(gg_lats),allocated(gg_lons) +allocate(gg_lats(latg)) +!print *,'aloocated lats' +allocate(gg_lons(lonf)) +!print *,'aloocated lons' +do k=1,latghf + gg_lats(k)=-1.0*colrad_a(latghf-k+1)*rad2deg + gg_lats(latg-k+1)=-1*gg_lats(k) +enddo +dx=360.0/lonf +!print*,'dx=',dx +do k=1,lonf + gg_lons(k)=dx*(k-1) +enddo +WLON=gg_lons(1)-(gg_lons(2)-gg_lons(1)) +RNLAT=gg_lats(1)*2-gg_lats(2) + +!print *,'done with init_stochastic_physics' + +if(Model%me == Model%master) print*,'do_skeb=',Model%do_skeb + +end subroutine stochastic_physics_init + + +!> \section arg_table_stochastic_physics_run Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_control_type | | in | F | +!! | Data | FV3-GFS_Data_type_all_blocks | Fortran DDT containing FV3-GFS data | DDT | 1 | GFS_data_type | | inout | F | +!! | nthreads | omp_threads | number of OpenMP threads available for physics schemes | count | 0 | integer | | in | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +subroutine stochastic_physics_run(Model, Data, nthreads, errmsg, errflg) +use stochy_internal_state_mod +use stochy_data_mod, only : nshum,rpattern_shum,rpattern_sppt,nsppt,rpattern_skeb,nskeb,& + rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_sppt,vfact_shum,vfact_skeb +use get_stochy_pattern_mod,only : get_random_pattern_fv3,get_random_pattern_fv3_vect,dump_patterns +use stochy_resol_def , only : latg,lonf +use stochy_namelist_def +use spectral_layout_mod,only:me +use GFS_typedefs, only: GFS_control_type, GFS_data_type +implicit none +type(GFS_control_type), intent(in) :: Model +type(GFS_data_type), intent(inout) :: Data(:) +integer, intent(in) :: nthreads +character(len=*), intent(out) :: errmsg +integer, intent(out) :: errflg + +real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:) +!D-grid +integer :: k +integer j,ierr,i +integer :: nblks, blk, len, maxlen +character*120 :: sfile +character*6 :: STRFH + +! Initialize CCPP error handling variables +errmsg = '' +errflg = 0 + +if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return + +! Update number of threads in shared variables in stochy_ccpp and set block-related variables +ompthreads = nthreads +nblks = size(Model%blksz) +maxlen = maxval(Model%blksz(:)) + +! check to see if it is time to write out random patterns +if (Model%phour .EQ. fhstoch) then + write(STRFH,FMT='(I6.6)') nint(Model%phour) + sfile='stoch_out.F'//trim(STRFH) + call dump_patterns(sfile) +endif +allocate(tmp_wts(nblks,maxlen)) +allocate(tmpu_wts(nblks,maxlen,Model%levs)) +allocate(tmpv_wts(nblks,maxlen,Model%levs)) +if (do_sppt) then + call get_random_pattern_fv3(rpattern_sppt,nsppt,gis_stochy,Model,Data(:)%Grid,nblks,maxlen,tmp_wts) + DO blk=1,nblks + len=size(Data(blk)%Grid%xlat,1) + DO k=1,Model%levs + Data(blk)%Coupling%sppt_wts(:,k)=tmp_wts(blk,1:len)*vfact_sppt(k) + ENDDO + if (sppt_logit) Data(blk)%Coupling%sppt_wts(:,:) = (2./(1.+exp(Data(blk)%Coupling%sppt_wts(:,:))))-1. + Data(blk)%Coupling%sppt_wts(:,:)= Data(blk)%Coupling%sppt_wts(:,:)+1.0 + ENDDO +endif +if (do_shum) then + call get_random_pattern_fv3(rpattern_shum,nshum,gis_stochy,Model,Data(:)%Grid,nblks,maxlen,tmp_wts) + DO blk=1,nblks + len=size(Data(blk)%Grid%xlat,1) + DO k=1,Model%levs + Data(blk)%Coupling%shum_wts(:,k)=tmp_wts(blk,1:len)*vfact_shum(k) + ENDDO + ENDDO +endif +if (do_skeb) then + call get_random_pattern_fv3_vect(rpattern_skeb,nskeb,gis_stochy,Model,Data(:)%Grid,nblks,maxlen,tmpu_wts,tmpv_wts) + DO blk=1,nblks + len=size(Data(blk)%Grid%xlat,1) + DO k=1,Model%levs + Data(blk)%Coupling%skebu_wts(:,k)=tmpu_wts(blk,1:len,k)*vfact_skeb(k) + Data(blk)%Coupling%skebv_wts(:,k)=tmpv_wts(blk,1:len,k)*vfact_skeb(k) + ENDDO + ENDDO +endif +deallocate(tmp_wts) +deallocate(tmpu_wts) +deallocate(tmpv_wts) + +end subroutine stochastic_physics_run + +!> \section arg_table_stochastic_physics_finalize Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +subroutine stochastic_physics_finalize(errmsg, errflg) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not.is_initialized) return + + is_initialized = .false. + +end subroutine stochastic_physics_finalize + +end module stochastic_physics + + +module stochastic_physics_sfc + +use stochy_ccpp, only : is_initialized, is_master + +implicit none + +private +public :: stochastic_physics_sfc_init, stochastic_physics_sfc_run, stochastic_physics_sfc_finalize + +contains + +!> \section arg_table_stochastic_physics_sfc_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_control_type | | in | F | +!! | Data | FV3-GFS_Data_type_all_blocks | Fortran DDT containing FV3-GFS data | DDT | 1 | GFS_data_type | | inout | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +subroutine stochastic_physics_sfc_init(Model, Data, errmsg, errflg) +use stochy_internal_state_mod +use stochy_data_mod, only : rad2deg,INTTYP,wlon,rnlat,gis_stochy, rpattern_sfc,npsfc ! mg, sfc-perts +use get_stochy_pattern_mod,only : get_random_pattern_sfc_fv3 ! mg, sfc-perts +use stochy_resol_def , only : latg,lonf +use stochy_namelist_def +use spectral_layout_mod,only:me +use GFS_typedefs, only: GFS_control_type, GFS_data_type +implicit none +type(GFS_control_type), intent(in) :: Model +type(GFS_data_type), intent(inout) :: Data(:) +character(len=*), intent(out) :: errmsg +integer, intent(out) :: errflg + +real,allocatable :: tmpsfc_wts(:,:,:) +!D-grid +integer :: k +integer j,ierr,i +integer :: nblks, blk, len, maxlen +character*120 :: sfile +character*6 :: STRFH + +! Initialize CCPP error handling variables +errmsg = '' +errflg = 0 + +if (.NOT. do_sfcperts) return + +! stochastic_physics_sfc_init depends on stochastic_physics_init being run first; +! in general, stochastic_physics_sfc can only be run with/after stochastic_physics +! check initialization status in stochy_ccpp to make sure this is true +if (.not.is_initialized) then + write(errmsg,'(*(a))') 'Logic error: stochastic_physics_init must be called before stochastic_physics_sfc_init' + errflg = 1 + return +end if + +! Set block-related variables +nblks = size(Model%blksz) +maxlen = maxval(Model%blksz(:)) + +allocate(tmpsfc_wts(nblks,maxlen,Model%nsfcpert)) ! mg, sfc-perts +if (is_master()) then + print*,'In stochastic_physics_sfc_init: do_sfcperts ',do_sfcperts +endif +call get_random_pattern_sfc_fv3(rpattern_sfc,npsfc,gis_stochy,Model,Data(:)%Grid,nblks,maxlen,tmpsfc_wts) +DO blk=1,nblks + len=size(Data(blk)%Grid%xlat,1) + DO k=1,Model%nsfcpert + Data(blk)%Coupling%sfc_wts(:,k)=tmpsfc_wts(blk,1:len,k) + ENDDO +ENDDO +if (is_master()) then + print*,'tmpsfc_wts(blk,1,:) =',tmpsfc_wts(1,1,1),tmpsfc_wts(1,1,2),tmpsfc_wts(1,1,3),tmpsfc_wts(1,1,4),tmpsfc_wts(1,1,5) + print*,'min(tmpsfc_wts(:,:,:)) =',minval(tmpsfc_wts(:,:,:)) +endif +deallocate(tmpsfc_wts) +end subroutine stochastic_physics_sfc_init + +subroutine stochastic_physics_sfc_run() +end subroutine stochastic_physics_sfc_run + +subroutine stochastic_physics_sfc_finalize() +end subroutine stochastic_physics_sfc_finalize + +end module stochastic_physics_sfc diff --git a/stochastic_physics/stochy_ccpp.F90 b/stochastic_physics/stochy_ccpp.F90 new file mode 100644 index 000000000..b6a84c49f --- /dev/null +++ b/stochastic_physics/stochy_ccpp.F90 @@ -0,0 +1,438 @@ +module stochy_ccpp + +#ifdef MPI + use mpi +#endif + + implicit none + + private + + public is_initialized + public is_master + public stochy_la2ga + public mp_bcst + public mp_reduce_sum + public mpicomm + public mpirank + public mpiroot + public mpisize + public mpp_alltoall + public ompthreads + + logical :: is_initialized = .false. + + interface mpp_alltoall + !module procedure mpp_alltoall_int4 + !module procedure mpp_alltoall_int8 + !module procedure mpp_alltoall_real4 + !module procedure mpp_alltoall_real8 + !module procedure mpp_alltoall_int4_v + !module procedure mpp_alltoall_int8_v + module procedure mpp_alltoall_real4_v + !module procedure mpp_alltoall_real8_v + end interface + + !> The interface 'mp_bcast contains routines that call SPMD broadcast + !! (one-to-many communication). + interface mp_bcst + module procedure mp_bcst_i + !module procedure mp_bcst_r4 + !module procedure mp_bcst_r8 + !module procedure mp_bcst_1d_r4 + module procedure mp_bcst_1d_r8 + !module procedure mp_bcst_2d_r4 + !module procedure mp_bcst_2d_r8 + !module procedure mp_bcst_3d_r4 + !module procedure mp_bcst_3d_r8 + !module procedure mp_bcst_4d_r4 + !module procedure mp_bcst_4d_r8 + module procedure mp_bcst_1d_i + !module procedure mp_bcst_2d_i + !module procedure mp_bcst_3d_i + !module procedure mp_bcst_4d_i + end interface + + !> The interface 'mp_reduce_sum' contains routines that call SPMD_REDUCE. + !! The routines compute the sums of values and place the net sum in a result. + interface mp_reduce_sum + !module procedure mp_reduce_sum_r4 + !module procedure mp_reduce_sum_r4_1d + !module procedure mp_reduce_sum_r4_1darr + !module procedure mp_reduce_sum_r4_2darr + !module procedure mp_reduce_sum_r8 + !module procedure mp_reduce_sum_r8_1d + module procedure mp_reduce_sum_r8_1darr + module procedure mp_reduce_sum_r8_2darr + end interface + +#ifdef MPI + integer, save, target :: mpicomm = MPI_COMM_WORLD +#else + integer, save, target :: mpicomm = 0 +#endif + integer, save, target :: mpirank = 0 + integer, save, target :: mpiroot = 0 + integer, save, target :: mpisize = 1 + + integer, save, target :: ompthreads = 1 + +contains + + + function is_master() result(match) + + logical :: match + + match = (mpirank==mpiroot) + + end function is_master + + + subroutine mpp_alltoall_real4_v(sbuf, ssize, sdispl, rbuf, rsize, rdispl) + + real(kind=4), intent(in) :: sbuf(:) + real(kind=4), intent(inout) :: rbuf(:) + + integer :: ierr + + integer, intent(in) :: ssize(:), rsize(:) + integer, intent(in) :: sdispl(:), rdispl(:) + +#ifdef MPI + call MPI_Alltoallv( sbuf, ssize, sdispl, MPI_REAL, & + rbuf, rsize, rdispl, MPI_REAL, & + mpicomm, ierr ) +#else + rbuf = sbuf +#endif + + end subroutine mpp_alltoall_real4_v + + + subroutine mp_bcst_i(q) + integer, intent(inout) :: q + integer :: ierr + +#ifdef MPI + call MPI_BCAST(q, 1, MPI_INTEGER, mpiroot, mpicomm, ierr) +#endif + + end subroutine mp_bcst_i + + + subroutine mp_bcst_1d_i(q, idim) + + integer, intent(in) :: idim + integer, intent(inout) :: q(idim) + + integer :: ierr + +#ifdef MPI + call MPI_BCAST(q, idim, MPI_INTEGER, mpiroot, mpicomm, ierr) +#endif + + end subroutine mp_bcst_1d_i + + + subroutine mp_bcst_1d_r8(q, idim) + + integer, intent(in) :: idim + real(kind=8), intent(inout) :: q(idim) + + integer :: ierr + +#ifdef MPI + call MPI_BCAST(q, idim, MPI_DOUBLE_PRECISION, mpiroot, mpicomm, ierr) +#endif + + end subroutine mp_bcst_1d_r8 + + + subroutine mp_reduce_sum_r8_1darr(mysum, npts) + + integer, intent(in) :: npts + real(kind=8), intent(inout) :: mysum(npts) + + real(kind=8) :: gsum(npts) + integer :: ierr + +#ifdef MPI + gsum = 0.0 + + call MPI_ALLREDUCE( mysum, gsum, npts, & + MPI_DOUBLE_PRECISION, MPI_SUM, & + mpicomm, ierr ) + + mysum = gsum +#endif + + end subroutine mp_reduce_sum_r8_1darr + + + subroutine mp_reduce_sum_r8_2darr(mysum,npts1,npts2) + + integer, intent(in) :: npts1,npts2 + real(kind=8), intent(inout) :: mysum(npts1,npts2) + + real(kind=8) :: gsum(npts1,npts2) + integer :: ierr + +#ifdef MPI + gsum = 0.0 + + call MPI_ALLREDUCE( mysum, gsum, npts1*npts2, & + MPI_DOUBLE_PRECISION, MPI_SUM, & + mpicomm, ierr ) + + mysum = gsum +#endif + + end subroutine mp_reduce_sum_r8_2darr + + + ! + ! interpolation from lat/lon or gaussian grid to other lat/lon grid + ! + subroutine stochy_la2ga(regin,imxin,jmxin,rinlon,rinlat,rlon,rlat, & + gauout,len,rslmsk, outlat, outlon) + use machine , only : kind_io8, kind_io4 + implicit none + ! interface variables + real (kind=kind_io8), intent(in) :: regin(imxin,jmxin) + integer, intent(in) :: imxin + integer, intent(in) :: jmxin + real (kind=kind_io8), intent(in) :: rinlon(imxin) + real (kind=kind_io8), intent(in) :: rinlat(jmxin) + real (kind=kind_io8), intent(in) :: rlon + real (kind=kind_io8), intent(in) :: rlat + real (kind=kind_io8), intent(out) :: gauout(len) + integer, intent(in) :: len + real (kind=kind_io8), intent(in) :: rslmsk(imxin,jmxin) + real (kind=kind_io8), intent(in) :: outlat(len) + real (kind=kind_io8), intent(in) :: outlon(len) + ! local variables + real (kind=kind_io8) :: wei4,wei3,wei2,sum2,sum1,sum3,wei1,sum4 + real (kind=kind_io8) :: wsum,wsumiv,sums,sumn,wi2j2,x,y,wi1j1 + real (kind=kind_io8) :: wi1j2,wi2j1,aphi,rnume,alamd,denom + integer :: jy,ix,i,j,jq,jx + integer :: j1,j2,ii,i1,i2,kmami,it + integer :: nx,kxs,kxt + integer :: iindx1(len) + integer :: iindx2(len) + integer :: jindx1(len) + integer :: jindx2(len) + real(kind=kind_io8) :: ddx(len) + real(kind=kind_io8) :: ddy(len) + real(kind=kind_io8) :: wrk(len) + integer :: len_thread_m + integer :: len_thread + integer :: i1_t + integer :: i2_t + integer, pointer :: num_threads => ompthreads +! + len_thread_m = (len+num_threads-1) / num_threads +! + !$omp parallel do default(none) & + !$omp private(i1_t,i2_t,len_thread,it,i,ii,i1,i2) & + !$omp private(j,j1,j2,jq,ix,jy,nx,kxs,kxt,kmami) & + !$omp private(alamd,denom,rnume,aphi,x,y,wsum,wsumiv,sum1,sum2) & + !$omp private(sum3,sum4,wi1j1,wi2j1,wi1j2,wi2j2,wei1,wei2,wei3,wei4) & + !$omp private(sumn,sums) & + !$omp shared(imxin,jmxin) & + !$omp shared(outlon,outlat,wrk,iindx1,rinlon,jindx1,rinlat,ddx,ddy) & + !$omp shared(rlon,rlat,regin,gauout) & + !$omp shared(num_threads,len_thread_m,len,iindx2,jindx2,rslmsk) + do it=1,num_threads ! start of threaded loop + i1_t = (it-1)*len_thread_m+1 + i2_t = min(i1_t+len_thread_m-1,len) + len_thread = i2_t-i1_t+1 +! +! find i-index for interpolation +! + do i=i1_t, i2_t + alamd = outlon(i) + if (alamd .lt. rlon) alamd = alamd + 360.0 + if (alamd .gt. 360.0+rlon) alamd = alamd - 360.0 + wrk(i) = alamd + iindx1(i) = imxin + enddo + do i=i1_t,i2_t + do ii=1,imxin + if(wrk(i) .ge. rinlon(ii)) iindx1(i) = ii + enddo + enddo + do i=i1_t,i2_t + i1 = iindx1(i) + if (i1 .lt. 1) i1 = imxin + i2 = i1 + 1 + if (i2 .gt. imxin) i2 = 1 + iindx1(i) = i1 + iindx2(i) = i2 + denom = rinlon(i2) - rinlon(i1) + if(denom.lt.0.) denom = denom + 360. + rnume = wrk(i) - rinlon(i1) + if(rnume.lt.0.) rnume = rnume + 360. + ddx(i) = rnume / denom + enddo +! +! find j-index for interplation +! + if(rlat.gt.0.) then + do j=i1_t,i2_t + jindx1(j)=0 + enddo + do jx=1,jmxin + do j=i1_t,i2_t + if(outlat(j).le.rinlat(jx)) jindx1(j) = jx + enddo + enddo + do j=i1_t,i2_t + jq = jindx1(j) + aphi=outlat(j) + if(jq.ge.1 .and. jq .lt. jmxin) then + j2=jq+1 + j1=jq + ddy(j)=(aphi-rinlat(j1))/(rinlat(j2)-rinlat(j1)) + elseif (jq .eq. 0) then + j2=1 + j1=1 + if(abs(90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + else + j2=jmxin + j1=jmxin + if(abs(-90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(-90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + endif + jindx1(j)=j1 + jindx2(j)=j2 + enddo + else + do j=i1_t,i2_t + jindx1(j) = jmxin+1 + enddo + do jx=jmxin,1,-1 + do j=i1_t,i2_t + if(outlat(j).le.rinlat(jx)) jindx1(j) = jx + enddo + enddo + do j=i1_t,i2_t + jq = jindx1(j) + aphi=outlat(j) + if(jq.gt.1 .and. jq .le. jmxin) then + j2=jq + j1=jq-1 + ddy(j)=(aphi-rinlat(j1))/(rinlat(j2)-rinlat(j1)) + elseif (jq .eq. 1) then + j2=1 + j1=1 + if(abs(-90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(-90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + else + j2=jmxin + j1=jmxin + if(abs(90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + endif + jindx1(j)=j1 + jindx2(j)=j2 + enddo + endif +! + sum1 = 0. + sum2 = 0. + sum3 = 0. + sum4 = 0. + do i=1,imxin + sum1 = sum1 + regin(i,1) + sum2 = sum2 + regin(i,jmxin) + enddo + sum1 = sum1 / imxin + sum2 = sum2 / imxin + sum3 = sum1 + sum4 = sum2 +! +! quasi-bilinear interpolation +! + do i=i1_t,i2_t + y = ddy(i) + j1 = jindx1(i) + j2 = jindx2(i) + x = ddx(i) + i1 = iindx1(i) + i2 = iindx2(i) +! + wi1j1 = (1.-x) * (1.-y) + wi2j1 = x *( 1.-y) + wi1j2 = (1.-x) * y + wi2j2 = x * y +! + wsum = wi1j1 + wi2j1 + wi1j2 + wi2j2 + wrk(i) = wsum + if(wsum.ne.0.) then + wsumiv = 1./wsum + if(j1.ne.j2) then + gauout(i) = (wi1j1*regin(i1,j1) + wi2j1*regin(i2,j1) + & + wi1j2*regin(i1,j2) + wi2j2*regin(i2,j2)) & + *wsumiv + else + if (rlat .gt. 0.0) then + sumn = sum3 + sums = sum4 + if( j1 .eq. 1) then + gauout(i) = (wi1j1*sumn +wi2j1*sumn + & + wi1j2*regin(i1,j2)+wi2j2*regin(i2,j2)) & + * wsumiv + elseif (j1 .eq. jmxin) then + gauout(i) = (wi1j1*regin(i1,j1)+wi2j1*regin(i2,j1)+ & + wi1j2*sums +wi2j2*sums ) & + * wsumiv + endif + else + sums = sum3 + sumn = sum4 + if( j1 .eq. 1) then + gauout(i) = (wi1j1*regin(i1,j1)+wi2j1*regin(i2,j1)+ & + wi1j2*sums +wi2j2*sums ) & + * wsumiv + elseif (j1 .eq. jmxin) then + gauout(i) = (wi1j1*sumn +wi2j1*sumn + & + wi1j2*regin(i1,j2)+wi2j2*regin(i2,j2)) & + * wsumiv + endif + endif + endif ! if j1 .ne. j2 + endif + enddo + do i=i1_t,i2_t + j1 = jindx1(i) + j2 = jindx2(i) + i1 = iindx1(i) + i2 = iindx2(i) + if(wrk(i) .eq. 0.0) then + write(6,*) ' la2ga: bad rslmsk given' + call sleep(2) + stop + endif + enddo + enddo ! end of threaded loop +!$omp end parallel do +! + return +! + end subroutine stochy_la2ga + +end module stochy_ccpp diff --git a/stochastic_physics/stochy_data_mod.F90 b/stochastic_physics/stochy_data_mod.F90 new file mode 100644 index 000000000..80b598f8b --- /dev/null +++ b/stochastic_physics/stochy_data_mod.F90 @@ -0,0 +1,393 @@ +module stochy_data_mod + +! set up and initialize stochastic random patterns. + + use spectral_layout_mod, only: len_trie_ls,len_trio_ls,ls_dim,ls_max_node + use stochy_resol_def, only : skeblevs,levs,jcap,lonf,latg + use stochy_namelist_def + use physcons, only : radius => con_rerth + use stochy_ccpp, only: is_master, & + mp_bcst, & + me => mpirank, & + nodes => mpisize + use stochy_patterngenerator_mod, only: random_pattern, patterngenerator_init,& + getnoise, patterngenerator_advance,ndimspec,chgres_pattern,computevarspec_r + use initialize_spectral_mod, only: initialize_spectral + use stochy_internal_state_mod +! use mersenne_twister_stochy, only : random_seed + use mersenne_twister, only : random_seed + use compns_stochy_mod, only : compns_stochy + + implicit none + private + public :: init_stochdata + + type(random_pattern), public, save, allocatable, dimension(:) :: & + rpattern_sppt,rpattern_shum,rpattern_skeb, rpattern_sfc + integer, public :: nsppt=0 + integer, public :: nshum=0 + integer, public :: nskeb=0 + integer, public :: npsfc=0 + real*8, public,allocatable :: sl(:) + + real(kind=kind_dbl_prec),public, allocatable :: vfact_sppt(:),vfact_shum(:),vfact_skeb(:) + real(kind=kind_dbl_prec),public, allocatable :: skeb_vwts(:,:),skeb_vpts(:,:) + real(kind=kind_dbl_prec),public, allocatable :: gg_lats(:),gg_lons(:) + real(kind=kind_dbl_prec),public :: wlon,rnlat,rad2deg + real(kind=kind_dbl_prec),public, allocatable :: skebu_save(:,:,:),skebv_save(:,:,:) + integer,public :: INTTYP + type(stochy_internal_state),public :: gis_stochy + + contains + subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) + +! initialize random patterns. A spinup period of spinup_efolds times the +! temporal time scale is run for each pattern. + integer, intent(in) :: nlunit,nlevs + character(len=*), intent(in) :: input_nml_file(:) + character(len=64), intent(in) :: fn_nml + real, intent(in) :: delt + integer, intent(out) :: iret + real :: pertsfc(1) + + real :: rnn1 + integer :: nn,nspinup,k,nm,spinup_efolds,stochlun,ierr,n + integer :: locl,indev,indod,indlsod,indlsev + integer :: l,jbasev,jbasod + real(kind_dbl_prec),allocatable :: noise_e(:,:),noise_o(:,:) + include 'function_indlsod' + include 'function_indlsev' + stochlun=99 + levs=nlevs + + iret=0 + if(is_master()) print*,'in init stochdata' + call compns_stochy (me,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) + if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return + if (nodes.GE.lat_s/2) then + lat_s=(int(nodes/12)+1)*24 + lon_s=lat_s*2 + ntrunc=lat_s-2 + if (is_master()) print*,'WARNING: spectral resolution is too low for number of mpi_tasks, resetting lon_s,lat_s,and ntrunc to',lon_s,lat_s,ntrunc + endif + call initialize_spectral(gis_stochy, iret) + if (iret/=0) return + allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) +! determine number of random patterns to be used for each scheme. + do n=1,size(sppt) + if (sppt(n) > 0) then + nsppt=nsppt+1 + else + exit + endif + enddo + if (is_master()) print *,'nsppt = ',nsppt + do n=1,size(shum) + if (shum(n) > 0) then + nshum=nshum+1 + else + exit + endif + enddo + if (is_master()) print *,'nshum = ',nshum + do n=1,size(skeb) + if (skeb(n) > 0) then + nskeb=nskeb+1 + else + exit + endif + enddo + if (is_master()) print *,'nskeb = ',nskeb + ! mg, sfc-perts + do n=1,size(pertz0) + if (pertz0(n) > 0 .or. pertzt(n)>0 .or. pertshc(n)>0 .or. & + pertvegf(n)>0 .or. pertlai(n)>0 .or. pertalb(n)>0) then + npsfc=npsfc+1 + else + exit + endif + enddo + if (is_master()) then + if (npsfc > 0) then + print *,' npsfc = ', npsfc + print *,' pertz0 = ', pertz0 + print *,' pertzt = ', pertzt + print *,' pertshc = ', pertshc + print *,' pertlai = ', pertlai + print *,' pertalb = ', pertalb + print *,' pertvegf = ', pertvegf + endif + endif + + if (nsppt > 0) allocate(rpattern_sppt(nsppt)) + if (nshum > 0) allocate(rpattern_shum(nshum)) + if (nskeb > 0) allocate(rpattern_skeb(nskeb)) + ! mg, sfc perts + if (npsfc > 0) allocate(rpattern_sfc(npsfc)) + +! if stochini is true, then read in pattern from a file + if (is_master()) then + if (stochini) then + print*,'opening stoch_ini' + OPEN(stochlun,file='stoch_ini',form='unformatted',iostat=ierr,status='old') + if (ierr .NE. 0) then + write(0,*) 'error opening stoch_ini' + iret = ierr + return + end if + endif + endif + ! no spinup needed if initial patterns are defined correctly. + spinup_efolds = 0 + if (nsppt > 0) then + if (is_master()) print *, 'Initialize random pattern for SPPT' + call patterngenerator_init(sppt_lscale,delt,sppt_tau,sppt,iseed_sppt,rpattern_sppt, & + lonf,latg,jcap,gis_stochy%ls_node,nsppt,1,0) + do n=1,nsppt + nspinup = spinup_efolds*sppt_tau(n)/delt + if (stochini) then + call read_pattern(rpattern_sppt(n),1,stochlun) + else + call getnoise(rpattern_sppt(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_sppt(n)%spec_e(nn,1,1)=noise_e(nn,1) + rpattern_sppt(n)%spec_e(nn,2,1)=noise_e(nn,2) + nm = rpattern_sppt(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_sppt(n)%spec_e(nn,1,1) = rpattern_sppt(n)%stdev*rpattern_sppt(n)%spec_e(nn,1,1)*rpattern_sppt(n)%varspectrum(nm) + rpattern_sppt(n)%spec_e(nn,2,1) = rpattern_sppt(n)%stdev*rpattern_sppt(n)%spec_e(nn,2,1)*rpattern_sppt(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_sppt(n)%spec_o(nn,1,1)=noise_o(nn,1) + rpattern_sppt(n)%spec_o(nn,2,1)=noise_o(nn,2) + nm = rpattern_sppt(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_sppt(n)%spec_o(nn,1,1) = rpattern_sppt(n)%stdev*rpattern_sppt(n)%spec_o(nn,1,1)*rpattern_sppt(n)%varspectrum(nm) + rpattern_sppt(n)%spec_o(nn,2,1) = rpattern_sppt(n)%stdev*rpattern_sppt(n)%spec_o(nn,2,1)*rpattern_sppt(n)%varspectrum(nm) + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_sppt(n),1,.false.) + enddo + endif + enddo + endif + if (nshum > 0) then + if (is_master()) print *, 'Initialize random pattern for SHUM' + call patterngenerator_init(shum_lscale,delt,shum_tau,shum,iseed_shum,rpattern_shum, & + lonf,latg,jcap,gis_stochy%ls_node,nshum,1,0) + do n=1,nshum + nspinup = spinup_efolds*shum_tau(n)/delt + if (stochini) then + call read_pattern(rpattern_shum(n),1,stochlun) + else + call getnoise(rpattern_shum(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_shum(n)%spec_e(nn,1,1)=noise_e(nn,1) + rpattern_shum(n)%spec_e(nn,2,1)=noise_e(nn,2) + nm = rpattern_shum(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_shum(n)%spec_e(nn,1,1) = rpattern_shum(n)%stdev*rpattern_shum(n)%spec_e(nn,1,1)*rpattern_shum(n)%varspectrum(nm) + rpattern_shum(n)%spec_e(nn,2,1) = rpattern_shum(n)%stdev*rpattern_shum(n)%spec_e(nn,2,1)*rpattern_shum(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_shum(n)%spec_o(nn,1,1)=noise_o(nn,1) + rpattern_shum(n)%spec_o(nn,2,1)=noise_o(nn,2) + nm = rpattern_shum(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_shum(n)%spec_o(nn,1,1) = rpattern_shum(n)%stdev*rpattern_shum(n)%spec_o(nn,1,1)*rpattern_shum(n)%varspectrum(nm) + rpattern_shum(n)%spec_o(nn,2,1) = rpattern_shum(n)%stdev*rpattern_shum(n)%spec_o(nn,2,1)*rpattern_shum(n)%varspectrum(nm) + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_shum(n),1,.false.) + enddo + endif + enddo + endif + + if (nskeb > 0) then + ! determine number of skeb levels to deal with temperoal/vertical correlations + skeblevs=nint(skeb_tau(1)/delt*skeb_vdof) +! backscatter noise. + if (is_master()) print *, 'Initialize random pattern for SKEB',skeblevs + call patterngenerator_init(skeb_lscale,delt,skeb_tau,skeb,iseed_skeb,rpattern_skeb, & + lonf,latg,jcap,gis_stochy%ls_node,nskeb,skeblevs,skeb_varspect_opt) + do n=1,nskeb + do k=1,skeblevs + nspinup = spinup_efolds*skeb_tau(n)/delt + if (stochini) then + call read_pattern(rpattern_skeb(n),k,stochlun) + if (is_master()) print *, 'skeb read',k,rpattern_skeb(n)%spec_o(5,1,k) + else + call getnoise(rpattern_skeb(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_skeb(n)%spec_e(nn,1,k)=noise_e(nn,1) + rpattern_skeb(n)%spec_e(nn,2,k)=noise_e(nn,2) + nm = rpattern_skeb(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_skeb(n)%spec_e(nn,1,k) = rpattern_skeb(n)%stdev*rpattern_skeb(n)%spec_e(nn,1,k)*rpattern_skeb(n)%varspectrum(nm) + rpattern_skeb(n)%spec_e(nn,2,k) = rpattern_skeb(n)%stdev*rpattern_skeb(n)%spec_e(nn,2,k)*rpattern_skeb(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_skeb(n)%spec_o(nn,1,k)=noise_o(nn,1) + rpattern_skeb(n)%spec_o(nn,2,k)=noise_o(nn,2) + nm = rpattern_skeb(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_skeb(n)%spec_o(nn,1,k) = rpattern_skeb(n)%stdev*rpattern_skeb(n)%spec_o(nn,1,k)*rpattern_skeb(n)%varspectrum(nm) + rpattern_skeb(n)%spec_o(nn,2,k) = rpattern_skeb(n)%stdev*rpattern_skeb(n)%spec_o(nn,2,k)*rpattern_skeb(n)%varspectrum(nm) + enddo + endif + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_skeb(n),skeblevs,.false.) + enddo + enddo + + gis_stochy%kenorm_e=1. + gis_stochy%kenorm_o=1. ! used to convert forcing pattern to wind field. +if (skebnorm==0) then + do locl=1,ls_max_node + l = gis_stochy%ls_node(locl) + jbasev = gis_stochy%ls_node(locl+ls_dim) + indev = indlsev(l,l) + jbasod = gis_stochy%ls_node(locl+2*ls_dim) + indod = indlsod(l+1,l) + do n=l,jcap,2 + rnn1 = n*(n+1.) + gis_stochy%kenorm_e(indev) = rnn1/radius**2 + indev = indev + 1 + enddo + do n=l+1,jcap,2 + rnn1 = n*(n+1.) + gis_stochy%kenorm_o(indod) = rnn1/radius**2 + indod = indod + 1 + enddo + enddo + if (is_master()) print*,'using streamfunction ',maxval(gis_stochy%kenorm_e(:)),minval(gis_stochy%kenorm_e(:)) +endif +if (skebnorm==1) then + do locl=1,ls_max_node + l = gis_stochy%ls_node(locl) + jbasev = gis_stochy%ls_node(locl+ls_dim) + indev = indlsev(l,l) + jbasod = gis_stochy%ls_node(locl+2*ls_dim) + indod = indlsod(l+1,l) + do n=l,jcap,2 + rnn1 = n*(n+1.) + gis_stochy%kenorm_e(indev) = sqrt(rnn1)/radius + indev = indev + 1 + enddo + do n=l+1,jcap,2 + rnn1 = n*(n+1.) + gis_stochy%kenorm_o(indod) = sqrt(rnn1)/radius + indod = indod + 1 + enddo + enddo + if (is_master()) print*,'using kenorm ',maxval(gis_stochy%kenorm_e(:)),minval(gis_stochy%kenorm_e(:)) +endif + ! set the even and odd (n-l) terms of the top row to zero +do locl=1,ls_max_node + l = gis_stochy%ls_node(locl) + jbasev = gis_stochy%ls_node(locl+ls_dim) + jbasod = gis_stochy%ls_node(locl+2*ls_dim) + if (mod(l,2) .eq. mod(jcap+1,2)) then + gis_stochy%kenorm_e(indlsev(jcap+1,l)) = 0. + endif + if (mod(l,2) .ne. mod(jcap+1,2)) then + gis_stochy%kenorm_o(indlsod(jcap+1,l)) = 0. + endif +enddo + + endif ! skeb > 0 +! mg, sfc-perts +if (npsfc > 0) then + pertsfc(1) = 1. + call patterngenerator_init(sfc_lscale,delt,sfc_tau,pertsfc,iseed_sfc,rpattern_sfc, & + lonf,latg,jcap,gis_stochy%ls_node,npsfc,nsfcpert,0) + do n=1,npsfc + if (is_master()) print *, 'Initialize random pattern for SFC-PERTS',n + do k=1,nsfcpert + nspinup = spinup_efolds*sfc_tau(n)/delt + call getnoise(rpattern_sfc(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_sfc(n)%spec_e(nn,1,k)=noise_e(nn,1) + rpattern_sfc(n)%spec_e(nn,2,k)=noise_e(nn,2) + nm = rpattern_sfc(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_sfc(n)%spec_e(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) + rpattern_sfc(n)%spec_e(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_sfc(n)%spec_o(nn,1,k)=noise_o(nn,1) + rpattern_sfc(n)%spec_o(nn,2,k)=noise_o(nn,2) + nm = rpattern_sfc(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_sfc(n)%spec_o(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) + rpattern_sfc(n)%spec_o(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_sfc(n),k,.false.) + enddo + if (is_master()) print *, 'Random pattern for SFC-PERTS: k, min, max ',k, minval(rpattern_sfc(1)%spec_o(:,:,k)), maxval(rpattern_sfc(1)%spec_o(:,:,k)) + enddo ! k, nsfcpert + enddo ! n, npsfc + endif ! npsfc > 0 + if (is_master() .and. stochini) CLOSE(stochlun) + deallocate(noise_e,noise_o) + end subroutine init_stochdata + +subroutine read_pattern(rpattern,k,lunptn) + type(random_pattern), intent(inout) :: rpattern + integer, intent(in) :: lunptn + real(kind_dbl_prec),allocatable :: pattern2d(:),pattern2din(:) + real(kind_dbl_prec) :: stdevin,varin + integer nm,nn,ierr,jcap,isize,k + integer, allocatable :: isave(:) + + allocate(pattern2d(2*ndimspec)) + pattern2d=0. + call random_seed(size=isize,stat=rpattern%rstate) ! get size of generator state seed array + allocate(isave(isize)) + ! read only on root process, and send to all tasks + if (is_master()) then + read(lunptn) jcap + read(lunptn) isave + allocate(pattern2din((jcap+1)*(jcap+2))) + print*,'reading in random pattern at ',jcap,ndimspec,size(pattern2din) + read(lunptn) pattern2din + print*,'reading in random pattern (min/max/size/seed)',& + minval(pattern2din),maxval(pattern2din),size(pattern2din),isave(1:4) + if (jcap .eq. ntrunc) then + pattern2d=pattern2din + else + call chgres_pattern(pattern2din,pattern2d,jcap,ntrunc) ! chgres of spectral files + ! change the standard deviation of the patterns for a resolution change + ! needed for SKEB & SHUM + call computevarspec_r(rpattern,pattern2d,varin) + print*,'stddev in and out..',sqrt(varin),rpattern%stdev + stdevin=rpattern%stdev/sqrt(varin) + pattern2d(:)=pattern2d(:)*stdevin + endif + deallocate(pattern2din) + endif + call mp_bcst(isave,isize) ! blast out seed + call mp_bcst(pattern2d,2*ndimspec) + call random_seed(put=isave,stat=rpattern%rstate) + ! subset + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + rpattern%spec_e(nn,1,k) = pattern2d(nm) + rpattern%spec_e(nn,2,k) = pattern2d(ndimspec+nm) + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + rpattern%spec_o(nn,1,k) = pattern2d(nm) + rpattern%spec_o(nn,2,k) = pattern2d(ndimspec+nm) + enddo + !print*,'after scatter...',me,maxval(pattern2d_e),maxval(pattern2d_o) & + ! ,minval(pattern2d_e),minval(pattern2d_o) + deallocate(pattern2d,isave) + end subroutine read_pattern + +end module stochy_data_mod diff --git a/stochastic_physics/stochy_gg_def.f b/stochastic_physics/stochy_gg_def.f new file mode 100644 index 000000000..9470c436e --- /dev/null +++ b/stochastic_physics/stochy_gg_def.f @@ -0,0 +1,9 @@ + module stochy_gg_def + use machine + implicit none + + real(kind=kind_dbl_prec), allocatable, dimension(:) :: colrad_a, + & wgt_a, wgtcs_a, rcs2_a, sinlat_a, coslat_a +! + integer ,allocatable, dimension(:) :: lats_nodes_h,global_lats_h + end module stochy_gg_def diff --git a/stochastic_physics/stochy_internal_state_mod.F90 b/stochastic_physics/stochy_internal_state_mod.F90 new file mode 100644 index 000000000..e62645bd1 --- /dev/null +++ b/stochastic_physics/stochy_internal_state_mod.F90 @@ -0,0 +1,136 @@ + +! +! !module: stochy_internal_state_mod +! --- internal state definition of the +! gridded component of the spectral random patterns +! +! !description: define the spectral internal state used to +! create the internal state. +!--------------------------------------------------------------------------- +! !revision history: +! +! Oct 11 2016 P Pegion port of gfs_dynamics_interal_state +! +! !interface: +! + + module stochy_internal_state_mod + +!!uses: +!------ + use spectral_layout_mod + use stochy_gg_def + use stochy_resol_def + + + implicit none + private + +! ----------------------------------------------- + type,public::stochy_internal_state ! start type define +! ----------------------------------------------- + + integer :: me, nodes + integer :: lnt2_s, llgg_s + integer :: lnt2 + integer :: grib_inp + +! + integer nxpt,nypt,jintmx + integer lonf,latg,lats_node_a_max + + integer npe_single_member + + character(16) :: cfhour1 +!jws + integer :: num_file + character(32) ,allocatable :: filename_base(:) + integer :: ipt_lats_node_a + integer :: lats_node_a +!jwe + + integer :: nblck,kdt +! real :: deltim + + integer ,allocatable :: lonsperlat (:) + integer ,allocatable :: ls_node (:) + integer ,allocatable :: ls_nodes (:, :) + integer ,allocatable :: max_ls_nodes (:) + + integer ,allocatable :: lats_nodes_a (:) + integer ,allocatable :: global_lats_a (:) + integer ,allocatable :: lats_nodes_ext (:) + integer ,allocatable :: global_lats_ext(:) + integer ,allocatable :: global_lats_h (:) + integer :: xhalo,yhalo + + integer ,allocatable :: lats_nodes_a_fix (:) + + real(kind=kind_dbl_prec) ,allocatable :: epse (:) + real(kind=kind_dbl_prec) ,allocatable :: epso (:) + real(kind=kind_dbl_prec) ,allocatable :: epsedn(:) + real(kind=kind_dbl_prec) ,allocatable :: epsodn(:) + real(kind=kind_dbl_prec) ,allocatable :: kenorm_e(:) + real(kind=kind_dbl_prec) ,allocatable :: kenorm_o(:) + + real(kind=kind_dbl_prec) ,allocatable :: snnp1ev(:) + real(kind=kind_dbl_prec) ,allocatable :: snnp1od(:) + + real(kind=kind_dbl_prec) ,allocatable :: plnev_a(:,:) + real(kind=kind_dbl_prec) ,allocatable :: plnod_a(:,:) + real(kind=kind_dbl_prec) ,allocatable :: pddev_a(:,:) + real(kind=kind_dbl_prec) ,allocatable :: pddod_a(:,:) + real(kind=kind_dbl_prec) ,allocatable :: plnew_a(:,:) + real(kind=kind_dbl_prec) ,allocatable :: plnow_a(:,:) + + + real(kind=kind_dbl_prec) ,allocatable :: trie_ls(:,:,:) + real(kind=kind_dbl_prec) ,allocatable :: trio_ls(:,:,:) + + INTEGER :: TRIEO_TOTAL_SIZE + INTEGER, ALLOCATABLE, DIMENSION(:) :: TRIE_LS_SIZE + INTEGER, ALLOCATABLE, DIMENSION(:) :: TRIO_LS_SIZE + INTEGER, ALLOCATABLE, DIMENSION(:) :: TRIEO_LS_SIZE + INTEGER, ALLOCATABLE, DIMENSION(:) :: LS_MAX_NODE_GLOBAL + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: LS_NODE_GLOBAL + + +! + +!! + integer init,jcount,jpt,node,ibmsign,lon_dim,ilat + + real(kind=kind_dbl_prec) colat1, rone, rlons_lat, scale_ibm + + integer lotls,lotgr,lots,lots_slg,lotd,lota,lotp + + integer ibrad,ifges,ihour,ini,j,jdt,ksout,maxstp + integer mdt,idt,timetot,timer,time0 + integer mods,n1,n2,ndgf,ndgi,nfiles,nflps + integer n1hyb, n2hyb,nlunit + integer nges,ngpken,niter,nnmod,nradf,nradr + integer nsfcf,nsfci,nsfcs,nsigi,nsigs,nstep + integer nznlf,nznli,nznls,id,iret,nsout,ndfi + + integer ierr,iprint,k,l,locl,n + integer lan,lat + integer spectral_loop + + + integer ikey,nrank_all,kcolor + + real(kind=kind_dbl_prec) cons0p5,cons1200,cons3600,cons0 + +! +! ----------------------------------------------------- + end type stochy_internal_state ! end type define +! ----------------------------------------------------- + +! this state is supported by c pointer not f90 pointer, thus +! need this wrap. +!----------------------------------------------------------- + type stochy_wrap ! begin type define + type (stochy_internal_state), pointer :: int_state + end type stochy_wrap ! end type define + + end module stochy_internal_state_mod diff --git a/stochastic_physics/stochy_layout_lag.f b/stochastic_physics/stochy_layout_lag.f new file mode 100644 index 000000000..b2ee13bab --- /dev/null +++ b/stochastic_physics/stochy_layout_lag.f @@ -0,0 +1,13 @@ + module stochy_layout_lag + use machine + implicit none + save +cc + integer lats_dim_h, + x lats_node_h, + x lats_node_h_max, + x ipt_lats_node_h, + x lon_dim_h +cc + INTEGER ,ALLOCATABLE :: lat1s_h(:) + end module stochy_layout_lag diff --git a/stochastic_physics/stochy_namelist_def.F90 b/stochastic_physics/stochy_namelist_def.F90 new file mode 100644 index 000000000..06fac4f48 --- /dev/null +++ b/stochastic_physics/stochy_namelist_def.F90 @@ -0,0 +1,37 @@ + module stochy_namelist_def +! +! program log +! 11 Oct 2016: Philip Pegion create standalone stochastic physics +! + use machine + implicit none + + public + integer nsskeb,lon_s,lat_s,ntrunc + +! pjp stochastic phyics + integer skeb_varspect_opt,skeb_npass + logical sppt_sfclimit + + real(kind=kind_dbl_prec) :: skeb_sigtop1,skeb_sigtop2, & + sppt_sigtop1,sppt_sigtop2,shum_sigefold, & + skeb_vdof + real(kind=kind_dbl_prec) fhstoch,skeb_diss_smooth,skebint,skebnorm + real(kind=kind_dbl_prec), dimension(5) :: skeb,skeb_lscale,skeb_tau + real(kind=kind_dbl_prec), dimension(5) :: sppt,sppt_lscale,sppt_tau + real(kind=kind_dbl_prec), dimension(5) :: shum,shum_lscale,shum_tau + integer,dimension(5) ::skeb_vfilt + integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb + logical stochini,sppt_logit + logical do_shum,do_sppt,do_skeb,use_zmtnblck + +! mg surface perturbations + real(kind=kind_dbl_prec), dimension(5) :: sfc_lscale,sfc_tau + real(kind=kind_dbl_prec), dimension(5) :: pertz0,pertshc,pertzt + real(kind=kind_dbl_prec), dimension(5) :: pertlai,pertvegf,pertalb + integer nsfcpert + integer(8),dimension(5) ::iseed_sfc + logical sppt_land + logical do_sfcperts + + end module stochy_namelist_def diff --git a/stochastic_physics/stochy_patterngenerator.F90 b/stochastic_physics/stochy_patterngenerator.F90 new file mode 100644 index 000000000..7b26e69d5 --- /dev/null +++ b/stochastic_physics/stochy_patterngenerator.F90 @@ -0,0 +1,357 @@ +module stochy_patterngenerator_mod + + ! generate random patterns with specified temporal and spatial auto-correlation + ! in spherical harmonic space. + use machine + use spectral_layout_mod, only: len_trie_ls, len_trio_ls, ls_dim, ls_max_node +! use mersenne_twister_stochy, only: random_setseed,random_gauss,random_stat + use mersenne_twister, only: random_setseed,random_gauss,random_stat + use stochy_ccpp, only: is_master, mp_bcst + implicit none + private + + public :: computevarspec, setvarspect,& + patterngenerator_init, patterngenerator_destroy, getnoise, & + patterngenerator_advance, random_pattern, ndimspec,& + chgres_pattern,computevarspec_r + + type random_pattern + real(kind_dbl_prec), public :: lengthscale + real(kind_dbl_prec), public :: tau + real(kind_dbl_prec), public :: dt + real(kind_dbl_prec), public :: phi + real(kind_dbl_prec), public :: stdev + real(kind_evod), allocatable, dimension(:), public :: varspectrum, varspectrum1d, lap + integer, allocatable, dimension(:), public ::& + degree,order,idx_e,idx_o + integer, allocatable, dimension(:,:), public :: idx + integer, public :: seed + real(kind_dbl_prec), allocatable, dimension(:,:,:), public :: spec_e,spec_o + type(random_stat), public :: rstate + end type random_pattern + + integer :: nlons,nlats,ntrunc,ndimspec + + contains + + subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& + nlon, nlat, jcap, ls_node, npatterns,& + nlevs, varspect_opt) + real(kind_dbl_prec), intent(in),dimension(npatterns) :: lscale,tscale,stdev + real, intent(in) :: delt + integer, intent(in) :: nlon,nlat,jcap,npatterns,varspect_opt + integer, intent(in) :: ls_node(ls_dim,3),nlevs + type(random_pattern), intent(out), dimension(npatterns) :: rpattern + integer(8), intent(inout) :: iseed(npatterns) + integer m,j,l,n,nm,nn,np,indev1,indev2,indod1,indod2 + integer(8) count, count_rate, count_max, count_trunc + integer(8) :: iscale = 10000000000 + integer count4, ierr +! integer member_id + integer indlsod,indlsev,jbasev,jbasod + include 'function_indlsod' + include 'function_indlsev' + nlons = nlon + nlats = nlat + ntrunc = jcap + ndimspec = (ntrunc+1)*(ntrunc+2)/2 +! propagate seed supplied from namelist to all patterns... + if (iseed(1) .NE. 0) then + do np=2,npatterns + if (iseed(np).EQ.0) then + iseed(np)=iseed(1)+np*100000000 + endif + enddo + endif + + do np=1,npatterns + allocate(rpattern(np)%idx(0:ntrunc,0:ntrunc)) + allocate(rpattern(np)%idx_e(len_trie_ls)) + allocate(rpattern(np)%idx_o(len_trio_ls)) + allocate(rpattern(np)%spec_e(len_trie_ls,2,nlevs)) + allocate(rpattern(np)%spec_o(len_trio_ls,2,nlevs)) + rpattern(np)%idx_e = 0; rpattern(np)%idx_o = 0; rpattern(np)%idx = 0 + rpattern(np)%spec_e(:,:,:)=0. + rpattern(np)%spec_o(:,:,:)=0. + nm = 0 + do m=0,ntrunc + do n=m,ntrunc + nm = nm + 1 + rpattern(np)%idx(m,n) = nm + enddo + enddo + do j = 1, ls_max_node + l=ls_node(j,1) ! zonal wavenumber + jbasev=ls_node(j,2) + jbasod=ls_node(j,3) + indev1 = indlsev(l,l) + indod1 = indlsod(l+1,l) + if (mod(l,2) .eq. mod(ntrunc+1,2)) then + indev2 = indlsev(ntrunc+1,l) + indod2 = indlsod(ntrunc ,l) + else + indev2 = indlsev(ntrunc ,l) + indod2 = indlsod(ntrunc+1,l) + endif + n = l ! degree + do nn=indev1,indev2 + if (n <= ntrunc .and. l <= ntrunc) then + nm = rpattern(np)%idx(l,n) + rpattern(np)%idx_e(nn) = nm + endif + n = n + 2 + enddo + n = l+1 + do nn=indod1,indod2 + if (n <= ntrunc .and. l <= ntrunc) then + nm = rpattern(np)%idx(l,n) + rpattern(np)%idx_o(nn) = nm + endif + n = n + 2 + enddo + enddo + allocate(rpattern(np)%degree(ndimspec),rpattern(np)%order(ndimspec),rpattern(np)%lap(ndimspec)) +#ifdef __GFORTRAN__ + j = 0 + do m=0,ntrunc + do n=m,ntrunc + j = j + 1 + rpattern(np)%degree(j) = n + rpattern(np)%order(j) = m + end do + end do +#else + rpattern(np)%degree = (/((n,n=m,ntrunc),m=0,ntrunc)/) + rpattern(np)%order = (/((m,n=m,ntrunc),m=0,ntrunc)/) +#endif + rpattern(np)%lap = -rpattern(np)%degree*(rpattern(np)%degree+1.0) + rpattern(np)%tau = tscale(np) + rpattern(np)%lengthscale = lscale(np) + rpattern(np)%dt = delt + rpattern(np)%phi = exp(-delt/tscale(np)) + rpattern(np)%stdev = stdev(np) + allocate(rpattern(np)%varspectrum(ndimspec)) + allocate(rpattern(np)%varspectrum1d(0:ntrunc)) + ! seed computed on root, then bcast to all tasks and set. + if (is_master()) then +! read(ens_nam(2:3),'(i2)') member_id +! print *,'ens_nam,member_id',trim(ens_nam),member_id + if (iseed(np) == 0) then + ! generate a random seed from system clock and ens member number + call system_clock(count, count_rate, count_max) + ! iseed is elapsed time since unix epoch began (secs) + ! truncate to 4 byte integer + count_trunc = iscale*(count/iscale) + count4 = count - count_trunc !+ member_id + print *,'using seed',count4 + else + !count4 = iseed(np) + member_id + ! don't rely on compiler to truncate integer(8) to integer(4) on + ! overflow, do wrap around explicitly. + !count4 = mod(iseed(np) + member_id + 2147483648, 4294967296) - 2147483648 + count4 = mod(iseed(np) + 2147483648, 4294967296) - 2147483648 + print *,'using seed',count4,iseed(np)!,member_id + endif + endif + ! broadcast seed to all tasks. + call mp_bcst(count4) + rpattern(np)%seed = count4 + ! set seed (to be the same) on all tasks. Save random state. + call random_setseed(rpattern(np)%seed,rpattern(np)%rstate) + if (varspect_opt .ne. 0 .and. varspect_opt .ne. 1) then + if (is_master()) then + print *,'WARNING: illegal value for varspect_opt (should be 0 or 1), using 0 (gaussian spectrum)...' + endif + call setvarspect(rpattern(np),0) + else + call setvarspect(rpattern(np),varspect_opt) + endif + enddo ! n=1,npatterns + end subroutine patterngenerator_init + + + + subroutine patterngenerator_destroy(rpattern,npatterns) + type(random_pattern), intent(inout) :: rpattern(npatterns) + integer, intent(in) :: npatterns + integer n + do n=1,npatterns + deallocate(rpattern(n)%varspectrum,rpattern(n)%varspectrum1d) + deallocate(rpattern(n)%degree,rpattern(n)%order,rpattern(n)%lap) + deallocate(rpattern(n)%idx,rpattern(n)%idx_e,rpattern(n)%idx_o) + enddo + end subroutine patterngenerator_destroy + + subroutine computevarspec(rpattern,dataspec,var) + ! compute globally integrated variance from spectral coefficients + complex(kind_evod), intent(in) :: dataspec(ndimspec) + real(kind_evod), intent(out) :: var + type(random_pattern), intent(in) :: rpattern + integer n + var = 0. + do n=1,ndimspec + if (rpattern%order(n) .ne. 0) then + var = var + dataspec(n)*conjg(dataspec(n)) + else + var = var + 0.5*dataspec(n)*conjg(dataspec(n)) + endif + enddo + end subroutine computevarspec + + subroutine computevarspec_r(rpattern,dataspec,var) + ! compute globally integrated variance from spectral coefficients + real(kind_dbl_prec), intent(in) :: dataspec(2*ndimspec) + real(kind_dbl_prec), intent(out) :: var + type(random_pattern), intent(in) :: rpattern + integer n + var = 0. + do n=1,ndimspec + if (rpattern%order(n) .ne. 0) then + var = var + dataspec(n)**2+dataspec(n+ndimspec)**2 + else + var = var + 0.5*(dataspec(n)**2+dataspec(n+ndimspec)**2) + endif + enddo + end subroutine computevarspec_r + + subroutine getnoise(rpattern,noise_e,noise_o) + real(kind_dbl_prec), intent(out) :: noise_e(len_trie_ls,2) + real(kind_dbl_prec), intent(out) :: noise_o(len_trio_ls,2) + ! generate white noise with unit variance in spectral space + type(random_pattern), intent(inout) :: rpattern + real :: noise(2*ndimspec) + integer nm,nn + call random_gauss(noise,rpattern%rstate) + noise(1) = 0.; noise(ndimspec+1) = 0. + noise = noise*sqrt(1./ntrunc) + noise_e = 0.; noise_o = 0. + ! subset + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + noise_e(nn,1) = noise(nm)/sqrt(2.*rpattern%degree(nm)+1) + noise_e(nn,2) = noise(ndimspec+nm)/sqrt(2.*rpattern%degree(nm)+1) + if (rpattern%order(nm) .eq. 0) then + noise_e(nn,1) = sqrt(2.)*noise_e(nn,1) + noise_e(nn,2) = 0. + endif + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + noise_o(nn,1) = noise(nm)/sqrt(2.*rpattern%degree(nm)+1) + noise_o(nn,2) = noise(ndimspec+nm)/sqrt(2.*rpattern%degree(nm)+1) + if (rpattern%order(nm) .eq. 0) then + noise_o(nn,1) = sqrt(2.)*noise_o(nn,1) + noise_o(nn,2) = 0. + endif + enddo + end subroutine getnoise + + subroutine patterngenerator_advance(rpattern,k,skeb_first_call) + ! advance 1st-order autoregressive process with + ! specified autocorrelation (phi) and variance spectrum (spectrum) + real(kind_dbl_prec) :: noise_e(len_trie_ls,2) + real(kind_dbl_prec) :: noise_o(len_trio_ls,2) + type(random_pattern), intent(inout) :: rpattern + logical, intent(in) :: skeb_first_call + integer j,l,n,nn,nm,k,k2 + call getnoise(rpattern,noise_e,noise_o) + if (k.GT.1.AND.skeb_first_call) then + k2=k-1 + else + k2=k + endif + do nn=1,len_trie_ls + nm = rpattern%idx_e(nn) + if (nm == 0) cycle + rpattern%spec_e(nn,1,k) = rpattern%phi*rpattern%spec_e(nn,1,k2) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*rpattern%varspectrum(nm)*noise_e(nn,1) + rpattern%spec_e(nn,2,k) = rpattern%phi*rpattern%spec_e(nn,2,k2) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*rpattern%varspectrum(nm)*noise_e(nn,2) + enddo + do nn=1,len_trio_ls + nm = rpattern%idx_o(nn) + if (nm == 0) cycle + rpattern%spec_o(nn,1,k) = rpattern%phi*rpattern%spec_o(nn,1,k2) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*rpattern%varspectrum(nm)*noise_o(nn,1) + rpattern%spec_o(nn,2,k) = rpattern%phi*rpattern%spec_o(nn,2,k2) + & + rpattern%stdev*sqrt(1.-rpattern%phi**2)*rpattern%varspectrum(nm)*noise_o(nn,2) + enddo + end subroutine patterngenerator_advance + + subroutine setvarspect(rpattern,varspect_opt) + ! define variance spectrum (isotropic covariance) + ! normalized to unit global variance + type(random_pattern), intent(inout) :: rpattern + integer, intent(in) :: varspect_opt + integer :: n + complex(kind_evod) noise(ndimspec) + real(kind_evod) var,rerth + rerth =6.3712e+6 ! radius of earth (m) + ! 1d variance spectrum (as a function of total wavenumber) + if (varspect_opt == 0) then ! gaussian + ! rpattern%lengthscale is interpreted as an efolding length + ! scale, in meters. + do n=0,ntrunc + rpattern%varspectrum1d(n) = exp(-rpattern%lengthscale**2*(float(n)*(float(n)+1.))/(4.*rerth**2)) + enddo + ! scaling factors for spectral coeffs of white noise pattern with unit variance + rpattern%varspectrum = sqrt(ntrunc*exp(rpattern%lengthscale**2*rpattern%lap/(4.*rerth**2))) + else if (varspect_opt == 1) then ! power law + ! rpattern%lengthscale is interpreted as a power, not a length. + do n=0,ntrunc + rpattern%varspectrum1d(n) = float(n)**(rpattern%lengthscale) + enddo + ! scaling factors for spectral coeffs of white noise pattern with unit variance + rpattern%varspectrum = sqrt(ntrunc*(rpattern%degree**(rpattern%lengthscale))) + endif + noise = 0. + do n=1,ndimspec + if (rpattern%order(n) .ne. 0.) then + noise(n) = cmplx(1.,1.)/sqrt(2.*rpattern%degree(n)+1) + else + noise(n) = sqrt(2.)/sqrt(2.*rpattern%degree(n)+1.) + endif + enddo + noise(1) = 0 ! no global mean. + ! make sure global mean variance is 1. + noise = noise*sqrt(1./ntrunc) + noise = rpattern%varspectrum*noise + call computevarspec(rpattern,noise,var) + rpattern%varspectrum = rpattern%varspectrum/sqrt(var) + rpattern%varspectrum1d = rpattern%varspectrum1d/var + + end subroutine setvarspect + + subroutine chgres_pattern(pattern2din,pattern2dout,ntruncin,ntruncout) + real(kind_dbl_prec), intent(in) :: pattern2din((ntruncin+1)*(ntruncin+2)) + real(kind_dbl_prec), intent(out) :: pattern2dout((ntruncout+1)*(ntruncout+2)) + integer, intent(in) :: ntruncin,ntruncout + integer :: m,n,nm,ndimsspecin,ndimsspecout + integer,allocatable, dimension(:,:):: idxin + allocate(idxin(0:ntruncin,0:ntruncin)) + ndimsspecin=(ntruncin+1)*(ntruncin+2)/2 + ndimsspecout=(ntruncout+1)*(ntruncout+2)/2 + nm = 0 + do m=0,ntruncin + do n=m,ntruncin + nm = nm + 1 + idxin(m,n) = nm + enddo + enddo + ! chgres + nm = 0 + do m=0,ntruncout + do n=m,ntruncout + nm = nm + 1 + if (m .le. ntruncin .and. n .le. ntruncin) then + pattern2dout(nm) = pattern2din(idxin(m,n)) + pattern2dout(ndimsspecout+nm) = pattern2din(ndimsspecin+idxin(m,n)) + endif + enddo + enddo + deallocate(idxin) +end subroutine chgres_pattern + +end module stochy_patterngenerator_mod diff --git a/stochastic_physics/stochy_resol_def.f b/stochastic_physics/stochy_resol_def.f new file mode 100644 index 000000000..708e31c84 --- /dev/null +++ b/stochastic_physics/stochy_resol_def.f @@ -0,0 +1,44 @@ + module stochy_resol_def + +! program log: +! 20110220: Henry Juang update index for MASS_DP and NDSLFV +! 20130202: Henry Juang revise reduced grid and add x number +! + implicit none + + integer jcap,jcap1,jcap2,latg,latg2 + integer levh,levm1,levp1,skeblevs,levs,lnt,lnt2,lnt22,levr + integer lnte,lnted,lnto,lntod,lnuv + integer lonf,lonfx,num_p2d,num_p3d + integer nxpt,nypt,jintmx,latgd + integer ntoz,ntcw,ncld,ntke,ixgr,ntiw,ntlnc,ntinc,nto,nto2 + integer ivsupa, ivsinp + integer nlunit, kdt_start + integer,target :: ntrac + integer,target :: ngrids_gg + integer,target :: thermodyn_id, sfcpress_id ! hmhj + logical,target :: adiabatic +! + INTEGER p_gz,p_lapgz,p_zslam,p_zsphi,p_dlam,p_dphi,p_uln,p_vln + INTEGER p_zem,p_dim,p_tem,p_rm,p_dpm,p_qm + INTEGER p_ze ,p_di ,p_te ,p_rq,p_dp ,p_q + INTEGER p_w ,p_x ,p_y ,p_rt,p_dpn,p_zq + INTEGER p_zz ,p_dpphi,p_dplam,p_zzphi,p_zzlam + INTEGER g_uum,g_vvm,g_ttm,g_rm ,g_dpm,g_qm,g_gz,g_zz + INTEGER g_uu ,g_vv ,g_tt ,g_rq ,g_dp ,g_q + INTEGER g_uup,g_vvp,g_ttp,g_rqp,g_dpp,g_zqp,g_rqtk + INTEGER g_u ,g_v ,g_t ,g_rt ,g_dpn,g_zq, g_p, g_dpdt + INTEGER lots,lots_slg,lotd,lota,lotp,lotls,lotgr,lotgr6 + + integer ksz, ksd, kst, ksr, ksdp, ksq, ksplam, kspphi + integer ksu, ksv, kzslam, kzsphi +! + integer kau, kav, kat, kar, kadp, kaps, kazs, kap2 +! + integer kdpphi, kzzphi, kdplam, kzzlam +! + integer kdtphi, kdrphi, kdtlam, kdrlam + integer kdulam, kdvlam, kduphi, kdvphi + + + end module stochy_resol_def diff --git a/stochastic_physics/sumfln_stochy.f b/stochastic_physics/sumfln_stochy.f new file mode 100644 index 000000000..973828e1b --- /dev/null +++ b/stochastic_physics/sumfln_stochy.f @@ -0,0 +1,294 @@ + module sumfln_stochy_mod + + implicit none + + contains + + subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, + & nvars,ls_node,latl2, + & workdim,nvarsdim,four_gr, + & ls_nodes,max_ls_nodes, + & lats_nodes,global_lats, + & lats_node,ipt_lats_node, + & lons_lat,londi,latl,nvars_0) +! + use stochy_resol_def , only : jcap,latgd + use spectral_layout_mod , only : len_trie_ls,len_trio_ls, + & ls_dim,ls_max_node,me,nodes + use machine + use stochy_ccpp, only : mpp_alltoall, + & num_parthds_stochy => ompthreads + + implicit none +! + external esmf_dgemm +! + integer lat1s(0:jcap),latl2 +! + integer nvars,nvars_0 + real(kind=kind_dbl_prec) flnev(len_trie_ls,2*nvars) + real(kind=kind_dbl_prec) flnod(len_trio_ls,2*nvars) +! + real(kind=kind_dbl_prec) plnev(len_trie_ls,latl2) + real(kind=kind_dbl_prec) plnod(len_trio_ls,latl2) +! + integer ls_node(ls_dim,3) +! +!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L +!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev +!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod +! +! local scalars +! ------------- +! + integer j, k, l, lat, lat1, n, kn, n2,indev,indod +! +! local arrays +! ------------ +! + real(kind=kind_dbl_prec), dimension(nvars*2,latl2) :: apev, apod + integer num_threads, nvar_thread_max, nvar_1, nvar_2 + &, thread +! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx +! + integer nvarsdim, latl, workdim, londi + &, lats_node, ipt_lats_node +! + real(kind=kind_dbl_prec) four_gr(londi,nvarsdim,workdim) +! + integer ls_nodes(ls_dim,nodes) + integer, dimension(nodes) :: max_ls_nodes, lats_nodes + integer, dimension(latl) :: global_lats, lons_lat + +!jfe integer global_lats(latg+2*jintmx+2*nypt*(nodes-1)) +! + real(kind=4),target,dimension(2,nvars,ls_dim*workdim,nodes):: + & workr,works +! real(kind=4),dimension(2*nvars*ls_dim*workdim*nodes):: +! & work1dr,work1ds + real(kind=4),pointer:: work1dr(:),work1ds(:) + integer, dimension(jcap+1) :: kpts, kptr, sendcounts, recvcounts, + & sdispls +! + integer ierr,ilat,ipt_ls, lmax,lval,i,jj,lonl,nv + integer node,nvar,arrsz + integer ilat_list(nodes) ! for OMP buffer copy +! +! statement functions +! ------------------- +! + integer indlsev, jbasev, indlsod, jbasod +! + include 'function_indlsev' + include 'function_indlsod' +! + real(kind=kind_dbl_prec), parameter :: cons0=0.0d0, cons1=1.0d0 +! +! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx +! + arrsz=2*nvars*ls_dim*workdim*nodes + num_threads = min(num_parthds_stochy,nvars) + nvar_thread_max = (nvars+num_threads-1)/num_threads + kpts = 0 +! write(0,*)' londi=',londi,'nvarsdim=',nvarsdim,'workdim=',workdim +! + do j = 1, ls_max_node ! start of do j loop ##################### +! + l = ls_node(j,1) + jbasev = ls_node(j,2) + jbasod = ls_node(j,3) + + indev = indlsev(l,l) + indod = indlsod(l+1,l) +! + lat1 = lat1s(l) + if ( kind_dbl_prec == 8 ) then !------------------------------------ + +!$omp parallel do private(thread,nvar_1,nvar_2,n2) + do thread=1,num_threads ! start of thread loop .............. + nvar_1 = (thread-1)*nvar_thread_max + 1 + nvar_2 = min(nvar_1+nvar_thread_max-1,nvars) + + if (nvar_2 >= nvar_1) then + n2 = 2*(nvar_2-nvar_1+1) + +! compute the even and odd components of the fourier coefficients +! +! compute the sum of the even real terms for each level +! compute the sum of the even imaginary terms for each level +! +! call dgemm('t','n',latl2-lat1+1, 2*(nvar_2-nvar_1+1), +! & (jcap+2-l)/2,cons1, !constant +! & plnev(indev,lat1), len_trio_ls, +! & flnev(indev,2*nvar_1-1),len_trio_ls,cons0, +! & apev(2*nvar_1-1,lat1),latl2) + call esmf_dgemm( + & 't', + & 'n', + & n2, + & latl2-lat1+1, + & (jcap+3-l)/2, + & cons1, + & flnev(indev,2*nvar_1-1), + & len_trie_ls, + & plnev(indev,lat1), + & len_trie_ls, + & cons0, + & apev(2*nvar_1-1,lat1), + & 2*nvars + & ) +! +! compute the sum of the odd real terms for each level +! compute the sum of the odd imaginary terms for each level +! +! call dgemm('t','n',latl2-lat1+1, 2*(nvar_2-nvar_1+1), +! & (jcap+2-l)/2,cons1, !constant +! & plnod(indod,lat1), len_trio_ls, +! & flnod(indod,2*nvar_1-1),len_trio_ls,cons0, +! & apod(2*nvar_1-1,lat1), latl2) + call esmf_dgemm( + & 't', + & 'n', + & n2, + & latl2-lat1+1, + & (jcap+2-l)/2, + & cons1, + & flnod(indod,2*nvar_1-1), + & len_trio_ls, + & plnod(indod,lat1), + & len_trio_ls, + & cons0, + & apod(2*nvar_1-1,lat1), + & 2*nvars + & ) +! + endif + enddo ! end of thread loop .................................. + else !------------------------------------------------------------ +!$omp parallel do private(thread,nvar_1,nvar_2) + do thread=1,num_threads ! start of thread loop .............. + nvar_1 = (thread-1)*nvar_thread_max + 1 + nvar_2 = min(nvar_1+nvar_thread_max-1,nvars) + enddo ! end of thread loop .................................. + endif !----------------------------------------------------------- +! +ccxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx +! +! compute the fourier coefficients for each level +! ----------------------------------------------- +! + ilat_list(1) = 0 + do node = 1, nodes - 1 + ilat_list(node+1) = ilat_list(node) + lats_nodes(node) + end do + +!$omp parallel do private(node,jj,ilat,lat,ipt_ls,nvar,kn,n2) + do node=1,nodes + do jj=1,lats_nodes(node) + ilat = ilat_list(node) + jj + lat = global_lats(ilat) + ipt_ls = min(lat,latl-lat+1) + if ( ipt_ls >= lat1s(ls_nodes(j,me+1)) ) then + kpts(node) = kpts(node) + 1 + kn = kpts(node) +! + if ( lat <= latl2 ) then +! northern hemisphere + do nvar=1,nvars + n2 = nvar + nvar + works(1,nvar,kn,node) = apev(n2-1,ipt_ls) + & + apod(n2-1,ipt_ls) + works(2,nvar,kn,node) = apev(n2, ipt_ls) + & + apod(n2, ipt_ls) + enddo + else +! southern hemisphere + do nvar=1,nvars + n2 = nvar + nvar + works(1,nvar,kn,node) = apev(n2-1,ipt_ls) + & - apod(n2-1,ipt_ls) + works(2,nvar,kn,node) = apev(n2, ipt_ls) + & - apod(n2, ipt_ls) + enddo + endif + endif + enddo + enddo +! + enddo ! end of do j loop ####################################### +! + kptr = 0 + do node=1,nodes + do l=1,max_ls_nodes(node) + lval = ls_nodes(l,node)+1 + do j=1,lats_node + lat = global_lats(ipt_lats_node-1+j) + if ( min(lat,latl-lat+1) >= lat1s(lval-1) ) then + kptr(node) = kptr(node) + 1 + endif + enddo + enddo + enddo +! +! + n2 = nvars + nvars +!$omp parallel do private(node) + do node=1,nodes + sendcounts(node) = kpts(node) * n2 + recvcounts(node) = kptr(node) * n2 + sdispls(node) = (node-1) * n2 * ls_dim * workdim + end do + work1dr(1:arrsz)=>workr + work1ds(1:arrsz)=>works + call mpp_alltoall(work1ds, sendcounts, sdispls, + & work1dr, recvcounts, sdispls) + nullify(work1dr) + nullify(work1ds) +!$omp parallel do private(j,lat,lmax,nvar,lval,n2,lonl,nv) + do j=1,lats_node + lat = global_lats(ipt_lats_node-1+j) + lonl = lons_lat(lat) + lmax = min(jcap,lonl/2) + n2 = lmax + lmax + 3 +! write(0,*)' j=',j,' lat=',lat,' lmax=',lmax,' n2=',n2 +! &,' nvars=',nvars,' lonl=',lonl + if ( n2 <= lonl+2 ) then + do nvar=1,nvars + nv = nvars_0 + nvar + do lval = n2, lonl+2 +! write(0,*)' lval=',lval,' nvar=',nvar,nvars_0 +! &,' n2=',n2,' lonl=',lonl,' nv=',nv,' j=',j +! &,'size=',size(four_gr,1),size(four_gr,2),size(four_gr,3) + four_gr(lval,nv,j) = cons0 + enddo + enddo + endif + enddo +! + kptr = 0 +! write(0,*)' kptr=',kptr(1) +!! +!$omp parallel do private(node,l,lval,j,lat,nvar,kn,n2) + do node=1,nodes + do l=1,max_ls_nodes(node) + lval = ls_nodes(l,node)+1 + n2 = lval + lval + do j=1,lats_node + lat = global_lats(ipt_lats_node-1+j) + if ( min(lat,latl-lat+1) >= lat1s(lval-1) ) then + kptr(node) = kptr(node) + 1 + kn = kptr(node) + + do nvar=1,nvars + four_gr(n2-1,nvars_0+nvar,j) = workr(1,nvar,kn,node) + four_gr(n2, nvars_0+nvar,j) = workr(2,nvar,kn,node) + enddo + endif + enddo + enddo + enddo +! + return + end + + end module sumfln_stochy_mod From 077b45558980b7b42a3287e2c0ab813307a96c4b Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 2 Aug 2018 15:52:37 -0600 Subject: [PATCH 4/4] physics/GFS_stochastics.F90: insert code from GFS_driver.F90 ('kludge for output', assigns stochastic weights in reverse vertical order from Coupling to Diag), add missing variables and their metadata --- physics/GFS_stochastics.F90 | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/physics/GFS_stochastics.F90 b/physics/GFS_stochastics.F90 index 6a4eba767..390bffd5a 100644 --- a/physics/GFS_stochastics.F90 +++ b/physics/GFS_stochastics.F90 @@ -22,10 +22,13 @@ end subroutine GFS_stochastics_finalize !! | do_skeb | flag_for_stochastic_skeb_option | flag for stochastic skeb option | flag | 0 | logical | | in | F | !! | zmtnblck | level_of_dividing_streamline | level of the dividing streamline | none | 1 | real | kind_phys | in | F | !! | sppt_wts | weights_for_stochastic_surface_physics_perturbation | weights for stochastic surface physics perturbation | none | 2 | real | kind_phys | inout | F | -!! | sppt_wts_inv | weights_for_stochastic_surface_physics_perturbation_flipped | weights for stochastic surface physics perturbation, flipped | none | 2 | real | kind_phys | out | F | !! | skebu_wts | weights_for_stochastic_skeb_perturbation_of_x_wind | weights for stochastic skeb perturbation of x wind | none | 2 | real | kind_phys | in | F | !! | skebv_wts | weights_for_stochastic_skeb_perturbation_of_y_wind | weights for stochastic skeb perturbation of y wind | none | 2 | real | kind_phys | in | F | !! | shum_wts | weights_for_stochastic_shum_perturbation | weights for stochastic shum perturbation | none | 2 | real | kind_phys | in | F | +!! | sppt_wts_inv | weights_for_stochastic_surface_physics_perturbation_flipped | weights for stochastic surface physics perturbation, flipped | none | 2 | real | kind_phys | inout | F | +!! | skebu_wts_inv | weights_for_stochastic_skeb_perturbation_of_x_wind_flipped | weights for stochastic skeb perturbation of x wind, flipped | none | 2 | real | kind_phys | inout | F | +!! | skebv_wts_inv | weights_for_stochastic_skeb_perturbation_of_y_wind_flipped | weights for stochastic skeb perturbation of y wind, flipped | none | 2 | real | kind_phys | inout | F | +!! | shum_wts_inv | weights_for_stochastic_shum_perturbation_flipped | weights for stochastic shum perturbation, flipped | none | 2 | real | kind_phys | inout | F | !! | diss_est | dissipation_estimate_of_air_temperature_at_model_layers | dissipation estimate model layer mean temperature | K | 2 | real | kind_phys | in | F | !! | ugrs | x_wind | zonal wind | m s-1 | 2 | real | kind_phys | in | F | !! | vgrs | y_wind | meridional wind | m s-1 | 2 | real | kind_phys | in | F | @@ -60,8 +63,9 @@ end subroutine GFS_stochastics_finalize ! 6) performs surface data cycling via the GFS gcycle routine !------------------------------------------------------------------------- subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, & - zmtnblck, sppt_wts, sppt_wts_inv, skebu_wts, & - skebv_wts, shum_wts, diss_est, & + zmtnblck, sppt_wts, skebu_wts, skebv_wts, shum_wts,& + sppt_wts_inv, skebu_wts_inv, skebv_wts_inv, & + shum_wts_inv, diss_est, & ugrs, vgrs, tgrs, qgrs, gu0, gv0, gt0, gq0, dtdtr, & rain, rainc, tprcp, totprcp, cnvprcp, & cplflx, rain_cpl, snow_cpl, drain_cpl, dsnow_cpl, & @@ -78,14 +82,18 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, logical, intent(in) :: do_shum logical, intent(in) :: do_skeb real(kind_phys), dimension(1:im), intent(in) :: zmtnblck - ! sppt_wts only allocated if do_sppt == .true. (sppt_wts_inv is always allocated) + ! sppt_wts only allocated if do_sppt == .true. real(kind_phys), dimension(:,:), intent(inout) :: sppt_wts - real(kind_phys), dimension(1:im,1:km), intent(out) :: sppt_wts_inv ! skebu_wts, skebv_wts only allocated if do_skeb == .true. real(kind_phys), dimension(:,:), intent(in) :: skebu_wts real(kind_phys), dimension(:,:), intent(in) :: skebv_wts ! shum_wts only allocated if do_shum == .true. real(kind_phys), dimension(:,:), intent(in) :: shum_wts + ! inverse/flipped weights are always allocated + real(kind_phys), dimension(1:im,1:km), intent(inout) :: sppt_wts_inv + real(kind_phys), dimension(1:im,1:km), intent(inout) :: skebu_wts_inv + real(kind_phys), dimension(1:im,1:km), intent(inout) :: skebv_wts_inv + real(kind_phys), dimension(1:im,1:km), intent(inout) :: shum_wts_inv real(kind_phys), dimension(1:im,1:km), intent(in) :: diss_est real(kind_phys), dimension(1:im,1:km), intent(in) :: ugrs real(kind_phys), dimension(1:im,1:km), intent(in) :: vgrs @@ -172,16 +180,20 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, endif endif - + if (do_shum) then - gq0(:,:) = gq0(:,:)*(1.0 + shum_wts(:,:)) + do k=1,km + gq0(:,k) = gq0(:,k)*(1.0 + shum_wts(:,k)) + shum_wts_inv(:,km-k+1) = shum_wts(:,k) + end do endif if (do_skeb) then do k=1,km - gu0(:,k) = gu0(:,k)+skebu_wts(:,k)*(diss_est(:,k)) - gv0(:,k) = gv0(:,k)+skebv_wts(:,k)*(diss_est(:,k)) - ! print*,'in do skeb',skebu_wts(1,k),diss_est(1,k) + gu0(:,k) = gu0(:,k)+skebu_wts(:,k)*(diss_est(:,k)) + gv0(:,k) = gv0(:,k)+skebv_wts(:,k)*(diss_est(:,k)) + skebu_wts_inv(:,km-k+1) = skebu_wts(:,k) + skebv_wts_inv(:,km-k+1) = skebv_wts(:,k) enddo endif