From 8e12246002c156b94b8e910d881db6dae1676db4 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 17 Oct 2022 16:33:31 -0600 Subject: [PATCH 01/24] initial implementation of physics grid TEM diagnostics new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm new file: src/physics/cam/phys_grid_ctem.F90 new file: src/utils/zonal_mean_mod.F90 modified: bld/namelist_files/namelist_definition.xml modified: src/control/cam_comp.F90 modified: src/control/runtime_opts.F90 modified: src/physics/cam/physpkg.F90 --- bld/namelist_files/namelist_definition.xml | 11 + .../cam/outfrq9s_physgrid_tem/shell_commands | 2 + .../cam/outfrq9s_physgrid_tem/user_nl_cam | 9 + .../cam/outfrq9s_physgrid_tem/user_nl_clm | 27 + src/control/cam_comp.F90 | 4 + src/control/runtime_opts.F90 | 2 + src/physics/cam/phys_grid_ctem.F90 | 341 +++ src/physics/cam/physpkg.F90 | 7 + src/utils/zonal_mean_mod.F90 | 2098 +++++++++++++++++ 9 files changed, 2501 insertions(+) create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm create mode 100644 src/physics/cam/phys_grid_ctem.F90 create mode 100644 src/utils/zonal_mean_mod.F90 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index d371f975f9..ef9e7afe31 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -422,6 +422,17 @@ 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 + + Turn on verbose output identifying columns that fail energy/water diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands new file mode 100644 index 0000000000..eb40ad83e0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_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/outfrq9s_physgrid_tem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam new file mode 100644 index 0000000000..378a76e449 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam @@ -0,0 +1,9 @@ +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_zm_nbas=16 +phys_grid_ctem_za_nlat=15 +fincl2 = 'VTHzm3d','WTHzm3d','UVzm3d','UWzm3d', + 'Uzm3d', 'Vzm3d', 'Wzm3d', 'THzm3d','THphys', +fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm new file mode 100644 index 0000000000..0d83b5367b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/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/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..d07e7db9ee --- /dev/null +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -0,0 +1,341 @@ +!---------------------------------------------------------------------------------- +! 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, pverp + use ref_pres, only: pref_edge + use interpolate_data, only: vertinterp + 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 + use namelist_utils,only: find_group_name + use spmd_utils, only: masterproc, mpi_integer, masterprocid, mpicom + + implicit none + + private + public :: phys_grid_ctem_readnl + public :: phys_grid_ctem_reg + public :: phys_grid_ctem_init + public :: phys_grid_ctem_diags + + type(ZonalMean_t) :: ZMobj + type(ZonalAverage_t) :: ZAobj + + integer :: nzalat = -huge(1) + integer :: nzmbas = -huge(1) + +contains + + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + subroutine phys_grid_ctem_readnl(nlfile) + character(len=*), intent(in) :: nlfile + integer :: ierr, unitn + + character(len=*), parameter :: prefix = 'phys_grid_ctem_readnl: ' + integer :: phys_grid_ctem_zm_nbas + integer :: phys_grid_ctem_za_nlat + + namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat + + ! 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) + call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr) + + if (masterproc) then + write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_zm_nbas: ',phys_grid_ctem_zm_nbas + write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_za_nlat: ',phys_grid_ctem_za_nlat + endif + + nzalat = phys_grid_ctem_za_nlat + nzmbas = phys_grid_ctem_zm_nbas + + 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 + + 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 + + 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 + + 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 + + if ( abs(1._r8-total_wght)>1.e-12_r8 .or. abs(fourpi-total_area)>1.e-12_r8 ) then + call endrun('zmean_phys_fields_reg: problem with area/wght calc') + end if + + call ZAobj%init(zalats,area,nzalat,GEN_GAUSSLATS=.false.) + call ZMobj%init(nzmbas) + + ! Zonal average grid + + 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)) + + 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 + + call addfld ('VTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('WTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('UVzaphys', (/'ilev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') + call addfld ('UWzaphys', (/'ilev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') + + call addfld ('VTHzm3d',(/'ilev' /), 'A','MK/S', 'Meridional Heat Flux: 3D zon. mean', gridname='physgrid' ) + call addfld ('WTHzm3d',(/'ilev' /), 'A','MK/S', 'Vertical Heat Flux: 3D zon. mean', gridname='physgrid' ) + call addfld ('UVzm3d', (/'ilev' /), 'A','M2/S2','Meridional Flux of Zonal Momentum: 3D zon. mean', gridname='physgrid' ) + call addfld ('UWzm3d', (/'ilev' /), 'A','M2/S2','Vertical Flux of Zonal Momentum: 3D zon. mean', gridname='physgrid' ) + + call addfld ('Uzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean zonal wind - defined on ilev', gridname='physgrid') + call addfld ('Vzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean meridional wind - defined on ilev', gridname='physgrid' ) + call addfld ('Wzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean vertical wind - defined on ilev', gridname='physgrid' ) + call addfld ('THzm3d', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) + call addfld ('THphys', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) + + end subroutine phys_grid_ctem_init + + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + subroutine phys_grid_ctem_diags(phys_state) + type(physics_state), intent(in) :: phys_state(begchunk:endchunk) + + real(r8) :: ui(pcols,pverp,begchunk:endchunk) + real(r8) :: vi(pcols,pverp,begchunk:endchunk) + real(r8) :: wi(pcols,pverp,begchunk:endchunk) + + real(r8) :: uzm(pcols,pverp,begchunk:endchunk) + real(r8) :: vzm(pcols,pverp,begchunk:endchunk) + real(r8) :: wzm(pcols,pverp,begchunk:endchunk) + + real(r8) :: ud(pcols,pverp,begchunk:endchunk) + real(r8) :: vd(pcols,pverp,begchunk:endchunk) + real(r8) :: wd(pcols,pverp,begchunk:endchunk) + real(r8) :: thd(pcols,pverp,begchunk:endchunk) + + real(r8) :: uvp(pcols,pverp,begchunk:endchunk) + real(r8) :: uwp(pcols,pverp,begchunk:endchunk) + real(r8) :: vthp(pcols,pverp,begchunk:endchunk) + real(r8) :: wthp(pcols,pverp,begchunk:endchunk) + + real(r8) :: uv(pcols,pverp,begchunk:endchunk) + real(r8) :: uw(pcols,pverp,begchunk:endchunk) + real(r8) :: vth(pcols,pverp,begchunk:endchunk) + real(r8) :: wth(pcols,pverp,begchunk:endchunk) + + integer :: lchnk, ncol, j, k + real(r8) :: fld_tmp(pcols,pverp) + + real(r8) :: theta(pcols,pver,begchunk:endchunk) ! potential temperature + real(r8) :: thi(pcols,pverp,begchunk:endchunk) + real(r8) :: thzm(pcols,pverp,begchunk:endchunk) + + real(r8) :: w(pcols,pver,begchunk:endchunk) + + real(r8) :: uvza(nzalat,pverp) + real(r8) :: uwza(nzalat,pverp) + real(r8) :: vthza(nzalat,pverp) + real(r8) :: wthza(nzalat,pverp) + + real(r8), parameter :: hscale = 7000._r8 ! pressure scale height + + ui(:,:,:) = 0._r8 + vi(:,:,:) = 0._r8 + wi(:,:,:) = 0._r8 + thi(:,:,:) = 0._r8 + + uzm(:,:,:) = 0._r8 + vzm(:,:,:) = 0._r8 + wzm(:,:,:) = 0._r8 + thzm(:,:,:) = 0._r8 + + ud(:,:,:) = 0._r8 + vd(:,:,:) = 0._r8 + uvp(:,:,:) = 0._r8 + + do lchnk = begchunk,endchunk + + ncol = phys_state(lchnk)%ncol + + theta(:ncol,:,lchnk) = phys_state(lchnk)%t(:ncol,:) * phys_state(lchnk)%exner(:ncol,:) + w(:ncol,:,lchnk) = - hscale * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) + + do k = 1,pverp + call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%u(:,:), ui(:,k,lchnk) ) + call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%v(:,:), vi(:,k,lchnk) ) + call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), theta(:,:,lchnk), thi(:,k,lchnk) ) + call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), w(:,:,lchnk), wi(:,k,lchnk) ) + end do + + end do + + ! these need to be evaluated on the physics grid (3D) + ! to be used in the deviations calculation below + uzm(:,:,:) = zmean_fld(ui(:,:,:)) + vzm(:,:,:) = zmean_fld(vi(:,:,:)) + wzm(:,:,:) = zmean_fld(wi(:,:,:)) + thzm(:,:,:) = zmean_fld(thi(:,:,:)) + + do lchnk = begchunk, endchunk + ncol = phys_state(lchnk)%ncol + + fld_tmp(:ncol,:) = thi(:ncol,:,lchnk) + call outfld( 'THphys', fld_tmp(:ncol,:), ncol, lchnk) + + fld_tmp(:ncol,:) = thzm(:ncol,:,lchnk) + call outfld( 'THzm3d', fld_tmp(:ncol,:), ncol, lchnk) + + fld_tmp(:ncol,:) = uzm(:ncol,:,lchnk) + call outfld( 'Uzm3d', fld_tmp(:ncol,:), ncol, lchnk) + fld_tmp(:ncol,:) = vzm(:ncol,:,lchnk) + call outfld( 'Vzm3d', fld_tmp(:ncol,:), ncol, lchnk) + fld_tmp(:ncol,:) = wzm(:ncol,:,lchnk) + call outfld( 'Wzm3d', fld_tmp(:ncol,:), ncol, lchnk) + end do + + + do lchnk = begchunk,endchunk + ncol = phys_state(lchnk)%ncol + do k = 1,pverp + ! zonal deviations + thd(:ncol,k,lchnk) = thi(:ncol,k,lchnk) - thzm(:ncol,k,lchnk) + ud(:ncol,k,lchnk) = ui(:ncol,k,lchnk) - uzm(:ncol,k,lchnk) + vd(:ncol,k,lchnk) = vi(:ncol,k,lchnk) - vzm(:ncol,k,lchnk) + wd(:ncol,k,lchnk) = wi(: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 these 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) + + 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 + + ! not needed --- only for sanity checks + uv(:,:,:) = zmean_fld(uvp(:,:,:)) + uw(:,:,:) = zmean_fld(uwp(:,:,:)) + vth(:,:,:) = zmean_fld(vthp(:,:,:)) + wth(:,:,:) = zmean_fld(wthp(:,:,:)) + + do lchnk = begchunk, endchunk + ncol = phys_state(lchnk)%ncol + fld_tmp(:ncol,:) = uv(:ncol,:,lchnk) + call outfld( 'UVzm3d', fld_tmp(:ncol,:), ncol, lchnk) + fld_tmp(:ncol,:) = uw(:ncol,:,lchnk) + call outfld( 'UWzm3d', fld_tmp(:ncol,:), ncol, lchnk) + fld_tmp(:ncol,:) = vth(:ncol,:,lchnk) + call outfld( 'VTHzm3d', fld_tmp(:ncol,:), ncol, lchnk) + fld_tmp(:ncol,:) = wth(:ncol,:,lchnk) + call outfld( 'WTHzm3d', fld_tmp(:ncol,:), ncol, lchnk) + end do + + contains + + !------------------------------------------------------------------------------ + ! utility function for evaluating 3D zonal mean fields + !------------------------------------------------------------------------------ + function zmean_fld( fld ) result(fldzm) + + real(r8), intent(in) :: fld(pcols,pverp,begchunk:endchunk) + + real(r8) :: fldzm(pcols,pverp,begchunk:endchunk) + + real(r8) :: Zonal_Bamp3d(nzmbas,pverp) + + call ZMobj%calc_amps(fld,Zonal_Bamp3d) + call ZMobj%eval_grid(Zonal_Bamp3d,fldzm) + + end function zmean_fld + + end subroutine phys_grid_ctem_diags + +end module phys_grid_ctem diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 01e05c48db..440e424a1d 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -770,6 +770,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(:) @@ -964,6 +965,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) @@ -2879,6 +2883,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 @@ -2952,6 +2957,8 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) !---------------------------------- if(Nudge_Model) call nudging_timestep_init(phys_state) + 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..a2acd602b7 --- /dev/null +++ b/src/utils/zonal_mean_mod.F90 @@ -0,0 +1,2098 @@ +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/synthsis of zonal mean values +! on a 2D grid of points distributed over the +! surface of a sphere. +! ZonalProfile_t: For the analysis/synthsis 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 gridpoints +! 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 +! the ZonalProfile_t grid and vise-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 a 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/NOT-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 gridpopints, initialize +! the ZonalAverage datastruture 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 + use spmd_utils, only: mpicom + + use phys_grid, only: ngcols_p => num_global_phys_cols + + 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 + 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 + 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 + end type ZonalAverage_t + + real(r8), parameter :: PI2 = 2._r8*atan(1._r8) ! pi/2 + +contains + !======================================================================= + subroutine init_ZonalMean(this,I_nbas) + ! + ! init_ZonalMean: Initialize the ZonalMean datastruture for the + ! given physgrid gridpoints. It is assumed that the domain + ! of these gridpoints spans the surface of the sphere. + ! The representation of basis functions 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 + + 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)) + allocate(this%basis(pcols,begchunk:endchunk,I_nbas)) + allocate(this%map (I_nbas,I_nbas)) + 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)) + allocate(Bcoef(I_nbas)) + allocate(Csum (nlcols, Cvec_len)) + allocate(Cvec (Cvec_len)) + allocate(Bsum (nlcols, I_nbas)) + allocate(Bnorm(I_nbas)) + allocate(Bcov (I_nbas,I_nbas)) + + Bsum(:,:) = 0._r8 + Csum(:,:) = 0._r8 + + ! Save a copy 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 + PI2 + end do + end do + + ! Add first basis for the mean values. + !------------------------------------------ + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + this%basis(cc,lchnk,1) = 1._r8/sqrt(8._r8*PI2) + end do + end do + + ! Loop over the remaining basis functions + !--------------------------------------- + do nn=2,this%nbas + nb = nn-1 + + ! Generate coefs for the basis + !------------------------------ + call dalfk(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 dlfpt(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn)) + end do + end do + end do ! nn=1,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) + do cc = 1,ncols + this%basis(cc,lchnk,nn) = this%basis(cc,lchnk,nn)/sqrt(Bnorm(nn)) + end do + 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,I_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 + + nlcols = get_nlcols_p() + + allocate(Gcov(this%nbas)) + allocate(Csum(nlcols, this%nbas)) + 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(:,:,:) + 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 + + integer :: nlev + + nlev = size(I_Gdata,dim=2) + + nlcols = get_nlcols_p() + allocate(Gcov(this%nbas)) + allocate(Csum(nlcols, this%nbas)) + + 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-begchunk+1)*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(:,:,:) + ! + ! 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-begchunk+1) = O_Gdata(cc,ll,lchnk-begchunk+1) + (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 datastruture 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 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 + logical :: generate_lats + + 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)) + allocate(this%basis(I_nlat,I_nbas)) + allocate(this%map (I_nbas,I_nbas)) + + allocate(Clats(I_nlat)) + allocate(Bcoef(I_nbas)) + allocate(Bcov (I_nbas,I_nbas)) + + ! Optionally create the Latitude Gridpoints + ! and their associated area weights. Otherwise it + ! is assumed that the user is supplying them. + !----------------------------------------------- + if(generate_lats) then + + ! Create a Gaussian grid from SP to NP + !-------------------------------------- + call dgaqd(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)/atan(1._r8)) - 90._r8 + IO_area(nn) = IO_area(nn)*(8._r8*atan(1._r8)) + end do + else + ! Convert Latitudes to SP->NP colatitudes in radians + !---------------------------------------------------- + do nn=1,I_nlat + Clats(nn) = (IO_lats(nn) + 90._r8)*atan(1._r8)/45._r8 + end do + endif + + ! Copy the area weights for each nlat + ! gridpoint to the datastructure + !--------------------------------------- + this%area(1:I_nlat) = IO_area(1:I_nlat) + + ! Add first basis for the mean values. + !------------------------------------------ + this%basis(:,1) = 1._r8/sqrt(16._r8*atan(1._r8)) + 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 dalfk(nb,0,Bcoef) + + ! Create an un-normalized basis for the + ! coefs at each nlat gridpoint + !--------------------------------------- + do ii=1,I_nlat + call dlfpt(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 + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + allocate(Gcov(this%nbas)) + 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,ll + + integer :: nlev + + nlev = size(I_Zdata,dim=2) + + ! Compute Covariance with input data and basis functions + !-------------------------------------------------------- + allocate(Gcov(this%nbas,nlev)) + do ll=1,nlev + do nn=1,this%nbas + Gcov(nn,ll) = 0._r8 + do ii=1,this%nlat + Gcov(nn,ll) = Gcov(nn,ll) + (I_Zdata(ii,ll)*this%basis(ii,nn)*this%area(ii)) + end do + end do + end do + + ! 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 + O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(n2,ll) + 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,ll + + 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 ll=1,nlev + do ii=1,this%nlat + O_Zdata(ii,ll) = O_Zdata(ii,ll) + (I_Bamp(nn,ll)*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) + ! + ! init_ZonalAverage: Initialize the ZonalAverage datastruture 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 + integer :: ncols,lchnk,cc,jlat + integer :: nlcols, count + logical :: generate_lats + + 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)) + allocate(this%a_norm (I_nlat)) + allocate(this%area_g (pcols,begchunk:endchunk)) + allocate(this%idx_map(pcols,begchunk:endchunk)) + + allocate(Clats (I_nlat)) + allocate(BinLat(I_nlat+1)) + allocate(Glats (pcols,begchunk:endchunk)) + allocate(Asum (nlcols,I_nlat)) + allocate(Anorm (I_nlat)) + + ! Optionally create the Latitude Gridpoints + ! and their associated area weights. Otherwise it + ! is assumed that the user is supplying them. + !----------------------------------------------- + if(generate_lats) then + + ! Create a Gaussin grid from SP to NP + !-------------------------------------- + call dgaqd(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)/atan(1._r8)) - 90._r8 + IO_area(nn) = IO_area(nn)*(8._r8*atan(1._r8)) + end do + else + ! Convert Latitudes to SP->NP colatitudes in radians + !---------------------------------------------------- + do nn=1,this%nlat + Clats(nn) = (IO_lats(nn) + 90._r8)*atan(1._r8)/45._r8 + end do + endif + + ! Copy the Lat grid area weights to the datastructure + !----------------------------------------------------- + this%area(1:this%nlat) = IO_area(1:this%nlat) + + ! Save a copy 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 + PI2 + end do + end do + + ! Set boundaroes for Latitude bins + !----------------------------------- + BinLat(1) = 0._r8 + BinLat(this%nlat+1) = 4._r8*atan(1._r8) + 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).le.BinLat(2)).and. & + (Glats(cc,lchnk).ge.BinLat(1)) ) then + jlat = 1 + elseif((Glats(cc,lchnk).ge.BinLat(this%nlat) ).and. & + (Glats(cc,lchnk).le.BinLat(this%nlat+1)) ) then + jlat = this%nlat + else + do jj=2,(this%nlat-1) + if((Glats(cc,lchnk).gt.BinLat(jj )).and. & + (Glats(cc,lchnk).le.BinLat(jj+1)) ) then + jlat = jj + exit + endif + end do + endif + if((jlat.lt.1).or.(jlat.gt.this%nlat)) 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 + print*, 'ZonalAverage init ERROR: Anorm; ',Anorm(:) + call endrun('init_ZonalAverage: error in Anorm') + 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_ZonalProfile_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 + + nlcols = get_nlcols_p() + + + ! Initialize Zonal profile + !--------------------------- + allocate(Asum(nlcols,this%nlat)) + 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 + + end subroutine calc_ZonalAverage_2DbinAvg + !======================================================================= + + + !======================================================================= + subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) + ! + ! calc_ZonalProfile_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(:,:,:) + real(r8),intent(out):: O_Zdata(:,:) + ! + ! Local Values + !-------------- + real(r8),allocatable:: Gsum(:) + real(r8),allocatable:: Asum(:,:) + integer:: nn,ncols,lchnk,cc,jlat + integer:: Nsum,ll,ns + + integer :: nlev + integer :: nlcols, count + + nlev = size(I_Gdata,dim=2) + nlcols = get_nlcols_p() + + ! Initialize Zonal profile + !--------------------------- + Nsum = this%nlat*nlev + allocate(Gsum(Nsum)) + allocate(Asum(nlcols,Nsum)) + Asum(:,:) = 0._r8 + + O_Zdata(1:this%nlat,1:nlev) = 0._r8 + + ! Compute area-weighted sums + !----------------------------- + do ll = 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 + (ll-1)*this%nlat + count=count+1 + Asum(count,ns) = I_Gdata(cc,ll,lchnk-begchunk+1)*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 ll = 1,nlev + do nn = 1,this%nlat + ns = nn + (ll-1)*this%nlat + O_Zdata(nn,ll) = Gsum(ns)/this%a_norm(nn) + end do + end do + + end subroutine calc_ZonalAverage_3DbinAvg + !======================================================================= + + + !======================================================================= + subroutine dalfk(nn,mm,cp) + ! + ! subroutine alfk (n,m,cp) + ! + ! dimension of real cp(n/2 + 1) + ! arguments + ! + ! purpose routine alfk computes single precision fourier + ! coefficients in the trigonometric series + ! representation of the normalized associated + ! legendre function pbar(n,m,theta) for use by + ! routines lfp and lfpt in calculating single + ! precision pbar(n,m,theta). + ! + ! first define the normalized associated + ! legendre functions + ! + ! pbar(m,n,theta) = sqrt((2*n+1)*factorial(n-m) + ! /(2*factorial(n+m)))*sin(theta)**m/(2**n* + ! factorial(n)) times the (n+m)th derivative of + ! (x**2-1)**n with respect to x=cos(theta) + ! + ! where theta is colatitude. + ! + ! then subroutine alfk computes the coefficients + ! cp(k) in the following trigonometric + ! expansion of pbar(m,n,theta). + ! + ! 1) for n even and m even, pbar(m,n,theta) = + ! .5*cp(1) plus the sum from k=1 to k=n/2 + ! of cp(k+1)*cos(2*k*th) + ! + ! 2) for n even and m odd, pbar(m,n,theta) = + ! the sum from k=1 to k=n/2 of + ! cp(k)*sin(2*k*th) + ! + ! 3) for n odd and m even, pbar(m,n,theta) = + ! the sum from k=1 to k=(n+1)/2 of + ! cp(k)*cos((2*k-1)*th) + ! + ! 4) for n odd and m odd, pbar(m,n,theta) = + ! the sum from k=1 to k=(n+1)/2 of + ! cp(k)*sin((2*k-1)*th) + ! + ! + ! usage call alfk(n,m,cp) + ! + ! arguments + ! + ! on input n + ! nonnegative integer specifying the degree of + ! pbar(n,m,theta) + ! + ! m + ! is the order of pbar(n,m,theta). m can be + ! any integer however cp is computed such that + ! pbar(n,m,theta) = 0 if abs(m) is greater + ! than n and pbar(n,m,theta) = (-1)**m* + ! pbar(n,-m,theta) for negative m. + ! + ! on output cp + ! single precision array of length (n/2)+1 + ! which contains the fourier coefficients in + ! the trigonometric series representation of + ! pbar(n,m,theta) + ! + ! + ! special conditions none + ! + ! precision single + ! + ! algorithm the highest order coefficient is determined in + ! closed form and the remainig coefficients are + ! determined as the solution of a backward + ! recurrence relation. + ! + ! accuracy comparison between routines alfk and double + ! precision dalfk on the cray1 indicates + ! greater accuracy for smaller values + ! of input parameter n. agreement to 14 + ! places was obtained for n=10 and to 13 + ! places for n=100. + ! + !===================================================================== + ! + ! 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,ll + + 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.gt.nn) return + + if((nn-1).lt.0) then + cp(1) = sqrt(2._r8) + return + elseif((nn-1).eq.0) then + if(ma.ne.0) then + cp(1) = sqrt(.75_r8) + if(mm.eq.-1) cp(1) = -cp(1) + else + cp(1) = sqrt(1.5_r8) + endif + return + else + if(mod(nn+ma,2).ne.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.ge.1) then + do ii = 1,nmms2 + t1 = fnum*t1/fden + if(t1.GT.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).ne.0) then + t1 = -t1/2._r8**(nn-1-nex) + else + t1 = t1/2._r8**(nn-1-nex) + endif + t2 = 1._r8 + if(ma.ne.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).eq.0).and.(mod(ma,2).eq.0)) then + ll = 1+(nn+1)/2 + else + ll = (nn+1)/2 + endif + + cp(ll) = cp2 + if(mm.lt.0) then + if(mod(ma,2).ne.0) cp(ll) = -cp(ll) + endif + if(ll.le.1) return + + fk = nn + a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 + b1 = 2._r8* (fk*fk-fnmsq) + cp(ll-1) = b1*cp(ll)/a1 + 30 continue + ll = ll - 1 + if(ll.le.1) return + 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(ll-1) = -(b1*cp(ll)+c1*cp(ll+1))/a1 + goto 30 + + end subroutine dalfk + !======================================================================= + + + !======================================================================= + subroutine dlfpt(nn,mm,theta,cp,pb) + ! + ! subroutine lfpt (n,m,theta,cp,pb) + ! + ! dimension of + ! arguments + ! cp((n/2)+1) + ! + ! purpose routine lfpt uses coefficients computed by + ! routine alfk to compute the single precision + ! normalized associated legendre function pbar(n, + ! m,theta) at colatitude theta. + ! + ! usage call lfpt(n,m,theta,cp,pb) + ! + ! arguments + ! + ! on input n + ! nonnegative integer specifying the degree of + ! pbar(n,m,theta) + ! m + ! is the order of pbar(n,m,theta). m can be + ! any integer however pbar(n,m,theta) = 0 + ! if abs(m) is greater than n and + ! pbar(n,m,theta) = (-1)**m*pbar(n,-m,theta) + ! for negative m. + ! + ! theta + ! single precision colatitude in radians + ! + ! cp + ! single precision array of length (n/2)+1 + ! containing coefficients computed by routine + ! alfk + ! + ! on output pb + ! single precision variable containing + ! pbar(n,m,theta) + ! + ! special conditions calls to routine lfpt must be preceded by an + ! appropriate call to routine alfk. + ! + ! precision single + ! + ! algorithm the trigonometric series formula used by + ! routine lfpt to calculate pbar(n,m,th) at + ! colatitude th depends on m and n as follows: + ! + ! 1) for n even and m even, the formula is + ! .5*cp(1) plus the sum from k=1 to k=n/2 + ! of cp(k)*cos(2*k*th) + ! 2) for n even and m odd. the formula is + ! the sum from k=1 to k=n/2 of + ! cp(k)*sin(2*k*th) + ! 3) for n odd and m even, the formula is + ! the sum from k=1 to k=(n+1)/2 of + ! cp(k)*cos((2*k-1)*th) + ! 4) for n odd and m odd, the formula is + ! the sum from k=1 to k=(n+1)/2 of + ! cp(k)*sin((2*k-1)*th) + ! + ! accuracy comparison between routines lfpt and double + ! precision dlfpt on the cray1 indicates greater + ! accuracy for greater values on input parameter + ! n. agreement to 13 places was obtained for + ! n=10 and to 12 places for n=100. + ! + ! timing time per call to routine lfpt is dependent on + ! the input parameter n. + ! + !===================================================================== + 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.gt.nn) return + + if(nn.le.0) then + if(ma.le.0) then + pb = sqrt(.5_r8) + goto 140 + endif + endif + + nmod = mod(nn,2) + mmod = mod(ma,2) + + if(nmod.le.0) then + if(mmod.le.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 + goto 140 + 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 + goto 140 + endif + + kdo = (nn+1)/2 + if(mmod.le.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 + goto 140 + 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 + +140 continue + + end subroutine dlfpt + !======================================================================= + + + !======================================================================= + 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,ll,i2,ii_max + + ! Allocate work space + !--------------------- + allocate(Mwrk(Nbas,Nbas)) + allocate(Rscl(Nbas)) + allocate(Indx(Nbas)) + + ! 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)).gt.Mmax) Mmax = abs(Mwrk(ii,jj)) + end do + if(Mmax.eq.0._r8) then + call endrun('Singular Matrix') + endif + Rscl(ii) = 1._r8/Mmax + end do + + ! Decompose the matrix + !----------------------- + do jj=1,Nbas + + if(jj.gt.1) then + do ii=1,(jj-1) + Sval = Mwrk(ii,jj) + if(ii.gt.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.gt.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.ge.Mmax) then + ii_max = ii + Mmax = Mval + endif + end do + + if(jj.ne.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.ne.Nbas) then + if(Mwrk(jj,jj).eq.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).eq.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 + ll = Indx(ii) + Sval = O_InvMat(ll,kk) + O_InvMat(ll,kk) = O_InvMat(ii,kk) + if(i2.ne.0) then + do jj=i2,(ii-1) + Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk) + end do + elseif(Sval.ne.0._r8) then + i2 = ii + endif + O_InvMat(ii,kk) = Sval + end do + + do ii=Nbas,1,-1 + Sval = O_InvMat(ii,kk) + if(ii.lt.Nbas) then + do jj=(ii+1),Nbas + Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk) + end do + endif + O_InvMat(ii,kk) = Sval/Mwrk(ii,ii) + end do + + end do ! kk=1,Nbas + + ! Clean up this mess + !--------------------- + deallocate(Mwrk) + deallocate(Rscl) + deallocate(Indx) + + end subroutine Invert_Matrix + !======================================================================= + + + !======================================================================= + subroutine dgaqd(nlat,theta,wts,ierr) + ! + ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + ! . . + ! . 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 that has been accepted for publication. + ! This routine is faster and more accurate than older program + ! with the same name. + ! + ! subroutine gaqd 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.le.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):: eps + real(r8):: sgnd + real(r8):: xx,pi,pis2,dtheta,dthalf + real(r8):: cmax,zprev,zlast,zero,zhold,pb,dpb,dcor,summ,cz + integer :: mnlat,ns2,nhalf,nix,it,ii + + ! check work space length + !------------------------ + if(nlat.le.0) then + return + ierr = 1 + endif + ierr = 0 + + ! compute weights and points analytically when nlat=1,2 + !------------------------------------------------------- + if(nlat.eq.1) then + theta(1) = acos(0._r8) + wts (1) = 2._r8 + return + elseif(nlat.eq.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 + !---------------------- + eps = sqrt(ddzeps(1._r8)) + eps = eps*sqrt(eps) + pis2 = 2._r8*atan(1._r8) + pi = pis2 + pis2 + mnlat = mod(nlat,2) + ns2 = nlat/2 + nhalf = (nlat+1)/2 + + call dcpdp(nlat,cz,theta(ns2+1),wts(ns2+1)) + + dtheta = pis2/nhalf + dthalf = dtheta/2._r8 + cmax = .2_r8*dtheta + + ! estimate first point next to theta = pi/2 + !------------------------------------------- + if(mnlat.ne.0) then + zero = pis2 - dtheta + zprev = pis2 + nix = nhalf - 1 + else + zero = pis2 - dthalf + nix = nhalf + endif + + 10 continue + it = 0 + 20 continue + it = it + 1 + zlast = zero + + ! newton iterations + !----------------------- + call dtpdp(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 + if(abs(zero-zlast).gt.eps*abs(zero)) goto 20 + + 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.eq.0) goto 30 + if(nix.eq.nhalf-1) zero = 3._r8*zero - pi + if(nix.lt.nhalf-1) zero = zero + zero - zprev + zprev = zhold + goto 10 + 30 continue + + ! extend points and weights via symmetries + !------------------------------------------- + if(mnlat.ne.0) then + theta(nhalf) = pis2 + call dtpdp(nlat,pis2,cz,theta(ns2+1),wts(ns2+1),pb,dpb) + wts(nhalf) = (nlat+nlat+1)/ (dpb*dpb) + endif + + do ii = 1,ns2 + wts (nlat-ii+1) = wts(ii) + theta(nlat-ii+1) = pi - theta(ii) + end do + + summ = 0._r8 + do ii = 1,nlat + summ = summ + wts(ii) + end do + do ii = 1,nlat + wts(ii) = 2._r8*wts(ii)/summ + end do + + end subroutine dgaqd + !======================================================================= + + + !======================================================================= + subroutine dcpdp(nn,cz,cp,dcp) + ! + ! computes the fourier coefficients of the legendre + ! polynomial p_n^0 and its derivative. + ! n is the degree and n/2 or (n+1)/2 + ! coefficients are returned in cp depending on whether + ! n is even or odd. The same number of coefficients + ! are returned in dcp. For n even the constant + ! coefficient is returned in cz. + !===================================================================== + ! + ! Passed variables + !----------------- + integer, intent(in) :: nn + real(r8), intent(out) :: cz + real(r8), intent(out) :: cp(nn/2+1) + real(r8), intent(out) :: dcp(nn/2+1) + ! + ! Local Values + !-------------- + real(r8):: t1,t2,t3,t4 + integer :: ncp,jj + + ncp = (nn+1)/2 + t1 = -1._r8 + t2 = nn + 1._r8 + t3 = 0._r8 + t4 = nn + nn + 1._r8 + if(mod(nn,2).eq.0) then + cp(ncp) = 1._r8 + do jj = ncp,2,-1 + t1 = t1 + 2._r8 + t2 = t2 - 1._r8 + t3 = t3 + 1._r8 + t4 = t4 - 2._r8 + cp(jj-1) = (t1*t2)/ (t3*t4)*cp(jj) + end do + t1 = t1 + 2._r8 + t2 = t2 - 1._r8 + t3 = t3 + 1._r8 + t4 = t4 - 2._r8 + cz = (t1*t2)/ (t3*t4)*cp(1) + do jj = 1,ncp + dcp(jj) = (jj+jj)*cp(jj) + end do + else + cp(ncp) = 1._r8 + do jj = ncp - 1,1,-1 + t1 = t1 + 2._r8 + t2 = t2 - 1._r8 + t3 = t3 + 1._r8 + t4 = t4 - 2._r8 + cp(jj) = (t1*t2)/ (t3*t4)*cp(jj+1) + end do + do jj = 1,ncp + dcp(jj) = (jj+jj-1)*cp(jj) + end do + endif + + end subroutine dcpdp + !======================================================================= + + + !======================================================================= + subroutine dtpdp(nn,theta,cz,cp,dcp,pb,dpb) + ! + ! computes pn(theta) and its derivative dpb(theta) with + ! respect to theta + !===================================================================== + ! + ! Passed variables + !------------------ + integer, intent(in) :: nn + real(r8), intent(in) :: theta + real(r8), intent(in) :: cz + real(r8), intent(in) :: cp (nn/2+1) + real(r8), intent(in) :: dcp(nn/2+1) + real(r8), intent(out) :: pb + real(r8), intent(out) :: dpb + ! + ! Local Values + !-------------- + real(r8):: cdt,sdt,cth,sth,chh + integer :: kdo,kk + + cdt = dcos(theta+theta) + sdt = dsin(theta+theta) + if(mod(nn,2).eq.0) then + ! n even + !---------- + kdo = nn/2 + pb = .5_r8*cz + dpb = 0._r8 + if(nn.gt.0) then + cth = cdt + sth = sdt + do kk = 1,kdo +! pb = pb+cp(k)*cos(2*k*theta) +! dpb = dpb-(k+k)*cp(k)*sin(2*k*theta) + 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(k)*cos((2*k-1)*theta) +! dpb = dpb-(k+k-1)*cp(k)*sin((2*k-1)*theta) + 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 dtpdp + !======================================================================= + + + !======================================================================= + real(r8) function ddzeps(xx) + ! + ! estimate unit roundoff in quantities of size x. + ! + ! this program should function properly on all systems + ! satisfying the following two assumptions, + ! 1. the base used in representing floating point + ! numbers is not a power of three. + ! 2. the quantity a in statement 10 is represented to + ! the accuracy used in floating point variables + ! that are stored in memory. + ! the statement number 10 and the go to 10 are intended to + ! force optimizing compilers to generate code satisfying + ! assumption 2. + ! under these assumptions, it should be true that, + ! a is not exactly equal to four-thirds, + ! b has a zero for its last bit or digit, + ! c is not exactly equal to one, + ! eps measures the separation of 1.0 from + ! the next larger floating point number. + ! the developers of eispack would appreciate being informed + ! about any systems where these assumptions do not hold. + ! + ! this version dated 4/6/83. + ! + !===================================================================== + ! + ! Passed variables + !----------------- + real(r8), intent(in) :: xx + ! + ! Local Values + !-------------- + real(r8):: bb,cc,eps + real(r8), parameter :: aa = 4._r8/3._r8 + + eps = 0.0_r8 + + do while(eps == 0.0_r8) + bb = aa - 1._r8 + cc = bb + bb + bb + eps = abs(cc-1._r8) + end do + + ddzeps = eps*abs(xx) + + end function ddzeps + !======================================================================= + +end module zonal_mean_mod From f059e4e801f9cc9032884e1dabd12f44eb7570e8 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 17 Oct 2022 17:51:42 -0600 Subject: [PATCH 02/24] add conditionals for doing the TEM daigs modified: src/physics/cam/phys_grid_ctem.F90 --- src/physics/cam/phys_grid_ctem.F90 | 31 ++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index d07e7db9ee..056e1bb9bf 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -15,6 +15,7 @@ module phys_grid_ctem use cam_abortutils,only: endrun 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 implicit none @@ -30,6 +31,11 @@ module phys_grid_ctem integer :: nzalat = -huge(1) integer :: nzmbas = -huge(1) + integer, parameter :: nhours = 6 ! number of hours bewteen TEM calculations + integer :: ntimesteps = -huge(1) ! number of time steps bewteen TEM calculations + + logical :: do_tem_diags = .false. + contains !------------------------------------------------------------------------------ @@ -61,9 +67,12 @@ subroutine phys_grid_ctem_readnl(nlfile) call MPI_bcast(phys_grid_ctem_zm_nbas, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr) + do_tem_diags = phys_grid_ctem_zm_nbas>0 .and. phys_grid_ctem_za_nlat>0 + if (masterproc) then write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_zm_nbas: ',phys_grid_ctem_zm_nbas write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_za_nlat: ',phys_grid_ctem_za_nlat + write(iulog,*) 'phys_grid_ctem_readnl... do_tem_diags: ', do_tem_diags endif nzalat = phys_grid_ctem_za_nlat @@ -95,6 +104,8 @@ subroutine phys_grid_ctem_reg 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) @@ -154,6 +165,10 @@ end subroutine phys_grid_ctem_reg !----------------------------------------------------------------------------- subroutine phys_grid_ctem_init + real(r8) :: dtime + + if (.not.do_tem_diags) return + call addfld ('VTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') call addfld ('WTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') call addfld ('UVzaphys', (/'ilev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') @@ -170,6 +185,9 @@ subroutine phys_grid_ctem_init call addfld ('THzm3d', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) call addfld ('THphys', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) + dtime = get_step_size() + ntimesteps = nint( nhours*3600._r8/dtime ) ! number of steps per nhours + end subroutine phys_grid_ctem_init !----------------------------------------------------------------------------- @@ -216,6 +234,8 @@ subroutine phys_grid_ctem_diags(phys_state) real(r8), parameter :: hscale = 7000._r8 ! pressure scale height + if (.not.do_calc()) return + ui(:,:,:) = 0._r8 vi(:,:,:) = 0._r8 wi(:,:,:) = 0._r8 @@ -336,6 +356,17 @@ function zmean_fld( fld ) result(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 end module phys_grid_ctem From 7683d547f7764093dea1a0ed0185e722bb08c32d Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 19 Oct 2022 17:22:31 -0600 Subject: [PATCH 03/24] calculate scale height; add comments modified: src/physics/cam/phys_grid_ctem.F90 modified: src/physics/cam/physpkg.F90 --- src/physics/cam/phys_grid_ctem.F90 | 40 +++++++++++++++++++----------- src/physics/cam/physpkg.F90 | 1 + 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 056e1bb9bf..64aa8ace8d 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -17,6 +17,10 @@ module phys_grid_ctem 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 + implicit none private @@ -117,6 +121,7 @@ subroutine phys_grid_ctem_reg 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 @@ -126,19 +131,19 @@ subroutine phys_grid_ctem_reg 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('zmean_phys_fields_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 + ! 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) + 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)) @@ -156,8 +161,8 @@ subroutine phys_grid_ctem_reg 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.) + 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 @@ -232,7 +237,7 @@ subroutine phys_grid_ctem_diags(phys_state) real(r8) :: vthza(nzalat,pverp) real(r8) :: wthza(nzalat,pverp) - real(r8), parameter :: hscale = 7000._r8 ! pressure scale height + real(r8) :: sheight(pcols,pver) ! pressure scale height (m) if (.not.do_calc()) return @@ -254,9 +259,16 @@ subroutine phys_grid_ctem_diags(phys_state) 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,:) - w(:ncol,:,lchnk) = - hscale * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) + ! vertical velocity + w(:ncol,:,lchnk) = -sheight(:ncol,:) * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) + + ! interpolate to layer interfaces do k = 1,pverp call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%u(:,:), ui(:,k,lchnk) ) call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%v(:,:), vi(:,k,lchnk) ) @@ -266,13 +278,13 @@ subroutine phys_grid_ctem_diags(phys_state) end do - ! these need to be evaluated on the physics grid (3D) - ! to be used in the deviations calculation below + ! zonal means evaluated on the physics grid (3D) to be used in the deviations calculation below uzm(:,:,:) = zmean_fld(ui(:,:,:)) vzm(:,:,:) = zmean_fld(vi(:,:,:)) wzm(:,:,:) = zmean_fld(wi(:,:,:)) thzm(:,:,:) = zmean_fld(thi(:,:,:)) + ! diagnostic output do lchnk = begchunk, endchunk ncol = phys_state(lchnk)%ncol @@ -290,7 +302,6 @@ subroutine phys_grid_ctem_diags(phys_state) call outfld( 'Wzm3d', fld_tmp(:ncol,:), ncol, lchnk) end do - do lchnk = begchunk,endchunk ncol = phys_state(lchnk)%ncol do k = 1,pverp @@ -307,12 +318,13 @@ subroutine phys_grid_ctem_diags(phys_state) end do end do - ! evaluate and output these on the zonal-average grid + ! 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) + ! diagnostic output do j = 1,nzalat call outfld('UVzaphys',uvza(j,:),1,j) call outfld('UWzaphys',uwza(j,:),1,j) @@ -320,7 +332,7 @@ subroutine phys_grid_ctem_diags(phys_state) call outfld('WTHzaphys',wthza(j,:),1,j) end do - ! not needed --- only for sanity checks + ! the following is not needed --- only for sanity checks uv(:,:,:) = zmean_fld(uvp(:,:,:)) uw(:,:,:) = zmean_fld(uwp(:,:,:)) vth(:,:,:) = zmean_fld(vthp(:,:,:)) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 440e424a1d..545451447c 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -2957,6 +2957,7 @@ 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 From 706ede398b7c21ccaa7c8deec4cce09fd938511d Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 20 Oct 2022 09:32:47 -0600 Subject: [PATCH 04/24] fix bug in do_tem_diags setting modified: src/physics/cam/phys_grid_ctem.F90 --- src/physics/cam/phys_grid_ctem.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 64aa8ace8d..0c7e5204db 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -54,6 +54,9 @@ subroutine phys_grid_ctem_readnl(nlfile) namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat + phys_grid_ctem_zm_nbas = -1 + phys_grid_ctem_za_nlat = -1 + ! Read in namelist values !------------------------ if(masterproc) then @@ -79,8 +82,10 @@ subroutine phys_grid_ctem_readnl(nlfile) write(iulog,*) 'phys_grid_ctem_readnl... do_tem_diags: ', do_tem_diags endif - nzalat = phys_grid_ctem_za_nlat - nzmbas = phys_grid_ctem_zm_nbas + 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 From e9caadc912a46bfa502ab39e673793608c5cbede Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 1 Nov 2022 08:51:49 -0600 Subject: [PATCH 05/24] declare lower bound begchunk in assumed shape array args modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index a2acd602b7..84b3ff665b 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -478,7 +478,7 @@ subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp) ! Passed Variables !------------------ class(ZonalMean_t) :: this - real(r8),intent(in ):: I_Gdata(:,:,:) + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) real(r8),intent(out):: O_Bamp (:,:) ! ! Local Values @@ -513,7 +513,7 @@ subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp) ncols = get_ncols_p(lchnk) do cc = 1,ncols count=count+1 - Csum(count,nn) = I_Gdata(cc,ll,lchnk-begchunk+1)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) + Csum(count,nn) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk) end do end do end do @@ -586,7 +586,7 @@ subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata) !------------------ class(ZonalMean_t) :: this real(r8),intent(in ):: I_Bamp (:,:) - real(r8),intent(out):: O_Gdata(:,:,:) + real(r8),intent(out):: O_Gdata(:,:,begchunk:) ! ! Local Values !-------------- @@ -606,7 +606,7 @@ subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata) do lchnk=begchunk,endchunk ncols = get_ncols_p(lchnk) do cc = 1,ncols - O_Gdata(cc,ll,lchnk-begchunk+1) = O_Gdata(cc,ll,lchnk-begchunk+1) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn)) + 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 @@ -1171,7 +1171,7 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) ! Passed Variables !------------------ class(ZonalAverage_t):: this - real(r8),intent(in ):: I_Gdata(:,:,:) + real(r8),intent(in ):: I_Gdata(:,:,begchunk:) real(r8),intent(out):: O_Zdata(:,:) ! ! Local Values @@ -1206,7 +1206,7 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) jlat = this%idx_map(cc,lchnk) ns = jlat + (ll-1)*this%nlat count=count+1 - Asum(count,ns) = I_Gdata(cc,ll,lchnk-begchunk+1)*this%area_g(cc,lchnk) + Asum(count,ns) = I_Gdata(cc,ll,lchnk)*this%area_g(cc,lchnk) end do end do end do From 145da7b18aa8034d725a8d5a2f62e22a0eb71452 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 3 Nov 2022 13:59:07 -0600 Subject: [PATCH 06/24] update default IC for waccmsc on ne30pg3 grid; new testmods dir modified: bld/namelist_files/namelist_defaults_cam.xml new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam new file: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm --- bld/namelist_files/namelist_defaults_cam.xml | 1 + .../outfrq9s_physgrid_tem_1deg/shell_commands | 2 ++ .../outfrq9s_physgrid_tem_1deg/user_nl_cam | 9 +++++++ .../outfrq9s_physgrid_tem_1deg/user_nl_clm | 27 +++++++++++++++++++ 4 files changed, 39 insertions(+) create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/shell_commands create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index b29e1e5be4..b2fcae85ee 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -241,6 +241,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/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..6d14fecd26 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam @@ -0,0 +1,9 @@ +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_zm_nbas=120 +phys_grid_ctem_za_nlat=90 +fincl2 = 'VTHzm3d','WTHzm3d','UVzm3d','UWzm3d', + 'Uzm3d', 'Vzm3d', 'Wzm3d', 'THzm3d','THphys', +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 + From 14ae0a6447c5ff9c47fab60403dee9c310a41505 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 7 Nov 2022 13:47:37 -0700 Subject: [PATCH 07/24] test with mpas dycore new file: cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands new file: cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam new file: cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm modified: bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml --- .../use_cases/waccm_sc_2000_cam6.xml | 6 ++--- .../shell_commands | 2 ++ .../user_nl_cam | 24 +++++++++++++++++++ .../user_nl_clm | 4 ++++ 4 files changed, 33 insertions(+), 3 deletions(-) create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/shell_commands create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_clm 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/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..83fd733ed5 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam @@ -0,0 +1,24 @@ +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_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' From b2f6a1418488cbad7a90cccec6acff612f395836 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 8 Nov 2022 10:41:04 -0700 Subject: [PATCH 08/24] remove debugging diagnostics (*zm3d history fields) modified: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam modified: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam modified: src/physics/cam/phys_grid_ctem.F90 --- .../cam/outfrq9s_physgrid_tem/user_nl_cam | 2 - .../outfrq9s_physgrid_tem_1deg/user_nl_cam | 2 - src/physics/cam/phys_grid_ctem.F90 | 41 ------------------- 3 files changed, 45 deletions(-) diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam index 378a76e449..d7332869c5 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam @@ -4,6 +4,4 @@ nhtfrq=9,9,9,9,9,9,9,9,9 inithist='ENDOFRUN' phys_grid_ctem_zm_nbas=16 phys_grid_ctem_za_nlat=15 -fincl2 = 'VTHzm3d','WTHzm3d','UVzm3d','UWzm3d', - 'Uzm3d', 'Vzm3d', 'Wzm3d', 'THzm3d','THphys', fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' 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 index 6d14fecd26..e6f681f955 100644 --- 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 @@ -4,6 +4,4 @@ nhtfrq=9,9,9,9,9,9,9,9,9 inithist='ENDOFRUN' phys_grid_ctem_zm_nbas=120 phys_grid_ctem_za_nlat=90 -fincl2 = 'VTHzm3d','WTHzm3d','UVzm3d','UWzm3d', - 'Uzm3d', 'Vzm3d', 'Wzm3d', 'THzm3d','THphys', fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 0c7e5204db..940742e19f 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -184,15 +184,6 @@ subroutine phys_grid_ctem_init call addfld ('UVzaphys', (/'ilev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') call addfld ('UWzaphys', (/'ilev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('VTHzm3d',(/'ilev' /), 'A','MK/S', 'Meridional Heat Flux: 3D zon. mean', gridname='physgrid' ) - call addfld ('WTHzm3d',(/'ilev' /), 'A','MK/S', 'Vertical Heat Flux: 3D zon. mean', gridname='physgrid' ) - call addfld ('UVzm3d', (/'ilev' /), 'A','M2/S2','Meridional Flux of Zonal Momentum: 3D zon. mean', gridname='physgrid' ) - call addfld ('UWzm3d', (/'ilev' /), 'A','M2/S2','Vertical Flux of Zonal Momentum: 3D zon. mean', gridname='physgrid' ) - - call addfld ('Uzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean zonal wind - defined on ilev', gridname='physgrid') - call addfld ('Vzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean meridional wind - defined on ilev', gridname='physgrid' ) - call addfld ('Wzm3d', (/'ilev' /), 'A','M/S', 'Zonal-Mean vertical wind - defined on ilev', gridname='physgrid' ) - call addfld ('THzm3d', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) call addfld ('THphys', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) dtime = get_step_size() @@ -223,11 +214,6 @@ subroutine phys_grid_ctem_diags(phys_state) real(r8) :: vthp(pcols,pverp,begchunk:endchunk) real(r8) :: wthp(pcols,pverp,begchunk:endchunk) - real(r8) :: uv(pcols,pverp,begchunk:endchunk) - real(r8) :: uw(pcols,pverp,begchunk:endchunk) - real(r8) :: vth(pcols,pverp,begchunk:endchunk) - real(r8) :: wth(pcols,pverp,begchunk:endchunk) - integer :: lchnk, ncol, j, k real(r8) :: fld_tmp(pcols,pverp) @@ -296,15 +282,6 @@ subroutine phys_grid_ctem_diags(phys_state) fld_tmp(:ncol,:) = thi(:ncol,:,lchnk) call outfld( 'THphys', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = thzm(:ncol,:,lchnk) - call outfld( 'THzm3d', fld_tmp(:ncol,:), ncol, lchnk) - - fld_tmp(:ncol,:) = uzm(:ncol,:,lchnk) - call outfld( 'Uzm3d', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = vzm(:ncol,:,lchnk) - call outfld( 'Vzm3d', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = wzm(:ncol,:,lchnk) - call outfld( 'Wzm3d', fld_tmp(:ncol,:), ncol, lchnk) end do do lchnk = begchunk,endchunk @@ -337,24 +314,6 @@ subroutine phys_grid_ctem_diags(phys_state) call outfld('WTHzaphys',wthza(j,:),1,j) end do - ! the following is not needed --- only for sanity checks - uv(:,:,:) = zmean_fld(uvp(:,:,:)) - uw(:,:,:) = zmean_fld(uwp(:,:,:)) - vth(:,:,:) = zmean_fld(vthp(:,:,:)) - wth(:,:,:) = zmean_fld(wthp(:,:,:)) - - do lchnk = begchunk, endchunk - ncol = phys_state(lchnk)%ncol - fld_tmp(:ncol,:) = uv(:ncol,:,lchnk) - call outfld( 'UVzm3d', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = uw(:ncol,:,lchnk) - call outfld( 'UWzm3d', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = vth(:ncol,:,lchnk) - call outfld( 'VTHzm3d', fld_tmp(:ncol,:), ncol, lchnk) - fld_tmp(:ncol,:) = wth(:ncol,:,lchnk) - call outfld( 'WTHzm3d', fld_tmp(:ncol,:), ncol, lchnk) - end do - contains !------------------------------------------------------------------------------ From 185d688ba66637275ef4887230aae06abc0258df Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 9 Nov 2022 16:45:41 -0700 Subject: [PATCH 09/24] Implement nfreq namelist variable; misc changes requested by code reviewer modified: bld/namelist_files/namelist_definition.xml modified: cime_config/testdefs/testmods_dirs/cam/outfrq1d_physgrid_tem_mpasa120_wcmsc/user_nl_cam modified: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam modified: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_cam modified: src/physics/cam/phys_grid_ctem.F90 --- bld/namelist_files/namelist_definition.xml | 7 +++ .../user_nl_cam | 1 + .../cam/outfrq9s_physgrid_tem/user_nl_cam | 1 + .../outfrq9s_physgrid_tem_1deg/user_nl_cam | 1 + src/physics/cam/phys_grid_ctem.F90 | 49 +++++++++++++------ 5 files changed, 44 insertions(+), 15 deletions(-) diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index ef9e7afe31..be4fb61972 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -433,6 +433,13 @@ 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/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 index 83fd733ed5..1212d35edd 100644 --- 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 @@ -10,6 +10,7 @@ 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 diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam index d7332869c5..2797dc0c3a 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam @@ -2,6 +2,7 @@ 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=16 phys_grid_ctem_za_nlat=15 fincl3 = 'VTHzaphys','WTHzaphys','UVzaphys','UWzaphys' 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 index e6f681f955..9cf4e0a97d 100644 --- 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 @@ -2,6 +2,7 @@ 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/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 940742e19f..e7fb045686 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -35,7 +35,6 @@ module phys_grid_ctem integer :: nzalat = -huge(1) integer :: nzmbas = -huge(1) - integer, parameter :: nhours = 6 ! number of hours bewteen TEM calculations integer :: ntimesteps = -huge(1) ! number of time steps bewteen TEM calculations logical :: do_tem_diags = .false. @@ -49,13 +48,17 @@ subroutine phys_grid_ctem_readnl(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 + namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat, phys_grid_ctem_nfreq - phys_grid_ctem_zm_nbas = -1 - phys_grid_ctem_za_nlat = -1 + phys_grid_ctem_zm_nbas = 0 + phys_grid_ctem_za_nlat = 0 + phys_grid_ctem_nfreq = 0 ! Read in namelist values !------------------------ @@ -73,13 +76,34 @@ subroutine phys_grid_ctem_readnl(nlfile) call MPI_bcast(phys_grid_ctem_zm_nbas, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr) + call MPI_bcast(phys_grid_ctem_nfreq, 1, mpi_integer, masterprocid, mpicom, ierr) - do_tem_diags = phys_grid_ctem_zm_nbas>0 .and. phys_grid_ctem_za_nlat>0 + 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') + 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') + end if + do_tem_diags = .true. + end if if (masterproc) then - write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_zm_nbas: ',phys_grid_ctem_zm_nbas - write(iulog,*) 'phys_grid_ctem_readnl... phys_grid_ctem_za_nlat: ',phys_grid_ctem_za_nlat - write(iulog,*) 'phys_grid_ctem_readnl... do_tem_diags: ', do_tem_diags + 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 @@ -175,20 +199,14 @@ end subroutine phys_grid_ctem_reg !----------------------------------------------------------------------------- subroutine phys_grid_ctem_init - real(r8) :: dtime - if (.not.do_tem_diags) return call addfld ('VTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') call addfld ('WTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') call addfld ('UVzaphys', (/'ilev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') call addfld ('UWzaphys', (/'ilev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('THphys', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) - dtime = get_step_size() - ntimesteps = nint( nhours*3600._r8/dtime ) ! number of steps per nhours - end subroutine phys_grid_ctem_init !----------------------------------------------------------------------------- @@ -217,7 +235,8 @@ subroutine phys_grid_ctem_diags(phys_state) integer :: lchnk, ncol, j, k real(r8) :: fld_tmp(pcols,pverp) - real(r8) :: theta(pcols,pver,begchunk:endchunk) ! potential temperature + ! potential temperature + real(r8) :: theta(pcols,pver,begchunk:endchunk) real(r8) :: thi(pcols,pverp,begchunk:endchunk) real(r8) :: thzm(pcols,pverp,begchunk:endchunk) From 5bbcc0e1feb1691599b0bdb685f549840a5768b1 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 9 Nov 2022 19:03:48 -0700 Subject: [PATCH 10/24] misc changes requested by code reviewer (round 2) modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 102 +++++++++++++++++------------------ 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 84b3ff665b..190e177bf2 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -6,14 +6,14 @@ module zonal_mean_mod ! 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/synthsis of zonal mean values +! 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/synthsis of zonal mean values +! 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 gridpoints +! area weighted bin-averaging of 2D grid points ! assigned to each latitude band. ! ! NOTE: The weighting of the Zonal Profiles values is scaled such @@ -747,13 +747,13 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! (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)) + 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 - end do ! Invert to get the basis amplitude map !-------------------------------------- @@ -830,7 +830,7 @@ subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp) ! Local Values !-------------- real(r8),allocatable:: Gcov(:,:) - integer:: ii,nn,n2,ll + integer:: ii,nn,n2,ilev integer :: nlev @@ -839,22 +839,22 @@ subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp) ! Compute Covariance with input data and basis functions !-------------------------------------------------------- allocate(Gcov(this%nbas,nlev)) - do ll=1,nlev + do ilev=1,nlev do nn=1,this%nbas - Gcov(nn,ll) = 0._r8 + Gcov(nn,ilev) = 0._r8 do ii=1,this%nlat - Gcov(nn,ll) = Gcov(nn,ll) + (I_Zdata(ii,ll)*this%basis(ii,nn)*this%area(ii)) + 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 ll=1,nlev + do ilev=1,nlev do nn=1,this%nbas - O_Bamp(nn,ll) = 0._r8 + O_Bamp(nn,ilev) = 0._r8 do n2=1,this%nbas - O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(n2,ll) + O_Bamp(nn,ilev) = O_Bamp(nn,ilev) + this%map(n2,nn)*Gcov(n2,ilev) end do end do end do @@ -885,9 +885,9 @@ subroutine eval_ZonalProfile_1Dgrid(this,I_Bamp,O_Zdata) !-------------------------------------------------- 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 + 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 @@ -909,7 +909,7 @@ subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata) ! ! Local Values !-------------- - integer:: ii,nn,ll + integer:: ii,nn,ilev integer :: nlev @@ -919,11 +919,11 @@ subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata) !-------------------------------------------------- O_Zdata(1:this%nlat,1:nlev) = 0._r8 do nn=1,this%nbas - do ll=1,nlev - do ii=1,this%nlat - O_Zdata(ii,ll) = O_Zdata(ii,ll) + (I_Bamp(nn,ll)*this%basis(ii,nn)) - end do - end do + 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 @@ -1037,7 +1037,7 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) end do end do - ! Set boundaroes for Latitude bins + ! Set boundaries for Latitude bins !----------------------------------- BinLat(1) = 0._r8 BinLat(this%nlat+1) = 4._r8*atan(1._r8) @@ -1179,7 +1179,7 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) real(r8),allocatable:: Gsum(:) real(r8),allocatable:: Asum(:,:) integer:: nn,ncols,lchnk,cc,jlat - integer:: Nsum,ll,ns + integer:: Nsum,ilev,ns integer :: nlev integer :: nlcols, count @@ -1198,15 +1198,15 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) ! Compute area-weighted sums !----------------------------- - do ll = 1,nlev + 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 + (ll-1)*this%nlat + ns = jlat + (ilev-1)*this%nlat count=count+1 - Asum(count,ns) = I_Gdata(cc,ll,lchnk)*this%area_g(cc,lchnk) + Asum(count,ns) = I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk) end do end do end do @@ -1215,10 +1215,10 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) ! Divide by area norm to get the averages !----------------------------------------- - do ll = 1,nlev + do ilev = 1,nlev do nn = 1,this%nlat - ns = nn + (ll-1)*this%nlat - O_Zdata(nn,ll) = Gsum(ns)/this%a_norm(nn) + ns = nn + (ilev-1)*this%nlat + O_Zdata(nn,ilev) = Gsum(ns)/this%a_norm(nn) end do end do @@ -1330,7 +1330,7 @@ subroutine dalfk(nn,mm,cp) real(r8):: fk real(r8):: a1,b1,C1 integer :: ma,nmms2,nex - integer :: ii,ll + integer :: ii,jj real(r8),parameter:: SC10=1024._r8 real(r8),parameter:: SC20=SC10*SC10 @@ -1398,29 +1398,29 @@ subroutine dalfk(nn,mm,cp) fnmsq = fnnp1 - 2._r8*ma*ma if((mod(nn,2).eq.0).and.(mod(ma,2).eq.0)) then - ll = 1+(nn+1)/2 + jj = 1+(nn+1)/2 else - ll = (nn+1)/2 + jj = (nn+1)/2 endif - cp(ll) = cp2 + cp(jj) = cp2 if(mm.lt.0) then - if(mod(ma,2).ne.0) cp(ll) = -cp(ll) + if(mod(ma,2).ne.0) cp(jj) = -cp(jj) endif - if(ll.le.1) return + if(jj.le.1) return fk = nn a1 = (fk-2._r8)*(fk-1._r8) - fnnp1 b1 = 2._r8* (fk*fk-fnmsq) - cp(ll-1) = b1*cp(ll)/a1 + cp(jj-1) = b1*cp(jj)/a1 30 continue - ll = ll - 1 - if(ll.le.1) return + jj = jj - 1 + if(jj.le.1) return 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(ll-1) = -(b1*cp(ll)+c1*cp(ll+1))/a1 + cp(jj-1) = -(b1*cp(jj)+c1*cp(jj+1))/a1 goto 30 end subroutine dalfk @@ -1616,7 +1616,7 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) real(r8),allocatable:: Mwrk(:,:),Rscl(:) integer ,allocatable:: Indx(:) real(r8):: Psgn,Mmax,Mval,Sval - integer :: ii,jj,kk,ll,i2,ii_max + integer :: ii,jj,kk,ndx,i2,ii_max ! Allocate work space !--------------------- @@ -1710,9 +1710,9 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) i2 = 0 do ii=11,Nbas - ll = Indx(ii) - Sval = O_InvMat(ll,kk) - O_InvMat(ll,kk) = O_InvMat(ii,kk) + ndx = Indx(ii) + Sval = O_InvMat(ndx,kk) + O_InvMat(ndx,kk) = O_InvMat(ii,kk) if(i2.ne.0) then do jj=i2,(ii-1) Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk) @@ -2048,22 +2048,22 @@ end subroutine dtpdp !======================================================================= real(r8) function ddzeps(xx) ! - ! estimate unit roundoff in quantities of size x. + ! estimate unit roundoff in quantities of size xx. ! ! this program should function properly on all systems ! satisfying the following two assumptions, ! 1. the base used in representing floating point ! numbers is not a power of three. - ! 2. the quantity a in statement 10 is represented to + ! 2. the quantity aa in statement 10 is represented to ! the accuracy used in floating point variables ! that are stored in memory. ! the statement number 10 and the go to 10 are intended to ! force optimizing compilers to generate code satisfying ! assumption 2. ! under these assumptions, it should be true that, - ! a is not exactly equal to four-thirds, - ! b has a zero for its last bit or digit, - ! c is not exactly equal to one, + ! aa is not exactly equal to four-thirds, + ! bb has a zero for its last bit or digit, + ! cc is not exactly equal to one, ! eps measures the separation of 1.0 from ! the next larger floating point number. ! the developers of eispack would appreciate being informed From f45f732af1f687f8664450f8d9deab550574c403 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 10 Nov 2022 08:20:53 -0700 Subject: [PATCH 11/24] implement final routines; mics changes modified: src/physics/cam/phys_grid_ctem.F90 modified: src/physics/cam/physpkg.F90 modified: src/utils/zonal_mean_mod.F90 --- src/physics/cam/phys_grid_ctem.F90 | 8 +++++ src/physics/cam/physpkg.F90 | 2 ++ src/utils/zonal_mean_mod.F90 | 54 +++++++++++++++++++++++------- 3 files changed, 52 insertions(+), 12 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index e7fb045686..e10d013919 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -28,6 +28,7 @@ module phys_grid_ctem 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 @@ -364,4 +365,11 @@ 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 545451447c..63b2a02728 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1305,6 +1305,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 !----------------------------------------------------------------------- ! @@ -1327,6 +1328,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 diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 190e177bf2..e4d564af68 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -187,7 +187,8 @@ module zonal_mean_mod procedure,private,pass:: calc_ZonalMean_3Damps procedure,private,pass:: eval_ZonalMean_2Dgrid procedure,private,pass:: eval_ZonalMean_3Dgrid - end type ZonalMean_t + procedure, pass :: final => final_ZonalMean + end type ZonalMean_t type ZonalProfile_t private @@ -206,6 +207,7 @@ module zonal_mean_mod 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 @@ -221,6 +223,7 @@ module zonal_mean_mod 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 :: PI2 = 2._r8*atan(1._r8) ! pi/2 @@ -615,6 +618,16 @@ subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata) 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) @@ -929,6 +942,16 @@ subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata) 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) @@ -1111,7 +1134,7 @@ end subroutine init_ZonalAverage !======================================================================= subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata) ! - ! calc_ZonalProfile_2DbinAvg: Given 2D data values for ncol gridpoints, + ! calc_ZonalAverage_2DbinAvg: Given 2D data values for ncol gridpoints, ! compute the nlat area weighted binAvg profile !===================================================================== ! @@ -1157,6 +1180,8 @@ subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata) O_Zdata(nn) = O_Zdata(nn)/this%a_norm(nn) end do + deallocate(Asum) + end subroutine calc_ZonalAverage_2DbinAvg !======================================================================= @@ -1164,7 +1189,7 @@ end subroutine calc_ZonalAverage_2DbinAvg !======================================================================= subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) ! - ! calc_ZonalProfile_3DbinAvg: Given 3D data values for ncol,nlev gridpoints, + ! calc_ZonalAverage_3DbinAvg: Given 3D data values for ncol,nlev gridpoints, ! compute the nlat,nlev area weighted binAvg profile !===================================================================== ! @@ -1222,9 +1247,23 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) 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 dalfk(nn,mm,cp) @@ -2013,8 +2052,6 @@ subroutine dtpdp(nn,theta,cz,cp,dcp,pb,dpb) cth = cdt sth = sdt do kk = 1,kdo -! pb = pb+cp(k)*cos(2*k*theta) -! dpb = dpb-(k+k)*cp(k)*sin(2*k*theta) pb = pb + cp(kk)*cth dpb = dpb - dcp(kk)*sth chh = cdt*cth - sdt*sth @@ -2031,8 +2068,6 @@ subroutine dtpdp(nn,theta,cz,cp,dcp,pb,dpb) cth = dcos(theta) sth = dsin(theta) do kk = 1,kdo -! pb = pb+cp(k)*cos((2*k-1)*theta) -! dpb = dpb-(k+k-1)*cp(k)*sin((2*k-1)*theta) pb = pb + cp(kk)*cth dpb = dpb - dcp(kk)*sth chh = cdt*cth - sdt*sth @@ -2057,9 +2092,6 @@ real(r8) function ddzeps(xx) ! 2. the quantity aa in statement 10 is represented to ! the accuracy used in floating point variables ! that are stored in memory. - ! the statement number 10 and the go to 10 are intended to - ! force optimizing compilers to generate code satisfying - ! assumption 2. ! under these assumptions, it should be true that, ! aa is not exactly equal to four-thirds, ! bb has a zero for its last bit or digit, @@ -2069,8 +2101,6 @@ real(r8) function ddzeps(xx) ! the developers of eispack would appreciate being informed ! about any systems where these assumptions do not hold. ! - ! this version dated 4/6/83. - ! !===================================================================== ! ! Passed variables From 5ef0f8b88e85347c605d5ad5976e3a584cf51b5b Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 10 Nov 2022 08:41:05 -0700 Subject: [PATCH 12/24] tweaks to some of the comments modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index e4d564af68..3c57f55c45 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -232,8 +232,8 @@ module zonal_mean_mod !======================================================================= subroutine init_ZonalMean(this,I_nbas) ! - ! init_ZonalMean: Initialize the ZonalMean datastruture for the - ! given physgrid gridpoints. It is assumed that the domain + ! init_ZonalMean: Initialize the ZonalMean datastrutures for the + ! physics grid. It is assumed that the domain ! of these gridpoints spans the surface of the sphere. ! The representation of basis functions functions is ! normalized w.r.t integration over the sphere. @@ -689,8 +689,8 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) allocate(Bcov (I_nbas,I_nbas)) ! Optionally create the Latitude Gridpoints - ! and their associated area weights. Otherwise it - ! is assumed that the user is supplying them. + ! and their associated area weights. Otherwise + ! they need to be supplied by the user. !----------------------------------------------- if(generate_lats) then @@ -1016,8 +1016,8 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) allocate(Anorm (I_nlat)) ! Optionally create the Latitude Gridpoints - ! and their associated area weights. Otherwise it - ! is assumed that the user is supplying them. + ! and their associated area weights. Otherwise + ! they need to be supplied by the user. !----------------------------------------------- if(generate_lats) then From 82e8203f3c9a6557ef035af8a3b0209157dab723 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Fri, 18 Nov 2022 14:46:24 -0700 Subject: [PATCH 13/24] check allocation status modified: src/physics/cam/phys_grid_ctem.F90 modified: src/utils/zonal_mean_mod.F90 --- src/physics/cam/phys_grid_ctem.F90 | 10 +- src/utils/zonal_mean_mod.F90 | 229 +++++++++++++++++------------ 2 files changed, 146 insertions(+), 93 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index e10d013919..df7a75822f 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -12,7 +12,7 @@ module phys_grid_ctem use zonal_mean_mod,only: ZonalAverage_t, ZonalMean_t use physconst, only: pi use cam_logfile, only: iulog - use cam_abortutils,only: endrun + 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 @@ -76,8 +76,11 @@ subroutine phys_grid_ctem_readnl(nlfile) 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 @@ -130,7 +133,7 @@ subroutine phys_grid_ctem_reg real(r8) :: dlatrad, dlatdeg, lat1, lat2 real(r8) :: total_area real(r8) :: total_wght - integer :: j + integer :: j, astat real(r8), parameter :: latdeg0 = -90._r8 real(r8), parameter :: latrad0 = -pi*0.5_r8 @@ -176,7 +179,8 @@ subroutine phys_grid_ctem_reg zalon_coord => horiz_coord_create('zalon', '', 1, 'longitude', 'degrees_east', 1, 1, zalons) ! grid decomposition map - allocate(grid_map(4,nzalat)) + 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 diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 3c57f55c45..67a5976261 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -157,7 +157,7 @@ module zonal_mean_mod 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 + use cam_abortutils, only: endrun, handle_allocate_error use spmd_utils, only: mpicom use phys_grid, only: ngcols_p => num_global_phys_cols @@ -258,7 +258,8 @@ subroutine init_ZonalMean(this,I_nbas) integer :: nn,n2,nb,lchnk,ncols,cc integer :: cnum,Cvec_len - integer :: nlcols, count + 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') @@ -271,29 +272,39 @@ subroutine init_ZonalMean(this,I_nbas) if(allocated(this%map )) deallocate(this%map) this%nbas = I_nbas - allocate(this%area (pcols,begchunk:endchunk)) - allocate(this%basis(pcols,begchunk:endchunk,I_nbas)) - allocate(this%map (I_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 + do n2=nn,this%nbas + Cvec_len = Cvec_len + 1 + end do end do nlcols = get_nlcols_p() - allocate(Clats(pcols,begchunk:endchunk)) - allocate(Bcoef(I_nbas)) - allocate(Csum (nlcols, Cvec_len)) - allocate(Cvec (Cvec_len)) - allocate(Bsum (nlcols, I_nbas)) - allocate(Bnorm(I_nbas)) - allocate(Bcov (I_nbas,I_nbas)) + allocate(Clats(pcols,begchunk:endchunk), stat=astat) + call handle_allocate_error(astat, subname, 'Clats') + allocate(Bcoef(I_nbas), 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 @@ -302,54 +313,54 @@ subroutine init_ZonalMean(this,I_nbas) ! 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 + PI2 - end do + 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 + PI2 + end do end do ! Add first basis for the mean values. !------------------------------------------ do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - this%basis(cc,lchnk,1) = 1._r8/sqrt(8._r8*PI2) - end do + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + this%basis(cc,lchnk,1) = 1._r8/sqrt(8._r8*PI2) + end do end do ! Loop over the remaining basis functions !--------------------------------------- do nn=2,this%nbas - nb = nn-1 + nb = nn-1 - ! Generate coefs for the basis - !------------------------------ - call dalfk(nb,0,Bcoef) + ! Generate coefs for the basis + !------------------------------ + call dalfk(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 dlfpt(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn)) - end do - end do - end do ! nn=1,this%nbas + ! Create basis for the coefs at each ncol gridpoint + !--------------------------------------------------- + do lchnk=begchunk,endchunk + ncols = get_ncols_p(lchnk) + do cc = 1,ncols + call dlfpt(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 + 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) @@ -387,11 +398,11 @@ subroutine init_ZonalMean(this,I_nbas) 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 + 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 @@ -430,12 +441,16 @@ subroutine calc_ZonalMean_2Damps(this,I_Gdata,O_Bamp) real(r8),allocatable :: Csum(:,:) real(r8),allocatable :: Gcov(:) integer :: nn,n2,ncols,lchnk,cc - integer :: nlcols, count + integer :: nlcols, count, astat + + character(len=*), parameter :: subname = 'calc_ZonalMean_2Damps' nlcols = get_nlcols_p() - allocate(Gcov(this%nbas)) - allocate(Csum(nlcols, this%nbas)) + 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 @@ -490,15 +505,18 @@ subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp) real(r8),allocatable:: Gcov (:) integer:: nn,n2,ncols,lchnk,cc integer:: Nsum,ns,ll - integer :: nlcols, count + 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)) - allocate(Csum(nlcols, this%nbas)) + 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 @@ -663,9 +681,11 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) real(r8),allocatable:: Bcoef(:) real(r8),allocatable:: Bcov (:,:) real(r8):: Bnorm - integer :: ii,nn,n2,nb,ierr + integer :: ii,nn,n2,nb,ierr, astat logical :: generate_lats + character(len=*), parameter :: subname = 'init_ZonalProfile' + generate_lats = .false. if (present(GEN_GAUSSLATS)) then @@ -680,13 +700,19 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) this%nlat = I_nlat this%nbas = I_nbas - allocate(this%area (I_nlat)) - allocate(this%basis(I_nlat,I_nbas)) - allocate(this%map (I_nbas,I_nbas)) - - allocate(Clats(I_nlat)) - allocate(Bcoef(I_nbas)) - allocate(Bcov (I_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), 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 @@ -799,11 +825,13 @@ subroutine calc_ZonalProfile_1Damps(this,I_Zdata,O_Bamp) ! Local Values !-------------- real(r8),allocatable:: Gcov(:) - integer:: ii,nn,n2 + integer:: ii,nn,n2, astat + character(len=*), parameter :: subname = 'calc_ZonalProfile_1Damps' ! Compute Covariance with input data and basis functions !-------------------------------------------------------- - allocate(Gcov(this%nbas)) + 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 @@ -845,13 +873,15 @@ subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp) real(r8),allocatable:: Gcov(:,:) integer:: ii,nn,n2,ilev - integer :: nlev + 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)) + 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 @@ -983,10 +1013,11 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) real(r8),allocatable:: Asum (:,:) real(r8),allocatable:: Anorm (:) real(r8):: area(pcols),rlat - integer :: nn,jj,ierr + 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. @@ -1004,16 +1035,25 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) if(allocated(this%idx_map)) deallocate(this%idx_map) this%nlat = I_nlat - allocate(this%area (I_nlat)) - allocate(this%a_norm (I_nlat)) - allocate(this%area_g (pcols,begchunk:endchunk)) - allocate(this%idx_map(pcols,begchunk:endchunk)) - - allocate(Clats (I_nlat)) - allocate(BinLat(I_nlat+1)) - allocate(Glats (pcols,begchunk:endchunk)) - allocate(Asum (nlcols,I_nlat)) - allocate(Anorm (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 @@ -1148,14 +1188,16 @@ subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata) !-------------- real(r8),allocatable:: Asum (:,:) integer:: nn,ncols,lchnk,cc,jlat - integer :: nlcols, count + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'calc_ZonalAverage_2DbinAvg' nlcols = get_nlcols_p() ! Initialize Zonal profile !--------------------------- - allocate(Asum(nlcols,this%nlat)) + allocate(Asum(nlcols,this%nlat), stat=astat) + call handle_allocate_error(astat, subname, 'Asum') Asum(:,:) = 0._r8 O_Zdata(1:this%nlat) = 0._r8 @@ -1207,7 +1249,8 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) integer:: Nsum,ilev,ns integer :: nlev - integer :: nlcols, count + integer :: nlcols, count, astat + character(len=*), parameter :: subname = 'calc_ZonalAverage_3DbinAvg' nlev = size(I_Gdata,dim=2) nlcols = get_nlcols_p() @@ -1215,8 +1258,10 @@ subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata) ! Initialize Zonal profile !--------------------------- Nsum = this%nlat*nlev - allocate(Gsum(Nsum)) - allocate(Asum(nlcols,Nsum)) + 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 @@ -1655,13 +1700,17 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) real(r8),allocatable:: Mwrk(:,:),Rscl(:) integer ,allocatable:: Indx(:) real(r8):: Psgn,Mmax,Mval,Sval - integer :: ii,jj,kk,ndx,i2,ii_max + integer :: ii,jj,kk,ndx,i2,ii_max, astat + character(len=*), parameter :: subname = 'Invert_Matrix' ! Allocate work space !--------------------- - allocate(Mwrk(Nbas,Nbas)) - allocate(Rscl(Nbas)) - allocate(Indx(Nbas)) + 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 !------------------------------------------------- From a01b6527a26225b84c74ffc7e7a9222ad1a05e03 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 21 Nov 2022 08:56:57 -0700 Subject: [PATCH 14/24] use physconst pi modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 67a5976261..e6310b330e 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -159,8 +159,8 @@ module zonal_mean_mod use shr_reprosum_mod,only: shr_reprosum_calc use cam_abortutils, only: endrun, handle_allocate_error use spmd_utils, only: mpicom - - use phys_grid, only: ngcols_p => num_global_phys_cols + use physconst, only: pi + use phys_grid, only: ngcols_p => num_global_phys_cols implicit none private @@ -226,7 +226,7 @@ module zonal_mean_mod procedure, pass :: final => final_ZonalAverage end type ZonalAverage_t - real(r8), parameter :: PI2 = 2._r8*atan(1._r8) ! pi/2 + real(r8), parameter :: PI2 = 0.5_r8*pi contains !======================================================================= From 3d9d37d8db6dc16a07a8043b70f23ff50ab26e53 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 21 Dec 2022 15:23:43 -0700 Subject: [PATCH 15/24] Remove vertical interpolation onto layer interfaces modified: src/physics/cam/phys_grid_ctem.F90 --- src/physics/cam/phys_grid_ctem.F90 | 120 ++++++++++++----------------- 1 file changed, 49 insertions(+), 71 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index df7a75822f..5303851c65 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -4,9 +4,7 @@ !---------------------------------------------------------------------------------- module phys_grid_ctem use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: begchunk, endchunk, pcols, pver, pverp - use ref_pres, only: pref_edge - use interpolate_data, only: vertinterp + 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 @@ -206,11 +204,11 @@ subroutine phys_grid_ctem_init if (.not.do_tem_diags) return - call addfld ('VTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') - call addfld ('WTHzaphys',(/'ilev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') - call addfld ('UVzaphys', (/'ilev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('UWzaphys', (/'ilev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('THphys', (/'ilev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) + call addfld ('VTHzaphys',(/'lev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('WTHzaphys',(/'lev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') + call addfld ('UVzaphys', (/'lev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') + call addfld ('UWzaphys', (/'lev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') + call addfld ('THphys', (/'lev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) end subroutine phys_grid_ctem_init @@ -219,57 +217,41 @@ end subroutine phys_grid_ctem_init subroutine phys_grid_ctem_diags(phys_state) type(physics_state), intent(in) :: phys_state(begchunk:endchunk) - real(r8) :: ui(pcols,pverp,begchunk:endchunk) - real(r8) :: vi(pcols,pverp,begchunk:endchunk) - real(r8) :: wi(pcols,pverp,begchunk:endchunk) + character(len=*), parameter :: prefix = 'phys_grid_ctem_diags: ' - real(r8) :: uzm(pcols,pverp,begchunk:endchunk) - real(r8) :: vzm(pcols,pverp,begchunk:endchunk) - real(r8) :: wzm(pcols,pverp,begchunk:endchunk) + 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,pverp,begchunk:endchunk) - real(r8) :: vd(pcols,pverp,begchunk:endchunk) - real(r8) :: wd(pcols,pverp,begchunk:endchunk) - real(r8) :: thd(pcols,pverp,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,pverp,begchunk:endchunk) - real(r8) :: uwp(pcols,pverp,begchunk:endchunk) - real(r8) :: vthp(pcols,pverp,begchunk:endchunk) - real(r8) :: wthp(pcols,pverp,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 - real(r8) :: fld_tmp(pcols,pverp) ! potential temperature real(r8) :: theta(pcols,pver,begchunk:endchunk) - real(r8) :: thi(pcols,pverp,begchunk:endchunk) - real(r8) :: thzm(pcols,pverp,begchunk:endchunk) - - real(r8) :: w(pcols,pver,begchunk:endchunk) + real(r8) :: thzm(pcols,pver,begchunk:endchunk) - real(r8) :: uvza(nzalat,pverp) - real(r8) :: uwza(nzalat,pverp) - real(r8) :: vthza(nzalat,pverp) - real(r8) :: wthza(nzalat,pverp) + 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 - ui(:,:,:) = 0._r8 - vi(:,:,:) = 0._r8 - wi(:,:,:) = 0._r8 - thi(:,:,:) = 0._r8 - - uzm(:,:,:) = 0._r8 - vzm(:,:,:) = 0._r8 - wzm(:,:,:) = 0._r8 - thzm(:,:,:) = 0._r8 - - ud(:,:,:) = 0._r8 - vd(:,:,:) = 0._r8 - uvp(:,:,:) = 0._r8 - do lchnk = begchunk,endchunk ncol = phys_state(lchnk)%ncol @@ -281,41 +263,32 @@ subroutine phys_grid_ctem_diags(phys_state) 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,:) - - ! interpolate to layer interfaces - do k = 1,pverp - call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%u(:,:), ui(:,k,lchnk) ) - call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), phys_state(lchnk)%v(:,:), vi(:,k,lchnk) ) - call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), theta(:,:,lchnk), thi(:,k,lchnk) ) - call vertinterp( ncol, pcols, pver, phys_state(lchnk)%pmid(:,:), pref_edge(k), w(:,:,lchnk), wi(:,k,lchnk) ) - end do + w(:ncol,:,lchnk) = -sheight(:ncol,:) * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) + + u(:ncol,:,lchnk) = phys_state(lchnk)%u(:,:) + v(:ncol,:,lchnk) = phys_state(lchnk)%v(:,:) end do ! zonal means evaluated on the physics grid (3D) to be used in the deviations calculation below - uzm(:,:,:) = zmean_fld(ui(:,:,:)) - vzm(:,:,:) = zmean_fld(vi(:,:,:)) - wzm(:,:,:) = zmean_fld(wi(:,:,:)) - thzm(:,:,:) = zmean_fld(thi(:,:,:)) + uzm(:,:,:) = zmean_fld(u(:,:,:)) + vzm(:,:,:) = zmean_fld(v(:,:,:)) + wzm(:,:,:) = zmean_fld(w(:,:,:)) + thzm(:,:,:) = zmean_fld(theta(:,:,:)) ! diagnostic output do lchnk = begchunk, endchunk - ncol = phys_state(lchnk)%ncol - - fld_tmp(:ncol,:) = thi(:ncol,:,lchnk) - call outfld( 'THphys', fld_tmp(:ncol,:), ncol, lchnk) - + call outfld( 'THphys', theta(:,:,lchnk), pcols, lchnk) end do do lchnk = begchunk,endchunk ncol = phys_state(lchnk)%ncol - do k = 1,pverp + do k = 1,pver ! zonal deviations - thd(:ncol,k,lchnk) = thi(:ncol,k,lchnk) - thzm(:ncol,k,lchnk) - ud(:ncol,k,lchnk) = ui(:ncol,k,lchnk) - uzm(:ncol,k,lchnk) - vd(:ncol,k,lchnk) = vi(:ncol,k,lchnk) - vzm(:ncol,k,lchnk) - wd(:ncol,k,lchnk) = wi(:ncol,k,lchnk) - wzm(:ncol,k,lchnk) + 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) @@ -330,6 +303,11 @@ subroutine phys_grid_ctem_diags(phys_state) 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) @@ -345,11 +323,11 @@ subroutine phys_grid_ctem_diags(phys_state) !------------------------------------------------------------------------------ function zmean_fld( fld ) result(fldzm) - real(r8), intent(in) :: fld(pcols,pverp,begchunk:endchunk) + real(r8), intent(in) :: fld(pcols,pver,begchunk:endchunk) - real(r8) :: fldzm(pcols,pverp,begchunk:endchunk) + real(r8) :: fldzm(pcols,pver,begchunk:endchunk) - real(r8) :: Zonal_Bamp3d(nzmbas,pverp) + real(r8) :: Zonal_Bamp3d(nzmbas,pver) call ZMobj%calc_amps(fld,Zonal_Bamp3d) call ZMobj%eval_grid(Zonal_Bamp3d,fldzm) From b39baa37b245e0a7d15477f1d973827b7a4ae581 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 3 Jan 2023 10:42:39 -0700 Subject: [PATCH 16/24] correction to array ncol extent modified: src/physics/cam/phys_grid_ctem.F90 --- src/physics/cam/phys_grid_ctem.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 5303851c65..95a1808f18 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -265,8 +265,8 @@ subroutine phys_grid_ctem_diags(phys_state) ! vertical velocity w(:ncol,:,lchnk) = -sheight(:ncol,:) * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:) - u(:ncol,:,lchnk) = phys_state(lchnk)%u(:,:) - v(:ncol,:,lchnk) = phys_state(lchnk)%v(:,:) + u(:ncol,:,lchnk) = phys_state(lchnk)%u(:ncol,:) + v(:ncol,:,lchnk) = phys_state(lchnk)%v(:ncol,:) end do From d29da6ef5ee496ebf821b059ea07f284079b0f1d Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 3 Jan 2023 14:44:38 -0700 Subject: [PATCH 17/24] the majority of Jesse's change requests modified: src/physics/cam/phys_grid_ctem.F90 modified: src/utils/zonal_mean_mod.F90 --- src/physics/cam/phys_grid_ctem.F90 | 9 +- src/utils/zonal_mean_mod.F90 | 338 ++++++++++++----------------- 2 files changed, 148 insertions(+), 199 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 95a1808f18..a0a8e84d8f 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -18,6 +18,7 @@ module phys_grid_ctem 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 @@ -83,7 +84,8 @@ subroutine phys_grid_ctem_readnl(nlfile) 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') + 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 @@ -92,7 +94,8 @@ subroutine phys_grid_ctem_readnl(nlfile) ntimesteps = nint( -phys_grid_ctem_nfreq*3600._r8/dtime ) end if if (ntimesteps<1) then - call endrun(prefix//'invalid ntimesteps') + 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 @@ -164,7 +167,7 @@ subroutine phys_grid_ctem_reg ! sanity check if ( abs(1._r8-total_wght)>1.e-12_r8 .or. abs(fourpi-total_area)>1.e-12_r8 ) then - call endrun('zmean_phys_fields_reg: problem with area/wght calc') + call endrun('phys_grid_ctem_reg: problem with area/wght calc') end if ! initialize zonal-average and zonal-mean utility objects diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index e6310b330e..695303b0b2 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -18,7 +18,7 @@ module zonal_mean_mod ! ! NOTE: The weighting of the Zonal Profiles values is scaled such ! that ZonalMean_t amplitudes can be used to evaluate values -! the ZonalProfile_t grid and vise-versa. +! 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 @@ -73,7 +73,7 @@ module zonal_mean_mod ! ------------------------------------------------------ ! - Initialize the data structure for the given number of ! latitudes. Either use the given Latitudes and weights, -! or OPTIONALLY create a profile gridpoints and associated +! 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 @@ -115,7 +115,7 @@ module zonal_mean_mod ! 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/NOT-ACCURATE) on Zonal profile grid +! (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) ! @@ -123,8 +123,8 @@ module zonal_mean_mod ! ========================================= ! call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.) ! -------------------------------------------------- -! - Given the latitude/area for the nlat meridional gridpopints, initialize -! the ZonalAverage datastruture for computing bin-averaging of physgrid +! - Given the latitude/area for the nlat meridional gridpoints, initialize +! the ZonalAverage data struture 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 @@ -161,6 +161,7 @@ module zonal_mean_mod 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 @@ -226,16 +227,19 @@ module zonal_mean_mod procedure, pass :: final => final_ZonalAverage end type ZonalAverage_t - real(r8), parameter :: PI2 = 0.5_r8*pi + 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 = .25_r8*pi contains !======================================================================= subroutine init_ZonalMean(this,I_nbas) ! - ! init_ZonalMean: Initialize the ZonalMean datastrutures for the + ! init_ZonalMean: Initialize the ZonalMean data strutures for the ! physics grid. It is assumed that the domain ! of these gridpoints spans the surface of the sphere. - ! The representation of basis functions functions is + ! The representation of basis functions is ! normalized w.r.t integration over the sphere. !===================================================================== ! @@ -293,7 +297,7 @@ subroutine init_ZonalMean(this,I_nbas) allocate(Clats(pcols,begchunk:endchunk), stat=astat) call handle_allocate_error(astat, subname, 'Clats') - allocate(Bcoef(I_nbas), stat=astat) + 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') @@ -309,7 +313,7 @@ subroutine init_ZonalMean(this,I_nbas) Bsum(:,:) = 0._r8 Csum(:,:) = 0._r8 - ! Save a copy the area weights for each ncol gridpoint + ! Save a copy of the area weights for each ncol gridpoint ! and convert Latitudes to SP->NP colatitudes in radians !------------------------------------------------------- do lchnk=begchunk,endchunk @@ -318,18 +322,13 @@ subroutine init_ZonalMean(this,I_nbas) do cc = 1,ncols rlat=get_rlat_p(lchnk,cc) this%area(cc,lchnk) = area(cc) - Clats (cc,lchnk) = rlat + PI2 + Clats (cc,lchnk) = rlat + halfPI end do end do ! Add first basis for the mean values. !------------------------------------------ - do lchnk=begchunk,endchunk - ncols = get_ncols_p(lchnk) - do cc = 1,ncols - this%basis(cc,lchnk,1) = 1._r8/sqrt(8._r8*PI2) - end do - end do + this%basis(:,begchunk:endchunk,1) = 1._r8/sqrt(fourPI) ! Loop over the remaining basis functions !--------------------------------------- @@ -366,12 +365,10 @@ subroutine init_ZonalMean(this,I_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) - do cc = 1,ncols - this%basis(cc,lchnk,nn) = this%basis(cc,lchnk,nn)/sqrt(Bnorm(nn)) - end do - end do + 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 @@ -379,7 +376,7 @@ subroutine init_ZonalMean(this,I_nbas) !--------------------------------------------------------------- cnum = 0 do nn= 1,this%nbas - do n2=nn,I_nbas + do n2=nn,this%nbas cnum = cnum + 1 count = 0 do lchnk=begchunk,endchunk @@ -387,7 +384,6 @@ subroutine init_ZonalMean(this,I_nbas) 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 @@ -650,7 +646,7 @@ end subroutine final_ZonalMean !======================================================================= subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! - ! init_ZonalProfile: Initialize the ZonalProfile datastruture for the + ! init_ZonalProfile: Initialize the ZonalProfile data struture for the ! given nlat gridpoints. It is assumed that the domain ! of these gridpoints of the profile span latitudes ! from SP to NP. @@ -658,7 +654,7 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! 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 for the ZonalMean_t class. + ! those for the ZonalMean_t class. ! ! The optional GEN_GAUSSLATS flag allows for the ! generation of Gaussian latitudes. The generated grid @@ -709,7 +705,7 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) allocate(Clats(I_nlat), stat=astat) call handle_allocate_error(astat, subname, 'Clats') - allocate(Bcoef(I_nbas), stat=astat) + 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') @@ -731,25 +727,25 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! to degrees and scale the area for global 2D integrals !----------------------------------------------------------- do nn=1,I_nlat - IO_lats(nn) = (45._r8*Clats(nn)/atan(1._r8)) - 90._r8 - IO_area(nn) = IO_area(nn)*(8._r8*atan(1._r8)) + IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 + IO_area(nn) = IO_area(nn)*pi end do else ! Convert Latitudes to SP->NP colatitudes in radians !---------------------------------------------------- do nn=1,I_nlat - Clats(nn) = (IO_lats(nn) + 90._r8)*atan(1._r8)/45._r8 + Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 end do endif ! Copy the area weights for each nlat - ! gridpoint to the datastructure + ! 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) = 1._r8/sqrt(16._r8*atan(1._r8)) + this%basis(:,1) = 1._r8/sqrt(fourPI) Bnorm = 0._r8 do ii=1,I_nlat Bnorm = Bnorm + (this%basis(ii,1)*this%basis(ii,1)*this%area(ii)) @@ -986,7 +982,7 @@ end subroutine final_ZonalProfile !======================================================================= subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) ! - ! init_ZonalAverage: Initialize the ZonalAverage datastruture for the + ! init_ZonalAverage: Initialize the ZonalAverage data struture for the ! given nlat gridpoints. It is assumed that the domain ! of these gridpoints of the profile span latitudes ! from SP to NP. @@ -1072,22 +1068,22 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) ! to degrees and scale the area for global 2D integrals !----------------------------------------------------------- do nn=1,this%nlat - IO_lats(nn) = (45._r8*Clats(nn)/atan(1._r8)) - 90._r8 - IO_area(nn) = IO_area(nn)*(8._r8*atan(1._r8)) + 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)*atan(1._r8)/45._r8 + Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8 end do endif - ! Copy the Lat grid area weights to the datastructure + ! Copy the Lat grid area weights to the data structure !----------------------------------------------------- this%area(1:this%nlat) = IO_area(1:this%nlat) - ! Save a copy the area weights for each 2D gridpoint + ! Save a copy of the area weights for each 2D gridpoint ! and convert Latitudes to SP->NP colatitudes in radians !------------------------------------------------------- do lchnk=begchunk,endchunk @@ -1096,14 +1092,14 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) do cc = 1,ncols rlat=get_rlat_p(lchnk,cc) this%area_g(cc,lchnk) = area(cc) - Glats (cc,lchnk) = rlat + PI2 + Glats (cc,lchnk) = rlat + halfPI end do end do ! Set boundaries for Latitude bins !----------------------------------- BinLat(1) = 0._r8 - BinLat(this%nlat+1) = 4._r8*atan(1._r8) + BinLat(this%nlat+1) = pi do nn=2,this%nlat BinLat(nn) = (Clats(nn-1)+Clats(nn))/2._r8 end do @@ -1114,22 +1110,22 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) ncols = get_ncols_p(lchnk) do cc = 1,ncols jlat = -1 - if((Glats(cc,lchnk).le.BinLat(2)).and. & - (Glats(cc,lchnk).ge.BinLat(1)) ) then + if((Glats(cc,lchnk)<=BinLat(2)).and. & + (Glats(cc,lchnk)>=BinLat(1)) ) then jlat = 1 - elseif((Glats(cc,lchnk).ge.BinLat(this%nlat) ).and. & - (Glats(cc,lchnk).le.BinLat(this%nlat+1)) ) then + 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).gt.BinLat(jj )).and. & - (Glats(cc,lchnk).le.BinLat(jj+1)) ) then + if((Glats(cc,lchnk)>BinLat(jj )).and. & + (Glats(cc,lchnk)<=BinLat(jj+1)) ) then jlat = jj exit endif end do endif - if((jlat.lt.1).or.(jlat.gt.this%nlat)) then + if (jlat<1) then call endrun('ZonalAverage init ERROR: jlat not in range') endif this%idx_map(cc,lchnk) = jlat @@ -1155,8 +1151,13 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) this%a_norm = Anorm if (.not.all(Anorm(:)>0._r8)) then - print*, 'ZonalAverage init ERROR: Anorm; ',Anorm(:) - call endrun('init_ZonalAverage: error in Anorm') + 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 @@ -1422,21 +1423,21 @@ subroutine dalfk(nn,mm,cp) cp(1) = 0._r8 ma = abs(mm) - if(ma.gt.nn) return + if(ma>nn) return - if((nn-1).lt.0) then + if((nn-1)<0) then cp(1) = sqrt(2._r8) return - elseif((nn-1).eq.0) then - if(ma.ne.0) then + elseif((nn-1)==0) then + if(ma/=0) then cp(1) = sqrt(.75_r8) - if(mm.eq.-1) cp(1) = -cp(1) + if(mm==-1) cp(1) = -cp(1) else cp(1) = sqrt(1.5_r8) endif return else - if(mod(nn+ma,2).ne.0) then + if(mod(nn+ma,2)/=0) then nmms2 = (nn-ma-1)/2 fnum = nn + ma + 2 fnmh = nn - ma + 2 @@ -1452,10 +1453,10 @@ subroutine dalfk(nn,mm,cp) t1 = 1._r8/SC20 nex = 20 fden = 2._r8 - if(nmms2.ge.1) then + if(nmms2>=1) then do ii = 1,nmms2 t1 = fnum*t1/fden - if(t1.GT.SC20) THEN + if(t1>SC20) THEN t1 = t1/SC40 nex = nex + 40 endif @@ -1464,13 +1465,13 @@ subroutine dalfk(nn,mm,cp) end do endif - if(mod(ma/2,2).ne.0) then + 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.ne.0) then + if(ma/=0) then do ii = 1,ma t2 = fnmh*t2/ (fnmh+pm1) fnmh = fnmh + 2._r8 @@ -1481,31 +1482,32 @@ subroutine dalfk(nn,mm,cp) fnnp1 = nn*(nn+1) fnmsq = fnnp1 - 2._r8*ma*ma - if((mod(nn,2).eq.0).and.(mod(ma,2).eq.0)) then + 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.lt.0) then - if(mod(ma,2).ne.0) cp(jj) = -cp(jj) + if(mm<0) then + if(mod(ma,2)/=0) cp(jj) = -cp(jj) endif - if(jj.le.1) return + 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 - 30 continue - jj = jj - 1 - if(jj.le.1) return + + 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 - goto 30 + jj = jj - 1 + end do end subroutine dalfk !======================================================================= @@ -1600,20 +1602,20 @@ subroutine dlfpt(nn,mm,theta,cp,pb) pb = 0._r8 ma = abs(mm) - if(ma.gt.nn) return + if(ma>nn) return - if(nn.le.0) then - if(ma.le.0) then + if(nn<=0) then + if(ma<=0) then pb = sqrt(.5_r8) - goto 140 + return endif endif nmod = mod(nn,2) mmod = mod(ma,2) - if(nmod.le.0) then - if(mmod.le.0) then + if(nmod<=0) then + if(mmod<=0) then kdo = nn/2 + 1 cdt = cos(theta+theta) sdt = sin(theta+theta) @@ -1627,7 +1629,7 @@ subroutine dlfpt(nn,mm,theta,cp,pb) summ = summ + cp(kp1)*ct end do pb = summ - goto 140 + return endif kdo = nn/2 cdt = cos(theta+theta) @@ -1642,11 +1644,11 @@ subroutine dlfpt(nn,mm,theta,cp,pb) summ = summ + cp(kk)*st end do pb = summ - goto 140 + return endif kdo = (nn+1)/2 - if(mmod.le.0) then + if(mmod<=0) then cdt = cos(theta+theta) sdt = sin(theta+theta) ct = cos(theta) @@ -1659,7 +1661,7 @@ subroutine dlfpt(nn,mm,theta,cp,pb) summ = summ + cp(kk)*ct end do pb = summ - goto 140 + return endif cdt = cos(theta+theta) @@ -1675,8 +1677,6 @@ subroutine dlfpt(nn,mm,theta,cp,pb) end do pb = summ -140 continue - end subroutine dlfpt !======================================================================= @@ -1722,10 +1722,10 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) do ii=1,Nbas Mmax = 0._r8 do jj=1,Nbas - if(abs(Mwrk(ii,jj)).gt.Mmax) Mmax = abs(Mwrk(ii,jj)) + if(abs(Mwrk(ii,jj))>Mmax) Mmax = abs(Mwrk(ii,jj)) end do - if(Mmax.eq.0._r8) then - call endrun('Singular Matrix') + if(Mmax==0._r8) then + call endrun('Invert_Matrix: Singular Matrix') endif Rscl(ii) = 1._r8/Mmax end do @@ -1734,10 +1734,10 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) !----------------------- do jj=1,Nbas - if(jj.gt.1) then + if(jj>1) then do ii=1,(jj-1) Sval = Mwrk(ii,jj) - if(ii.gt.1) then + if(ii>1) then do kk=1,(ii-1) Sval = Sval - Mwrk(ii,kk)*Mwrk(kk,jj) end do @@ -1749,20 +1749,20 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) Mmax = 0._r8 do ii=jj,Nbas Sval = Mwrk(ii,jj) - if(jj.gt.1) then + 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.ge.Mmax) then + if(Mval>=Mmax) then ii_max = ii Mmax = Mval endif end do - if(jj.ne.ii_max) then + if(jj/=ii_max) then do kk=1,Nbas Mval = Mwrk(ii_max,kk) Mwrk(ii_max,kk) = Mwrk(jj,kk) @@ -1773,8 +1773,8 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) endif Indx(jj) = ii_max - if(jj.ne.Nbas) then - if(Mwrk(jj,jj).eq.0._r8) Mwrk(jj,jj) = TINY + 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 @@ -1783,7 +1783,7 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) end do ! jj=1,Nbas - if(Mwrk(Nbas,Nbas).eq.0._r8) Mwrk(Nbas,Nbas) = TINY + if(Mwrk(Nbas,Nbas)==0._r8) Mwrk(Nbas,Nbas) = TINY ! Initialize the inverse array with the identity matrix !------------------------------------------------------- @@ -1801,11 +1801,11 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) ndx = Indx(ii) Sval = O_InvMat(ndx,kk) O_InvMat(ndx,kk) = O_InvMat(ii,kk) - if(i2.ne.0) then + if(i2/=0) then do jj=i2,(ii-1) Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk) end do - elseif(Sval.ne.0._r8) then + elseif(Sval/=0._r8) then i2 = ii endif O_InvMat(ii,kk) = Sval @@ -1813,7 +1813,7 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) do ii=Nbas,1,-1 Sval = O_InvMat(ii,kk) - if(ii.lt.Nbas) then + if(ii 2 !---------------------- - eps = sqrt(ddzeps(1._r8)) - eps = eps*sqrt(eps) - pis2 = 2._r8*atan(1._r8) - pi = pis2 + pis2 mnlat = mod(nlat,2) ns2 = nlat/2 nhalf = (nlat+1)/2 call dcpdp(nlat,cz,theta(ns2+1),wts(ns2+1)) - dtheta = pis2/nhalf + dtheta = halfPI/nhalf dthalf = dtheta/2._r8 cmax = .2_r8*dtheta ! estimate first point next to theta = pi/2 !------------------------------------------- - if(mnlat.ne.0) then - zero = pis2 - dtheta - zprev = pis2 + if(mnlat/=0) then + zero = halfPI - dtheta + zprev = halfPI nix = nhalf - 1 else - zero = pis2 - dthalf + zero = halfPI - dthalf nix = nhalf endif - 10 continue - it = 0 - 20 continue - it = it + 1 - zlast = zero - - ! newton iterations - !----------------------- - call dtpdp(nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb) + do while(nix/=0) + it = 0 + do while (abs(dcor) > eps*abs(zero)) + it = it + 1 + ! newton iterations + !----------------------- + call dtpdp(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 - 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 - if(abs(zero-zlast).gt.eps*abs(zero)) goto 20 - - 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.eq.0) goto 30 - if(nix.eq.nhalf-1) zero = 3._r8*zero - pi - if(nix.lt.nhalf-1) zero = zero + zero - zprev - zprev = zhold - goto 10 - 30 continue + 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 @@ -2128,50 +2120,4 @@ subroutine dtpdp(nn,theta,cz,cp,dcp,pb,dpb) end subroutine dtpdp !======================================================================= - - !======================================================================= - real(r8) function ddzeps(xx) - ! - ! estimate unit roundoff in quantities of size xx. - ! - ! this program should function properly on all systems - ! satisfying the following two assumptions, - ! 1. the base used in representing floating point - ! numbers is not a power of three. - ! 2. the quantity aa in statement 10 is represented to - ! the accuracy used in floating point variables - ! that are stored in memory. - ! under these assumptions, it should be true that, - ! aa is not exactly equal to four-thirds, - ! bb has a zero for its last bit or digit, - ! cc is not exactly equal to one, - ! eps measures the separation of 1.0 from - ! the next larger floating point number. - ! the developers of eispack would appreciate being informed - ! about any systems where these assumptions do not hold. - ! - !===================================================================== - ! - ! Passed variables - !----------------- - real(r8), intent(in) :: xx - ! - ! Local Values - !-------------- - real(r8):: bb,cc,eps - real(r8), parameter :: aa = 4._r8/3._r8 - - eps = 0.0_r8 - - do while(eps == 0.0_r8) - bb = aa - 1._r8 - cc = bb + bb + bb - eps = abs(cc-1._r8) - end do - - ddzeps = eps*abs(xx) - - end function ddzeps - !======================================================================= - end module zonal_mean_mod From a9a24d7f83d8b5459be90865a384d3dfc359cf16 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 4 Jan 2023 10:34:30 -0700 Subject: [PATCH 18/24] more misc changes suggested by Jesse modified: src/physics/cam/phys_grid_ctem.F90 modified: src/utils/zonal_mean_mod.F90 --- src/physics/cam/phys_grid_ctem.F90 | 2 +- src/utils/zonal_mean_mod.F90 | 512 ++++++++++++++--------------- 2 files changed, 244 insertions(+), 270 deletions(-) diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index a0a8e84d8f..6e2ccf92a9 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -211,7 +211,7 @@ subroutine phys_grid_ctem_init call addfld ('WTHzaphys',(/'lev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') call addfld ('UVzaphys', (/'lev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') call addfld ('UWzaphys', (/'lev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('THphys', (/'lev' /), 'A', 'K', 'Zonal-Mean potential temp - defined on ilev', gridname='physgrid' ) + call addfld ('THphys', (/'lev'/), 'A', 'K', 'Potential temp - defined on ilev', gridname='physgrid' ) end subroutine phys_grid_ctem_init diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 695303b0b2..28511a0c5f 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -337,14 +337,14 @@ subroutine init_ZonalMean(this,I_nbas) ! Generate coefs for the basis !------------------------------ - call dalfk(nb,0,Bcoef) + 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 dlfpt(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn)) + 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 @@ -718,7 +718,7 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! Create a Gaussian grid from SP to NP !-------------------------------------- - call dgaqd(I_nlat,Clats,IO_area,ierr) + 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 @@ -759,13 +759,13 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! Generate coefs for the basis !------------------------------ - call dalfk(nb,0,Bcoef) + 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 dlfpt(nb,0,Clats(ii),Bcoef,this%basis(ii,nn)) + call sh_create_basis(nb,0,Clats(ii),Bcoef,this%basis(ii,nn)) end do ! Numerically normalize the basis funnction @@ -1059,7 +1059,7 @@ subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) ! Create a Gaussin grid from SP to NP !-------------------------------------- - call dgaqd(this%nlat,Clats,IO_area,ierr) + 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 @@ -1311,90 +1311,229 @@ subroutine final_ZonalAverage(this) end subroutine final_ZonalAverage !======================================================================= + !======================================================================= - subroutine dalfk(nn,mm,cp) + 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 ! - ! subroutine alfk (n,m,cp) + ! Passed Variables + !------------------ + real(r8),intent(in ):: I_Mat (:,:) + integer ,intent(in ):: Nbas + real(r8),intent(out):: O_InvMat(:,:) ! - ! dimension of real cp(n/2 + 1) + ! 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(ii=1) then do ii = 1,nmms2 t1 = fnum*t1/fden - if(t1>SC20) THEN + if (t1>SC20) then t1 = t1/SC40 nex = nex + 40 endif @@ -1509,80 +1648,65 @@ subroutine dalfk(nn,mm,cp) jj = jj - 1 end do - end subroutine dalfk + end subroutine sh_gen_basis_coefs !======================================================================= - !======================================================================= - subroutine dlfpt(nn,mm,theta,cp,pb) + subroutine sh_create_basis(nn,mm,theta,cp,pb) ! - ! subroutine lfpt (n,m,theta,cp,pb) + ! spherepack lfpt ! ! dimension of ! arguments - ! cp((n/2)+1) + ! cp((nn/2)+1) ! - ! purpose routine lfpt uses coefficients computed by - ! routine alfk to compute the single precision - ! normalized associated legendre function pbar(n, - ! m,theta) at colatitude theta. - ! - ! usage call lfpt(n,m,theta,cp,pb) + ! 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 n + ! on input nn ! nonnegative integer specifying the degree of - ! pbar(n,m,theta) - ! m - ! is the order of pbar(n,m,theta). m can be - ! any integer however pbar(n,m,theta) = 0 - ! if abs(m) is greater than n and - ! pbar(n,m,theta) = (-1)**m*pbar(n,-m,theta) - ! for negative m. + ! 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 - ! single precision colatitude in radians + ! colatitude in radians ! ! cp - ! single precision array of length (n/2)+1 + ! array of length (nn/2)+1 ! containing coefficients computed by routine - ! alfk + ! sh_gen_basis_coefs ! ! on output pb - ! single precision variable containing - ! pbar(n,m,theta) - ! - ! special conditions calls to routine lfpt must be preceded by an - ! appropriate call to routine alfk. + ! variable containing pbar(n,m,theta) ! - ! precision single + ! 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 lfpt to calculate pbar(n,m,th) at - ! colatitude th depends on m and n as follows: + ! routine sh_create_basis to calculate pbar(nn,mm,theta) at + ! colatitude theta depends on mm and nn as follows: ! - ! 1) for n even and m even, the formula is + ! 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*th) - ! 2) for n even and m odd. the formula is - ! the sum from k=1 to k=n/2 of - ! cp(k)*sin(2*k*th) - ! 3) for n odd and m even, the formula is - ! the sum from k=1 to k=(n+1)/2 of - ! cp(k)*cos((2*k-1)*th) - ! 4) for n odd and m odd, the formula is - ! the sum from k=1 to k=(n+1)/2 of - ! cp(k)*sin((2*k-1)*th) - ! - ! accuracy comparison between routines lfpt and double - ! precision dlfpt on the cray1 indicates greater - ! accuracy for greater values on input parameter - ! n. agreement to 13 places was obtained for - ! n=10 and to 12 places for n=100. - ! - ! timing time per call to routine lfpt is dependent on - ! the input parameter n. + ! 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 @@ -1677,165 +1801,13 @@ subroutine dlfpt(nn,mm,theta,cp,pb) end do pb = summ - end subroutine dlfpt + end subroutine sh_create_basis !======================================================================= - !======================================================================= - 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(ii Date: Wed, 4 Jan 2023 14:14:43 -0700 Subject: [PATCH 19/24] add regression tests modified: cime_config/testdefs/testlist_cam.xml modified: cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam modified: cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm renamed: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands -> cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands renamed: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam -> cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam renamed: cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm -> cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm --- cime_config/testdefs/testlist_cam.xml | 35 +++++++++++++++++++ .../shell_commands | 0 .../user_nl_cam | 4 +-- .../user_nl_clm | 9 +++-- 4 files changed, 41 insertions(+), 7 deletions(-) rename cime_config/testdefs/testmods_dirs/cam/{outfrq9s_physgrid_tem => outfrq3s_physgrid_tem}/shell_commands (100%) rename cime_config/testdefs/testmods_dirs/cam/{outfrq9s_physgrid_tem => outfrq3s_physgrid_tem}/user_nl_cam (78%) rename cime_config/testdefs/testmods_dirs/cam/{outfrq9s_physgrid_tem => outfrq3s_physgrid_tem}/user_nl_clm (94%) diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index fffe69590c..a01ba25b20 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/outfrq9s_physgrid_tem/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands similarity index 100% rename from cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/shell_commands rename to cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/shell_commands diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam similarity index 78% rename from cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam rename to cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam index 2797dc0c3a..55e353983a 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_cam @@ -1,8 +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 +nhtfrq=3,3,3,3,3,3,3,3,3 inithist='ENDOFRUN' -phys_grid_ctem_nfreq=2 +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/outfrq9s_physgrid_tem/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm similarity index 94% rename from cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm rename to cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm index 0d83b5367b..dfd79c6314 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq3s_physgrid_tem/user_nl_clm @@ -1,11 +1,11 @@ !---------------------------------------------------------------------------------- -! Users should add all user specific namelist changes below in the form of -! namelist_var = new_namelist_value +! 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: +! 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 @@ -21,7 +21,6 @@ ! Set maxpatch_glcmec with GLC_NEC option ! Set glc_do_dynglacier with GLC_TWO_WAY_COUPLING env variable !---------------------------------------------------------------------------------- -hist_nhtfrq = 9 +hist_nhtfrq = 3 hist_mfilt = 1 hist_ndens = 1 - From 8927292dc8c5338162393ef6fb6416aeaf4435d3 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Fri, 6 Jan 2023 15:59:41 -0700 Subject: [PATCH 20/24] update to history fields; move mpas tem test to prealpha modified: cime_config/testdefs/testlist_cam.xml modified: src/physics/cam/phys_grid_ctem.F90 --- cime_config/testdefs/testlist_cam.xml | 2 +- src/physics/cam/phys_grid_ctem.F90 | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index a01ba25b20..1a99f5be60 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -754,7 +754,7 @@ - + diff --git a/src/physics/cam/phys_grid_ctem.F90 b/src/physics/cam/phys_grid_ctem.F90 index 6e2ccf92a9..e2d0cb5c71 100644 --- a/src/physics/cam/phys_grid_ctem.F90 +++ b/src/physics/cam/phys_grid_ctem.F90 @@ -207,11 +207,11 @@ subroutine phys_grid_ctem_init if (.not.do_tem_diags) return - call addfld ('VTHzaphys',(/'lev'/), 'A', 'MK/S', 'Meridional Heat Flux:', gridname='ctem_zavg_phys') - call addfld ('WTHzaphys',(/'lev'/), 'A', 'MK/S', 'Vertical Heat Flux:', gridname='ctem_zavg_phys') - call addfld ('UVzaphys', (/'lev'/), 'A', 'M2/S2','Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('UWzaphys', (/'lev'/), 'A', 'M2/S2','Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys') - call addfld ('THphys', (/'lev'/), 'A', 'K', 'Potential temp - defined on ilev', gridname='physgrid' ) + 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 From b2901b0af5098f815487fc3fb4374027568325b3 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 9 Jan 2023 13:33:49 -0700 Subject: [PATCH 21/24] fix 2pi issue; misc clean up modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 28511a0c5f..285b05bcfd 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -124,7 +124,7 @@ module zonal_mean_mod ! call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.) ! -------------------------------------------------- ! - Given the latitude/area for the nlat meridional gridpoints, initialize -! the ZonalAverage data struture for computing bin-averaging of physgrid +! 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 @@ -230,13 +230,14 @@ module zonal_mean_mod 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 = .25_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 strutures for the + ! 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 @@ -328,7 +329,7 @@ subroutine init_ZonalMean(this,I_nbas) ! Add first basis for the mean values. !------------------------------------------ - this%basis(:,begchunk:endchunk,1) = 1._r8/sqrt(fourPI) + this%basis(:,begchunk:endchunk,1) = invSqrt4pi ! Loop over the remaining basis functions !--------------------------------------- @@ -646,7 +647,7 @@ end subroutine final_ZonalMean !======================================================================= subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! - ! init_ZonalProfile: Initialize the ZonalProfile data struture for the + ! 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. @@ -728,7 +729,7 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) !----------------------------------------------------------- do nn=1,I_nlat IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8 - IO_area(nn) = IO_area(nn)*pi + IO_area(nn) = IO_area(nn)*twoPI end do else ! Convert Latitudes to SP->NP colatitudes in radians @@ -745,7 +746,7 @@ subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS) ! Add first basis for the mean values. !------------------------------------------ - this%basis(:,1) = 1._r8/sqrt(fourPI) + 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)) @@ -982,7 +983,7 @@ end subroutine final_ZonalProfile !======================================================================= subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS) ! - ! init_ZonalAverage: Initialize the ZonalAverage data struture for the + ! 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. From 1457f8ed480dd8675e45f24360f5cc4cf7c54231 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Mon, 9 Jan 2023 13:42:26 -0700 Subject: [PATCH 22/24] initialize dcor in sh_create_gaus_grid modified: src/utils/zonal_mean_mod.F90 --- src/utils/zonal_mean_mod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 285b05bcfd..3c41a4b4ec 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -1914,6 +1914,7 @@ subroutine sh_create_gaus_grid(nlat,theta,wts,ierr) endif do while(nix/=0) + dcor = huge(1._r8) it = 0 do while (abs(dcor) > eps*abs(zero)) it = it + 1 From 903c2df0ab26097cfd2a389a759a07cc1673a251 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 10 Jan 2023 09:56:33 -0700 Subject: [PATCH 23/24] ChangeLog draft --- doc/ChangeLog | 119 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index 6342032df6..7776043007 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,124 @@ =============================================================== +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 +M bld/namelist_files/use_cases/waccm_sc_2000_cam6.xml +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: + +izumi/nag/aux_cam: + +izumi/gnu/aux_cam: + +CAM tag used for the baseline comparison tests if different than previous +tag: + +Summarize any changes to answers, i.e., +- what code configurations: +- what platforms/compilers: +- nature of change (roundoff; larger than roundoff but same climate; new + climate): + +If bitwise differences were observed, how did you show they were no worse +than roundoff? + +If this tag changes climate describe the run(s) done to evaluate the new +climate in enough detail that it(they) could be reproduced, i.e., +- source tag (all code used must be in the repository): +- platform/compilers: +- configure commandline: +- build-namelist command (or complete namelist): +- MSS location of output: + +MSS location of control simulations used to validate new climate: + +URL for AMWG diagnostics output used to validate new climate: + +=============================================================== +=============================================================== + Tag name: cam6_3_088 Originator(s): cacraig Date: Dec 21, 2022 From 479bf6062f8a27413a23ef7e849c138b539f5be5 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 10 Jan 2023 11:40:06 -0700 Subject: [PATCH 24/24] ChangeLog updates --- doc/ChangeLog | 38 ++++++++++++++------------------------ 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/doc/ChangeLog b/doc/ChangeLog index 7776043007..ea503ba1ce 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -65,7 +65,14 @@ A cime_config/testdefs/testmods_dirs/cam/outfrq9s_physgrid_tem_1deg/user_nl_clm 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 @@ -86,35 +93,18 @@ 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: +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 -izumi/gnu/aux_cam: - -CAM tag used for the baseline comparison tests if different than previous -tag: - -Summarize any changes to answers, i.e., -- what code configurations: -- what platforms/compilers: -- nature of change (roundoff; larger than roundoff but same climate; new - climate): - -If bitwise differences were observed, how did you show they were no worse -than roundoff? - -If this tag changes climate describe the run(s) done to evaluate the new -climate in enough detail that it(they) could be reproduced, i.e., -- source tag (all code used must be in the repository): -- platform/compilers: -- configure commandline: -- build-namelist command (or complete namelist): -- MSS location of output: + SMS_D_Ln6_Vnuopc.ne5_ne5_mg37.QPWmaC4.izumi_nag.cam-outfrq3s_physgrid_tem (Overall: DIFF) details: + - new test -MSS location of control simulations used to validate new climate: +izumi/gnu/aux_cam: All PASS -URL for AMWG diagnostics output used to validate new climate: +Summarize any changes to answers: bit-for-bit unchanged =============================================================== ===============================================================