diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 681200fd41..65ab4229fb 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -244,6 +244,7 @@ atm/waccm/ic/wa3_ne5np4_1950_spinup.cam2.i.1960-01-01-00000_c150810.nc atm/waccm/ic/FW2000_ne30_L70_01-01-0001_c200602.nc +atm/waccm/ic/FWsc2000climo_ne30pg3_L70_0002-01-01_c221103.nc atm/waccm/ic/FW2000_ne30pg3_L70_01-01-0001_c200602.nc atm/waccm/ic/FWsc2000_ne30pg3_L110_01-01-0001_c200521.nc diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index b122a0d1b7..afc9923da5 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -422,6 +422,24 @@ Default: .false., unless it is overridden (WACCM with interactive chemistry and configurations do this) + +Number of zonal mean basis functions (number of m=0 spherical harmonics) used in +Transformed Eulerian Mean (TEM) diagnostics + + + +Number of latitude grid points for zonal average TEM diagnostics history fields + + + + Frequency of TEM diagnostics calucation. + If > 0, frequency is specified as number of timesteps. + If < 0, frequency is specified as number of hours. + + Turn on verbose output identifying columns that fail energy/water diff --git a/bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml b/bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml index ff0b781607..fb1183f373 100644 --- a/bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml +++ b/bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml @@ -11,7 +11,7 @@ atm/cam/solar/SolarForcing1995-2005avg_c160929.nc -cesm2_init/f.e21.FWsc2000climo.f09_f09_mg17.cesm2.1-exp011.001/0003-01-01/f.e21.FWsc2000climo.f09_f09_mg17.cesm2.1-exp011.001.cam.i.0003-01-01-00000.nc +cesm2_init/f.e21.FWsc2000climo.f09_f09_mg17.cesm2.1-exp011.001_v2/0003-01-01/f.e21.FWsc2000climo.f09_f09_mg17.cesm2.1-exp011.001_v2.cam.i.0003-01-01-00000.nc atm/waccm/ic/f2000.waccm-mam3_1.9x2.5_L70.cam2.i.0017-01-01.c120410.nc @@ -91,7 +91,7 @@ -.true. +.true. .true. 1, 30, 120, 240, 240, 480, 365, 73, 30 @@ -109,7 +109,7 @@ - + 'BTAUN', 'BTAUS', 'BTAUE', 'BTAUW', 'BTAUNET', 'BUTEND1', 'BUTEND2', 'BUTEND3', 'BUTEND4', 'BUTEND5', 'BVTGWSPEC', 'MAXQ0', 'HDEPTH', 'NETDT', 'TAUN', 'TAUS', 'TAUE', 'TAUW', 'TAUGWX', 'TAUGWY', 'UTEND1', 'UTEND2', 'UTEND3', 'UTEND4', 'UTEND5', 'FRONTGF', 'FRONTGFA', 'EKGW', 'QNO', 'QRLNLTE', 'QRL_TOT', 'DUV', 'DVV', 'TTPXMLC' diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index fa19e6e00e..805c037249 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -738,6 +738,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands new file mode 100644 index 0000000000..bbc9a3e804 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands @@ -0,0 +1,2 @@ +./xmlchange ROF_NCPL=\$ATM_NCPL +./xmlchange LND_DOMAIN_FILE="domain.lnd.mpasa120_gx1v7.201215.nc" diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam new file mode 100644 index 0000000000..1212d35edd --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam @@ -0,0 +1,25 @@ +ncdata = '$DIN_LOC_ROOT/atm/waccm/ic/mpasa120km.waccm_fulltopo_c220818.nc' + +mpas_cam_coef = 0.2D0 +mpas_rayleigh_damp_u_timescale_days = 5.D0 +mpas_zd = 80000.0D0 +mpas_apvm_upwinding = 0.0D0 +mpas_dt = 600.D0 +mpas_dynamics_split_steps = 4 +mpas_epssm = 0.5D0 + +use_gw_front = .false. + +phys_grid_ctem_nfreq = -12 +phys_grid_ctem_zm_nbas = 120 +phys_grid_ctem_za_nlat = 90 + +fincl1 = ' ' +fexcl1 = ' ' + +fincl2 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' + +mfilt=1,1,1,1,1,1,1,1,1,1 +ndens=1,1,1,1,1,1,1,1,1,1 +nhtfrq=-24,-24,-24,-24,-24,-24,-24,-24,-24,-24 +inithist='ENDOFRUN' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm new file mode 100644 index 0000000000..a8dec89f15 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm @@ -0,0 +1,4 @@ +hist_nhtfrq = -24 +hist_mfilt = 1 +hist_ndens = 1 +fsurdat = '${DIN_LOC_ROOT}/lnd/clm2/surfdata_map/release-clm5.0.34/surfdata_mpasa120_hist_78pfts_CMIP6_simyr2000_c201215.nc' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands new file mode 100644 index 0000000000..eb40ad83e0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands @@ -0,0 +1,2 @@ +./xmlchange ROF_NCPL=\$ATM_NCPL +./xmlchange GLC_NCPL=\$ATM_NCPL diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam new file mode 100644 index 0000000000..55e353983a --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam @@ -0,0 +1,8 @@ +mfilt=1,1,1,1,1,1,1,1,1 +ndens=1,1,1,1,1,1,1,1,1 +nhtfrq=3,3,3,3,3,3,3,3,3 +inithist='ENDOFRUN' +phys_grid_ctem_nfreq=3 +phys_grid_ctem_zm_nbas=16 +phys_grid_ctem_za_nlat=15 +fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm new file mode 100644 index 0000000000..dfd79c6314 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm @@ -0,0 +1,26 @@ +!---------------------------------------------------------------------------------- +! Users should add all user specific namelist changes below in the form of +! namelist_var = new_namelist_value +! +! Include namelist variables for drv_flds_in ONLY if -megan and/or -drydep options +! are set in the CLM_NAMELIST_OPTS env variable. +! +! EXCEPTIONS: +! Set use_cndv by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting +! Set use_vichydro by the compset you use and the CLM_BLDNML_OPTS -vichydro setting +! Set use_cn by the compset you use and CLM_BLDNML_OPTS -bgc setting +! Set use_crop by the compset you use and CLM_BLDNML_OPTS -crop setting +! Set spinup_state by the CLM_BLDNML_OPTS -bgc_spinup setting +! Set irrigate by the CLM_BLDNML_OPTS -irrig setting +! Set dtime with L_NCPL option +! Set fatmlndfrc with LND_DOMAIN_PATH/LND_DOMAIN_FILE options +! Set finidat with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases +! (includes $inst_string for multi-ensemble cases) +! Set glc_grid with CISM_GRID option +! Set glc_smb with GLC_SMB option +! Set maxpatch_glcmec with GLC_NEC option +! Set glc_do_dynglacier with GLC_TWO_WAY_COUPLING env variable +!---------------------------------------------------------------------------------- +hist_nhtfrq = 3 +hist_mfilt = 1 +hist_ndens = 1 diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands new file mode 100644 index 0000000000..eb40ad83e0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands @@ -0,0 +1,2 @@ +./xmlchange ROF_NCPL=\$ATM_NCPL +./xmlchange GLC_NCPL=\$ATM_NCPL diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam new file mode 100644 index 0000000000..9cf4e0a97d --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam @@ -0,0 +1,8 @@ +mfilt=1,1,1,1,1,1,1,1,1 +ndens=1,1,1,1,1,1,1,1,1 +nhtfrq=9,9,9,9,9,9,9,9,9 +inithist='ENDOFRUN' +phys_grid_ctem_nfreq=-2 +phys_grid_ctem_zm_nbas=120 +phys_grid_ctem_za_nlat=90 +fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm new file mode 100644 index 0000000000..0d83b5367b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm @@ -0,0 +1,27 @@ +!---------------------------------------------------------------------------------- +! Users should add all user specific namelist changes below in the form of +! namelist_var = new_namelist_value +! +! Include namelist variables for drv_flds_in ONLY if -megan and/or -drydep options +! are set in the CLM_NAMELIST_OPTS env variable. +! +! EXCEPTIONS: +! Set use_cndv by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting +! Set use_vichydro by the compset you use and the CLM_BLDNML_OPTS -vichydro setting +! Set use_cn by the compset you use and CLM_BLDNML_OPTS -bgc setting +! Set use_crop by the compset you use and CLM_BLDNML_OPTS -crop setting +! Set spinup_state by the CLM_BLDNML_OPTS -bgc_spinup setting +! Set irrigate by the CLM_BLDNML_OPTS -irrig setting +! Set dtime with L_NCPL option +! Set fatmlndfrc with LND_DOMAIN_PATH/LND_DOMAIN_FILE options +! Set finidat with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases +! (includes $inst_string for multi-ensemble cases) +! Set glc_grid with CISM_GRID option +! Set glc_smb with GLC_SMB option +! Set maxpatch_glcmec with GLC_NEC option +! Set glc_do_dynglacier with GLC_TWO_WAY_COUPLING env variable +!---------------------------------------------------------------------------------- +hist_nhtfrq = 9 +hist_mfilt = 1 +hist_ndens = 1 + diff --git a/doc/ChangeLog b/doc/ChangeLog index 6342032df6..ea503ba1ce 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,114 @@ =============================================================== +Tag name: cam6_3_089 +Originator(s): fvitt, patc +Date: 10 Jan 2013 +One-line Summary: Introduce zonal mean and Transformed Eulerian Mean diagnostics capabilities on arbitrary grids +Github PR URL: https://github.com/ESCOMP/CAM/pull/677 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + + Add capability to calculate zonal mean of physics fields on arbatrary grids (structured or unstructured). + The synthesis of zonal mean values is based on m=0 spherical harmonics. (patc) + + Implement TEM circulation diagnostics in physics #653 + Need to add zonal mean grid to history for fields on arbitrary grids #652 + +Describe any changes made to build system: N/A + +Describe any changes made to the namelist: + + Namelist options added: + + phys_grid_ctem_zm_nbas: + Number of zonal mean basis functions (number of m=0 spherical harmonics) used in + Transformed Eulerian Mean (TEM) diagnostics + + phys_grid_ctem_za_nlat: + Number of latitude grid points for zonal average TEM diagnostics history fields + + phys_grid_ctem_nfreq: + Frequency of TEM diagnostics calucation. + If > 0, frequency is specified as number of timesteps. + If < 0, frequency is specified as number of hours. + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: + + The physics grid TEM diagnostics is expected have a small impact (on the order of a few percent) + when activated. + +Code reviewed by: peverwhee, cacraigucar, nusbaume + +List all files eliminated: N/A + +List all files added and what they do: + +A src/physics/cam/phys_grid_ctem.F90 + - computes terms of the TEM diags on the physics grid + +A src/utils/zonal_mean_mod.F90 + - computes zonal means on arbitrary physics grids based on m=0 order spherical harmonics + +A cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands +A cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam +A cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm +A cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands +A cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam +A cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam +A cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm + - for testing TEM diags on physics grids + +List all existing files that have been modified, and describe the changes: + +M bld/namelist_files/namelist_definition.xml + - new phys_grid_ctem options (see above) + +M bld/namelist_files/namelist_defaults_cam.xml + - update IC file for WACCM-SC on ne30pg3 grid (spun up stable IC file) + +M bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml + - correct IC file path for 0.9x1.25 grid + +M cime_config/testdefs/testlist_cam.xml + - add tests for TEM diags on physics grids + +M src/control/cam_comp.F90 + - invoke phys_grid_ctem_reg + +M src/control/runtime_opts.F90 + - invoke phys_grid_ctem_readnl + +M src/physics/cam/physpkg.F90 + - invoke : + phys_grid_ctem_init + phys_grid_ctem_final + phys_grid_ctem_diags + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +cheyenne/intel/aux_cam: All PASS + +izumi/nag/aux_cam: + DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details: + - pre-existing failure + + SMS_D_Ln6_Vnuopc.ne5_ne5_mg37.QPWmaC4.izumi_nag.cam-outfrq3s_physgrid_tem (Overall: DIFF) details: + - new test + +izumi/gnu/aux_cam: All PASS + +Summarize any changes to answers: bit-for-bit unchanged + +=============================================================== +=============================================================== + Tag name: cam6_3_088 Originator(s): cacraig Date: Dec 21, 2022 diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index f28cf66568..835fa3e452 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -92,6 +92,7 @@ subroutine cam_init( & #if (defined BFB_CAM_SCAM_IOP) use history_defaults, only: initialize_iop_history #endif + use phys_grid_ctem, only: phys_grid_ctem_reg ! Arguments character(len=cl), intent(in) :: caseid ! case ID @@ -168,6 +169,9 @@ subroutine cam_init( & ! Initialize physics grid decomposition call phys_grid_init() + ! Register zonal average grid for phys TEM diagnostics + call phys_grid_ctem_reg() + ! Register advected tracers and physics buffer fields call phys_register () diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 50ee489e71..356ec0f6d4 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -97,6 +97,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use qneg_module, only: qneg_readnl use lunar_tides, only: lunar_tides_readnl use upper_bc, only: ubc_readnl + use phys_grid_ctem, only: phys_grid_ctem_readnl !---------------------------Arguments----------------------------------- @@ -195,6 +196,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call dyn_readnl(nlfilename) call ionosphere_readnl(nlfilename) call qneg_readnl(nlfilename) + call phys_grid_ctem_readnl(nlfilename) end subroutine read_namelist diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 new file mode 100644 index 0000000000..e2d0cb5c71 --- /dev/null +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -0,0 +1,360 @@ +!---------------------------------------------------------------------------------- +! circulation diagnostics -- terms of the Transformed Eulerian Mean (TEM) equation +! +!---------------------------------------------------------------------------------- +module phys_grid_ctem + use shr_kind_mod, only: r8 => shr_kind_r8 + 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 physconst, only: pi + use cam_logfile, only: iulog + use cam_abortutils,only: endrun, handle_allocate_error + use namelist_utils,only: find_group_name + use spmd_utils, only: masterproc, mpi_integer, masterprocid, mpicom + use time_manager, only: get_step_size, get_nstep + + use shr_const_mod, only: rgas => shr_const_rgas ! J/K/kmole + use shr_const_mod, only: grav => shr_const_g ! m/s2 + use air_composition, only: mbarv ! g/mole + use string_utils, only: int2str + + implicit none + + private + public :: phys_grid_ctem_readnl + public :: phys_grid_ctem_reg + public :: phys_grid_ctem_init + public :: phys_grid_ctem_diags + public :: phys_grid_ctem_final + + type(ZonalMean_t) :: ZMobj + type(ZonalAverage_t) :: ZAobj + + integer :: nzalat = -huge(1) + integer :: nzmbas = -huge(1) + + integer :: ntimesteps = -huge(1) ! number of time steps bewteen TEM calculations + + logical :: do_tem_diags = .false. + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine phys_grid_ctem_readnl(nlfile) + character(len=*), intent(in) :: nlfile + integer :: ierr, unitn + + character(len=*), parameter :: prefix = 'phys_grid_ctem_readnl: ' + real(r8) :: dtime + + integer :: phys_grid_ctem_zm_nbas + integer :: phys_grid_ctem_za_nlat + integer :: phys_grid_ctem_nfreq + + namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat, phys_grid_ctem_nfreq + + phys_grid_ctem_zm_nbas = 0 + phys_grid_ctem_za_nlat = 0 + phys_grid_ctem_nfreq = 0 + + ! Read in namelist values + !------------------------ + if(masterproc) then + open(newunit=unitn, file=trim(nlfile), status='old') + call find_group_name(unitn, 'phys_grid_ctem_opts', status=ierr) + if(ierr == 0) then + read(unitn,phys_grid_ctem_opts,iostat=ierr) + if(ierr /= 0) then + call endrun(prefix//'ERROR reading namelist') + end if + end if + close(unitn) + end if + + call MPI_bcast(phys_grid_ctem_zm_nbas, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_zm_nbas') + call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_za_nlat') + call MPI_bcast(phys_grid_ctem_nfreq, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_nfreq') + + do_tem_diags = .false. + if (phys_grid_ctem_nfreq/=0) then + if (.not.(phys_grid_ctem_zm_nbas>0 .and. phys_grid_ctem_za_nlat>0)) then + call endrun(prefix//'inconsistent phys_grid_ctem namelist settings -- phys_grid_ctem_zm_nbas=' & + //int2str(phys_grid_ctem_zm_nbas)//', phys_grid_ctem_za_nlat='//int2str(phys_grid_ctem_za_nlat)) + end if + if (phys_grid_ctem_nfreq>0) then + ntimesteps = phys_grid_ctem_nfreq + else + dtime = get_step_size() + ntimesteps = nint( -phys_grid_ctem_nfreq*3600._r8/dtime ) + end if + if (ntimesteps<1) then + call endrun(prefix//'invalid ntimesteps -- phys_grid_ctem_nfreq needs to be a larger negative value ' & + //'or the model time step needs to be shorter') + end if + do_tem_diags = .true. + end if + + if (masterproc) then + if (do_tem_diags) then + write(iulog,*) 'TEM diagnostics will be calculated every ',ntimesteps,' time steps' + write(iulog,*) ' phys_grid_ctem_zm_nbas = ', phys_grid_ctem_zm_nbas + write(iulog,*) ' phys_grid_ctem_za_nlat = ', phys_grid_ctem_za_nlat + write(iulog,*) ' phys_grid_ctem_nfreq = ', phys_grid_ctem_nfreq + else + write(iulog,*) 'TEM diagnostics will not be performed' + end if + endif + + if (do_tem_diags) then + nzalat = phys_grid_ctem_za_nlat + nzmbas = phys_grid_ctem_zm_nbas + end if + + end subroutine phys_grid_ctem_readnl + + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + subroutine phys_grid_ctem_reg + + use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap, cam_grid_register + + type(horiz_coord_t), pointer :: zalon_coord + type(horiz_coord_t), pointer :: zalat_coord + integer(iMap), pointer :: grid_map(:,:) + + real(r8) :: zalats(nzalat) + real(r8) :: area(nzalat) + real(r8) :: zalons(1) + real(r8) :: dlatrad, dlatdeg, lat1, lat2 + real(r8) :: total_area + real(r8) :: total_wght + integer :: j, astat + + real(r8), parameter :: latdeg0 = -90._r8 + real(r8), parameter :: latrad0 = -pi*0.5_r8 + real(r8), parameter :: fourpi = pi*4._r8 + + integer, parameter :: ctem_zavg_phys_decomp = 201 ! Must be unique within CAM + + if (.not.do_tem_diags) return + + nullify(zalat_coord) + nullify(zalon_coord) + nullify(grid_map) + + zalons(1) = 0._r8 + + dlatrad = pi/real(nzalat,kind=r8) + dlatdeg = 180._r8/real(nzalat,kind=r8) + total_area = 0._r8 + total_wght = 0._r8 + + ! calculate latitudes and areas of zonal average grid boxes + do j = 1,nzalat + zalats(j) = latdeg0 + (real(j,kind=r8)-0.5_r8)*dlatdeg + lat1 = latrad0 + real(j-1,kind=r8)*dlatrad + lat2 = latrad0 + real(j ,kind=r8)*dlatrad + area(j) = 2._r8*pi*(sin(lat2)-sin(lat1)) + total_area = total_area + area(j) + total_wght = total_wght + 0.5_r8*(sin(lat2)-sin(lat1)) + end do + + ! sanity check + if ( abs(1._r8-total_wght)>1.e-12_r8 .or. abs(fourpi-total_area)>1.e-12_r8 ) then + call endrun('phys_grid_ctem_reg: problem with area/wght calc') + end if + + ! initialize zonal-average and zonal-mean utility objects + call ZAobj%init(zalats,area,nzalat,GEN_GAUSSLATS=.false.) + call ZMobj%init(nzmbas) + + ! Zonal average grid for history fields + + zalat_coord => horiz_coord_create('zalat', '', nzalat, 'latitude', 'degrees_north', 1, nzalat, zalats) + zalon_coord => horiz_coord_create('zalon', '', 1, 'longitude', 'degrees_east', 1, 1, zalons) + + ! grid decomposition map + allocate(grid_map(4,nzalat), stat=astat) + call handle_allocate_error(astat, 'phys_grid_ctem_reg', 'grid_map') + + do j = 1,nzalat + grid_map(1,j) = 1 + grid_map(2,j) = j + if (masterproc) then + grid_map(3,j) = 1 + grid_map(4,j) = j + else + grid_map(3,j) = 0 + grid_map(4,j) = 0 + end if + end do + + ! register the zonal average grid + call cam_grid_register('ctem_zavg_phys', ctem_zavg_phys_decomp, zalat_coord, zalon_coord, grid_map, & + unstruct=.false., zonal_grid=.true.) + + end subroutine phys_grid_ctem_reg + + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + subroutine phys_grid_ctem_init + + if (.not.do_tem_diags) return + + call addfld ('VTHzaphys',(/'lev'/), 'A', 'K m s-1','Meridional Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('WTHzaphys',(/'lev'/), 'A', 'K m s-1','Vertical Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('UVzaphys', (/'lev'/), 'A', 'm2 s-2', 'Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') + call addfld ('UWzaphys', (/'lev'/), 'A', 'm2 s-2', 'Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') + call addfld ('THphys', (/'lev'/), 'A', 'K', 'Potential temp', gridname='physgrid' ) + + end subroutine phys_grid_ctem_init + + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + subroutine phys_grid_ctem_diags(phys_state) + type(physics_state), intent(in) :: phys_state(begchunk:endchunk) + + character(len=*), parameter :: prefix = 'phys_grid_ctem_diags: ' + + real(r8) :: u(pcols,pver,begchunk:endchunk) + real(r8) :: v(pcols,pver,begchunk:endchunk) + real(r8) :: w(pcols,pver,begchunk:endchunk) + + real(r8) :: uzm(pcols,pver,begchunk:endchunk) + real(r8) :: vzm(pcols,pver,begchunk:endchunk) + real(r8) :: wzm(pcols,pver,begchunk:endchunk) + + real(r8) :: ud(pcols,pver,begchunk:endchunk) + real(r8) :: vd(pcols,pver,begchunk:endchunk) + real(r8) :: wd(pcols,pver,begchunk:endchunk) + real(r8) :: thd(pcols,pver,begchunk:endchunk) + + real(r8) :: uvp(pcols,pver,begchunk:endchunk) + real(r8) :: uwp(pcols,pver,begchunk:endchunk) + real(r8) :: vthp(pcols,pver,begchunk:endchunk) + real(r8) :: wthp(pcols,pver,begchunk:endchunk) + + integer :: lchnk, ncol, j, k + + ! potential temperature + real(r8) :: theta(pcols,pver,begchunk:endchunk) + real(r8) :: thzm(pcols,pver,begchunk:endchunk) + + real(r8) :: uvza(nzalat,pver) + real(r8) :: uwza(nzalat,pver) + real(r8) :: vthza(nzalat,pver) + real(r8) :: wthza(nzalat,pver) + + real(r8) :: sheight(pcols,pver) ! pressure scale height (m) + + if (.not.do_calc()) return + + do lchnk = begchunk,endchunk + + ncol = phys_state(lchnk)%ncol + + ! scale height + sheight(:ncol,:) = phys_state(lchnk)%t(:ncol,:) * rgas / ( mbarv(:ncol,:,lchnk) * grav ) ! meters + + ! potential temperature + theta(:ncol,:,lchnk) = phys_state(lchnk)%t(:ncol,:) * phys_state(lchnk)%exner(:ncol,:) + + ! vertical velocity + w(:ncol,:,lchnk) = -sheight(:ncol,:) * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) + + u(:ncol,:,lchnk) = phys_state(lchnk)%u(:ncol,:) + v(:ncol,:,lchnk) = phys_state(lchnk)%v(:ncol,:) + + end do + + ! zonal means evaluated on the physics grid (3D) to be used in the deviations calculation below + uzm(:,:,:) = zmean_fld(u(:,:,:)) + vzm(:,:,:) = zmean_fld(v(:,:,:)) + wzm(:,:,:) = zmean_fld(w(:,:,:)) + thzm(:,:,:) = zmean_fld(theta(:,:,:)) + + ! diagnostic output + do lchnk = begchunk, endchunk + call outfld( 'THphys', theta(:,:,lchnk), pcols, lchnk) + end do + + do lchnk = begchunk,endchunk + ncol = phys_state(lchnk)%ncol + do k = 1,pver + ! zonal deviations + thd(:ncol,k,lchnk) = theta(:ncol,k,lchnk) - thzm(:ncol,k,lchnk) + ud(:ncol,k,lchnk) = u(:ncol,k,lchnk) - uzm(:ncol,k,lchnk) + vd(:ncol,k,lchnk) = v(:ncol,k,lchnk) - vzm(:ncol,k,lchnk) + wd(:ncol,k,lchnk) = w(:ncol,k,lchnk) - wzm(:ncol,k,lchnk) + ! fluxes + uvp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * vd(:ncol,k,lchnk) + uwp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * wd(:ncol,k,lchnk) + vthp(:ncol,k,lchnk) = vd(:ncol,k,lchnk) * thd(:ncol,k,lchnk) + wthp(:ncol,k,lchnk) = wd(:ncol,k,lchnk) * thd(:ncol,k,lchnk) + end do + end do + + ! evaluate and output fluxes on the zonal-average grid + call ZAobj%binAvg(uvp, uvza) + call ZAobj%binAvg(uwp, uwza) + call ZAobj%binAvg(vthp, vthza) + call ZAobj%binAvg(wthp, wthza) + + if (any(abs(uvza)>1.e20_r8)) call endrun(prefix//'bad values in uvza') + if (any(abs(uwza)>1.e20_r8)) call endrun(prefix//'bad values in uwza') + if (any(abs(vthza)>1.e20_r8)) call endrun(prefix//'bad values in vthza') + if (any(abs(wthza)>1.e20_r8)) call endrun(prefix//'bad values in wthza') + + ! diagnostic output + do j = 1,nzalat + call outfld('UVzaphys',uvza(j,:),1,j) + call outfld('UWzaphys',uwza(j,:),1,j) + call outfld('VTHzaphys',vthza(j,:),1,j) + call outfld('WTHzaphys',wthza(j,:),1,j) + end do + + contains + + !------------------------------------------------------------------------------ + ! utility function for evaluating 3D zonal mean fields + !------------------------------------------------------------------------------ + function zmean_fld( fld ) result(fldzm) + + real(r8), intent(in) :: fld(pcols,pver,begchunk:endchunk) + + real(r8) :: fldzm(pcols,pver,begchunk:endchunk) + + real(r8) :: Zonal_Bamp3d(nzmbas,pver) + + call ZMobj%calc_amps(fld,Zonal_Bamp3d) + call ZMobj%eval_grid(Zonal_Bamp3d,fldzm) + + end function zmean_fld + + !------------------------------------------------------------------------------ + ! utility function returns TRUE when time to update TEM diags + !------------------------------------------------------------------------------ + logical function do_calc() + + integer :: nstep + nstep = get_nstep() + do_calc = do_tem_diags .and. mod(nstep,ntimesteps) == 0 + + end function do_calc + + end subroutine phys_grid_ctem_diags + + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + subroutine phys_grid_ctem_final + call ZAobj%final() + call ZMobj%final() + end subroutine phys_grid_ctem_final + +end module phys_grid_ctem diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 4f8cc56a6d..3003c9f42c 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -772,6 +772,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cam_snapshot, only: cam_snapshot_init use cam_history, only: addfld, register_vector_field, add_default use phys_control, only: phys_getopts + use phys_grid_ctem, only: phys_grid_ctem_init ! Input/output arguments type(physics_state), pointer :: phys_state(:) @@ -966,6 +967,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Initialize qneg3 and qneg4 call qneg_init() + ! Initialize phys TEM diagnostics + call phys_grid_ctem_init() + ! Initialize the snapshot capability call cam_snapshot_init(cam_in, cam_out, pbuf2d, begchunk) @@ -1303,6 +1307,7 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d ) use carma_intr, only : carma_final use wv_saturation, only : wv_sat_final use microp_aero, only : microp_aero_final + use phys_grid_ctem, only : phys_grid_ctem_final !----------------------------------------------------------------------- ! @@ -1325,6 +1330,7 @@ subroutine phys_final( phys_state, phys_tend, pbuf2d ) call carma_final call wv_sat_final call microp_aero_final() + call phys_grid_ctem_final() end subroutine phys_final @@ -2881,6 +2887,7 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) use iop_forcing, only: scam_use_iop_srf use nudging, only: Nudge_Model, nudging_timestep_init use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_timestep_init + use phys_grid_ctem, only: phys_grid_ctem_diags implicit none @@ -2954,6 +2961,9 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) !---------------------------------- if(Nudge_Model) call nudging_timestep_init(phys_state) + ! Update TEM diagnostics + call phys_grid_ctem_diags(phys_state) + end subroutine phys_timestep_init end module physpkg diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 new file mode 100644 index 0000000000..3c41a4b4ec --- /dev/null +++ b/src/utils/zonal_mean_mod.F90 @@ -0,0 +1,2099 @@ +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. + !==================================================================== + real(r8),parameter:: TINY = 1.e-20_r8 + ! + ! Passed Variables + !------------------ + real(r8),intent(in ):: I_Mat (:,:) + integer ,intent(in ):: Nbas + real(r8),intent(out):: O_InvMat(:,:) + ! + ! Local Values + !------------- + real(r8),allocatable:: Mwrk(:,:),Rscl(:) + integer ,allocatable:: Indx(:) + real(r8):: Psgn,Mmax,Mval,Sval + integer :: ii,jj,kk,ndx,i2,ii_max, astat + character(len=*), parameter :: subname = 'Invert_Matrix' + + ! Allocate work space + !--------------------- + allocate(Mwrk(Nbas,Nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Mwrk') + allocate(Rscl(Nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Rscl') + allocate(Indx(Nbas), stat=astat) + call handle_allocate_error(astat, subname, 'Indx') + + ! Copy the Input matrix so it can be decomposed + !------------------------------------------------- + Mwrk(1:Nbas,1:Nbas) = I_Mat(1:Nbas,1:Nbas) + + ! Initailize Row scales + !---------------------- + Psgn = 1._r8 + do ii=1,Nbas + Mmax = 0._r8 + do jj=1,Nbas + if(abs(Mwrk(ii,jj))>Mmax) Mmax = abs(Mwrk(ii,jj)) + end do + if(Mmax==0._r8) then + call endrun('Invert_Matrix: Singular Matrix') + endif + Rscl(ii) = 1._r8/Mmax + end do + + ! Decompose the matrix + !----------------------- + do jj=1,Nbas + + if(jj>1) then + do ii=1,(jj-1) + Sval = Mwrk(ii,jj) + if(ii>1) then + do kk=1,(ii-1) + Sval = Sval - Mwrk(ii,kk)*Mwrk(kk,jj) + end do + Mwrk(ii,jj) = Sval + endif + end do + endif + + Mmax = 0._r8 + do ii=jj,Nbas + Sval = Mwrk(ii,jj) + if(jj>1) then + do kk=1,(jj-1) + Sval = Sval - Mwrk(ii,kk)*Mwrk(kk,jj) + end do + Mwrk(ii,jj) = Sval + endif + Mval = Rscl(ii)*abs(Sval) + if(Mval>=Mmax) then + ii_max = ii + Mmax = Mval + endif + end do + + if(jj/=ii_max) then + do kk=1,Nbas + Mval = Mwrk(ii_max,kk) + Mwrk(ii_max,kk) = Mwrk(jj,kk) + Mwrk(jj,kk) = Mval + end do + Psgn = -Psgn + Rscl(ii_max) = Rscl(jj) + endif + + Indx(jj) = ii_max + if(jj/=Nbas) then + if(Mwrk(jj,jj)==0._r8) Mwrk(jj,jj) = TINY + Mval = 1._r8/Mwrk(jj,jj) + do ii=(jj+1),Nbas + Mwrk(ii,jj) = Mwrk(ii,jj)*Mval + end do + endif + + end do ! jj=1,Nbas + + if(Mwrk(Nbas,Nbas)==0._r8) Mwrk(Nbas,Nbas) = TINY + + ! Initialize the inverse array with the identity matrix + !------------------------------------------------------- + O_InvMat(:,:) = 0._r8 + do ii=1,Nbas + O_InvMat(ii,ii) = 1._r8 + end do + + ! Back substitution to construct the inverse + !--------------------------------------------- + do kk=1,Nbas + + i2 = 0 + do ii=11,Nbas + ndx = Indx(ii) + Sval = O_InvMat(ndx,kk) + O_InvMat(ndx,kk) = O_InvMat(ii,kk) + if(i2/=0) then + do jj=i2,(ii-1) + Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk) + end do + elseif(Sval/=0._r8) then + i2 = ii + endif + O_InvMat(ii,kk) = Sval + end do + + do ii=Nbas,1,-1 + Sval = O_InvMat(ii,kk) + if(iinn) 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