diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index ca66966f56..3eeaaf3499 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -220,6 +220,53 @@ Default: FALSE + + Toggle Spectral Filtering of Nudging Tendencies ON/OFF. + + Nudge_SpectralFilter - LOGICAL Option to apply spherical harminic filtering to + the model state and target data so that nudging + tendencies are only applied to scales larger than + the specified truncation. + Default: FALSE + + + + Set Horizonal Scale for Spectral Filtering + + Nudge_SpectralNtrunc - INT The number of meridional spherical harmonic modes used + for spectral filtering. The nominal horizontal scale of + the filtering can be estimated as: + + Hscale = PI*6350/Nudge_SpectralNtrunc + + i.e. Nudge_SpectralNtrunc=40 corresponds to a horizontal + nudging scale Hscale~500km. + Default: -1 + + + + Set size of local-area grid used to suppress Sampling Errors. + + Nudge_SpectralNring - INT The number of sampling rings used for local area averaging + of spherical harmonic modes, to suppress sampling errors. + When initializing each basis, a local average of SH values + is computed for the area associated with each grid point. + SpectralNring set the number of rings of equal-area points + in this sampling domain. + Each ring (kk) contains 8*(kk-1) sample points. + + Nudge_SpectralNring Number of Sampling Points + ------------------- ------------------------- + 1 1 (DEFAULT SampleGrid NOT used) + 2 9 + 3 25 + 4 49 + Default: 1 + + Full pathname of analyses data to use for nudging. @@ -307,6 +354,9 @@ 0 == OFF (No Nudging of this variable) 1 == CONSTANT (Spatially Uniform Nudging) 2 == HEAVISIDE WINDOW FUNCTION + 3 == HEAVISIDE WINDOW FUNCTION IN VERTICAL + and USER SPECIFIED HORIZONAL DOMAIN + (Specified by user via Nudge_Bwindow_File) Default: 0 @@ -323,6 +373,9 @@ 0 == OFF (No Nudging of this variable) 1 == CONSTANT (Spatially Uniform Nudging) 2 == HEAVISIDE WINDOW FUNCTION + 3 == HEAVISIDE WINDOW FUNCTION IN VERTICAL + and USER SPECIFIED HORIZONAL DOMAIN + (Specified by user via Nudge_Bwindow_File) Default: 0 @@ -339,6 +392,9 @@ 0 == OFF (No Nudging of this variable) 1 == CONSTANT (Spatially Uniform Nudging) 2 == HEAVISIDE WINDOW FUNCTION + 3 == HEAVISIDE WINDOW FUNCTION IN VERTICAL + and USER SPECIFIED HORIZONAL DOMAIN + (Specified by user via Nudge_Bwindow_File) Default: 0 @@ -355,6 +411,9 @@ 0 == OFF (No Nudging of this variable) 1 == CONSTANT (Spatially Uniform Nudging) 2 == HEAVISIDE WINDOW FUNCTION + 3 == HEAVISIDE WINDOW FUNCTION IN VERTICAL + and USER SPECIFIED HORIZONAL DOMAIN + (Specified by user via Nudge_Bwindow_File) Default: 0 @@ -367,10 +426,12 @@ - Profile index for PS nudging. + Profile index for PS nudging.(Horizontal Only) 0 == OFF (No Nudging of this variable) 1 == CONSTANT (Spatially Uniform Nudging) 2 == HEAVISIDE WINDOW FUNCTION + 3 == USER SPECIFIED HORIZONAL DOMAIN + (Specified by user via Nudge_Bwindow_File) Default: 0 @@ -381,6 +442,14 @@ Default: 0. + + Coeffcient controling the vertical influence function. + Nudge_PSscal < 1.0 : P'(n) = [ B(n) ]*(Target_Ps - Model_Ps) + Nudge_PSscal > 1.0 : P'(n) = [(Pref/P0)**scl]*(Target_Ps - Model_Ps) + Default: 0. + + LAT0 center of Horizontal Window in degrees [-90.,90.]. @@ -461,6 +530,13 @@ Default: FALSE + + Filename for the Horizonal Boundary Window to use for nudging with profile option 3. + (e.g. '/$DIN_LOC_ROOT/atm/cam/nudging/BoundaryWindowFile.nc') + Default: 'BoundaryWindowFile.nc' + + Switch to turn on zonal mean filtering nudging. If TRUE, the nudging scheme diff --git a/cime_config/testdefs/testmods_dirs/cam/nudging/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/nudging/user_nl_cam index ce798ca005..081248b587 100644 --- a/cime_config/testdefs/testmods_dirs/cam/nudging/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/nudging/user_nl_cam @@ -20,6 +20,11 @@ Nudge_TimeScale_Opt = 0 Nudge_Times_Per_Day=4 Model_Times_Per_Day=48 + Nudge_ZonalFilter=.false. + Nudge_ZonalNbasis=-1 + Nudge_SpectralFilter=.false. + Nudge_SpectralNtrunc=-1 + Nudge_SpectralNring= 1 Nudge_Uprof =1 Nudge_Ucoef =1.00 Nudge_Vprof =1 @@ -30,6 +35,7 @@ Nudge_Qcoef =0.00 Nudge_PSprof =0 Nudge_PScoef =0.00 + Nudge_PSscal =0.00 Nudge_Beg_Year =2008 Nudge_Beg_Month=12 Nudge_Beg_Day =16 @@ -48,5 +54,6 @@ Nudge_Vwin_Lindex =5. Nudge_Vwin_Ldelta =1.0 Nudge_Vwin_Invert =.false. + Nudge_Bwindow_File='BoundaryWindowFile.nc' / diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_f10_L26/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_f10_L26/user_nl_cam index 05a64cd2a2..9726ff6f94 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_f10_L26/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_f10_L26/user_nl_cam @@ -13,6 +13,9 @@ Nudge_TimeScale_Opt = 0 Nudge_Times_Per_Day=4 Model_Times_Per_Day=48 + Nudge_SpectralFilter=.false. + Nudge_SpectralNtrunc=-1 + Nudge_SpectralNring= 1 Nudge_Uprof =1 Nudge_Ucoef =1.00 Nudge_Vprof =1 @@ -21,8 +24,9 @@ Nudge_Tcoef =1.00 Nudge_Qprof =1 Nudge_Qcoef =1.00 - Nudge_PSprof =1 - Nudge_PScoef =1.00 + Nudge_PSprof =0 + Nudge_PScoef =0.00 + Nudge_PSscal =0.00 Nudge_Beg_Year =0001 Nudge_Beg_Month=01 Nudge_Beg_Day =01 diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_ne5_L26/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_ne5_L26/user_nl_cam index 4b17143322..44cc74f9d5 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_ne5_L26/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_ne5_L26/user_nl_cam @@ -13,6 +13,9 @@ Nudge_TimeScale_Opt = 0 Nudge_Times_Per_Day=4 Model_Times_Per_Day=48 + Nudge_SpectralFilter=.false. + Nudge_SpectralNtrunc=-1 + Nudge_SpectralNring= 1 Nudge_Uprof =1 Nudge_Ucoef =1.00 Nudge_Vprof =1 @@ -21,8 +24,9 @@ Nudge_Tcoef =1.00 Nudge_Qprof =1 Nudge_Qcoef =1.00 - Nudge_PSprof =1 - Nudge_PScoef =1.00 + Nudge_PSprof =0 + Nudge_PScoef =0.00 + Nudge_PSscal =0.00 Nudge_Beg_Year =0001 Nudge_Beg_Month=01 Nudge_Beg_Day =01 diff --git a/doc/ChangeLog b/doc/ChangeLog index 4ef0649b98..aea1f981bc 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,101 @@ +=============================================================== + +Tag name: cam6_4_191 +Originator(s): patc, jet, fvitt +Date: Aug 3, 2026 +One-line Add Spectral Nudging Option and Fix Nudging Behavior When Files Are Missing- #1402 +Github PR URL: https://github.com/ESCOMP/CAM/pull/1402 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): +This pull request addresses issues #1393 and #1273 + +For Issue #1273 an error exit with a message is invoked if a missing nudging file is encountered. + +For Issue #1393 a module for computing spherical harmonics on the physics grid is implemented. The nudging module has additional functionality that uses this module to optionally filter nudging tendencies to scales larger then a specified truncation scale. + +This functionality is controlled by 2 additional variables in the &nudging_nl namelist. +Nudge_SpectralFilter- LOGICAL Option to apply spherical harminic filtering to +the model state and target data so that nudging +tendencies are only applied to scales larger than +the specified truncation. + +Nudge_SpectralNtrunc - INT The number of meridional spherical harmonic modes used +for spectral filtering. The nominal horizontal scale of +the filtering can be estimated as: + Hscale = PI*6350/Nudge_SpectralNtrunc + i.e. Nudge_SpectralNtrunc=40 corresponds to a horizontal nudging scale Hscale~500km. + +New Functionality was also added for PS nudging and user defined nudging windows. + - PS nudging is now available via namelist variables nudge_pscoef, nudge_psprof, + and nudge_psscal + - User-defined nudging window can be specified via a netcdf file defining the + horizontal boundaries and using a profile index of 3 for nudge_uprof, + nudge_vprof, nudge_tprof, nudge_qprof or nudge_psprof + +Describe any changes made to build system: N/A + +Describe any changes made to the namelist: + + New Namelist variables added - see nudging.F90 for more information. + - nudge_bwindow_file - CHAR path to the Horizonal Boundary Window file. + - nudge_psscal - REAL Coeffcient controling the vertical influence function. + - nudge_spectralfilter - LOGICAL Option to apply spherical harminic filtering to model state + - nudge_spectralnring - INT number of sampling rings used for local area avg of spherical harmonic modes + - nudge_spectralntrunc - INT number of meridional spherical harmonic modes + - nudge_zonalfilter - LOGICAL apply zonal mean filtering + - nudge_zonalnbasis - INT number of meridional modes(Legendre Polynomials) used for filtering + + New Namelist Defaults + - nudge_pscoef 1.00 => 0.00 + - nudge_psprof 1 => 0 + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: jet, fvitt + +List all files eliminated: +D src/utils/zonal_mean_mod.F90 + +List all files added and what they do: +A src/utils/ug_spectralmethods_mod.F90 + - Refactored spectral methods which added spectral functionality and replace zonal_mean_mod.F90 + +List all existing files that have been modified, and describe the changes: +M bld/namelist_files/namelist_definition.xml + - new spectral namelist variables +M cime_config/testdefs/testmods_dirs/cam/nudging/user_nl_cam +M cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_f10_L26/user_nl_cam +M cime_config/testdefs/testmods_dirs/cam/outfrq3s_nudging_ne5_L26/user_nl_cam + - new user_nl_cam files for spectral nudge testing +M src/physics/cam/nudging.F90 + - refactor of nudging routine. +M src/physics/cam/phys_grid_ctem.F90 + - new location for zonal_mean procedures, moved from zonal_mean_mod.F90 +M tools/nudging/user_nl_cam-NUDGING_TEMPLATE + - nudgeting template +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +derecho/intel/aux_cam: + FAIL SMS_Ln9.f09_f09_mg17.F2010climo.derecho_intel.cam-nudging NLCOMP + - expected due to addition of new nl variables + +derecho/nvhpc/aux_cam: all BFB + +izumi/nag/aux_cam: + FAIL ERC_D_Ln9.f10_f10_mt232.FHIST_C5.izumi_nag.cam-outfrq3s_subcol COMPARE_base_rest + - pre-existing failure -- see https://github.com/ESCOMP/CAM/issues/1514 + +izumi/gnu/aux_cam: + NLFAIL ERC_D_Ln9.ne5_ne5_mg37.QPC4.izumi_gnu.cam-outfrq3s_nudging_ne5_L26 NLCOMP + - expected due to addition of new nl variables + +=============================================================== =============================================================== Tag name: cam6_4_190 @@ -1101,8 +1198,8 @@ but without RTTOV output), RTTOV will not be compiled. Describe any changes made to the namelist: Trigger use of COSP-RTTOV with cosp_lrttov_sim, cosp_rttov_Ninstruments, cosp_rttov_instrument_namelists -Allows COSP output to apply simple swathing using “cosp_N_SWATHS_[INSTRUMENT_NAME]”, “cosp_SWATH_LOCALTIMES_[INSTRUMENT_NAME]”, -“cosp_SWATH_WIDTHS_[INSTRUMENT_NAME]" +Allows COSP output to apply simple swathing using ?cosp_N_SWATHS_[INSTRUMENT_NAME]?, ?cosp_SWATH_LOCALTIMES_[INSTRUMENT_NAME]?, +?cosp_SWATH_WIDTHS_[INSTRUMENT_NAME]" List any changes to the defaults for the boundary datasets: N/A @@ -13312,7 +13409,7 @@ Github PR URL: https://github.com/ESCOMP/CAM/pull/1117 Purpose of changes (include the issue number and title text for each relevant GitHub issue): - Continuing development of gravity wave parameterization (//github.com/ESCOMP/CAM/issues/1115) - Added vorticity calculation to SE dycore. Vorticity is passed to the gravity wave (GW) scheme in - model physics to provide a possible source for ‘moving mountain’ GW, i.e., low-phase speed GW forced + model physics to provide a possible source for ?moving mountain? GW, i.e., low-phase speed GW forced by atmospheric circulations. This provides another forcing option, in addition to boundary layer momentum flux implemented earlier. Vorticity anomalies as sources for GW have been proposed by other researchers in published papers. @@ -35744,7 +35841,7 @@ Purpose of changes (include the issue number and title text for each relevant Gi - add value of zmconv_parcel_pbl in the namelist #498 - Update externals to cesm2_3_beta07 #499 - default pecount for ne120xx needs to be increased to run "out of the box" #500 - - Update dyn_grid to read from the INI grid if the initial conditions file has dimension ‘ncol.’ Always write using the GLL grid. (no issue opened for this bug) + - Update dyn_grid to read from the INI grid if the initial conditions file has dimension ?ncol.? Always write using the GLL grid. (no issue opened for this bug) Describe any changes made to build system: @@ -62420,7 +62517,7 @@ yellowstone/intel/aux_cam: SMS_D_Ln9.f19_f19.FXHIST.yellowstone_intel.cam-outfrq9s (Overall: FAIL) details: FAIL SMS_D_Ln9.f19_f19.FXHIST.yellowstone_intel.cam-outfrq9s NLCOMP FAIL SMS_D_Ln9.f19_f19.FXHIST.yellowstone_intel.cam-outfrq9s RUN time=118 - SMS_Ld1.f09_f09.FW1850.yellowstone_intel.cam-reduced_hist1d (Overall: NLFAI≈ßL) details: + SMS_Ld1.f09_f09.FW1850.yellowstone_intel.cam-reduced_hist1d (Overall: NLFAI?L) details: FAIL SMS_Ld1.f09_f09.FW1850.yellowstone_intel.cam-reduced_hist1d NLCOMP Run failures for SPCAM and WACCMX are pre-existing due to CLM version @@ -64087,8 +64184,8 @@ Describe any changes made to build system: Describe any changes made to the namelist: - Added the following namelist variables for gw_rdg gw_rdg_do_divstream -- If .true. use separate dividing streamlines for downslope wind and - flow splitting regimes (“DS” configuration). If .false. use single - dividing streamline as in Scinocca&McFarlane 2000 (“SM” configuration). + flow splitting regimes (?DS? configuration). If .false. use single + dividing streamline as in Scinocca&McFarlane 2000 (?SM? configuration). gw_rdg_C_BetaMax_DS -- Enhancement factor for downslope wind stress in DS configuration gw_rdg_C_GammaMax -- Enhancement factor for depth of downslope wind regime in DS configuration gw_rdg_Frx0 -- Lower inverse Froude number limits on linear ramp terminating downslope wind @@ -65069,7 +65166,7 @@ Purpose of changes: Modify the chemical tropopause definition poleward of 50 degrees from using a fixed 300 hPa level (see note for 5_4_89) to use the lapse rate tropopause at P > 125 hPa and the climatology if the lapse rate tropopause is found at - P < 125 hPa or is not found. This eliminates “false” tropopause detection in the + P < 125 hPa or is not found. This eliminates ?false? tropopause detection in the 75 hPa < P < 125 hPa range that are not physical and have a detrimental effect on heterogeneous chemistry in the polar stratosphere, provides a diagnostic that can respond to some climate change, and gives a smoother and more physical tropopause @@ -89667,7 +89764,7 @@ M models/atm/cam/test/system/tests_chem_mpi - updates to chemistry regression tests M models/atm/cam/src/chemistry/utils/tracer_data.F90 - -- Bugzilla – Bug 1824 -- don't use cyc_ndx_beg before it is defined + -- Bugzilla ? Bug 1824 -- don't use cyc_ndx_beg before it is defined If there were any failures reported from running test_driver.sh on any test platform, and checkin with these failures has been OK'd by the gatekeeper, diff --git a/src/physics/cam/nudging.F90 b/src/physics/cam/nudging.F90 index ced2ef57d2..4eeed18a02 100644 --- a/src/physics/cam/nudging.F90 +++ b/src/physics/cam/nudging.F90 @@ -19,13 +19,11 @@ module nudging ! forcing discontinues. ! ! Some analyses products can have gaps in the available data, where values -! are missing for some interval of time. When files are missing, the nudging -! force is switched off for that interval of time, so we effectively 'coast' -! thru the gap. -! -! Currently, the nudging module is set up to accomodate nudging of PS -! values, however that functionality requires forcing that is applied in -! the selected dycore and is not yet implemented. +! are missing for some interval of time. The default behavior is now for the +! model to error exit if there is a gap. Users with known gaps in their nuding +! data can manually change the gap behavior to accomodate their needs. +! When files are missing, the nudging force can be switched off for that interval +! of time, in order to effectively 'coast' thru the gap. ! ! The nudging of the model toward the analyses data is controlled by ! the 'nudging_nl' namelist in 'user_nl_cam'; whose variables control the @@ -43,13 +41,41 @@ module nudging ! ! F_nudge = Alpha*((Target-Model(t_curr))/TimeScale ! +! SURFACE PRESSURE NUDGING: +! ------------------------ +! Surface pressure is not a prognostic variable, so nudging cannot be implemented directly. +! Hydrostatic balance provides a constraint on the values of pressure, temperature, and model +! layer thicknesses(geopotential heights). Since any change in thickness values will translate +! into a momentum tendency via the gradient of height, we require that dZ/dn must remain constant +! under changes to pressures so that the result does not indirectly impose a momentum nudging. +! The nudging of surface pressure values then has a corresponding nudging of temperature values. +! Surface pressure nudging is just an alternate form of temperature nudging, thus both temperature +! and surface pressure cannot be applied at the same time. ! +! Vertical Influence Function: +! To convey the effect of the 2D surface pressure changes to the 3D temperature values, there are +! two choices for the vertical influence function. The first, most natural choice, is to just +! convey the surface pressures using the (Bn) hybrid pressure coeffcients. However, due to the +! inherent structure of the hybrid coeffcients this results in a vertical profile that is small +! at the surface increases to a maximum at the upper levels, and then abruptly drops to 0 for the +! pure-pressure levels at the mode top. Very un-physical behavior, but that is what the mathematics +! dictates for the current (An,Bn) values. An alternative vertical influence function is implemented +! based on the reference profile of the atmosphere. +! (Pref/P0)**scl = exp[-scl*(Ztop/H)*(1-n)]. (H=7km) +! For scl=1.0, the vertical profile is similar to the Bn profile, but with a smaller amplitude. +! For larger values of scl, the profile transitions to one in which the vertical influence function +! falls off with height from a maximum at the surface layer. There is typically a discontinuity in +! the profile at the transition to pure pressure levels, which gets smaller with increasing scl values. +! ! WINDOWING: ! ---------- ! The region of applied nudging can be limited using Horizontal/Vertical ! window functions that are constructed using a parameterization of the ! Heaviside step function. ! +! ------------------------------------------------------------------------------- +! HEAVISIDE WINDOW: (Nudge_Xprof = 2) where X=[U,V,T,Q,PS] +! ------------------------------------------------------------------------------- ! The Heaviside window function is the product of separate horizonal and vertical ! windows that are controled via 12 parameters: ! @@ -93,16 +119,43 @@ module nudging ! running the model. Lookat_NudgeWindow.ncl is a script avalable in the tools directory ! which will read in the values for a given namelist and display the resulting window. ! -! The module is currently configured for only 1 window function. It can readily be -! extended for multiple windows if the need arises. +! ------------------------------------------------------------------------------- +! USER-DEFINED WINDOW: (Nudge_Xprof = 3) where X=[U,V,T,Q,PS] +! ------------------------------------------------------------------------------- +! The horizontal window domain can be customized by providing a netCDF file containing +! a lat/lon grid of window coeffcients [0.,1.], where 1.0 represents the nudged domain +! and the values taper to 0. for un-nudged gridpoints. The Horizonal window values at +! model gridpoints are linearly interpolated from the given reclilinear grid. For this +! case the user specifies the Nudge_Bwindow_File file in the namelist and setting +! (Nudge_Hwin_Invert=.true.) will invert the window given in the file, but all of the +! other Nudge_Hwin_* namelist values are ignored. The vertical domain of the window is +! still specified via the Nudge_Vwin_* namelist values. +! +! The format for the Nudge_Bwindow_File can be either: +! +! dimensions: ! dimensions: +! latitude; ! latitude; +! longitude; ! longitude; +! variables: ! variables: +! double latitude(latitude); ! double latitude(latitude); +! double longitude(longitude); ! double longitude(longitude); +! double refineMap(latitude, longitude); ! double boundaryMap(latitude, longitude); ! +! The left format is the REFMAP file output that is created from the variable mesh VRM editor +! program. The default resolution for these files is [720,360]. There is no restriction on the +! horizonal dimensions, so higher resoultion domains can be provided as lonag as the grid points +! span the lat/lon domain. +! ! -! Input/Output Values: +! INPUT/OUTPUT VALUES: +! --------------------- ! Forcing contributions are available for history file output by -! the names: {'Nudge_U','Nudge_V','Nudge_T',and 'Nudge_Q'} +! the names: {'Nudge_U','Nudge_V','Nudge_T','Nudge_Q', and 'Nudge_PS'} ! The target values that the model state is nudged toward are available for history -! file output via the variables: {'Target_U','Target_V','Target_T',and 'Target_Q'} +! file output via the variables: {'Target_U','Target_V','Target_T','Target_Q',and 'Target_PS'} ! +! NAMELIST SPECIFICATION: +! --------------------- ! &nudging_nl ! Nudge_Model - LOGICAL toggle to activate nudging. ! TRUE -> Nudging is on. @@ -149,26 +202,69 @@ module nudging ! 0 --> TimeScale = 1/Tdlt_Anal [DEFAULT] ! 1 --> TimeScale = 1/(t'_next - t_curr ) ! -! Nudge_Uprof - INT index of profile structure to use for U. [0,1,2] -! Nudge_Vprof - INT index of profile structure to use for V. [0,1,2] -! Nudge_Tprof - INT index of profile structure to use for T. [0,1,2] -! Nudge_Qprof - INT index of profile structure to use for Q. [0,1,2] -! Nudge_PSprof - INT index of profile structure to use for PS. [0,N/A] +! Nudge_ZonalFilter - LOGICAL Option to apply zonal mean filtering to the +! model state and target data. +! +! Nudge_ZonalNbasis - INT The number of meridional modes(Legendre Polynomials) +! used for zonal filtering. +! +! Nudge_SpectralFilter - LOGICAL Option to apply spherical harminic filtering to +! the model state and target data so that nudging +! tendencies are only applied to scales larger than +! the specified truncation. +! +! Nudge_SpectralNtrunc - INT The number of meridional spherical harmonic modes used +! for spectral filtering. The nominal horizontal scale (km) of +! the filtering can be estimated as: +! +! Hscale = PI*6350/Nudge_SpectralNtrunc +! +! i.e. Nudge_SpectralNtrunc=40 corresponds to a horizontal +! nudging scale Hscale~500km. +! +! Nudge_SpectralNring - INT The number of sampling rings used for local area averaging +! of spherical harmonic modes, to suppress sampling errors. +! When initializing each basis, a local average of SH values +! is computed for the area associated with each grid point. +! SpectralNring set the number of rings of equal-area points +! in this sampling domain. +! Each ring (kk) contains 8*(kk-1) sample points. +! +! Nudge_SpectralNring Number of Samping Points +! ------------------- ------------------------- +! 1 1 (DEFAULT SampleGrid NOT used) +! 2 9 +! 3 25 +! 4 49 +! +! Nudge_Uprof - INT index of profile structure to use for U. [0,1,2,3] +! Nudge_Vprof - INT index of profile structure to use for V. [0,1,2,3] +! Nudge_Tprof - INT index of profile structure to use for T. [0,1,2,3] +! Nudge_Qprof - INT index of profile structure to use for Q. [0,1,2,3] ! ! The spatial distribution is specified with a profile index. ! Where: 0 == OFF (No Nudging of this variable) ! 1 == CONSTANT (Spatially Uniform Nudging) ! 2 == HEAVISIDE WINDOW FUNCTION +! 3 == HEAVISIDE WINDOW FUNCTION IN VERTICAL +! and USER SPECIFIED HORIZONAL DOMAIN +! (Specified by user via Nudge_Bwindow_File) ! ! Nudge_Ucoef - REAL fractional nudging coeffcient for U. ! Nudge_Vcoef - REAL fractional nudging coeffcient for V. ! Nudge_Tcoef - REAL fractional nudging coeffcient for T. ! Nudge_Qcoef - REAL fractional nudging coeffcient for Q. -! Nudge_PScoef - REAL fractional nudging coeffcient for PS. ! ! The strength of the nudging is specified as a fractional ! coeffcient between [0,1]. ! +! Nudge_PSprof - INT index of (Horizontal Only) window structure to use for PS. [0,1,2,3] +! Nudge_PScoef - REAL fractional nudging coeffcient for PS. +! Nudge_PSscal - REAL Coeffcient controling the vertical influence function. +! for: Nudge_PSscal < 1.0 : P'(n) = [ B(n) ]*(Target_Ps - Model_Ps) +! Nudge_PSscal > 1.0 : P'(n) = [(Pref/P0)**scl]*(Target_Ps - Model_Ps) +! (default: Nudge_PSscal = 0.0) +! ! Nudge_Hwin_lat0 - REAL latitudinal center of window in degrees. ! Nudge_Hwin_lon0 - REAL longitudinal center of window in degrees. ! Nudge_Hwin_latWidth - REAL latitudinal width of window in degrees. @@ -183,14 +279,10 @@ module nudging ! Nudge_Vwin_Hdelta - REAL HI transition length ! Nudge_Vwin_Invert - LOGICAL FALSE= value=1 inside the specified window, 0 outside ! TRUE = value=0 inside the specified window, 1 outside +! Nudge_Bwindow_File - CHAR path to the Horizonal Boundary Window file. +! (e.g. '/glade/scratch/USER/inputdata/nudging/BoundaryWindowFile.nc') ! / ! -!================ -! -! TO DO: -! ----------- -! ** Implement Ps Nudging???? -! !===================================================================== ! Useful modules !------------------ @@ -201,7 +293,7 @@ module nudging use spmd_utils, only: masterproc, mstrid=>masterprocid, mpicom, mpi_success use spmd_utils, only: mpi_integer, mpi_real8, mpi_logical, mpi_character use cam_logfile, only: iulog - use zonal_mean_mod, only: ZonalMean_t + use ug_spectralmethods_mod, only: SphericalHarmonic_GS_t, ZonalMean_t ! Set all Global values and routines to private by default ! and then explicitly set their exposure. @@ -215,10 +307,12 @@ module nudging public :: nudging_timestep_init public :: nudging_timestep_tend private :: nudging_update_analyses - private :: nudging_set_PSprofile private :: nudging_set_profile - private :: calc_DryStaticEnergy + private :: nudging_set_PSprofile + private :: nudging_set_Vwindow public :: nudging_final + private :: calc_DryStaticEnergy + private :: interp_Bwin ! Nudging Parameters !-------------------- @@ -237,6 +331,7 @@ module nudging real(r8) :: Nudge_Qcoef,Nudge_Tcoef integer :: Nudge_Qprof,Nudge_Tprof real(r8) :: Nudge_PScoef + real(r8) :: Nudge_PSscal integer :: Nudge_PSprof integer :: Nudge_Beg_Year ,Nudge_Beg_Month integer :: Nudge_Beg_Day ,Nudge_Beg_Sec @@ -272,15 +367,24 @@ module nudging real(r8) :: Nudge_Hwin_lonWidthH real(r8) :: Nudge_Hwin_max real(r8) :: Nudge_Hwin_min + character(len=cl) :: Nudge_Bwindow_File - ! Nudging Zonal Filter variables - !--------------------------------- + ! Nudging Zonal/Spectral Filter variables + !----------------------------------------- logical :: Nudge_ZonalFilter =.false. integer :: Nudge_ZonalNbasis = -1 type(ZonalMean_t) :: ZM real(r8),allocatable:: Zonal_Bamp2d(:) real(r8),allocatable:: Zonal_Bamp3d(:,:) + logical :: Nudge_SpectralFilter =.false. + integer :: Nudge_SpectralNtrunc = -1 + integer :: Nudge_SpectralNbasis = -1 + integer :: Nudge_SpectralNring = 1 + type(SphericalHarmonic_GS_t):: SH + real(r8),allocatable:: Spectral_Bamp2d(:) + real(r8),allocatable:: Spectral_Bamp3d(:,:) + ! Nudging State Arrays !----------------------- integer :: Nudge_nlon,Nudge_nlat,Nudge_ncol,Nudge_nlev @@ -295,7 +399,8 @@ module nudging real(r8),allocatable:: Model_T (:,:,:) !(pcols,pver,begchunk:endchunk) real(r8),allocatable:: Model_S (:,:,:) !(pcols,pver,begchunk:endchunk) real(r8),allocatable:: Model_Q (:,:,:) !(pcols,pver,begchunk:endchunk) - real(r8),allocatable:: Model_PS (:,:) !(pcols,begchunk:endchunk) + real(r8),allocatable:: Model_PS (:,:) !(pcols,begchunk:endchunk) + real(r8),allocatable:: Model_PSfilt(:,:) !(pcols,begchunk:endchunk) real(r8),allocatable:: Nudge_Utau (:,:,:) !(pcols,pver,begchunk:endchunk) real(r8),allocatable:: Nudge_Vtau (:,:,:) !(pcols,pver,begchunk:endchunk) real(r8),allocatable:: Nudge_Stau (:,:,:) !(pcols,pver,begchunk:endchunk) @@ -343,11 +448,13 @@ subroutine nudging_readnl(nlfile) Nudge_File_Template, Nudge_Force_Opt, & Nudge_TimeScale_Opt, & Nudge_Times_Per_Day, Model_Times_Per_Day, & + Nudge_SpectralFilter, Nudge_SpectralNtrunc, & + Nudge_SpectralNring, & Nudge_Ucoef , Nudge_Uprof, & Nudge_Vcoef , Nudge_Vprof, & Nudge_Qcoef , Nudge_Qprof, & Nudge_Tcoef , Nudge_Tprof, & - Nudge_PScoef, Nudge_PSprof, & + Nudge_PScoef, Nudge_PSprof, Nudge_PSscal, & Nudge_Beg_Year, Nudge_Beg_Month, Nudge_Beg_Day, & Nudge_End_Year, Nudge_End_Month, Nudge_End_Day, & Nudge_Hwin_lat0, Nudge_Hwin_lon0, & @@ -356,7 +463,7 @@ subroutine nudging_readnl(nlfile) Nudge_Hwin_Invert, & Nudge_Vwin_Lindex, Nudge_Vwin_Hindex, & Nudge_Vwin_Ldelta, Nudge_Vwin_Hdelta, & - Nudge_Vwin_Invert + Nudge_Vwin_Invert, Nudge_Bwindow_File ! For Zonal Mean Filtering namelist /nudging_nl/ Nudge_ZonalFilter, Nudge_ZonalNbasis @@ -383,11 +490,12 @@ subroutine nudging_readnl(nlfile) Nudge_Vcoef = 0._r8 Nudge_Qcoef = 0._r8 Nudge_Tcoef = 0._r8 - Nudge_PScoef = 0._r8 Nudge_Uprof = 0 Nudge_Vprof = 0 Nudge_Qprof = 0 Nudge_Tprof = 0 + Nudge_PScoef = 0._r8 + Nudge_PSscal = 0._r8 Nudge_PSprof = 0 Nudge_Beg_Year = 2008 Nudge_Beg_Month = 5 @@ -411,6 +519,13 @@ subroutine nudging_readnl(nlfile) Nudge_Vwin_Invert = .false. Nudge_Vwin_lo = 0.0_r8 Nudge_Vwin_hi = 1.0_r8 + Nudge_Bwindow_File = 'BoundaryWindowFile.nc' + Nudge_ZonalFilter =.false. + Nudge_ZonalNbasis = -1 + Nudge_SpectralFilter =.false. + Nudge_SpectralNtrunc = -1 + Nudge_SpectralNbasis = -1 + Nudge_SpectralNring = 1 ! Read in namelist values !------------------------ @@ -458,7 +573,7 @@ subroutine nudging_readnl(nlfile) call endrun('nudging_readnl:: ERROR in namelist') endif - if((Nudge_Vwin_Lindex > Nudge_Vwin_Hindex) .or. & + if((Nudge_Vwin_Lindex > Nudge_Vwin_Hindex) .or. & (Nudge_Vwin_Hindex > float(pver+1)) .or. (Nudge_Vwin_Hindex < 0._r8) .or. & (Nudge_Vwin_Lindex > float(pver+1)) .or. (Nudge_Vwin_Lindex < 0._r8) ) then write(iulog,*) 'NUDGING: Window Lindex must be in [0,pver+1]' @@ -487,6 +602,27 @@ subroutine nudging_readnl(nlfile) call endrun('nudging_readnl:: ERROR in namelist') endif + if((Nudge_ZonalFilter).and.(Nudge_SpectralFilter)) then + write(iulog,*) 'NUDGING: Zonal Nudging and Spectral Nudging cannot both be active' + write(iulog,*) 'NUDGING: Nudge_SpectralFilter=',Nudge_SpectralFilter + write(iulog,*) 'NUDGING: Nudge_ZonalFilter=',Nudge_ZonalFilter + call endrun('nudging_readnl:: ERROR in namelist') + endif + + if((Nudge_ZonalFilter).and.(Nudge_ZonalNbasis.le.0)) then + write(iulog,*) 'NUDGING: Zonal Nudging requires that (Nudge_ZonalNbasis > 0)' + write(iulog,*) 'NUDGING: Nudge_ZonalFilter=',Nudge_ZonalFilter + write(iulog,*) 'NUDGING: Nudge_ZonalNbasis=',Nudge_ZonalNbasis + call endrun('nudging_readnl:: ERROR in namelist') + endif + + if((Nudge_SpectralFilter).and.(Nudge_SpectralNtrunc.le.0)) then + write(iulog,*) 'NUDGING: Spectral Nudging requires that (Nudge_SpectralNtrunc > 0)' + write(iulog,*) 'NUDGING: Nudge_SpectralFilter=',Nudge_SpectralFilter + write(iulog,*) 'NUDGING: Nudge_SpectralNtrunc=',Nudge_SpectralNtrunc + call endrun('nudging_readnl:: ERROR in namelist') + endif + ! Broadcast namelist variables !------------------------------ call MPI_bcast(Nudge_Path , len(Nudge_Path), & @@ -521,6 +657,8 @@ subroutine nudging_readnl(nlfile) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Qcoef') call MPI_bcast(Nudge_PScoef , 1, mpi_real8 , mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_PScoef') + call MPI_bcast(Nudge_PSscal , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_PSscal') call MPI_bcast(Nudge_Uprof , 1, mpi_integer, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Uprof') call MPI_bcast(Nudge_Vprof , 1, mpi_integer, mstrid, mpicom, ierr) @@ -579,10 +717,18 @@ subroutine nudging_readnl(nlfile) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Vwin_Ldelta') call MPI_bcast(Nudge_Vwin_Invert, 1, mpi_logical, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Vwin_Invert') + call MPI_bcast(Nudge_Bwindow_File, len(Nudge_Bwindow_File), mpi_character, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Bwindow_File') call MPI_bcast(Nudge_ZonalFilter, 1, mpi_logical, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_ZonalFilter') call MPI_bcast(Nudge_ZonalNbasis, 1, mpi_integer, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_ZonalNbasis') + call MPI_bcast(Nudge_SpectralFilter, 1, mpi_logical, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_SpectralFilter') + call MPI_bcast(Nudge_SpectralNtrunc, 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_SpectralNtrunc') + call MPI_bcast(Nudge_SpectralNring, 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_SpectralNring') ! End Routine !------------ @@ -601,9 +747,10 @@ subroutine nudging_init use dycore ,only: dycore_is use dyn_grid ,only: get_horiz_grid_dim_d use phys_grid ,only: get_rlat_p,get_rlon_p,get_ncols_p - use cam_history ,only: addfld + use cam_history ,only: addfld,horiz_only use shr_const_mod ,only: SHR_CONST_PI use filenames ,only: interpret_filename_spec + use netcdf ! Local values !---------------- @@ -622,6 +769,13 @@ subroutine nudging_init real(r8) :: Val1_n,Val2_n,Val3_n,Val4_n integer :: nn + real(r8),allocatable:: Bwindow(:,:) + real(r8),allocatable:: B_lon (:) + real(r8),allocatable:: B_lat (:) + integer :: Bnlon,Bnlat,ncid,varid + real(r8):: Hcoef + real(r8):: Vwindow(pver) + character(len=*), parameter :: prefix = 'nudging_init: ' ! Get the time step size @@ -655,6 +809,8 @@ subroutine nudging_init call alloc_err(istat,'nudging_init','Model_Q',pcols*pver*((endchunk-begchunk)+1)) allocate(Model_PS(pcols,begchunk:endchunk),stat=istat) call alloc_err(istat,'nudging_init','Model_PS',pcols*((endchunk-begchunk)+1)) + allocate(Model_PSfilt(pcols,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_init','Model_PSfilt',pcols*((endchunk-begchunk)+1)) ! Allocate Space for spatial dependence of ! Nudging Coefs and Nudging Forcing. @@ -690,7 +846,10 @@ subroutine nudging_init call addfld('Target_U',(/ 'lev' /),'A','m/s' ,'U Nudging Target' ) call addfld('Target_V',(/ 'lev' /),'A','m/s' ,'V Nudging Target' ) call addfld('Target_T',(/ 'lev' /),'A','K' ,'T Nudging Target' ) - call addfld('Target_Q',(/ 'lev' /),'A','kg/kg' ,'Q Nudging Target ') + call addfld('Target_Q',(/ 'lev' /),'A','kg/kg' ,'Q Nudging Target' ) + call addfld( 'Nudge_PS',horiz_only,'A','hPa/s' ,'PS Nudging Tendency') + call addfld('Target_PS',horiz_only,'A','hPa' ,'PS Nudging Target' ) + call addfld( 'Model_PS',horiz_only,'A','hPa' ,'PS Model Surface P' ) !----------------------------------------- ! Values initialized only by masterproc @@ -848,11 +1007,15 @@ subroutine nudging_init write(iulog,*) 'NUDGING: Model_Step=',Model_Step write(iulog,*) 'NUDGING: Nudge_ZonalFilter=',Nudge_ZonalFilter write(iulog,*) 'NUDGING: Nudge_ZonalNbasis=',Nudge_ZonalNbasis + write(iulog,*) 'NUDGING: Nudge_SpectralFilter=',Nudge_SpectralFilter + write(iulog,*) 'NUDGING: Nudge_SpectralNtrunc=',Nudge_SpectralNtrunc + write(iulog,*) 'NUDGING: Nudge_SpectralNring=',Nudge_SpectralNring write(iulog,*) 'NUDGING: Nudge_Ucoef =',Nudge_Ucoef write(iulog,*) 'NUDGING: Nudge_Vcoef =',Nudge_Vcoef write(iulog,*) 'NUDGING: Nudge_Qcoef =',Nudge_Qcoef write(iulog,*) 'NUDGING: Nudge_Tcoef =',Nudge_Tcoef write(iulog,*) 'NUDGING: Nudge_PScoef =',Nudge_PScoef + write(iulog,*) 'NUDGING: Nudge_PSscal =',Nudge_PSscal write(iulog,*) 'NUDGING: Nudge_Uprof =',Nudge_Uprof write(iulog,*) 'NUDGING: Nudge_Vprof =',Nudge_Vprof write(iulog,*) 'NUDGING: Nudge_Qprof =',Nudge_Qprof @@ -884,11 +1047,50 @@ subroutine nudging_init write(iulog,*) 'NUDGING: Nudge_Hwin_lonWidthH=',Nudge_Hwin_lonWidthH write(iulog,*) 'NUDGING: Nudge_Hwin_max =',Nudge_Hwin_max write(iulog,*) 'NUDGING: Nudge_Hwin_min =',Nudge_Hwin_min + write(iulog,*) 'NUDGING: Nudge_Bwindow_File =',trim(Nudge_Bwindow_File) write(iulog,*) 'NUDGING: Nudge_Initialized =',Nudge_Initialized write(iulog,*) ' ' write(iulog,*) 'NUDGING: Nudge_NumObs=',Nudge_NumObs write(iulog,*) ' ' + ! Error Check TSmode usage + !--------------------------- + if((Nudge_ZonalFilter).or.(Nudge_SpectralFilter).or.(Nudge_PSprof.ne.0)) then + if(Nudge_TSmode.ne.0) then + write(iulog,*) 'NUDGING: Nudge_TSmode must be set to 0 when nudging Ps or ' + write(iulog,*) 'NUDGING: when spectral filtering is active. Until this module ' + write(iulog,*) 'NUDGING: is restructured to fix it, the computation DSE gives ' + write(iulog,*) 'NUDGING: inconsistent results for the Model and Target states.' + call endrun('NUDGING: Eror with TSmode nudging option') + endif + endif + + ! Error Check and set up for Ps Nudging option + !--------------------------------------------- + if(Nudge_PSprof.ne.0) then + ! Ps Nudging is activated. + ! Make sure that Direct nudging of T is disbled + !------------------------------------------------ + if(Nudge_Tprof.ne.0) then + write(iulog,*) 'NUDGING: ***ERROR STOP*** ' + write(iulog,*) 'NUDGING: Nudge_Tprof and Nudge_PSprof cannot both be non-zero' + write(iulog,*) 'NUDGING: at the same time. The effects of PS nudging are implemented ' + write(iulog,*) 'NUDGING: via temperature tendenies in the near surface layers. ' + write(iulog,*) 'NUDGING: Set Nudge_Tprof=0 if you wish to apply Ps nudging. ' + call endrun('NUDGING: Eror with PS nudging option') + endif + + ! Internally turn on T nudging and set the + ! nudging coef with the Nudge_PScoef value + !------------------------------------------ + Nudge_Tprof = Nudge_PSprof + Nudge_Tcoef = Nudge_PScoef + + write(iulog,*) 'NUDGING: Activating Surface Pressure (Ps) nudging ' + write(iulog,*) 'NUDGING: Nudge_Tprof =',Nudge_Tprof + write(iulog,*) 'NUDGING: Nudge_Tcoef =',Nudge_Tcoef + endif + endif ! (masterproc) then ! Broadcast other variables that have changed @@ -937,6 +1139,10 @@ subroutine nudging_init if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Hwin_latWidthH') call MPI_bcast(Nudge_NumObs , 1, mpi_integer, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_NumObs') + call MPI_bcast(Nudge_Tprof , 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Tprof') + call MPI_bcast(Nudge_Tcoef , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(prefix//'FATAL: mpi_bcast: Nudge_Tcoef') ! All non-masterproc processes also need to allocate space ! before the broadcast of Nudge_NumObs dependent data. @@ -985,8 +1191,8 @@ subroutine nudging_init endif !!DIAG - ! Initialize the Zonal Mean type if needed - !------------------------------------------ + ! Initialize the Zonal Mean Spectral type if needed + !-------------------------------------------------- if(Nudge_ZonalFilter) then call ZM%init(Nudge_ZonalNbasis) allocate(Zonal_Bamp2d(Nudge_ZonalNbasis),stat=istat) @@ -995,6 +1201,18 @@ subroutine nudging_init call alloc_err(istat,'nudging_init','Zonal_Bamp3d',Nudge_ZonalNbasis*pver) endif + if(Nudge_SpectralFilter) then + write(iulog,*) 'NUDGING: calling SH%init() Nudge_SpectralNtrunc =',Nudge_SpectralNtrunc + write(iulog,*) 'NUDGING: calling SH%init() Nudge_SpectralNring =',Nudge_SpectralNring + call SH%init(Nudge_SpectralNtrunc,Nudge_SpectralNbasis,SAMPLE_NRING=Nudge_SpectralNring) + write(iulog,*) 'NUDGING: done SH%init() Nudge_SpectralNbasis =',Nudge_SpectralNbasis + allocate(Spectral_Bamp2d(Nudge_SpectralNbasis),stat=istat) + call alloc_err(istat,'nudging_init','Spectral_Bamp2d',Nudge_SpectralNbasis) + allocate(Spectral_Bamp3d(Nudge_SpectralNbasis,pver),stat=istat) + call alloc_err(istat,'nudging_init','Spectral_Bamp3d',Nudge_SpectralNbasis*pver) + write(iulog,*) 'NUDGING: SH% Arrays allocated' + endif + ! Initialize the analysis filename at the NEXT time for startup. !--------------------------------------------------------------- Nudge_File=interpret_filename_spec(Nudge_File_Template , & @@ -1012,6 +1230,110 @@ subroutine nudging_init !---------------------------------------------------------- call nudging_update_analyses (trim(Nudge_Path)//trim(Nudge_File)) + ! Prepare needed variables if a custom window + ! (Nudge_Xprof==3) is used + !------------------------------------------------- + if((Nudge_Uprof.eq.3).or.(Nudge_Vprof.eq.3).or. & + (Nudge_Tprof.eq.3).or.(Nudge_Qprof.eq.3).or.(Nudge_PSprof.eq.3)) then + + ! Initialize values that specify the + ! custom horizonal Boundary window + !------------------------------------- + istat = nf90_open(trim(Nudge_Bwindow_File),NF90_NOWRITE,ncid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) 'NF90_OPEN: failed for file:',trim(Nudge_Bwindow_File) + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + else + if(masterproc) then + write(iulog,*) 'NUDGING: Opened Boundary Window File:',trim(Nudge_Bwindow_File) + endif + endif + + istat=nf90_inq_dimid(ncid,'longitude',varid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat=nf90_inquire_dimension(ncid,varid,len=Bnlon) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + + istat=nf90_inq_dimid(ncid,'latitude',varid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat=nf90_inquire_dimension(ncid,varid,len=Bnlat) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + + allocate(Bwindow(Bnlon,Bnlat)) + allocate(B_lon (Bnlon)) + allocate(B_lat (Bnlat)) + + istat=nf90_inq_varid(ncid,'longitude',varid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat=nf90_get_var(ncid,varid,B_lon) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + + istat=nf90_inq_varid(ncid,'latitude',varid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat=nf90_get_var(ncid,varid,B_lat) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + + istat=nf90_inq_varid(ncid,'refineMap',varid) + if(istat.ne.NF90_NOERR) then + istat=nf90_inq_varid(ncid,'boundaryMap',varid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat=nf90_get_var(ncid,varid,Bwindow) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + endif + istat=nf90_get_var(ncid,varid,Bwindow) + if(istat.ne.NF90_NOERR) then + write(iulog,*) nf90_strerror(istat) + call endrun('subroutine nudging_init') + endif + istat = nf90_close(ncid) + + ! Check Bwindow values are in range [0,1] + !------------------------------------------ + if((MAXVAL(Bwindow).gt.1._r8).or.(MINVAL(Bwindow).lt.0._r8)) then + write(iulog,*) 'NUDGING: ERROR Boundary Window values are not in the range [0,1]' + write(iulog,*) 'NUDGING: ERROR Boundary Window MAX =',MAXVAL(Bwindow) + write(iulog,*) 'NUDGING: ERROR Boundary Window MIN =',MINVAL(Bwindow) + call endrun('subroutine nudging_init') + endif + + ! Optionally: Invert the Boudnary window + !---------------------------------------- + if(Nudge_Hwin_Invert) then + Bwindow(:,:) = 1._r8 - Bwindow(:,:) + endif + endif + ! Initialize Nudging Coeffcient profiles in local arrays ! Load zeros into nudging arrays !------------------------------------------------------ @@ -1021,17 +1343,81 @@ subroutine nudging_init rlat=get_rlat_p(lchnk,icol)*180._r8/SHR_CONST_PI rlon=get_rlon_p(lchnk,icol)*180._r8/SHR_CONST_PI - call nudging_set_profile(rlat,rlon,Nudge_Uprof,Wprof,pver) - Nudge_Utau(icol,:,lchnk)=Wprof(:) - call nudging_set_profile(rlat,rlon,Nudge_Vprof,Wprof,pver) - Nudge_Vtau(icol,:,lchnk)=Wprof(:) - call nudging_set_profile(rlat,rlon,Nudge_Tprof,Wprof,pver) - Nudge_Stau(icol,:,lchnk)=Wprof(:) - call nudging_set_profile(rlat,rlon,Nudge_Qprof,Wprof,pver) - Nudge_Qtau(icol,:,lchnk)=Wprof(:) + if(Nudge_Uprof.eq.3) then + ! Interp Horizontal window values from user input + ! Initialize the vertical profile from namelist parameters + !------------------------------------------------ + call nudging_set_Vwindow(Vwindow,pver) + Hcoef = interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,Bnlon,Bnlat) + Nudge_Utau(icol,:,lchnk)=Hcoef*Vwindow(:) + else + ! Set prodfile for Nudge_Uprof = [0,1,2] + !--------------------------------------- + call nudging_set_profile(rlat,rlon,Nudge_Uprof,Wprof,pver) + Nudge_Utau(icol,:,lchnk)=Wprof(:) + endif + + if(Nudge_Vprof.eq.3) then + ! Interp Horizontal window values from user input + ! Initialize the vertical profile from namelist parameters + !------------------------------------------------ + call nudging_set_Vwindow(Vwindow,pver) + Hcoef = interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,Bnlon,Bnlat) + Nudge_Vtau(icol,:,lchnk)=Hcoef*Vwindow(:) + else + ! Set prodfile for Nudge_Uprof = [0,1,2] + !--------------------------------------- + call nudging_set_profile(rlat,rlon,Nudge_Vprof,Wprof,pver) + Nudge_Vtau(icol,:,lchnk)=Wprof(:) + endif + + if(Nudge_Tprof.eq.3) then + ! Interp Horizontal window values from user input + ! Initialize the vertical profile from namelist parameters + !------------------------------------------------ + call nudging_set_Vwindow(Vwindow,pver) + Hcoef = interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,Bnlon,Bnlat) + Nudge_Stau(icol,:,lchnk)=Hcoef*Vwindow(:) + else + ! Set prodfile for Nudge_Uprof = [0,1,2] + !--------------------------------------- + call nudging_set_profile(rlat,rlon,Nudge_Tprof,Wprof,pver) + Nudge_Stau(icol,:,lchnk)=Wprof(:) + endif + + if(Nudge_Qprof.eq.3) then + ! Interp Horizontal window values from user input + ! Initialize the vertical profile from namelist parameters + !------------------------------------------------ + call nudging_set_Vwindow(Vwindow,pver) + Hcoef = interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,Bnlon,Bnlat) + Nudge_Qtau(icol,:,lchnk)=Hcoef*Vwindow(:) + else + ! Set prodfile for Nudge_Uprof = [0,1,2] + !--------------------------------------- + call nudging_set_profile(rlat,rlon,Nudge_Qprof,Wprof,pver) + Nudge_Qtau(icol,:,lchnk)=Wprof(:) + endif + + if(Nudge_PSprof.eq.3) then + ! Interp Horizontal window values from user input + !------------------------------------------------ + Hcoef = interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,Bnlon,Bnlat) + Nudge_PStau(icol,lchnk)=Hcoef + else + ! Set prodfile for Nudge_Uprof = [0,1,2] + !--------------------------------------- + Nudge_PStau(icol,lchnk)=nudging_set_PSprofile(rlat,rlon,Nudge_PSprof) + endif - Nudge_PStau(icol,lchnk)=nudging_set_PSprofile(rlat,rlon,Nudge_PSprof) end do + end do + + ! Initialize windowed Coeffcient profiles + ! Load zeros into nudging arrays + !------------------------------------------------------ + do lchnk=begchunk,endchunk + ncol=get_ncols_p(lchnk) Nudge_Utau(:ncol,:pver,lchnk) = & Nudge_Utau(:ncol,:pver,lchnk) * Nudge_Ucoef/float(Nudge_Step) Nudge_Vtau(:ncol,:pver,lchnk) = & @@ -1058,7 +1444,9 @@ subroutine nudging_init ! End Routine !------------ - + if(allocated(Bwindow)) deallocate(Bwindow) + if(allocated(B_lon )) deallocate(B_lon) + if(allocated(B_lat )) deallocate(B_lat) end subroutine nudging_init !================================================================ @@ -1071,12 +1459,14 @@ subroutine nudging_timestep_init(phys_state) ! arrays when necessary. Toggle the Nudging flag ! when the time is withing the nudging window. !=============================================================== - use physconst ,only: cpair - use physics_types,only: physics_state - use constituents ,only: cnst_get_ind - use dycore ,only: dycore_is - use ppgrid ,only: pver,pcols,begchunk,endchunk - use filenames ,only: interpret_filename_spec + use hycoef ,only: hyai, hybi, ps0, hyam, hybm + use error_messages,only: alloc_err + use physconst ,only: cpair + use physics_types ,only: physics_state + use constituents ,only: cnst_get_ind + use dycore ,only: dycore_is + use ppgrid ,only: pver,pverp,pcols,begchunk,endchunk + use filenames ,only: interpret_filename_spec use ESMF ! Arguments @@ -1102,12 +1492,32 @@ subroutine nudging_timestep_init(phys_state) real(r8) :: Sbar,Qbar,Wsum integer :: dtime + real(r8),allocatable:: hyai_w(:,:,:) + real(r8),allocatable:: hybi_w(:,:,:) + real(r8):: hyai_m(pcols,pverp) + real(r8):: hybi_m(pcols,pverp) + real(r8):: WETsum(pcols,pver) + real(r8):: PSpert(pcols) + real(r8):: Dlt_PS(pcols) + real(r8):: lnP0 (pcols,pverp) + real(r8):: lnP1 (pcols,pverp) + integer :: istat + ! Check if Nudging is initialized !--------------------------------- if(.not.Nudge_Initialized) then call endrun('nudging_timestep_init:: Nudging NOT Initialized') endif + ! PS nudging Needs some workspace + !------------------------------------- + if(Nudge_PSprof.ne.0) then + allocate(hyai_w(pcols,pverp,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','hyai_w',pcols*pverp*((endchunk-begchunk)+1)) + allocate(hybi_w(pcols,pverp,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','hybi_w',pcols*pverp*((endchunk-begchunk)+1)) + endif + ! Get time step size !-------------------- dtime = get_step_size() @@ -1176,7 +1586,8 @@ subroutine nudging_timestep_init(phys_state) Model_V(:ncol,:pver,lchnk)=phys_state(lchnk)%v(:ncol,:pver) Model_T(:ncol,:pver,lchnk)=phys_state(lchnk)%t(:ncol,:pver) Model_Q(:ncol,:pver,lchnk)=phys_state(lchnk)%q(:ncol,:pver,indw) - Model_PS(:ncol,lchnk)=phys_state(lchnk)%ps(:ncol) + Model_PS (:ncol,lchnk)=phys_state(lchnk)%ps(:ncol) + Model_PSfilt(:ncol,lchnk)=Model_PS(:ncol,lchnk) end do ! Load Dry Static Energy values for Model @@ -1199,27 +1610,93 @@ subroutine nudging_timestep_init(phys_state) end do endif - ! Optionally: Apply Zonal Filtering to Model state data - !------------------------------------------------------- - if(Nudge_ZonalFilter) then - call ZM%calc_amps(Model_U,Zonal_Bamp3d) - call ZM%eval_grid(Zonal_Bamp3d,Model_U) - - call ZM%calc_amps(Model_V,Zonal_Bamp3d) - call ZM%eval_grid(Zonal_Bamp3d,Model_V) - - call ZM%calc_amps(Model_T,Zonal_Bamp3d) - call ZM%eval_grid(Zonal_Bamp3d,Model_T) - - call ZM%calc_amps(Model_S,Zonal_Bamp3d) - call ZM%eval_grid(Zonal_Bamp3d,Model_S) + ! Process values differently if PS nudging is active + !---------------------------------------------------- + if(Nudge_PSprof.ne.0) then + ! Ps Nudging is active: + ! Compute hybrid coef for wet pressures on dry pressure surfaces. + !---------------------------------------------------------------- + do lchnk=begchunk,endchunk + ncol=phys_state(lchnk)%ncol - call ZM%calc_amps(Model_Q,Zonal_Bamp3d) - call ZM%eval_grid(Zonal_Bamp3d,Model_Q) + ! Forumate WET factor + !----------------------- + WETsum(:ncol,:pver) = 1._r8 + Model_Q(:ncol,:pver,lchnk) + + ! Integrate downward to compute moist A/B values + !---------------------------------------------- + hyai_m(:ncol,1) = hyai(1) + hybi_m(:ncol,1) = hybi(1) + do kk=1,pver + hyai_m(:ncol,kk+1) = hyai_m(:ncol,kk) + WETsum(:ncol,kk)*(hyai(kk+1)-hyai(kk)) + hybi_m(:ncol,kk+1) = hybi_m(:ncol,kk) + WETsum(:ncol,kk)*(hybi(kk+1)-hybi(kk)) + end do + + ! Now formulate the A/B coefs for wet surface pressures + !------------------------------------------------------ + do kk=1,pverp + hybi_w(:ncol,kk,lchnk) = hybi_m(:ncol,kk)/hybi_m(:ncol,pverp) + hyai_w(:ncol,kk,lchnk) = (hyai_m(:ncol, 1)+hyai_m(:ncol, kk)) & + -(hyai_m(:ncol, 1)+hyai_m(:ncol,pverp))*hybi_w(:ncol,kk,lchnk) + end do + end do - call ZM%calc_amps(Model_PS,Zonal_Bamp2d) - call ZM%eval_grid(Zonal_Bamp2d,Model_PS) + ! Optionally: Apply Zonal/Spectral Filtering to only to U,V,Q,Ps Model state data + !--------------------------------------------------------------------------------- + if(Nudge_ZonalFilter) then + call ZM%calc_amps(Model_U,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_U) + call ZM%calc_amps(Model_V,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_V) + call ZM%calc_amps(Model_Q,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_Q) + call ZM%calc_amps(Model_PSfilt,Zonal_Bamp2d) + call ZM%eval_grid(Zonal_Bamp2d,Model_PSfilt) + endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Model_U,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_U) + call SH%calc_amps(Model_V,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_V) + call SH%calc_amps(Model_Q,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_Q) + call SH%calc_amps(Model_PSfilt,Spectral_Bamp2d) + call SH%eval_grid(Spectral_Bamp2d,Model_PSfilt) + endif + else + ! Ps Nudging is not active: + ! Optionally: Apply Zonal/Spectral Filtering to all Model state data + !-------------------------------------------------------------------- + if(Nudge_ZonalFilter) then + call ZM%calc_amps(Model_U,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_U) + call ZM%calc_amps(Model_V,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_V) + call ZM%calc_amps(Model_T,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_T) + call ZM%calc_amps(Model_S,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_S) + call ZM%calc_amps(Model_Q,Zonal_Bamp3d) + call ZM%eval_grid(Zonal_Bamp3d,Model_Q) + call ZM%calc_amps(Model_PSfilt,Zonal_Bamp2d) + call ZM%eval_grid(Zonal_Bamp2d,Model_PSfilt) + endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Model_U,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_U) + call SH%calc_amps(Model_V,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_V) + call SH%calc_amps(Model_T,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_T) + call SH%calc_amps(Model_S,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_S) + call SH%calc_amps(Model_Q,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Model_Q) + call SH%calc_amps(Model_PSfilt,Spectral_Bamp2d) + call SH%eval_grid(Spectral_Bamp2d,Model_PSfilt) + endif endif + endif ! ((Before_End) .and. (Update_Model)) then !---------------------------------------------------------------- @@ -1286,9 +1763,11 @@ subroutine nudging_timestep_init(phys_state) endif if(.not.Nudge_ON) then if(masterproc) then - write(iulog,*) 'NUDGING: WARNING - analyses file NOT FOUND. Switching ' - write(iulog,*) 'NUDGING: nudging OFF to coast thru the gap. ' + write(iulog,*) 'NUDGING: WARNING - analyses file NOT FOUND. You can switch nudging ' + write(iulog,*) 'NUDGING: OFF to coast thru a known gap in your files ' + write(iulog,*) 'NUDGING: by commenting out the following endrun command.' endif + call endrun('nudging_timestep_init:: ERROR Missing Nudging File') endif else Nudge_ON=.false. @@ -1344,24 +1823,85 @@ subroutine nudging_timestep_init(phys_state) call endrun('nudging_timestep_init:: ERROR unknown Nudging_Force_Opt') endif - ! Now load Dry Static Energy values for Target - !--------------------------------------------- - if(Nudge_TSmode == 0) then - ! DSE tendencies from Temperature only - !--------------------------------------- - do lchnk=begchunk,endchunk - ncol=phys_state(lchnk)%ncol - Target_S(:ncol,:pver,lchnk)=cpair*Target_T(:ncol,:pver,lchnk) - end do - elseif(Nudge_TSmode == 1) then - ! Caluculate DSE tendencies from Temperature, Water Vapor, and Surface Pressure - !------------------------------------------------------------------------------ - do lchnk=begchunk,endchunk - ncol=phys_state(lchnk)%ncol - call calc_DryStaticEnergy(Target_T(:,:,lchnk), Target_Q(:,:,lchnk), & - phys_state(lchnk)%phis, Target_PS(:,lchnk), & - Target_S(:,:,lchnk), ncol) - end do + ! Set the Target_S() values used for Temperature nudging + !-------------------------------------------------------- + if(Nudge_PSprof.ne.0) then + ! Ps Nudging is active + !---------------------- + if(Nudge_PSscal < 1._r8) then + ! Use hybrid Bn coeffcients for the vertical + ! influence of Ps perturbations + !----------------------------------------------------- + do lchnk=begchunk,endchunk + ncol=phys_state(lchnk)%ncol + + ! Compute moist pressures for hybi vertical influence + !----------------------------------------------------- + Dlt_PS(:ncol)=(Target_PS(:ncol,lchnk) - Model_PSfilt(:ncol,lchnk)) + do kk=1,pverp + lnP0(:ncol,kk) = ps0*hyai_w(:ncol,kk,lchnk) + hybi_w(:ncol,kk,lchnk)*Model_PS(:ncol,lchnk) + lnP1(:ncol,kk) = lnP0(:ncol,kk) + hybi_w(:ncol,kk,lchnk)*Dlt_PS(:ncol) + lnP0(:ncol,kk) = log(lnP0(:ncol,kk)) + lnP1(:ncol,kk) = log(lnP1(:ncol,kk)) + end do + + ! Compute Target_T/Target_S values for Ps perturbation + !------------------------------------------------------ + do kk=1,pver + PSpert(:ncol)= (lnP0(:ncol,kk+1)-lnP0(:ncol,kk))/(lnP1(:ncol,kk+1)-lnP1(:ncol,kk)) + Target_T(:ncol,kk,lchnk)=Model_T(:ncol,kk,lchnk)*PSpert(:ncol) + Target_S(:ncol,kk,lchnk)=cpair*Target_T(:ncol,kk,lchnk) + end do ! kk=1,pver + end do + else! (Nudge_PSscal >= 1._r8) + ! Use the value of Nudge_PSscal to set the scale + ! height for the vertical influence function: + ! (Pref/P0)**Nudge_PSscal + !----------------------------------------------------- + do lchnk=begchunk,endchunk + ncol=phys_state(lchnk)%ncol + + ! Compute moist pressures for (Pref/P0) vertical influence + !--------------------------------------------------------- + Dlt_PS(:ncol)=(Target_PS(:ncol,lchnk) - Model_PSfilt(:ncol,lchnk)) + do kk=1,pverp + lnP0(:ncol,kk) = ps0*hyai_w(:ncol,kk,lchnk) + hybi_w(:ncol,kk,lchnk)*Model_PS(:ncol,lchnk) + lnP1(:ncol,kk) = lnP0(:ncol,kk) & + +Dlt_PS(:ncol)*((hyai_w(:ncol,kk,lchnk)+hybi_w(:ncol,kk,lchnk))**Nudge_PSscal) + lnP0(:ncol,kk) = log(lnP0(:ncol,kk)) + lnP1(:ncol,kk) = log(lnP1(:ncol,kk)) + end do + + ! Compute Target_T/Target_S values for Ps perturbation + !------------------------------------------------------ + do kk=1,pver + PSpert(:ncol)= (lnP0(:ncol,kk+1)-lnP0(:ncol,kk))/(lnP1(:ncol,kk+1)-lnP1(:ncol,kk)) + Target_T(:ncol,kk,lchnk)=Model_T(:ncol,kk,lchnk)*PSpert(:ncol) + Target_S(:ncol,kk,lchnk)=cpair*Target_T(:ncol,kk,lchnk) + end do ! kk=1,pver + end do + endif + else + ! Ps Nudging is NOT active + ! Now load Dry Static Energy values for Target + !--------------------------------------------- + if(Nudge_TSmode == 0) then + ! DSE tendencies from Temperature only + !--------------------------------------- + do lchnk=begchunk,endchunk + ncol=phys_state(lchnk)%ncol + Target_S(:ncol,:pver,lchnk)=cpair*Target_T(:ncol,:pver,lchnk) + end do + elseif(Nudge_TSmode == 1) then + ! Caluculate DSE tendencies from Temperature, Water Vapor, and Surface Pressure + !------------------------------------------------------------------------------ + do lchnk=begchunk,endchunk + ncol=phys_state(lchnk)%ncol + call calc_DryStaticEnergy(Target_T(:,:,lchnk), Target_Q(:,:,lchnk), & + phys_state(lchnk)%phis, Target_PS(:,lchnk), & + Target_S(:,:,lchnk), ncol) + end do + endif endif ! Set Tscale for the specified Forcing Option @@ -1410,7 +1950,8 @@ subroutine nudging_timestep_init(phys_state) ! write(iulog,*) 'PFC: Target_S(1,:pver,begchunk)=',Target_S(1,:pver,begchunk) ! write(iulog,*) 'PFC: Model_S(1,:pver,begchunk)=',Model_S(1,:pver,begchunk) ! write(iulog,*) 'PFC: Target_PS(1,begchunk)=',Target_PS(1,begchunk) -! write(iulog,*) 'PFC: Model_PS(1,begchunk)=',Model_PS(1,begchunk) +! write(iulog,*) 'PFC: Model_PS(1,begchunk)=',Model_PS (1,begchunk) +! write(iulog,*) 'PFC: Model_PSfilt(1,begchunk)=',Model_PSfilt(1,begchunk) ! write(iulog,*) 'PFC: Nudge_Sstep(1,:pver,begchunk)=',Nudge_Sstep(1,:pver,begchunk) ! write(iulog,*) 'PFC: Nudge_Xstep arrays updated:' ! endif @@ -1418,7 +1959,8 @@ subroutine nudging_timestep_init(phys_state) ! End Routine !------------ - + if(allocated(hyai_w)) deallocate(hyai_w) + if(allocated(hybi_w)) deallocate(hybi_w) end subroutine nudging_timestep_init !================================================================ @@ -1468,6 +2010,9 @@ subroutine nudging_timestep_tend(phys_state,phys_tend) call outfld('Target_V',Target_V(:,:,lchnk),pcols,lchnk) call outfld('Target_T',Target_T(:,:,lchnk),pcols,lchnk) call outfld('Target_Q',Target_Q(:,:,lchnk),pcols,lchnk) + call outfld( 'Nudge_PS',Nudge_PSstep(:,lchnk),pcols,lchnk) + call outfld('Target_PS', Target_PS(:,lchnk),pcols,lchnk) + call outfld( 'Model_PS',Model_PSfilt(:,lchnk),pcols,lchnk) endif ! End Routine @@ -1556,6 +2101,10 @@ subroutine nudging_update_analyses(anal_file) call ZM%calc_amps(Tmp3D,Zonal_Bamp3d) call ZM%eval_grid(Zonal_Bamp3d,Tmp3D) endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Tmp3D,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Tmp3D) + endif Nobs_U(:,:,begchunk:endchunk,Nudge_ObsInd(1)) = Tmp3D(:,:,begchunk:endchunk) else call endrun('Variable "U" is missing in '//trim(anal_file)) @@ -1569,6 +2118,10 @@ subroutine nudging_update_analyses(anal_file) call ZM%calc_amps(Tmp3D,Zonal_Bamp3d) call ZM%eval_grid(Zonal_Bamp3d,Tmp3D) endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Tmp3D,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Tmp3D) + endif Nobs_V(:,:,begchunk:endchunk,Nudge_ObsInd(1)) = Tmp3D(:,:,begchunk:endchunk) else call endrun('Variable "V" is missing in '//trim(anal_file)) @@ -1582,6 +2135,10 @@ subroutine nudging_update_analyses(anal_file) call ZM%calc_amps(Tmp3D,Zonal_Bamp3d) call ZM%eval_grid(Zonal_Bamp3d,Tmp3D) endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Tmp3D,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Tmp3D) + endif Nobs_T(:,:,begchunk:endchunk,Nudge_ObsInd(1)) = Tmp3D(:,:,begchunk:endchunk) else call endrun('Variable "T" is missing in '//trim(anal_file)) @@ -1595,6 +2152,10 @@ subroutine nudging_update_analyses(anal_file) call ZM%calc_amps(Tmp3D,Zonal_Bamp3d) call ZM%eval_grid(Zonal_Bamp3d,Tmp3D) endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Tmp3D,Spectral_Bamp3d) + call SH%eval_grid(Spectral_Bamp3d,Tmp3D) + endif Nobs_Q(:,:,begchunk:endchunk,Nudge_ObsInd(1)) = Tmp3D(:,:,begchunk:endchunk) else call endrun('Variable "Q" is missing in '//trim(anal_file)) @@ -1608,6 +2169,10 @@ subroutine nudging_update_analyses(anal_file) call ZM%calc_amps(Tmp2D,Zonal_Bamp2d) call ZM%eval_grid(Zonal_Bamp2d,Tmp2D) endif + if(Nudge_SpectralFilter) then + call SH%calc_amps(Tmp2D,Spectral_Bamp2d) + call SH%eval_grid(Spectral_Bamp2d,Tmp2D) + endif Nobs_PS(:,begchunk:endchunk,Nudge_ObsInd(1)) = Tmp2D(:,begchunk:endchunk) else call endrun('Variable "PS" is missing in '//trim(anal_file)) @@ -1675,7 +2240,7 @@ subroutine nudging_set_profile(rlat,rlon,Nudge_prof,Wprof,nlev) !------------------------------------------ latx=rlat-Nudge_Hwin_lat0 lonx=rlon-Nudge_Hwin_lon0 - if(lonx > 180._r8) lonx=lonx-360._r8 + if(lonx > 180._r8) lonx=lonx-360._r8 if(lonx <= -180._r8) lonx=lonx+360._r8 ! Calcualte RAW window value @@ -1733,6 +2298,131 @@ subroutine nudging_set_profile(rlat,rlon,Nudge_prof,Wprof,nlev) end subroutine nudging_set_profile !================================================================ + + !================================================================ + real(r8) function nudging_set_PSprofile(rlat,rlon,Nudge_PSprof) + ! + ! NUDGING_SET_PSPROFILE: for the given lat and lon set the surface + ! pressure profile value for the specified index. + ! Values range from 0. to 1. to affect spatial + ! variations on nudging strength. + !=============================================================== + + ! Arguments + !-------------- + real(r8) :: rlat,rlon + integer :: Nudge_PSprof + + ! Local values + !---------------- + real(r8) :: Hcoef,latx,lonx + real(r8) :: lon_lo,lon_hi,lat_lo,lat_hi,lev_lo,lev_hi + + !--------------- + ! set coeffcient + !--------------- + if(Nudge_PSprof == 0) then + ! No Nudging + !------------- + nudging_set_PSprofile=0.0_r8 + elseif(Nudge_PSprof == 1) then + ! Uniform Nudging + !----------------- + nudging_set_PSprofile=1.0_r8 + elseif(Nudge_PSprof == 2) then + ! Localized Nudging with specified Heaviside window function + !------------------------------------------------------------ + if(Nudge_Hwin_max <= Nudge_Hwin_min) then + ! For a constant Horizontal window function, + ! just set Hcoef to the maximum of Hlo/Hhi. + !-------------------------------------------- + Hcoef=max(Nudge_Hwin_lo,Nudge_Hwin_hi) + else + ! get lat/lon relative to window center + !------------------------------------------ + latx=rlat-Nudge_Hwin_lat0 + lonx=rlon-Nudge_Hwin_lon0 + if(lonx > 180._r8) lonx=lonx-360._r8 + if(lonx <= -180._r8) lonx=lonx+360._r8 + + ! Calcualte RAW window value + !------------------------------- + lon_lo=(Nudge_Hwin_lonWidthH+lonx)/Nudge_Hwin_lonDelta + lon_hi=(Nudge_Hwin_lonWidthH-lonx)/Nudge_Hwin_lonDelta + lat_lo=(Nudge_Hwin_latWidthH+latx)/Nudge_Hwin_latDelta + lat_hi=(Nudge_Hwin_latWidthH-latx)/Nudge_Hwin_latDelta + Hcoef=((1._r8+tanh(lon_lo))/2._r8)*((1._r8+tanh(lon_hi))/2._r8) & + *((1._r8+tanh(lat_lo))/2._r8)*((1._r8+tanh(lat_hi))/2._r8) + + ! Scale the horizontal window coef for specfied range of values. + !-------------------------------------------------------- + Hcoef=(Hcoef-Nudge_Hwin_min)/(Nudge_Hwin_max-Nudge_Hwin_min) + Hcoef=(1._r8-Hcoef)*Nudge_Hwin_lo + Hcoef*Nudge_Hwin_hi + endif + nudging_set_PSprofile=Hcoef + else + call endrun('nudging_set_PSprofile:: Unknown Nudge_prof value') + endif + + ! End Routine + !------------ + + end function nudging_set_PSprofile + !================================================================ + + + !================================================================ + subroutine nudging_set_Vwindow(Wprof,nlev) + ! + ! NUDGING_SET_VWINDOW: for the currently set namelist values + ! return the verical profile of window coeffcients. + ! Values range from 0. to 1. to affect spatial + ! variations on nudging strength. + !=============================================================== + + ! Arguments + !-------------- + integer :: nlev + real(r8) :: Wprof(nlev) + + ! Local values + !---------------- + real(r8) :: latx,lonx,Vmax,Vmin + real(r8) :: lev_lo,lev_hi + integer :: ilev + + ! Load the RAW vertical window + !------------------------------ + do ilev=1,nlev + lev_lo=(float(ilev)-Nudge_Vwin_Lindex)/Nudge_Vwin_Ldelta + lev_hi=(Nudge_Vwin_Hindex-float(ilev))/Nudge_Vwin_Hdelta + Wprof(ilev)=((1._r8+tanh(lev_lo))/2._r8)*((1._r8+tanh(lev_hi))/2._r8) + end do + + ! Scale the Window function to span the values between Vlo and Vhi: + !----------------------------------------------------------------- + Vmax=maxval(Wprof) + Vmin=minval(Wprof) + if((Vmax <= Vmin) .or. ((Nudge_Vwin_Hindex >= (nlev+1)) .and. & + (Nudge_Vwin_Lindex <= 0 ) )) then + ! For a constant Vertical window function, + ! load maximum of Vlo/Vhi into Wprof() + !-------------------------------------------- + Vmax=max(Nudge_Vwin_lo,Nudge_Vwin_hi) + Wprof(:)=Vmax + else + ! Scale the RAW vertical window for specfied range of values. + !-------------------------------------------------------- + Wprof(:)=(Wprof(:)-Vmin)/(Vmax-Vmin) + Wprof(:)=Nudge_Vwin_lo + Wprof(:)*(Nudge_Vwin_hi-Nudge_Vwin_lo) + endif + + ! End Routine + !------------ + end subroutine nudging_set_Vwindow + !================================================================ + + !================================================================ subroutine nudging_final @@ -1748,6 +2438,7 @@ subroutine nudging_final if (allocated(Model_S)) deallocate(Model_S) if (allocated(Model_Q)) deallocate(Model_Q) if (allocated(Model_PS)) deallocate(Model_PS) + if (allocated(Model_PSfilt)) deallocate(Model_PSfilt) if (allocated(Nudge_Utau)) deallocate(Nudge_Utau) if (allocated(Nudge_Vtau)) deallocate(Nudge_Vtau) if (allocated(Nudge_Stau)) deallocate(Nudge_Stau) @@ -1768,50 +2459,15 @@ subroutine nudging_final if (allocated(Nobs_PS)) deallocate(Nobs_PS) if (allocated(Zonal_Bamp2d)) deallocate(Zonal_Bamp2d) if (allocated(Zonal_Bamp3d)) deallocate(Zonal_Bamp3d) + if (allocated(Spectral_Bamp2d)) deallocate(Spectral_Bamp2d) + if (allocated(Spectral_Bamp3d)) deallocate(Spectral_Bamp3d) call ZM%final() + call SH%final() end subroutine nudging_final !================================================================ - !================================================================ - real(r8) function nudging_set_PSprofile(rlat,rlon,Nudge_PSprof) - ! - ! NUDGING_SET_PSPROFILE: for the given lat and lon set the surface - ! pressure profile value for the specified index. - ! Values range from 0. to 1. to affect spatial - ! variations on nudging strength. - !=============================================================== - - ! Arguments - !-------------- - real(r8) :: rlat,rlon - integer :: Nudge_PSprof - - ! Local values - !---------------- - - !--------------- - ! set coeffcient - !--------------- - if(Nudge_PSprof == 0) then - ! No Nudging - !------------- - nudging_set_PSprofile=0.0_r8 - elseif(Nudge_PSprof == 1) then - ! Uniform Nudging - !----------------- - nudging_set_PSprofile=1.0_r8 - else - call endrun('nudging_set_PSprofile:: Unknown Nudge_prof value') - endif - - ! End Routine - !------------ - - end function nudging_set_PSprofile - !================================================================ - !================================================================ subroutine calc_DryStaticEnergy(t, q, phis, ps, dse, ncol) @@ -1906,4 +2562,110 @@ subroutine calc_DryStaticEnergy(t, q, phis, ps, dse, ncol) end subroutine calc_DryStaticEnergy !================================================================ + + !================================================================ + real(r8) function interp_Bwin(rlon,rlat,B_lon,B_lat,Bwindow,nlon,nlat) + ! + ! interp_Bwin: for the given lat and lon return the interpolated + ! Bounary Nudinging coef from the given lat,lon + ! array of coeffcients. + !=============================================================== + + ! Arguments + !-------------- + integer nlat,nlon + real(r8) rlat,rlon + real(r8) B_lon (nlon) + real(r8) B_lat (nlat) + real(r8) Bwindow(nlon,nlat) + + ! Local values + !---------------- + integer ii + integer LonLo,LonHi,LatLo,LatHi + integer LonLx,LonHx,LatLx,LatHx + real(r8) LonFrac,LatFrac + + ! Get the Bounding indices for longitude + !---------------------------------------- + if(rlon.lt.B_lon(1)) then + LonLo = nlon + LonHi = 1 + LonFrac = ( rlon +360._r8 - B_lon(nlon)) & + /(B_lon(1)+360._r8 - B_lon(nlon)) + elseif(rlon.gt.B_lon(nlon)) then + LonLo = nlon + LonHi = 1 + LonFrac = ( rlon - B_lon(nlon)) & + /(B_lon(1)+360._r8 - B_lon(nlon)) + else + do ii=1,(nlon-1) + if((rlon.ge.B_lon(ii)).and.(rlon.lt.B_lon(ii+1))) then + LonLo = ii + LonHi = ii+1 + LonFrac = ( rlon - B_lon(ii)) & + /(B_lon(ii+1) - B_lon(ii)) + exit + endif + end do + endif + + ! Get the Bounding indices for latitude + !--------------------------------------- + if(rlat.lt.B_lat(1)) then + LatLo = 0 + LatHi = 1 + LatFrac = (rlat + B_lat(1))/(2._r8*B_lat(1)) + elseif(rlat.gt.B_lat(nlat)) then + LatLo = nlat + LatHi = nlat + 1 + LatFrac = ( rlat - B_lat(nlat)) & + /(180._r8 - 2._r8*B_lat(nlat)) + else + do ii=1,(nlat-1) + if((rlat.ge.B_lat(ii)).and.(rlat.lt.B_lat(ii+1))) then + LatLo = ii + LatHi = ii+1 + LatFrac = ( rlat - B_lat(ii)) & + /(B_lat(ii+1) - B_lat(ii)) + exit + endif + end do + endif + + ! Interpolate + !------------- + if(LatHi.gt.nlat) then + LatHx = LatLo + LonLx = LonHi + (nlon/2) + LonHx = LonLo + (nlon/2) + if(LonLx.gt.nlon) LonLx = LonLx - nlon + if(LonHx.gt.nlon) LonHx = LonHx - nlon + interp_Bwin = Bwindow(LonLo,LatLo)*(1._r8-LonFrac)*(1._r8-LatFrac) & + +Bwindow(LonHi,LatLo)* LonFrac *(1._r8-LatFrac) & + +Bwindow(LonHx,LatHx)* LonFrac * LatFrac & + +Bwindow(LonLx,LatHx)*(1._r8-LonFrac)* LatFrac + elseif(LatLo.lt. 1 ) then + LatLx = LatHi + LonLx = LonHi + (nlon/2) + LonHx = LonLo + (nlon/2) + if(LonLx.gt.nlon) LonLx = LonLx - nlon + if(LonHx.gt.nlon) LonHx = LonHx - nlon + interp_Bwin = Bwindow(LonLx,LatLx)*(1._r8-LonFrac)*(1._r8-LatFrac) & + +Bwindow(LonHx,LatLx)* LonFrac *(1._r8-LatFrac) & + +Bwindow(LonHi,LatHi)* LonFrac * LatFrac & + +Bwindow(LonLo,LatHi)*(1._r8-LonFrac)* LatFrac + else + interp_Bwin = Bwindow(LonLo,LatLo)*(1._r8-LonFrac)*(1._r8-LatFrac) & + +Bwindow(LonHi,LatLo)* LonFrac *(1._r8-LatFrac) & + +Bwindow(LonHi,LatHi)* LonFrac * LatFrac & + +Bwindow(LonLo,LatHi)*(1._r8-LonFrac)* LatFrac + endif + + ! End Routine + !------------ + return + end function ! interp_Bwin + !================================================================ + end module nudging diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 4d38ed5b93..22aab85a2f 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -7,7 +7,7 @@ module phys_grid_ctem use ppgrid, only: begchunk, endchunk, pcols, pver use physics_types, only: physics_state use cam_history, only: addfld, outfld - use zonal_mean_mod,only: ZonalAverage_t, ZonalMean_t + use ug_spectralmethods_mod,only: ZonalAverage_t, ZonalMean_t use physconst, only: pi use cam_logfile, only: iulog use cam_abortutils,only: endrun, handle_allocate_error diff --git a/src/utils/ug_spectralmethods_mod.F90 b/src/utils/ug_spectralmethods_mod.F90 new file mode 100644 index 0000000000..0dc06b76db --- /dev/null +++ b/src/utils/ug_spectralmethods_mod.F90 @@ -0,0 +1,3869 @@ +module ug_spectralmethods_mod +!====================================================================== +! +! Purpose: Implement Spherical Harmonic Analysis/Synthesis methods +! for Unstructured Grids(UG). +! +! Spherical Harmonic Analysis/Synthesis: +! ====================================== +! This module implements 2 classes for the spectral analysis and +! synthesis of spherical harmonic functions. For each, the global +! integrals to compute basis amlitudes are implemented with an area +! weighted sum of the UG gridpoints on the surface of the sphere. +! +! One computes ampitudes via a least-squares fit in the typical mannor. +! This approach distributes numerical errors in the basis representation +! based on covariances, but it can be very costly for the MPI environment. +! The other approach orthogonalizes the spherical harmonic basis functions +! via the Gram-Schmidt method. While this is much less computationally +! expensive, it propagates numerical errors in the basis representations +! from the gravest scales down to the higher resolution modes. +! +! SphericalHarmonic_COV_t: For the analysis/synthesis of spherical harmonic +! functions and amplitudes on a 2D/3D grids of points +! distributed over the surface of a sphere. +! Amplitudes are computed based on COV values between +! basis functions. +! +! SphericalHarmonic_GS_t: For the analysis/synthesis of spherical harmonic +! functions and amplitudes on a 2D/3D grids of points +! distributed over the surface of a sphere. +! Orthonormal basis functions are constructed via +! Gram-Schmidt. +! +! Zonal Mean Analysis/Synthesis: +! ====================================== +! This module implements 3 classes for the spectral analysis and +! synthesis of zonal mean values based on m=0 spherical harmonics. +! +! ZonalMean_t: For the analysis/synthesis of zonal mean values +! on a 2D/3D grid of points distributed over the +! surface of a sphere. +! ZonalProfile_t: For the analysis/synthesis of zonal mean values +! on a meridional grid that spans the latitudes +! from SP to NP +! ZonalAverage_t: To calculate zonal mean values via a simple +! area weighted bin-averaging of 2D/3D grid points +! assigned to each latitude band. +! +! NOTE: The weighting of the Zonal Profiles values is scaled such +! that ZonalMean_t amplitudes can be used to evaluate values +! on the ZonalProfile_t grid and vice-versa. +! +! The ZonalMean_t computes global integrals to compute basis +! amplitudes. For distributed environments the cost of these +! can be reduced using the The ZonalAverage_t data structures. +! +! ====================== +! BASIS SAMPLING ERRORS: +! ====================== +! The variable values on the unstructured gridpoints represent the +! average value over the associated gridpoint's area, but the basis +! values used at these gridpoints are just the spherical harmonics +! evaluated at that point. These sampling errors can be significant +! particularly in the polar regions where the range of lat/lon values +! can vary greatly over the domain of a finite area. +! +! To suppress these errors, during the initailization process, each +! spherical harmonic basis functions is averaged over a local circular +! domain centered at each gridpoint. The radius of this domain is set +! to match the area associated with the given gridpoint. The number of +! equal-area gridpoints within this circular domain is specified by the +! user via an optional argument to the initialization routine. While the +! representation of the area is approximate, the goal here is for the +! gridpoint values of each basis function to also represent the local +! area average over the associated gridpoint's area. +! +! The gridpoints in the circular domain are constructed in the the form +! of concentric rings such that each point has equal area weighting. The +! users controls the density of samping points by setting the number of +! rings. +! +! There is always 1 point at center, each additional +! ring (kk) contains 8*(kk-1) sample points. The total +! number of sample points (Nsamp) is thus: +! +! SAMPLE_NRING=1 Nsamp=1 (DEFAULT SampleGrid NOT used) +! SAMPLE_NRING=2 Nsamp=9 +! SAMPLE_NRING=3 Nsamp=25 +! SAMPLE_NRING=4 Nsamp=49 +! +! The SampleGrid is only used during initialization, after that there +! is no additional computational cost. +! +! +! USAGE: +! +! Spherical Harmonic Analysis/Synthesis: +! ====================================== +! Compute Spherical Harmonic amplitudes and synthesize values on 2D/3D physgrid. +! All of the differences between these classes are internal, so methods/interfaces +! are the same for the user: +! +! Usage: type(SphericalHarmonic_COV_t):: SH +! Usage: type(SphericalHarmonic_GS_t ):: SH +! ========================================= +! call SH%init(nmax,nbas,SAMPLE_NRING) +! ------------------ +! - Initialize the data structure for the given spherical +! truncation 'nmax' and return 'nbas',the number of spherical +! harmonic basis functions. +! +! Arguments: +! integer,intent(in ) :: nmax -Number of meridional modes +! integer,intent(out) :: nbas -Total number spherical harmonic functions +! integer,intent(in ),optional:: SAMPLE_NRING (Default=1) +! +! call SH%calc_amps(Gdata,Bamp) +! ----------------------------- +! - For the initialized SphericalHarmonic_t; Given Gdata() values on +! the physgrid, compute the harmonic basis amplitudes Bamp(). +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk) +! real(r8),intent(out):: Bamp (nbas) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk) +! real(r8),intent(out):: Bamp (nbas,pver) +! +! call SH%eval_grid(Bamp,Gdata) +! ----------------------------- +! - For the initialized SphericalHarmonic_t; Given Bamp() spherical +! harmonic basis amplitudes, compute the Gdata() values on the physgrid. +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Bamp (nbas) +! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Bamp (nbas,pver) +! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk) +! +! USAGE: +! +! Zonal Mean Analysis/Synthesis: +! ====================================== +! (1) Compute Zonal mean amplitudes and synthesize values on 2D/3D physgrid +! +! Usage: type(ZonalMean_t):: ZM +! ========================================= +! call ZM%init(nbas,SAMPLE_NRING) +! --------------------------------- +! - Initialize the data structure with 'nbas' basis functions +! for the given physgrid latitudes and areas. +! +! Arguments: +! integer,intent(in) :: nbas -Number of m=0 spherical harmonics +! integer,intent(in),optional:: SAMPLE_NRING (Default=1) +! +! call ZM%calc_amps(Gdata,Bamp) +! ----------------------------- +! - For the initialized ZonalMean_t; Given Gdata() values on the physgrid, +! compute the zonal mean basis amplitudes Bamp(). +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk) +! real(r8),intent(out):: Bamp (nbas) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk) +! real(r8),intent(out):: Bamp (nbas,pver) +! +! call ZM%eval_grid(Bamp,Gdata) +! ----------------------------- +! - For the initialized ZonalMean_t; Given Bamp() zonal mean basis +! amplitudes, compute the Gdata() values on the physgrid. +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Bamp (nbas) +! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Bamp (nbas,pver) +! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk) +! +! +! (2) Compute Zonal mean amplitudes and synthesize values on Zonal profile grid +! +! Usage: type(ZonalProfile_t):: ZP +! ========================================= +! call ZP%init(lats,area,nlat,nbas,GEN_GAUSSLATS=.true.) +! ------------------------------------------------------ +! - Initialize the data structure for the given number of +! latitudes. Either use the given Latitudes and weights, +! or OPTIONALLY create profile gridpoints and associated +! area weights from SP to NP. Then initialize 'nbas' basis +! functions for the profile gridpoints. +! If the user supplies the lats/area values, the area values must +! be correctly scaled such that the global area adds up to 4PI. +! Otherwise, the ampitudes between ZonalProfile_t and ZonalMean_t +! are not interchangable. +! +! Arguments: +! real(r8),intent(inout):: lats(:) - Latitudes of meridional grid. +! real(r8),intent(inout):: area(:) - Area of each meridional gridpoint. +! integer ,intent(in) :: nlat - Number of meridional gridpoints. +! integer ,intent(in) :: nbas - Number of m=0 spherical harmonics +! logical ,intent(in),optional:: GEN_GAUSLATS - Flag to generate +! lats/areas values. +! +! call ZP%calc_amps(Zdata,Bamp) +! ----------------------------- +! - Given Zdata() on the Zonal profile grid, compute the +! zonal basis amplitudes Bamp(). +! +! Interface: 1D data on (nlat) grid +! real(r8),intent(in ):: Zdata(nlat) - Meridional Profile data +! real(r8),intent(out):: Bamp (nbas) - Zonal Basis Amplitudes +! +! Interface: 2D data on (nlat,pver) grid +! real(r8),intent(in ):: Zdata(nlat,pver) - Meridional Profile data +! real(r8),intent(out):: Bamp (nbas,pver) - Zonal Basis Amplitudes +! +! call ZP%eval_grid(Bamp,Zdata) +! ----------------------------- +! - Given Bamp() zonal basis amplitudes, evaluate the Zdata() +! values on the Zonal profile grid. +! +! Interface: 1D data on (nlat) grid +! real(r8),intent(in ):: Bamp (nbas) - Zonal Basis Amplitudes +! real(r8),intent(out):: Zdata(nlat) - Meridional Profile data +! +! Interface: 2D data on (nlat,pver) grid +! real(r8),intent(in ):: Bamp (nbas,pver) - Zonal Basis Amplitudes +! real(r8),intent(out):: Zdata(nlat,pver) - Meridional Profile data +! +! (3) Compute Zonal mean averages (FASTER/LESS-ACCURATE) on Zonal profile grid +! (For the created zonal profile, just bin average area weighted +! 2D/3D physgrid grid values) +! +! Usage: type(ZonalAverage_t):: ZA +! ========================================= +! call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.,USE_LINEARWGTS=.true.) +! -------------------------------------------------- +! - Given the latitude/area for the nlat meridional gridpoints, initialize +! the ZonalAverage data structure for computing bin-averaging of physgrid +! values. It is assumed that the domain of these gridpoints of the +! profile span latitudes from SP to NP. +! +! The optional GEN_GAUSSLATS flag allows for the generation of Gaussian +! latitude gridpoints. The generated grid over-writes the given values +! lats and area passed by the user. +! +! The optional USE_LINEARWGTS flag allows for the linear weighting of +! gridpoint values as they are distributed to latitude bins, or to assign +! grid point values to a single bin with equal weighting. +! (default=TRUE for linear weighting) +! +! Arguments: +! real(r8),intent(inout) :: lats(nlat) - Latitudes of meridional grid. +! real(r8),intent(inout) :: area(nlat) - Area of meridional gridpoints. +! integer,intent(in) :: nlat - Number of meridional gridpoints +! logical,intent(in),optional:: GEN_GAUSLATS - Flag to generate lats/areas values. +! logical,intent(in),optional:: USE_LINEARWGTS - Flag for weighting option. +! +! call ZA%binAvg(Gdata,Zdata) +! --------------------------- +! - For the initialized ZonalAverage_t; Given Gdata() on the physgrid, +! compute bin averages and return Zdata() on the Zonal profile grid. +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk) +! real(r8),intent(out):: Zdata(nlat) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk) +! real(r8),intent(out):: Zdata(nlat,pver) +! +! call ZA%setVal(Zdata,Gdata) +! --------------------------- +! - For the initialized ZonalAverage_t; Given Zdata for the profile bins, +! set the 2D physgrid values in each bin to the corresponding value. +! +! Interface: 2D data on the physgrid +! real(r8),intent(in ):: Zdata(nlat) +! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk) +! +! Interface: 3D data on the physgrid +! real(r8),intent(in ):: Zdata(nlat,pver) +! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk) +! +!====================================================================== + use shr_kind_mod, only: r8=>SHR_KIND_R8 + use phys_grid, only: get_ncols_p, get_rlat_p, get_rlon_p, get_wght_all_p, get_nlcols_p + use ppgrid, only: begchunk, endchunk, pcols + use shr_reprosum_mod,only: shr_reprosum_calc + use cam_abortutils, only: endrun, handle_allocate_error + use spmd_utils, only: mpicom + use physconst, only: pi + use phys_grid, only: ngcols_p => num_global_phys_cols + use cam_logfile, only: iulog + + implicit none + private + + public :: SphericalHarmonic_COV_t + public :: SphericalHarmonic_GS_t + public :: ZonalMean_t + public :: ZonalProfile_t + public :: ZonalAverage_t + + ! Type definitions + !------------------- + type SphericalHarmonic_COV_t + private + integer :: nmax + integer :: nbas + integer :: nring + real(r8),allocatable:: area (:,:) + real(r8),allocatable:: basis(:,:,:) + real(r8),allocatable:: map (:,:) + integer ,allocatable:: bindx(:,:) + contains + procedure,pass:: init => init_SphericalHarmonic_COV + procedure,pass:: final => final_SphericalHarmonic_COV + generic,public:: calc_amps => calc_SphericalHarmonic_2Damps_COV, & + calc_SphericalHarmonic_3Damps_COV + generic,public:: eval_grid => eval_SphericalHarmonic_2Dgrid_COV, & + eval_SphericalHarmonic_3Dgrid_COV + procedure,private,pass:: calc_SphericalHarmonic_2Damps_COV + procedure,private,pass:: calc_SphericalHarmonic_3Damps_COV + procedure,private,pass:: eval_SphericalHarmonic_2Dgrid_COV + procedure,private,pass:: eval_SphericalHarmonic_3Dgrid_COV + end type SphericalHarmonic_COV_t + + type SphericalHarmonic_GS_t + private + integer :: nmax + integer :: nbas + integer :: nring + real(r8),allocatable:: area (:,:) + real(r8),allocatable:: basis(:,:,:) + contains + procedure,pass:: init => init_SphericalHarmonic_GS + procedure,pass:: final => final_SphericalHarmonic_GS + generic,public:: calc_amps => calc_SphericalHarmonic_2Damps_GS, & + calc_SphericalHarmonic_3Damps_GS + generic,public:: eval_grid => eval_SphericalHarmonic_2Dgrid_GS, & + eval_SphericalHarmonic_3Dgrid_GS + procedure,private,pass:: calc_SphericalHarmonic_2Damps_GS + procedure,private,pass:: calc_SphericalHarmonic_3Damps_GS + procedure,private,pass:: eval_SphericalHarmonic_2Dgrid_GS + procedure,private,pass:: eval_SphericalHarmonic_3Dgrid_GS + end type SphericalHarmonic_GS_t + + type ZonalMean_t + private + integer :: nbas + integer :: nring + real(r8),allocatable:: area (:,:) + real(r8),allocatable:: basis(:,:,:) + real(r8),allocatable:: map (:,:) + contains + procedure,pass:: init => init_ZonalMean + procedure,pass:: final => final_ZonalMean + generic,public:: calc_amps => calc_ZonalMean_2Damps, & + calc_ZonalMean_3Damps + generic,public:: eval_grid => eval_ZonalMean_2Dgrid, & + eval_ZonalMean_3Dgrid + procedure,private,pass:: calc_ZonalMean_2Damps + procedure,private,pass:: calc_ZonalMean_3Damps + procedure,private,pass:: eval_ZonalMean_2Dgrid + procedure,private,pass:: eval_ZonalMean_3Dgrid + end type ZonalMean_t + + type ZonalProfile_t + private + integer :: nlat + integer :: nbas + real(r8),allocatable:: area (:) + real(r8),allocatable:: basis(:,:) + real(r8),allocatable:: map (:,:) + contains + procedure,pass:: init => init_ZonalProfile + procedure,pass:: final => final_ZonalProfile + generic,public:: calc_amps => calc_ZonalProfile_1Damps, & + calc_ZonalProfile_2Damps + generic,public:: eval_grid => eval_ZonalProfile_1Dgrid, & + eval_ZonalProfile_2Dgrid + procedure,private,pass:: calc_ZonalProfile_1Damps + procedure,private,pass:: calc_ZonalProfile_2Damps + procedure,private,pass:: eval_ZonalProfile_1Dgrid + procedure,private,pass:: eval_ZonalProfile_2Dgrid + end type ZonalProfile_t + + type ZonalAverage_t + private + logical :: LINEARWGTS = .false. ! DEFAULT: Should be .true. ? + integer :: nlat + real(r8),allocatable:: area (:) + real(r8),allocatable:: a_norm (:) + real(r8),allocatable:: area_g (:,:) + integer ,allocatable:: idx_map(:,:) + integer ,allocatable:: idx_wgt(:,:) + contains + procedure,pass:: init => init_ZonalAverage + procedure,pass:: final => final_ZonalAverage + generic,public:: binAvg => calc_ZonalAverage_2DbinAvg, & + calc_ZonalAverage_3DbinAvg + generic,public:: setVal => set_ZonalAverage_2Dgrid, & + set_ZonalAverage_3Dgrid + procedure,private,pass:: calc_ZonalAverage_2DbinAvg + procedure,private,pass:: calc_ZonalAverage_3DbinAvg + procedure,private,pass:: set_ZonalAverage_2Dgrid + procedure,private,pass:: set_ZonalAverage_3Dgrid + end type ZonalAverage_t + + real(r8), parameter :: halfPI = 0.5_r8*pi + real(r8), parameter :: twoPI = 2.0_r8*pi + real(r8), parameter :: fourPI = 4.0_r8*pi + real(r8), parameter :: qrtrPI = 0.25_r8*pi + real(r8), parameter :: invSqrt4pi = 1._r8/sqrt(fourPI) + +contains + !======================================================================= + subroutine init_SphericalHarmonic_COV(this,I_nmax,O_nbas,SAMPLE_NRING) + ! + ! init_SphericalHarmonic: Initialize the SphericalHarmonic data structure + ! for the physics grid. It is assumed that the domain + ! of these gridpoints spans the surface of the sphere. + ! The representation of basis functions is + ! normalized w.r.t integration over the sphere. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_COV_t) :: this + integer ,intent(in ):: I_nmax + integer ,intent(out):: O_nbas + integer ,intent(in ),optional:: SAMPLE_NRING + ! + ! Local Values + !-------------- + real(r8),allocatable:: Clons(:,:) + real(r8),allocatable:: Clats(:,:) + real(r8),allocatable:: Bcoef(:) + real(r8),allocatable:: Csum (:,:) + real(r8),allocatable:: Cvec (:) + real(r8),allocatable:: Bsum (:,:) + real(r8),allocatable:: Bnorm(:) + real(r8),allocatable:: Bcov (:,:) + real(r8):: area(pcols),rlat,rlon + + real(r8),allocatable:: Nr (:) + real(r8),allocatable:: Nr_sum(:) + real(r8),allocatable:: Rs (:) + real(r8),allocatable:: As (:) + real(r8):: Rc,Dth,Ws + real(r8):: Slat0,Slon0,X0,Y0,Z0 + real(r8):: Slon, Slat,Xs,Ys,Zs + real(r8):: Rot(3,3),Rse,Bavg,Pnm + integer :: kk,ns,Nsamp + + integer :: mm,nn,n2,nb,lchnk,ncols,cc + integer :: cnum,Cvec_len + + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'init_SphericalHarmonic_COV' + + if (I_nmax<1) then + call endrun('SphericalHarmonic_COV%init: ERROR I_nmax must be greater than 0') + end if + + ! Allocate space + !----------------- + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + if(allocated(this%bindx)) deallocate(this%bindx) + + O_nbas = (I_nmax+1)**2 + this%nmax = I_nmax + this%nbas = O_nbas + allocate(this%area (pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%area') + allocate(this%basis(pcols,begchunk:endchunk,O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%basis') + allocate(this%map (O_nbas,O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%map') + allocate(this%bindx(0:O_nbas,-O_nbas:O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%bindx') + this%area (:,:) = 0._r8 + this%basis(:,:,:) = 0._r8 + this%map (:,:) = 0._r8 + this%bindx(:,:) = -1 + + + Cvec_len = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + Cvec_len = Cvec_len + 1 + end do + end do + + nlcols = get_nlcols_p() + + allocate(Clons(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clons') + allocate(Clats(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(Bcoef(O_nbas/2+1), stat=astat) + call handle_allocate_error(astat, subname, 'Bcoef') + allocate(Csum (nlcols, Cvec_len), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + allocate(Cvec (Cvec_len), stat=astat) + call handle_allocate_error(astat, subname, 'Cvec') + allocate(Bsum (nlcols, O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bsum') + allocate(Bnorm(O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bnorm') + allocate(Bcov (O_nbas,O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bcov') + + Bsum(:,:) = 0._r8 + Csum(:,:) = 0._r8 + + ! Save a copy of the area weights for each ncol gridpoint + ! and convert Latitudes to SP->NP colatitudes in radians + !------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + call get_wght_all_p(lchnk, ncols, area) + do cc = 1,ncols + rlat=get_rlat_p(lchnk,cc) + rlon=get_rlon_p(lchnk,cc) + this%area(cc,lchnk) = area(cc) + Clons (cc,lchnk) = rlon + Clats (cc,lchnk) = rlat + end do + end do + + ! Initialize Samping Grid Size + !----------------------------- + if(present(SAMPLE_NRING)) then + this%nring = SAMPLE_NRING + else + this%nring = 1 + endif + + if(this%nring < 1) then + call endrun('SphericalHarmonic_COV%init: ERROR nring be 1 or more') + endif + + ! Set the number of gridpoints in each concentric ring + !----------------------------------------------------- + allocate(Nr (this%nring)) + allocate(Nr_sum(this%nring)) + Nr (1) = 1 + Nr_sum(1) = 1 + do kk=2,this%nring + Nr (kk) = 8*(kk-1) + Nr_sum(kk) = Nr_sum(kk-1) + Nr(kk) + end do + Nsamp = Nr_sum(this%nring) + + ! Init the equal area grid with the center point, + ! then add the polar coordinate gridpoints for + ! each ring for a reference domain radius=1. + !------------------------------------------------- + allocate(Rs(Nsamp)) + allocate(As(Nsamp)) + Rs(1) = 0._r8 + As(1) = 0._r8 + ns = 1 + do kk =2,this%nring + Rc = (sqrt(real(Nr_sum(kk)))+sqrt(real(Nr_sum(kk-1))))/(2._r8*sqrt(real(Nsamp))) + Dth = twoPI/real(Nr(kk)) + do nn=1,Nr(kk) + ns = ns + 1 + Rs(ns) = Rc + As(ns) = Dth*(nn-1) + end do + end do + Ws = 1._r8/real(Nsamp) + + if(ns.ne.Nsamp) then + call endrun('SphericalHarmonic_COV%init: ERROR mismatch in sample point number') + endif + + ! Samping INFO/DIAG Output + !-------------------------- + IF((.TRUE.).and.(Nsamp>1)) THEN + write(iulog,*) ' ' + write(iulog,*) 'SphericalHarmonic_COV%init: SAMPLE GRID: Nsamp=',Nsamp + do ns=1,Nsamp + write(iulog,*) ' ns=',ns,' Grid: Rs=',Rs(ns),' As=',As(ns),As(ns)*360./twoPI + end do + write(iulog,*) ' ' + ENDIF + + ! Add first basis for the mean values. + !------------------------------------------ + this%nbas = 1 + this%bindx(0,0) = 1 + this%basis(:,begchunk:endchunk,1) = invSqrt4pi + + ! Loop over the remaining meridional modes + !------------------------------------------ + do nb=2,(this%nmax+1) + nn = nb-1 + + ! Loop over zonal modes for each nn + !----------------------------------- + do mm=-nn,nn + this%nbas = this%nbas + 1 + this%bindx(nn,mm) = this%nbas + + ! Generate coefs for the basis + !------------------------------ + call sh_gen_basis_coefs(nn,iabs(mm),Bcoef) + + if(this%nring.eq.1) then + ! Create basis for the coefs at each ncol gridpoint + ! Lats shifted by PI/2 for SH's SP origin + !--------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + Slat = Clats(cc,lchnk) + halfPI + call sh_create_basis(nn,iabs(mm),Slat,Bcoef,this%basis(cc,lchnk,this%nbas)) + if(mm < 0) then + this%basis(cc,lchnk,this%nbas) = this%basis(cc,lchnk,this%nbas)*sin(mm*Clons(cc,lchnk)) + elseif(mm > 0) then + this%basis(cc,lchnk,this%nbas) = this%basis(cc,lchnk,this%nbas)*cos(mm*Clons(cc,lchnk)) + endif + end do + end do + else + ! Optionally Compute local area averge of + ! basis values for each gridpoiont + !--------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + + ! Set values to rotate/scale the sample + ! grid to the current SE gridpoint + !----------------------------------------- + Rse = sqrt(this%area(cc,lchnk)/pi) + Rot(1,1) = -sin(Clons(cc,lchnk)) + Rot(2,1) = -cos(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,1) = cos(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,2) = cos(Clons(cc,lchnk)) + Rot(2,2) = -sin(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,2) = sin(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,3) = 0._r8 + Rot(2,3) = cos(Clats(cc,lchnk)) + Rot(3,3) = sin(Clats(cc,lchnk)) + + ! For each SE grid point, we rotate/scale the reference + ! grid for a local sampling domain with the elements area. + ! The computed basis values are then averaged over this domain + !------------------------------------------------------------ + Bavg = 0._r8 + do ns=1,Nsamp + ! scale polar coordinates to lat/lon circle centered at NP + !--------------------------------------------------------- + Slat0 = halfPI - Rse*Rs(ns) + Slon0 = As(ns) + + ! Compute cartesian coordinates + !------------------------------- + X0 = cos(Slon0)*cos(Slat0) + Y0 = sin(Slon0)*cos(Slat0) + Z0 = sin(Slat0) + + ! Apply the Rotation from the NP + !-------------------------------- + Xs = X0*Rot(1,1) + Y0*Rot(2,1) + Z0*Rot(3,1) + Ys = X0*Rot(1,2) + Y0*Rot(2,2) + Z0*Rot(3,2) + Zs = X0*Rot(1,3) + Y0*Rot(2,3) + Z0*Rot(3,3) + + ! Compute resulting lat/lon gridpoint + !------------------------------------- + Slon = atan2(Ys,Xs) + Slat = halfPI + asin(Zs) + + ! Compute basis value and add the result to the average + ! Slat shifted by PI/2 for SH's SP origin + !------------------------------------------------------- + Slat = Slat + halfPI + call sh_create_basis(nn,iabs(mm),Slat,Bcoef,Pnm) + if(mm < 0) then + Bavg = Bavg + Ws*Pnm*sin(mm*Slon) + elseif(mm > 0) then + Bavg = Bavg + Ws*Pnm*cos(mm*Slon) + endif + end do + this%basis(cc,lchnk,this%nbas) = Bavg + end do + end do + endif + + end do ! mm=-nn,nn + end do ! nn=2,this%nbas + O_nbas = this%nbas + + ! Numerically normalize the basis funnctions + !-------------------------------------------------------------- + do nn=1,this%nbas + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Bsum(count,nn) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do ! nn=1,this%nbas + + call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) + + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(nn)) + end do + end do ! nn=1,this%nbas + + ! Compute covariance matrix for basis functions + ! (Yes, they are theoretically orthonormal, but lets make sure) + !--------------------------------------------------------------- + cnum = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + cnum = cnum + 1 + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,cnum) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk) + end do + end do + + end do + end do + + call shr_reprosum_calc(Csum, Cvec, count, nlcols, Cvec_len, gbl_count=ngcols_p, commid=mpicom) + + cnum = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + cnum = cnum + 1 + Bcov(nn,n2) = Cvec(cnum) + Bcov(n2,nn) = Cvec(cnum) + end do + end do + + ! Invert to get the basis amplitude map + !-------------------------------------- + call Invert_Matrix(Bcov,this%nbas,this%map) + + ! End Routine + !------------ + deallocate(Nr) + deallocate(Nr_sum) + deallocate(Rs) + deallocate(As) + deallocate(Clons) + deallocate(Clats) + deallocate(Bcoef) + deallocate(Csum ) + deallocate(Cvec ) + deallocate(Bsum ) + deallocate(Bnorm) + deallocate(Bcov ) + end subroutine init_SphericalHarmonic_COV + !======================================================================= + + + !======================================================================= + subroutine final_SphericalHarmonic_COV(this) + class(SphericalHarmonic_COV_t) :: this + + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + + end subroutine final_SphericalHarmonic_COV + !======================================================================= + + + !======================================================================= + subroutine calc_SphericalHarmonic_2Damps_COV(this,I_Gdata,O_Bamp) + ! + ! calc_SphericalHarmonic_2Damps: Given 2D data values for the ncol gridpoints, + ! compute the spherical harmonic basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_COV_t) :: this + real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk) + real(r8),intent(out) :: O_Bamp(:) + ! + ! Local Values + !-------------- + real(r8),allocatable :: Csum(:,:) + real(r8),allocatable :: Gcov(:) + integer :: nn,n2,ncols,lchnk,cc + integer :: nlcols, count, astat + + character(len=*), parameter :: subname = 'calc_SphericalHarmonic_2Damps_COV' + + nlcols = get_nlcols_p() + + allocate(Gcov(this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + allocate(Csum(nlcols, this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + Csum(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do nn= 1,this%nbas + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + + call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do nn=1,this%nbas + O_Bamp(nn) = 0._r8 + do n2=1,this%nbas + O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2) + end do + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Gcov) + + end subroutine calc_SphericalHarmonic_2Damps_COV + !======================================================================= + + + !======================================================================= + subroutine calc_SphericalHarmonic_3Damps_COV(this,I_Gdata,O_Bamp) + ! + ! calc_SphericalHarmonic_3Damps: Given 3D data values for the ncol,nlev gridpoints, + ! compute the zonal mean basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_COV_t) :: this + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) + real(r8),intent(out):: O_Bamp (:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Csum (:,:) + real(r8),allocatable:: Gcov (:) + integer:: nn,n2,ncols,lchnk,cc + integer:: Nsum,ns,ll,idx + integer :: nlcols, count, astat + + integer :: nlev + integer :: nbas_x_nlev + character(len=*), parameter :: subname = 'calc_SphericalHarmonic_3Damps_COV' + + nlev = size(I_Gdata,dim=2) + nlcols = get_nlcols_p() + + ! Flatten the vertical and basis dimensions to + ! allow a single MPI reduction + !----------------------------------------------- + nbas_x_nlev = this%nbas*nlev + allocate(Gcov(nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + allocate(Csum(nlcols, nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + + Csum(:,:) = 0._r8 + O_Bamp(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do ll= 1,nlev + do nn= 1,this%nbas + idx = (ll-1)*this%nbas + nn + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,idx) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + end do + + call shr_reprosum_calc(Csum, Gcov, count, nlcols, nbas_x_nlev, gbl_count=ngcols_p, commid=mpicom) + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do ll= 1,nlev + do nn=1,this%nbas + O_Bamp(nn,ll) = 0._r8 + do n2=1,this%nbas + idx = (ll-1)*this%nbas + n2 + O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(idx) + end do + end do + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Gcov) + + end subroutine calc_SphericalHarmonic_3Damps_COV + !======================================================================= + + + !======================================================================= + subroutine eval_SphericalHarmonic_2Dgrid_COV(this,I_Bamp,O_Gdata) + ! + ! eval_SphericalHarmonic_2Dgrid: Given the zonal mean basis amplitudes, + ! compute 2D data values for the ncol gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_COV_t) :: this + real(r8),intent(in ):: I_Bamp (:) + real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + + O_Gdata(:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn)) + end do + end do + end do + + end subroutine eval_SphericalHarmonic_2Dgrid_COV + !======================================================================= + + + !======================================================================= + subroutine eval_SphericalHarmonic_3Dgrid_COV(this,I_Bamp,O_Gdata) + ! + ! eval_SphericalHarmonic_3Dgrid: Given the zonal mean basis amplitudes, + ! compute 3D data values for the ncol,nlev gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_COV_t) :: this + real(r8),intent(in ):: I_Bamp (:,:) + real(r8),intent(out):: O_Gdata(:,:,begchunk:) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + integer:: ll + + integer :: nlev + nlev = size(O_Gdata,dim=2) + + O_Gdata(:,:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do ll = 1,nlev + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn)) + end do + end do + end do + end do + + end subroutine eval_SphericalHarmonic_3Dgrid_COV + !======================================================================= + + + !======================================================================= + subroutine init_SphericalHarmonic_GS(this,I_nmax,O_nbas,SAMPLE_NRING) + ! + ! init_SphericalHarmonic: Initialize the SphericalHarmonic data structure + ! for the physics grid. It is assumed that the domain + ! of these gridpoints spans the surface of the sphere. + ! The representation of basis functions is + ! normalized w.r.t integration over the sphere. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_GS_t) :: this + integer ,intent(in ):: I_nmax + integer ,intent(out):: O_nbas + integer ,intent(in ),optional:: SAMPLE_NRING + ! + ! Local Values + !-------------- + real(r8),allocatable:: Clons(:,:) + real(r8),allocatable:: Clats(:,:) + real(r8),allocatable:: Bcoef(:) + real(r8),allocatable:: Bsum (:,:) + real(r8),allocatable:: Bnorm (:) + real(r8),allocatable:: Bamp (:) + real(r8):: area(pcols),rlat,rlon + real(r8):: Pnm + + real(r8),allocatable:: Nr (:) + real(r8),allocatable:: Nr_sum(:) + real(r8),allocatable:: Rs (:) + real(r8),allocatable:: As (:) + real(r8):: Rc,Dth,Ws + real(r8):: Slat0,Slon0,X0,Y0,Z0 + real(r8):: Slon, Slat,Xs,Ys,Zs + real(r8):: Rot(3,3),Rse,Bavg,Bavg_e,Bavg_o + integer :: kk,ns,Nsamp + + integer :: nbas_e,nbas_o + integer :: mm,nn,n2,nb,lchnk,ncols,cc + + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'init_SphericalHarmonic_GS' + + if (I_nmax<1) then + call endrun('SphericalHarmonic_GS%init: ERROR I_nmax must be greater than 0') + end if + + ! Allocate space + !----------------- + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + + O_nbas = (I_nmax+1)**2 + this%nmax = I_nmax + this%nbas = O_nbas + allocate(this%area (pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%area') + allocate(this%basis(pcols,begchunk:endchunk,O_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%basis') + this%area (:,:) = 0._r8 + this%basis(:,:,:) = 0._r8 + + nlcols = get_nlcols_p() + + allocate(Clons(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clons') + allocate(Clats(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(Bcoef(O_nbas/2+1), stat=astat) + call handle_allocate_error(astat, subname, 'Bcoef') + allocate(Bsum (nlcols,1), stat=astat) + call handle_allocate_error(astat, subname, 'Bsum') + allocate(Bamp (1), stat=astat) + call handle_allocate_error(astat, subname, 'Bamp') + allocate(Bnorm(1), stat=astat) + call handle_allocate_error(astat, subname, 'Bnorm') + + ! Save a copy of the area weights for each ncol gridpoint + ! and convert Latitudes to SP->NP colatitudes in radians + !------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + call get_wght_all_p(lchnk, ncols, area) + do cc = 1,ncols + rlat=get_rlat_p(lchnk,cc) + rlon=get_rlon_p(lchnk,cc) + this%area(cc,lchnk) = area(cc) + Clons (cc,lchnk) = rlon + Clats (cc,lchnk) = rlat + end do + end do + + ! Initialize Samping Grid Size + !----------------------------- + if(present(SAMPLE_NRING)) then + this%nring = SAMPLE_NRING + else + this%nring = 1 + endif + + if(this%nring < 1) then + call endrun('SphericalHarmonic_GS%init: ERROR nring be 1 or more') + endif + + ! Set the number of gridpoints in each concentric ring + !----------------------------------------------------- + allocate(Nr (this%nring)) + allocate(Nr_sum(this%nring)) + Nr (1) = 1 + Nr_sum(1) = 1 + do kk=2,this%nring + Nr (kk) = 8*(kk-1) + Nr_sum(kk) = Nr_sum(kk-1) + Nr(kk) + end do + Nsamp = Nr_sum(this%nring) + + ! Init the equal area grid with the center point, + ! then add the polar coordinate gridpoints for + ! each ring for a reference domain radius=1. + !------------------------------------------------- + allocate(Rs(Nsamp)) + allocate(As(Nsamp)) + Rs(1) = 0._r8 + As(1) = 0._r8 + ns = 1 + do kk =2,this%nring + Rc = (sqrt(real(Nr_sum(kk)))+sqrt(real(Nr_sum(kk-1))))/(2._r8*sqrt(real(Nsamp))) + Dth = twoPI/real(Nr(kk)) + do nn=1,Nr(kk) + ns = ns + 1 + Rs(ns) = Rc + As(ns) = Dth*(nn-1) + end do + end do + Ws = 1._r8/real(Nsamp) + + if(ns.ne.Nsamp) then + call endrun('SphericalHarmonic_GS%init: ERROR mismatch in sample point number') + endif + + ! Samping INFO/DIAG Output + !-------------------------- + IF((.TRUE.).and.(Nsamp>1)) THEN + write(iulog,*) ' ' + write(iulog,*) 'SphericalHarmonic_GS%init: SAMPLE GRID: Nsamp=',Nsamp + do ns=1,Nsamp + write(iulog,*) ' ns=',ns,' Grid: Rs=',Rs(ns),' As=',As(ns),As(ns)*360./twoPI + end do + write(iulog,*) ' ' + ENDIF + + ! Add first basis for the mean values. + !------------------------------------------ + this%nbas = 1 + this%basis(:,begchunk:endchunk,1) = invSqrt4pi + + ! Loop over the remaining meridional modes, + ! The ordering of the basis functions here is + ! set to accomodate Gram-Schmidt + !------------------------------------------------ + do nb=2,(this%nmax+1) + nn = nb-1 + + ! Add the m=0 mode first + !------------------------ + this%nbas = this%nbas + 1 + call sh_gen_basis_coefs(nn,0,Bcoef) + + if(this%nring.eq.1) then + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + ! Lats shifted by PI/2 for SH's SP origin + !-------------------------------------------- + Slat = Clats(cc,lchnk) + halfPI + call sh_create_basis(nn,0,Slat,Bcoef,this%basis(cc,lchnk,this%nbas)) + end do + end do + else + ! Optionally Compute local area averge of + ! basis values for each gridpoiont + !--------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + + ! Set values to rotate/scale the sample + ! grid to the current SE gridpoint + !----------------------------------------- + Rse = sqrt(this%area(cc,lchnk)/pi) + Rot(1,1) = -sin(Clons(cc,lchnk)) + Rot(2,1) = -cos(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,1) = cos(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,2) = cos(Clons(cc,lchnk)) + Rot(2,2) = -sin(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,2) = sin(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,3) = 0._r8 + Rot(2,3) = cos(Clats(cc,lchnk)) + Rot(3,3) = sin(Clats(cc,lchnk)) + + ! For each SE grid point, we rotate/scale the reference + ! grid for a local sampling domain with the elements area. + ! The computed basis values are then averaged over this domain + !------------------------------------------------------------ + Bavg = 0._r8 + do ns=1,Nsamp + ! scale polar coordinates to lat/lon circle centered at NP + !--------------------------------------------------------- + Slat0 = halfPI - Rse*Rs(ns) + Slon0 = As(ns) + + ! Compute cartesian coordinates + !------------------------------- + X0 = cos(Slon0)*cos(Slat0) + Y0 = sin(Slon0)*cos(Slat0) + Z0 = sin(Slat0) + + ! Apply the Rotation from the NP + !-------------------------------- + Xs = X0*Rot(1,1) + Y0*Rot(2,1) + Z0*Rot(3,1) + Ys = X0*Rot(1,2) + Y0*Rot(2,2) + Z0*Rot(3,2) + Zs = X0*Rot(1,3) + Y0*Rot(2,3) + Z0*Rot(3,3) + + ! Compute resulting lat/lon gridpoint + !------------------------------------- + Slon = atan2(Ys,Xs) + Slat = halfPI + asin(Zs) + + ! Compute basis value and add the result to the average + ! Slat shifted by PI/2 for SH's SP origin + !------------------------------------------------------- + Slat = Slat + halfPI + call sh_create_basis(nn,0,Slat,Bcoef,Pnm) + Bavg = Bavg + Ws*Pnm + end do + this%basis(cc,lchnk,this%nbas) = Bavg + end do + end do + endif + + ! Now loop over zonal modes mm=1,nn + ! and add even/odd basis functions + !---------------------------------- + do mm=1,nn + nbas_e = this%nbas + 1 + nbas_o = this%nbas + 2 + this%nbas = this%nbas + 2 + call sh_gen_basis_coefs(nn,mm,Bcoef) + + if(this%nring.eq.1) then + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + ! Lats shifted by PI/2 for SH's SP origin + !-------------------------------------------- + Slat = Clats(cc,lchnk) + halfPI + call sh_create_basis(nn,mm,Slat,Bcoef,Pnm) + this%basis(cc,lchnk,nbas_e) = Pnm*cos(mm*Clons(cc,lchnk)) + this%basis(cc,lchnk,nbas_o) = Pnm*sin(mm*Clons(cc,lchnk)) + end do + end do + else + ! Optionally Compute local area averge of + ! basis values for each gridpoiont + !--------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + + ! Set values to rotate/scale the sample + ! grid to the current SE gridpoint + !----------------------------------------- + Rse = sqrt(this%area(cc,lchnk)/pi) + Rot(1,1) = -sin(Clons(cc,lchnk)) + Rot(2,1) = -cos(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,1) = cos(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,2) = cos(Clons(cc,lchnk)) + Rot(2,2) = -sin(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,2) = sin(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,3) = 0._r8 + Rot(2,3) = cos(Clats(cc,lchnk)) + Rot(3,3) = sin(Clats(cc,lchnk)) + + ! For each SE grid point, we rotate/scale the reference + ! grid for a local sampling domain with the elements area. + ! The computed basis values are then averaged over this domain + !------------------------------------------------------------ + Bavg_e = 0._r8 + Bavg_o = 0._r8 + do ns=1,Nsamp + ! scale polar coordinates to lat/lon circle centered at NP + !--------------------------------------------------------- + Slat0 = halfPI - Rse*Rs(ns) + Slon0 = As(ns) + + ! Compute cartesian coordinates + !------------------------------- + X0 = cos(Slon0)*cos(Slat0) + Y0 = sin(Slon0)*cos(Slat0) + Z0 = sin(Slat0) + + ! Apply the Rotation from the NP + !-------------------------------- + Xs = X0*Rot(1,1) + Y0*Rot(2,1) + Z0*Rot(3,1) + Ys = X0*Rot(1,2) + Y0*Rot(2,2) + Z0*Rot(3,2) + Zs = X0*Rot(1,3) + Y0*Rot(2,3) + Z0*Rot(3,3) + + ! Compute resulting lat/lon gridpoint + !------------------------------------- + Slon = atan2(Ys,Xs) + Slat = halfPI + asin(Zs) + + ! Compute basis value and add the result to the average + ! Slat shifted by PI/2 for SH's SP origin + !------------------------------------------------------- + Slat = Slat + halfPI + call sh_create_basis(nn,mm,Slat,Bcoef,Pnm) + Bavg_o = Bavg_o + Ws*Pnm*sin(mm*Slon) + Bavg_e = Bavg_e + Ws*Pnm*cos(mm*Slon) + end do + this%basis(cc,lchnk,nbas_e) = Bavg_e + this%basis(cc,lchnk,nbas_o) = Bavg_o + end do + end do + endif + + end do + end do ! nn=2,this%nbas + O_nbas = this%nbas + + !------------------------------------------------------------- + ! The Discrete basis representation needs to be orthogonalized + ! Grahm-Schmidt Orthogonization + !------------------------------------------------------------- + + ! Numerically normalize the gravest functon + ! as the first orthonormal basis function. + !------------------------------------------- + count = 0 + Bsum(:,:) = 0._r8 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count = count+1 + Bsum(count,1) = this%basis(cc,lchnk,1)*this%basis(cc,lchnk,1)*this%area(cc,lchnk) + end do + end do + + call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, 1, gbl_count=ngcols_p, commid=mpicom) + + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + this%basis(:ncols,lchnk,1) = this%basis(:ncols,lchnk,1)/sqrt(Bnorm(1)) + end do + + ! Loop over the remaining basis functions + !------------------------------------------- + do nn=2,this%nbas + + ! Remove contributions from exisiting set of orthonormal functions + !------------------------------------------------------------------ + do n2=1,(nn-1) + count = 0 + Bsum(:,:) = 0._r8 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count = count+1 + Bsum(count,1) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk) + end do + end do + + call shr_reprosum_calc(Bsum, Bamp, count, nlcols, 1, gbl_count=ngcols_p, commid=mpicom) + + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn) - Bamp(1)*this%basis(:ncols,lchnk,n2) + end do + end do ! n2=1,(nn-1) + + ! Normalize the result for the newest member of the orthonomal set + !-------------------------------------------------------------------- + count = 0 + Bsum(:,:) = 0._r8 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count = count+1 + Bsum(count,1) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + + call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, 1, gbl_count=ngcols_p, commid=mpicom) + + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(1)) + end do + end do ! nn=2,this%nbas + + !DIAG: Check for blatent orthogonality errors + !----------------------------------------------- +! IF(.TRUE.) THEN +! write(iulog,*) 'PFC: ORTHONORM CHECK:' +! do nn=2,this%nbas +! do n2=1,(nn-1) +! count = 0 +! Bsum(:,:) = 0._r8 +! do lchnk=begchunk,endchunk +! ncols = get_ncols_p(lchnk) +! do cc = 1,ncols +! count = count+1 +! Bsum(count,1) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk) +! end do +! end do +! call shr_reprosum_calc(Bsum, Bamp, count, nlcols, 1, gbl_count=ngcols_p, commid=mpicom) +! if(abs(Bamp(1)).gt.1.d-5) then +! write(iulog,*) 'PFC: *** nn=',nn,' n2=',n2,' Bamp=',Bamp(1) +! endif +! end do +! end do ! nn=2,this%nbas +! write(iulog,*) 'PFC: ORTHONORM CHECK: done' +! ENDIF + + ! End Routine + !------------ + deallocate(Nr) + deallocate(Nr_sum) + deallocate(Rs) + deallocate(As) + deallocate(Clons) + deallocate(Clats) + deallocate(Bcoef) + deallocate(Bsum ) + deallocate(Bamp ) + deallocate(Bnorm) + end subroutine init_SphericalHarmonic_GS + !======================================================================= + + + !======================================================================= + subroutine final_SphericalHarmonic_GS(this) + class(SphericalHarmonic_GS_t) :: this + + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + + end subroutine final_SphericalHarmonic_GS + !======================================================================= + + + !======================================================================= + subroutine calc_SphericalHarmonic_2Damps_GS(this,I_Gdata,O_Bamp) + ! + ! calc_SphericalHarmonic_2Damps: Given 2D data values for the ncol gridpoints, + ! compute the spherical harmonic basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_GS_t) :: this + real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk) + real(r8),intent(out) :: O_Bamp(:) + ! + ! Local Values + !-------------- + real(r8),allocatable :: Csum(:,:) + real(r8),allocatable :: Bamp(:) + integer :: nn,n2,ncols,lchnk,cc + integer :: nlcols, count, astat + + character(len=*), parameter :: subname = 'calc_SphericalHarmonic_2Damps_GS' + + nlcols = get_nlcols_p() + + allocate(Bamp(this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bamp') + allocate(Csum(nlcols, this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + Csum(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do nn= 1,this%nbas + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + + call shr_reprosum_calc(Csum, Bamp, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) + + ! Output the amplitudes + !-------------------------- + do nn=1,this%nbas + O_Bamp(nn) = Bamp(nn) + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Bamp) + + end subroutine calc_SphericalHarmonic_2Damps_GS + !======================================================================= + + + !======================================================================= + subroutine calc_SphericalHarmonic_3Damps_GS(this,I_Gdata,O_Bamp) + ! + ! calc_SphericalHarmonic_3Damps: Given 3D data values for the ncol,nlev gridpoints, + ! compute the zonal mean basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_GS_t) :: this + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) + real(r8),intent(out):: O_Bamp (:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Csum (:,:) + real(r8),allocatable:: Bamp (:) + integer:: nn,n2,ncols,lchnk,cc + integer:: Nsum,ns,ll,idx + integer :: nlcols, count, astat + + integer :: nlev + integer :: nbas_x_nlev + character(len=*), parameter :: subname = 'calc_SphericalHarmonic_3Damps_GS' + + nlev = size(I_Gdata,dim=2) + nlcols = get_nlcols_p() + + ! Flatten the vertical and basis dimensions to + ! allow a single MPI reduction + !----------------------------------------------- + nbas_x_nlev = this%nbas*nlev + allocate(Bamp(nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Bamp') + allocate(Csum(nlcols, nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + + Csum(:,:) = 0._r8 + O_Bamp(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do ll= 1,nlev + do nn= 1,this%nbas + idx = (ll-1)*this%nbas + nn + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,idx) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + end do + + call shr_reprosum_calc(Csum, Bamp, count, nlcols, nbas_x_nlev, gbl_count=ngcols_p, commid=mpicom) + + ! Unpack the amplitudes into the Output array + !-------------------------------------------- + do ll= 1,nlev + do nn=1,this%nbas + idx = (ll-1)*this%nbas + nn + O_Bamp(nn,ll) = Bamp(idx) + end do + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Bamp) + + end subroutine calc_SphericalHarmonic_3Damps_GS + !======================================================================= + + + !======================================================================= + subroutine eval_SphericalHarmonic_2Dgrid_GS(this,I_Bamp,O_Gdata) + ! + ! eval_SphericalHarmonic_2Dgrid: Given the zonal mean basis amplitudes, + ! compute 2D data values for the ncol gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_GS_t) :: this + real(r8),intent(in ):: I_Bamp (:) + real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + + O_Gdata(:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn)) + end do + end do + end do + + end subroutine eval_SphericalHarmonic_2Dgrid_GS + !======================================================================= + + + !======================================================================= + subroutine eval_SphericalHarmonic_3Dgrid_GS(this,I_Bamp,O_Gdata) + ! + ! eval_SphericalHarmonic_3Dgrid: Given the zonal mean basis amplitudes, + ! compute 3D data values for the ncol,nlev gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(SphericalHarmonic_GS_t) :: this + real(r8),intent(in ):: I_Bamp (:,:) + real(r8),intent(out):: O_Gdata(:,:,begchunk:) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + integer:: ll + + integer :: nlev + nlev = size(O_Gdata,dim=2) + + O_Gdata(:,:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do ll = 1,nlev + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn)) + end do + end do + end do + end do + + end subroutine eval_SphericalHarmonic_3Dgrid_GS + !======================================================================= + + + !======================================================================= + subroutine init_ZonalMean(this,I_nbas,SAMPLE_NRING) + ! + ! init_ZonalMean: Initialize the ZonalMean data structures for the + ! physics grid. It is assumed that the domain + ! of these gridpoints spans the surface of the sphere. + ! The representation of basis functions is + ! normalized w.r.t integration over the sphere. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalMean_t) :: this + integer ,intent(in):: I_nbas + integer ,intent(in),optional:: SAMPLE_NRING + ! + ! Local Values + !-------------- + real(r8),allocatable:: Clons(:,:) + real(r8),allocatable:: Clats(:,:) + real(r8),allocatable:: Bcoef(:) + real(r8),allocatable:: Csum (:,:) + real(r8),allocatable:: Cvec (:) + real(r8),allocatable:: Bsum (:,:) + real(r8),allocatable:: Bnorm(:) + real(r8),allocatable:: Bcov (:,:) + real(r8):: area(pcols),rlat,rlon + + real(r8),allocatable:: Nr (:) + real(r8),allocatable:: Nr_sum(:) + real(r8),allocatable:: Rs (:) + real(r8),allocatable:: As (:) + real(r8):: Rc,Dth,Ws + real(r8):: Slat0,Slon0,X0,Y0,Z0 + real(r8):: Slon, Slat,Xs,Ys,Zs + real(r8):: Rot(3,3),Rse,Bavg,Pnm + integer :: kk,ns,mm,Nsamp + + integer :: nn,n2,nb,lchnk,ncols,cc + integer :: cnum,Cvec_len + + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'init_ZonalMean' + + if (I_nbas<1) then + call endrun('ZonalMean%init: ERROR I_nbas must be greater than 0') + end if + + ! Allocate space + !----------------- + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + + this%nbas = I_nbas + allocate(this%area (pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%area') + allocate(this%basis(pcols,begchunk:endchunk,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%basis') + allocate(this%map (I_nbas,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%map') + this%area (:,:) = 0._r8 + this%basis(:,:,:) = 0._r8 + this%map (:,:) = 0._r8 + + Cvec_len = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + Cvec_len = Cvec_len + 1 + end do + end do + + nlcols = get_nlcols_p() + + allocate(Clons(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clons') + allocate(Clats(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(Bcoef(I_nbas/2+1), stat=astat) + call handle_allocate_error(astat, subname, 'Bcoef') + allocate(Csum (nlcols, Cvec_len), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + allocate(Cvec (Cvec_len), stat=astat) + call handle_allocate_error(astat, subname, 'Cvec') + allocate(Bsum (nlcols, I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bsum') + allocate(Bnorm(I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bnorm') + allocate(Bcov (I_nbas,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bcov') + + Bsum(:,:) = 0._r8 + Csum(:,:) = 0._r8 + + ! Save a copy of the area weights for each ncol gridpoint + ! and convert Latitudes to SP->NP colatitudes in radians + !------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + call get_wght_all_p(lchnk, ncols, area) + do cc = 1,ncols + rlat=get_rlat_p(lchnk,cc) + rlon=get_rlon_p(lchnk,cc) + this%area(cc,lchnk) = area(cc) + Clats (cc,lchnk) = rlat + Clons (cc,lchnk) = rlon + end do + end do + + ! Initialize Samping Grid Size + !----------------------------- + if(present(SAMPLE_NRING)) then + this%nring = SAMPLE_NRING + else + this%nring = 1 + endif + + if(this%nring < 1) then + call endrun('ZonalMean%init: ERROR nring be 1 or more') + endif + + ! Set the number of gridpoints in each concentric ring + !----------------------------------------------------- + allocate(Nr (this%nring)) + allocate(Nr_sum(this%nring)) + Nr (1) = 1 + Nr_sum(1) = 1 + do kk=2,this%nring + Nr (kk) = 8*(kk-1) + Nr_sum(kk) = Nr_sum(kk-1) + Nr(kk) + end do + Nsamp = Nr_sum(this%nring) + + ! Init the equal area grid with the center point, + ! then add the polar coordinate gridpoints for + ! each ring for a reference domain radius=1. + !------------------------------------------------- + allocate(Rs(Nsamp)) + allocate(As(Nsamp)) + Rs(1) = 0._r8 + As(1) = 0._r8 + ns = 1 + do kk =2,this%nring + Rc = (sqrt(real(Nr_sum(kk)))+sqrt(real(Nr_sum(kk-1))))/(2._r8*sqrt(real(Nsamp))) + Dth = twoPI/real(Nr(kk)) + do nn=1,Nr(kk) + ns = ns + 1 + Rs(ns) = Rc + As(ns) = Dth*(nn-1) + end do + end do + Ws = 1._r8/real(Nsamp) + + if(ns.ne.Nsamp) then + call endrun('ZonalMean%init: ERROR mismatch in sample point number') + endif + + ! Samping INFO/DIAG Output + !-------------------------- + IF((.TRUE.).and.(Nsamp>1)) THEN + write(iulog,*) ' ' + write(iulog,*) 'ZonalMean%init: SAMPLE GRID: Nsamp=',Nsamp + do ns=1,Nsamp + write(iulog,*) ' ns=',ns,' Grid: Rs=',Rs(ns),' As=',As(ns),As(ns)*360./twoPI + end do + write(iulog,*) ' ' + ENDIF + + ! Add first basis for the mean values. + !------------------------------------------ + this%basis(:,begchunk:endchunk,1) = invSqrt4pi + + ! Loop over the remaining basis functions + !--------------------------------------- + do nn=2,this%nbas + nb = nn-1 + + ! Generate coefs for the basis + !------------------------------ + call sh_gen_basis_coefs(nb,0,Bcoef) + + if(this%nring.eq.1) then + ! Create basis for the coefs at each ncol gridpoint + ! Lats shifted by PI/2 for SH's SP origin + !--------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + Slat = Clats(cc,lchnk) + halfPI + call sh_create_basis(nb,0,Slat,Bcoef,this%basis(cc,lchnk,nn)) + end do + end do + else + ! Optionally Compute local area averge of + ! basis values for each gridpoiont + !--------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + + ! Set values to rotate/scale the sample + ! grid to the current SE gridpoint + !----------------------------------------- + Rse = sqrt(this%area(cc,lchnk)/pi) + Rot(1,1) = -sin(Clons(cc,lchnk)) + Rot(2,1) = -cos(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,1) = cos(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,2) = cos(Clons(cc,lchnk)) + Rot(2,2) = -sin(Clons(cc,lchnk))*sin(Clats(cc,lchnk)) + Rot(3,2) = sin(Clons(cc,lchnk))*cos(Clats(cc,lchnk)) + Rot(1,3) = 0._r8 + Rot(2,3) = cos(Clats(cc,lchnk)) + Rot(3,3) = sin(Clats(cc,lchnk)) + + ! For each SE grid point, we rotate/scale the reference + ! grid for a local sampling domain with the elements area. + ! The computed basis values are then averaged over this domain + !------------------------------------------------------------ + Bavg = 0._r8 + do ns=1,Nsamp + ! scale polar coordinates to lat/lon circle centered at NP + !--------------------------------------------------------- + Slat0 = halfPI - Rse*Rs(ns) + Slon0 = As(ns) + + ! Compute cartesian coordinates + !------------------------------- + X0 = cos(Slon0)*cos(Slat0) + Y0 = sin(Slon0)*cos(Slat0) + Z0 = sin(Slat0) + + ! Apply the Rotation from the NP + !-------------------------------- + Xs = X0*Rot(1,1) + Y0*Rot(2,1) + Z0*Rot(3,1) + Ys = X0*Rot(1,2) + Y0*Rot(2,2) + Z0*Rot(3,2) + Zs = X0*Rot(1,3) + Y0*Rot(2,3) + Z0*Rot(3,3) + + ! Compute resulting lat/lon gridpoint + !------------------------------------- + Slon = atan2(Ys,Xs) + Slat = halfPI + asin(Zs) + + ! Compute basis value and add the result to the average + ! Slat shifted by PI/2 for SH's SP origin + !------------------------------------------------------- + Slat = Slat + halfPI + call sh_create_basis(nn,0,Slat,Bcoef,Pnm) + if(mm < 0) then + Bavg = Bavg + Ws*Pnm*sin(mm*Slon) + elseif(mm > 0) then + Bavg = Bavg + Ws*Pnm*cos(mm*Slon) + endif + end do + this%basis(cc,lchnk,this%nbas) = Bavg + end do + end do + endif + + end do ! nn=2,this%nbas + + ! Numerically normalize the basis funnctions + !-------------------------------------------------------------- + do nn=1,this%nbas + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Bsum(count,nn) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do ! nn=1,this%nbas + + call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) + + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(nn)) + end do + end do ! nn=1,this%nbas + + ! Compute covariance matrix for basis functions + ! (Yes, they are theoretically orthonormal, but lets make sure) + !--------------------------------------------------------------- + cnum = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + cnum = cnum + 1 + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,cnum) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk) + end do + end do + + end do + end do + + call shr_reprosum_calc(Csum, Cvec, count, nlcols, Cvec_len, gbl_count=ngcols_p, commid=mpicom) + + cnum = 0 + do nn= 1,this%nbas + do n2=nn,this%nbas + cnum = cnum + 1 + Bcov(nn,n2) = Cvec(cnum) + Bcov(n2,nn) = Cvec(cnum) + end do + end do + + ! Invert to get the basis amplitude map + !-------------------------------------- + call Invert_Matrix(Bcov,this%nbas,this%map) + + ! End Routine + !------------ + deallocate(Nr) + deallocate(Nr_sum) + deallocate(Rs) + deallocate(As) + deallocate(Clats) + deallocate(Bcoef) + deallocate(Csum ) + deallocate(Cvec ) + deallocate(Bsum ) + deallocate(Bnorm) + deallocate(Bcov ) + + end subroutine init_ZonalMean + !======================================================================= + + + !======================================================================= + subroutine final_ZonalMean(this) + class(ZonalMean_t) :: this + + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + + end subroutine final_ZonalMean + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalMean_2Damps(this,I_Gdata,O_Bamp) + ! + ! calc_ZonalMean_2Damps: Given 2D data values for the ncol gridpoints, + ! compute the zonal mean basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalMean_t) :: this + real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk) + real(r8),intent(out) :: O_Bamp(:) + ! + ! Local Values + !-------------- + real(r8),allocatable :: Csum(:,:) + real(r8),allocatable :: Gcov(:) + integer :: nn,n2,ncols,lchnk,cc + integer :: nlcols, count, astat + + character(len=*), parameter :: subname = 'calc_ZonalMean_2Damps' + + nlcols = get_nlcols_p() + + allocate(Gcov(this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + allocate(Csum(nlcols, this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + Csum(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do nn= 1,this%nbas + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + + call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do nn=1,this%nbas + O_Bamp(nn) = 0._r8 + do n2=1,this%nbas + O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2) + end do + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Gcov) + + end subroutine calc_ZonalMean_2Damps + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp) + ! + ! calc_ZonalMean_3Damps: Given 3D data values for the ncol,nlev gridpoints, + ! compute the zonal mean basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalMean_t) :: this + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) + real(r8),intent(out):: O_Bamp (:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Csum (:,:) + real(r8),allocatable:: Gcov (:) + integer:: nn,n2,ncols,lchnk,cc + integer:: Nsum,ns,ll,idx + integer :: nlcols, count, astat + + integer :: nlev + integer :: nbas_x_nlev + character(len=*), parameter :: subname = 'calc_ZonalMean_3Damps' + + nlev = size(I_Gdata,dim=2) + nlcols = get_nlcols_p() + + ! Flatten the vertical and basis dimensions to + ! allow a single MPI reduction + !----------------------------------------------- + nbas_x_nlev = this%nbas*nlev + allocate(Gcov(nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + allocate(Csum(nlcols, nbas_x_nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Csum') + + Csum(:,:) = 0._r8 + O_Bamp(:,:) = 0._r8 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + do ll= 1,nlev + do nn= 1,this%nbas + idx = (ll-1)*this%nbas + nn + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + count=count+1 + Csum(count,idx) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + end do + end do + end do + end do + + call shr_reprosum_calc(Csum, Gcov, count, nlcols, nbas_x_nlev, gbl_count=ngcols_p, commid=mpicom) + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do ll= 1,nlev + do nn=1,this%nbas + O_Bamp(nn,ll) = 0._r8 + do n2=1,this%nbas + idx = (ll-1)*this%nbas + n2 + O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(idx) + end do + end do + end do + + ! End Routine + !------------ + deallocate(Csum) + deallocate(Gcov) + + end subroutine calc_ZonalMean_3Damps + !======================================================================= + + + !======================================================================= + subroutine eval_ZonalMean_2Dgrid(this,I_Bamp,O_Gdata) + ! + ! eval_ZonalMean_2Dgrid: Given the zonal mean basis amplitudes, + ! compute 2D data values for the ncol gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalMean_t) :: this + real(r8),intent(in ):: I_Bamp (:) + real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + + O_Gdata(:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn)) + end do + end do + end do + + end subroutine eval_ZonalMean_2Dgrid + !======================================================================= + + + !======================================================================= + subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata) + ! + ! eval_ZonalMean_3Dgrid: Given the zonal mean basis amplitudes, + ! compute 3D data values for the ncol,nlev gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalMean_t) :: this + real(r8),intent(in ):: I_Bamp (:,:) + real(r8),intent(out):: O_Gdata(:,:,begchunk:) + ! + ! Local Values + !-------------- + integer:: nn,ncols,lchnk,cc + integer:: ll + + integer :: nlev + nlev = size(O_Gdata,dim=2) + + O_Gdata(:,:,:) = 0._r8 + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + + do ll = 1,nlev + do nn=1,this%nbas + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn)) + end do + end do + end do + end do + + end subroutine eval_ZonalMean_3Dgrid + !======================================================================= + + + !======================================================================= + subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) + ! + ! init_ZonalProfile: Initialize the ZonalProfile data structure for the + ! given nlat gridpoints. It is assumed that the domain + ! of these gridpoints of the profile span latitudes + ! from SP to NP. + ! The representation of basis functions functions is + ! normalized w.r.t integration over the sphere so that + ! when configured for tha same number of basis functions, + ! the calculated amplitudes are interchangable with + ! those for the ZonalMean_t class. + ! + ! The optional GEN_GAUSSLATS flag allows for the + ! generation of Gaussian latitudes. The generated grid + ! over-writes the values of IO_lats/IO_area passed by + ! the user. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalProfile_t) :: this + real(r8) ,intent(inout):: IO_lats(:) + real(r8) ,intent(inout):: IO_area(:) + integer ,intent(in):: I_nlat + integer ,intent(in):: I_nbas + logical,optional,intent(in):: GEN_GAUSSLATS + ! + ! Local Values + !-------------- + real(r8),allocatable:: Clats(:) + real(r8),allocatable:: Bcoef(:) + real(r8),allocatable:: Bcov (:,:) + real(r8):: Bnorm + integer :: ii,nn,n2,nb,ierr, astat + logical :: generate_lats + + character(len=*), parameter :: subname = 'init_ZonalProfile' + + generate_lats = .false. + + if (present(GEN_GAUSSLATS)) then + generate_lats = GEN_GAUSSLATS + end if + + ! Allocate space + !----------------- + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + + this%nlat = I_nlat + this%nbas = I_nbas + allocate(this%area (I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'this%area') + allocate(this%basis(I_nlat,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%basis') + allocate(this%map (I_nbas,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'this%map') + + allocate(Clats(I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(Bcoef(I_nbas/2+1), stat=astat) + call handle_allocate_error(astat, subname, 'Bcoef') + allocate(Bcov (I_nbas,I_nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Bcov') + + ! Optionally create the Latitude Gridpoints + ! and their associated area weights. Otherwise + ! they need to be supplied by the user. + !----------------------------------------------- + if(generate_lats) then + + ! Create a Gaussian grid from SP to NP + !-------------------------------------- + call sh_create_gaus_grid(I_nlat,Clats,IO_area,ierr) + if (ierr/=0) then + call endrun('init_ZonalProfile: Error creating Gaussian grid') + end if + + ! Convert generated colatitudes SP->NP to Lats and convert + ! to degrees and scale the area for global 2D integrals + !----------------------------------------------------------- + do nn=1,I_nlat + IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 + IO_area(nn) = IO_area(nn)*twoPI + end do + else + ! Convert Latitudes to SP->NP colatitudes in radians + !---------------------------------------------------- + do nn=1,I_nlat + Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 + end do + endif + + ! Copy the area weights for each nlat + ! gridpoint to the data structure + !--------------------------------------- + this%area(1:I_nlat) = IO_area(1:I_nlat) + + ! Add first basis for the mean values. + !------------------------------------------ + this%basis(:,1) = invSqrt4pi + Bnorm = 0._r8 + do ii=1,I_nlat + Bnorm = Bnorm + (this%basis(ii,1)*this%basis(ii,1)*this%area(ii)) + end do + this%basis(:,1) = this%basis(:,1)/sqrt(Bnorm) + + ! Loop over the remaining basis functions + !--------------------------------------- + do nn=2,I_nbas + nb = nn-1 + + ! Generate coefs for the basis + !------------------------------ + call sh_gen_basis_coefs(nb,0,Bcoef) + + ! Create an un-normalized basis for the + ! coefs at each nlat gridpoint + !--------------------------------------- + do ii=1,I_nlat + call sh_create_basis(nb,0,Clats(ii),Bcoef,this%basis(ii,nn)) + end do + + ! Numerically normalize the basis funnction + !-------------------------------------------------------------- + Bnorm = 0._r8 + do ii=1,I_nlat + Bnorm = Bnorm + (this%basis(ii,nn)*this%basis(ii,nn)*this%area(ii)) + end do + this%basis(:,nn) = this%basis(:,nn)/sqrt(Bnorm) + + end do ! nn=1,I_nbas + + ! Compute covariance matrix for basis functions + ! (Yes, they are theoretically orthonormal, but lets make sure) + !-------------------------------------------------------------- + do nn=1,I_nbas + do n2=1,I_nbas + Bcov(nn,n2) = 0._r8 + do ii=1,I_nlat + Bcov(nn,n2) = Bcov(nn,n2) + (this%basis(ii,nn)*this%basis(ii,n2)*this%area(ii)) + end do + end do + end do + + ! Invert to get the basis amplitude map + !-------------------------------------- + call Invert_Matrix(Bcov,I_nbas,this%map) + + ! End Routine + !------------ + deallocate(Clats) + deallocate(Bcoef) + deallocate(Bcov ) + + end subroutine init_ZonalProfile + !======================================================================= + + + !======================================================================= + subroutine final_ZonalProfile(this) + class(ZonalProfile_t) :: this + + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%basis)) deallocate(this%basis) + if(allocated(this%map )) deallocate(this%map) + + end subroutine final_ZonalProfile + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalProfile_1Damps(this,I_Zdata,O_Bamp) + ! + ! calc_ZonalProfile_1Damps: Given 1D data values for the nlat zonal + ! profiles gridpoints, compute the zonal + ! profile basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalProfile_t):: this + real(r8),intent(in ):: I_Zdata(:) + real(r8),intent(out):: O_Bamp (:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Gcov(:) + integer:: ii,nn,n2, astat + character(len=*), parameter :: subname = 'calc_ZonalProfile_1Damps' + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + allocate(Gcov(this%nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + do nn=1,this%nbas + Gcov(nn) = 0._r8 + do ii=1,this%nlat + Gcov(nn) = Gcov(nn) + (I_Zdata(ii)*this%basis(ii,nn)*this%area(ii)) + end do + end do + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do nn=1,this%nbas + O_Bamp(nn) = 0._r8 + do n2=1,this%nbas + O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2) + end do + end do + + deallocate(Gcov) + + end subroutine calc_ZonalProfile_1Damps + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp) + ! + ! calc_ZonalProfile_2Damps: Given 2D data values for the nlat,nlev zonal + ! profiles gridpoints, compute the zonal + ! profile basis amplitudes. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalProfile_t):: this + real(r8),intent(in ):: I_Zdata(:,:) + real(r8),intent(out):: O_Bamp (:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Gcov(:,:) + integer:: ii,nn,n2,ilev + + integer :: nlev, astat + character(len=*), parameter :: subname = 'calc_ZonalProfile_2Damps' + + nlev = size(I_Zdata,dim=2) + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + allocate(Gcov(this%nbas,nlev), stat=astat) + call handle_allocate_error(astat, subname, 'Gcov') + do ilev=1,nlev + do nn=1,this%nbas + Gcov(nn,ilev) = 0._r8 + do ii=1,this%nlat + Gcov(nn,ilev) = Gcov(nn,ilev) + (I_Zdata(ii,ilev)*this%basis(ii,nn)*this%area(ii)) + end do + end do + end do + + ! Multiply by map to get the amplitudes + !------------------------------------------- + do ilev=1,nlev + do nn=1,this%nbas + O_Bamp(nn,ilev) = 0._r8 + do n2=1,this%nbas + O_Bamp(nn,ilev) = O_Bamp(nn,ilev) + this%map(n2,nn)*Gcov(n2,ilev) + end do + end do + end do + deallocate(Gcov) + + end subroutine calc_ZonalProfile_2Damps + !======================================================================= + + + !======================================================================= + subroutine eval_ZonalProfile_1Dgrid(this,I_Bamp,O_Zdata) + ! + ! eval_ZonalProfile_1Dgrid: Given the zonal profile basis amplitudes, + ! compute 1D data values for the nlat gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalProfile_t):: this + real(r8),intent(in ):: I_Bamp (:) + real(r8),intent(out):: O_Zdata(:) + ! + ! Local Values + !-------------- + integer:: ii,nn + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + O_Zdata(1:this%nlat) = 0._r8 + do nn=1,this%nbas + do ii=1,this%nlat + O_Zdata(ii) = O_Zdata(ii) + (I_Bamp(nn)*this%basis(ii,nn)) + end do + end do + + end subroutine eval_ZonalProfile_1Dgrid + !======================================================================= + + + !======================================================================= + subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata) + ! + ! eval_ZonalProfile_2Dgrid: Given the zonal profile basis amplitudes, + ! compute 2D data values for the nlat,nlev gridpoints. + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalProfile_t):: this + real(r8),intent(in ):: I_Bamp (:,:) + real(r8),intent(out):: O_Zdata(:,:) + ! + ! Local Values + !-------------- + integer:: ii,nn,ilev + + integer :: nlev + + nlev = size(I_Bamp,dim=2) + + ! Construct grid values from basis amplitudes. + !-------------------------------------------------- + O_Zdata(1:this%nlat,1:nlev) = 0._r8 + do nn=1,this%nbas + do ilev=1,nlev + do ii=1,this%nlat + O_Zdata(ii,ilev) = O_Zdata(ii,ilev) + (I_Bamp(nn,ilev)*this%basis(ii,nn)) + end do + end do + end do + + end subroutine eval_ZonalProfile_2Dgrid + !======================================================================= + + + !======================================================================= + subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS,USE_LINEARWGTS) + ! + ! init_ZonalAverage: Initialize the ZonalAverage data structure for the + ! given nlat gridpoints. It is assumed that the domain + ! of these gridpoints of the profile span latitudes + ! from SP to NP. + ! + ! The optional GEN_GAUSSLATS flag allows for the + ! generation of Gaussian latitudes. The generated grid + ! over-writes the values of IO_lats/IO_area passed by + ! the user. + ! + ! The optional USE_LINEARWGTS flag allows for the + ! linear weighting of gridpoint values as they are + ! distributed to latitude bins, or to assign grid point + ! values to a single bin with equal weighting. + ! (default=TRUE for linear weighting) + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalAverage_t) :: this + real(r8) ,intent(inout):: IO_lats(:) + real(r8) ,intent(inout):: IO_area(:) + integer ,intent(in):: I_nlat + logical,optional,intent(in):: GEN_GAUSSLATS + logical,optional,intent(in):: USE_LINEARWGTS + ! + ! Local Values + !-------------- + real(r8),allocatable:: Clats (:) + real(r8),allocatable:: Glats (:,:) + real(r8),allocatable:: BinLat(:) + real(r8),allocatable:: Asum (:,:) + real(r8),allocatable:: Anorm (:) + real(r8):: area(pcols),rlat + integer :: nn,jj,ierr, astat + integer :: ncols,lchnk,cc,jlat + integer :: nlcols, count + logical :: generate_lats + character(len=*), parameter :: subname = 'init_ZonalAverage' + + generate_lats = .false. + + if(present(GEN_GAUSSLATS)) then + generate_lats = GEN_GAUSSLATS + endif + this%LINEARWGTS = .false. ! DEFAULT: Should be .true. ? + if(present(USE_LINEARWGTS)) then + this%LINEARWGTS = USE_LINEARWGTS + endif + + nlcols = get_nlcols_p() + + ! Allocate space + !----------------- + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%a_norm )) deallocate(this%a_norm) + if(allocated(this%area_g )) deallocate(this%area_g) + if(allocated(this%idx_map)) deallocate(this%idx_map) + if(allocated(this%idx_wgt)) deallocate(this%idx_wgt) + + this%nlat = I_nlat + allocate(this%area (I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'this%area') + allocate(this%a_norm (I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'this%a_norm') + allocate(this%area_g (pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%area_g') + allocate(this%idx_map(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%idx_map') + allocate(this%idx_wgt(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'this%idx_wgt') + + allocate(Clats (I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(BinLat(I_nlat+1), stat=astat) + call handle_allocate_error(astat, subname, 'BinLat') + allocate(Glats (pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Glats') + allocate(Asum (nlcols,I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Asum') + allocate(Anorm (I_nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Anorm') + + ! Optionally create the Latitude Gridpoints + ! and their associated area weights. Otherwise + ! they need to be supplied by the user. + !----------------------------------------------- + if(generate_lats) then + + ! Create a Gaussin grid from SP to NP + !-------------------------------------- + call sh_create_gaus_grid(this%nlat,Clats,IO_area,ierr) + if (ierr/=0) then + call endrun('init_ZonalAverage: Error creating Gaussian grid') + end if + + ! Convert generated colatitudes SP->NP to Lats and convert + ! to degrees and scale the area for global 2D integrals + !----------------------------------------------------------- + do nn=1,this%nlat + IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 + IO_area(nn) = IO_area(nn)*twoPI + end do + else + ! Convert Latitudes to SP->NP colatitudes in radians + !---------------------------------------------------- + do nn=1,this%nlat + Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 + end do + endif + + ! Copy the Lat grid area weights to the data structure + !----------------------------------------------------- + this%area(1:this%nlat) = IO_area(1:this%nlat) + + ! Save a copy of the area weights for each 2D gridpoint + ! and convert Latitudes to SP->NP colatitudes in radians + !------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + call get_wght_all_p(lchnk, ncols, area) + do cc = 1,ncols + rlat=get_rlat_p(lchnk,cc) + this%area_g(cc,lchnk) = area(cc) + Glats (cc,lchnk) = rlat + halfPI + end do + end do + + ! Set boundaries for Latitude bins + !----------------------------------- + BinLat(1) = 0._r8 + BinLat(this%nlat+1) = pi + do nn=2,this%nlat + BinLat(nn) = (Clats(nn-1)+Clats(nn))/2._r8 + end do + + if(this%LINEARWGTS) then + ! Each Gridpoint value is distributed with LINEAR weighting + ! to the two neighboring Latitude bins + !------------------------------------------------------ + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + if(Glats(cc,lchnk).lt.Clats(1)) then + this%idx_map(cc,lchnk) = 1 + this%idx_wgt(cc,lchnk) = (Glats(cc,lchnk)-Clats(2))/(Clats(2)-Clats(1)) + elseif(Glats(cc,lchnk).gt.Clats(I_nlat)) then + this%idx_map(cc,lchnk) = I_nlat-1 + this%idx_wgt(cc,lchnk) = (Glats(cc,lchnk)-Clats(I_nlat-1))/(Clats(I_nlat)-Clats(I_nlat-1)) + else + do jj=1,(I_nlat-1) + if((Glats(cc,lchnk).ge.Clats(jj )).and. & + (Glats(cc,lchnk).le.Clats(jj+1)) ) then + this%idx_map(cc,lchnk) = jj + this%idx_wgt(cc,lchnk) = (Clats(jj+1)-Glats(cc,lchnk))/(Clats(jj+1)-Clats(jj)) + exit + endif + end do + endif + end do + end do + + ! Compute normalization for weighted sums + !----------------------------------------- + Asum(:,:) = 0._r8 + Anorm(:) = 0._r8 + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + count=count+1 + if(this%idx_wgt(cc,lchnk) < -1._r8) then + ! South Pole point + !------------------ + Asum(count,jlat) = Asum(count,jlat) + this%area_g(cc,lchnk) + elseif(this%idx_wgt(cc,lchnk) > 1._r8) then + ! North Pole point + !------------------ + Asum(count,jlat+1) = Asum(count,jlat+1) + this%area_g(cc,lchnk) + else + ! Distribute area weight to bins + !-------------------------------- + Asum(count,jlat ) = Asum(count,jlat ) + this%area_g(cc,lchnk)* this%idx_wgt(cc,lchnk) + Asum(count,jlat+1) = Asum(count,jlat+1) + this%area_g(cc,lchnk)*(1._r8-this%idx_wgt(cc,lchnk)) + endif + end do + end do + call shr_reprosum_calc(Asum, Anorm, count, nlcols, I_nlat, gbl_count=ngcols_p, commid=mpicom) + + else + ! Each Gridpoint value is assigned to its Latitude bin + ! with equal weighting. + !------------------------------------------------------ + this%idx_wgt(:,:) = 1.0_r8 + + ! Loop over 2D gridpoints and determine its lat bin index + !--------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = -1 + if((Glats(cc,lchnk)<=BinLat(2)).and. & + (Glats(cc,lchnk)>=BinLat(1)) ) then + jlat = 1 + elseif((Glats(cc,lchnk)>=BinLat(this%nlat) ).and. & + (Glats(cc,lchnk)<=BinLat(this%nlat+1)) ) then + jlat = this%nlat + else + do jj=2,(this%nlat-1) + if((Glats(cc,lchnk)>BinLat(jj )).and. & + (Glats(cc,lchnk)<=BinLat(jj+1)) ) then + jlat = jj + exit + endif + end do + endif + if (jlat<1) then + call endrun('ZonalAverage init ERROR: jlat not in range') + endif + this%idx_map(cc,lchnk) = jlat + end do + end do + + ! Compute normalization for weighted sums + !-------------------------------------- + Asum(:,:) = 0._r8 + Anorm(:) = 0._r8 + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + count=count+1 + Asum(count,jlat) = this%area_g(cc,lchnk) + end do + end do + call shr_reprosum_calc(Asum, Anorm, count, nlcols, I_nlat, gbl_count=ngcols_p, commid=mpicom) + + endif + + ! Save normalization vallues for future + ! weghted-average calculations + !----------------------------------------- + this%a_norm = Anorm + if(.not.all(Anorm(:)>0._r8)) then + write(iulog,*) 'init_ZonalAverage -- ERROR in Anorm values: ' + do jlat = 1,I_nlat + if (.not.Anorm(jlat)>0._r8) then + write(iulog,*) ' Anorm(',jlat,'): ', Anorm(jlat) + endif + end do + call endrun('init_ZonalAverage -- ERROR in Anorm values') + endif + + ! End Routine + !------------ + deallocate(Clats) + deallocate(BinLat) + deallocate(Glats) + deallocate(Asum) + deallocate(Anorm) + + end subroutine init_ZonalAverage + !======================================================================= + + + !======================================================================= + subroutine final_ZonalAverage(this) + class(ZonalAverage_t) :: this + + if(allocated(this%area )) deallocate(this%area) + if(allocated(this%a_norm )) deallocate(this%a_norm) + if(allocated(this%area_g )) deallocate(this%area_g) + if(allocated(this%idx_map)) deallocate(this%idx_map) + + end subroutine final_ZonalAverage + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata) + ! + ! calc_ZonalAverage_2DbinAvg: Given 2D data values for ncol gridpoints, + ! compute the nlat area weighted binAvg profile + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalAverage_t):: this + real(r8),intent(in ):: I_Gdata(pcols,begchunk:endchunk) + real(r8),intent(out):: O_Zdata(:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Asum (:,:) + integer:: nn,ncols,lchnk,cc,jlat + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'calc_ZonalAverage_2DbinAvg' + + nlcols = get_nlcols_p() + + ! Initialize Zonal profile + !--------------------------- + allocate(Asum(nlcols,this%nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Asum') + + Asum (:,:) = 0._r8 + O_Zdata(1:this%nlat) = 0._r8 + + ! Compute area-weighted sums + !----------------------------- + if(this%LINEARWGTS) then + ! Each Gridpoint value is distributed with LINEAR weighting + ! to the two neighboring Latitude bins + !------------------------------------------------------ + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + count=count+1 + if(this%idx_wgt(cc,lchnk) < -1._r8) then + ! South Pole point + !------------------ + Asum(count,jlat) = Asum(count,jlat) + I_Gdata(cc,lchnk)*this%area_g(cc,lchnk) + elseif(this%idx_wgt(cc,lchnk) > 1._r8) then + ! North Pole point + !------------------ + Asum(count,jlat+1) = Asum(count,jlat+1) + I_Gdata(cc,lchnk)*this%area_g(cc,lchnk) + else + ! Distribute area weight to bins + !-------------------------------- + Asum(count,jlat ) = Asum(count,jlat ) & + + I_Gdata(cc,lchnk)*this%area_g(cc,lchnk)* this%idx_wgt(cc,lchnk) + Asum(count,jlat+1) = Asum(count,jlat+1) & + + I_Gdata(cc,lchnk)*this%area_g(cc,lchnk)*(1._r8-this%idx_wgt(cc,lchnk)) + endif + end do + end do + call shr_reprosum_calc(Asum,O_Zdata,count, nlcols, this%nlat,gbl_count=ngcols_p, commid=mpicom) + + else + ! Each Gridpoint value is assigned to its Latitude bin + ! with equal weighting. + !------------------------------------------------------ + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + count=count+1 + Asum(count,jlat) = I_Gdata(cc,lchnk)*this%area_g(cc,lchnk) + end do + end do + call shr_reprosum_calc(Asum,O_Zdata,count, nlcols, this%nlat,gbl_count=ngcols_p, commid=mpicom) + + endif + + ! Divide by area norm to get the averages + !----------------------------------------- + do nn=1,this%nlat + O_Zdata(nn) = O_Zdata(nn)/this%a_norm(nn) + end do + + deallocate(Asum) + + end subroutine calc_ZonalAverage_2DbinAvg + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) + ! + ! calc_ZonalAverage_3DbinAvg: Given 3D data values for ncol,nlev gridpoints, + ! compute the nlat,nlev area weighted binAvg profile + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalAverage_t):: this + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) + real(r8),intent(out):: O_Zdata(:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Gsum(:) + real(r8),allocatable:: Asum(:,:) + integer:: nn,ncols,lchnk,cc,jlat + integer:: Nsum,ilev,ns + + integer :: nlev + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'calc_ZonalAverage_3DbinAvg' + + nlev = size(I_Gdata,dim=2) + nlcols = get_nlcols_p() + + ! Initialize Zonal profile + !--------------------------- + Nsum = this%nlat*nlev + allocate(Gsum(Nsum), stat=astat) + call handle_allocate_error(astat, subname, 'Gsum') + allocate(Asum(nlcols,Nsum), stat=astat) + call handle_allocate_error(astat, subname, 'Asum') + + Asum (:,:) = 0._r8 + O_Zdata(1:this%nlat,1:nlev) = 0._r8 + + ! Compute area-weighted sums + !----------------------------- + if(this%LINEARWGTS) then + ! Each Gridpoint value is distributed with LINEAR weighting + ! to the two neighboring Latitude bins + !------------------------------------------------------ + do ilev = 1,nlev + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + ns = jlat + (ilev-1)*this%nlat + count=count+1 + if(this%idx_wgt(cc,lchnk) < -1._r8) then + ! South Pole point + !------------------ + Asum(count,ns) = Asum(count,ns) + I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk) + elseif(this%idx_wgt(cc,lchnk) > 1._r8) then + ! North Pole point + !------------------ + Asum(count,ns+1) = Asum(count,ns+1) + I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk) + else + ! Distribute area weight to bins + !-------------------------------- + Asum(count,ns ) = Asum(count,ns ) & + + I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk)* this%idx_wgt(cc,lchnk) + Asum(count,ns+1) = Asum(count,ns+1) & + + I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk)*(1._r8-this%idx_wgt(cc,lchnk)) + endif + end do + end do + end do + call shr_reprosum_calc(Asum,Gsum, count, nlcols, Nsum, gbl_count=ngcols_p, commid=mpicom) + + else + ! Each Gridpoint value is assigned to its Latitude bin + ! with equal weighting. + !------------------------------------------------------ + do ilev = 1,nlev + count = 0 + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + jlat = this%idx_map(cc,lchnk) + ns = jlat + (ilev-1)*this%nlat + count=count+1 + Asum(count,ns) = I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk) + end do + end do + end do + call shr_reprosum_calc(Asum,Gsum, count, nlcols, Nsum, gbl_count=ngcols_p, commid=mpicom) + + endif + + ! Divide by area norm to get the averages + !----------------------------------------- + do ilev = 1,nlev + do nn = 1,this%nlat + ns = nn + (ilev-1)*this%nlat + O_Zdata(nn,ilev) = Gsum(ns)/this%a_norm(nn) + end do + end do + + deallocate(Gsum) + deallocate(Asum) + + end subroutine calc_ZonalAverage_3DbinAvg + !======================================================================= + + + !======================================================================= + subroutine set_ZonalAverage_2Dgrid(this,I_Zdata,O_Gdata) + ! + ! set_ZonalAverage_2Dgrid: Set the 2D data values given + ! the 1D zonal Average profile values, + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalAverage_t) :: this + real(r8),intent(in ):: I_Zdata(:) + real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk) + ! + ! Local Values + !-------------- + integer:: ncols,lchnk,cc + + O_Gdata(:,:) = 0._r8 + + ! Construct 2D values from bin averages. + !---------------------------------------- + if(this%LINEARWGTS) then + ! Use linear weights for neighboring bins + !---------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + if(this%idx_wgt(cc,lchnk).lt.-1._r8) then + ! South Pole point + !------------------ + O_Gdata(cc,lchnk) = -this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk) ) & + +(1._r8+this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk)+1) + elseif(this%idx_wgt(cc,lchnk).gt.1._r8) then + ! North Pole point + !------------------ + O_Gdata(cc,lchnk) = this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk)+1) & + +(1._r8-this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk) ) + else + O_Gdata(cc,lchnk) = this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk) ) & + +(1._r8-this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk)+1) + endif + end do + end do + else + ! Set 2D gridpoint values in each bin equal to the given value + !--------------------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,lchnk) = I_Zdata(this%idx_map(cc,lchnk)) + end do + end do + endif + + end subroutine set_ZonalAverage_2Dgrid + !======================================================================= + + + !======================================================================= + subroutine set_ZonalAverage_3Dgrid(this,I_Zdata,O_Gdata) + ! + ! set_ZonalAverage_3Dgrid: Set the 3D data values given + ! the 2D zonal Average profile values, + !===================================================================== + ! + ! Passed Variables + !------------------ + class(ZonalAverage_t) :: this + real(r8),intent(in ):: I_Zdata(:,:) + real(r8),intent(out):: O_Gdata(:,:,begchunk:) + ! + ! Local Values + !-------------- + integer:: ncols,lchnk,cc + integer:: ll + + integer :: nlev + nlev = size(O_Gdata,dim=2) + + O_Gdata(:,:,:) = 0._r8 + + ! Construct 3D values from bin averages. + !---------------------------------------- + if(this%LINEARWGTS) then + ! Use linear weights for neighboring bins + !---------------------------------------- + do ll = 1,nlev + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + if(this%idx_wgt(cc,lchnk).lt.-1._r8) then + ! South Pole point + !------------------ + O_Gdata(cc,ll,lchnk) = -this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk) ,ll) & + +(1._r8+this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk)+1,ll) + elseif(this%idx_wgt(cc,lchnk).gt.1._r8) then + ! North Pole point + !------------------ + O_Gdata(cc,ll,lchnk) = this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk)+1,ll) & + +(1._r8-this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk) ,ll) + else + O_Gdata(cc,ll,lchnk) = this%idx_wgt(cc,lchnk) *I_Zdata(this%idx_map(cc,lchnk) ,ll) & + +(1._r8-this%idx_wgt(cc,lchnk))*I_Zdata(this%idx_map(cc,lchnk)+1,ll) + endif + end do + end do + end do + else + ! Set 2D gridpoint values in each bin equal to the given value + !--------------------------------------------------------------- + do ll = 1,nlev + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + O_Gdata(cc,ll,lchnk) = I_Zdata(this%idx_map(cc,lchnk),ll) + end do + end do + end do + endif + + end subroutine set_ZonalAverage_3Dgrid + !======================================================================= + + + !======================================================================= + subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) + ! + ! Invert_Matrix: Given the NbasxNbas matrix, calculate and return + ! the inverse of the matrix. + ! + ! Implemented with the LAPACK DGESV routine. + ! + !==================================================================== + ! + ! Passed Variables + !------------------ + real(r8), intent(inout) :: I_Mat(:,:) ! input matrix contains P*L*U + ! decomposition on output + integer, intent(in) :: Nbas + real(r8), intent(out) :: O_InvMat(:,:) + ! + ! Local Values + !------------- + integer, allocatable :: Indx(:) ! pivot indices + integer :: astat, ii + character(len=*), parameter :: subname = 'Invert_Matrix' + character(len=80) :: msg + + external DGESV + + ! Allocate work space + !--------------------- + allocate(Indx(Nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Indx') + + ! Initialize the inverse array with the identity matrix + !------------------------------------------------------- + O_InvMat(:,:) = 0._r8 + do ii=1,Nbas + O_InvMat(ii,ii) = 1._r8 + end do + + call DGESV(Nbas, Nbas, I_Mat, Nbas, Indx, O_InvMat, Nbas, astat) + + if (astat < 0) then + write(msg, '(a, i1, a)') 'argument # ', abs(astat), ' has an illegal value' + call endrun(subname//': DGESV error return: '//msg) + else if (astat > 0) then + call endrun(subname//': DGESV error return: matrix is singular') + end if + + deallocate(Indx) + + end subroutine Invert_Matrix + !======================================================================= + + !======================================================================= + ! legacy spherepack routines + !======================================================================= + subroutine sh_gen_basis_coefs(nn,mm,cp) + ! + ! spherepack alfk + ! + ! dimension of real cp(nn/2 + 1) + ! arguments + ! + ! purpose computes fourier coefficients in the trigonometric series + ! representation of the normalized associated + ! legendre function pbar(nn,mm,theta) for use by + ! sh_gen_basis_coefs in calculating pbar(nn,mm,theta). + ! + ! first define the normalized associated + ! legendre functions + ! + ! pbar(mm,nn,theta) = sqrt((2*nn+1)*factorial(nn-mm) + ! /(2*factorial(nn+mm)))*sin(theta)**mm/(2**nn* + ! factorial(nn)) times the (nn+mm)th derivative of + ! (x**2-1)**nn with respect to x=cos(theta) + ! + ! where theta is colatitude. + ! + ! then subroutine sh_gen_basis_coefs computes the coefficients + ! cp(k) in the following trigonometric + ! expansion of pbar(m,n,theta). + ! + ! 1) for n even and m even, pbar(mm,nn,theta) = + ! .5*cp(1) plus the sum from k=1 to k=nn/2 + ! of cp(k+1)*cos(2*k*th) + ! + ! 2) for nn even and mm odd, pbar(mm,nn,theta) = + ! the sum from k=1 to k=nn/2 of + ! cp(k)*sin(2*k*th) + ! + ! 3) for n odd and m even, pbar(mm,nn,theta) = + ! the sum from k=1 to k=(nn+1)/2 of + ! cp(k)*cos((2*k-1)*th) + ! + ! 4) for nn odd and mm odd, pbar(mm,nn,theta) = + ! the sum from k=1 to k=(nn+1)/2 of + ! cp(k)*sin((2*k-1)*th) + ! + ! arguments + ! + ! on input nn + ! nonnegative integer specifying the degree of + ! pbar(nn,mm,theta) + ! + ! mm + ! is the order of pbar(nn,mm,theta). mm can be + ! any integer however cp is computed such that + ! pbar(nn,mm,theta) = 0 if abs(m) is greater + ! than nn and pbar(nn,mm,theta) = (-1)**mm* + ! pbar(nn,-mm,theta) for negative mm. + ! + ! on output cp + ! array of length (nn/2)+1 + ! which contains the fourier coefficients in + ! the trigonometric series representation of + ! pbar(nn,mm,theta) + ! + ! special conditions none + ! + ! algorithm the highest order coefficient is determined in + ! closed form and the remainig coefficients are + ! determined as the solution of a backward + ! recurrence relation. + ! + !===================================================================== + ! + ! Passed Variables + !------------------ + integer ,intent(in ):: nn + integer ,intent(in ):: mm + real(r8),intent(out):: cp(nn/2+1) + ! + ! Local Values + !---------------- + real(r8):: fnum,fnmh + real(r8):: pm1 + real(r8):: t1,t2 + real(r8):: fden + real(r8):: cp2 + real(r8):: fnnp1 + real(r8):: fnmsq + real(r8):: fk + real(r8):: a1,b1,C1 + integer :: ma,nmms2,nex + integer :: ii,jj + + real(r8),parameter:: SC10=1024._r8 + real(r8),parameter:: SC20=SC10*SC10 + real(r8),parameter:: SC40=SC20*SC20 + + cp(1) = 0._r8 + ma = abs(mm) + if(ma>nn) return + + if((nn-1)<0) then + cp(1) = sqrt(2._r8) + return + elseif((nn-1)==0) then + if(ma/=0) then + cp(1) = sqrt(.75_r8) + if(mm==-1) cp(1) = -cp(1) + else + cp(1) = sqrt(1.5_r8) + endif + return + else + if(mod(nn+ma,2)/=0) then + nmms2 = (nn-ma-1)/2 + fnum = nn + ma + 2 + fnmh = nn - ma + 2 + pm1 = -1._r8 + else + nmms2 = (nn-ma)/2 + fnum = nn + ma + 1 + fnmh = nn - ma + 1 + pm1 = 1._r8 + endif + endif + + t1 = 1._r8/SC20 + nex = 20 + fden = 2._r8 + if(nmms2>=1) then + do ii = 1,nmms2 + t1 = fnum*t1/fden + if (t1>SC20) then + t1 = t1/SC40 + nex = nex + 40 + endif + fnum = fnum + 2._r8 + fden = fden + 2._r8 + end do + endif + + if(mod(ma/2,2)/=0) then + t1 = -t1/2._r8**(nn-1-nex) + else + t1 = t1/2._r8**(nn-1-nex) + endif + t2 = 1._r8 + if(ma/=0) then + do ii = 1,ma + t2 = fnmh*t2/ (fnmh+pm1) + fnmh = fnmh + 2._r8 + end do + endif + + cp2 = t1*sqrt((nn+.5_r8)*t2) + fnnp1 = nn*(nn+1) + fnmsq = fnnp1 - 2._r8*ma*ma + + if((mod(nn,2)==0).and.(mod(ma,2)==0)) then + jj = 1+(nn+1)/2 + else + jj = (nn+1)/2 + endif + + cp(jj) = cp2 + if(mm<0) then + if(mod(ma,2)/=0) cp(jj) = -cp(jj) + endif + if(jj<=1) return + + fk = nn + a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 + b1 = 2._r8* (fk*fk-fnmsq) + cp(jj-1) = b1*cp(jj)/a1 + + jj = jj - 1 + do while(jj>1) + fk = fk - 2._r8 + a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 + b1 = -2._r8*(fk*fk-fnmsq) + c1 = (fk+1._r8)*(fk+2._r8) - fnnp1 + cp(jj-1) = -(b1*cp(jj)+c1*cp(jj+1))/a1 + jj = jj - 1 + end do + + end subroutine sh_gen_basis_coefs + !======================================================================= + + !======================================================================= + subroutine sh_create_basis(nn,mm,theta,cp,pb) + ! + ! spherepack lfpt + ! + ! dimension of + ! arguments + ! cp((nn/2)+1) + ! + ! purpose routine sh_create_basis uses coefficients computed by + ! routine sh_gen_basis_coefs to compute the + ! normalized associated legendre function pbar(nn,mm,theta) + ! at colatitude theta. + ! + ! arguments + ! + ! on input nn + ! nonnegative integer specifying the degree of + ! pbar(nn,mm,theta) + ! mm + ! is the order of pbar(nn,mm,theta). mm can be + ! any integer however pbar(nn,mm,theta) = 0 + ! if abs(mm) is greater than nn and + ! pbar(nn,mm,theta) = (-1)**mm*pbar(nn,-mm,theta) + ! for negative mm. + ! + ! theta + ! colatitude in radians + ! + ! cp + ! array of length (nn/2)+1 + ! containing coefficients computed by routine + ! sh_gen_basis_coefs + ! + ! on output pb + ! variable containing pbar(n,m,theta) + ! + ! special conditions calls to routine sh_create_basis must be preceded by an + ! appropriate call to routine sh_gen_basis_coefs. + ! + ! algorithm the trigonometric series formula used by + ! routine sh_create_basis to calculate pbar(nn,mm,theta) at + ! colatitude theta depends on mm and nn as follows: + ! + ! 1) for nn even and mm even, the formula is + ! .5*cp(1) plus the sum from k=1 to k=n/2 + ! of cp(k)*cos(2*k*theta) + ! 2) for nn even and mm odd. the formula is + ! the sum from k=1 to k=nn/2 of + ! cp(k)*sin(2*k*theta) + ! 3) for nn odd and mm even, the formula is + ! the sum from k=1 to k=(nn+1)/2 of + ! cp(k)*cos((2*k-1)*theta) + ! 4) for nn odd and mm odd, the formula is + ! the sum from k=1 to k=(nn+1)/2 of + ! cp(k)*sin((2*k-1)*theta) + ! + !===================================================================== + integer, intent(in) :: nn,mm + real(r8), intent(in) :: theta + real(r8), intent(in) :: cp(:) + real(r8), intent(out) :: pb + + real(r8) :: cdt + real(r8) :: sdt + real(r8) :: ct + real(r8) :: st + real(r8) :: summ + real(r8) :: cth + + integer:: ma,nmod,mmod,kdo + integer:: kp1,kk + + pb = 0._r8 + ma = abs(mm) + if(ma>nn) return + + if(nn<=0) then + if(ma<=0) then + pb = sqrt(.5_r8) + return + endif + endif + + nmod = mod(nn,2) + mmod = mod(ma,2) + + if(nmod<=0) then + if(mmod<=0) then + kdo = nn/2 + 1 + cdt = cos(theta+theta) + sdt = sin(theta+theta) + ct = 1._r8 + st = 0._r8 + summ = .5_r8*cp(1) + do kp1 = 2,kdo + cth = cdt*ct - sdt*st + st = sdt*ct + cdt*st + ct = cth + summ = summ + cp(kp1)*ct + end do + pb = summ + return + endif + kdo = nn/2 + cdt = cos(theta+theta) + sdt = sin(theta+theta) + ct = 1._r8 + st = 0._r8 + summ = 0._r8 + do kk = 1,kdo + cth = cdt*ct - sdt*st + st = sdt*ct + cdt*st + ct = cth + summ = summ + cp(kk)*st + end do + pb = summ + return + endif + + kdo = (nn+1)/2 + if(mmod<=0) then + cdt = cos(theta+theta) + sdt = sin(theta+theta) + ct = cos(theta) + st = -sin(theta) + summ = 0._r8 + do kk = 1,kdo + cth = cdt*ct - sdt*st + st = sdt*ct + cdt*st + ct = cth + summ = summ + cp(kk)*ct + end do + pb = summ + return + endif + + cdt = cos(theta+theta) + sdt = sin(theta+theta) + ct = cos(theta) + st = -sin(theta) + summ = 0._r8 + do kk = 1,kdo + cth = cdt*ct - sdt*st + st = sdt*ct + cdt*st + ct = cth + summ = summ + cp(kk)*st + end do + pb = summ + + end subroutine sh_create_basis + !======================================================================= + + !======================================================================= + subroutine sh_create_gaus_grid(nlat,theta,wts,ierr) + ! + ! spherepack gaqd + ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + ! . . + ! . copyright (c) 2001 by ucar . + ! . . + ! . university corporation for atmospheric research . + ! . . + ! . all rights reserved . + ! . . + ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + ! + ! February 2002 + ! + ! gauss points and weights are computed using the fourier-newton + ! described in "on computing the points and weights for + ! gauss-legendre quadrature", paul n. swarztrauber, siam journal + ! on scientific computing (DOI 10.1137/S1064827500379690). + ! This routine is faster and more accurate than older program + ! with the same name. + ! + ! computes the nlat gaussian colatitudes and weights + ! in double precision. the colatitudes are in radians and lie in the + ! in the interval (0,pi). + ! + ! input parameters + ! + ! nlat the number of gaussian colatitudes in the interval (0,pi) + ! (between the two poles). nlat must be greater than zero. + ! + ! output parameters + ! + ! theta a double precision array with length nlat + ! containing the gaussian colatitudes in + ! increasing radians on the interval (0,pi). + ! + ! wts a double precision array with lenght nlat + ! containing the gaussian weights. + ! + ! ierror = 0 no errors + ! = 1 if nlat<=0 + ! + !=================================================================== + ! + ! Passed variables + !----------------- + integer ,intent(in ) :: nlat + real(r8),intent(out) :: theta(nlat) + real(r8),intent(out) :: wts(nlat) + integer ,intent(out) :: ierr + ! + ! Local Values + !------------- + real(r8):: sgnd + real(r8):: xx,dtheta,dthalf + real(r8):: cmax,zprev,zlast,zero,zhold,pb,dpb,dcor,summ,cz + integer :: mnlat,ns2,nhalf,nix,it,ii + + real(r8), parameter :: eps = epsilon(1._r8) + + ! check work space length + !------------------------ + if(nlat<=0) then + ierr = 1 + return + endif + ierr = 0 + + ! compute weights and points analytically when nlat=1,2 + !------------------------------------------------------- + if(nlat==1) then + theta(1) = acos(0._r8) + wts (1) = 2._r8 + return + elseif(nlat==2) then + xx = sqrt(1._r8/3._r8) + theta(1) = acos( xx) + theta(2) = acos(-xx) + wts (1) = 1._r8 + wts (2) = 1._r8 + return + endif + + ! Proceed for nlat > 2 + !---------------------- + mnlat = mod(nlat,2) + ns2 = nlat/2 + nhalf = (nlat+1)/2 + + call sh_fourier_coefs_dp(nlat,cz,theta(ns2+1),wts(ns2+1)) + + dtheta = halfPI/nhalf + dthalf = dtheta/2._r8 + cmax = .2_r8*dtheta + + ! estimate first point next to theta = pi/2 + !------------------------------------------- + if(mnlat/=0) then + zero = halfPI - dtheta + zprev = halfPI + nix = nhalf - 1 + else + zero = halfPI - dthalf + nix = nhalf + endif + + do while(nix/=0) + dcor = huge(1._r8) + it = 0 + do while (abs(dcor) > eps*abs(zero)) + it = it + 1 + ! newton iterations + !----------------------- + call sh_legp_dlegp_theta(nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb) + dcor = pb/dpb + if(dcor.ne.0._r8) then + sgnd = dcor/abs(dcor) + else + sgnd = 1._r8 + endif + dcor = sgnd*min(abs(dcor),cmax) + zero = zero - dcor + end do + + theta(nix) = zero + zhold = zero + + ! wts(nix) = (nlat+nlat+1)/(dpb*dpb) + ! yakimiw's formula permits using old pb and dpb + !-------------------------------------------------- + wts(nix) = (nlat+nlat+1)/ (dpb+pb*dcos(zlast)/dsin(zlast))**2 + nix = nix - 1 + if(nix==nhalf-1) zero = 3._r8*zero - pi + if(nix0) then + cth = cdt + sth = sdt + do kk = 1,kdo + pb = pb + cp(kk)*cth + dpb = dpb - dcp(kk)*sth + chh = cdt*cth - sdt*sth + sth = sdt*cth + cdt*sth + cth = chh + end do + endif + else + ! n odd + !----------- + kdo = (nn+1)/2 + pb = 0._r8 + dpb = 0._r8 + cth = dcos(theta) + sth = dsin(theta) + do kk = 1,kdo + pb = pb + cp(kk)*cth + dpb = dpb - dcp(kk)*sth + chh = cdt*cth - sdt*sth + sth = sdt*cth + cdt*sth + cth = chh + end do + endif + + end subroutine sh_legp_dlegp_theta + !======================================================================= + +end module ug_spectralmethods_mod diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 deleted file mode 100644 index 25e3f8564a..0000000000 --- a/src/utils/zonal_mean_mod.F90 +++ /dev/null @@ -1,2000 +0,0 @@ -module zonal_mean_mod -!====================================================================== -! -! Purpose: Compute and make use of Zonal Mean values on physgrid -! -! This module implements 3 data structures for the spectral analysis -! and synthesis of zonal mean values based on m=0 spherical harmonics. -! -! ZonalMean_t: For the analysis/synthesis of zonal mean values -! on a 2D grid of points distributed over the -! surface of a sphere. -! ZonalProfile_t: For the analysis/synthesis of zonal mean values -! on a meridional grid that spans the latitudes -! from SP to NP -! ZonalAverage_t: To calculate zonal mean values via a simple -! area weighted bin-averaging of 2D grid points -! assigned to each latitude band. -! -! NOTE: The weighting of the Zonal Profiles values is scaled such -! that ZonalMean_t amplitudes can be used to evaluate values -! on the ZonalProfile_t grid and vice-versa. -! -! The ZonalMean_t computes global integrals to compute basis -! amplitudes. For distributed environments the cost of these -! can be reduced using the The ZonalAverage_t data structures. -! -! USAGE: -! -! (1) Compute Zonal mean amplitudes and synthesize values on 2D/3D physgrid -! -! Usage: type(ZonalMean_t):: ZM -! ========================================= -! call ZM%init(nbas) -! ------------------ -! - Initialize the data structure with 'nbas' basis functions -! for the given physgrid latitudes and areas. -! -! Arguments: -! integer ,intent(in):: nbas -Number of m=0 spherical harmonics -! -! call ZM%calc_amps(Gdata,Bamp) -! ----------------------------- -! - For the initialized ZonalMean_t; Given Gdata() values on the physgrid, -! compute the zonal mean basis amplitudes Bamp(). -! -! Interface: 2D data on the physgrid -! real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk) -! real(r8),intent(out):: Bamp (nbas) -! -! Interface: 3D data on the physgrid -! real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk) -! real(r8),intent(out):: Bamp (nbas,pver) -! -! call ZM%eval_grid(Bamp,Gdata) -! ----------------------------- -! - For the initialized ZonalMean_t; Given Bamp() zonal mean basis -! amplitudes, compute the Gdata() values on the physgrid. -! -! Interface: 2D data on the physgrid -! real(r8),intent(in ):: Bamp (nbas) -! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk) -! -! Interface: 3D data on the physgrid -! real(r8),intent(in ):: Bamp (nbas,pver) -! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk) -! -! -! (2) Compute Zonal mean amplitudes and synthesize values on Zonal profile grid -! -! Usage: type(ZonalProfile_t):: ZP -! ========================================= -! call ZP%init(lats,area,nlat,nbas,GEN_GAUSSLATS=.true.) -! ------------------------------------------------------ -! - Initialize the data structure for the given number of -! latitudes. Either use the given Latitudes and weights, -! or OPTIONALLY create profile gridpoints and associated -! area weights from SP to NP. Then initialize 'nbas' basis -! functions for the profile gridpoints. -! If the user supplies the lats/area values, the area values must -! be correctly scaled such that the global area adds up to 4PI. -! Otherwise, the ampitudes between ZonalProfile_t and ZonalMean_t -! are not interchangable. -! -! Arguments: -! real(r8),intent(inout):: lats(:) - Latitudes of meridional grid. -! real(r8),intent(inout):: area(:) - Area of each meridional gridpoint. -! integer ,intent(in) :: nlat - Number of meridional gridpoints. -! integer ,intent(in) :: nbas - Number of m=0 spherical harmonics -! logical ,intent(in),optional:: GEN_GAUSLATS - Flag to generate -! lats/areas values. -! -! call ZP%calc_amps(Zdata,Bamp) -! ----------------------------- -! - Given Zdata() on the Zonal profile grid, compute the -! zonal basis amplitudes Bamp(). -! -! Interface: 1D data on (nlat) grid -! real(r8),intent(in ):: Zdata(nlat) - Meridional Profile data -! real(r8),intent(out):: Bamp (nbas) - Zonal Basis Amplitudes -! -! Interface: 2D data on (nlat,pver) grid -! real(r8),intent(in ):: Zdata(nlat,pver) - Meridional Profile data -! real(r8),intent(out):: Bamp (nbas,pver) - Zonal Basis Amplitudes -! -! call ZP%eval_grid(Bamp,Zdata) -! ----------------------------- -! - Given Bamp() zonal basis amplitudes, evaluate the Zdata() -! values on the Zonal profile grid. -! -! Interface: 1D data on (nlat) grid -! real(r8),intent(in ):: Bamp (nbas) - Zonal Basis Amplitudes -! real(r8),intent(out):: Zdata(nlat) - Meridional Profile data -! -! Interface: 2D data on (nlat,pver) grid -! real(r8),intent(in ):: Bamp (nbas,pver) - Zonal Basis Amplitudes -! real(r8),intent(out):: Zdata(nlat,pver) - Meridional Profile data -! -! (3) Compute Zonal mean averages (FASTER/LESS-ACCURATE) on Zonal profile grid -! (For the created zonal profile, just bin average area weighted -! 2D/3D physgrid grid values) -! -! Usage: type(ZonalAverage_t):: ZA -! ========================================= -! call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.) -! -------------------------------------------------- -! - Given the latitude/area for the nlat meridional gridpoints, initialize -! the ZonalAverage data structure for computing bin-averaging of physgrid -! values. It is assumed that the domain of these gridpoints of the -! profile span latitudes from SP to NP. -! The optional GEN_GAUSSLATS flag allows for the generation of Gaussian -! latitude gridpoints. The generated grid over-writes the given values -! lats and area passed by the user. -! -! Arguments: -! real(r8),intent(inout):: lats(nlat) - Latitudes of meridional grid. -! real(r8),intent(inout):: area(nlat) - Area of meridional gridpoints. -! integer ,intent(in):: nlat - Number of meridional gridpoints -! logical,intent(in),optional:: GEN_GAUSLATS - Flag to generate -! lats/areas values. -! -! call ZA%binAvg(Gdata,Zdata) -! --------------------------- -! - For the initialized ZonalAverage_t; Given Gdata() on the physgrid, -! compute bin averages and return Zdata() on the Zonal profile grid. -! -! Interface: 2D data on the physgrid -! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk) -! real(r8),intent(out):: Zdata(nlat) -! -! Interface: 3D data on the physgrid -! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk) -! real(r8),intent(out):: Zdata(nlat,pver) -! -!====================================================================== - - use shr_kind_mod, only: r8=>SHR_KIND_R8 - use phys_grid, only: get_ncols_p, get_rlat_p, get_wght_all_p, get_nlcols_p - use ppgrid, only: begchunk, endchunk, pcols - use shr_reprosum_mod,only: shr_reprosum_calc - use cam_abortutils, only: endrun, handle_allocate_error - use spmd_utils, only: mpicom - use physconst, only: pi - use phys_grid, only: ngcols_p => num_global_phys_cols - use cam_logfile, only: iulog - - implicit none - private - - public :: ZonalMean_t - public :: ZonalProfile_t - public :: ZonalAverage_t - - ! Type definitions - !------------------- - type ZonalMean_t - private - integer :: nbas - real(r8),allocatable:: area (:,:) - real(r8),allocatable:: basis(:,:,:) - real(r8),allocatable:: map (:,:) - contains - procedure,pass:: init => init_ZonalMean - generic,public:: calc_amps => calc_ZonalMean_2Damps, & - calc_ZonalMean_3Damps - generic,public:: eval_grid => eval_ZonalMean_2Dgrid, & - eval_ZonalMean_3Dgrid - procedure,private,pass:: calc_ZonalMean_2Damps - procedure,private,pass:: calc_ZonalMean_3Damps - procedure,private,pass:: eval_ZonalMean_2Dgrid - procedure,private,pass:: eval_ZonalMean_3Dgrid - procedure, pass :: final => final_ZonalMean - end type ZonalMean_t - - type ZonalProfile_t - private - integer :: nlat - integer :: nbas - real(r8),allocatable:: area (:) - real(r8),allocatable:: basis(:,:) - real(r8),allocatable:: map (:,:) - contains - procedure,pass:: init => init_ZonalProfile - generic,public:: calc_amps => calc_ZonalProfile_1Damps, & - calc_ZonalProfile_2Damps - generic,public:: eval_grid => eval_ZonalProfile_1Dgrid, & - eval_ZonalProfile_2Dgrid - procedure,private,pass:: calc_ZonalProfile_1Damps - procedure,private,pass:: calc_ZonalProfile_2Damps - procedure,private,pass:: eval_ZonalProfile_1Dgrid - procedure,private,pass:: eval_ZonalProfile_2Dgrid - procedure, pass :: final => final_ZonalProfile - end type ZonalProfile_t - - type ZonalAverage_t - private - integer :: nlat - real(r8),allocatable:: area (:) - real(r8),allocatable:: a_norm (:) - real(r8),allocatable:: area_g (:,:) - integer ,allocatable:: idx_map(:,:) - contains - procedure,pass:: init => init_ZonalAverage - generic,public:: binAvg => calc_ZonalAverage_2DbinAvg, & - calc_ZonalAverage_3DbinAvg - procedure,private,pass:: calc_ZonalAverage_2DbinAvg - procedure,private,pass:: calc_ZonalAverage_3DbinAvg - procedure, pass :: final => final_ZonalAverage - end type ZonalAverage_t - - real(r8), parameter :: halfPI = 0.5_r8*pi - real(r8), parameter :: twoPI = 2.0_r8*pi - real(r8), parameter :: fourPI = 4.0_r8*pi - real(r8), parameter :: qrtrPI = 0.25_r8*pi - real(r8), parameter :: invSqrt4pi = 1._r8/sqrt(fourPI) - -contains - !======================================================================= - subroutine init_ZonalMean(this,I_nbas) - ! - ! init_ZonalMean: Initialize the ZonalMean data structures for the - ! physics grid. It is assumed that the domain - ! of these gridpoints spans the surface of the sphere. - ! The representation of basis functions is - ! normalized w.r.t integration over the sphere. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalMean_t) :: this - integer ,intent(in):: I_nbas - ! - ! Local Values - !-------------- - real(r8),allocatable:: Clats(:,:) - real(r8),allocatable:: Bcoef(:) - real(r8),allocatable:: Csum (:,:) - real(r8),allocatable:: Cvec (:) - real(r8),allocatable:: Bsum (:,:) - real(r8),allocatable:: Bnorm(:) - real(r8),allocatable:: Bcov (:,:) - real(r8):: area(pcols),rlat - - integer :: nn,n2,nb,lchnk,ncols,cc - integer :: cnum,Cvec_len - - integer :: nlcols, count, astat - character(len=*), parameter :: subname = 'init_ZonalMean' - - if (I_nbas<1) then - call endrun('ZonalMean%init: ERROR I_nbas must be greater than 0') - end if - - ! Allocate space - !----------------- - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%basis)) deallocate(this%basis) - if(allocated(this%map )) deallocate(this%map) - - this%nbas = I_nbas - allocate(this%area (pcols,begchunk:endchunk), stat=astat) - call handle_allocate_error(astat, subname, 'this%area') - allocate(this%basis(pcols,begchunk:endchunk,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'this%basis') - allocate(this%map (I_nbas,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'this%map') - this%area (:,:) = 0._r8 - this%basis(:,:,:) = 0._r8 - this%map (:,:) = 0._r8 - - Cvec_len = 0 - do nn= 1,this%nbas - do n2=nn,this%nbas - Cvec_len = Cvec_len + 1 - end do - end do - - nlcols = get_nlcols_p() - - allocate(Clats(pcols,begchunk:endchunk), stat=astat) - call handle_allocate_error(astat, subname, 'Clats') - allocate(Bcoef(I_nbas/2+1), stat=astat) - call handle_allocate_error(astat, subname, 'Bcoef') - allocate(Csum (nlcols, Cvec_len), stat=astat) - call handle_allocate_error(astat, subname, 'Csum') - allocate(Cvec (Cvec_len), stat=astat) - call handle_allocate_error(astat, subname, 'Cvec') - allocate(Bsum (nlcols, I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Bsum') - allocate(Bnorm(I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Bnorm') - allocate(Bcov (I_nbas,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Bcov') - - Bsum(:,:) = 0._r8 - Csum(:,:) = 0._r8 - - ! Save a copy of the area weights for each ncol gridpoint - ! and convert Latitudes to SP->NP colatitudes in radians - !------------------------------------------------------- - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - call get_wght_all_p(lchnk, ncols, area) - do cc = 1,ncols - rlat=get_rlat_p(lchnk,cc) - this%area(cc,lchnk) = area(cc) - Clats (cc,lchnk) = rlat + halfPI - end do - end do - - ! Add first basis for the mean values. - !------------------------------------------ - this%basis(:,begchunk:endchunk,1) = invSqrt4pi - - ! Loop over the remaining basis functions - !--------------------------------------- - do nn=2,this%nbas - nb = nn-1 - - ! Generate coefs for the basis - !------------------------------ - call sh_gen_basis_coefs(nb,0,Bcoef) - - ! Create basis for the coefs at each ncol gridpoint - !--------------------------------------------------- - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - call sh_create_basis(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn)) - end do - end do - end do ! nn=2,this%nbas - - ! Numerically normalize the basis funnctions - !-------------------------------------------------------------- - do nn=1,this%nbas - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - count=count+1 - Bsum(count,nn) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) - end do - end do - end do ! nn=1,this%nbas - - call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) - - do nn=1,this%nbas - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(nn)) - end do - end do ! nn=1,this%nbas - - ! Compute covariance matrix for basis functions - ! (Yes, they are theoretically orthonormal, but lets make sure) - !--------------------------------------------------------------- - cnum = 0 - do nn= 1,this%nbas - do n2=nn,this%nbas - cnum = cnum + 1 - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - count=count+1 - Csum(count,cnum) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk) - end do - end do - - end do - end do - - call shr_reprosum_calc(Csum, Cvec, count, nlcols, Cvec_len, gbl_count=ngcols_p, commid=mpicom) - - cnum = 0 - do nn= 1,this%nbas - do n2=nn,this%nbas - cnum = cnum + 1 - Bcov(nn,n2) = Cvec(cnum) - Bcov(n2,nn) = Cvec(cnum) - end do - end do - - ! Invert to get the basis amplitude map - !-------------------------------------- - call Invert_Matrix(Bcov,this%nbas,this%map) - - ! End Routine - !------------ - deallocate(Clats) - deallocate(Bcoef) - deallocate(Csum ) - deallocate(Cvec ) - deallocate(Bsum ) - deallocate(Bnorm) - deallocate(Bcov ) - - end subroutine init_ZonalMean - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalMean_2Damps(this,I_Gdata,O_Bamp) - ! - ! calc_ZonalMean_2Damps: Given 2D data values for the ncol gridpoints, - ! compute the zonal mean basis amplitudes. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalMean_t) :: this - real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk) - real(r8),intent(out) :: O_Bamp(:) - ! - ! Local Values - !-------------- - real(r8),allocatable :: Csum(:,:) - real(r8),allocatable :: Gcov(:) - integer :: nn,n2,ncols,lchnk,cc - integer :: nlcols, count, astat - - character(len=*), parameter :: subname = 'calc_ZonalMean_2Damps' - - nlcols = get_nlcols_p() - - allocate(Gcov(this%nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Gcov') - allocate(Csum(nlcols, this%nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Csum') - Csum(:,:) = 0._r8 - - ! Compute Covariance with input data and basis functions - !-------------------------------------------------------- - do nn= 1,this%nbas - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - count=count+1 - Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) - end do - end do - end do - - call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) - - ! Multiply by map to get the amplitudes - !------------------------------------------- - do nn=1,this%nbas - O_Bamp(nn) = 0._r8 - do n2=1,this%nbas - O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2) - end do - end do - - ! End Routine - !------------ - deallocate(Csum) - deallocate(Gcov) - - end subroutine calc_ZonalMean_2Damps - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp) - ! - ! calc_ZonalMean_3Damps: Given 3D data values for the ncol,nlev gridpoints, - ! compute the zonal mean basis amplitudes. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalMean_t) :: this - real(r8),intent(in ):: I_Gdata(:,:,begchunk:) - real(r8),intent(out):: O_Bamp (:,:) - ! - ! Local Values - !-------------- - real(r8),allocatable:: Csum (:,:) - real(r8),allocatable:: Gcov (:) - integer:: nn,n2,ncols,lchnk,cc - integer:: Nsum,ns,ll - integer :: nlcols, count, astat - - integer :: nlev - character(len=*), parameter :: subname = 'calc_ZonalMean_3Damps' - - nlev = size(I_Gdata,dim=2) - - nlcols = get_nlcols_p() - allocate(Gcov(this%nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Gcov') - allocate(Csum(nlcols, this%nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Csum') - - Csum(:,:) = 0._r8 - O_Bamp(:,:) = 0._r8 - - ! Compute Covariance with input data and basis functions - !-------------------------------------------------------- - do ll= 1,nlev - - Csum(:,:) = 0._r8 - Gcov(:) = 0._r8 - - do nn= 1,this%nbas - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - count=count+1 - Csum(count,nn) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) - end do - end do - end do - - call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom) - - ! Multiply by map to get the amplitudes - !------------------------------------------- - do nn=1,this%nbas - O_Bamp(nn,ll) = 0._r8 - do n2=1,this%nbas - O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(n2) - end do - end do - - end do - - ! End Routine - !------------ - deallocate(Csum) - deallocate(Gcov) - - end subroutine calc_ZonalMean_3Damps - !======================================================================= - - - !======================================================================= - subroutine eval_ZonalMean_2Dgrid(this,I_Bamp,O_Gdata) - ! - ! eval_ZonalMean_2Dgrid: Given the zonal mean basis amplitudes, - ! compute 2D data values for the ncol gridpoints. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalMean_t) :: this - real(r8),intent(in ):: I_Bamp (:) - real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk) - ! - ! Local Values - !-------------- - integer:: nn,ncols,lchnk,cc - - O_Gdata(:,:) = 0._r8 - - ! Construct grid values from basis amplitudes. - !-------------------------------------------------- - - do nn=1,this%nbas - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn)) - end do - end do - end do - - end subroutine eval_ZonalMean_2Dgrid - !======================================================================= - - - !======================================================================= - subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata) - ! - ! eval_ZonalMean_3Dgrid: Given the zonal mean basis amplitudes, - ! compute 3D data values for the ncol,nlev gridpoints. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalMean_t) :: this - real(r8),intent(in ):: I_Bamp (:,:) - real(r8),intent(out):: O_Gdata(:,:,begchunk:) - ! - ! Local Values - !-------------- - integer:: nn,ncols,lchnk,cc - integer:: ll - - integer :: nlev - nlev = size(O_Gdata,dim=2) - - O_Gdata(:,:,:) = 0._r8 - - ! Construct grid values from basis amplitudes. - !-------------------------------------------------- - - do ll = 1,nlev - do nn=1,this%nbas - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn)) - end do - end do - end do - end do - - end subroutine eval_ZonalMean_3Dgrid - !======================================================================= - - !======================================================================= - subroutine final_ZonalMean(this) - class(ZonalMean_t) :: this - - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%basis)) deallocate(this%basis) - if(allocated(this%map )) deallocate(this%map) - - end subroutine final_ZonalMean - !======================================================================= - - !======================================================================= - subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) - ! - ! init_ZonalProfile: Initialize the ZonalProfile data structure for the - ! given nlat gridpoints. It is assumed that the domain - ! of these gridpoints of the profile span latitudes - ! from SP to NP. - ! The representation of basis functions functions is - ! normalized w.r.t integration over the sphere so that - ! when configured for tha same number of basis functions, - ! the calculated amplitudes are interchangable with - ! those for the ZonalMean_t class. - ! - ! The optional GEN_GAUSSLATS flag allows for the - ! generation of Gaussian latitudes. The generated grid - ! over-writes the values of IO_lats/IO_area passed by - ! the user. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalProfile_t) :: this - real(r8) ,intent(inout):: IO_lats(:) - real(r8) ,intent(inout):: IO_area(:) - integer ,intent(in):: I_nlat - integer ,intent(in):: I_nbas - logical,optional,intent(in):: GEN_GAUSSLATS - ! - ! Local Values - !-------------- - real(r8),allocatable:: Clats(:) - real(r8),allocatable:: Bcoef(:) - real(r8),allocatable:: Bcov (:,:) - real(r8):: Bnorm - integer :: ii,nn,n2,nb,ierr, astat - logical :: generate_lats - - character(len=*), parameter :: subname = 'init_ZonalProfile' - - generate_lats = .false. - - if (present(GEN_GAUSSLATS)) then - generate_lats = GEN_GAUSSLATS - end if - - ! Allocate space - !----------------- - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%basis)) deallocate(this%basis) - if(allocated(this%map )) deallocate(this%map) - - this%nlat = I_nlat - this%nbas = I_nbas - allocate(this%area (I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'this%area') - allocate(this%basis(I_nlat,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'this%basis') - allocate(this%map (I_nbas,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'this%map') - - allocate(Clats(I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'Clats') - allocate(Bcoef(I_nbas/2+1), stat=astat) - call handle_allocate_error(astat, subname, 'Bcoef') - allocate(Bcov (I_nbas,I_nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Bcov') - - ! Optionally create the Latitude Gridpoints - ! and their associated area weights. Otherwise - ! they need to be supplied by the user. - !----------------------------------------------- - if(generate_lats) then - - ! Create a Gaussian grid from SP to NP - !-------------------------------------- - call sh_create_gaus_grid(I_nlat,Clats,IO_area,ierr) - if (ierr/=0) then - call endrun('init_ZonalProfile: Error creating Gaussian grid') - end if - - ! Convert generated colatitudes SP->NP to Lats and convert - ! to degrees and scale the area for global 2D integrals - !----------------------------------------------------------- - do nn=1,I_nlat - IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 - IO_area(nn) = IO_area(nn)*twoPI - end do - else - ! Convert Latitudes to SP->NP colatitudes in radians - !---------------------------------------------------- - do nn=1,I_nlat - Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 - end do - endif - - ! Copy the area weights for each nlat - ! gridpoint to the data structure - !--------------------------------------- - this%area(1:I_nlat) = IO_area(1:I_nlat) - - ! Add first basis for the mean values. - !------------------------------------------ - this%basis(:,1) = invSqrt4pi - Bnorm = 0._r8 - do ii=1,I_nlat - Bnorm = Bnorm + (this%basis(ii,1)*this%basis(ii,1)*this%area(ii)) - end do - this%basis(:,1) = this%basis(:,1)/sqrt(Bnorm) - - ! Loop over the remaining basis functions - !--------------------------------------- - do nn=2,I_nbas - nb = nn-1 - - ! Generate coefs for the basis - !------------------------------ - call sh_gen_basis_coefs(nb,0,Bcoef) - - ! Create an un-normalized basis for the - ! coefs at each nlat gridpoint - !--------------------------------------- - do ii=1,I_nlat - call sh_create_basis(nb,0,Clats(ii),Bcoef,this%basis(ii,nn)) - end do - - ! Numerically normalize the basis funnction - !-------------------------------------------------------------- - Bnorm = 0._r8 - do ii=1,I_nlat - Bnorm = Bnorm + (this%basis(ii,nn)*this%basis(ii,nn)*this%area(ii)) - end do - this%basis(:,nn) = this%basis(:,nn)/sqrt(Bnorm) - - end do ! nn=1,I_nbas - - ! Compute covariance matrix for basis functions - ! (Yes, they are theoretically orthonormal, but lets make sure) - !-------------------------------------------------------------- - do nn=1,I_nbas - do n2=1,I_nbas - Bcov(nn,n2) = 0._r8 - do ii=1,I_nlat - Bcov(nn,n2) = Bcov(nn,n2) + (this%basis(ii,nn)*this%basis(ii,n2)*this%area(ii)) - end do - end do - end do - - ! Invert to get the basis amplitude map - !-------------------------------------- - call Invert_Matrix(Bcov,I_nbas,this%map) - - ! End Routine - !------------ - deallocate(Clats) - deallocate(Bcoef) - deallocate(Bcov ) - - end subroutine init_ZonalProfile - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalProfile_1Damps(this,I_Zdata,O_Bamp) - ! - ! calc_ZonalProfile_1Damps: Given 1D data values for the nlat zonal - ! profiles gridpoints, compute the zonal - ! profile basis amplitudes. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalProfile_t):: this - real(r8),intent(in ):: I_Zdata(:) - real(r8),intent(out):: O_Bamp (:) - ! - ! Local Values - !-------------- - real(r8),allocatable:: Gcov(:) - integer:: ii,nn,n2, astat - character(len=*), parameter :: subname = 'calc_ZonalProfile_1Damps' - - ! Compute Covariance with input data and basis functions - !-------------------------------------------------------- - allocate(Gcov(this%nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Gcov') - do nn=1,this%nbas - Gcov(nn) = 0._r8 - do ii=1,this%nlat - Gcov(nn) = Gcov(nn) + (I_Zdata(ii)*this%basis(ii,nn)*this%area(ii)) - end do - end do - - ! Multiply by map to get the amplitudes - !------------------------------------------- - do nn=1,this%nbas - O_Bamp(nn) = 0._r8 - do n2=1,this%nbas - O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2) - end do - end do - - deallocate(Gcov) - - end subroutine calc_ZonalProfile_1Damps - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp) - ! - ! calc_ZonalProfile_2Damps: Given 2D data values for the nlat,nlev zonal - ! profiles gridpoints, compute the zonal - ! profile basis amplitudes. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalProfile_t):: this - real(r8),intent(in ):: I_Zdata(:,:) - real(r8),intent(out):: O_Bamp (:,:) - ! - ! Local Values - !-------------- - real(r8),allocatable:: Gcov(:,:) - integer:: ii,nn,n2,ilev - - integer :: nlev, astat - character(len=*), parameter :: subname = 'calc_ZonalProfile_2Damps' - - nlev = size(I_Zdata,dim=2) - - ! Compute Covariance with input data and basis functions - !-------------------------------------------------------- - allocate(Gcov(this%nbas,nlev), stat=astat) - call handle_allocate_error(astat, subname, 'Gcov') - do ilev=1,nlev - do nn=1,this%nbas - Gcov(nn,ilev) = 0._r8 - do ii=1,this%nlat - Gcov(nn,ilev) = Gcov(nn,ilev) + (I_Zdata(ii,ilev)*this%basis(ii,nn)*this%area(ii)) - end do - end do - end do - - ! Multiply by map to get the amplitudes - !------------------------------------------- - do ilev=1,nlev - do nn=1,this%nbas - O_Bamp(nn,ilev) = 0._r8 - do n2=1,this%nbas - O_Bamp(nn,ilev) = O_Bamp(nn,ilev) + this%map(n2,nn)*Gcov(n2,ilev) - end do - end do - end do - deallocate(Gcov) - - end subroutine calc_ZonalProfile_2Damps - !======================================================================= - - - !======================================================================= - subroutine eval_ZonalProfile_1Dgrid(this,I_Bamp,O_Zdata) - ! - ! eval_ZonalProfile_1Dgrid: Given the zonal profile basis amplitudes, - ! compute 1D data values for the nlat gridpoints. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalProfile_t):: this - real(r8),intent(in ):: I_Bamp (:) - real(r8),intent(out):: O_Zdata(:) - ! - ! Local Values - !-------------- - integer:: ii,nn - - ! Construct grid values from basis amplitudes. - !-------------------------------------------------- - O_Zdata(1:this%nlat) = 0._r8 - do nn=1,this%nbas - do ii=1,this%nlat - O_Zdata(ii) = O_Zdata(ii) + (I_Bamp(nn)*this%basis(ii,nn)) - end do - end do - - end subroutine eval_ZonalProfile_1Dgrid - !======================================================================= - - - !======================================================================= - subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata) - ! - ! eval_ZonalProfile_2Dgrid: Given the zonal profile basis amplitudes, - ! compute 2D data values for the nlat,nlev gridpoints. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalProfile_t):: this - real(r8),intent(in ):: I_Bamp (:,:) - real(r8),intent(out):: O_Zdata(:,:) - ! - ! Local Values - !-------------- - integer:: ii,nn,ilev - - integer :: nlev - - nlev = size(I_Bamp,dim=2) - - ! Construct grid values from basis amplitudes. - !-------------------------------------------------- - O_Zdata(1:this%nlat,1:nlev) = 0._r8 - do nn=1,this%nbas - do ilev=1,nlev - do ii=1,this%nlat - O_Zdata(ii,ilev) = O_Zdata(ii,ilev) + (I_Bamp(nn,ilev)*this%basis(ii,nn)) - end do - end do - end do - - end subroutine eval_ZonalProfile_2Dgrid - !======================================================================= - - !======================================================================= - subroutine final_ZonalProfile(this) - class(ZonalProfile_t) :: this - - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%basis)) deallocate(this%basis) - if(allocated(this%map )) deallocate(this%map) - - end subroutine final_ZonalProfile - !======================================================================= - - !======================================================================= - subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) - ! - ! init_ZonalAverage: Initialize the ZonalAverage data structure for the - ! given nlat gridpoints. It is assumed that the domain - ! of these gridpoints of the profile span latitudes - ! from SP to NP. - ! - ! The optional GEN_GAUSSLATS flag allows for the - ! generation of Gaussian latitudes. The generated grid - ! over-writes the values of IO_lats/IO_area passed by - ! the user. - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalAverage_t) :: this - real(r8) ,intent(inout):: IO_lats(:) - real(r8) ,intent(inout):: IO_area(:) - integer ,intent(in):: I_nlat - logical,optional,intent(in):: GEN_GAUSSLATS - ! - ! Local Values - !-------------- - real(r8),allocatable:: Clats (:) - real(r8),allocatable:: Glats (:,:) - real(r8),allocatable:: BinLat(:) - real(r8),allocatable:: Asum (:,:) - real(r8),allocatable:: Anorm (:) - real(r8):: area(pcols),rlat - integer :: nn,jj,ierr, astat - integer :: ncols,lchnk,cc,jlat - integer :: nlcols, count - logical :: generate_lats - character(len=*), parameter :: subname = 'init_ZonalAverage' - - generate_lats = .false. - - if (present(GEN_GAUSSLATS)) then - generate_lats = GEN_GAUSSLATS - end if - - nlcols = get_nlcols_p() - - ! Allocate space - !----------------- - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%a_norm )) deallocate(this%a_norm) - if(allocated(this%area_g )) deallocate(this%area_g) - if(allocated(this%idx_map)) deallocate(this%idx_map) - - this%nlat = I_nlat - allocate(this%area (I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'this%area') - allocate(this%a_norm (I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'this%a_norm') - allocate(this%area_g (pcols,begchunk:endchunk), stat=astat) - call handle_allocate_error(astat, subname, 'this%area_g') - allocate(this%idx_map(pcols,begchunk:endchunk), stat=astat) - call handle_allocate_error(astat, subname, 'this%idx_map') - - allocate(Clats (I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'Clats') - allocate(BinLat(I_nlat+1), stat=astat) - call handle_allocate_error(astat, subname, 'BinLat') - allocate(Glats (pcols,begchunk:endchunk), stat=astat) - call handle_allocate_error(astat, subname, 'Glats') - allocate(Asum (nlcols,I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'Asum') - allocate(Anorm (I_nlat), stat=astat) - call handle_allocate_error(astat, subname, 'Anorm') - - ! Optionally create the Latitude Gridpoints - ! and their associated area weights. Otherwise - ! they need to be supplied by the user. - !----------------------------------------------- - if(generate_lats) then - - ! Create a Gaussin grid from SP to NP - !-------------------------------------- - call sh_create_gaus_grid(this%nlat,Clats,IO_area,ierr) - if (ierr/=0) then - call endrun('init_ZonalAverage: Error creating Gaussian grid') - end if - - ! Convert generated colatitudes SP->NP to Lats and convert - ! to degrees and scale the area for global 2D integrals - !----------------------------------------------------------- - do nn=1,this%nlat - IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 - IO_area(nn) = IO_area(nn)*twoPI - end do - else - ! Convert Latitudes to SP->NP colatitudes in radians - !---------------------------------------------------- - do nn=1,this%nlat - Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 - end do - endif - - ! Copy the Lat grid area weights to the data structure - !----------------------------------------------------- - this%area(1:this%nlat) = IO_area(1:this%nlat) - - ! Save a copy of the area weights for each 2D gridpoint - ! and convert Latitudes to SP->NP colatitudes in radians - !------------------------------------------------------- - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - call get_wght_all_p(lchnk, ncols, area) - do cc = 1,ncols - rlat=get_rlat_p(lchnk,cc) - this%area_g(cc,lchnk) = area(cc) - Glats (cc,lchnk) = rlat + halfPI - end do - end do - - ! Set boundaries for Latitude bins - !----------------------------------- - BinLat(1) = 0._r8 - BinLat(this%nlat+1) = pi - do nn=2,this%nlat - BinLat(nn) = (Clats(nn-1)+Clats(nn))/2._r8 - end do - - ! Loop over 2D gridpoints and determine its lat bin index - !--------------------------------------------------------- - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - jlat = -1 - if((Glats(cc,lchnk)<=BinLat(2)).and. & - (Glats(cc,lchnk)>=BinLat(1)) ) then - jlat = 1 - elseif((Glats(cc,lchnk)>=BinLat(this%nlat) ).and. & - (Glats(cc,lchnk)<=BinLat(this%nlat+1)) ) then - jlat = this%nlat - else - do jj=2,(this%nlat-1) - if((Glats(cc,lchnk)>BinLat(jj )).and. & - (Glats(cc,lchnk)<=BinLat(jj+1)) ) then - jlat = jj - exit - endif - end do - endif - if (jlat<1) then - call endrun('ZonalAverage init ERROR: jlat not in range') - endif - this%idx_map(cc,lchnk) = jlat - end do - end do - - ! Initialize 2D Area sums for each bin - !-------------------------------------- - Asum(:,:) = 0._r8 - Anorm(:) = 0._r8 - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - jlat = this%idx_map(cc,lchnk) - count=count+1 - Asum(count,jlat) = this%area_g(cc,lchnk) - end do - end do - - call shr_reprosum_calc(Asum, Anorm, count, nlcols, I_nlat, gbl_count=ngcols_p, commid=mpicom) - - this%a_norm = Anorm - - if (.not.all(Anorm(:)>0._r8)) then - write(iulog,*) 'init_ZonalAverage -- ERROR in Anorm values: ' - do jlat = 1,I_nlat - if (.not.Anorm(jlat)>0._r8) then - write(iulog,*) ' Anorm(',jlat,'): ', Anorm(jlat) - endif - end do - call endrun('init_ZonalAverage -- ERROR in Anorm values') - end if - - ! End Routine - !------------ - deallocate(Clats) - deallocate(BinLat) - deallocate(Glats) - deallocate(Asum) - deallocate(Anorm) - - end subroutine init_ZonalAverage - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata) - ! - ! calc_ZonalAverage_2DbinAvg: Given 2D data values for ncol gridpoints, - ! compute the nlat area weighted binAvg profile - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalAverage_t):: this - real(r8),intent(in ):: I_Gdata(pcols,begchunk:endchunk) - real(r8),intent(out):: O_Zdata(:) - ! - ! Local Values - !-------------- - real(r8),allocatable:: Asum (:,:) - integer:: nn,ncols,lchnk,cc,jlat - integer :: nlcols, count, astat - character(len=*), parameter :: subname = 'calc_ZonalAverage_2DbinAvg' - - nlcols = get_nlcols_p() - - - ! Initialize Zonal profile - !--------------------------- - allocate(Asum(nlcols,this%nlat), stat=astat) - call handle_allocate_error(astat, subname, 'Asum') - Asum(:,:) = 0._r8 - - O_Zdata(1:this%nlat) = 0._r8 - - ! Compute area-weighted sums - !----------------------------- - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - jlat = this%idx_map(cc,lchnk) - count=count+1 - Asum(count,jlat) = I_Gdata(cc,lchnk)*this%area_g(cc,lchnk) - end do - end do - - call shr_reprosum_calc(Asum,O_Zdata,count, nlcols, this%nlat,gbl_count=ngcols_p, commid=mpicom) - - ! Divide by area norm to get the averages - !----------------------------------------- - do nn=1,this%nlat - O_Zdata(nn) = O_Zdata(nn)/this%a_norm(nn) - end do - - deallocate(Asum) - - end subroutine calc_ZonalAverage_2DbinAvg - !======================================================================= - - - !======================================================================= - subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) - ! - ! calc_ZonalAverage_3DbinAvg: Given 3D data values for ncol,nlev gridpoints, - ! compute the nlat,nlev area weighted binAvg profile - !===================================================================== - ! - ! Passed Variables - !------------------ - class(ZonalAverage_t):: this - real(r8),intent(in ):: I_Gdata(:,:,begchunk:) - real(r8),intent(out):: O_Zdata(:,:) - ! - ! Local Values - !-------------- - real(r8),allocatable:: Gsum(:) - real(r8),allocatable:: Asum(:,:) - integer:: nn,ncols,lchnk,cc,jlat - integer:: Nsum,ilev,ns - - integer :: nlev - integer :: nlcols, count, astat - character(len=*), parameter :: subname = 'calc_ZonalAverage_3DbinAvg' - - nlev = size(I_Gdata,dim=2) - nlcols = get_nlcols_p() - - ! Initialize Zonal profile - !--------------------------- - Nsum = this%nlat*nlev - allocate(Gsum(Nsum), stat=astat) - call handle_allocate_error(astat, subname, 'Gsum') - allocate(Asum(nlcols,Nsum), stat=astat) - call handle_allocate_error(astat, subname, 'Asum') - Asum(:,:) = 0._r8 - - O_Zdata(1:this%nlat,1:nlev) = 0._r8 - - ! Compute area-weighted sums - !----------------------------- - do ilev = 1,nlev - count = 0 - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - jlat = this%idx_map(cc,lchnk) - ns = jlat + (ilev-1)*this%nlat - count=count+1 - Asum(count,ns) = I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk) - end do - end do - end do - - call shr_reprosum_calc(Asum,Gsum, count, nlcols, Nsum, gbl_count=ngcols_p, commid=mpicom) - - ! Divide by area norm to get the averages - !----------------------------------------- - do ilev = 1,nlev - do nn = 1,this%nlat - ns = nn + (ilev-1)*this%nlat - O_Zdata(nn,ilev) = Gsum(ns)/this%a_norm(nn) - end do - end do - - deallocate(Gsum) - deallocate(Asum) - - end subroutine calc_ZonalAverage_3DbinAvg - !======================================================================= - - !======================================================================= - subroutine final_ZonalAverage(this) - class(ZonalAverage_t) :: this - - if(allocated(this%area )) deallocate(this%area) - if(allocated(this%a_norm )) deallocate(this%a_norm) - if(allocated(this%area_g )) deallocate(this%area_g) - if(allocated(this%idx_map)) deallocate(this%idx_map) - - end subroutine final_ZonalAverage - !======================================================================= - - - !======================================================================= - subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) - ! - ! Invert_Matrix: Given the NbasxNbas matrix, calculate and return - ! the inverse of the matrix. - ! - ! Implemented with the LAPACK DGESV routine. - ! - !==================================================================== - ! - ! Passed Variables - !------------------ - real(r8), intent(inout) :: I_Mat(:,:) ! input matrix contains P*L*U - ! decomposition on output - integer, intent(in) :: Nbas - real(r8), intent(out) :: O_InvMat(:,:) - ! - ! Local Values - !------------- - integer, allocatable :: Indx(:) ! pivot indices - integer :: astat, ii - character(len=*), parameter :: subname = 'Invert_Matrix' - character(len=80) :: msg - - external DGESV - - ! Allocate work space - !--------------------- - allocate(Indx(Nbas), stat=astat) - call handle_allocate_error(astat, subname, 'Indx') - - ! Initialize the inverse array with the identity matrix - !------------------------------------------------------- - O_InvMat(:,:) = 0._r8 - do ii=1,Nbas - O_InvMat(ii,ii) = 1._r8 - end do - - call DGESV(Nbas, Nbas, I_Mat, Nbas, Indx, O_InvMat, Nbas, astat) - - if (astat < 0) then - write(msg, '(a, i1, a)') 'argument # ', abs(astat), ' has an illegal value' - call endrun(subname//': DGESV error return: '//msg) - else if (astat > 0) then - call endrun(subname//': DGESV error return: matrix is singular') - end if - - deallocate(Indx) - - end subroutine Invert_Matrix - !======================================================================= - - !======================================================================= - ! legacy spherepack routines - !======================================================================= - subroutine sh_gen_basis_coefs(nn,mm,cp) - ! - ! spherepack alfk - ! - ! dimension of real cp(nn/2 + 1) - ! arguments - ! - ! purpose computes fourier coefficients in the trigonometric series - ! representation of the normalized associated - ! legendre function pbar(nn,mm,theta) for use by - ! sh_gen_basis_coefs in calculating pbar(nn,mm,theta). - ! - ! first define the normalized associated - ! legendre functions - ! - ! pbar(mm,nn,theta) = sqrt((2*nn+1)*factorial(nn-mm) - ! /(2*factorial(nn+mm)))*sin(theta)**mm/(2**nn* - ! factorial(nn)) times the (nn+mm)th derivative of - ! (x**2-1)**nn with respect to x=cos(theta) - ! - ! where theta is colatitude. - ! - ! then subroutine sh_gen_basis_coefs computes the coefficients - ! cp(k) in the following trigonometric - ! expansion of pbar(m,n,theta). - ! - ! 1) for n even and m even, pbar(mm,nn,theta) = - ! .5*cp(1) plus the sum from k=1 to k=nn/2 - ! of cp(k+1)*cos(2*k*th) - ! - ! 2) for nn even and mm odd, pbar(mm,nn,theta) = - ! the sum from k=1 to k=nn/2 of - ! cp(k)*sin(2*k*th) - ! - ! 3) for n odd and m even, pbar(mm,nn,theta) = - ! the sum from k=1 to k=(nn+1)/2 of - ! cp(k)*cos((2*k-1)*th) - ! - ! 4) for nn odd and mm odd, pbar(mm,nn,theta) = - ! the sum from k=1 to k=(nn+1)/2 of - ! cp(k)*sin((2*k-1)*th) - ! - ! arguments - ! - ! on input nn - ! nonnegative integer specifying the degree of - ! pbar(nn,mm,theta) - ! - ! mm - ! is the order of pbar(nn,mm,theta). mm can be - ! any integer however cp is computed such that - ! pbar(nn,mm,theta) = 0 if abs(m) is greater - ! than nn and pbar(nn,mm,theta) = (-1)**mm* - ! pbar(nn,-mm,theta) for negative mm. - ! - ! on output cp - ! array of length (nn/2)+1 - ! which contains the fourier coefficients in - ! the trigonometric series representation of - ! pbar(nn,mm,theta) - ! - ! special conditions none - ! - ! algorithm the highest order coefficient is determined in - ! closed form and the remainig coefficients are - ! determined as the solution of a backward - ! recurrence relation. - ! - !===================================================================== - ! - ! Passed Variables - !------------------ - integer ,intent(in ):: nn - integer ,intent(in ):: mm - real(r8),intent(out):: cp(nn/2+1) - ! - ! Local Values - !---------------- - real(r8):: fnum,fnmh - real(r8):: pm1 - real(r8):: t1,t2 - real(r8):: fden - real(r8):: cp2 - real(r8):: fnnp1 - real(r8):: fnmsq - real(r8):: fk - real(r8):: a1,b1,C1 - integer :: ma,nmms2,nex - integer :: ii,jj - - real(r8),parameter:: SC10=1024._r8 - real(r8),parameter:: SC20=SC10*SC10 - real(r8),parameter:: SC40=SC20*SC20 - - cp(1) = 0._r8 - ma = abs(mm) - if(ma>nn) return - - if((nn-1)<0) then - cp(1) = sqrt(2._r8) - return - elseif((nn-1)==0) then - if(ma/=0) then - cp(1) = sqrt(.75_r8) - if(mm==-1) cp(1) = -cp(1) - else - cp(1) = sqrt(1.5_r8) - endif - return - else - if(mod(nn+ma,2)/=0) then - nmms2 = (nn-ma-1)/2 - fnum = nn + ma + 2 - fnmh = nn - ma + 2 - pm1 = -1._r8 - else - nmms2 = (nn-ma)/2 - fnum = nn + ma + 1 - fnmh = nn - ma + 1 - pm1 = 1._r8 - endif - endif - - t1 = 1._r8/SC20 - nex = 20 - fden = 2._r8 - if(nmms2>=1) then - do ii = 1,nmms2 - t1 = fnum*t1/fden - if (t1>SC20) then - t1 = t1/SC40 - nex = nex + 40 - endif - fnum = fnum + 2._r8 - fden = fden + 2._r8 - end do - endif - - if(mod(ma/2,2)/=0) then - t1 = -t1/2._r8**(nn-1-nex) - else - t1 = t1/2._r8**(nn-1-nex) - endif - t2 = 1._r8 - if(ma/=0) then - do ii = 1,ma - t2 = fnmh*t2/ (fnmh+pm1) - fnmh = fnmh + 2._r8 - end do - endif - - cp2 = t1*sqrt((nn+.5_r8)*t2) - fnnp1 = nn*(nn+1) - fnmsq = fnnp1 - 2._r8*ma*ma - - if((mod(nn,2)==0).and.(mod(ma,2)==0)) then - jj = 1+(nn+1)/2 - else - jj = (nn+1)/2 - endif - - cp(jj) = cp2 - if(mm<0) then - if(mod(ma,2)/=0) cp(jj) = -cp(jj) - endif - if(jj<=1) return - - fk = nn - a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 - b1 = 2._r8* (fk*fk-fnmsq) - cp(jj-1) = b1*cp(jj)/a1 - - jj = jj - 1 - do while(jj>1) - fk = fk - 2._r8 - a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 - b1 = -2._r8*(fk*fk-fnmsq) - c1 = (fk+1._r8)*(fk+2._r8) - fnnp1 - cp(jj-1) = -(b1*cp(jj)+c1*cp(jj+1))/a1 - jj = jj - 1 - end do - - end subroutine sh_gen_basis_coefs - !======================================================================= - - !======================================================================= - subroutine sh_create_basis(nn,mm,theta,cp,pb) - ! - ! spherepack lfpt - ! - ! dimension of - ! arguments - ! cp((nn/2)+1) - ! - ! purpose routine sh_create_basis uses coefficients computed by - ! routine sh_gen_basis_coefs to compute the - ! normalized associated legendre function pbar(nn,mm,theta) - ! at colatitude theta. - ! - ! arguments - ! - ! on input nn - ! nonnegative integer specifying the degree of - ! pbar(nn,mm,theta) - ! mm - ! is the order of pbar(nn,mm,theta). mm can be - ! any integer however pbar(nn,mm,theta) = 0 - ! if abs(mm) is greater than nn and - ! pbar(nn,mm,theta) = (-1)**mm*pbar(nn,-mm,theta) - ! for negative mm. - ! - ! theta - ! colatitude in radians - ! - ! cp - ! array of length (nn/2)+1 - ! containing coefficients computed by routine - ! sh_gen_basis_coefs - ! - ! on output pb - ! variable containing pbar(n,m,theta) - ! - ! special conditions calls to routine sh_create_basis must be preceded by an - ! appropriate call to routine sh_gen_basis_coefs. - ! - ! algorithm the trigonometric series formula used by - ! routine sh_create_basis to calculate pbar(nn,mm,theta) at - ! colatitude theta depends on mm and nn as follows: - ! - ! 1) for nn even and mm even, the formula is - ! .5*cp(1) plus the sum from k=1 to k=n/2 - ! of cp(k)*cos(2*k*theta) - ! 2) for nn even and mm odd. the formula is - ! the sum from k=1 to k=nn/2 of - ! cp(k)*sin(2*k*theta) - ! 3) for nn odd and mm even, the formula is - ! the sum from k=1 to k=(nn+1)/2 of - ! cp(k)*cos((2*k-1)*theta) - ! 4) for nn odd and mm odd, the formula is - ! the sum from k=1 to k=(nn+1)/2 of - ! cp(k)*sin((2*k-1)*theta) - ! - !===================================================================== - integer, intent(in) :: nn,mm - real(r8), intent(in) :: theta - real(r8), intent(in) :: cp(:) - real(r8), intent(out) :: pb - - real(r8) :: cdt - real(r8) :: sdt - real(r8) :: ct - real(r8) :: st - real(r8) :: summ - real(r8) :: cth - - integer:: ma,nmod,mmod,kdo - integer:: kp1,kk - - pb = 0._r8 - ma = abs(mm) - if(ma>nn) return - - if(nn<=0) then - if(ma<=0) then - pb = sqrt(.5_r8) - return - endif - endif - - nmod = mod(nn,2) - mmod = mod(ma,2) - - if(nmod<=0) then - if(mmod<=0) then - kdo = nn/2 + 1 - cdt = cos(theta+theta) - sdt = sin(theta+theta) - ct = 1._r8 - st = 0._r8 - summ = .5_r8*cp(1) - do kp1 = 2,kdo - cth = cdt*ct - sdt*st - st = sdt*ct + cdt*st - ct = cth - summ = summ + cp(kp1)*ct - end do - pb = summ - return - endif - kdo = nn/2 - cdt = cos(theta+theta) - sdt = sin(theta+theta) - ct = 1._r8 - st = 0._r8 - summ = 0._r8 - do kk = 1,kdo - cth = cdt*ct - sdt*st - st = sdt*ct + cdt*st - ct = cth - summ = summ + cp(kk)*st - end do - pb = summ - return - endif - - kdo = (nn+1)/2 - if(mmod<=0) then - cdt = cos(theta+theta) - sdt = sin(theta+theta) - ct = cos(theta) - st = -sin(theta) - summ = 0._r8 - do kk = 1,kdo - cth = cdt*ct - sdt*st - st = sdt*ct + cdt*st - ct = cth - summ = summ + cp(kk)*ct - end do - pb = summ - return - endif - - cdt = cos(theta+theta) - sdt = sin(theta+theta) - ct = cos(theta) - st = -sin(theta) - summ = 0._r8 - do kk = 1,kdo - cth = cdt*ct - sdt*st - st = sdt*ct + cdt*st - ct = cth - summ = summ + cp(kk)*st - end do - pb = summ - - end subroutine sh_create_basis - !======================================================================= - - !======================================================================= - subroutine sh_create_gaus_grid(nlat,theta,wts,ierr) - ! - ! spherepack gaqd - ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - ! . . - ! . copyright (c) 2001 by ucar . - ! . . - ! . university corporation for atmospheric research . - ! . . - ! . all rights reserved . - ! . . - ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - ! - ! February 2002 - ! - ! gauss points and weights are computed using the fourier-newton - ! described in "on computing the points and weights for - ! gauss-legendre quadrature", paul n. swarztrauber, siam journal - ! on scientific computing (DOI 10.1137/S1064827500379690). - ! This routine is faster and more accurate than older program - ! with the same name. - ! - ! computes the nlat gaussian colatitudes and weights - ! in double precision. the colatitudes are in radians and lie in the - ! in the interval (0,pi). - ! - ! input parameters - ! - ! nlat the number of gaussian colatitudes in the interval (0,pi) - ! (between the two poles). nlat must be greater than zero. - ! - ! output parameters - ! - ! theta a double precision array with length nlat - ! containing the gaussian colatitudes in - ! increasing radians on the interval (0,pi). - ! - ! wts a double precision array with lenght nlat - ! containing the gaussian weights. - ! - ! ierror = 0 no errors - ! = 1 if nlat<=0 - ! - !=================================================================== - ! - ! Passed variables - !----------------- - integer ,intent(in ) :: nlat - real(r8),intent(out) :: theta(nlat) - real(r8),intent(out) :: wts(nlat) - integer ,intent(out) :: ierr - ! - ! Local Values - !------------- - real(r8):: sgnd - real(r8):: xx,dtheta,dthalf - real(r8):: cmax,zprev,zlast,zero,zhold,pb,dpb,dcor,summ,cz - integer :: mnlat,ns2,nhalf,nix,it,ii - - real(r8), parameter :: eps = epsilon(1._r8) - - ! check work space length - !------------------------ - if(nlat<=0) then - ierr = 1 - return - endif - ierr = 0 - - ! compute weights and points analytically when nlat=1,2 - !------------------------------------------------------- - if(nlat==1) then - theta(1) = acos(0._r8) - wts (1) = 2._r8 - return - elseif(nlat==2) then - xx = sqrt(1._r8/3._r8) - theta(1) = acos( xx) - theta(2) = acos(-xx) - wts (1) = 1._r8 - wts (2) = 1._r8 - return - endif - - ! Proceed for nlat > 2 - !---------------------- - mnlat = mod(nlat,2) - ns2 = nlat/2 - nhalf = (nlat+1)/2 - - call sh_fourier_coefs_dp(nlat,cz,theta(ns2+1),wts(ns2+1)) - - dtheta = halfPI/nhalf - dthalf = dtheta/2._r8 - cmax = .2_r8*dtheta - - ! estimate first point next to theta = pi/2 - !------------------------------------------- - if(mnlat/=0) then - zero = halfPI - dtheta - zprev = halfPI - nix = nhalf - 1 - else - zero = halfPI - dthalf - nix = nhalf - endif - - do while(nix/=0) - dcor = huge(1._r8) - it = 0 - do while (abs(dcor) > eps*abs(zero)) - it = it + 1 - ! newton iterations - !----------------------- - call sh_legp_dlegp_theta(nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb) - dcor = pb/dpb - if(dcor.ne.0._r8) then - sgnd = dcor/abs(dcor) - else - sgnd = 1._r8 - endif - dcor = sgnd*min(abs(dcor),cmax) - zero = zero - dcor - end do - - theta(nix) = zero - zhold = zero - - ! wts(nix) = (nlat+nlat+1)/(dpb*dpb) - ! yakimiw's formula permits using old pb and dpb - !-------------------------------------------------- - wts(nix) = (nlat+nlat+1)/ (dpb+pb*dcos(zlast)/dsin(zlast))**2 - nix = nix - 1 - if(nix==nhalf-1) zero = 3._r8*zero - pi - if(nix0) then - cth = cdt - sth = sdt - do kk = 1,kdo - pb = pb + cp(kk)*cth - dpb = dpb - dcp(kk)*sth - chh = cdt*cth - sdt*sth - sth = sdt*cth + cdt*sth - cth = chh - end do - endif - else - ! n odd - !----------- - kdo = (nn+1)/2 - pb = 0._r8 - dpb = 0._r8 - cth = dcos(theta) - sth = dsin(theta) - do kk = 1,kdo - pb = pb + cp(kk)*cth - dpb = dpb - dcp(kk)*sth - chh = cdt*cth - sdt*sth - sth = sdt*cth + cdt*sth - cth = chh - end do - endif - - end subroutine sh_legp_dlegp_theta - !======================================================================= - -end module zonal_mean_mod diff --git a/tools/nudging/user_nl_cam-NUDGING_TEMPLATE b/tools/nudging/user_nl_cam-NUDGING_TEMPLATE index a8fa5980ef..39f7193982 100644 --- a/tools/nudging/user_nl_cam-NUDGING_TEMPLATE +++ b/tools/nudging/user_nl_cam-NUDGING_TEMPLATE @@ -24,11 +24,42 @@ ! current model timestep and shorter than the analyses ! timestep. As this number is increased, the nudging ! force has the form of newtonian cooling. -! Nudge_Uprof - INT index of profile structure to use for U. [0=OFF,1=ON,2=WINDOW] -! Nudge_Vprof - INT index of profile structure to use for V. [0=OFF,1=ON,2=WINDOW] -! Nudge_Tprof - INT index of profile structure to use for T. [0=OFF,1=ON,2=WINDOW] -! Nudge_Qprof - INT index of profile structure to use for Q. [0=OFF,1=ON,2=WINDOW] -! Nudge_PSprof - INT index of profile structure to use for PS. [0=OFF,N/A] +! Nudge_ZonalFilter - LOGICAL Option to apply zonal mean filtering to the +! model state and target data. +! Nudge_ZonalNbasis - INT The number of meridional modes(Legendre Polynomials) +! used for zonal filtering. +! Nudge_SpectralFilter - LOGICAL Option to apply spherical harminic filtering to +! the model state and target data so that nudging +! tendencies are only applied to scales larger than +! the specified truncation. +! Nudge_SpectralNtrunc - INT The number of meridional spherical harmonic modes used +! for spectral filtering. The nominal horizontal scale of +! the filtering can be estimated as: +! +! Hscale = PI*6350/Nudge_SpectralNtrunc +! +! i.e. Nudge_SpectralNtrunc=40 corresponds to a horizontal +! nudging scale Hscale~500km. +! +! Nudge_SpectralNring - INT The number of sampling rings used for local area averaging +! of spherical harmonic modes, to suppress sampling errors. +! When initializing each basis, a local average of SH values +! is computed for the area associated with each grid point. +! SpectralNring set the number of rings of equal-area points +! in this sampling domain. +! Each ring (kk) contains 8*(kk-1) sample points. +! +! Nudge_SpectralNring Number of Samping Points +! ------------------- ------------------------- +! 1 1 (DEFAULT SampleGrid NOT used) +! 2 9 +! 3 25 +! 4 49 +! +! Nudge_Uprof - INT index of profile structure to use for U. [0=OFF,1=ON,2=WINDOW,3=USER-WINDOW] +! Nudge_Vprof - INT index of profile structure to use for V. [0=OFF,1=ON,2=WINDOW,3=USER-WINDOW] +! Nudge_Tprof - INT index of profile structure to use for T. [0=OFF,1=ON,2=WINDOW,3=USER-WINDOW] +! Nudge_Qprof - INT index of profile structure to use for Q. [0=OFF,1=ON,2=WINDOW,3=USER-WINDOW] ! Nudge_Ucoef - REAL fractional nudging coeffcient for U. ! Utau=(Nudge_Ucoef/analyses_timestep) ! Nudge_Vcoef - REAL fractional nudging coeffcient for V. @@ -37,8 +68,13 @@ ! Ttau=(Nudge_Tcoef/analyses_timestep) ! Nudge_Qcoef - REAL fractional nudging coeffcient for Q. ! Qtau=(Nudge_Qcoef/analyses_timestep) +! Nudge_PSprof - INT index of profile structure to use for PS. [0=OFF,1=ON,2=WINDOW,3=USER-WINDOW] ! Nudge_PScoef - REAL fractional nudging coeffcient for PS. ! PStau=(Nudge_PScoef/analyses_timestep) +! Nudge_PSscal - REAL Coeffcient controling the vertical influence function. +! for: Nudge_PSscal < 1.0 : P'(n) = [ B(n) ]*(Target_Ps - Model_Ps) +! Nudge_PSscal > 1.0 : P'(n) = [(Pref/P0)**scl]*(Target_Ps - Model_Ps) +! (default: Nudge_PSscal = 0.0) ! Nudge_Beg_Year - INT nudging begining year. ! Nudge_Beg_Month - INT nudging begining month. ! Nudge_Beg_Day - INT nudging begining day. @@ -61,6 +97,7 @@ ! Nudge_Vwin_Invert - LOGICAL Invert Vertical Window Function to its Compliment. ! TRUE = value=0 inside the specified window, 1 outside ! FALSE = value=1 inside the specified window, 0 outside +! Nudge_Bwindow_File - CHAR path to the Horizonal Boundary Window file. ! / !====================================================== &nudging_nl @@ -71,6 +108,11 @@ Nudge_TimeScale_Opt= 0 Nudge_Times_Per_Day= 4 Model_Times_Per_Day= 48 + Nudge_ZonalFilter = .false. + Nudge_ZonalNbasis = -1 + Nudge_SpectralFilter = .false. + Nudge_SpectralNtrunc = -1 + Nudge_SpectralNring = 1 Nudge_Uprof =1 Nudge_Ucoef =1.00 Nudge_Vprof =1 @@ -81,6 +123,7 @@ Nudge_Qcoef =1.00 Nudge_PSprof =0 Nudge_PScoef =0.00 + Nudge_PSscal =0.00 Nudge_Beg_Year =1979 Nudge_Beg_Month=1 Nudge_Beg_Day =2