diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index a55f3ab13f..2c32abb0bd 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/src/physics/cam/nudging.F90 b/src/physics/cam/nudging.F90 index ced2ef57d2..158f0f192f 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,52 @@ subroutine nudging_timestep_init(phys_state) real(r8) :: Sbar,Qbar,Wsum integer :: dtime + real(r8),allocatable:: hyam_w(:,:,:) + real(r8),allocatable:: hybm_w(:,:,:) + real(r8),allocatable:: dadn_w(:,:,:) + real(r8),allocatable:: dbdn_w(:,:,:) + real(r8):: Astar (pcols,pver) + real(r8):: Bstar (pcols,pver) + real(r8):: Asum (pcols) + real(r8):: Bsum (pcols) + real(r8):: PSvert1(pcols) + real(r8):: PSvert2(pcols) + integer :: istat + + 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) + ! 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)) + + allocate(hyam_w(pcols,pver,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','hyam_w',pcols*pver*((endchunk-begchunk)+1)) + allocate(hybm_w(pcols,pver,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','hybm_w',pcols*pver*((endchunk-begchunk)+1)) + allocate(dadn_w(pcols,pver,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','dadn_w',pcols*pver*((endchunk-begchunk)+1)) + allocate(dbdn_w(pcols,pver,begchunk:endchunk),stat=istat) + call alloc_err(istat,'nudging_timestep_init','dbdn_w',pcols*pver*((endchunk-begchunk)+1)) + endif + ! Get time step size !-------------------- dtime = get_step_size() @@ -1176,7 +1606,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 +1630,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 +1783,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 +1843,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 +1970,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 +1979,10 @@ subroutine nudging_timestep_init(phys_state) ! End Routine !------------ - + if(allocated(hyam_w)) deallocate(hyam_w) + if(allocated(hybm_w)) deallocate(hybm_w) + if(allocated(dadn_w)) deallocate(dadn_w) + if(allocated(dbdn_w)) deallocate(dbdn_w) end subroutine nudging_timestep_init !================================================================ @@ -1468,6 +2032,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 +2123,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 +2140,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 +2157,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 +2174,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 +2191,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 +2262,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 +2320,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 +2460,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 +2481,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 +2584,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 6863799864..7edd6054ce 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..bd1cdd4c3e --- /dev/null +++ b/src/utils/ug_spectralmethods_mod.F90 @@ -0,0 +1,3767 @@ +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 = .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 = .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 + 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 + + ! 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') + endif + 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 + !----------------------------- + 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 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